(git:0de0cc2)
topology_generate_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 !> Teodor Laino 09.2006 - Major rewriting with linear scaling routines
12 ! **************************************************************************************************
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
17  USE cell_types, ONLY: pbc
20  cp_logger_type,&
21  cp_to_string
24  USE cp_units, ONLY: cp_unit_to_cp2k
26  fist_neighbor_type
28  USE input_constants, ONLY: do_add,&
31  do_conn_off,&
32  do_conn_user,&
33  do_remove
36  section_vals_type,&
38  USE kinds, ONLY: default_string_length,&
39  dp
40  USE memory_utilities, ONLY: reallocate
41  USE message_passing, ONLY: mp_para_env_type
44  particle_type
46  USE qmmm_types_low, ONLY: qmmm_env_mm_type
47  USE string_table, ONLY: id2str,&
48  s2s,&
49  str2id
51  uppercase
52  USE topology_types, ONLY: atom_info_type,&
55  USE topology_util, ONLY: array1_list_type,&
60  reorder_structure
61  USE util, ONLY: find_boundary,&
62  sort
63 #include "./base/base_uses.f90"
64 
65  IMPLICIT NONE
66 
67  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_generate_util'
68 
69  PRIVATE
70  LOGICAL, PARAMETER :: debug_this_module = .false.
71 
72  PUBLIC :: topology_generate_bend, &
80 
81 CONTAINS
82 
83 ! **************************************************************************************************
84 !> \brief Generates molnames: useful when the connectivity on file does not
85 !> provide them
86 !> \param conn_info ...
87 !> \param natom ...
88 !> \param natom_prev ...
89 !> \param nbond_prev ...
90 !> \param id_molname ...
91 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
92 ! **************************************************************************************************
93  SUBROUTINE topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, &
94  id_molname)
95  TYPE(connectivity_info_type), POINTER :: conn_info
96  INTEGER, INTENT(IN) :: natom, natom_prev, nbond_prev
97  INTEGER, DIMENSION(:), INTENT(INOUT) :: id_molname
98 
99  CHARACTER(LEN=default_string_length), PARAMETER :: basename = "MOL"
100 
101  CHARACTER(LEN=default_string_length) :: molname
102  INTEGER :: i, n, nmol
103  LOGICAL :: check
104  TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
105 
106 ! convert a simple list of bonds to a list of bonds per atom
107 ! (each bond is present in the forward and backward direction
108 
109  ALLOCATE (atom_bond_list(natom))
110  DO i = 1, natom
111  ALLOCATE (atom_bond_list(i)%array1(0))
112  END DO
113  n = 0
114  IF (ASSOCIATED(conn_info%bond_a)) n = SIZE(conn_info%bond_a) - nbond_prev
115  CALL reorder_structure(atom_bond_list, conn_info%bond_a(nbond_prev + 1:) - natom_prev, &
116  conn_info%bond_b(nbond_prev + 1:) - natom_prev, n)
117 
118  nmol = 0
119  check = all(id_molname == str2id(s2s("__UNDEF__"))) .OR. all(id_molname /= str2id(s2s("__UNDEF__")))
120  cpassert(check)
121  DO i = 1, natom
122  IF (id2str(id_molname(i)) == "__UNDEF__") THEN
123  molname = trim(basename)//adjustl(cp_to_string(nmol))
124  CALL generate_molname_low(i, atom_bond_list, molname, id_molname)
125  nmol = nmol + 1
126  END IF
127  END DO
128  DO i = 1, natom
129  DEALLOCATE (atom_bond_list(i)%array1)
130  END DO
131  DEALLOCATE (atom_bond_list)
132 
133  END SUBROUTINE topology_generate_molname
134 
135 ! **************************************************************************************************
136 !> \brief Generates molnames: useful when the connectivity on file does not
137 !> provide them
138 !> \param i ...
139 !> \param atom_bond_list ...
140 !> \param molname ...
141 !> \param id_molname ...
142 !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
143 ! **************************************************************************************************
144  RECURSIVE SUBROUTINE generate_molname_low(i, atom_bond_list, molname, id_molname)
145  INTEGER, INTENT(IN) :: i
146  TYPE(array1_list_type), DIMENSION(:) :: atom_bond_list
147  CHARACTER(LEN=default_string_length), INTENT(IN) :: molname
148  INTEGER, DIMENSION(:), INTENT(INOUT) :: id_molname
149 
150  INTEGER :: j, k
151 
152  IF (debug_this_module) THEN
153  WRITE (*, *) "Entered with :", i
154  WRITE (*, *) trim(molname)//": entering with i:", i, " full series to test:: ", atom_bond_list(i)%array1
155  IF ((trim(id2str(id_molname(i))) /= "__UNDEF__") .AND. &
156  (trim(id2str(id_molname(i))) /= trim(molname))) THEN
157  WRITE (*, *) "Atom (", i, ") has already a molecular name assigned ! ("//trim(id2str(id_molname(i)))//")."
158  WRITE (*, *) "New molecular name would be: ("//trim(molname)//")."
159  WRITE (*, *) "Detecting something wrong in the molecular setup!"
160  cpabort("")
161  END IF
162  END IF
163  id_molname(i) = str2id(molname)
164  DO j = 1, SIZE(atom_bond_list(i)%array1)
165  k = atom_bond_list(i)%array1(j)
166  IF (debug_this_module) WRITE (*, *) "entering with i:", i, "testing :", k
167  IF (k == -1) cycle
168  atom_bond_list(i)%array1(j) = -1
169  WHERE (atom_bond_list(k)%array1 == i) atom_bond_list(k)%array1 = -1
170  CALL generate_molname_low(k, atom_bond_list, molname, id_molname)
171  END DO
172  END SUBROUTINE generate_molname_low
173 
174 ! **************************************************************************************************
175 !> \brief Use information from bond list to generate molecule. (ie clustering)
176 !> \param topology ...
177 !> \param qmmm ...
178 !> \param qmmm_env ...
179 !> \param subsys_section ...
180 ! **************************************************************************************************
181  SUBROUTINE topology_generate_molecule(topology, qmmm, qmmm_env, subsys_section)
182  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
183  LOGICAL, INTENT(in), OPTIONAL :: qmmm
184  TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
185  TYPE(section_vals_type), POINTER :: subsys_section
186 
187  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_molecule'
188  INTEGER, PARAMETER :: nblock = 100
189 
190  INTEGER :: atom_in_kind, atom_in_mol, first, handle, handle2, i, iatm, iatom, iend, ifirst, &
191  ilast, inum, istart, itype, iw, j, jump1, jump2, last, max_mol_num, mol_num, mol_res, &
192  mol_typ, myind, n, natom, nlocl, ntype, resid
193  INTEGER, DIMENSION(:), POINTER :: qm_atom_index, wrk1, wrk2
194  LOGICAL :: do_again, found, my_qmmm
195  TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
196  TYPE(atom_info_type), POINTER :: atom_info
197  TYPE(connectivity_info_type), POINTER :: conn_info
198  TYPE(cp_logger_type), POINTER :: logger
199 
200  NULLIFY (logger)
201  logger => cp_get_default_logger()
202  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
203  extension=".subsysLog")
204  CALL timeset(routinen, handle)
205  NULLIFY (qm_atom_index)
206  NULLIFY (wrk1)
207  NULLIFY (wrk2)
208 
209  atom_info => topology%atom_info
210  conn_info => topology%conn_info
211  !
212  ! QM/MM coordinate_control
213  !
214  my_qmmm = .false.
215  IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
216 
217  natom = topology%natoms
218  IF (ASSOCIATED(atom_info%map_mol_typ)) DEALLOCATE (atom_info%map_mol_typ)
219  ALLOCATE (atom_info%map_mol_typ(natom))
220 
221  IF (ASSOCIATED(atom_info%map_mol_num)) DEALLOCATE (atom_info%map_mol_num)
222  ALLOCATE (atom_info%map_mol_num(natom))
223 
224  IF (ASSOCIATED(atom_info%map_mol_res)) DEALLOCATE (atom_info%map_mol_res)
225  ALLOCATE (atom_info%map_mol_res(natom))
226 
227  ! Initialisation
228  atom_info%map_mol_typ(:) = 0
229  atom_info%map_mol_num(:) = -1
230  atom_info%map_mol_res(:) = 1
231 
232  ! Parse the atom list to find the different molecule types and residues
233  ntype = 1
234  atom_info%map_mol_typ(1) = 1
235  resid = 1
236  CALL reallocate(wrk1, 1, nblock)
237  wrk1(1) = atom_info%id_molname(1)
238  DO iatom = 2, natom
239  IF (topology%conn_type == do_conn_off) THEN
240  ! No connectivity: each atom becomes a molecule of its own molecule kind
241  ntype = ntype + 1
242  atom_info%map_mol_typ(iatom) = ntype
243  ELSE IF (topology%conn_type == do_conn_user) THEN
244  ! User-defined connectivity: 5th column of COORD section or molecule
245  ! or residue name in the case of PDB files
246  IF (atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1)) THEN
247  atom_info%map_mol_typ(iatom) = atom_info%map_mol_typ(iatom - 1)
248  IF (atom_info%id_resname(iatom) == atom_info%id_resname(iatom - 1)) THEN
249  atom_info%map_mol_res(iatom) = atom_info%map_mol_res(iatom - 1)
250  ELSE
251  resid = resid + 1
252  atom_info%map_mol_res(iatom) = resid
253  END IF
254  ELSE
255  ! Check if the type is already known
256  found = .false.
257  DO itype = 1, ntype
258  IF (atom_info%id_molname(iatom) == wrk1(itype)) THEN
259  atom_info%map_mol_typ(iatom) = itype
260  found = .true.
261  EXIT
262  END IF
263  END DO
264  IF (.NOT. found) THEN
265  ntype = ntype + 1
266  atom_info%map_mol_typ(iatom) = ntype
267  IF (ntype > SIZE(wrk1)) CALL reallocate(wrk1, 1, 2*SIZE(wrk1))
268  wrk1(ntype) = atom_info%id_molname(iatom)
269  END IF
270  resid = resid + 1
271  atom_info%map_mol_res(iatom) = resid
272  END IF
273  ELSE
274  IF (atom_info%id_molname(iatom - 1) == atom_info%id_molname(iatom)) THEN
275  atom_info%map_mol_typ(iatom) = ntype
276  ELSE
277  ntype = ntype + 1
278  atom_info%map_mol_typ(iatom) = ntype
279  END IF
280  END IF
281  END DO
282  DEALLOCATE (wrk1)
283 
284  IF (iw > 0) WRITE (iw, '(/,T2,A)') "Start of molecule generation"
285 
286  ! convert a simple list of bonds to a list of bonds per atom
287  ! (each bond is present in the forward and backward direction
288  ALLOCATE (atom_bond_list(natom))
289  DO i = 1, natom
290  ALLOCATE (atom_bond_list(i)%array1(0))
291  END DO
292  n = 0
293  IF (ASSOCIATED(conn_info%bond_a)) n = SIZE(conn_info%bond_a)
294  CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, n)
295  CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
296  DO i = 1, natom
297  DEALLOCATE (atom_bond_list(i)%array1)
298  END DO
299  DEALLOCATE (atom_bond_list)
300  IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of molecule generation"
301 
302  ! Modify according map_mol_typ the array map_mol_num
303  IF (iw > 0) WRITE (iw, '(/,T2,A)') "Checking for non-continuous generated molecules"
304  ! Check molecule number
305  ALLOCATE (wrk1(natom))
306  ALLOCATE (wrk2(natom))
307  wrk1 = atom_info%map_mol_num
308 
309  IF (debug_this_module) THEN
310  DO i = 1, natom
311  WRITE (*, '(2I10)') i, atom_info%map_mol_num(i)
312  END DO
313  END IF
314 
315  CALL sort(wrk1, natom, wrk2)
316  istart = 1
317  mol_typ = wrk1(istart)
318  DO i = 2, natom
319  IF (mol_typ /= wrk1(i)) THEN
320  iend = i - 1
321  first = minval(wrk2(istart:iend))
322  last = maxval(wrk2(istart:iend))
323  nlocl = last - first + 1
324  IF (iend - istart + 1 /= nlocl) THEN
325  IF (debug_this_module) WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl
326  CALL cp_abort(__location__, &
327  "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// &
328  "In particular a molecule defined from index ("//cp_to_string(first)//") to ("// &
329  cp_to_string(last)//") contains other molecules, not connected! "// &
330  "Too late at this stage everything should be already ordered! "// &
331  "If you have not yet employed the REORDERING keyword, please do so. "// &
332  "It may help to fix this issue.")
333  END IF
334  istart = i
335  mol_typ = wrk1(istart)
336  END IF
337  END DO
338  iend = i - 1
339  first = minval(wrk2(istart:iend))
340  last = maxval(wrk2(istart:iend))
341  nlocl = last - first + 1
342  IF (iend - istart + 1 /= nlocl) THEN
343  IF (debug_this_module) WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl
344  CALL cp_abort(__location__, &
345  "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// &
346  "In particular a molecule defined from index ("//cp_to_string(first)//") to ("// &
347  cp_to_string(last)//") contains other molecules, not connected! "// &
348  "Too late at this stage everything should be already ordered! "// &
349  "If you have not yet employed the REORDERING keyword, please do so. "// &
350  "It may help to fix this issue.")
351  END IF
352  DEALLOCATE (wrk1)
353  DEALLOCATE (wrk2)
354  IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of check"
355 
356  IF (iw > 0) WRITE (unit=iw, fmt="(/,T2,A)") "Start of renumbering molecules"
357  IF (topology%conn_type == do_conn_user) THEN
358  mol_num = 1
359  atom_info%map_mol_num(1) = 1
360  DO iatom = 2, natom
361  IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN
362  mol_num = 1
363  ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN
364  mol_num = mol_num + 1
365  END IF
366  atom_info%map_mol_num(iatom) = mol_num
367  END DO
368  ELSE
369  mol_typ = atom_info%map_mol_typ(1)
370  mol_num = atom_info%map_mol_num(1)
371  DO i = 2, natom
372  IF (atom_info%map_mol_typ(i) /= mol_typ) THEN
373  myind = atom_info%map_mol_num(i) - mol_num + 1
374  cpassert(myind /= atom_info%map_mol_num(i - 1))
375  mol_typ = atom_info%map_mol_typ(i)
376  mol_num = atom_info%map_mol_num(i)
377  END IF
378  atom_info%map_mol_num(i) = atom_info%map_mol_num(i) - mol_num + 1
379  END DO
380  END IF
381  IF (iw > 0) WRITE (unit=iw, fmt="(/,T2,A)") "End of renumbering molecules"
382 
383  ! Optionally, use the residues as molecules
384  CALL timeset(routinen//"_PARA_RES", handle2)
385  IF (iw > 0) WRITE (unit=iw, fmt="(/,T2,A,L2)") "Starting PARA_RES: ", topology%para_res
386  IF (topology%para_res) THEN
387  IF (topology%conn_type == do_conn_user) THEN
388  atom_info%id_molname(:) = atom_info%id_resname(:)
389  ntype = 1
390  atom_info%map_mol_typ(1) = 1
391  mol_num = 1
392  atom_info%map_mol_num(1) = 1
393  DO iatom = 2, natom
394  IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN
395  ntype = ntype + 1
396  mol_num = 1
397  ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN
398  mol_num = mol_num + 1
399  END IF
400  atom_info%map_mol_typ(iatom) = ntype
401  atom_info%map_mol_num(iatom) = mol_num
402  END DO
403  ELSE
404  mol_res = 1
405  mol_typ = atom_info%map_mol_typ(1)
406  mol_num = atom_info%map_mol_num(1)
407  atom_info%map_mol_res(1) = mol_res
408  DO i = 2, natom
409  IF ((atom_info%resid(i - 1) /= atom_info%resid(i)) .OR. &
410  (atom_info%id_resname(i - 1) /= atom_info%id_resname(i))) THEN
411  mol_res = mol_res + 1
412  END IF
413  IF ((atom_info%map_mol_typ(i) /= mol_typ) .OR. &
414  (atom_info%map_mol_num(i) /= mol_num)) THEN
415  mol_typ = atom_info%map_mol_typ(i)
416  mol_num = atom_info%map_mol_num(i)
417  mol_res = 1
418  END IF
419  atom_info%map_mol_res(i) = mol_res
420  END DO
421  END IF
422  END IF
423  IF (iw > 0) WRITE (unit=iw, fmt="(/,T2,A)") "End of PARA_RES"
424  CALL timestop(handle2)
425 
426  IF (iw > 0) THEN
427  DO iatom = 1, natom
428  WRITE (iw, '(4(1X,A,":",I0),2(1X,A,1X,A))') "iatom", iatom, &
429  "map_mol_typ", atom_info%map_mol_typ(iatom), &
430  "map_mol_num", atom_info%map_mol_num(iatom), &
431  "map_mol_res", atom_info%map_mol_res(iatom), &
432  "mol_name:", trim(id2str(atom_info%id_molname(iatom))), &
433  "res_name:", trim(id2str(atom_info%id_resname(iatom)))
434  END DO
435  END IF
436 
437  IF (my_qmmm) THEN
438  do_again = .false.
439  IF (iw > 0) WRITE (iw, *) "MAP_MOL_NUM ", atom_info%map_mol_num
440  IF (iw > 0) WRITE (iw, *) "MAP_MOL_TYP ", atom_info%map_mol_typ
441  IF (iw > 0) WRITE (iw, *) "MAP_MOL_RES ", atom_info%map_mol_res
442  ALLOCATE (qm_atom_index(SIZE(qmmm_env%qm_atom_index)))
443  qm_atom_index = qmmm_env%qm_atom_index
444  cpassert(all(qm_atom_index /= 0))
445  DO myind = 1, SIZE(qm_atom_index)
446  IF (qm_atom_index(myind) == 0) cycle
447  CALL find_boundary(atom_info%map_mol_typ, natom, ifirst, ilast, &
448  atom_info%map_mol_typ(qm_atom_index(myind)))
449  CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, &
450  atom_info%map_mol_typ(qm_atom_index(myind)), atom_info%map_mol_num(qm_atom_index(myind)))
451  IF (iw > 0) WRITE (iw, *) "qm fragment:: ifirst, ilast", ifirst, ilast
452  cpassert(((ifirst /= 0) .OR. (ilast /= natom)))
453  DO iatm = ifirst, ilast
454  atom_info%id_molname(iatm) = str2id(s2s("_QM_"// &
455  trim(id2str(atom_info%id_molname(iatm)))))
456  IF (iw > 0) WRITE (iw, *) "QM Molecule name :: ", id2str(atom_info%id_molname(iatm))
457  WHERE (qm_atom_index == iatm) qm_atom_index = 0
458  END DO
459  DO iatm = 1, ifirst - 1
460  IF (any(qm_atom_index == iatm)) do_again = .true.
461  END DO
462  DO iatm = ilast + 1, natom
463  IF (any(qm_atom_index == iatm)) do_again = .true.
464  END DO
465  IF (iw > 0) WRITE (iw, *) " Another QM fragment? :: ", do_again
466  IF (ifirst /= 1) THEN
467  jump1 = atom_info%map_mol_typ(ifirst) - atom_info%map_mol_typ(ifirst - 1)
468  cpassert(jump1 <= 1 .AND. jump1 >= 0)
469  jump1 = abs(jump1 - 1)
470  ELSE
471  jump1 = 0
472  END IF
473  IF (ilast /= natom) THEN
474  jump2 = atom_info%map_mol_typ(ilast + 1) - atom_info%map_mol_typ(ilast)
475  cpassert(jump2 <= 1 .AND. jump2 >= 0)
476  jump2 = abs(jump2 - 1)
477  ELSE
478  jump2 = 0
479  END IF
480 
481  ! Changing mol_type consistently
482  DO iatm = ifirst, natom
483  atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump1
484  END DO
485  DO iatm = ilast + 1, natom
486  atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump2
487  END DO
488  IF (jump1 == 1) THEN
489  DO iatm = ifirst, ilast
490  atom_info%map_mol_num(iatm) = 1
491  END DO
492  END IF
493 
494  IF (jump2 == 1) THEN
495  CALL find_boundary(atom_info%map_mol_typ, natom, first, last, atom_info%map_mol_typ(ilast + 1))
496  CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, &
497  atom_info%map_mol_typ(ilast + 1), atom_info%map_mol_num(ilast + 1))
498  atom_in_mol = ilast - ifirst + 1
499  inum = 1
500  DO iatm = first, last, atom_in_mol
501  atom_info%map_mol_num(iatm:iatm + atom_in_mol - 1) = inum
502  inum = inum + 1
503  END DO
504  END IF
505 
506  IF (.NOT. do_again) EXIT
507  END DO
508  DEALLOCATE (qm_atom_index)
509 
510  IF (iw > 0) THEN
511  WRITE (iw, *) "After the QM/MM Setup:"
512  DO iatom = 1, natom
513  WRITE (iw, *) " iatom,map_mol_typ,map_mol_num ", iatom, &
514  atom_info%map_mol_typ(iatom), atom_info%map_mol_num(iatom)
515  END DO
516  END IF
517  END IF
518  !
519  ! Further check : see if the number of atoms belonging to same molecule kinds
520  ! are equal
521  IF (iw > 0) THEN
522  WRITE (iw, *) "SUMMARY:: Number of molecule kinds found:", ntype
523  ntype = maxval(atom_info%map_mol_typ)
524  DO i = 1, ntype
525  atom_in_kind = count(atom_info%map_mol_typ == i)
526  WRITE (iw, *) "Molecule kind:", i, " contains", atom_in_kind, " atoms"
527  IF (atom_in_kind <= 1) cycle
528  CALL find_boundary(atom_info%map_mol_typ, natom, first, last, i)
529  WRITE (iw, *) "Boundary atoms:", first, last
530  cpassert(last - first + 1 == atom_in_kind)
531  max_mol_num = maxval(atom_info%map_mol_num(first:last))
532  WRITE (iw, *) "Number of molecules of kind", i, "is ::", max_mol_num
533  atom_in_mol = atom_in_kind/max_mol_num
534  WRITE (iw, *) "Number of atoms per each molecule:", atom_in_mol
535  WRITE (iw, *) "MAP_MOL_TYP::", atom_info%map_mol_typ(first:last)
536  WRITE (iw, *) "MAP_MOL_NUM::", atom_info%map_mol_num(first:last)
537  WRITE (iw, *) "MAP_MOL_RES::", atom_info%map_mol_res(first:last)
538  !
539  DO j = 1, max_mol_num
540  IF (count(atom_info%map_mol_num(first:last) == j) /= atom_in_mol) THEN
541  WRITE (iw, *) "molecule type:", i, "molecule num:", j, " has ", &
542  count(atom_info%map_mol_num(first:last) == j), &
543  " atoms instead of ", atom_in_mol, " ."
544  CALL cp_abort(__location__, &
545  "Two molecules of the same kind have "// &
546  "been created with different numbers of atoms!")
547  END IF
548  END DO
549  END DO
550  END IF
551  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
552  "PRINT%TOPOLOGY_INFO/UTIL_INFO")
553  CALL timestop(handle)
554  END SUBROUTINE topology_generate_molecule
555 
556 ! **************************************************************************************************
557 !> \brief Use info from periodic table and assumptions to generate bonds
558 !> \param topology ...
559 !> \param para_env ...
560 !> \param subsys_section ...
561 !> \author Teodoro Laino 09.2006
562 ! **************************************************************************************************
563  SUBROUTINE topology_generate_bond(topology, para_env, subsys_section)
564  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
565  TYPE(mp_para_env_type), POINTER :: para_env
566  TYPE(section_vals_type), POINTER :: subsys_section
567 
568  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_bond'
569 
570  CHARACTER(LEN=2) :: upper_sym_1
571  INTEGER :: cbond, handle, handle2, i, iatm1, iatm2, iatom, ibond, idim, iw, j, jatom, k, &
572  n_bonds, n_heavy_bonds, n_hydr_bonds, n_rep, natom, npairs, output_unit
573  INTEGER, ALLOCATABLE, DIMENSION(:) :: bond_a, bond_b, list, map_nb
574  INTEGER, DIMENSION(:), POINTER :: isolated_atoms, tmp_v
575  LOGICAL :: connectivity_ok, explicit, print_info
576  LOGICAL, ALLOCATABLE, DIMENSION(:) :: h_list
577  REAL(kind=dp) :: bondparm_factor, cell_v(3), dr(3), &
578  ksign, my_maxrad, r2, r2_min, rbond, &
579  rbond2, tmp
580  REAL(kind=dp), DIMENSION(1, 1) :: r_max, r_minsq
581  REAL(kind=dp), DIMENSION(:), POINTER :: radius
582  REAL(kind=dp), DIMENSION(:, :), POINTER :: pbc_coord
583  TYPE(array2_list_type), DIMENSION(:), POINTER :: bond_list
584  TYPE(atom_info_type), POINTER :: atom_info
585  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
586  TYPE(atomic_kind_type), POINTER :: atomic_kind
587  TYPE(connectivity_info_type), POINTER :: conn_info
588  TYPE(cp_logger_type), POINTER :: logger
589  TYPE(fist_neighbor_type), POINTER :: nonbonded
590  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
591  TYPE(section_vals_type), POINTER :: bond_section, generate_section, &
592  isolated_section
593 
594  NULLIFY (logger, particle_set, atomic_kind_set, nonbonded, bond_section, generate_section)
595  NULLIFY (isolated_atoms, tmp_v)
596  CALL timeset(routinen, handle)
597  logger => cp_get_default_logger()
598  output_unit = cp_logger_get_default_io_unit(logger)
599  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
600  extension=".subsysLog")
601  ! Get atoms that one considers isolated (like ions in solution)
602  ALLOCATE (isolated_atoms(0))
603  generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
604  isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
605  CALL section_vals_get(isolated_section, explicit=explicit)
606  IF (explicit) THEN
607  CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
608  DO i = 1, n_rep
609  CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
610  CALL reallocate(isolated_atoms, 1, SIZE(isolated_atoms) + SIZE(tmp_v))
611  isolated_atoms(SIZE(isolated_atoms) - SIZE(tmp_v) + 1:SIZE(isolated_atoms)) = tmp_v
612  END DO
613  END IF
614  atom_info => topology%atom_info
615  conn_info => topology%conn_info
616  bondparm_factor = topology%bondparm_factor
617  cbond = 0
618  natom = topology%natoms
619  NULLIFY (radius)
620  ! Allocate temporary arrays
621  ALLOCATE (radius(natom))
622  ALLOCATE (list(natom))
623  ALLOCATE (h_list(natom))
624  ALLOCATE (pbc_coord(3, natom))
625  h_list = .false.
626  CALL timeset(trim(routinen)//"_1", handle2)
627  DO iatom = 1, natom
628  list(iatom) = iatom
629  upper_sym_1 = id2str(atom_info%id_element(iatom))
630  IF (topology%bondparm_type == do_bondparm_covalent) THEN
631  CALL get_ptable_info(symbol=upper_sym_1, covalent_radius=radius(iatom))
632  ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN
633  CALL get_ptable_info(symbol=upper_sym_1, vdw_radius=radius(iatom))
634  ELSE
635  cpabort("Illegal bondparm_type")
636  END IF
637  IF (upper_sym_1 == "H ") h_list(iatom) = .true.
638  ! isolated atoms? put the radius to 0.0_dp
639  IF (any(isolated_atoms == iatom)) radius(iatom) = 0.0_dp
640  radius(iatom) = cp_unit_to_cp2k(radius(iatom), "angstrom")
641  IF (iw > 0) WRITE (iw, '(T2,"GENERATE|",5X,A,T50,A5,T60,A,T69,F12.6)') &
642  "In topology_generate_bond :: iatom = ", upper_sym_1, &
643  "radius:", radius(iatom)
644  END DO
645  CALL timestop(handle2)
646  CALL timeset(trim(routinen)//"_2", handle2)
647  ! Initialize fake particle_set and atomic_kinds to generate the bond list
648  ! using the neighboring list routine
649  ALLOCATE (atomic_kind_set(1))
650  CALL allocate_particle_set(particle_set, natom)
651  !
652  my_maxrad = maxval(radius)*2.0_dp
653  atomic_kind => atomic_kind_set(1)
654  CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=1, &
655  name="XXX", element_symbol="XXX", mass=0.0_dp, atom_list=list)
656  CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MAX", r_val=tmp)
657  r_max = tmp
658  IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN
659  IF (output_unit > 0) THEN
660  WRITE (output_unit, '(T2,"GENERATE|",A)') &
661  " ERROR in connectivity generation!", &
662  " The THRESHOLD to select possible bonds is larger than the max. bondlength", &
663  " used to build the neighbors lists. Increase the BONDLENGTH_MAX parameter"
664  WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') &
665  " Present THRESHOLD (", my_maxrad*bondparm_factor, " )."// &
666  " Present BONDLENGTH_MAX (", r_max(1, 1), " )"
667  END IF
668  cpabort("")
669  END IF
670  DO i = 1, natom
671  particle_set(i)%atomic_kind => atomic_kind_set(1)
672  particle_set(i)%r(1) = atom_info%r(1, i)
673  particle_set(i)%r(2) = atom_info%r(2, i)
674  particle_set(i)%r(3) = atom_info%r(3, i)
675  pbc_coord(:, i) = pbc(atom_info%r(:, i), topology%cell)
676  END DO
677  CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MIN", r_val=tmp)
678  r_minsq = tmp*tmp
679  CALL timestop(handle2)
680  CALL timeset(trim(routinen)//"_3", handle2)
681  CALL build_fist_neighbor_lists(atomic_kind_set, particle_set, &
682  cell=topology%cell, r_max=r_max, r_minsq=r_minsq, &
683  ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nonbonded, &
684  para_env=para_env, build_from_scratch=.true., geo_check=.true., &
685  mm_section=generate_section)
686  IF (iw > 0) THEN
687  WRITE (iw, '(T2,"GENERATE| Number of prescreened bonds (neighbors):",T71,I10)') &
688  nonbonded%neighbor_kind_pairs(1)%npairs
689  END IF
690  npairs = 0
691  DO i = 1, SIZE(nonbonded%neighbor_kind_pairs)
692  npairs = npairs + nonbonded%neighbor_kind_pairs(i)%npairs
693  END DO
694  ALLOCATE (bond_a(npairs))
695  ALLOCATE (bond_b(npairs))
696  ALLOCATE (map_nb(npairs))
697  idim = 0
698  DO j = 1, SIZE(nonbonded%neighbor_kind_pairs)
699  DO i = 1, nonbonded%neighbor_kind_pairs(j)%npairs
700  idim = idim + 1
701  bond_a(idim) = nonbonded%neighbor_kind_pairs(j)%list(1, i)
702  bond_b(idim) = nonbonded%neighbor_kind_pairs(j)%list(2, i)
703  map_nb(idim) = j
704  END DO
705  END DO
706  CALL timestop(handle2)
707  CALL timeset(trim(routinen)//"_4", handle2)
708  ! We have a list of neighbors let's order the list w.r.t. the particle number
709  ALLOCATE (bond_list(natom))
710  DO i = 1, natom
711  ALLOCATE (bond_list(i)%array1(0))
712  ALLOCATE (bond_list(i)%array2(0))
713  END DO
714  CALL reorder_structure(bond_list, bond_a, bond_b, map_nb, SIZE(bond_a))
715  DEALLOCATE (bond_a)
716  DEALLOCATE (bond_b)
717  DEALLOCATE (map_nb)
718  ! Find the Real bonds in the system
719  ! Let's start with heavy atoms.. hydrogens will be treated only later on...
720  ! Heavy atoms loop
721  CALL reallocate(conn_info%bond_a, 1, 1)
722  CALL reallocate(conn_info%bond_b, 1, 1)
723  connectivity_ok = .false.
724  ! No need to check consistency between provided molecule name and
725  ! generated connectivity since we overrided the molecule definition.
726  IF (topology%create_molecules) THEN
727  atom_info%id_molname = str2id(s2s("TO_DEFINE_LATER"))
728  ! A real name assignment will then be performed in the reorder module..
729  END IF
730  ! It may happen that the connectivity created is fault for the missing
731  ! of one bond.. this external loop ensures that everything was created
732  ! fits exactly with the definition of molecules..
733  DO WHILE (.NOT. connectivity_ok)
734  n_heavy_bonds = 0
735  n_bonds = 0
736  DO iatm1 = 1, natom
737  IF (h_list(iatm1)) cycle
738  DO j = 1, SIZE(bond_list(iatm1)%array1)
739  iatm2 = bond_list(iatm1)%array1(j)
740  IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) cycle
741  IF (h_list(iatm2) .OR. (iatm2 <= iatm1)) cycle
742  k = bond_list(iatm1)%array2(j)
743  ksign = sign(1.0_dp, real(k, kind=dp))
744  k = abs(k)
745  cell_v = matmul(topology%cell%hmat, &
746  REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, kind=dp))
747  dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
748  r2 = dot_product(dr, dr)
749  IF (r2 <= r_minsq(1, 1)) THEN
750  CALL cp_abort(__location__, &
751  "bond distance between atoms less then the smallest distance provided "// &
752  "in input "//cp_to_string(tmp)//" [bohr]")
753  END IF
754  ! Screen isolated atoms
755  IF ((any(isolated_atoms == iatm1)) .OR. (any(isolated_atoms == iatm2))) cycle
756 
757  ! Screen neighbors
758  IF (topology%bondparm_type == do_bondparm_covalent) THEN
759  rbond = radius(iatm1) + radius(iatm2)
760  ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN
761  rbond = max(radius(iatm1), radius(iatm2))
762  END IF
763  rbond2 = rbond*rbond
764  rbond2 = rbond2*(bondparm_factor)**2
765  !Test the distance to the sum of the covalent radius
766  IF (r2 <= rbond2) THEN
767  n_heavy_bonds = n_heavy_bonds + 1
768  CALL add_bonds_list(conn_info, iatm1, iatm2, n_heavy_bonds)
769  END IF
770  END DO
771  END DO
772  n_hydr_bonds = 0
773  n_bonds = n_heavy_bonds
774  ! Now check bonds formed by hydrogens...
775  ! The hydrogen valence is 1 so we can choose the closest atom..
776  IF (output_unit > 0) WRITE (output_unit, *)
777  DO iatm1 = 1, natom
778  IF (.NOT. h_list(iatm1)) cycle
779  r2_min = huge(0.0_dp)
780  ibond = -1
781  print_info = .true.
782  DO j = 1, SIZE(bond_list(iatm1)%array1)
783  iatm2 = bond_list(iatm1)%array1(j)
784  print_info = .false.
785  IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) cycle
786  IF (h_list(iatm2) .AND. (iatm2 <= iatm1)) cycle
787  ! Screen isolated atoms
788  IF ((any(isolated_atoms == iatm1)) .OR. (any(isolated_atoms == iatm2))) cycle
789 
790  k = bond_list(iatm1)%array2(j)
791  ksign = sign(1.0_dp, real(k, kind=dp))
792  k = abs(k)
793  cell_v = matmul(topology%cell%hmat, &
794  REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, kind=dp))
795  dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
796  r2 = dot_product(dr, dr)
797  IF (r2 <= r_minsq(1, 1)) THEN
798  CALL cp_abort(__location__, &
799  "bond distance between atoms less then the smallest distance provided "// &
800  "in input "//cp_to_string(tmp)//" [bohr]")
801  END IF
802  IF (r2 <= r2_min) THEN
803  r2_min = r2
804  ibond = iatm2
805  END IF
806  END DO
807  IF (ibond == -1) THEN
808  IF (output_unit > 0 .AND. print_info) THEN
809  WRITE (output_unit, '(T2,"GENERATE|",1X,A,I10,A)') &
810  "WARNING:: No connections detected for Hydrogen - Atom Nr:", iatm1, " !"
811  END IF
812  ELSE
813  n_hydr_bonds = n_hydr_bonds + 1
814  n_bonds = n_bonds + 1
815  CALL add_bonds_list(conn_info, min(iatm1, ibond), max(iatm1, ibond), n_bonds)
816  END IF
817  END DO
818  IF (output_unit > 0) THEN
819  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') &
820  " Preliminary Number of Bonds generated:", n_bonds
821  END IF
822  ! External defined bonds (useful for complex connectivity)
823  bond_section => section_vals_get_subs_vals(generate_section, "BOND")
824  CALL connectivity_external_control(section=bond_section, &
825  iarray1=conn_info%bond_a, &
826  iarray2=conn_info%bond_b, &
827  nvar=n_bonds, &
828  topology=topology, &
829  output_unit=output_unit)
830  ! Resize arrays to their proper size..
831  CALL reallocate(conn_info%bond_a, 1, n_bonds)
832  CALL reallocate(conn_info%bond_b, 1, n_bonds)
833  IF (topology%create_molecules) THEN
834  ! Since we create molecule names we're not sure that all atoms are contiguous
835  ! so we need to reorder them on the basis of the generated name
836  IF (.NOT. topology%reorder_atom) THEN
837  topology%reorder_atom = .true.
838  IF (output_unit > 0) WRITE (output_unit, '(T2,"GENERATE|",A)') &
839  " Molecules names have been generated. Now reordering particle set in order to have ", &
840  " atoms belonging to the same molecule in a sequential order."
841  END IF
842  connectivity_ok = .true.
843  ELSE
844  ! Check created connectivity and possibly give the OK to proceed
845  connectivity_ok = check_generate_mol(conn_info%bond_a, conn_info%bond_b, &
846  atom_info, bondparm_factor, output_unit)
847  END IF
848  IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN
849  IF (output_unit > 0) THEN
850  WRITE (output_unit, '(T2,"GENERATE|",A)') &
851  " ERROR in connectivity generation!", &
852  " The THRESHOLD to select possible bonds is bigger than the MAX bondlength", &
853  " used to build the neighbors lists. Increase the BONDLENGTH_MAX patameter"
854  WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') &
855  " Present THRESHOLD (", my_maxrad*bondparm_factor, " )."// &
856  " Present BONDLENGTH_MAX (", r_max(1, 1), " )"
857  END IF
858  cpabort("")
859  END IF
860  END DO
861  IF (connectivity_ok .AND. (output_unit > 0)) THEN
862  WRITE (output_unit, '(T2,"GENERATE|",A)') &
863  " Achieved consistency in connectivity generation."
864  END IF
865  CALL fist_neighbor_deallocate(nonbonded)
866  CALL timestop(handle2)
867  CALL timeset(trim(routinen)//"_6", handle2)
868  ! Deallocate temporary working arrays
869  DO i = 1, natom
870  DEALLOCATE (bond_list(i)%array1)
871  DEALLOCATE (bond_list(i)%array2)
872  END DO
873  DEALLOCATE (bond_list)
874  DEALLOCATE (pbc_coord)
875  DEALLOCATE (radius)
876  DEALLOCATE (list)
877  CALL deallocate_particle_set(particle_set)
878  CALL deallocate_atomic_kind_set(atomic_kind_set)
879  !
880  CALL timestop(handle2)
881  IF (output_unit > 0 .AND. n_bonds > 0) THEN
882  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bonds generated:", &
883  n_bonds
884  END IF
885  CALL timeset(trim(routinen)//"_7", handle2)
886  ! If PARA_RES then activate RESIDUES
887  CALL reallocate(conn_info%c_bond_a, 1, 0)
888  CALL reallocate(conn_info%c_bond_b, 1, 0)
889  IF (topology%para_res) THEN
890  DO ibond = 1, SIZE(conn_info%bond_a)
891  iatom = conn_info%bond_a(ibond)
892  jatom = conn_info%bond_b(ibond)
893  IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
894  (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
895  (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN
896  IF (iw > 0) WRITE (iw, *) " PARA_RES, bond between molecules atom ", &
897  iatom, jatom
898  cbond = cbond + 1
899  CALL reallocate(conn_info%c_bond_a, 1, cbond)
900  CALL reallocate(conn_info%c_bond_b, 1, cbond)
901  conn_info%c_bond_a(cbond) = iatom
902  conn_info%c_bond_b(cbond) = jatom
903  ELSE
904  IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN
905  cpabort("Bonds between different molecule types?")
906  END IF
907  END IF
908  END DO
909  END IF
910  CALL timestop(handle2)
911  DEALLOCATE (isolated_atoms)
912  CALL timestop(handle)
913  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
914  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
915  END SUBROUTINE topology_generate_bond
916 
917 ! **************************************************************************************************
918 !> \brief Performs a check on the generated connectivity
919 !> \param bond_a ...
920 !> \param bond_b ...
921 !> \param atom_info ...
922 !> \param bondparm_factor ...
923 !> \param output_unit ...
924 !> \return ...
925 !> \author Teodoro Laino 09.2006
926 ! **************************************************************************************************
927  FUNCTION check_generate_mol(bond_a, bond_b, atom_info, bondparm_factor, output_unit) &
928  result(conn_ok)
929  INTEGER, DIMENSION(:), POINTER :: bond_a, bond_b
930  TYPE(atom_info_type), POINTER :: atom_info
931  REAL(kind=dp), INTENT(INOUT) :: bondparm_factor
932  INTEGER, INTENT(IN) :: output_unit
933  LOGICAL :: conn_ok
934 
935  CHARACTER(len=*), PARAMETER :: routinen = 'check_generate_mol'
936 
937  CHARACTER(LEN=10) :: ctmp1, ctmp2, ctmp3
938  INTEGER :: handle, i, idim, itype, j, mol_natom, &
939  natom, nsize
940  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mol_info_tmp
941  INTEGER, DIMENSION(:), POINTER :: mol_map, mol_map_o, wrk
942  INTEGER, DIMENSION(:, :), POINTER :: mol_info
943  LOGICAL, DIMENSION(:), POINTER :: icheck
944  TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
945 
946  CALL timeset(routinen, handle)
947  conn_ok = .true.
948  natom = SIZE(atom_info%id_atmname)
949  ALLOCATE (bond_list(natom))
950  DO i = 1, natom
951  ALLOCATE (bond_list(i)%array1(0))
952  END DO
953  CALL reorder_structure(bond_list, bond_a, bond_b, SIZE(bond_a))
954  ALLOCATE (mol_map(natom))
955  ALLOCATE (mol_map_o(natom))
956  ALLOCATE (wrk(natom))
957 
958  DO i = 1, natom
959  mol_map(i) = atom_info%id_molname(i)
960  END DO
961  mol_map_o = mol_map
962 
963  CALL sort(mol_map, natom, wrk)
964  !
965  ! mol(i,1) : stores id of the molecule
966  ! mol(i,2) : stores the total number of atoms forming that kind of molecule
967  ! mol(i,3) : contains the number of molecules generated for that kind
968  ! mol(i,4) : contains the number of atoms forming one molecule of that kind
969  ! Connectivity will be considered correct only if for each i:
970  !
971  ! mol(i,2) = mol(i,3)*mol(i,4)
972  !
973  ! If not, very probably, a bond is missing increase bondparm by 10% and let's
974  ! check if the newest connectivity is bug free..
975  !
976 
977  ALLOCATE (mol_info_tmp(natom, 2))
978 
979  itype = mol_map(1)
980  nsize = 1
981  idim = 1
982  mol_info_tmp(1, 1) = itype
983  DO i = 2, natom
984  IF (mol_map(i) /= itype) THEN
985  nsize = nsize + 1
986  itype = mol_map(i)
987  mol_info_tmp(nsize, 1) = itype
988  mol_info_tmp(nsize - 1, 2) = idim
989  idim = 1
990  ELSE
991  idim = idim + 1
992  END IF
993  END DO
994  mol_info_tmp(nsize, 2) = idim
995 
996  ALLOCATE (mol_info(nsize, 4))
997  mol_info(1:nsize, 1:2) = mol_info_tmp(1:nsize, 1:2)
998  DEALLOCATE (mol_info_tmp)
999 
1000  DO i = 1, nsize
1001  mol_info(i, 3) = 0
1002  mol_info(i, 4) = 0
1003  END DO
1004  !
1005  ALLOCATE (icheck(natom))
1006  icheck = .false.
1007  DO i = 1, natom
1008  IF (icheck(i)) cycle
1009  itype = mol_map_o(i)
1010  mol_natom = 0
1011  CALL give_back_molecule(icheck, bond_list, i, mol_natom, mol_map_o, mol_map_o(i))
1012  DO j = 1, SIZE(mol_info)
1013  IF (itype == mol_info(j, 1)) EXIT
1014  END DO
1015  mol_info(j, 3) = mol_info(j, 3) + 1
1016  IF (mol_info(j, 4) == 0) mol_info(j, 4) = mol_natom
1017  IF (mol_info(j, 4) /= mol_natom) THEN
1018  ! Two same molecules have been found with different number
1019  ! of atoms. This usually indicates a missing bond in the
1020  ! generated connectivity
1021  ! Set connectivity to .false. EXIT and increase bondparm_factor by 1.05
1022  conn_ok = .false.
1023  bondparm_factor = bondparm_factor*1.05_dp
1024  IF (output_unit < 0) EXIT
1025  WRITE (output_unit, '(/,T2,"GENERATE|",A)') " WARNING in connectivity generation!"
1026  WRITE (output_unit, '(T2,"GENERATE|",A)') &
1027  ' Two molecules/residues named ('//trim(id2str(itype))//') have different '// &
1028  ' number of atoms.'
1029  CALL integer_to_string(i, ctmp1)
1030  CALL integer_to_string(mol_natom, ctmp2)
1031  CALL integer_to_string(mol_info(j, 4), ctmp3)
1032  WRITE (output_unit, '(T2,"GENERATE|",A)') ' Molecule starting at position ('// &
1033  trim(ctmp1)//') has Nr. <'//trim(ctmp2)// &
1034  '> of atoms.', ' while the other same molecules have Nr. <'// &
1035  trim(ctmp3)//'> of atoms!'
1036  WRITE (output_unit, '(T2,"GENERATE|",A)') &
1037  ' Increasing bondparm_factor by 1.05.. An error was found in the generated', &
1038  ' connectivity. Retry...'
1039  WRITE (output_unit, '(T2,"GENERATE|",A,F11.6,A,/)') &
1040  " Present value of BONDPARM_FACTOR (", bondparm_factor, " )."
1041  EXIT
1042  END IF
1043  END DO
1044 
1045  DEALLOCATE (icheck)
1046  DEALLOCATE (mol_info)
1047  DEALLOCATE (mol_map)
1048  DEALLOCATE (mol_map_o)
1049  DEALLOCATE (wrk)
1050  DO i = 1, natom
1051  DEALLOCATE (bond_list(i)%array1)
1052  END DO
1053  DEALLOCATE (bond_list)
1054  CALL timestop(handle)
1055  END FUNCTION check_generate_mol
1056 
1057 ! **************************************************************************************************
1058 !> \brief Add/Remove a bond to the generated list
1059 !> Particularly useful for system with complex connectivity
1060 !> \param section ...
1061 !> \param Iarray1 ...
1062 !> \param Iarray2 ...
1063 !> \param Iarray3 ...
1064 !> \param Iarray4 ...
1065 !> \param nvar ...
1066 !> \param topology ...
1067 !> \param output_unit ...
1068 !> \param is_impr ...
1069 !> \author Teodoro Laino 09.2006
1070 ! **************************************************************************************************
1071  SUBROUTINE connectivity_external_control(section, Iarray1, Iarray2, Iarray3, Iarray4, nvar, &
1072  topology, output_unit, is_impr)
1073  TYPE(section_vals_type), POINTER :: section
1074  INTEGER, DIMENSION(:), POINTER :: iarray1, iarray2
1075  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: iarray3, iarray4
1076  INTEGER, INTENT(INOUT) :: nvar
1077  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1078  INTEGER, INTENT(IN) :: output_unit
1079  LOGICAL, INTENT(IN), OPTIONAL :: is_impr
1080 
1081  CHARACTER(LEN=8) :: fmt
1082  INTEGER :: do_action, do_it, i, j, k, n_rep, &
1083  n_rep_val, natom, new_size, nsize
1084  INTEGER, DIMENSION(:), POINTER :: atlist, ilist1, ilist2, ilist3, ilist4
1085  LOGICAL :: explicit, ip3, ip4
1086 
1087  natom = topology%natoms
1088  ! Preliminary sort of arrays
1089  ip3 = PRESENT(iarray3)
1090  ip4 = PRESENT(iarray4)
1091  nsize = 2
1092  IF (ip3) nsize = nsize + 1
1093  IF (ip3 .AND. ip4) nsize = nsize + 1
1094  ! Put the lists always in the canonical order
1095  CALL reorder_list_array(iarray1, iarray2, iarray3, iarray4, nsize, nvar)
1096  ! Go on with external control
1097  CALL section_vals_get(section, explicit=explicit, n_repetition=n_rep)
1098  IF (explicit) THEN
1099  NULLIFY (ilist1, ilist2, ilist3, ilist4, atlist)
1100  ALLOCATE (ilist1(nvar))
1101  ALLOCATE (ilist2(nvar))
1102  ilist1 = iarray1(1:nvar)
1103  ilist2 = iarray2(1:nvar)
1104  SELECT CASE (nsize)
1105  CASE (2) !do nothing
1106  CASE (3)
1107  ALLOCATE (ilist3(nvar))
1108  ilist3 = iarray3(1:nvar)
1109  CASE (4)
1110  ALLOCATE (ilist3(nvar))
1111  ALLOCATE (ilist4(nvar))
1112  ilist3 = iarray3(1:nvar)
1113  ilist4 = iarray4(1:nvar)
1114  CASE DEFAULT
1115  ! Should never reach this point
1116  cpabort("")
1117  END SELECT
1118  CALL list_canonical_order(ilist1, ilist2, ilist3, ilist4, nsize, is_impr)
1119  !
1120  DO i = 1, n_rep
1121  CALL section_vals_val_get(section, "ATOMS", i_rep_section=i, n_rep_val=n_rep_val)
1122  CALL section_vals_val_get(section, "_SECTION_PARAMETERS_", i_rep_section=i, &
1123  i_val=do_action)
1124  !
1125  DO j = 1, n_rep_val
1126  CALL section_vals_val_get(section, "ATOMS", i_rep_section=i, i_rep_val=j, &
1127  i_vals=atlist)
1128  cpassert(SIZE(atlist) == nsize)
1129  CALL integer_to_string(nsize - 1, fmt)
1130  CALL check_element_list(do_it, do_action, atlist, ilist1, ilist2, ilist3, ilist4, &
1131  is_impr)
1132  IF (do_action == do_add) THEN
1133  ! Add to the element to the list
1134  IF (do_it > 0) THEN
1135  nvar = nvar + 1
1136  IF (output_unit > 0) THEN
1137  WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//trim(fmt)//'(A,I6),A,T64,A,I6)') &
1138  "element (", &
1139  atlist(1), (",", atlist(k), k=2, nsize), ") added.", " NEW size::", nvar
1140  END IF
1141  IF (nvar > SIZE(iarray1)) THEN
1142  new_size = int(5 + 1.2*nvar)
1143  CALL reallocate(iarray1, 1, new_size)
1144  CALL reallocate(iarray2, 1, new_size)
1145  SELECT CASE (nsize)
1146  CASE (3)
1147  CALL reallocate(iarray3, 1, new_size)
1148  CASE (4)
1149  CALL reallocate(iarray3, 1, new_size)
1150  CALL reallocate(iarray4, 1, new_size)
1151  END SELECT
1152  END IF
1153  ! Using Ilist instead of atlist the canonical order is preserved..
1154  iarray1(do_it + 1:nvar) = iarray1(do_it:nvar - 1)
1155  iarray2(do_it + 1:nvar) = iarray2(do_it:nvar - 1)
1156  iarray1(do_it) = ilist1(do_it)
1157  iarray2(do_it) = ilist2(do_it)
1158  SELECT CASE (nsize)
1159  CASE (3)
1160  iarray3(do_it + 1:nvar) = iarray3(do_it:nvar - 1)
1161  iarray3(do_it) = ilist3(do_it)
1162  CASE (4)
1163  iarray3(do_it + 1:nvar) = iarray3(do_it:nvar - 1)
1164  iarray4(do_it + 1:nvar) = iarray4(do_it:nvar - 1)
1165  iarray3(do_it) = ilist3(do_it)
1166  iarray4(do_it) = ilist4(do_it)
1167  END SELECT
1168  ELSE
1169  IF (output_unit > 0) THEN
1170  WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//trim(fmt)//'(A,I6),A,T80,A)') &
1171  "element (", &
1172  atlist(1), (",", atlist(k), k=2, nsize), ") already found.", "X"
1173  END IF
1174  END IF
1175  ELSE
1176  ! Remove element from the list
1177  IF (do_it > 0) THEN
1178  nvar = nvar - 1
1179  IF (output_unit > 0) THEN
1180  WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//trim(fmt)//'(A,I6),A,T64,A,I6)') &
1181  "element (", &
1182  atlist(1), (",", atlist(k), k=2, nsize), ") removed.", " NEW size::", nvar
1183  END IF
1184  iarray1(do_it:nvar) = iarray1(do_it + 1:nvar + 1)
1185  iarray2(do_it:nvar) = iarray2(do_it + 1:nvar + 1)
1186  iarray1(nvar + 1) = -huge(0)
1187  iarray2(nvar + 1) = -huge(0)
1188  SELECT CASE (nsize)
1189  CASE (3)
1190  iarray3(do_it:nvar) = iarray3(do_it + 1:nvar + 1)
1191  iarray3(nvar + 1) = -huge(0)
1192  CASE (4)
1193  iarray3(do_it:nvar) = iarray3(do_it + 1:nvar + 1)
1194  iarray4(do_it:nvar) = iarray4(do_it + 1:nvar + 1)
1195  iarray3(nvar + 1) = -huge(0)
1196  iarray4(nvar + 1) = -huge(0)
1197  END SELECT
1198  ELSE
1199  IF (output_unit > 0) THEN
1200  WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//trim(fmt)//'(A,I6),A,T80,A)') &
1201  "element (", &
1202  atlist(1), (",", atlist(k), k=2, nsize), ") not found.", "X"
1203  END IF
1204  END IF
1205  END IF
1206 
1207  END DO
1208  END DO
1209  DEALLOCATE (ilist1)
1210  DEALLOCATE (ilist2)
1211  SELECT CASE (nsize)
1212  CASE (2) ! do nothing
1213  CASE (3)
1214  DEALLOCATE (ilist3)
1215  CASE (4)
1216  DEALLOCATE (ilist3)
1217  DEALLOCATE (ilist4)
1218  CASE DEFAULT
1219  ! Should never reach this point
1220  cpabort("")
1221  END SELECT
1222  END IF
1223  END SUBROUTINE connectivity_external_control
1224 
1225 ! **************************************************************************************************
1226 !> \brief Orders list in the canonical order: the extrema of the list are such
1227 !> that the first extrema is always smaller or equal to the last extrema.
1228 !> \param Ilist1 ...
1229 !> \param Ilist2 ...
1230 !> \param Ilist3 ...
1231 !> \param Ilist4 ...
1232 !> \param nsize ...
1233 !> \param is_impr ...
1234 !> \author Teodoro Laino 09.2006
1235 ! **************************************************************************************************
1236  SUBROUTINE list_canonical_order(Ilist1, Ilist2, Ilist3, Ilist4, nsize, is_impr)
1237  INTEGER, DIMENSION(:), POINTER :: ilist1, ilist2
1238  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ilist3, ilist4
1239  INTEGER, INTENT(IN) :: nsize
1240  LOGICAL, INTENT(IN), OPTIONAL :: is_impr
1241 
1242  INTEGER :: i, ss(3), tmp1, tmp2, tmp3, tt(3)
1243  LOGICAL :: do_impr
1244 
1245  do_impr = .false.
1246  IF (PRESENT(is_impr)) do_impr = is_impr
1247  SELECT CASE (nsize)
1248  CASE (2)
1249  DO i = 1, SIZE(ilist1)
1250  tmp1 = ilist1(i)
1251  tmp2 = ilist2(i)
1252  ilist1(i) = min(tmp1, tmp2)
1253  ilist2(i) = max(tmp1, tmp2)
1254  END DO
1255  CASE (3)
1256  DO i = 1, SIZE(ilist1)
1257  tmp1 = ilist1(i)
1258  tmp2 = ilist3(i)
1259  ilist1(i) = min(tmp1, tmp2)
1260  ilist3(i) = max(tmp1, tmp2)
1261  END DO
1262  CASE (4)
1263  DO i = 1, SIZE(ilist1)
1264  IF (.NOT. do_impr) THEN
1265  tmp1 = ilist1(i)
1266  tmp2 = ilist4(i)
1267  ilist1(i) = min(tmp1, tmp2)
1268  IF (ilist1(i) == tmp2) THEN
1269  tmp3 = ilist3(i)
1270  ilist3(i) = ilist2(i)
1271  ilist2(i) = tmp3
1272  END IF
1273  ilist4(i) = max(tmp1, tmp2)
1274  ELSE
1275  tt(1) = ilist2(i)
1276  tt(2) = ilist3(i)
1277  tt(3) = ilist4(i)
1278  CALL sort(tt, 3, ss)
1279  ilist2(i) = tt(1)
1280  ilist3(i) = tt(2)
1281  ilist4(i) = tt(3)
1282  END IF
1283  END DO
1284  END SELECT
1285 
1286  END SUBROUTINE list_canonical_order
1287 
1288 ! **************************************************************************************************
1289 !> \brief finds an element in the ordered list
1290 !> \param do_it ...
1291 !> \param do_action ...
1292 !> \param atlist ...
1293 !> \param Ilist1 ...
1294 !> \param Ilist2 ...
1295 !> \param Ilist3 ...
1296 !> \param Ilist4 ...
1297 !> \param is_impr ...
1298 !> \author Teodoro Laino 09.2006
1299 ! **************************************************************************************************
1300  SUBROUTINE check_element_list(do_it, do_action, atlist, Ilist1, Ilist2, Ilist3, Ilist4, &
1301  is_impr)
1302  INTEGER, INTENT(OUT) :: do_it
1303  INTEGER, INTENT(IN) :: do_action
1304  INTEGER, DIMENSION(:), POINTER :: atlist, ilist1, ilist2
1305  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ilist3, ilist4
1306  LOGICAL, INTENT(IN), OPTIONAL :: is_impr
1307 
1308  INTEGER :: i, iend, istart, ndim, new_size, nsize, &
1309  ss(3), tmp1, tmp2, tmp3, tt(3)
1310  INTEGER, DIMENSION(4) :: tmp
1311  LOGICAL :: do_impr, found
1312 
1313  do_impr = .false.
1314  IF (PRESENT(is_impr)) do_impr = is_impr
1315  found = .false.
1316  nsize = SIZE(atlist)
1317  ndim = SIZE(ilist1)
1318  DO i = 1, nsize
1319  tmp(i) = atlist(i)
1320  END DO
1321  SELECT CASE (nsize)
1322  CASE (2)
1323  tmp1 = tmp(1)
1324  tmp2 = tmp(2)
1325  tmp(1) = min(tmp1, tmp2)
1326  tmp(2) = max(tmp1, tmp2)
1327  CASE (3)
1328  tmp1 = tmp(1)
1329  tmp2 = tmp(3)
1330  tmp(1) = min(tmp1, tmp2)
1331  tmp(3) = max(tmp1, tmp2)
1332  CASE (4)
1333  IF (.NOT. do_impr) THEN
1334  tmp1 = tmp(1)
1335  tmp2 = tmp(4)
1336  tmp(1) = min(tmp1, tmp2)
1337  IF (tmp(1) == tmp2) THEN
1338  tmp3 = tmp(3)
1339  tmp(3) = tmp(2)
1340  tmp(2) = tmp3
1341  END IF
1342  tmp(4) = max(tmp1, tmp2)
1343  ELSE
1344  tt(1) = tmp(2)
1345  tt(2) = tmp(3)
1346  tt(3) = tmp(4)
1347  CALL sort(tt, 3, ss)
1348  tmp(2) = tt(1)
1349  tmp(3) = tt(2)
1350  tmp(4) = tt(3)
1351  END IF
1352  END SELECT
1353  ! boundary to search
1354  DO istart = 1, ndim
1355  IF (ilist1(istart) >= tmp(1)) EXIT
1356  END DO
1357  ! if nothing there stay within bounds
1358  IF (istart <= ndim) THEN
1359  IF (ilist1(istart) > tmp(1) .AND. (istart /= 1)) istart = istart - 1
1360  END IF
1361  DO iend = istart, ndim
1362  IF (ilist1(iend) /= tmp(1)) EXIT
1363  END DO
1364  IF (iend == ndim + 1) iend = ndim
1365  ! Final search in array
1366  SELECT CASE (nsize)
1367  CASE (2)
1368  DO i = istart, iend
1369  IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2))) EXIT
1370  IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2))) THEN
1371  found = .true.
1372  EXIT
1373  END IF
1374  END DO
1375  CASE (3)
1376  DO i = istart, iend
1377  IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2)) .OR. (ilist3(i) > tmp(3))) EXIT
1378  IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2)) .AND. (ilist3(i) == tmp(3))) THEN
1379  found = .true.
1380  EXIT
1381  END IF
1382  END DO
1383  CASE (4)
1384  DO i = istart, iend
1385  IF ((ilist1(i) > tmp(1)) .OR. (ilist2(i) > tmp(2)) .OR. (ilist3(i) > tmp(3)) .OR. (ilist4(i) > tmp(4))) EXIT
1386  IF ((ilist1(i) == tmp(1)) .AND. (ilist2(i) == tmp(2)) &
1387  .AND. (ilist3(i) == tmp(3)) .AND. (ilist4(i) == tmp(4))) THEN
1388  found = .true.
1389  EXIT
1390  END IF
1391  END DO
1392  END SELECT
1393  SELECT CASE (do_action)
1394  CASE (do_add)
1395  IF (found) THEN
1396  do_it = -i
1397  ! Nothing to modify. Element already present
1398  ! in this case ABS(do_it) gives the exact location of the element
1399  ! in the list
1400  ELSE
1401  ! Let's add the element in the right place of the list.. so that we can keep the
1402  ! canonical order
1403  ! in this case do_it gives the index of the list with indexes bigger than
1404  ! the one we're searching for
1405  ! At the end do_it gives the exact location of the element in the canonical list
1406  do_it = i
1407  new_size = ndim + 1
1408  CALL reallocate(ilist1, 1, new_size)
1409  CALL reallocate(ilist2, 1, new_size)
1410  ilist1(i + 1:new_size) = ilist1(i:ndim)
1411  ilist2(i + 1:new_size) = ilist2(i:ndim)
1412  ilist1(i) = tmp(1)
1413  ilist2(i) = tmp(2)
1414  SELECT CASE (nsize)
1415  CASE (3)
1416  CALL reallocate(ilist3, 1, new_size)
1417  ilist3(i + 1:new_size) = ilist3(i:ndim)
1418  ilist3(i) = tmp(3)
1419  CASE (4)
1420  CALL reallocate(ilist3, 1, new_size)
1421  CALL reallocate(ilist4, 1, new_size)
1422  ilist3(i + 1:new_size) = ilist3(i:ndim)
1423  ilist4(i + 1:new_size) = ilist4(i:ndim)
1424  ilist3(i) = tmp(3)
1425  ilist4(i) = tmp(4)
1426  END SELECT
1427  END IF
1428  CASE (do_remove)
1429  IF (found) THEN
1430  do_it = i
1431  ! Let's delete the element in position do_it
1432  new_size = ndim - 1
1433  ilist1(i:new_size) = ilist1(i + 1:ndim)
1434  ilist2(i:new_size) = ilist2(i + 1:ndim)
1435  CALL reallocate(ilist1, 1, new_size)
1436  CALL reallocate(ilist2, 1, new_size)
1437  SELECT CASE (nsize)
1438  CASE (3)
1439  ilist3(i:new_size) = ilist3(i + 1:ndim)
1440  CALL reallocate(ilist3, 1, new_size)
1441  CASE (4)
1442  ilist3(i:new_size) = ilist3(i + 1:ndim)
1443  ilist4(i:new_size) = ilist4(i + 1:ndim)
1444  CALL reallocate(ilist3, 1, new_size)
1445  CALL reallocate(ilist4, 1, new_size)
1446  END SELECT
1447  ELSE
1448  do_it = -i
1449  ! Nothing to modify. Element not present in the list
1450  ! in this case ABS(do_it) gives the exact location of the element
1451  ! in the list
1452  END IF
1453  END SELECT
1454  END SUBROUTINE check_element_list
1455 
1456 ! **************************************************************************************************
1457 !> \brief Adds a bond to the generated bond list
1458 !> \param conn_info ...
1459 !> \param atm1 ...
1460 !> \param atm2 ...
1461 !> \param n_bonds ...
1462 !> \author Teodoro Laino 09.2006
1463 ! **************************************************************************************************
1464  SUBROUTINE add_bonds_list(conn_info, atm1, atm2, n_bonds)
1465  TYPE(connectivity_info_type), POINTER :: conn_info
1466  INTEGER, INTENT(IN) :: atm1, atm2, n_bonds
1467 
1468  INTEGER :: new_size, old_size
1469 
1470  old_size = SIZE(conn_info%bond_a)
1471  IF (n_bonds > old_size) THEN
1472  new_size = int(5 + 1.2*old_size)
1473  CALL reallocate(conn_info%bond_a, 1, new_size)
1474  CALL reallocate(conn_info%bond_b, 1, new_size)
1475  END IF
1476  conn_info%bond_a(n_bonds) = atm1
1477  conn_info%bond_b(n_bonds) = atm2
1478  END SUBROUTINE add_bonds_list
1479 
1480 ! **************************************************************************************************
1481 !> \brief Using a list of bonds, generate a list of bends
1482 !> \param topology ...
1483 !> \param subsys_section ...
1484 !> \author Teodoro Laino 09.2006
1485 ! **************************************************************************************************
1486  SUBROUTINE topology_generate_bend(topology, subsys_section)
1487  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1488  TYPE(section_vals_type), POINTER :: subsys_section
1489 
1490  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_bend'
1491 
1492  INTEGER :: handle, handle2, i, iw, natom, nbond, &
1493  nsize, ntheta, output_unit
1494  TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
1495  TYPE(connectivity_info_type), POINTER :: conn_info
1496  TYPE(cp_logger_type), POINTER :: logger
1497  TYPE(section_vals_type), POINTER :: bend_section
1498 
1499  NULLIFY (logger)
1500  logger => cp_get_default_logger()
1501  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1502  extension=".subsysLog")
1503  CALL timeset(routinen, handle)
1504  output_unit = cp_logger_get_default_io_unit(logger)
1505  conn_info => topology%conn_info
1506  nbond = 0
1507  ntheta = 0
1508  natom = topology%natoms
1509  ! This call is for connectivity off
1510  IF (ASSOCIATED(conn_info%bond_a)) THEN
1511  nbond = SIZE(conn_info%bond_a)
1512  ELSE
1513  CALL reallocate(conn_info%bond_a, 1, nbond)
1514  CALL reallocate(conn_info%bond_b, 1, nbond)
1515  END IF
1516  IF (nbond /= 0) THEN
1517  nsize = int(5 + 1.2*ntheta)
1518  CALL reallocate(conn_info%theta_a, 1, nsize)
1519  CALL reallocate(conn_info%theta_b, 1, nsize)
1520  CALL reallocate(conn_info%theta_c, 1, nsize)
1521  ! Get list of bonds to pre-process theta
1522  ALLOCATE (bond_list(natom))
1523  DO i = 1, natom
1524  ALLOCATE (bond_list(i)%array1(0))
1525  END DO
1526  CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1527  ! All the dirty job is handled by this routine.. for bends it_levl is equal 3
1528  CALL timeset(routinen//"_1", handle2)
1529  CALL match_iterative_path(iarray1=bond_list, &
1530  iarray2=bond_list, &
1531  max_levl=3, &
1532  nvar=ntheta, &
1533  oarray1=conn_info%theta_a, &
1534  oarray2=conn_info%theta_b, &
1535  oarray3=conn_info%theta_c)
1536  CALL timestop(handle2)
1537  DO i = 1, natom
1538  DEALLOCATE (bond_list(i)%array1)
1539  END DO
1540  DEALLOCATE (bond_list)
1541  IF (output_unit > 0) THEN
1542  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Bends generated:", &
1543  ntheta
1544  END IF
1545  ! External defined bends (useful for complex connectivity)
1546  bend_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%ANGLE")
1547  CALL connectivity_external_control(section=bend_section, &
1548  iarray1=conn_info%theta_a, &
1549  iarray2=conn_info%theta_b, &
1550  iarray3=conn_info%theta_c, &
1551  nvar=ntheta, &
1552  topology=topology, &
1553  output_unit=output_unit)
1554  END IF
1555  ! Resize arrays to their proper size..
1556  CALL reallocate(conn_info%theta_a, 1, ntheta)
1557  CALL reallocate(conn_info%theta_b, 1, ntheta)
1558  CALL reallocate(conn_info%theta_c, 1, ntheta)
1559  IF (output_unit > 0 .AND. ntheta > 0) THEN
1560  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bends generated:", &
1561  ntheta
1562  END IF
1563  CALL timestop(handle)
1564  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1565  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1566  END SUBROUTINE topology_generate_bend
1567 
1568 !
1569 
1570 ! **************************************************************************************************
1571 !> \brief Routine matching iteratively along a graph
1572 !> \param Iarray1 ...
1573 !> \param Iarray2 ...
1574 !> \param Iarray3 ...
1575 !> \param max_levl ...
1576 !> \param Oarray1 ...
1577 !> \param Oarray2 ...
1578 !> \param Oarray3 ...
1579 !> \param Oarray4 ...
1580 !> \param Ilist ...
1581 !> \param it_levl ...
1582 !> \param nvar ...
1583 !> \author Teodoro Laino 09.2006
1584 ! **************************************************************************************************
1585  RECURSIVE SUBROUTINE match_iterative_path(Iarray1, Iarray2, Iarray3, &
1586  max_levl, Oarray1, Oarray2, Oarray3, Oarray4, Ilist, it_levl, nvar)
1587  TYPE(array1_list_type), DIMENSION(:), POINTER :: iarray1
1588  TYPE(array1_list_type), DIMENSION(:), OPTIONAL, &
1589  POINTER :: iarray2, iarray3
1590  INTEGER, INTENT(IN) :: max_levl
1591  INTEGER, DIMENSION(:), POINTER :: oarray1, oarray2
1592  INTEGER, DIMENSION(:), OPTIONAL, POINTER :: oarray3, oarray4
1593  INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: ilist
1594  INTEGER, INTENT(IN), OPTIONAL :: it_levl
1595  INTEGER, INTENT(INOUT) :: nvar
1596 
1597  INTEGER :: i, ind, j, my_levl, natom
1598  INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
1599  LOGICAL :: check
1600  TYPE(array1_list_type), DIMENSION(:), POINTER :: wrk
1601 
1602  check = max_levl >= 2 .AND. max_levl <= 4
1603  cpassert(check)
1604  IF (.NOT. PRESENT(ilist)) THEN
1605  SELECT CASE (max_levl)
1606  CASE (2)
1607  cpassert(.NOT. PRESENT(iarray2))
1608  cpassert(.NOT. PRESENT(iarray3))
1609  cpassert(.NOT. PRESENT(oarray3))
1610  cpassert(.NOT. PRESENT(oarray4))
1611  CASE (3)
1612  cpassert(PRESENT(iarray2))
1613  cpassert(.NOT. PRESENT(iarray3))
1614  cpassert(PRESENT(oarray3))
1615  cpassert(.NOT. PRESENT(oarray4))
1616  CASE (4)
1617  cpassert(PRESENT(iarray2))
1618  cpassert(PRESENT(iarray3))
1619  cpassert(PRESENT(oarray3))
1620  cpassert(PRESENT(oarray4))
1621  END SELECT
1622  END IF
1623  natom = SIZE(iarray1)
1624  IF (.NOT. PRESENT(ilist)) THEN
1625  ! Start a new loop.. Only the first time the routine is called
1626  ALLOCATE (my_list(max_levl))
1627  DO i = 1, natom
1628  my_levl = 1
1629  my_list = -1
1630  my_list(my_levl) = i
1631  CALL match_iterative_path(iarray1=iarray1, &
1632  iarray2=iarray2, &
1633  iarray3=iarray3, &
1634  it_levl=my_levl + 1, &
1635  max_levl=max_levl, &
1636  oarray1=oarray1, &
1637  oarray2=oarray2, &
1638  oarray3=oarray3, &
1639  oarray4=oarray4, &
1640  nvar=nvar, &
1641  ilist=my_list)
1642  END DO
1643  DEALLOCATE (my_list)
1644  ELSE
1645  SELECT CASE (it_levl)
1646  CASE (2)
1647  wrk => iarray1
1648  CASE (3)
1649  wrk => iarray2
1650  CASE (4)
1651  wrk => iarray3
1652  END SELECT
1653  i = ilist(it_levl - 1)
1654  DO j = 1, SIZE(iarray1(i)%array1)
1655  ind = wrk(i)%array1(j)
1656  IF (any(ilist == ind)) cycle
1657  IF (it_levl < max_levl) THEN
1658  ilist(it_levl) = ind
1659  CALL match_iterative_path(iarray1=iarray1, &
1660  iarray2=iarray2, &
1661  iarray3=iarray3, &
1662  it_levl=it_levl + 1, &
1663  max_levl=max_levl, &
1664  oarray1=oarray1, &
1665  oarray2=oarray2, &
1666  oarray3=oarray3, &
1667  oarray4=oarray4, &
1668  nvar=nvar, &
1669  ilist=ilist)
1670  ilist(it_levl) = -1
1671  ELSEIF (it_levl == max_levl) THEN
1672  IF (ilist(1) > ind) cycle
1673  ilist(it_levl) = ind
1674  nvar = nvar + 1
1675  SELECT CASE (it_levl)
1676  CASE (2)
1677  IF (nvar > SIZE(oarray1)) THEN
1678  CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1679  CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1680  END IF
1681  oarray1(nvar) = ilist(1)
1682  oarray2(nvar) = ilist(2)
1683  CASE (3)
1684  IF (nvar > SIZE(oarray1)) THEN
1685  CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1686  CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1687  CALL reallocate(oarray3, 1, int(5 + 1.2*nvar))
1688  END IF
1689  oarray1(nvar) = ilist(1)
1690  oarray2(nvar) = ilist(2)
1691  oarray3(nvar) = ilist(3)
1692  CASE (4)
1693  IF (nvar > SIZE(oarray1)) THEN
1694  CALL reallocate(oarray1, 1, int(5 + 1.2*nvar))
1695  CALL reallocate(oarray2, 1, int(5 + 1.2*nvar))
1696  CALL reallocate(oarray3, 1, int(5 + 1.2*nvar))
1697  CALL reallocate(oarray4, 1, int(5 + 1.2*nvar))
1698  END IF
1699  oarray1(nvar) = ilist(1)
1700  oarray2(nvar) = ilist(2)
1701  oarray3(nvar) = ilist(3)
1702  oarray4(nvar) = ilist(4)
1703  CASE DEFAULT
1704  !should never reach this point
1705  cpabort("")
1706  END SELECT
1707  ilist(it_levl) = -1
1708  ELSE
1709  !should never reach this point
1710  cpabort("")
1711  END IF
1712  END DO
1713  END IF
1714  END SUBROUTINE match_iterative_path
1715 
1716 !
1717 
1718 ! **************************************************************************************************
1719 !> \brief The list of Urey-Bradley is equal to the list of bends
1720 !> \param topology ...
1721 !> \param subsys_section ...
1722 ! **************************************************************************************************
1723  SUBROUTINE topology_generate_ub(topology, subsys_section)
1724  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1725  TYPE(section_vals_type), POINTER :: subsys_section
1726 
1727  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_ub'
1728 
1729  INTEGER :: handle, itheta, iw, ntheta, output_unit
1730  TYPE(connectivity_info_type), POINTER :: conn_info
1731  TYPE(cp_logger_type), POINTER :: logger
1732 
1733  NULLIFY (logger)
1734  logger => cp_get_default_logger()
1735  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1736  extension=".subsysLog")
1737  output_unit = cp_logger_get_default_io_unit(logger)
1738  CALL timeset(routinen, handle)
1739  conn_info => topology%conn_info
1740  ntheta = SIZE(conn_info%theta_a)
1741  CALL reallocate(conn_info%ub_a, 1, ntheta)
1742  CALL reallocate(conn_info%ub_b, 1, ntheta)
1743  CALL reallocate(conn_info%ub_c, 1, ntheta)
1744 
1745  DO itheta = 1, ntheta
1746  conn_info%ub_a(itheta) = conn_info%theta_a(itheta)
1747  conn_info%ub_b(itheta) = conn_info%theta_b(itheta)
1748  conn_info%ub_c(itheta) = conn_info%theta_c(itheta)
1749  END DO
1750  IF (output_unit > 0 .AND. ntheta > 0) THEN
1751  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of UB generated:", &
1752  ntheta
1753  END IF
1754  CALL timestop(handle)
1755  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1756  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1757 
1758  END SUBROUTINE topology_generate_ub
1759 
1760 ! **************************************************************************************************
1761 !> \brief Generate a list of torsions from bonds
1762 !> \param topology ...
1763 !> \param subsys_section ...
1764 !> \author Teodoro Laino 09.2006
1765 ! **************************************************************************************************
1766  SUBROUTINE topology_generate_dihe(topology, subsys_section)
1767  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1768  TYPE(section_vals_type), POINTER :: subsys_section
1769 
1770  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_dihe'
1771 
1772  INTEGER :: handle, i, iw, natom, nbond, nphi, &
1773  nsize, output_unit
1774  TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
1775  TYPE(connectivity_info_type), POINTER :: conn_info
1776  TYPE(cp_logger_type), POINTER :: logger
1777  TYPE(section_vals_type), POINTER :: torsion_section
1778 
1779  NULLIFY (logger)
1780  logger => cp_get_default_logger()
1781  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1782  extension=".subsysLog")
1783  output_unit = cp_logger_get_default_io_unit(logger)
1784  CALL timeset(routinen, handle)
1785  conn_info => topology%conn_info
1786  nphi = 0
1787  nbond = SIZE(conn_info%bond_a)
1788  IF (nbond /= 0) THEN
1789  nsize = int(5 + 1.2*nphi)
1790  CALL reallocate(conn_info%phi_a, 1, nsize)
1791  CALL reallocate(conn_info%phi_b, 1, nsize)
1792  CALL reallocate(conn_info%phi_c, 1, nsize)
1793  CALL reallocate(conn_info%phi_d, 1, nsize)
1794  ! Get list of bonds to pre-process phi
1795  natom = topology%natoms
1796  ALLOCATE (bond_list(natom))
1797  DO i = 1, natom
1798  ALLOCATE (bond_list(i)%array1(0))
1799  END DO
1800  CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1801  ! All the dirty job is handled by this routine.. for torsions it_levl is equal 4
1802  CALL match_iterative_path(iarray1=bond_list, &
1803  iarray2=bond_list, &
1804  iarray3=bond_list, &
1805  max_levl=4, &
1806  nvar=nphi, &
1807  oarray1=conn_info%phi_a, &
1808  oarray2=conn_info%phi_b, &
1809  oarray3=conn_info%phi_c, &
1810  oarray4=conn_info%phi_d)
1811  DO i = 1, natom
1812  DEALLOCATE (bond_list(i)%array1)
1813  END DO
1814  DEALLOCATE (bond_list)
1815  IF (output_unit > 0) THEN
1816  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Torsions generated:", &
1817  nphi
1818  END IF
1819  ! External defined torsions (useful for complex connectivity)
1820  torsion_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%TORSION")
1821  CALL connectivity_external_control(section=torsion_section, &
1822  iarray1=conn_info%phi_a, &
1823  iarray2=conn_info%phi_b, &
1824  iarray3=conn_info%phi_c, &
1825  iarray4=conn_info%phi_d, &
1826  nvar=nphi, &
1827  topology=topology, &
1828  output_unit=output_unit)
1829  END IF
1830  ! Resize arrays to their proper size..
1831  CALL reallocate(conn_info%phi_a, 1, nphi)
1832  CALL reallocate(conn_info%phi_b, 1, nphi)
1833  CALL reallocate(conn_info%phi_c, 1, nphi)
1834  CALL reallocate(conn_info%phi_d, 1, nphi)
1835  IF (output_unit > 0 .AND. nphi > 0) THEN
1836  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Torsions generated:", &
1837  nphi
1838  END IF
1839  CALL timestop(handle)
1840  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1841  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1842 
1843  END SUBROUTINE topology_generate_dihe
1844 
1845 ! **************************************************************************************************
1846 !> \brief Using a list of bends, generate a list of impr
1847 !> \param topology ...
1848 !> \param subsys_section ...
1849 !> \author Teodoro Laino 09.2006
1850 ! **************************************************************************************************
1851  SUBROUTINE topology_generate_impr(topology, subsys_section)
1852  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1853  TYPE(section_vals_type), POINTER :: subsys_section
1854 
1855  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_impr'
1856 
1857  CHARACTER(LEN=2) :: atm_symbol
1858  INTEGER :: handle, i, ind, iw, j, natom, nbond, &
1859  nimpr, nsize, output_unit
1860  LOGICAL :: accept_impr
1861  TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
1862  TYPE(atom_info_type), POINTER :: atom_info
1863  TYPE(connectivity_info_type), POINTER :: conn_info
1864  TYPE(cp_logger_type), POINTER :: logger
1865  TYPE(section_vals_type), POINTER :: impr_section
1866 
1867  NULLIFY (logger)
1868  logger => cp_get_default_logger()
1869  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1870  extension=".subsysLog")
1871  output_unit = cp_logger_get_default_io_unit(logger)
1872  CALL timeset(routinen, handle)
1873  atom_info => topology%atom_info
1874  conn_info => topology%conn_info
1875  natom = topology%natoms
1876  nimpr = 0
1877  nbond = SIZE(conn_info%bond_a)
1878  IF (nbond /= 0) THEN
1879  nsize = int(5 + 1.2*nimpr)
1880  CALL reallocate(conn_info%impr_a, 1, nsize)
1881  CALL reallocate(conn_info%impr_b, 1, nsize)
1882  CALL reallocate(conn_info%impr_c, 1, nsize)
1883  CALL reallocate(conn_info%impr_d, 1, nsize)
1884  ! Get list of bonds to pre-process phi
1885  ALLOCATE (bond_list(natom))
1886  DO i = 1, natom
1887  ALLOCATE (bond_list(i)%array1(0))
1888  END DO
1889  CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1890  DO i = 1, natom
1891  ! Count all atoms with three bonds
1892  IF (SIZE(bond_list(i)%array1) == 3) THEN
1893  ! Problematic cases::
1894  ! Nitrogen
1895  accept_impr = .true.
1896  atm_symbol = id2str(atom_info%id_element(i))
1897  CALL uppercase(atm_symbol)
1898  IF (atm_symbol == "N ") THEN
1899  accept_impr = .false.
1900  ! Impropers on Nitrogen only when there is another atom close to it
1901  ! with other 3 bonds
1902  DO j = 1, 3
1903  ind = bond_list(i)%array1(j)
1904  IF (SIZE(bond_list(ind)%array1) == 3) accept_impr = .true.
1905  END DO
1906  END IF
1907  IF (.NOT. accept_impr) cycle
1908  nimpr = nimpr + 1
1909  IF (nimpr > SIZE(conn_info%impr_a)) THEN
1910  nsize = int(5 + 1.2*nimpr)
1911  CALL reallocate(conn_info%impr_a, 1, nsize)
1912  CALL reallocate(conn_info%impr_b, 1, nsize)
1913  CALL reallocate(conn_info%impr_c, 1, nsize)
1914  CALL reallocate(conn_info%impr_d, 1, nsize)
1915  END IF
1916  conn_info%impr_a(nimpr) = i
1917  conn_info%impr_b(nimpr) = bond_list(i)%array1(1)
1918  conn_info%impr_c(nimpr) = bond_list(i)%array1(2)
1919  conn_info%impr_d(nimpr) = bond_list(i)%array1(3)
1920  END IF
1921  END DO
1922  DO i = 1, natom
1923  DEALLOCATE (bond_list(i)%array1)
1924  END DO
1925  DEALLOCATE (bond_list)
1926  ! External defined impropers (useful for complex connectivity)
1927  impr_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%IMPROPER")
1928  CALL connectivity_external_control(section=impr_section, &
1929  iarray1=conn_info%impr_a, &
1930  iarray2=conn_info%impr_b, &
1931  iarray3=conn_info%impr_c, &
1932  iarray4=conn_info%impr_d, &
1933  nvar=nimpr, &
1934  topology=topology, &
1935  output_unit=output_unit, &
1936  is_impr=.true.)
1937  END IF
1938  ! Resize arrays to their proper size..
1939  CALL reallocate(conn_info%impr_a, 1, nimpr)
1940  CALL reallocate(conn_info%impr_b, 1, nimpr)
1941  CALL reallocate(conn_info%impr_c, 1, nimpr)
1942  CALL reallocate(conn_info%impr_d, 1, nimpr)
1943  IF (output_unit > 0 .AND. nimpr > 0) THEN
1944  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Impropers generated:", &
1945  nimpr
1946  END IF
1947  CALL timestop(handle)
1948  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1949  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1950 
1951  END SUBROUTINE topology_generate_impr
1952 
1953 ! **************************************************************************************************
1954 !> \brief Using a list of torsion, generate a list of onfo
1955 !> \param topology ...
1956 !> \param subsys_section ...
1957 ! **************************************************************************************************
1958  SUBROUTINE topology_generate_onfo(topology, subsys_section)
1959  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1960  TYPE(section_vals_type), POINTER :: subsys_section
1961 
1962  CHARACTER(len=*), PARAMETER :: routinen = 'topology_generate_onfo'
1963 
1964  INTEGER :: atom_a, atom_b, handle, i, ionfo, iw, &
1965  natom, nbond, nphi, ntheta, output_unit
1966  TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list, phi_list, theta_list
1967  TYPE(connectivity_info_type), POINTER :: conn_info
1968  TYPE(cp_logger_type), POINTER :: logger
1969 
1970  NULLIFY (logger)
1971  logger => cp_get_default_logger()
1972  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1973  extension=".subsysLog")
1974  output_unit = cp_logger_get_default_io_unit(logger)
1975  CALL timeset(routinen, handle)
1976 
1977  conn_info => topology%conn_info
1978  natom = topology%natoms
1979 
1980  ! Get list of bonds (sic). Get a list of bonded neighbors for every atom.
1981  ALLOCATE (bond_list(natom))
1982  DO i = 1, natom
1983  ALLOCATE (bond_list(i)%array1(0))
1984  END DO
1985  nbond = SIZE(conn_info%bond_a)
1986  CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1987 
1988  ! Get a list of next nearest neighbors for every atom.
1989  ALLOCATE (theta_list(natom))
1990  DO i = 1, natom
1991  ALLOCATE (theta_list(i)%array1(0))
1992  END DO
1993  ntheta = SIZE(conn_info%theta_a)
1994  CALL reorder_structure(theta_list, conn_info%theta_a, conn_info%theta_c, ntheta)
1995 
1996  ! Get a list of next next nearest neighbors for every atom.
1997  ALLOCATE (phi_list(natom))
1998  DO i = 1, natom
1999  ALLOCATE (phi_list(i)%array1(0))
2000  END DO
2001  nphi = SIZE(conn_info%phi_a)
2002  CALL reorder_structure(phi_list, conn_info%phi_a, conn_info%phi_d, nphi)
2003 
2004  ! Allocate enough (possible too much)
2005  CALL reallocate(conn_info%onfo_a, 1, nphi)
2006  CALL reallocate(conn_info%onfo_b, 1, nphi)
2007 
2008  ionfo = 0
2009  DO atom_a = 1, natom
2010  DO i = 1, SIZE(phi_list(atom_a)%array1)
2011  atom_b = phi_list(atom_a)%array1(i)
2012  ! Avoid trivial duplicates.
2013  IF (atom_a > atom_b) cycle
2014  ! Avoid onfo's in 4-rings.
2015  IF (any(atom_b == bond_list(atom_a)%array1)) cycle
2016  ! Avoid onfo's in 5-rings.
2017  IF (any(atom_b == theta_list(atom_a)%array1)) cycle
2018  ! Avoid onfo's in 6-rings.
2019  IF (any(atom_b == phi_list(atom_a)%array1(:i - 1))) cycle
2020  ionfo = ionfo + 1
2021  conn_info%onfo_a(ionfo) = atom_a
2022  conn_info%onfo_b(ionfo) = atom_b
2023  END DO
2024  END DO
2025 
2026  ! Reallocate such that just enough memory is used.
2027  CALL reallocate(conn_info%onfo_a, 1, ionfo)
2028  CALL reallocate(conn_info%onfo_b, 1, ionfo)
2029 
2030  ! Deallocate bond_list
2031  DO i = 1, natom
2032  DEALLOCATE (bond_list(i)%array1)
2033  END DO
2034  DEALLOCATE (bond_list)
2035  ! Deallocate theta_list
2036  DO i = 1, natom
2037  DEALLOCATE (theta_list(i)%array1)
2038  END DO
2039  DEALLOCATE (theta_list)
2040  ! Deallocate phi_list
2041  DO i = 1, natom
2042  DEALLOCATE (phi_list(i)%array1)
2043  END DO
2044  DEALLOCATE (phi_list)
2045 
2046  ! Final output
2047  IF (output_unit > 0 .AND. ionfo > 0) THEN
2048  WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of 1-4 interactions generated:", &
2049  ionfo
2050  END IF
2051  CALL timestop(handle)
2052  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
2053  "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
2054 
2055  END SUBROUTINE topology_generate_onfo
2056 
2057 END MODULE topology_generate_util
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
subroutine, public set_atomic_kind(atomic_kind, element_symbol, name, mass, kind_number, natom, atom_list, fist_potential, shell, shell_active, damping)
Set the components of an atomic kind data set.
subroutine, public deallocate_atomic_kind_set(atomic_kind_set)
Destructor routine for a set of atomic kinds.
Handles all functions related to the CELL.
Definition: cell_types.F:15
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
unit conversion facility
Definition: cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
Define the neighbor list data types and the corresponding functionality.
subroutine, public fist_neighbor_deallocate(fist_neighbor)
...
Generate the atomic neighbor lists for FIST.
subroutine, public build_fist_neighbor_lists(atomic_kind_set, particle_set, local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, nonbonded, para_env, build_from_scratch, geo_check, mm_section, full_nl, exclusions)
...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_remove
integer, parameter, public do_bondparm_covalent
integer, parameter, public do_conn_off
integer, parameter, public do_bondparm_vdw
integer, parameter, public do_conn_user
integer, parameter, public do_add
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
Utility routines for the memory handling.
Interface to the message passing library MPI.
Define the data structure for the particle information.
subroutine, public deallocate_particle_set(particle_set)
Deallocate a particle set.
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Periodic Table related data definitions.
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.
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
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
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition: string_table.F:115
Utilities for string manipulations.
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Collection of subroutine needed for topology related things.
subroutine, public topology_generate_impr(topology, subsys_section)
Using a list of bends, generate a list of impr.
subroutine, public topology_generate_onfo(topology, subsys_section)
Using a list of torsion, generate a list of onfo.
subroutine, public topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, id_molname)
Generates molnames: useful when the connectivity on file does not provide them.
subroutine, public topology_generate_bend(topology, subsys_section)
Using a list of bonds, generate a list of bends.
subroutine, public topology_generate_molecule(topology, qmmm, qmmm_env, subsys_section)
Use information from bond list to generate molecule. (ie clustering)
subroutine, public topology_generate_dihe(topology, subsys_section)
Generate a list of torsions from bonds.
subroutine, public topology_generate_ub(topology, subsys_section)
The list of Urey-Bradley is equal to the list of bends.
subroutine, public topology_generate_bond(topology, para_env, subsys_section)
Use info from periodic table and assumptions to generate bonds.
Collection of subroutine needed for topology related things.
Definition: topology_util.F:13
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...
recursive subroutine, public give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
...
recursive subroutine, public reorder_list_array(Ilist1, Ilist2, Ilist3, Ilist4, nsize, ndim)
Order arrays of lists.
Control for reading in different topologies and coordinates.
Definition: topology.F:13
All kind of helpful little routines.
Definition: util.F:14