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