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 29208 : 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 9736 : NULLIFY (logger)
155 19472 : logger => cp_get_default_logger()
156 9736 : CALL timeset(routineN, handle)
157 9736 : NULLIFY (cell_section, constraint_section, topology_section)
158 :
159 9736 : cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
160 9736 : IF (use_motion_section) THEN
161 9417 : constraint_section => section_vals_get_subs_vals(root_section, "MOTION%CONSTRAINT")
162 : END IF
163 9736 : topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
164 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
165 9736 : extension=".mmLog")
166 9736 : my_qmmm = .FALSE.
167 9736 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
168 :
169 9736 : IF (PRESENT(elkind)) THEN
170 6900 : CALL section_vals_val_get(topology_section, "USE_ELEMENT_AS_KIND", explicit=explicit)
171 6900 : IF (explicit) THEN
172 0 : CALL section_vals_val_get(topology_section, "USE_ELEMENT_AS_KIND", l_val=el_as_kind)
173 : ELSE
174 6900 : el_as_kind = elkind
175 : END IF
176 : ELSE
177 2836 : 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 9736 : 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 9736 : para_env=para_env)
186 9736 : CALL write_cell(topology%cell, subsys_section, tag="CELL_TOP")
187 9736 : 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 9736 : CALL read_topology_section(topology, topology_section)
191 :
192 : ! 4. Read in the constraints section
193 9736 : 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 9736 : binary_coord_read)
198 9736 : IF (.NOT. binary_coord_read) THEN
199 9690 : 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 9736 : force_env_section)
205 :
206 9736 : IF (el_as_kind) THEN
207 : ! redefine atom names with the name of the element
208 24778 : 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 9736 : 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 9736 : 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 9736 : 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 9342 : force_env_section=force_env_section, exclusions=exclusions)
233 : END IF
234 :
235 : !10. Post-Process colvar definitions (if needed)
236 9736 : CALL topology_post_proc_colvar(colvar_p, particle_set)
237 :
238 : !11. Deal with the constraint stuff if requested
239 9736 : 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 9342 : 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 9736 : file_status="REPLACE", extension=".pdb")
252 9736 : 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 9736 : "TOPOLOGY%DUMP_PDB")
257 : iw2 = cp_print_key_unit_nr(logger, subsys_section, "TOPOLOGY%DUMP_PSF", &
258 9736 : file_status="REPLACE", extension=".psf")
259 9736 : 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 9736 : "TOPOLOGY%DUMP_PSF")
264 :
265 : !13. Cleanup the topology structure type
266 9736 : CALL deallocate_topology(topology)
267 9736 : CALL timestop(handle)
268 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
269 9736 : "PRINT%TOPOLOGY_INFO")
270 9736 : 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 10274 : 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 10274 : 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 10274 : NULLIFY (logger)
305 10274 : logger => cp_get_default_logger()
306 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
307 10274 : extension=".mmLog")
308 10274 : CALL timeset(routineN, handle)
309 :
310 10274 : my_qmmm = .FALSE.
311 10274 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
312 :
313 : ! 1. Read in the connectivity information (if this is the case)
314 11071 : 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 10274 : 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 10274 : IF (topology%aa_element) THEN
335 9046 : check = SIZE(topology%atom_info%id_element) == SIZE(topology%atom_info%id_atmname)
336 9046 : CPASSERT(check)
337 568364 : topology%atom_info%id_element = topology%atom_info%id_atmname
338 : END IF
339 :
340 : ! 3. Check for the element name..
341 10274 : 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 11071 : 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 10274 : use_mm_map_first = .FALSE.
350 : END SELECT
351 10274 : CALL create_ff_map("AMBER")
352 10274 : CALL create_ff_map("CHARMM")
353 10274 : CALL create_ff_map("GROMOS")
354 :
355 30822 : ALLOCATE (elements(SIZE(topology%atom_info%id_element)))
356 733701 : DO i = 1, SIZE(elements)
357 733701 : elements(i) = id2str(topology%atom_info%id_element(i))
358 : END DO
359 :
360 733701 : DO i = 1, topology%natoms
361 723427 : IF (elements(i) == "__DEF__") CYCLE
362 : ! If present an underscore let's skip all that over the underscore
363 23755 : id = INDEX(elements(i), "_") - 1
364 23755 : 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 28211 : DO j = id, 1, -1
370 28199 : itmp = ICHAR(elements(i) (j:j))
371 28211 : IF ((itmp < map0) .OR. (itmp > map9)) EXIT
372 : END DO
373 23755 : 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 23755 : subsys_section, use_mm_map_first)
377 : ! Earn time fixing same element labels for same atoms
378 23755 : element0 = elements(i)
379 25124492 : DO k = i, topology%natoms
380 25813890 : IF (element0 == id2str(topology%atom_info%id_element(k))) THEN
381 733109 : topology%atom_info%id_element(k) = str2id(s2s(my_element))
382 25823572 : elements(k) = "__DEF__"
383 : END IF
384 : END DO
385 : END DO
386 10274 : DEALLOCATE (elements)
387 10274 : CALL destroy_ff_map("GROMOS")
388 10274 : CALL destroy_ff_map("CHARMM")
389 10274 : CALL destroy_ff_map("AMBER")
390 10274 : CALL timestop(handle2)
391 :
392 : ! 4. Generate the connectivity information otherwise
393 18103 : SELECT CASE (topology%conn_type)
394 : CASE (do_conn_generate)
395 7829 : CALL topology_set_atm_mass(topology, subsys_section)
396 7829 : CALL topology_generate_bond(topology, para_env, subsys_section)
397 7829 : 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 7829 : CALL topology_generate_bend(topology, subsys_section)
407 7829 : CALL topology_generate_ub(topology, subsys_section)
408 7829 : CALL topology_generate_dihe(topology, subsys_section)
409 7829 : CALL topology_generate_impr(topology, subsys_section)
410 7829 : 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 11922 : CALL topology_generate_onfo(topology, subsys_section)
418 : END SELECT
419 :
420 : ! 5. Handle multiple unit_cell - Update atoms_info
421 10274 : CALL topology_muc(topology, subsys_section)
422 :
423 : ! 6. Handle multiple unit_cell - Update conn_info
424 10274 : CALL topology_conn_multiple(topology, subsys_section)
425 :
426 : ! 7. Generate Molecules
427 10274 : CALL topology_generate_molecule(topology, my_qmmm, qmmm_env, subsys_section)
428 10274 : IF (topology%molecules_check) CALL topology_molecules_check(topology, subsys_section)
429 :
430 : ! 8. Modify for QM/MM
431 10274 : IF (my_qmmm) THEN
432 394 : CALL qmmm_coordinate_control(topology, qmmm_env, subsys_section)
433 : END IF
434 10274 : CALL timestop(handle)
435 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
436 10274 : "PRINT%TOPOLOGY_INFO")
437 :
438 20548 : 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 38760 : 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 9690 : NULLIFY (logger)
533 9690 : logger => cp_get_default_logger()
534 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO", &
535 9690 : extension=".mmLog")
536 9690 : CALL timeset(routineN, handle)
537 :
538 9690 : NULLIFY (global_section)
539 9690 : global_section => section_vals_get_subs_vals(root_section, "GLOBAL")
540 9690 : CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem)
541 : !-----------------------------------------------------------------------------
542 : ! 1. If reading in from external file, make sure its there first
543 : !-----------------------------------------------------------------------------
544 9690 : IF (topology%coordinate) THEN
545 1883 : INQUIRE (FILE=topology%coord_file_name, EXIST=found, IOSTAT=istat)
546 1883 : IF (istat /= 0) THEN
547 : WRITE (UNIT=message, FMT="(A,I0,A)") &
548 : "An error occurred inquiring the file <"// &
549 0 : TRIM(topology%coord_file_name)//"> (IOSTAT = ", istat, ")"
550 0 : CPABORT(TRIM(message))
551 : END IF
552 1883 : IF (.NOT. found) THEN
553 : CALL cp_abort(__LOCATION__, &
554 : "Coordinate file <"//TRIM(topology%coord_file_name)// &
555 0 : "> not found.")
556 : END IF
557 : END IF
558 : !-----------------------------------------------------------------------------
559 : ! 2. Read in the coordinates from the corresponding locations
560 : !-----------------------------------------------------------------------------
561 9690 : CALL timeset(routineN//"_READ_COORDINATE", handle2)
562 9704 : SELECT CASE (topology%coord_type)
563 : CASE (do_coord_off)
564 : ! Do nothing.. we will parse later from the &COORD section..
565 : CASE (do_coord_g96)
566 14 : CALL read_coordinate_g96(topology, para_env, subsys_section)
567 : CASE (do_coord_crd)
568 26 : CALL read_coordinate_crd(topology, para_env, subsys_section)
569 : CASE (do_coord_pdb)
570 804 : CALL read_coordinate_pdb(topology, para_env, subsys_section)
571 : CASE (do_coord_xyz)
572 1011 : CALL read_coordinate_xyz(topology, para_env, subsys_section)
573 : CASE (do_coord_cif)
574 14 : CALL read_coordinate_cif(topology, para_env, subsys_section)
575 : CASE (do_coord_xtl)
576 2 : CALL read_coordinate_xtl(topology, para_env, subsys_section)
577 : CASE (do_coord_cp2k)
578 12 : CALL read_coordinate_cp2k(topology, para_env, subsys_section)
579 : CASE DEFAULT
580 : ! We should never reach this point..
581 9690 : CPABORT("")
582 : END SELECT
583 :
584 : ! Parse &COORD section and in case overwrite
585 9690 : IF (topology%coord_type /= do_coord_cp2k) THEN
586 : CALL read_atoms_input(topology, overwrite=(topology%coord_type /= do_coord_off), &
587 9678 : subsys_section=subsys_section, save_mem=save_mem)
588 : END IF
589 : CALL section_vals_val_set(subsys_section, "TOPOLOGY%NUMBER_OF_ATOMS", &
590 9690 : i_val=topology%natoms)
591 9690 : CALL timestop(handle2)
592 : ! Check on atom numbers
593 9690 : IF (topology%natoms <= 0) &
594 0 : CPABORT("No atomic coordinates have been found! ")
595 9690 : CALL timestop(handle)
596 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
597 9690 : "PRINT%TOPOLOGY_INFO")
598 9690 : END SUBROUTINE coordinate_control
599 :
600 : ! **************************************************************************************************
601 : !> \brief ...
602 : !> \param colvar_p ...
603 : !> \param particle_set ...
604 : !> \par History
605 : !> none
606 : !> \author Teodoro Laino [tlaino] - 07.2007
607 : ! **************************************************************************************************
608 9736 : SUBROUTINE topology_post_proc_colvar(colvar_p, particle_set)
609 :
610 : TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p
611 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
612 :
613 : INTEGER :: i, j
614 :
615 9736 : IF (ASSOCIATED(colvar_p)) THEN
616 10197 : DO i = 1, SIZE(colvar_p)
617 10197 : IF (colvar_p(i)%colvar%type_id == combine_colvar_id) THEN
618 30 : DO j = 1, SIZE(colvar_p(i)%colvar%combine_cvs_param%colvar_p)
619 30 : CALL post_process_colvar(colvar_p(i)%colvar%combine_cvs_param%colvar_p(j)%colvar, particle_set)
620 : END DO
621 10 : CALL colvar_setup(colvar_p(i)%colvar)
622 : ELSE
623 454 : CALL post_process_colvar(colvar_p(i)%colvar, particle_set)
624 : END IF
625 : END DO
626 : END IF
627 9736 : END SUBROUTINE topology_post_proc_colvar
628 :
629 : ! **************************************************************************************************
630 : !> \brief Setup the cell used for handling properly the multiple_unit_cell option
631 : !> \param cell_muc ...
632 : !> \param cell ...
633 : !> \param subsys_section ...
634 : !> \author Teodoro Laino [tlaino] - 06.2009
635 : ! **************************************************************************************************
636 9736 : SUBROUTINE setup_cell_muc(cell_muc, cell, subsys_section)
637 :
638 : TYPE(cell_type), POINTER :: cell_muc, cell
639 : TYPE(section_vals_type), POINTER :: subsys_section
640 :
641 9736 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
642 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_ref
643 :
644 0 : CPASSERT(.NOT. ASSOCIATED(cell_muc))
645 :
646 : CALL section_vals_val_get(subsys_section, "CELL%MULTIPLE_UNIT_CELL", &
647 9736 : i_vals=multiple_unit_cell)
648 38510 : IF (ANY(multiple_unit_cell /= 1)) THEN
649 : ! Restore the original cell
650 592 : hmat_ref(:, 1) = cell%hmat(:, 1)/multiple_unit_cell(1)
651 592 : hmat_ref(:, 2) = cell%hmat(:, 2)/multiple_unit_cell(2)
652 592 : hmat_ref(:, 3) = cell%hmat(:, 3)/multiple_unit_cell(3)
653 : ! Create the MUC cell
654 148 : CALL cell_create(cell_muc, hmat=hmat_ref, periodic=cell%perd, tag="CELL_UC")
655 148 : CALL write_cell(cell_muc, subsys_section)
656 : ELSE
657 : ! If a multiple_unit_cell was not requested just point to the original cell
658 9588 : CALL cell_retain(cell)
659 9588 : cell_muc => cell
660 : END IF
661 :
662 9736 : END SUBROUTINE setup_cell_muc
663 :
664 : END MODULE topology
|
