(git:0446f17)
Loading...
Searching...
No Matches
topology_util.F
Go to the documentation of this file.
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 Collection of subroutine needed for topology related things
10!> \par History
11!> jgh (23-05-2004) Last atom of molecule information added
12! **************************************************************************************************
13MODULE topology_util
14 USE cp_log_handling, ONLY: cp_get_default_logger,&
15 cp_logger_get_default_io_unit,&
16 cp_logger_type,&
17 cp_to_string
18 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 cp_print_key_unit_nr
20 USE graphcon, ONLY: graph_type,&
21 hash_molecule,&
22 reorder_graph,&
23 vertex
24 USE input_section_types, ONLY: section_vals_get,&
25 section_vals_get_subs_vals,&
26 section_vals_type,&
27 section_vals_val_get,&
28 section_vals_val_set
29 USE kinds, ONLY: default_string_length,&
30 dp
31 USE mm_mapping_library, ONLY: amber_map,&
32 charmm_map,&
33 gromos_map
34 USE periodic_table, ONLY: get_ptable_info
35 USE qmmm_types_low, ONLY: qmmm_env_mm_type
36 USE string_table, ONLY: id2str,&
37 str2id
38 USE string_utilities, ONLY: uppercase
39 USE topology_types, ONLY: atom_info_type,&
40 connectivity_info_type,&
41 topology_parameters_type
42 USE util, ONLY: sort
43#include "./base/base_uses.f90"
44
45 IMPLICIT NONE
46
47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_util'
48
49 TYPE array1_list_type
50 INTEGER, POINTER, DIMENSION(:) :: array1 => null()
51 END TYPE array1_list_type
52
53 TYPE array2_list_type
54 INTEGER, POINTER, DIMENSION(:) :: array1 => null()
55 INTEGER, POINTER, DIMENSION(:) :: array2 => null()
56 END TYPE array2_list_type
57
58 PRIVATE
59 PUBLIC :: topology_set_atm_mass, &
60 topology_reorder_atoms, &
61 topology_molecules_check, &
62 check_subsys_element, &
63 reorder_structure, &
64 find_molecule, &
65 array1_list_type, &
66 array2_list_type, &
67 give_back_molecule, &
68 reorder_list_array, &
69 tag_molecule
70
71 INTERFACE reorder_structure
72 MODULE PROCEDURE reorder_structure1d, reorder_structure2d
73 END INTERFACE
74
75CONTAINS
76
77! **************************************************************************************************
78!> \brief ...
79!> \param topology ...
80!> \param qmmm ...
81!> \param qmmm_env_mm ...
82!> \param subsys_section ...
83!> \param force_env_section ...
84!> \par History
85!> Teodoro Laino 09.2006 - Rewritten with a graph matching algorithm
86! **************************************************************************************************
87 SUBROUTINE topology_reorder_atoms(topology, qmmm, qmmm_env_mm, subsys_section, force_env_section)
88 TYPE(topology_parameters_type), INTENT(INOUT) :: topology
89 LOGICAL, INTENT(in), OPTIONAL :: qmmm
90 TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env_mm
91 TYPE(section_vals_type), POINTER :: subsys_section, force_env_section
92
93 CHARACTER(len=*), PARAMETER :: routinen = 'topology_reorder_atoms'
94
95 CHARACTER(LEN=default_string_length) :: mol_id
96 CHARACTER(LEN=default_string_length), POINTER :: molname(:), telement(:), &
97 tlabel_atmname(:), tlabel_molname(:), &
98 tlabel_resname(:)
99 INTEGER :: handle, i, iatm, iindex, ikind, imol, imol_ref, iref, iund, iw, j, k, location, &
100 max_mol_num, mm_index, n, n_rep, n_var, natom, natom_loc, nkind, nlinks, old_hash, &
101 old_mol, output_unit, qm_index, unique_mol
102 INTEGER, DIMENSION(:), POINTER :: mm_link_atoms, qm_atom_index
103 INTEGER, POINTER :: atm_map1(:), atm_map2(:), atm_map3(:), map_atom_type(:), &
104 map_mol_hash(:), mm_indexes_n(:), mm_indexes_v(:), mol_bnd(:, :), mol_hash(:), &
105 mol_num(:), new_position(:), order(:), tmp_n(:), tmp_v(:), wrk(:)
106 LOGICAL :: explicit, matches, my_qmmm
107 REAL(kind=dp), DIMENSION(:, :), POINTER :: tr
108 REAL(kind=dp), POINTER :: tatm_charge(:), tatm_mass(:)
109 TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
110 TYPE(atom_info_type), POINTER :: atom_info
111 TYPE(connectivity_info_type), POINTER :: conn_info
112 TYPE(cp_logger_type), POINTER :: logger
113 TYPE(graph_type), DIMENSION(:), POINTER :: reference_set
114 TYPE(section_vals_type), POINTER :: generate_section, isolated_section, &
115 qm_kinds, qmmm_link_section, &
116 qmmm_section
117 TYPE(vertex), DIMENSION(:), POINTER :: reference, unordered
118
119 NULLIFY (logger, generate_section, isolated_section, tmp_v, tmp_n)
120 logger => cp_get_default_logger()
121 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
122 extension=".subsysLog")
123 CALL timeset(routinen, handle)
124 output_unit = cp_logger_get_default_io_unit(logger)
125 IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
126
127 my_qmmm = .false.
128 IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env_mm)) my_qmmm = qmmm
129
130 atom_info => topology%atom_info
131 conn_info => topology%conn_info
132 natom = topology%natoms
133
134 NULLIFY (new_position, reference_set)
135 NULLIFY (tlabel_atmname, telement, mol_num, tlabel_molname, tlabel_resname)
136 NULLIFY (tr, tatm_charge, tatm_mass, atm_map1, atm_map2, atm_map3)
137 ! This routine can be called only at a very high level where these structures are still
138 ! not even taken into account...
139 cpassert(.NOT. ASSOCIATED(atom_info%map_mol_num))
140 cpassert(.NOT. ASSOCIATED(atom_info%map_mol_typ))
141 cpassert(.NOT. ASSOCIATED(atom_info%map_mol_res))
142 !ALLOCATE all the temporary arrays needed for this routine
143 ALLOCATE (new_position(natom))
144 ALLOCATE (mol_num(natom))
145 ALLOCATE (molname(natom))
146 ALLOCATE (tlabel_atmname(natom))
147 ALLOCATE (tlabel_molname(natom))
148 ALLOCATE (tlabel_resname(natom))
149 ALLOCATE (tr(3, natom))
150 ALLOCATE (tatm_charge(natom))
151 ALLOCATE (tatm_mass(natom))
152 ALLOCATE (telement(natom))
153 ALLOCATE (atm_map1(natom))
154 ALLOCATE (atm_map2(natom))
155
156 ! The only information we have at this level is the connectivity of the system.
157 ! 0. Build a list of mapping atom types
158 ALLOCATE (map_atom_type(natom))
159 ! 1. Build a list of bonds
160 ALLOCATE (atom_bond_list(natom))
161 DO i = 1, natom
162 map_atom_type(i) = atom_info%id_atmname(i)
163 ALLOCATE (atom_bond_list(i)%array1(0))
164 END DO
165 n = 0
166 IF (ASSOCIATED(conn_info%bond_a)) n = SIZE(conn_info%bond_a)
167 CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, n)
168 ALLOCATE (atom_info%map_mol_num(natom))
169 atom_info%map_mol_num = -1
170 CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
171 max_mol_num = maxval(atom_info%map_mol_num)
172 ! In atom_info%map_mol_num have already been mapped molecules
173 ALLOCATE (mol_bnd(2, max_mol_num))
174 ALLOCATE (mol_hash(max_mol_num))
175 ALLOCATE (map_mol_hash(max_mol_num))
176 ! 2. Sort the map_mol_num array.. atm_map1 contains now the mapped index
177 ! of the reordered array
178 CALL sort(atom_info%map_mol_num, natom, atm_map1)
179 old_mol = 0
180 iindex = 0
181 imol = 0
182 DO i = 1, natom
183 IF (old_mol .NE. atom_info%map_mol_num(i)) THEN
184 old_mol = atom_info%map_mol_num(i)
185 iindex = 0
186 IF (imol > 0) THEN
187 mol_bnd(2, imol) = i - 1
188 END IF
189 imol = imol + 1
190 mol_bnd(1, imol) = i
191 END IF
192 iindex = iindex + 1
193 atm_map2(atm_map1(i)) = iindex
194 END DO
195 mol_bnd(2, imol) = natom
196 ! Indexes of the two molecules to check
197 iref = 1
198 iund = max_mol_num/2 + 1
199 ! Allocate reference and unordered
200 NULLIFY (reference, unordered)
201 ! This is the real matching of graphs
202 DO j = 1, max_mol_num
203 CALL setup_graph(j, reference, map_atom_type, &
204 atom_bond_list, mol_bnd, atm_map1, atm_map2)
205
206 ALLOCATE (order(SIZE(reference)))
207 CALL hash_molecule(reference, order, mol_hash(j))
208
209 DEALLOCATE (order)
210 DO i = 1, SIZE(reference)
211 DEALLOCATE (reference(i)%bonds)
212 END DO
213 DEALLOCATE (reference)
214 END DO
215 ! Reorder molecules hashes
216 CALL sort(mol_hash, max_mol_num, map_mol_hash)
217 ! Now find unique molecules and reorder atoms too (if molecules match)
218 old_hash = -1
219 unique_mol = 0
220 natom_loc = 0
221 IF (output_unit > 0) THEN
222 WRITE (output_unit, '(T2,"REORDER | ",A)') &
223 "Reordering Molecules. The Reordering of molecules will override all", &
224 "information regarding molecule names and residue names.", &
225 "New ones will be provided on the basis of the connectivity!"
226 END IF
227 DO j = 1, max_mol_num
228 IF (mol_hash(j) .NE. old_hash) THEN
229 unique_mol = unique_mol + 1
230 old_hash = mol_hash(j)
231 CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
232 map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
233 atm_map3)
234 ! Reorder Last added reference
235 mol_id = trim(adjustl(cp_to_string(unique_mol)))
236 DO i = 1, SIZE(atm_map3)
237 natom_loc = natom_loc + 1
238 new_position(natom_loc) = atm_map3(i)
239 molname(natom_loc) = mol_id
240 mol_num(natom_loc) = unique_mol
241 END DO
242 DEALLOCATE (atm_map3)
243 ELSE
244 CALL setup_graph(map_mol_hash(j), unordered, map_atom_type, &
245 atom_bond_list, mol_bnd, atm_map1, atm_map2, atm_map3)
246 DO imol_ref = 1, unique_mol
247 !
248 reference => reference_set(imol_ref)%graph
249 ALLOCATE (order(SIZE(reference)))
250 CALL reorder_graph(reference, unordered, order, matches)
251 IF (matches) EXIT
252 DEALLOCATE (order)
253 END DO
254 IF (matches) THEN
255 ! Reorder according the correct index
256 ALLOCATE (wrk(SIZE(order)))
257 CALL sort(order, SIZE(order), wrk)
258 DO i = 1, SIZE(order)
259 natom_loc = natom_loc + 1
260 new_position(natom_loc) = atm_map3(wrk(i))
261 molname(natom_loc) = mol_id
262 mol_num(natom_loc) = unique_mol
263 END DO
264 !
265 DEALLOCATE (order)
266 DEALLOCATE (wrk)
267 ELSE
268 unique_mol = unique_mol + 1
269 CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
270 map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
271 atm_map3)
272 ! Reorder Last added reference
273 mol_id = trim(adjustl(cp_to_string(unique_mol)))
274 DO i = 1, SIZE(atm_map3)
275 natom_loc = natom_loc + 1
276 new_position(natom_loc) = atm_map3(i)
277 molname(natom_loc) = mol_id
278 mol_num(natom_loc) = unique_mol
279 END DO
280 DEALLOCATE (atm_map3)
281 END IF
282 DO i = 1, SIZE(unordered)
283 DEALLOCATE (unordered(i)%bonds)
284 END DO
285 DEALLOCATE (unordered)
286 DEALLOCATE (atm_map3)
287 END IF
288 END DO
289 IF (output_unit > 0) THEN
290 WRITE (output_unit, '(T2,"REORDER | ",A,I7,A)') "Number of unique molecules found:", unique_mol, "."
291 END IF
292 cpassert(natom_loc == natom)
293 DEALLOCATE (map_atom_type)
294 DEALLOCATE (atm_map1)
295 DEALLOCATE (atm_map2)
296 DEALLOCATE (mol_bnd)
297 DEALLOCATE (mol_hash)
298 DEALLOCATE (map_mol_hash)
299 ! Deallocate working arrays
300 DO i = 1, natom
301 DEALLOCATE (atom_bond_list(i)%array1)
302 END DO
303 DEALLOCATE (atom_bond_list)
304 DEALLOCATE (atom_info%map_mol_num)
305 ! Deallocate reference_set
306 DO i = 1, SIZE(reference_set)
307 DO j = 1, SIZE(reference_set(i)%graph)
308 DEALLOCATE (reference_set(i)%graph(j)%bonds)
309 END DO
310 DEALLOCATE (reference_set(i)%graph)
311 END DO
312 DEALLOCATE (reference_set)
313 !Lets swap the atoms now
314 DO iatm = 1, natom
315 location = new_position(iatm)
316 tlabel_molname(iatm) = id2str(atom_info%id_molname(location))
317 tlabel_resname(iatm) = id2str(atom_info%id_resname(location))
318 tlabel_atmname(iatm) = id2str(atom_info%id_atmname(location))
319 telement(iatm) = id2str(atom_info%id_element(location))
320 tr(1, iatm) = atom_info%r(1, location)
321 tr(2, iatm) = atom_info%r(2, location)
322 tr(3, iatm) = atom_info%r(3, location)
323 tatm_charge(iatm) = atom_info%atm_charge(location)
324 tatm_mass(iatm) = atom_info%atm_mass(location)
325 END DO
326 IF (topology%create_molecules) THEN
327 DO iatm = 1, natom
328 tlabel_molname(iatm) = "MOL"//trim(molname(iatm))
329 tlabel_resname(iatm) = "R"//trim(molname(iatm))
330 END DO
331 topology%create_molecules = .false.
332 END IF
333 DO iatm = 1, natom
334 atom_info%id_molname(iatm) = str2id(tlabel_molname(iatm))
335 atom_info%id_resname(iatm) = str2id(tlabel_resname(iatm))
336 atom_info%id_atmname(iatm) = str2id(tlabel_atmname(iatm))
337 atom_info%id_element(iatm) = str2id(telement(iatm))
338 atom_info%resid(iatm) = mol_num(iatm)
339 atom_info%r(1, iatm) = tr(1, iatm)
340 atom_info%r(2, iatm) = tr(2, iatm)
341 atom_info%r(3, iatm) = tr(3, iatm)
342 atom_info%atm_charge(iatm) = tatm_charge(iatm)
343 atom_info%atm_mass(iatm) = tatm_mass(iatm)
344 END DO
345
346 ! Let's reorder all the list provided in the input..
347 ALLOCATE (wrk(SIZE(new_position)))
348 CALL sort(new_position, SIZE(new_position), wrk)
349
350 ! NOTE: In the future the whole list of possible integers should be updated at this level..
351 ! Let's update all the integer lists in the qmmm_env_mm section and in the input sections
352 ! from where qmmm_env_qm will be read.
353 IF (my_qmmm) THEN
354 ! Update the qmmm_env_mm
355 cpassert(SIZE(qmmm_env_mm%qm_atom_index) /= 0)
356 cpassert(all(qmmm_env_mm%qm_atom_index /= 0))
357 ALLOCATE (qm_atom_index(SIZE(qmmm_env_mm%qm_atom_index)))
358 DO iatm = 1, SIZE(qmmm_env_mm%qm_atom_index)
359 qm_atom_index(iatm) = wrk(qmmm_env_mm%qm_atom_index(iatm))
360 END DO
361 qmmm_env_mm%qm_atom_index = qm_atom_index
362 cpassert(all(qmmm_env_mm%qm_atom_index /= 0))
363 DEALLOCATE (qm_atom_index)
364 ! Update the qmmm_section: MM_INDEX of the QM_KIND
365 qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
366 qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
367 CALL section_vals_get(qm_kinds, n_repetition=nkind)
368 DO ikind = 1, nkind
369 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
370 DO k = 1, n_var
371 CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
372 i_vals=mm_indexes_v)
373 ALLOCATE (mm_indexes_n(SIZE(mm_indexes_v)))
374 DO j = 1, SIZE(mm_indexes_v)
375 mm_indexes_n(j) = wrk(mm_indexes_v(j))
376 END DO
377 CALL section_vals_val_set(qm_kinds, "MM_INDEX", i_rep_section=ikind, &
378 i_vals_ptr=mm_indexes_n, i_rep_val=k)
379 END DO
380 END DO
381 ! Handle the link atoms
382 IF (qmmm_env_mm%qmmm_link) THEN
383 ! Update the qmmm_env_mm
384 cpassert(SIZE(qmmm_env_mm%mm_link_atoms) > 0)
385 ALLOCATE (mm_link_atoms(SIZE(qmmm_env_mm%mm_link_atoms)))
386 DO iatm = 1, SIZE(qmmm_env_mm%mm_link_atoms)
387 mm_link_atoms(iatm) = wrk(qmmm_env_mm%mm_link_atoms(iatm))
388 END DO
389 qmmm_env_mm%mm_link_atoms = mm_link_atoms
390 cpassert(all(qmmm_env_mm%mm_link_atoms /= 0))
391 DEALLOCATE (mm_link_atoms)
392 ! Update the qmmm_section: MM_INDEX,QM_INDEX of the LINK atom list
393 qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
394 CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
395 cpassert(nlinks /= 0)
396 DO ikind = 1, nlinks
397 CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
398 CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
399 mm_index = wrk(mm_index)
400 qm_index = wrk(qm_index)
401 CALL section_vals_val_set(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
402 CALL section_vals_val_set(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
403 END DO
404 END IF
405 END IF
406 !
407 !LIST of ISOLATED atoms
408 !
409 generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
410 isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
411 CALL section_vals_get(isolated_section, explicit=explicit)
412 IF (explicit) THEN
413 CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
414 DO i = 1, n_rep
415 CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
416 ALLOCATE (tmp_n(SIZE(tmp_v)))
417 DO j = 1, SIZE(tmp_v)
418 tmp_n(j) = wrk(tmp_v(j))
419 END DO
420 CALL section_vals_val_set(isolated_section, "LIST", i_vals_ptr=tmp_n, i_rep_val=i)
421 END DO
422 END IF
423 DEALLOCATE (wrk)
424 !DEALLOCATE all the temporary arrays needed for this routine
425 DEALLOCATE (new_position)
426 DEALLOCATE (tlabel_atmname)
427 DEALLOCATE (tlabel_molname)
428 DEALLOCATE (tlabel_resname)
429 DEALLOCATE (telement)
430 DEALLOCATE (tr)
431 DEALLOCATE (tatm_charge)
432 DEALLOCATE (tatm_mass)
433 DEALLOCATE (molname)
434 DEALLOCATE (mol_num)
435
436 ! DEALLOCATE the bond structures in the connectivity
437 DEALLOCATE (conn_info%bond_a)
438 DEALLOCATE (conn_info%bond_b)
439 IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
440 CALL timestop(handle)
441 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
442 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
443 END SUBROUTINE topology_reorder_atoms
444
445! **************************************************************************************************
446!> \brief Set up a SET of graph kind
447!> \param graph_set ...
448!> \param idim ...
449!> \param ind ...
450!> \param array2 ...
451!> \param atom_bond_list ...
452!> \param map_mol ...
453!> \param atm_map1 ...
454!> \param atm_map2 ...
455!> \param atm_map3 ...
456!> \author Teodoro Laino 10.2006
457! **************************************************************************************************
458 SUBROUTINE setup_graph_set(graph_set, idim, ind, array2, atom_bond_list, map_mol, &
459 atm_map1, atm_map2, atm_map3)
460 TYPE(graph_type), DIMENSION(:), POINTER :: graph_set
461 INTEGER, INTENT(IN) :: idim, ind
462 INTEGER, DIMENSION(:), POINTER :: array2
463 TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
464 INTEGER, DIMENSION(:, :), POINTER :: map_mol
465 INTEGER, DIMENSION(:), POINTER :: atm_map1, atm_map2, atm_map3
466
467 INTEGER :: ldim
468 TYPE(graph_type), DIMENSION(:), POINTER :: tmp_graph_set
469
470 ldim = 0
471 NULLIFY (tmp_graph_set)
472 IF (ASSOCIATED(graph_set)) THEN
473 ldim = SIZE(graph_set)
474 cpassert(ldim + 1 == idim)
475 mark_used(idim)
476 NULLIFY (tmp_graph_set)
477 CALL allocate_graph_set(graph_set, tmp_graph_set)
478 END IF
479 CALL allocate_graph_set(tmp_graph_set, graph_set, ldim, ldim + 1)
480 CALL setup_graph(ind, graph_set(ldim + 1)%graph, array2, &
481 atom_bond_list, map_mol, atm_map1, atm_map2, atm_map3)
482
483 END SUBROUTINE setup_graph_set
484
485! **************************************************************************************************
486!> \brief Allocate a new graph_set deallocating an old one..
487!> \param graph_set_in ...
488!> \param graph_set_out ...
489!> \param ldim_in ...
490!> \param ldim_out ...
491!> \author Teodoro Laino 10.2006
492! **************************************************************************************************
493 SUBROUTINE allocate_graph_set(graph_set_in, graph_set_out, ldim_in, ldim_out)
494 TYPE(graph_type), DIMENSION(:), POINTER :: graph_set_in, graph_set_out
495 INTEGER, INTENT(IN), OPTIONAL :: ldim_in, ldim_out
496
497 INTEGER :: b_dim, i, j, mydim_in, mydim_out, v_dim
498
499 cpassert(.NOT. ASSOCIATED(graph_set_out))
500 mydim_in = 0
501 mydim_out = 0
502 IF (ASSOCIATED(graph_set_in)) THEN
503 mydim_in = SIZE(graph_set_in)
504 mydim_out = SIZE(graph_set_in)
505 END IF
506 IF (PRESENT(ldim_in)) mydim_in = ldim_in
507 IF (PRESENT(ldim_out)) mydim_out = ldim_out
508 ALLOCATE (graph_set_out(mydim_out))
509 DO i = 1, mydim_out
510 NULLIFY (graph_set_out(i)%graph)
511 END DO
512 ! Copy graph structure into the temporary array
513 DO i = 1, mydim_in
514 v_dim = SIZE(graph_set_in(i)%graph)
515 ALLOCATE (graph_set_out(i)%graph(v_dim))
516 DO j = 1, v_dim
517 graph_set_out(i)%graph(j)%kind = graph_set_in(i)%graph(j)%kind
518 b_dim = SIZE(graph_set_in(i)%graph(j)%bonds)
519 ALLOCATE (graph_set_out(i)%graph(j)%bonds(b_dim))
520 graph_set_out(i)%graph(j)%bonds = graph_set_in(i)%graph(j)%bonds
521 DEALLOCATE (graph_set_in(i)%graph(j)%bonds)
522 END DO
523 DEALLOCATE (graph_set_in(i)%graph)
524 END DO
525 IF (ASSOCIATED(graph_set_in)) THEN
526 DEALLOCATE (graph_set_in)
527 END IF
528
529 END SUBROUTINE allocate_graph_set
530
531! **************************************************************************************************
532!> \brief Set up a graph kind
533!> \param ind ...
534!> \param graph ...
535!> \param array2 ...
536!> \param atom_bond_list ...
537!> \param map_mol ...
538!> \param atm_map1 ...
539!> \param atm_map2 ...
540!> \param atm_map3 ...
541!> \author Teodoro Laino 09.2006
542! **************************************************************************************************
543 SUBROUTINE setup_graph(ind, graph, array2, atom_bond_list, map_mol, &
544 atm_map1, atm_map2, atm_map3)
545 INTEGER, INTENT(IN) :: ind
546 TYPE(vertex), DIMENSION(:), POINTER :: graph
547 INTEGER, DIMENSION(:), POINTER :: array2
548 TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
549 INTEGER, DIMENSION(:, :), POINTER :: map_mol
550 INTEGER, DIMENSION(:), POINTER :: atm_map1, atm_map2
551 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: atm_map3
552
553 INTEGER :: i, idim, ifirst, ilast, j, nbonds, &
554 nelement
555
556 IF (PRESENT(atm_map3)) THEN
557 cpassert(.NOT. ASSOCIATED(atm_map3))
558 END IF
559 cpassert(.NOT. ASSOCIATED(graph))
560 ! Setup reference graph
561 idim = 0
562 ifirst = map_mol(1, ind)
563 ilast = map_mol(2, ind)
564 nelement = ilast - ifirst + 1
565 ALLOCATE (graph(nelement))
566 IF (PRESENT(atm_map3)) THEN
567 ALLOCATE (atm_map3(nelement))
568 END IF
569 DO i = ifirst, ilast
570 idim = idim + 1
571 graph(idim)%kind = array2(atm_map1(i))
572 nbonds = SIZE(atom_bond_list(atm_map1(i))%array1)
573 ALLOCATE (graph(idim)%bonds(nbonds))
574 DO j = 1, nbonds
575 graph(idim)%bonds(j) = atm_map2(atom_bond_list(atm_map1(i))%array1(j))
576 END DO
577 IF (PRESENT(atm_map3)) THEN
578 atm_map3(idim) = atm_map1(i)
579 END IF
580 END DO
581
582 END SUBROUTINE setup_graph
583
584! **************************************************************************************************
585!> \brief Order arrays of lists
586!> \param Ilist1 ...
587!> \param Ilist2 ...
588!> \param Ilist3 ...
589!> \param Ilist4 ...
590!> \param nsize ...
591!> \param ndim ...
592!> \author Teodoro Laino 09.2006
593! **************************************************************************************************
594 RECURSIVE SUBROUTINE reorder_list_array(Ilist1, Ilist2, Ilist3, Ilist4, nsize, ndim)
595 INTEGER, DIMENSION(:), POINTER :: ilist1
596 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ilist2, ilist3, ilist4
597 INTEGER, INTENT(IN) :: nsize, ndim
598
599 INTEGER :: i, iend, istart, ldim
600 INTEGER, DIMENSION(:), POINTER :: tmp, wrk
601
602 cpassert(nsize > 0)
603 ALLOCATE (wrk(ndim))
604 CALL sort(ilist1, ndim, wrk)
605 IF (nsize /= 1) THEN
606 ALLOCATE (tmp(ndim))
607 tmp = ilist2(1:ndim)
608 DO i = 1, ndim
609 ilist2(i) = tmp(wrk(i))
610 END DO
611 SELECT CASE (nsize)
612 CASE (3)
613 tmp = ilist3(1:ndim)
614 DO i = 1, ndim
615 ilist3(i) = tmp(wrk(i))
616 END DO
617 CASE (4)
618 tmp = ilist3(1:ndim)
619 DO i = 1, ndim
620 ilist3(i) = tmp(wrk(i))
621 END DO
622 tmp = ilist4(1:ndim)
623 DO i = 1, ndim
624 ilist4(i) = tmp(wrk(i))
625 END DO
626 END SELECT
627 DEALLOCATE (tmp)
628 istart = 1
629 DO i = 1, ndim
630 IF (ilist1(i) /= ilist1(istart)) THEN
631 iend = i - 1
632 ldim = iend - istart + 1
633 CALL reorder_list_array_low(ilist2, ilist3, ilist4, nsize, &
634 ldim, istart, iend)
635 istart = i
636 END IF
637 END DO
638 ! Last term to sort
639 iend = ndim
640 ldim = iend - istart + 1
641 CALL reorder_list_array_low(ilist2, ilist3, ilist4, nsize, &
642 ldim, istart, iend)
643 END IF
644 DEALLOCATE (wrk)
645 END SUBROUTINE reorder_list_array
646
647! **************************************************************************************************
648!> \brief Low level routine for ordering arrays of lists
649!> \param Ilist2 ...
650!> \param Ilist3 ...
651!> \param Ilist4 ...
652!> \param nsize ...
653!> \param ldim ...
654!> \param istart ...
655!> \param iend ...
656!> \author Teodoro Laino 09.2006
657! **************************************************************************************************
658 RECURSIVE SUBROUTINE reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
659 ldim, istart, iend)
660 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ilist2, ilist3, ilist4
661 INTEGER, INTENT(IN) :: nsize, ldim, istart, iend
662
663 INTEGER, DIMENSION(:), POINTER :: tmp_2, tmp_3, tmp_4
664
665 SELECT CASE (nsize)
666 CASE (2)
667 ALLOCATE (tmp_2(ldim))
668 tmp_2(:) = ilist2(istart:iend)
669 CALL reorder_list_array(tmp_2, nsize=nsize - 1, ndim=ldim)
670 ilist2(istart:iend) = tmp_2(:)
671 DEALLOCATE (tmp_2)
672 CASE (3)
673 ALLOCATE (tmp_2(ldim))
674 ALLOCATE (tmp_3(ldim))
675 tmp_2(:) = ilist2(istart:iend)
676 tmp_3(:) = ilist3(istart:iend)
677 CALL reorder_list_array(tmp_2, tmp_3, nsize=nsize - 1, ndim=ldim)
678 ilist2(istart:iend) = tmp_2(:)
679 ilist3(istart:iend) = tmp_3(:)
680 DEALLOCATE (tmp_2)
681 DEALLOCATE (tmp_3)
682 CASE (4)
683 ALLOCATE (tmp_2(ldim))
684 ALLOCATE (tmp_3(ldim))
685 ALLOCATE (tmp_4(ldim))
686 tmp_2(:) = ilist2(istart:iend)
687 tmp_3(:) = ilist3(istart:iend)
688 tmp_4(:) = ilist4(istart:iend)
689 CALL reorder_list_array(tmp_2, tmp_3, tmp_4, nsize=nsize - 1, ndim=ldim)
690 ilist2(istart:iend) = tmp_2(:)
691 ilist3(istart:iend) = tmp_3(:)
692 ilist4(istart:iend) = tmp_4(:)
693 DEALLOCATE (tmp_2)
694 DEALLOCATE (tmp_3)
695 DEALLOCATE (tmp_4)
696 END SELECT
697
698 END SUBROUTINE reorder_list_array_low
699
700! **************************************************************************************************
701!> \brief ...
702!> \param icheck ...
703!> \param bond_list ...
704!> \param i ...
705!> \param mol_natom ...
706!> \param mol_map ...
707!> \param my_mol ...
708!> \author Teodoro Laino 09.2006
709! **************************************************************************************************
710 RECURSIVE SUBROUTINE give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
711 LOGICAL, DIMENSION(:), POINTER :: icheck
712 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
713 INTEGER, INTENT(IN) :: i
714 INTEGER, INTENT(INOUT) :: mol_natom
715 INTEGER, DIMENSION(:), POINTER :: mol_map
716 INTEGER, INTENT(IN) :: my_mol
717
718 INTEGER :: j, k
719
720 IF (mol_map(i) == my_mol) THEN
721 icheck(i) = .true.
722 DO j = 1, SIZE(bond_list(i)%array1)
723 k = bond_list(i)%array1(j)
724 IF (icheck(k)) cycle
725 mol_natom = mol_natom + 1
726 CALL give_back_molecule(icheck, bond_list, k, mol_natom, mol_map, my_mol)
727 END DO
728 ELSE
729 ! Do nothing means only that bonds were found between molecules
730 ! as we defined them.. This could easily be a bond detected but not
731 ! physically present.. so just skip this part and go on..
732 END IF
733 END SUBROUTINE give_back_molecule
734
735! **************************************************************************************************
736!> \brief gives back a mapping of molecules.. icheck needs to be initialized with -1
737!> \param icheck ...
738!> \param bond_list ...
739!> \param i ...
740!> \param my_mol ...
741!> \author Teodoro Laino 04.2007 - Zurich University
742! **************************************************************************************************
743 RECURSIVE SUBROUTINE tag_molecule(icheck, bond_list, i, my_mol)
744 INTEGER, DIMENSION(:), POINTER :: icheck
745 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
746 INTEGER, INTENT(IN) :: i, my_mol
747
748 INTEGER :: j, k
749
750 icheck(i) = my_mol
751 DO j = 1, SIZE(bond_list(i)%array1)
752 k = bond_list(i)%array1(j)
753 IF (k <= i) cycle
754 CALL tag_molecule(icheck, bond_list, k, my_mol)
755 END DO
756
757 END SUBROUTINE tag_molecule
758
759! **************************************************************************************************
760!> \brief Given two lists of corresponding element a complex type is built in
761!> which each element of the type has a list of mapping elements
762!> \param work ...
763!> \param list1 ...
764!> \param list2 ...
765!> \param N ...
766!> \author Teodoro Laino 08.2006
767! **************************************************************************************************
768 SUBROUTINE reorder_structure1d(work, list1, list2, N)
769 TYPE(array1_list_type), DIMENSION(:), &
770 INTENT(INOUT) :: work
771 INTEGER, DIMENSION(:), INTENT(IN) :: list1, list2
772 INTEGER, INTENT(IN) :: N
773
774 INTEGER :: I, index1, index2, Nsize
775 INTEGER, DIMENSION(:), POINTER :: wrk_tmp
776
777 DO i = 1, n
778 index1 = list1(i)
779 index2 = list2(i)
780
781 wrk_tmp => work(index1)%array1
782 nsize = SIZE(wrk_tmp)
783 ALLOCATE (work(index1)%array1(nsize + 1))
784 work(index1)%array1(1:nsize) = wrk_tmp
785 work(index1)%array1(nsize + 1) = index2
786 DEALLOCATE (wrk_tmp)
787
788 wrk_tmp => work(index2)%array1
789 nsize = SIZE(wrk_tmp)
790 ALLOCATE (work(index2)%array1(nsize + 1))
791 work(index2)%array1(1:nsize) = wrk_tmp
792 work(index2)%array1(nsize + 1) = index1
793 DEALLOCATE (wrk_tmp)
794 END DO
795
796 END SUBROUTINE reorder_structure1d
797
798! **************************************************************************************************
799!> \brief Given two lists of corresponding element a complex type is built in
800!> which each element of the type has a list of mapping elements
801!> \param work ...
802!> \param list1 ...
803!> \param list2 ...
804!> \param list3 ...
805!> \param N ...
806!> \author Teodoro Laino 09.2006
807! **************************************************************************************************
808 SUBROUTINE reorder_structure2d(work, list1, list2, list3, N)
809 TYPE(array2_list_type), DIMENSION(:), &
810 INTENT(INOUT) :: work
811 INTEGER, DIMENSION(:), INTENT(IN) :: list1, list2, list3
812 INTEGER, INTENT(IN) :: N
813
814 INTEGER :: I, index1, index2, index3, Nsize
815 INTEGER, DIMENSION(:), POINTER :: wrk_tmp
816
817 DO i = 1, n
818 index1 = list1(i)
819 index2 = list2(i)
820 index3 = list3(i)
821
822 wrk_tmp => work(index1)%array1
823 nsize = SIZE(wrk_tmp)
824 ALLOCATE (work(index1)%array1(nsize + 1))
825 work(index1)%array1(1:nsize) = wrk_tmp
826 work(index1)%array1(nsize + 1) = index2
827 DEALLOCATE (wrk_tmp)
828
829 wrk_tmp => work(index2)%array1
830 nsize = SIZE(wrk_tmp)
831 ALLOCATE (work(index2)%array1(nsize + 1))
832 work(index2)%array1(1:nsize) = wrk_tmp
833 work(index2)%array1(nsize + 1) = index1
834 DEALLOCATE (wrk_tmp)
835
836 wrk_tmp => work(index1)%array2
837 nsize = SIZE(wrk_tmp)
838 ALLOCATE (work(index1)%array2(nsize + 1))
839 work(index1)%array2(1:nsize) = wrk_tmp
840 work(index1)%array2(nsize + 1) = index3
841 DEALLOCATE (wrk_tmp)
842
843 wrk_tmp => work(index2)%array2
844 nsize = SIZE(wrk_tmp)
845 ALLOCATE (work(index2)%array2(nsize + 1))
846 work(index2)%array2(1:nsize) = wrk_tmp
847 work(index2)%array2(nsize + 1) = -index3
848 DEALLOCATE (wrk_tmp)
849 END DO
850
851 END SUBROUTINE reorder_structure2d
852
853! **************************************************************************************************
854!> \brief each atom will be assigned a molecule number based on bonded fragments
855!> The array mol_info should be initialized with -1 before calling the
856!> find_molecule routine
857!> \param atom_bond_list ...
858!> \param mol_info ...
859!> \param mol_name ...
860!> \author Joost 05.2006
861! **************************************************************************************************
862 SUBROUTINE find_molecule(atom_bond_list, mol_info, mol_name)
863 TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
864 INTEGER, DIMENSION(:), POINTER :: mol_info, mol_name
865
866 INTEGER :: i, my_mol_name, n, nmol
867
868 n = SIZE(atom_bond_list)
869 nmol = 0
870 DO i = 1, n
871 IF (mol_info(i) == -1) THEN
872 nmol = nmol + 1
873 my_mol_name = mol_name(i)
874 CALL spread_mol(atom_bond_list, mol_info, i, nmol, my_mol_name, &
875 mol_name)
876 END IF
877 END DO
878 END SUBROUTINE find_molecule
879
880! **************************************************************************************************
881!> \brief spreads the molnumber over the bonded list
882!> \param atom_bond_list ...
883!> \param mol_info ...
884!> \param iatom ...
885!> \param imol ...
886!> \param my_mol_name ...
887!> \param mol_name ...
888!> \author Joost 05.2006
889! **************************************************************************************************
890 RECURSIVE SUBROUTINE spread_mol(atom_bond_list, mol_info, iatom, imol, &
891 my_mol_name, mol_name)
892 TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
893 INTEGER, DIMENSION(:), POINTER :: mol_info
894 INTEGER, INTENT(IN) :: iatom, imol, my_mol_name
895 INTEGER, DIMENSION(:), POINTER :: mol_name
896
897 INTEGER :: atom_b, i
898
899 mol_info(iatom) = imol
900 DO i = 1, SIZE(atom_bond_list(iatom)%array1)
901 atom_b = atom_bond_list(iatom)%array1(i)
902 ! In this way we're really sure that all atoms belong to the same
903 ! molecule. This should correct possible errors in the generation of
904 ! the bond list..
905 IF (mol_info(atom_b) == -1 .AND. my_mol_name == mol_name(atom_b)) THEN
906 CALL spread_mol(atom_bond_list, mol_info, atom_b, imol, my_mol_name, mol_name)
907 IF (mol_info(atom_b) /= imol) cpabort("internal error")
908 END IF
909 END DO
910 END SUBROUTINE spread_mol
911
912! **************************************************************************************************
913!> \brief Use info from periodic table and set atm_mass
914!> \param topology ...
915!> \param subsys_section ...
916! **************************************************************************************************
917 SUBROUTINE topology_set_atm_mass(topology, subsys_section)
918 TYPE(topology_parameters_type), INTENT(INOUT) :: topology
919 TYPE(section_vals_type), POINTER :: subsys_section
920
921 CHARACTER(len=*), PARAMETER :: routinen = 'topology_set_atm_mass'
922
923 CHARACTER(LEN=2) :: upper_sym_1
924 CHARACTER(LEN=default_string_length) :: atmname_upper
925 CHARACTER(LEN=default_string_length), &
926 DIMENSION(:), POINTER :: keyword
927 INTEGER :: handle, i, i_rep, iatom, ielem_found, &
928 iw, n_rep, natom
929 LOGICAL :: user_defined
930 REAL(kind=dp), DIMENSION(:), POINTER :: mass
931 TYPE(atom_info_type), POINTER :: atom_info
932 TYPE(cp_logger_type), POINTER :: logger
933 TYPE(section_vals_type), POINTER :: kind_section
934
935 NULLIFY (logger)
936 logger => cp_get_default_logger()
937 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
938 extension=".subsysLog")
939 CALL timeset(routinen, handle)
940
941 atom_info => topology%atom_info
942 natom = topology%natoms
943
944 ! Available external info
945 kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
946 CALL section_vals_get(kind_section, n_repetition=n_rep)
947 ALLOCATE (keyword(n_rep))
948 ALLOCATE (mass(n_rep))
949 mass = huge(0.0_dp)
950 DO i_rep = 1, n_rep
951 CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
952 c_val=keyword(i_rep), i_rep_section=i_rep)
953 CALL uppercase(keyword(i_rep))
954 CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
955 keyword_name="MASS", n_rep_val=i)
956 IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
957 keyword_name="MASS", r_val=mass(i_rep))
958 END DO
959 !
960 DO iatom = 1, natom
961 !If we reach this point then we've definitely identified the element..
962 !Let's look if an external mass has been defined..
963 user_defined = .false.
964 DO i = 1, SIZE(keyword)
965 atmname_upper = id2str(atom_info%id_atmname(iatom))
966 CALL uppercase(atmname_upper)
967 IF (trim(atmname_upper) == trim(keyword(i)) .AND. mass(i) /= huge(0.0_dp)) THEN
968 atom_info%atm_mass(iatom) = mass(i)
969 user_defined = .true.
970 EXIT
971 END IF
972 END DO
973 ! If name didn't match let's try with the element
974 IF (.NOT. user_defined) THEN
975 upper_sym_1 = trim(id2str(atom_info%id_element(iatom)))
976 CALL get_ptable_info(symbol=upper_sym_1, ielement=ielem_found, amass=atom_info%atm_mass(iatom))
977 END IF
978 IF (iw > 0) WRITE (iw, '(7X,A,A5,A,F12.5)') "In topology_set_atm_mass :: element = ", &
979 id2str(atom_info%id_element(iatom)), " a_mass ", atom_info%atm_mass(iatom)
980 END DO
981 DEALLOCATE (keyword)
982 DEALLOCATE (mass)
983
984 CALL timestop(handle)
985 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
986 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
987
988 END SUBROUTINE topology_set_atm_mass
989
990! **************************************************************************************************
991!> \brief Check and verify that all molecules of the same kind are bonded the same
992!> \param topology ...
993!> \param subsys_section ...
994! **************************************************************************************************
995 SUBROUTINE topology_molecules_check(topology, subsys_section)
996 TYPE(topology_parameters_type), INTENT(INOUT) :: topology
997 TYPE(section_vals_type), POINTER :: subsys_section
998
999 CHARACTER(len=*), PARAMETER :: routinen = 'topology_molecules_check'
1000
1001 INTEGER :: counter, first, first_loc, handle, i, &
1002 iatom, iw, k, loc_counter, mol_num, &
1003 mol_typ, n, natom
1004 LOGICAL :: icheck_num, icheck_typ
1005 TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
1006 TYPE(atom_info_type), POINTER :: atom_info
1007 TYPE(connectivity_info_type), POINTER :: conn_info
1008 TYPE(cp_logger_type), POINTER :: logger
1009
1010 CALL timeset(routinen, handle)
1011
1012 NULLIFY (logger)
1013 logger => cp_get_default_logger()
1014 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
1015 extension=".subsysLog")
1016
1017 atom_info => topology%atom_info
1018 conn_info => topology%conn_info
1019 natom = topology%natoms
1020
1021 IF (iw > 0) THEN
1022 WRITE (unit=iw, fmt="(/,T2,A)") &
1023 "FORCEFIELD| Checking consistency of the generated molecules"
1024 END IF
1025
1026 ALLOCATE (atom_bond_list(natom))
1027
1028 DO i = 1, natom
1029 ALLOCATE (atom_bond_list(i)%array1(0))
1030 END DO
1031 n = 0
1032 IF (ASSOCIATED(conn_info%bond_a)) n = SIZE(conn_info%bond_a)
1033 CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, n)
1034
1035 mol_typ = atom_info%map_mol_typ(1)
1036 mol_num = atom_info%map_mol_num(1)
1037 counter = 1
1038 loc_counter = 1
1039 first = 1
1040 first_loc = 1
1041 DO iatom = 2, natom
1042 icheck_num = (atom_info%map_mol_num(iatom) == mol_num)
1043 icheck_typ = (atom_info%map_mol_typ(iatom) == mol_typ)
1044 IF ((icheck_typ .AND. (.NOT. icheck_num)) .OR. (.NOT. icheck_typ)) THEN
1045 !-----------------------------------------------------------------------------
1046 ! 1. Check each molecule have the same number of atoms
1047 !-----------------------------------------------------------------------------
1048 IF (counter /= loc_counter) THEN
1049 CALL cp_abort(__location__, &
1050 "Different number of atoms for the same molecule kind found "// &
1051 "(molecule type = "//trim(adjustl(cp_to_string(mol_typ)))// &
1052 ", molecule number = "//trim(adjustl(cp_to_string(mol_num)))// &
1053 ", expected number of atoms = "// &
1054 trim(adjustl(cp_to_string(counter)))// &
1055 ", number of atoms found = "// &
1056 trim(adjustl(cp_to_string(loc_counter)))//").")
1057 END IF
1058 END IF
1059 IF (.NOT. icheck_typ) THEN
1060 first = iatom
1061 first_loc = iatom
1062 counter = 1
1063 loc_counter = 1
1064 mol_typ = atom_info%map_mol_typ(iatom)
1065 END IF
1066 IF (icheck_num) THEN
1067 IF (icheck_typ) loc_counter = loc_counter + 1
1068 !-----------------------------------------------------------------------------
1069 ! 2. Check that each molecule has the same atom sequences
1070 !-----------------------------------------------------------------------------
1071 IF (atom_info%id_atmname(iatom) /= &
1072 atom_info%id_atmname(first + loc_counter - 1)) THEN
1073 CALL cp_abort(__location__, &
1074 "Different atom names for the same molecule kind found "// &
1075 "(atom number = "//trim(adjustl(cp_to_string(iatom)))// &
1076 ", molecule type = "//trim(adjustl(cp_to_string(mol_typ)))// &
1077 ", molecule number = "//trim(adjustl(cp_to_string(mol_num)))// &
1078 ", expected atom name = "// &
1079 trim(adjustl(id2str(atom_info%id_atmname(first + loc_counter - 1))))// &
1080 ", atom name found = "// &
1081 trim(adjustl(id2str(atom_info%id_atmname(iatom))))//").")
1082 END IF
1083 !-----------------------------------------------------------------------------
1084 ! 3. Check that each molecule have the same bond sequences
1085 !-----------------------------------------------------------------------------
1086 IF (SIZE(atom_bond_list(iatom)%array1) /= SIZE(atom_bond_list(first + loc_counter - 1)%array1)) THEN
1087 CALL cp_abort(__location__, &
1088 "Different number of bonds for the same molecule kind found "// &
1089 "(molecule type = "//trim(adjustl(cp_to_string(mol_typ)))// &
1090 ", molecule number = "//trim(adjustl(cp_to_string(mol_num)))// &
1091 ", expected number of bonds = "// &
1092 trim(adjustl(cp_to_string(SIZE(atom_bond_list(first + loc_counter - 1)%array1))))// &
1093 ", number of bonds found = "// &
1094 trim(adjustl(cp_to_string(SIZE(atom_bond_list(iatom)%array1))))// &
1095 "). Check the connectivity of your system.")
1096 END IF
1097 DO k = 1, SIZE(atom_bond_list(iatom)%array1)
1098 IF (all(atom_bond_list(first + loc_counter - 1)%array1 - first /= &
1099 atom_bond_list(iatom)%array1(k) - first_loc)) THEN
1100 CALL cp_abort(__location__, &
1101 "Different sequence of bonds for the same molecule kind found "// &
1102 "(molecule type = "//trim(adjustl(cp_to_string(mol_typ)))// &
1103 ", molecule number = "//trim(adjustl(cp_to_string(mol_num)))// &
1104 "). Check the connectivity of your system.")
1105 END IF
1106 END DO
1107 ELSE
1108 mol_num = atom_info%map_mol_num(iatom)
1109 loc_counter = 1
1110 first_loc = iatom
1111 END IF
1112 IF (mol_num == 1 .AND. icheck_typ) counter = counter + 1
1113 END DO
1114
1115 IF (iw > 0) THEN
1116 WRITE (unit=iw, fmt="(T2,A)") &
1117 "FORCEFIELD| Consistency check for molecules completed"
1118 END IF
1119
1120 DO i = 1, natom
1121 DEALLOCATE (atom_bond_list(i)%array1)
1122 END DO
1123 DEALLOCATE (atom_bond_list)
1124
1125 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1126 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
1127
1128 CALL timestop(handle)
1129
1130 END SUBROUTINE topology_molecules_check
1131
1132! **************************************************************************************************
1133!> \brief Check and returns the ELEMENT label
1134!> \param element_in ...
1135!> \param atom_name_in ...
1136!> \param element_out ...
1137!> \param subsys_section ...
1138!> \param use_mm_map_first ...
1139!> \par History
1140!> 12.2005 created [teo]
1141!> \author Teodoro Laino
1142! **************************************************************************************************
1143 SUBROUTINE check_subsys_element(element_in, atom_name_in, element_out, subsys_section, use_mm_map_first)
1144 CHARACTER(len=*), INTENT(IN) :: element_in, atom_name_in
1145 CHARACTER(len=default_string_length), INTENT(OUT) :: element_out
1146 TYPE(section_vals_type), POINTER :: subsys_section
1147 LOGICAL :: use_mm_map_first
1148
1149 CHARACTER(len=default_string_length) :: atom_name, element_symbol, keyword
1150 INTEGER :: i, i_rep, n_rep
1151 LOGICAL :: defined_kind_section, found, in_ptable
1152 TYPE(section_vals_type), POINTER :: kind_section
1153
1154 found = .false.
1155 element_symbol = element_in
1156 atom_name = atom_name_in
1157 element_out = ""
1158 defined_kind_section = .false.
1159
1160 ! First check if a KIND section is overriding the element
1161 ! definition
1162 CALL uppercase(atom_name)
1163 kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
1164 CALL section_vals_get(kind_section, n_repetition=n_rep)
1165 DO i_rep = 1, n_rep
1166 CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1167 c_val=keyword, i_rep_section=i_rep)
1168 CALL uppercase(keyword)
1169 IF (trim(keyword) == trim(atom_name)) THEN
1170 CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1171 keyword_name="ELEMENT", n_rep_val=i)
1172 IF (i > 0) THEN
1173 CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1174 keyword_name="ELEMENT", c_val=element_symbol)
1175 defined_kind_section = .true.
1176 EXIT
1177 ELSE
1178 element_symbol = element_in
1179 defined_kind_section = .true.
1180 END IF
1181 END IF
1182 END DO
1183
1184 ! Let's check the validity of the element so far stored..
1185 ! if we are not having a connectivity file, we must first match against the ptable.
1186 ! this helps to resolve Ca/CA (calcium and Calpha) or Cn/CN7 (Coppernicum (112) CN) conflicts
1187 ! so, in the presence of connectivity CA will be 'C', while in the absence of connectivity CA will be 'Ca'
1188 IF (defined_kind_section .OR. .NOT. use_mm_map_first) THEN
1189 ! lengths larger than 2 should not match, because 'trailing' characters are ignored.
1190 in_ptable = .false.
1191 IF (len_trim(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1192 IF (in_ptable) THEN
1193 element_out = trim(element_symbol)
1194 found = .true.
1195 END IF
1196 END IF
1197
1198 ! This is clearly a user error
1199 IF (.NOT. found .AND. defined_kind_section) &
1200 CALL cp_abort(__location__, "Element <"//trim(element_symbol)// &
1201 "> provided for KIND <"//trim(atom_name_in)//"> "// &
1202 "which cannot be mapped with any standard element label. Please correct your "// &
1203 "input file!")
1204
1205 ! Last chance.. are these atom_kinds of AMBER or CHARMM or GROMOS FF ?
1206 CALL uppercase(element_symbol)
1207 IF ((.NOT. found) .AND. (ASSOCIATED(amber_map))) THEN
1208 ! First we go through the AMBER library
1209 DO i = 1, SIZE(amber_map%kind)
1210 IF (element_symbol == amber_map%kind(i)) THEN
1211 found = .true.
1212 EXIT
1213 END IF
1214 END DO
1215 IF (found) THEN
1216 element_out = amber_map%element(i)
1217 END IF
1218 END IF
1219 IF ((.NOT. found) .AND. (ASSOCIATED(charmm_map))) THEN
1220 ! Then we go through the CHARMM library
1221 DO i = 1, SIZE(charmm_map%kind)
1222 IF (element_symbol == charmm_map%kind(i)) THEN
1223 found = .true.
1224 EXIT
1225 END IF
1226 END DO
1227 IF (found) THEN
1228 element_out = charmm_map%element(i)
1229 END IF
1230 END IF
1231 IF ((.NOT. found) .AND. (ASSOCIATED(gromos_map))) THEN
1232 ! Then we go through the GROMOS library
1233 DO i = 1, SIZE(gromos_map%kind)
1234 IF (element_symbol == gromos_map%kind(i)) THEN
1235 found = .true.
1236 EXIT
1237 END IF
1238 END DO
1239 IF (found) THEN
1240 element_out = gromos_map%element(i)
1241 END IF
1242 END IF
1243
1244 ! final check. We has a connectivity, so we first tried to match against the ff_maps, but the element was not there.
1245 ! Try again the periodic table.
1246 IF (.NOT. found) THEN
1247 in_ptable = .false.
1248 IF (len_trim(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1249 IF (in_ptable) THEN
1250 element_out = trim(element_symbol)
1251 found = .true.
1252 END IF
1253 END IF
1254
1255 ! If no element is found the job stops here.
1256 IF (.NOT. found) THEN
1257 CALL cp_abort(__location__, &
1258 " Unknown element for KIND <"//trim(atom_name_in)//">."// &
1259 " This problem can be fixed specifying properly elements in PDB"// &
1260 " or specifying a KIND section or getting in touch with one of"// &
1261 " the developers!")
1262 END IF
1263
1264 END SUBROUTINE check_subsys_element
1265
1266END MODULE topology_util
graph
program graph
Program to Map on grid the hills spawned during a metadynamics run.
Definition graph.F:19
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
topology_util::reorder_structure
Definition topology_util.F:71
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
graphcon
uses a combination of graphs and hashing to determine if two molecules are topologically equivalent,...
Definition graphcon.F:23
graphcon::hash_molecule
subroutine, public hash_molecule(reference, kind_ref, hash)
hashes a molecule to a number. Molecules that are the (topologically) the same have the same hash....
Definition graphcon.F:80
graphcon::reorder_graph
subroutine, public reorder_graph(reference, unordered, order, matches)
If two molecules are topologically the same, finds the ordering that maps the unordered one on the or...
Definition graphcon.F:134
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1223
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mm_mapping_library
Contains the mapping ATOM_KIND -> ELEMENT for the most common cases in CHARMM and AMBER This should a...
Definition mm_mapping_library.F:18
mm_mapping_library::charmm_map
type(ff_map_type), pointer, public charmm_map
Definition mm_mapping_library.F:34
mm_mapping_library::gromos_map
type(ff_map_type), pointer, public gromos_map
Definition mm_mapping_library.F:34
mm_mapping_library::amber_map
type(ff_map_type), pointer, public amber_map
Definition mm_mapping_library.F:34
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
qmmm_types_low
Definition qmmm_types_low.F:13
string_table
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition string_table.F:22
string_table::str2id
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition string_table.F:72
string_table::id2str
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition string_table.F:115
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
topology_types
Definition topology_types.F:18
topology_util
Collection of subroutine needed for topology related things.
Definition topology_util.F:13
topology_util::topology_molecules_check
subroutine, public topology_molecules_check(topology, subsys_section)
Check and verify that all molecules of the same kind are bonded the same.
Definition topology_util.F:996
topology_util::find_molecule
subroutine, public find_molecule(atom_bond_list, mol_info, mol_name)
each atom will be assigned a molecule number based on bonded fragments The array mol_info should be i...
Definition topology_util.F:863
topology_util::give_back_molecule
recursive subroutine, public give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
...
Definition topology_util.F:711
topology_util::topology_reorder_atoms
subroutine, public topology_reorder_atoms(topology, qmmm, qmmm_env_mm, subsys_section, force_env_section)
...
Definition topology_util.F:88
topology_util::reorder_list_array
recursive subroutine, public reorder_list_array(ilist1, ilist2, ilist3, ilist4, nsize, ndim)
Order arrays of lists.
Definition topology_util.F:595
topology_util::topology_set_atm_mass
subroutine, public topology_set_atm_mass(topology, subsys_section)
Use info from periodic table and set atm_mass.
Definition topology_util.F:918
topology_util::tag_molecule
recursive subroutine, public tag_molecule(icheck, bond_list, i, my_mol)
gives back a mapping of molecules.. icheck needs to be initialized with -1
Definition topology_util.F:744
topology_util::check_subsys_element
subroutine, public check_subsys_element(element_in, atom_name_in, element_out, subsys_section, use_mm_map_first)
Check and returns the ELEMENT label.
Definition topology_util.F:1144
topology
Control for reading in different topologies and coordinates.
Definition topology.F:13
util
All kind of helpful little routines.
Definition util.F:14
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
graphcon::graph_type
Definition graphcon.F:37
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
qmmm_types_low::qmmm_env_mm_type
...
Definition qmmm_types_low.F:238
topology_types::atom_info_type
Definition topology_types.F:35
topology_types::connectivity_info_type
Definition topology_types.F:53
topology_types::topology_parameters_type
Definition topology_types.F:143
topology_util::array1_list_type
Definition topology_util.F:49
topology_util::array2_list_type
Definition topology_util.F:53