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