Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Control for reading in different topologies and coordinates
10 : !> \par History
11 : !> none
12 : ! **************************************************************************************************
13 : MODULE topology
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE atoms_input, ONLY: read_atoms_input
16 : USE cell_methods, ONLY: cell_create,&
17 : read_cell,&
18 : write_cell
19 : USE cell_types, ONLY: cell_clone,&
20 : cell_release,&
21 : cell_retain,&
22 : cell_transform_input_cartesian,&
23 : cell_type
24 : USE colvar_types, ONLY: colvar_p_type,&
25 : colvar_setup,&
26 : combine_colvar_id
27 : USE colvar_utils, ONLY: post_process_colvar
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type
30 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31 : cp_print_key_unit_nr
32 : USE cp_subsys_types, ONLY: cp_subsys_set,&
33 : cp_subsys_type
34 : USE exclusion_types, ONLY: exclusion_type
35 : USE input_constants, ONLY: &
36 : do_conn_amb7, do_conn_g87, do_conn_g96, do_conn_generate, do_conn_mol_set, do_conn_off, &
37 : do_conn_psf, do_conn_psf_u, do_conn_user, do_coord_cif, do_coord_cp2k, do_coord_crd, &
38 : do_coord_g96, do_coord_off, do_coord_pdb, do_coord_xtl, do_coord_xyz
39 : USE input_cp2k_binary_restarts, ONLY: read_binary_coordinates
40 : USE input_section_types, ONLY: section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get,&
44 : section_vals_val_set
45 : USE kinds, ONLY: default_path_length,&
46 : default_string_length,&
47 : dp
48 : USE message_passing, ONLY: mp_para_env_type
49 : USE mm_mapping_library, ONLY: create_ff_map,&
50 : destroy_ff_map
51 : USE molecule_kind_types, ONLY: molecule_kind_type
52 : USE molecule_types, ONLY: global_constraint_type,&
53 : molecule_type
54 : USE particle_types, ONLY: particle_type
55 : USE qmmm_topology_util, ONLY: qmmm_connectivity_control,&
56 : qmmm_coordinate_control
57 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
58 : USE string_table, ONLY: id2str,&
59 : s2s,&
60 : str2id
61 : USE topology_amber, ONLY: read_connectivity_amber,&
62 : read_coordinate_crd
63 : USE topology_cif, ONLY: read_coordinate_cif
64 : USE topology_connectivity_util, ONLY: topology_conn_multiple,&
65 : topology_connectivity_pack
66 : USE topology_constraint_util, ONLY: topology_constraint_pack
67 : USE topology_coordinate_util, ONLY: topology_coordinate_pack
68 : USE topology_cp2k, ONLY: read_coordinate_cp2k
69 : USE topology_generate_util, ONLY: topology_generate_bend,&
70 : topology_generate_bond,&
71 : topology_generate_dihe,&
72 : topology_generate_impr,&
73 : topology_generate_molecule,&
74 : topology_generate_onfo,&
75 : topology_generate_ub
76 : USE topology_gromos, ONLY: read_coordinate_g96,&
77 : read_topology_gromos
78 : USE topology_input, ONLY: read_constraints_section,&
79 : read_topology_section
80 : USE topology_multiple_unit_cell, ONLY: topology_muc
81 : USE topology_pdb, ONLY: read_coordinate_pdb,&
82 : write_coordinate_pdb
83 : USE topology_psf, ONLY: idm_psf,&
84 : psf_post_process,&
85 : read_topology_psf,&
86 : write_topology_psf
87 : USE topology_types, ONLY: deallocate_topology,&
88 : init_topology,&
89 : pre_read_topology,&
90 : topology_parameters_type
91 : USE topology_util, ONLY: check_subsys_element,&
92 : topology_molecules_check,&
93 : topology_reorder_atoms,&
94 : topology_set_atm_mass
95 : USE topology_xtl, ONLY: read_coordinate_xtl
96 : USE topology_xyz, ONLY: read_coordinate_xyz
97 : #include "./base/base_uses.f90"
98 :
99 : IMPLICIT NONE
100 :
101 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology'
102 :
103 : PRIVATE
104 :
105 : ! Public parameters
106 : PUBLIC :: topology_control, &
107 : connectivity_control
108 :
109 : CONTAINS
110 :
111 : ! **************************************************************************************************
112 : !> \brief ...
113 : !> \param atomic_kind_set ...
114 : !> \param particle_set ...
115 : !> \param molecule_kind_set ...
116 : !> \param molecule_set ...
117 : !> \param colvar_p ...
118 : !> \param gci ...
119 : !> \param root_section ...
120 : !> \param para_env ...
121 : !> \param qmmm ...
122 : !> \param qmmm_env ...
123 : !> \param force_env_section ...
124 : !> \param subsys_section ...
125 : !> \param use_motion_section ...
126 : !> \param exclusions ...
127 : !> \param elkind ...
128 : !> \param subsys ...
129 : ! **************************************************************************************************
130 53630 : SUBROUTINE topology_control(atomic_kind_set, particle_set, molecule_kind_set, &
131 : molecule_set, colvar_p, gci, root_section, para_env, qmmm, qmmm_env, &
132 : force_env_section, subsys_section, use_motion_section, &
133 : exclusions, elkind, subsys)
134 :
135 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
136 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
137 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
138 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
139 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
140 : TYPE(global_constraint_type), POINTER :: gci
141 : TYPE(section_vals_type), POINTER :: root_section
142 : TYPE(mp_para_env_type), POINTER :: para_env
143 : LOGICAL, INTENT(IN), OPTIONAL :: qmmm
144 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
145 : TYPE(section_vals_type), POINTER :: force_env_section, subsys_section
146 : LOGICAL, INTENT(IN) :: use_motion_section
147 : TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
148 : POINTER :: exclusions
149 : LOGICAL, INTENT(IN), OPTIONAL :: elkind
150 : TYPE(cp_subsys_type), OPTIONAL, POINTER :: subsys
151 :
152 : CHARACTER(LEN=*), PARAMETER :: routineN = 'topology_control'
153 :
154 : INTEGER :: handle, iw, iw2
155 : LOGICAL :: binary_coord_read, el_as_kind, explicit, &
156 : my_qmmm, use_ref_cell
157 : TYPE(cell_type), POINTER :: my_cell, my_cell_ref
158 : TYPE(cp_logger_type), POINTER :: logger
159 : TYPE(section_vals_type), POINTER :: cell_ref_section, cell_section, &
160 : constraint_section, fragment_section, &
161 : topology_section
162 : TYPE(topology_parameters_type) :: topology
163 :
164 10726 : NULLIFY (logger)
165 21452 : logger => cp_get_default_logger()
166 10726 : CALL timeset(routineN, handle)
167 10726 : NULLIFY (cell_section, cell_ref_section, constraint_section, topology_section, fragment_section)
168 :
169 10726 : cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
170 10726 : IF (use_motion_section) THEN
171 10391 : constraint_section => section_vals_get_subs_vals(root_section, "MOTION%CONSTRAINT")
172 : END IF
173 10726 : topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
174 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
175 10726 : extension=".mmLog")
176 10726 : my_qmmm = .FALSE.
177 10726 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
178 :
179 10726 : IF (PRESENT(elkind)) THEN
180 7886 : CALL section_vals_val_get(topology_section, "USE_ELEMENT_AS_KIND", explicit=explicit)
181 7886 : IF (explicit) THEN
182 0 : CALL section_vals_val_get(topology_section, "USE_ELEMENT_AS_KIND", l_val=el_as_kind)
183 : ELSE
184 7886 : el_as_kind = elkind
185 : END IF
186 : ELSE
187 2840 : CALL section_vals_val_get(topology_section, "USE_ELEMENT_AS_KIND", l_val=el_as_kind)
188 : END IF
189 :
190 : ! 1. Initialize the topology structure type
191 10726 : CALL init_topology(topology)
192 :
193 : ! 2. Get the cell info
194 10726 : NULLIFY (my_cell, my_cell_ref)
195 : CALL read_cell(my_cell, my_cell_ref, cell_section=cell_section, &
196 10726 : topology_section=topology_section, para_env=para_env)
197 : use_ref_cell = .FALSE.
198 10726 : cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
199 10726 : CALL section_vals_get(cell_ref_section, explicit=use_ref_cell)
200 10726 : IF (.NOT. ASSOCIATED(topology%cell)) CALL cell_create(topology%cell)
201 10726 : IF (.NOT. ASSOCIATED(topology%cell_ref)) CALL cell_create(topology%cell_ref)
202 10726 : CALL cell_clone(my_cell, topology%cell)
203 10726 : CALL cell_clone(my_cell_ref, topology%cell_ref)
204 10726 : CALL setup_cell_muc(topology%cell_muc, topology%cell, subsys_section)
205 10726 : IF (PRESENT(subsys)) THEN ! subsys is passed from cp_subsys_create
206 : CALL cp_subsys_set(subsys, cell=my_cell, cell_ref=my_cell_ref, &
207 10723 : use_ref_cell=use_ref_cell)
208 : END IF
209 10726 : CALL write_cell(my_cell, subsys_section, tag="CELL")
210 10726 : CALL write_cell(my_cell_ref, subsys_section, tag="CELL_REF")
211 10726 : CALL cell_release(my_cell)
212 10726 : CALL cell_release(my_cell_ref)
213 :
214 : ! 3. Read in the topology section in the input file if any
215 10726 : CALL read_topology_section(topology, topology_section)
216 :
217 : ! 4. Read in the constraints section
218 10726 : CALL read_constraints_section(topology, colvar_p, constraint_section)
219 :
220 : ! 5. Read in the coordinates
221 : CALL read_binary_coordinates(topology, root_section, para_env, subsys_section, &
222 10726 : binary_coord_read)
223 10726 : IF (.NOT. binary_coord_read) THEN
224 10680 : CALL coordinate_control(topology, root_section, para_env, subsys_section)
225 : END IF
226 :
227 : ! 6. Read in or generate the molecular connectivity
228 : CALL connectivity_control(topology, para_env, my_qmmm, qmmm_env, subsys_section, &
229 10726 : force_env_section)
230 :
231 10726 : IF (el_as_kind) THEN
232 : ! redefine atom names with the name of the element
233 26284 : topology%atom_info%id_atmname(:) = topology%atom_info%id_element(:)
234 : END IF
235 :
236 : ! 7. Pack everything into the molecular types
237 : CALL topology_connectivity_pack(molecule_kind_set, molecule_set, &
238 10726 : topology, subsys_section)
239 :
240 : ! 8. Set up the QM/MM linkage (if any)
241 : ! This part takes care of the molecule in which QM atoms were defined.
242 : ! Preliminary setup for QM/MM link region
243 10726 : IF (my_qmmm) THEN
244 394 : CALL qmmm_connectivity_control(molecule_set, qmmm_env, subsys_section)
245 : END IF
246 :
247 : ! 9. Pack everything into the atomic types
248 10726 : IF (my_qmmm) THEN
249 : CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
250 : molecule_kind_set, molecule_set, &
251 : topology, my_qmmm, qmmm_env, subsys_section, &
252 394 : force_env_section=force_env_section, exclusions=exclusions)
253 : ELSE
254 : CALL topology_coordinate_pack(particle_set, atomic_kind_set, &
255 : molecule_kind_set, molecule_set, &
256 : topology, subsys_section=subsys_section, &
257 10332 : force_env_section=force_env_section, exclusions=exclusions)
258 : END IF
259 :
260 : !10. Post-Process colvar definitions (if needed)
261 10726 : CALL topology_post_proc_colvar(colvar_p, particle_set)
262 :
263 : !11. Deal with the constraint stuff if requested
264 10726 : IF (my_qmmm) THEN
265 : CALL topology_constraint_pack(molecule_kind_set, molecule_set, &
266 : topology, qmmm_env, particle_set, root_section, subsys_section, &
267 394 : gci)
268 : ELSE
269 : CALL topology_constraint_pack(molecule_kind_set, molecule_set, &
270 : topology, particle_set=particle_set, input_file=root_section, &
271 10332 : subsys_section=subsys_section, gci=gci)
272 : END IF
273 :
274 : !12. Dump the topology informations
275 : iw2 = cp_print_key_unit_nr(logger, subsys_section, "TOPOLOGY%DUMP_PDB", &
276 10726 : file_status="REPLACE", extension=".pdb")
277 10726 : IF (iw2 > 0) THEN
278 54 : CALL write_coordinate_pdb(iw2, topology, subsys_section)
279 : END IF
280 : CALL cp_print_key_finished_output(iw2, logger, subsys_section, &
281 10726 : "TOPOLOGY%DUMP_PDB")
282 : iw2 = cp_print_key_unit_nr(logger, subsys_section, "TOPOLOGY%DUMP_PSF", &
283 10726 : file_status="REPLACE", extension=".psf")
284 10726 : IF (iw2 > 0) THEN
285 56 : CALL write_topology_psf(iw2, topology, subsys_section, force_env_section)
286 : END IF
287 : CALL cp_print_key_finished_output(iw2, logger, subsys_section, &
288 10726 : "TOPOLOGY%DUMP_PSF")
289 : !13. Read fragmentation
290 10726 : fragment_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%FRAGMENTS")
291 10726 : CALL section_vals_get(fragment_section, explicit=explicit)
292 10726 : IF (explicit) THEN
293 2 : CALL read_fragments(topology, subsys_section, particle_set=particle_set)
294 : END IF
295 :
296 : !13. Cleanup the topology structure type
297 10726 : CALL deallocate_topology(topology)
298 10726 : CALL timestop(handle)
299 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
300 10726 : "PRINT%TOPOLOGY_INFO")
301 :
302 10726 : END SUBROUTINE topology_control
303 :
304 : ! **************************************************************************************************
305 : !> \brief 1. If reading in from external file, make sure its there first
306 : !> 2. Generate the connectivity if no information to be read in
307 : !> \param topology ...
308 : !> \param para_env ...
309 : !> \param qmmm ...
310 : !> \param qmmm_env ...
311 : !> \param subsys_section ...
312 : !> \param force_env_section ...
313 : !> \par History
314 : !> none
315 : !> \author IKUO 08.01.2003
316 : ! **************************************************************************************************
317 11264 : SUBROUTINE connectivity_control(topology, para_env, qmmm, qmmm_env, subsys_section, &
318 : force_env_section)
319 :
320 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
321 : TYPE(mp_para_env_type), POINTER :: para_env
322 : LOGICAL, INTENT(in), OPTIONAL :: qmmm
323 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
324 : TYPE(section_vals_type), POINTER :: subsys_section, force_env_section
325 :
326 : CHARACTER(len=*), PARAMETER :: routineN = 'connectivity_control'
327 : INTEGER, PARAMETER :: map0 = ICHAR("0"), map9 = ICHAR("9")
328 :
329 : CHARACTER(len=default_string_length) :: element0, my_element
330 : CHARACTER(len=default_string_length), &
331 11264 : ALLOCATABLE, DIMENSION(:) :: elements
332 : INTEGER :: handle, handle2, i, id, itmp, iw, j, k
333 : LOGICAL :: check, my_qmmm, use_mm_map_first
334 : TYPE(cp_logger_type), POINTER :: logger
335 :
336 11264 : NULLIFY (logger)
337 11264 : logger => cp_get_default_logger()
338 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
339 11264 : extension=".mmLog")
340 11264 : CALL timeset(routineN, handle)
341 :
342 11264 : my_qmmm = .FALSE.
343 11264 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
344 :
345 : ! 1. Read in the connectivity information (if this is the case)
346 12061 : SELECT CASE (topology%conn_type)
347 : CASE (do_conn_generate, do_conn_off, do_conn_user)
348 : ! Do nothing for the time being.. after we check element and proceed with the workflow..
349 : CASE DEFAULT
350 : ! Prepare arrays
351 797 : CALL pre_read_topology(topology)
352 :
353 : ! Read connectivity from file
354 : CALL read_topology_conn(topology, topology%conn_type, topology%conn_file_name, &
355 797 : para_env, subsys_section)
356 :
357 : ! Post process of PSF and AMBER information
358 11264 : SELECT CASE (topology%conn_type)
359 : CASE (do_conn_mol_set, do_conn_psf, do_conn_psf_u, do_conn_amb7)
360 797 : CALL psf_post_process(topology, subsys_section)
361 : END SELECT
362 : END SELECT
363 :
364 : ! 2. In case element was autoassigned let's keep up2date the element name
365 : ! with the atom_name
366 11264 : IF (topology%aa_element) THEN
367 10072 : check = SIZE(topology%atom_info%id_element) == SIZE(topology%atom_info%id_atmname)
368 10072 : CPASSERT(check)
369 595008 : topology%atom_info%id_element = topology%atom_info%id_atmname
370 : END IF
371 :
372 : ! 3. Check for the element name..
373 11264 : CALL timeset(routineN//"_check_element_name", handle2)
374 : ! Fix element name
375 :
376 : ! we will only translate names if we actually have a connectivity file given
377 12061 : SELECT CASE (topology%conn_type)
378 : CASE (do_conn_mol_set, do_conn_psf, do_conn_psf_u, do_conn_g96, do_conn_g87, do_conn_amb7)
379 797 : use_mm_map_first = .TRUE.
380 : CASE DEFAULT
381 11264 : use_mm_map_first = .FALSE.
382 : END SELECT
383 11264 : CALL create_ff_map("AMBER")
384 11264 : CALL create_ff_map("CHARMM")
385 11264 : CALL create_ff_map("GROMOS")
386 :
387 33792 : ALLOCATE (elements(SIZE(topology%atom_info%id_element)))
388 745069 : DO i = 1, SIZE(elements)
389 745069 : elements(i) = id2str(topology%atom_info%id_element(i))
390 : END DO
391 :
392 745069 : DO i = 1, topology%natoms
393 733805 : IF (elements(i) == "__DEF__") CYCLE
394 : ! If present an underscore let's skip all that over the underscore
395 25391 : id = INDEX(elements(i), "_") - 1
396 25391 : IF (id == -1) id = LEN_TRIM(elements(i))
397 : ! Many atomic kind have been defined as ELEMENT+LETTER+NUMBER
398 : ! the number at the end can vary arbitrarily..
399 : ! Let's check all ELEMENT+LETTER skipping the number.. we should
400 : ! identify in any case the element
401 29859 : DO j = id, 1, -1
402 29847 : itmp = ICHAR(elements(i) (j:j))
403 29859 : IF ((itmp < map0) .OR. (itmp > map9)) EXIT
404 : END DO
405 25391 : element0 = elements(i) (1:j)
406 : ! ALWAYS check for elements..
407 : CALL check_subsys_element(element0, id2str(topology%atom_info%id_atmname(i)), my_element, &
408 25391 : subsys_section, use_mm_map_first)
409 : ! Earn time fixing same element labels for same atoms
410 25391 : element0 = elements(i)
411 25270860 : DO k = i, topology%natoms
412 25968010 : IF (element0 == id2str(topology%atom_info%id_element(k))) THEN
413 743491 : topology%atom_info%id_element(k) = str2id(s2s(my_element))
414 25977696 : elements(k) = "__DEF__"
415 : END IF
416 : END DO
417 : END DO
418 11264 : DEALLOCATE (elements)
419 11264 : CALL destroy_ff_map("GROMOS")
420 11264 : CALL destroy_ff_map("CHARMM")
421 11264 : CALL destroy_ff_map("AMBER")
422 11264 : CALL timestop(handle2)
423 :
424 : ! 4. Generate the connectivity information otherwise
425 20029 : SELECT CASE (topology%conn_type)
426 : CASE (do_conn_generate)
427 8765 : CALL topology_set_atm_mass(topology, subsys_section)
428 8765 : CALL topology_generate_bond(topology, para_env, subsys_section)
429 8765 : IF (topology%reorder_atom) THEN
430 : ! If we generate connectivity we can save memory reordering the molecules
431 : ! in this case once a first connectivity has been created we match according
432 : ! molecule names provided in the PDB and reorder the connectivity according to that.
433 : CALL topology_reorder_atoms(topology, qmmm, qmmm_env, subsys_section, &
434 316 : force_env_section)
435 316 : CALL topology_set_atm_mass(topology, subsys_section)
436 316 : CALL topology_generate_bond(topology, para_env, subsys_section)
437 : END IF
438 8765 : CALL topology_generate_bend(topology, subsys_section)
439 8765 : CALL topology_generate_ub(topology, subsys_section)
440 8765 : CALL topology_generate_dihe(topology, subsys_section)
441 8765 : CALL topology_generate_impr(topology, subsys_section)
442 8765 : CALL topology_generate_onfo(topology, subsys_section)
443 : CASE (do_conn_off, do_conn_user)
444 1702 : CALL topology_set_atm_mass(topology, subsys_section)
445 1702 : CALL topology_generate_bend(topology, subsys_section)
446 1702 : CALL topology_generate_ub(topology, subsys_section)
447 1702 : CALL topology_generate_dihe(topology, subsys_section)
448 1702 : CALL topology_generate_impr(topology, subsys_section)
449 12966 : CALL topology_generate_onfo(topology, subsys_section)
450 : END SELECT
451 :
452 : ! 5. Handle multiple unit_cell - Update atoms_info
453 11264 : CALL topology_muc(topology, subsys_section)
454 :
455 : ! 6. Handle multiple unit_cell - Update conn_info
456 11264 : CALL topology_conn_multiple(topology, subsys_section)
457 :
458 : ! 7. Generate Molecules
459 11264 : CALL topology_generate_molecule(topology, my_qmmm, qmmm_env, subsys_section)
460 11264 : IF (topology%molecules_check) CALL topology_molecules_check(topology, subsys_section)
461 :
462 : ! 8. Modify for QM/MM
463 11264 : IF (my_qmmm) THEN
464 394 : CALL qmmm_coordinate_control(topology, qmmm_env, subsys_section)
465 : END IF
466 11264 : CALL timestop(handle)
467 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
468 11264 : "PRINT%TOPOLOGY_INFO")
469 :
470 22528 : END SUBROUTINE connectivity_control
471 : ! **************************************************************************************************
472 : !> \brief Reads fragments from input file
473 : !> \param topology ...
474 : !> \param subsys_section ...
475 : !> \param particle_set ...
476 : ! **************************************************************************************************
477 2 : SUBROUTINE read_fragments(topology, subsys_section, particle_set)
478 :
479 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
480 : TYPE(section_vals_type), POINTER :: subsys_section
481 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
482 :
483 : INTEGER :: i_rep, ifrag, n_rep
484 2 : INTEGER, DIMENSION(:), POINTER :: frag_range
485 : TYPE(section_vals_type), POINTER :: section
486 :
487 2 : NULLIFY (section)
488 2 : NULLIFY (frag_range)
489 :
490 4 : section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%FRAGMENTS")
491 2 : section => section_vals_get_subs_vals(section, "FRAGMENT")
492 2 : CALL section_vals_get(section, n_repetition=n_rep)
493 10 : ALLOCATE (topology%fragments(n_rep))
494 :
495 6 : DO i_rep = 1, n_rep
496 4 : CALL section_vals_val_get(section, "FRAG_RANGE", i_vals=frag_range, i_rep_section=i_rep)
497 4 : topology%fragments(i_rep)%frag_start = frag_range(1)
498 4 : topology%fragments(i_rep)%frag_end = frag_range(2)
499 20 : DO ifrag = topology%fragments(i_rep)%frag_start, topology%fragments(i_rep)%frag_end
500 18 : particle_set(ifrag)%fragment_index = i_rep
501 : END DO
502 : END DO
503 :
504 4 : END SUBROUTINE read_fragments
505 :
506 : ! **************************************************************************************************
507 : !> \brief Reads connectivity from file
508 : !> \param topology ...
509 : !> \param conn_type ...
510 : !> \param conn_file_name ...
511 : !> \param para_env ...
512 : !> \param subsys_section ...
513 : !> \author Teodoro Laino [tlaino] - 10.2009
514 : ! **************************************************************************************************
515 8629 : RECURSIVE SUBROUTINE read_topology_conn(topology, conn_type, conn_file_name, para_env, &
516 : subsys_section)
517 :
518 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
519 : INTEGER, INTENT(IN) :: conn_type
520 : CHARACTER(LEN=default_path_length), INTENT(IN) :: conn_file_name
521 : TYPE(mp_para_env_type), POINTER :: para_env
522 : TYPE(section_vals_type), POINTER :: subsys_section
523 :
524 : CHARACTER(len=default_path_length) :: filename
525 : INTEGER :: i_rep, imol, loc_conn_type, n_rep, nmol
526 : TYPE(section_vals_type), POINTER :: section
527 :
528 8629 : NULLIFY (section)
529 :
530 8902 : SELECT CASE (conn_type)
531 : CASE (do_conn_mol_set)
532 273 : section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%MOL_SET")
533 273 : section => section_vals_get_subs_vals(section, "MOLECULE")
534 273 : CALL section_vals_get(section, n_repetition=n_rep)
535 776 : DO i_rep = 1, n_rep
536 503 : CALL section_vals_val_get(section, "NMOL", i_val=nmol, i_rep_section=i_rep)
537 503 : CALL section_vals_val_get(section, "CONN_FILE_NAME", c_val=filename, i_rep_section=i_rep)
538 503 : CALL section_vals_val_get(section, "CONN_FILE_FORMAT", i_val=loc_conn_type, i_rep_section=i_rep)
539 :
540 1279 : SELECT CASE (loc_conn_type)
541 : CASE (do_conn_psf, do_conn_psf_u, do_conn_g96, do_conn_g87, do_conn_amb7)
542 8335 : DO imol = 1, nmol
543 8335 : CALL read_topology_conn(topology, loc_conn_type, filename, para_env, subsys_section)
544 : END DO
545 : CASE DEFAULT
546 : CALL cp_abort(__LOCATION__, &
547 : "MOL_SET feature implemented only for PSF/UPSF, G87/G96 and AMBER "// &
548 503 : "connectivity type.")
549 : END SELECT
550 : END DO
551 273 : IF (SIZE(topology%atom_info%id_molname) /= topology%natoms) &
552 : CALL cp_abort(__LOCATION__, &
553 : "Number of atoms in connectivity control is larger than the "// &
554 : "number of atoms in coordinate control. check coordinates and "// &
555 0 : "connectivity. ")
556 :
557 : ! Merge defined structures
558 273 : section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%MOL_SET%MERGE_MOLECULES")
559 273 : CALL idm_psf(topology, section, subsys_section)
560 :
561 : CASE (do_conn_g96, do_conn_g87)
562 14 : CALL read_topology_gromos(conn_file_name, topology, para_env, subsys_section)
563 : CASE (do_conn_psf, do_conn_psf_u)
564 8320 : CALL read_topology_psf(conn_file_name, topology, para_env, subsys_section, conn_type)
565 : CASE (do_conn_amb7)
566 8902 : CALL read_connectivity_amber(conn_file_name, topology, para_env, subsys_section)
567 : END SELECT
568 :
569 8629 : END SUBROUTINE read_topology_conn
570 :
571 : ! **************************************************************************************************
572 : !> \brief 1. If reading in from external file, make sure its there first
573 : !> 2. Read in the coordinates from the corresponding locations
574 : !> \param topology ...
575 : !> \param root_section ...
576 : !> \param para_env ...
577 : !> \param subsys_section ...
578 : !> \par History
579 : !> - Teodoro Laino [tlaino] - University of Zurich 10.2008
580 : !> adding support for AMBER coordinates
581 : !> \author IKUO 08.11.2003
582 : ! **************************************************************************************************
583 42720 : SUBROUTINE coordinate_control(topology, root_section, para_env, subsys_section)
584 :
585 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
586 : TYPE(section_vals_type), POINTER :: root_section
587 : TYPE(mp_para_env_type), POINTER :: para_env
588 : TYPE(section_vals_type), POINTER :: subsys_section
589 :
590 : CHARACTER(len=*), PARAMETER :: routineN = 'coordinate_control'
591 :
592 : CHARACTER(LEN=default_string_length) :: message
593 : INTEGER :: handle, handle2, istat, iw
594 : LOGICAL :: found, save_mem
595 : TYPE(cp_logger_type), POINTER :: logger
596 : TYPE(section_vals_type), POINTER :: global_section
597 :
598 10680 : NULLIFY (logger)
599 10680 : logger => cp_get_default_logger()
600 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
601 10680 : extension=".mmLog")
602 10680 : CALL timeset(routineN, handle)
603 :
604 10680 : NULLIFY (global_section)
605 10680 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
606 10680 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
607 : !-----------------------------------------------------------------------------
608 : ! 1. If reading in from external file, make sure its there first
609 : !-----------------------------------------------------------------------------
610 10680 : IF (topology%coordinate) THEN
611 1921 : INQUIRE (FILE=topology%coord_file_name, EXIST=found, IOSTAT=istat)
612 1921 : IF (istat /= 0) THEN
613 : WRITE (UNIT=message, FMT="(A,I0,A)") &
614 : "An error occurred inquiring the file <"// &
615 0 : TRIM(topology%coord_file_name)//"> (IOSTAT = ", istat, ")"
616 0 : CPABORT(TRIM(message))
617 : END IF
618 1921 : IF (.NOT. found) THEN
619 : CALL cp_abort(__LOCATION__, &
620 : "Coordinate file <"//TRIM(topology%coord_file_name)// &
621 0 : "> not found.")
622 : END IF
623 : END IF
624 : !-----------------------------------------------------------------------------
625 : ! 2. Read in the coordinates from the corresponding locations
626 : !-----------------------------------------------------------------------------
627 10680 : CALL timeset(routineN//"_READ_COORDINATE", handle2)
628 10694 : SELECT CASE (topology%coord_type)
629 : CASE (do_coord_off)
630 : ! Do nothing.. we will parse later from the &COORD section..
631 : CASE (do_coord_g96)
632 14 : CALL read_coordinate_g96(topology, para_env, subsys_section)
633 14 : CALL transform_external_cartesian_coordinates(topology)
634 : CASE (do_coord_crd)
635 26 : CALL read_coordinate_crd(topology, para_env, subsys_section)
636 26 : CALL transform_external_cartesian_coordinates(topology)
637 : CASE (do_coord_pdb)
638 768 : CALL read_coordinate_pdb(topology, para_env, subsys_section)
639 768 : CALL transform_external_cartesian_coordinates(topology)
640 : CASE (do_coord_xyz)
641 1085 : CALL read_coordinate_xyz(topology, para_env, subsys_section)
642 1085 : CALL transform_external_cartesian_coordinates(topology)
643 : CASE (do_coord_cif)
644 14 : CALL read_coordinate_cif(topology, para_env, subsys_section)
645 : CASE (do_coord_xtl)
646 2 : CALL read_coordinate_xtl(topology, para_env, subsys_section)
647 : CASE (do_coord_cp2k)
648 12 : CALL read_coordinate_cp2k(topology, para_env, subsys_section)
649 : CASE DEFAULT
650 : ! We should never reach this point..
651 10680 : CPABORT("")
652 : END SELECT
653 :
654 : ! Parse &COORD section and in case overwrite
655 10680 : IF (topology%coord_type /= do_coord_cp2k) THEN
656 : CALL read_atoms_input(topology, overwrite=(topology%coord_type /= do_coord_off), &
657 10668 : subsys_section=subsys_section, save_mem=save_mem)
658 : END IF
659 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", &
660 10680 : i_val=topology%natoms)
661 10680 : CALL timestop(handle2)
662 : ! Check on atom numbers
663 10680 : IF (topology%natoms <= 0) &
664 0 : CPABORT("No atomic coordinates have been found! ")
665 10680 : CALL timestop(handle)
666 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
667 10680 : "PRINT%TOPOLOGY_INFO")
668 10680 : END SUBROUTINE coordinate_control
669 :
670 : ! **************************************************************************************************
671 : !> \brief Transform Cartesian coordinates read from external files to the internal cell frame.
672 : !> \param topology ...
673 : ! **************************************************************************************************
674 1893 : SUBROUTINE transform_external_cartesian_coordinates(topology)
675 :
676 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
677 :
678 : INTEGER :: iatom
679 :
680 1893 : IF (.NOT. ASSOCIATED(topology%cell_muc)) RETURN
681 :
682 505964 : DO iatom = 1, topology%natoms
683 505964 : CALL cell_transform_input_cartesian(topology%cell_muc, topology%atom_info%r(:, iatom))
684 : END DO
685 :
686 : END SUBROUTINE transform_external_cartesian_coordinates
687 :
688 : ! **************************************************************************************************
689 : !> \brief ...
690 : !> \param colvar_p ...
691 : !> \param particle_set ...
692 : !> \par History
693 : !> none
694 : !> \author Teodoro Laino [tlaino] - 07.2007
695 : ! **************************************************************************************************
696 10726 : SUBROUTINE topology_post_proc_colvar(colvar_p, particle_set)
697 :
698 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
699 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
700 :
701 : INTEGER :: i, j
702 :
703 10726 : IF (ASSOCIATED(colvar_p)) THEN
704 11189 : DO i = 1, SIZE(colvar_p)
705 11189 : IF (colvar_p(i)%colvar%type_id == combine_colvar_id) THEN
706 36 : DO j = 1, SIZE(colvar_p(i)%colvar%combine_cvs_param%colvar_p)
707 36 : CALL post_process_colvar(colvar_p(i)%colvar%combine_cvs_param%colvar_p(j)%colvar, particle_set)
708 : END DO
709 12 : CALL colvar_setup(colvar_p(i)%colvar)
710 : ELSE
711 454 : CALL post_process_colvar(colvar_p(i)%colvar, particle_set)
712 : END IF
713 : END DO
714 : END IF
715 10726 : END SUBROUTINE topology_post_proc_colvar
716 :
717 : ! **************************************************************************************************
718 : !> \brief Setup the cell used for handling properly the multiple_unit_cell option
719 : !> \param cell_muc ...
720 : !> \param cell ...
721 : !> \param subsys_section ...
722 : !> \author Teodoro Laino [tlaino] - 06.2009
723 : ! **************************************************************************************************
724 10726 : SUBROUTINE setup_cell_muc(cell_muc, cell, subsys_section)
725 :
726 : TYPE(cell_type), POINTER :: cell_muc, cell
727 : TYPE(section_vals_type), POINTER :: subsys_section
728 :
729 10726 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
730 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_ref
731 :
732 0 : CPASSERT(.NOT. ASSOCIATED(cell_muc))
733 :
734 : CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
735 10726 : i_vals=multiple_unit_cell)
736 42476 : IF (ANY(multiple_unit_cell /= 1)) THEN
737 : ! Restore the original cell
738 584 : hmat_ref(:, 1) = cell%hmat(:, 1)/multiple_unit_cell(1)
739 584 : hmat_ref(:, 2) = cell%hmat(:, 2)/multiple_unit_cell(2)
740 584 : hmat_ref(:, 3) = cell%hmat(:, 3)/multiple_unit_cell(3)
741 : ! Create the MUC cell
742 146 : CALL cell_create(cell_muc, hmat=hmat_ref, periodic=cell%perd, tag="CELL_UC")
743 146 : CALL write_cell(cell_muc, subsys_section)
744 : ELSE
745 : ! If a multiple_unit_cell was not requested just point to the original cell
746 10580 : CALL cell_retain(cell)
747 10580 : cell_muc => cell
748 : END IF
749 :
750 10726 : END SUBROUTINE setup_cell_muc
751 :
752 : END MODULE topology
|
