(git:15c1bfc)
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, 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 = trim(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 = trim(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
2057END MODULE topology_generate_util
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:73
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:1852
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:1959
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:1487
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:182
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:1767
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:1724
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:564
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:865
topology_util::give_back_molecule
recursive subroutine, public give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
...
Definition topology_util.F:713
topology_util::reorder_list_array
recursive subroutine, public reorder_list_array(ilist1, ilist2, ilist3, ilist4, nsize, ndim)
Order arrays of lists.
Definition topology_util.F:597
topology
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:50
topology_util::array2_list_type
Definition topology_util.F:55