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