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