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