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