Line data Source code
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 : ! **************************************************************************************************
13 : MODULE 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 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Generates molnames: useful when the connectivity on file does not
85 : !> provide them
86 : !> \param conn_info ...
87 : !> \param natom ...
88 : !> \param natom_prev ...
89 : !> \param nbond_prev ...
90 : !> \param id_molname ...
91 : !> \author Teodoro Laino [tlaino] - University of Zurich 10.2008
92 : ! **************************************************************************************************
93 22 : SUBROUTINE topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, &
94 22 : 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 22 : 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 78904 : ALLOCATE (atom_bond_list(natom))
110 78860 : DO i = 1, natom
111 78860 : ALLOCATE (atom_bond_list(i)%array1(0))
112 : END DO
113 22 : n = 0
114 22 : 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 114378 : conn_info%bond_b(nbond_prev + 1:) - natom_prev, n)
117 :
118 22 : nmol = 0
119 22 : id_undef = str2id(s2s("__UNDEF__"))
120 78882 : check = ALL(id_molname == id_undef) .OR. ALL(id_molname /= id_undef)
121 22 : CPASSERT(check)
122 78860 : DO i = 1, natom
123 78860 : IF (id_molname(i) == id_undef) THEN
124 21954 : molname = TRIM(basename)//ADJUSTL(cp_to_string(nmol))
125 21954 : CALL generate_molname_low(i, atom_bond_list, molname, id_molname)
126 21954 : nmol = nmol + 1
127 : END IF
128 : END DO
129 78860 : DO i = 1, natom
130 78860 : DEALLOCATE (atom_bond_list(i)%array1)
131 : END DO
132 22 : DEALLOCATE (atom_bond_list)
133 :
134 22 : 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 79132 : 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 79132 : id_molname(i) = str2id(molname)
164 194494 : DO j = 1, SIZE(atom_bond_list(i)%array1)
165 115362 : k = atom_bond_list(i)%array1(j)
166 : IF (debug_this_module) WRITE (*, *) "entering with i:", i, "testing :", k
167 115362 : IF (k == -1) CYCLE
168 57178 : atom_bond_list(i)%array1(j) = -1
169 128560 : WHERE (atom_bond_list(k)%array1 == i) atom_bond_list(k)%array1 = -1
170 194494 : CALL generate_molname_low(k, atom_bond_list, molname, id_molname)
171 : END DO
172 79132 : 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 11404 : 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 11404 : INTEGER, DIMENSION(:), POINTER :: qm_atom_index, wrk1, wrk2
194 : LOGICAL :: do_again, found, my_qmmm
195 11404 : 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 11404 : NULLIFY (logger)
201 22808 : logger => cp_get_default_logger()
202 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
203 11404 : extension=".subsysLog")
204 11404 : CALL timeset(routineN, handle)
205 11404 : NULLIFY (qm_atom_index)
206 11404 : NULLIFY (wrk1)
207 11404 : NULLIFY (wrk2)
208 :
209 11404 : atom_info => topology%atom_info
210 11404 : conn_info => topology%conn_info
211 : !
212 : ! QM/MM coordinate_control
213 : !
214 11404 : my_qmmm = .FALSE.
215 11404 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
216 :
217 11404 : natom = topology%natoms
218 11404 : IF (ASSOCIATED(atom_info%map_mol_typ)) DEALLOCATE (atom_info%map_mol_typ)
219 34212 : ALLOCATE (atom_info%map_mol_typ(natom))
220 :
221 11404 : IF (ASSOCIATED(atom_info%map_mol_num)) DEALLOCATE (atom_info%map_mol_num)
222 22808 : ALLOCATE (atom_info%map_mol_num(natom))
223 :
224 11404 : IF (ASSOCIATED(atom_info%map_mol_res)) DEALLOCATE (atom_info%map_mol_res)
225 22808 : ALLOCATE (atom_info%map_mol_res(natom))
226 :
227 : ! Initialisation
228 776609 : atom_info%map_mol_typ(:) = 0
229 776609 : atom_info%map_mol_num(:) = -1
230 776609 : atom_info%map_mol_res(:) = 1
231 :
232 : ! Parse the atom list to find the different molecule types and residues
233 11404 : ntype = 1
234 11404 : atom_info%map_mol_typ(1) = 1
235 11404 : resid = 1
236 11404 : CALL reallocate(wrk1, 1, nblock)
237 11404 : wrk1(1) = atom_info%id_molname(1)
238 765205 : DO iatom = 2, natom
239 765205 : IF (topology%conn_type == do_conn_off) THEN
240 : ! No connectivity: each atom becomes a molecule of its own molecule kind
241 44894 : ntype = ntype + 1
242 44894 : atom_info%map_mol_typ(iatom) = ntype
243 708907 : 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 29570 : IF ((atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1)) .AND. &
247 : (.NOT. MODULO(iatom, topology%natom_muc) == 1)) THEN
248 28186 : atom_info%map_mol_typ(iatom) = atom_info%map_mol_typ(iatom - 1)
249 28186 : IF (atom_info%id_resname(iatom) == atom_info%id_resname(iatom - 1)) THEN
250 26940 : atom_info%map_mol_res(iatom) = atom_info%map_mol_res(iatom - 1)
251 : ELSE
252 1246 : resid = resid + 1
253 1246 : atom_info%map_mol_res(iatom) = resid
254 : END IF
255 : ELSE
256 : ! Check if the type is already known
257 1384 : found = .FALSE.
258 19630 : DO itype = 1, ntype
259 19630 : IF (atom_info%id_molname(iatom) == wrk1(itype)) THEN
260 998 : atom_info%map_mol_typ(iatom) = itype
261 : found = .TRUE.
262 : EXIT
263 : END IF
264 : END DO
265 : IF (.NOT. found) THEN
266 386 : ntype = ntype + 1
267 386 : atom_info%map_mol_typ(iatom) = ntype
268 386 : IF (ntype > SIZE(wrk1)) CALL reallocate(wrk1, 1, 2*SIZE(wrk1))
269 386 : wrk1(ntype) = atom_info%id_molname(iatom)
270 : END IF
271 1384 : resid = resid + 1
272 1384 : atom_info%map_mol_res(iatom) = resid
273 : END IF
274 : ELSE
275 679337 : IF (atom_info%id_molname(iatom - 1) == atom_info%id_molname(iatom)) THEN
276 593316 : atom_info%map_mol_typ(iatom) = ntype
277 : ELSE
278 86021 : ntype = ntype + 1
279 86021 : atom_info%map_mol_typ(iatom) = ntype
280 : END IF
281 : END IF
282 : END DO
283 11404 : DEALLOCATE (wrk1)
284 :
285 11404 : 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 799417 : ALLOCATE (atom_bond_list(natom))
290 776609 : DO I = 1, natom
291 776609 : ALLOCATE (atom_bond_list(I)%array1(0))
292 : END DO
293 11404 : N = 0
294 11404 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
295 11404 : CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
296 11404 : CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
297 776609 : DO I = 1, natom
298 776609 : DEALLOCATE (atom_bond_list(I)%array1)
299 : END DO
300 11404 : DEALLOCATE (atom_bond_list)
301 11404 : IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of molecule generation"
302 :
303 : ! Modify according map_mol_typ the array map_mol_num
304 11404 : IF (iw > 0) WRITE (iw, '(/,T2,A)') "Checking for non-continuous generated molecules"
305 : ! Check molecule number
306 22808 : ALLOCATE (wrk1(natom))
307 22808 : ALLOCATE (wrk2(natom))
308 1553218 : 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 11404 : CALL sort(wrk1, natom, wrk2)
317 11404 : istart = 1
318 11404 : mol_typ = wrk1(istart)
319 765205 : DO i = 2, natom
320 765205 : IF (mol_typ /= wrk1(i)) THEN
321 327790 : iend = i - 1
322 1063655 : first = MINVAL(wrk2(istart:iend))
323 1063655 : last = MAXVAL(wrk2(istart:iend))
324 327790 : nlocl = last - first + 1
325 327790 : 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 0 : "It may help to fix this issue.")
334 : END IF
335 327790 : istart = i
336 327790 : mol_typ = wrk1(istart)
337 : END IF
338 : END DO
339 11404 : iend = i - 1
340 40744 : first = MINVAL(wrk2(istart:iend))
341 40744 : last = MAXVAL(wrk2(istart:iend))
342 11404 : nlocl = last - first + 1
343 11404 : 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 0 : "It may help to fix this issue.")
352 : END IF
353 11404 : DEALLOCATE (wrk1)
354 11404 : DEALLOCATE (wrk2)
355 11404 : IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of check"
356 :
357 11404 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "Start of renumbering molecules"
358 11404 : IF (topology%conn_type == do_conn_user) THEN
359 170 : mol_num = 1
360 170 : atom_info%map_mol_num(1) = 1
361 29740 : DO iatom = 2, natom
362 29570 : IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN
363 : mol_num = 1
364 29184 : ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN
365 2244 : mol_num = mol_num + 1
366 : END IF
367 29740 : atom_info%map_mol_num(iatom) = mol_num
368 : END DO
369 : ELSE
370 11234 : mol_typ = atom_info%map_mol_typ(1)
371 11234 : mol_num = atom_info%map_mol_num(1)
372 735465 : DO i = 2, natom
373 724231 : IF (atom_info%map_mol_typ(i) /= mol_typ) THEN
374 130915 : myind = atom_info%map_mol_num(i) - mol_num + 1
375 130915 : CPASSERT(myind /= atom_info%map_mol_num(i - 1))
376 130915 : mol_typ = atom_info%map_mol_typ(i)
377 130915 : mol_num = atom_info%map_mol_num(i)
378 : END IF
379 735465 : atom_info%map_mol_num(i) = atom_info%map_mol_num(i) - mol_num + 1
380 : END DO
381 : END IF
382 11404 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "End of renumbering molecules"
383 :
384 : ! Optionally, use the residues as molecules
385 11404 : CALL timeset(routineN//"_PARA_RES", handle2)
386 11404 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A,L2)") "Starting PARA_RES: ", topology%para_res
387 11404 : IF (topology%para_res) THEN
388 10672 : IF (topology%conn_type == do_conn_user) THEN
389 34 : atom_info%id_molname(:) = atom_info%id_resname(:)
390 6 : ntype = 1
391 6 : atom_info%map_mol_typ(1) = 1
392 6 : mol_num = 1
393 6 : atom_info%map_mol_num(1) = 1
394 28 : DO iatom = 2, natom
395 22 : IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN
396 6 : ntype = ntype + 1
397 6 : mol_num = 1
398 16 : ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN
399 0 : mol_num = mol_num + 1
400 : END IF
401 22 : atom_info%map_mol_typ(iatom) = ntype
402 28 : atom_info%map_mol_num(iatom) = mol_num
403 : END DO
404 : ELSE
405 10666 : mol_res = 1
406 10666 : mol_typ = atom_info%map_mol_typ(1)
407 10666 : mol_num = atom_info%map_mol_num(1)
408 10666 : atom_info%map_mol_res(1) = mol_res
409 731469 : DO i = 2, natom
410 720803 : 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 226039 : mol_res = mol_res + 1
413 : END IF
414 720803 : IF ((atom_info%map_mol_typ(i) /= mol_typ) .OR. &
415 : (atom_info%map_mol_num(i) /= mol_num)) THEN
416 294920 : mol_typ = atom_info%map_mol_typ(i)
417 294920 : mol_num = atom_info%map_mol_num(i)
418 294920 : mol_res = 1
419 : END IF
420 731469 : atom_info%map_mol_res(i) = mol_res
421 : END DO
422 : END IF
423 : END IF
424 11404 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "End of PARA_RES"
425 11404 : CALL timestop(handle2)
426 :
427 11404 : IF (iw > 0) THEN
428 1559 : DO iatom = 1, natom
429 1532 : WRITE (iw, '(4(1X,A,":",I0),2(1X,A,1X,A))') "iatom", iatom, &
430 1532 : "map_mol_typ", atom_info%map_mol_typ(iatom), &
431 1532 : "map_mol_num", atom_info%map_mol_num(iatom), &
432 1532 : "map_mol_res", atom_info%map_mol_res(iatom), &
433 1532 : "mol_name:", TRIM(id2str(atom_info%id_molname(iatom))), &
434 3091 : "res_name:", TRIM(id2str(atom_info%id_resname(iatom)))
435 : END DO
436 : END IF
437 :
438 11404 : IF (my_qmmm) THEN
439 394 : do_again = .FALSE.
440 394 : IF (iw > 0) WRITE (iw, *) "MAP_MOL_NUM ", atom_info%map_mol_num
441 394 : IF (iw > 0) WRITE (iw, *) "MAP_MOL_TYP ", atom_info%map_mol_typ
442 394 : IF (iw > 0) WRITE (iw, *) "MAP_MOL_RES ", atom_info%map_mol_res
443 1182 : ALLOCATE (qm_atom_index(SIZE(qmmm_env%qm_atom_index)))
444 6572 : qm_atom_index = qmmm_env%qm_atom_index
445 3286 : CPASSERT(ALL(qm_atom_index /= 0))
446 1958 : DO myind = 1, SIZE(qm_atom_index)
447 1854 : IF (qm_atom_index(myind) == 0) CYCLE
448 : CALL find_boundary(atom_info%map_mol_typ, natom, ifirst, ilast, &
449 974 : 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 974 : atom_info%map_mol_typ(qm_atom_index(myind)), atom_info%map_mol_num(qm_atom_index(myind)))
452 974 : IF (iw > 0) WRITE (iw, *) "qm fragment:: ifirst, ilast", ifirst, ilast
453 974 : CPASSERT(((ifirst /= 0) .OR. (ilast /= natom)))
454 16314 : DO iatm = ifirst, ilast
455 : atom_info%id_molname(iatm) = str2id(s2s("_QM_"// &
456 15340 : TRIM(id2str(atom_info%id_molname(iatm)))))
457 15340 : IF (iw > 0) WRITE (iw, *) "QM Molecule name :: ", id2str(atom_info%id_molname(iatm))
458 786952 : WHERE (qm_atom_index == iatm) qm_atom_index = 0
459 : END DO
460 466890 : DO iatm = 1, ifirst - 1
461 59902382 : IF (ANY(qm_atom_index == iatm)) do_again = .TRUE.
462 : END DO
463 626780 : DO iatm = ilast + 1, natom
464 62275092 : IF (ANY(qm_atom_index == iatm)) do_again = .TRUE.
465 : END DO
466 974 : IF (iw > 0) WRITE (iw, *) " Another QM fragment? :: ", do_again
467 974 : IF (ifirst /= 1) THEN
468 656 : jump1 = atom_info%map_mol_typ(ifirst) - atom_info%map_mol_typ(ifirst - 1)
469 656 : CPASSERT(jump1 <= 1 .AND. jump1 >= 0)
470 656 : jump1 = ABS(jump1 - 1)
471 : ELSE
472 : jump1 = 0
473 : END IF
474 974 : IF (ilast /= natom) THEN
475 878 : jump2 = atom_info%map_mol_typ(ilast + 1) - atom_info%map_mol_typ(ilast)
476 878 : CPASSERT(jump2 <= 1 .AND. jump2 >= 0)
477 878 : jump2 = ABS(jump2 - 1)
478 : ELSE
479 : jump2 = 0
480 : END IF
481 :
482 : ! Changing mol_type consistently
483 642120 : DO iatm = ifirst, natom
484 642120 : atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump1
485 : END DO
486 626780 : DO iatm = ilast + 1, natom
487 626780 : atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump2
488 : END DO
489 974 : IF (jump1 == 1) THEN
490 608 : DO iatm = ifirst, ilast
491 608 : atom_info%map_mol_num(iatm) = 1
492 : END DO
493 : END IF
494 :
495 974 : IF (jump2 == 1) THEN
496 254 : 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 254 : atom_info%map_mol_typ(ilast + 1), atom_info%map_mol_num(ilast + 1))
499 254 : atom_in_mol = ilast - ifirst + 1
500 254 : inum = 1
501 254 : DO iatm = first, last, atom_in_mol
502 167580 : atom_info%map_mol_num(iatm:iatm + atom_in_mol - 1) = inum
503 42224 : inum = inum + 1
504 : END DO
505 : END IF
506 :
507 2052 : IF (.NOT. do_again) EXIT
508 : END DO
509 394 : DEALLOCATE (qm_atom_index)
510 :
511 394 : IF (iw > 0) THEN
512 0 : WRITE (iw, *) "After the QM/MM Setup:"
513 0 : DO iatom = 1, natom
514 0 : WRITE (iw, *) " iatom,map_mol_typ,map_mol_num ", iatom, &
515 0 : 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 11404 : IF (iw > 0) THEN
523 27 : WRITE (iw, *) "SUMMARY:: Number of molecule kinds found:", ntype
524 1559 : ntype = MAXVAL(atom_info%map_mol_typ)
525 460 : DO i = 1, ntype
526 154989 : atom_in_kind = COUNT(atom_info%map_mol_typ == i)
527 433 : WRITE (iw, *) "Molecule kind:", i, " contains", atom_in_kind, " atoms"
528 433 : IF (atom_in_kind <= 1) CYCLE
529 24 : CALL find_boundary(atom_info%map_mol_typ, natom, first, last, i)
530 24 : WRITE (iw, *) "Boundary atoms:", first, last
531 24 : CPASSERT(last - first + 1 == atom_in_kind)
532 1147 : max_mol_num = MAXVAL(atom_info%map_mol_num(first:last))
533 24 : WRITE (iw, *) "Number of molecules of kind", i, "is ::", max_mol_num
534 24 : atom_in_mol = atom_in_kind/max_mol_num
535 24 : WRITE (iw, *) "Number of atoms per each molecule:", atom_in_mol
536 1147 : WRITE (iw, *) "MAP_MOL_TYP::", atom_info%map_mol_typ(first:last)
537 1147 : WRITE (iw, *) "MAP_MOL_NUM::", atom_info%map_mol_num(first:last)
538 1147 : WRITE (iw, *) "MAP_MOL_RES::", atom_info%map_mol_res(first:last)
539 : !
540 195 : DO j = 1, max_mol_num
541 31706 : IF (COUNT(atom_info%map_mol_num(first:last) == j) /= atom_in_mol) THEN
542 0 : WRITE (iw, *) "molecule type:", i, "molecule num:", j, " has ", &
543 0 : COUNT(atom_info%map_mol_num(first:last) == j), &
544 0 : " atoms instead of ", atom_in_mol, " ."
545 : CALL cp_abort(__LOCATION__, &
546 : "Two molecules of the same kind have "// &
547 0 : "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 11404 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
554 11404 : CALL timestop(handle)
555 45616 : 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 9197 : 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 9197 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bond_a, bond_b, list, map_nb
575 9197 : INTEGER, DIMENSION(:), POINTER :: isolated_atoms, tmp_v
576 : LOGICAL :: connectivity_ok, explicit, print_info
577 9197 : 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 9197 : REAL(KIND=dp), DIMENSION(:), POINTER :: radius
583 9197 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pbc_coord
584 9197 : TYPE(array2_list_type), DIMENSION(:), POINTER :: bond_list
585 : TYPE(atom_info_type), POINTER :: atom_info
586 9197 : 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 9197 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
592 : TYPE(section_vals_type), POINTER :: bond_section, generate_section, &
593 : isolated_section
594 :
595 9197 : NULLIFY (logger, particle_set, atomic_kind_set, nonbonded, bond_section, generate_section)
596 9197 : NULLIFY (isolated_atoms, tmp_v)
597 9197 : CALL timeset(routineN, handle)
598 9197 : logger => cp_get_default_logger()
599 9197 : 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 9197 : extension=".subsysLog")
602 : ! Get atoms that one considers isolated (like ions in solution)
603 9197 : ALLOCATE (isolated_atoms(0))
604 9197 : generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
605 9197 : isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
606 9197 : CALL section_vals_get(isolated_section, explicit=explicit)
607 9197 : IF (explicit) THEN
608 8 : CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
609 20 : DO i = 1, n_rep
610 12 : CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
611 12 : CALL reallocate(isolated_atoms, 1, SIZE(isolated_atoms) + SIZE(tmp_v))
612 196 : isolated_atoms(SIZE(isolated_atoms) - SIZE(tmp_v) + 1:SIZE(isolated_atoms)) = tmp_v
613 : END DO
614 : END IF
615 9197 : atom_info => topology%atom_info
616 9197 : conn_info => topology%conn_info
617 9197 : bondparm_factor = topology%bondparm_factor
618 9197 : cbond = 0
619 9197 : natom = topology%natoms
620 9197 : NULLIFY (radius)
621 : ! Allocate temporary arrays
622 27591 : ALLOCATE (radius(natom))
623 27591 : ALLOCATE (list(natom))
624 18394 : ALLOCATE (h_list(natom))
625 27591 : ALLOCATE (pbc_coord(3, natom))
626 9197 : h_list = .FALSE.
627 9197 : CALL timeset(TRIM(routineN)//"_1", handle2)
628 236993 : DO iatom = 1, natom
629 227796 : list(iatom) = iatom
630 227796 : upper_sym_1 = TRIM(id2str(atom_info%id_element(iatom)))
631 227796 : IF (topology%bondparm_type == do_bondparm_covalent) THEN
632 227796 : CALL get_ptable_info(symbol=upper_sym_1, covalent_radius=radius(iatom))
633 0 : ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN
634 0 : CALL get_ptable_info(symbol=upper_sym_1, vdw_radius=radius(iatom))
635 : ELSE
636 0 : CPABORT("Illegal bondparm_type")
637 : END IF
638 227796 : IF (upper_sym_1 == "H ") h_list(iatom) = .TRUE.
639 : ! isolated atoms? put the radius to 0.0_dp
640 353672 : IF (ANY(isolated_atoms == iatom)) radius(iatom) = 0.0_dp
641 227796 : radius(iatom) = cp_unit_to_cp2k(radius(iatom), "angstrom")
642 227796 : IF (iw > 0) WRITE (iw, '(T2,"GENERATE|",5X,A,T50,A5,T60,A,T69,F12.6)') &
643 3186 : "In topology_generate_bond :: iatom = ", upper_sym_1, &
644 15569 : "radius:", radius(iatom)
645 : END DO
646 9197 : CALL timestop(handle2)
647 9197 : 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 18394 : ALLOCATE (atomic_kind_set(1))
651 9197 : CALL allocate_particle_set(particle_set, natom)
652 : !
653 236993 : my_maxrad = MAXVAL(radius)*2.0_dp
654 9197 : atomic_kind => atomic_kind_set(1)
655 : CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=1, &
656 9197 : name="XXX", element_symbol="XXX", mass=0.0_dp, atom_list=list)
657 9197 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MAX", r_val=tmp)
658 27591 : r_max = tmp
659 9197 : IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN
660 0 : IF (output_unit > 0) THEN
661 : WRITE (output_unit, '(T2,"GENERATE|",A)') &
662 0 : " ERROR in connectivity generation!", &
663 0 : " The THRESHOLD to select possible bonds is larger than the max. bondlength", &
664 0 : " used to build the neighbors lists. Increase the BONDLENGTH_MAX parameter"
665 : WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') &
666 0 : " Present THRESHOLD (", my_maxrad*bondparm_factor, " )."// &
667 0 : " Present BONDLENGTH_MAX (", r_max(1, 1), " )"
668 : END IF
669 0 : CPABORT("Unable to generate connectivity")
670 : END IF
671 236993 : DO i = 1, natom
672 227796 : particle_set(i)%atomic_kind => atomic_kind_set(1)
673 227796 : particle_set(i)%r(1) = atom_info%r(1, i)
674 227796 : particle_set(i)%r(2) = atom_info%r(2, i)
675 227796 : particle_set(i)%r(3) = atom_info%r(3, i)
676 920381 : pbc_coord(:, i) = pbc(atom_info%r(:, i), topology%cell)
677 : END DO
678 9197 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MIN", r_val=tmp)
679 27591 : r_minsq = tmp*tmp
680 9197 : CALL timestop(handle2)
681 9197 : 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 9197 : mm_section=generate_section)
687 9197 : IF (iw > 0) THEN
688 : WRITE (iw, '(T2,"GENERATE| Number of prescreened bonds (neighbors):",T71,I10)') &
689 8 : nonbonded%neighbor_kind_pairs(1)%npairs
690 : END IF
691 9197 : npairs = 0
692 170668 : DO i = 1, SIZE(nonbonded%neighbor_kind_pairs)
693 170668 : npairs = npairs + nonbonded%neighbor_kind_pairs(i)%npairs
694 : END DO
695 26757 : ALLOCATE (bond_a(npairs))
696 17560 : ALLOCATE (bond_b(npairs))
697 17560 : ALLOCATE (map_nb(npairs))
698 9197 : idim = 0
699 170668 : DO j = 1, SIZE(nonbonded%neighbor_kind_pairs)
700 1345077 : DO i = 1, nonbonded%neighbor_kind_pairs(j)%npairs
701 1174409 : idim = idim + 1
702 1174409 : bond_a(idim) = nonbonded%neighbor_kind_pairs(j)%list(1, i)
703 1174409 : bond_b(idim) = nonbonded%neighbor_kind_pairs(j)%list(2, i)
704 1335880 : map_nb(idim) = j
705 : END DO
706 : END DO
707 9197 : CALL timestop(handle2)
708 9197 : 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 255387 : ALLOCATE (bond_list(natom))
711 236993 : DO I = 1, natom
712 227796 : ALLOCATE (bond_list(I)%array1(0))
713 236993 : ALLOCATE (bond_list(I)%array2(0))
714 : END DO
715 9197 : CALL reorder_structure(bond_list, bond_a, bond_b, map_nb, SIZE(bond_a))
716 9197 : DEALLOCATE (bond_a)
717 9197 : DEALLOCATE (bond_b)
718 9197 : 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 9197 : CALL reallocate(conn_info%bond_a, 1, 1)
723 9197 : CALL reallocate(conn_info%bond_b, 1, 1)
724 9197 : connectivity_ok = .FALSE.
725 : ! No need to check consistency between provided molecule name and
726 : ! generated connectivity since we overrided the molecule definition.
727 9197 : IF (topology%create_molecules) THEN
728 9436 : 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 18396 : DO WHILE (.NOT. connectivity_ok)
735 9199 : n_heavy_bonds = 0
736 : n_bonds = 0
737 238225 : DO iatm1 = 1, natom
738 229026 : IF (h_list(iatm1)) CYCLE
739 1177381 : DO j = 1, SIZE(bond_list(iatm1)%array1)
740 1053768 : iatm2 = bond_list(iatm1)%array1(j)
741 1053768 : IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) CYCLE
742 683562 : IF (h_list(iatm2) .OR. (iatm2 <= iatm1)) CYCLE
743 139888 : k = bond_list(iatm1)%array2(j)
744 139888 : ksign = SIGN(1.0_dp, REAL(k, KIND=dp))
745 139888 : k = ABS(k)
746 : cell_v = MATMUL(topology%cell%hmat, &
747 2238208 : REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, KIND=dp))
748 559552 : dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
749 559552 : r2 = DOT_PRODUCT(dr, dr)
750 139888 : IF (r2 <= r_minsq(1, 1)) THEN
751 : CALL cp_abort(__LOCATION__, &
752 : "bond distance between atoms less then the smallest distance provided "// &
753 0 : "in input "//cp_to_string(tmp)//" [bohr]")
754 : END IF
755 : ! Screen isolated atoms
756 1617502 : IF ((ANY(isolated_atoms == iatm1)) .OR. (ANY(isolated_atoms == iatm2))) CYCLE
757 :
758 : ! Screen neighbors
759 138424 : IF (topology%bondparm_type == do_bondparm_covalent) THEN
760 138424 : rbond = radius(iatm1) + radius(iatm2)
761 0 : ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN
762 0 : rbond = MAX(radius(iatm1), radius(iatm2))
763 : END IF
764 138424 : rbond2 = rbond*rbond
765 138424 : rbond2 = rbond2*(bondparm_factor)**2
766 : !Test the distance to the sum of the covalent radius
767 367450 : IF (r2 <= rbond2) THEN
768 18044 : n_heavy_bonds = n_heavy_bonds + 1
769 18044 : CALL add_bonds_list(conn_info, iatm1, iatm2, n_heavy_bonds)
770 : END IF
771 : END DO
772 : END DO
773 9199 : n_hydr_bonds = 0
774 9199 : 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 9199 : IF (output_unit > 0) WRITE (output_unit, *)
778 238225 : DO iatm1 = 1, natom
779 229026 : IF (.NOT. h_list(iatm1)) CYCLE
780 114612 : r2_min = HUGE(0.0_dp)
781 114612 : ibond = -1
782 114612 : print_info = .TRUE.
783 1416226 : DO j = 1, SIZE(bond_list(iatm1)%array1)
784 1301614 : iatm2 = bond_list(iatm1)%array1(j)
785 1301614 : print_info = .FALSE.
786 1301614 : IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) CYCLE
787 1202826 : IF (h_list(iatm2) .AND. (iatm2 <= iatm1)) CYCLE
788 : ! Screen isolated atoms
789 12228642 : IF ((ANY(isolated_atoms == iatm1)) .OR. (ANY(isolated_atoms == iatm2))) CYCLE
790 :
791 799594 : k = bond_list(iatm1)%array2(j)
792 799594 : ksign = SIGN(1.0_dp, REAL(k, KIND=dp))
793 799594 : k = ABS(k)
794 : cell_v = MATMUL(topology%cell%hmat, &
795 12793504 : REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, KIND=dp))
796 3198376 : dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v
797 3198376 : r2 = DOT_PRODUCT(dr, dr)
798 799594 : IF (r2 <= r_minsq(1, 1)) THEN
799 : CALL cp_abort(__LOCATION__, &
800 : "bond distance between atoms less then the smallest distance provided "// &
801 0 : "in input "//cp_to_string(tmp)//" [bohr]")
802 : END IF
803 914206 : IF (r2 <= r2_min) THEN
804 227442 : r2_min = r2
805 227442 : ibond = iatm2
806 : END IF
807 : END DO
808 123811 : IF (ibond == -1) THEN
809 15570 : IF (output_unit > 0 .AND. print_info) THEN
810 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,I10,A)') &
811 135 : "WARNING:: No connections detected for Hydrogen - Atom Nr:", iatm1, " !"
812 : END IF
813 : ELSE
814 99042 : n_hydr_bonds = n_hydr_bonds + 1
815 99042 : n_bonds = n_bonds + 1
816 99042 : CALL add_bonds_list(conn_info, MIN(iatm1, ibond), MAX(iatm1, ibond), n_bonds)
817 : END IF
818 : END DO
819 9199 : IF (output_unit > 0) THEN
820 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') &
821 4701 : " Preliminary Number of Bonds generated:", n_bonds
822 : END IF
823 : ! External defined bonds (useful for complex connectivity)
824 9199 : 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 9199 : output_unit=output_unit)
831 : ! Resize arrays to their proper size..
832 9199 : CALL reallocate(conn_info%bond_a, 1, n_bonds)
833 9199 : CALL reallocate(conn_info%bond_b, 1, n_bonds)
834 9199 : 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 314 : IF (.NOT. topology%reorder_atom) THEN
838 304 : topology%reorder_atom = .TRUE.
839 304 : IF (output_unit > 0) WRITE (output_unit, '(T2,"GENERATE|",A)') &
840 152 : " Molecules names have been generated. Now reordering particle set in order to have ", &
841 304 : " 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 8885 : atom_info, bondparm_factor, output_unit)
848 : END IF
849 18396 : IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN
850 0 : IF (output_unit > 0) THEN
851 : WRITE (output_unit, '(T2,"GENERATE|",A)') &
852 0 : " ERROR in connectivity generation!", &
853 0 : " The THRESHOLD to select possible bonds is bigger than the MAX bondlength", &
854 0 : " used to build the neighbors lists. Increase the BONDLENGTH_MAX patameter"
855 : WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') &
856 0 : " Present THRESHOLD (", my_maxrad*bondparm_factor, " )."// &
857 0 : " Present BONDLENGTH_MAX (", r_max(1, 1), " )"
858 : END IF
859 0 : CPABORT("Unable to generate connectivity")
860 : END IF
861 : END DO
862 9197 : IF (connectivity_ok .AND. (output_unit > 0)) THEN
863 : WRITE (output_unit, '(T2,"GENERATE|",A)') &
864 4700 : " Achieved consistency in connectivity generation."
865 : END IF
866 9197 : CALL fist_neighbor_deallocate(nonbonded)
867 9197 : CALL timestop(handle2)
868 9197 : CALL timeset(TRIM(routineN)//"_6", handle2)
869 : ! Deallocate temporary working arrays
870 236993 : DO I = 1, natom
871 227796 : DEALLOCATE (bond_list(I)%array1)
872 236993 : DEALLOCATE (bond_list(I)%array2)
873 : END DO
874 9197 : DEALLOCATE (bond_list)
875 9197 : DEALLOCATE (pbc_coord)
876 9197 : DEALLOCATE (radius)
877 9197 : DEALLOCATE (list)
878 9197 : CALL deallocate_particle_set(particle_set)
879 9197 : CALL deallocate_atomic_kind_set(atomic_kind_set)
880 : !
881 9197 : CALL timestop(handle2)
882 9197 : IF (output_unit > 0 .AND. n_bonds > 0) THEN
883 1182 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bonds generated:", &
884 2364 : n_bonds
885 : END IF
886 9197 : CALL timeset(TRIM(routineN)//"_7", handle2)
887 : ! If PARA_RES then activate RESIDUES
888 9197 : CALL reallocate(conn_info%c_bond_a, 1, 0)
889 9197 : CALL reallocate(conn_info%c_bond_b, 1, 0)
890 9197 : IF (topology%para_res) THEN
891 125249 : DO ibond = 1, SIZE(conn_info%bond_a)
892 116052 : iatom = conn_info%bond_a(ibond)
893 116052 : jatom = conn_info%bond_b(ibond)
894 : IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. &
895 116052 : (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. &
896 9197 : (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN
897 7214 : IF (iw > 0) WRITE (iw, *) " PARA_RES, bond between molecules atom ", &
898 4 : iatom, jatom
899 7212 : cbond = cbond + 1
900 7212 : CALL reallocate(conn_info%c_bond_a, 1, cbond)
901 7212 : CALL reallocate(conn_info%c_bond_b, 1, cbond)
902 7212 : conn_info%c_bond_a(cbond) = iatom
903 7212 : 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 9197 : CALL timestop(handle2)
912 9197 : DEALLOCATE (isolated_atoms)
913 9197 : CALL timestop(handle)
914 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
915 9197 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
916 101167 : 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 8885 : 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 8885 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mol_info_tmp
942 8885 : INTEGER, DIMENSION(:), POINTER :: mol_map, mol_map_o, wrk
943 8885 : INTEGER, DIMENSION(:, :), POINTER :: mol_info
944 8885 : LOGICAL, DIMENSION(:), POINTER :: icheck
945 8885 : TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
946 :
947 8885 : CALL timeset(routineN, handle)
948 8885 : conn_ok = .TRUE.
949 8885 : natom = SIZE(atom_info%id_atmname)
950 246559 : ALLOCATE (bond_list(natom))
951 228789 : DO I = 1, natom
952 228789 : ALLOCATE (bond_list(I)%array1(0))
953 : END DO
954 8885 : CALL reorder_structure(bond_list, bond_a, bond_b, SIZE(bond_a))
955 26655 : ALLOCATE (mol_map(natom))
956 17770 : ALLOCATE (mol_map_o(natom))
957 17770 : ALLOCATE (wrk(natom))
958 :
959 228789 : DO i = 1, natom
960 228789 : mol_map(i) = atom_info%id_molname(i)
961 : END DO
962 448693 : mol_map_o = mol_map
963 :
964 8885 : 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 26655 : ALLOCATE (mol_info_tmp(natom, 2))
979 :
980 8885 : itype = mol_map(1)
981 8885 : nsize = 1
982 8885 : idim = 1
983 8885 : mol_info_tmp(1, 1) = itype
984 219904 : DO i = 2, natom
985 219904 : IF (mol_map(i) /= itype) THEN
986 60147 : nsize = nsize + 1
987 60147 : itype = mol_map(i)
988 60147 : mol_info_tmp(nsize, 1) = itype
989 60147 : mol_info_tmp(nsize - 1, 2) = idim
990 60147 : idim = 1
991 : ELSE
992 150872 : idim = idim + 1
993 : END IF
994 : END DO
995 8885 : mol_info_tmp(nsize, 2) = idim
996 :
997 26655 : ALLOCATE (mol_info(nsize, 4))
998 164719 : mol_info(1:nsize, 1:2) = mol_info_tmp(1:nsize, 1:2)
999 8885 : DEALLOCATE (mol_info_tmp)
1000 :
1001 77917 : DO i = 1, nsize
1002 69032 : mol_info(i, 3) = 0
1003 77917 : mol_info(i, 4) = 0
1004 : END DO
1005 : !
1006 17770 : ALLOCATE (icheck(natom))
1007 228789 : icheck = .FALSE.
1008 228723 : DO i = 1, natom
1009 219840 : IF (icheck(i)) CYCLE
1010 112794 : itype = mol_map_o(i)
1011 112794 : mol_natom = 0
1012 112794 : CALL give_back_molecule(icheck, bond_list, i, mol_natom, mol_map_o, mol_map_o(i))
1013 13781513 : DO j = 1, SIZE(mol_info)
1014 13555925 : IF (itype == mol_info(j, 1)) EXIT
1015 : END DO
1016 112794 : mol_info(j, 3) = mol_info(j, 3) + 1
1017 112794 : IF (mol_info(j, 4) == 0) mol_info(j, 4) = mol_natom
1018 121678 : 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 2 : conn_ok = .FALSE.
1024 2 : bondparm_factor = bondparm_factor*1.05_dp
1025 2 : IF (output_unit < 0) EXIT
1026 1 : 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 1 : ' number of atoms.'
1030 1 : CALL integer_to_string(i, ctmp1)
1031 1 : CALL integer_to_string(mol_natom, ctmp2)
1032 1 : 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 1 : '> of atoms.', ' while the other same molecules have Nr. <'// &
1036 2 : TRIM(ctmp3)//'> of atoms!'
1037 : WRITE (output_unit, '(T2,"GENERATE|",A)') &
1038 1 : ' Increasing bondparm_factor by 1.05.. An error was found in the generated', &
1039 2 : ' connectivity. Retry...'
1040 : WRITE (output_unit, '(T2,"GENERATE|",A,F11.6,A,/)') &
1041 1 : " Present value of BONDPARM_FACTOR (", bondparm_factor, " )."
1042 1 : EXIT
1043 : END IF
1044 : END DO
1045 :
1046 8885 : DEALLOCATE (icheck)
1047 8885 : DEALLOCATE (mol_info)
1048 8885 : DEALLOCATE (mol_map)
1049 8885 : DEALLOCATE (mol_map_o)
1050 8885 : DEALLOCATE (wrk)
1051 228789 : DO I = 1, natom
1052 228789 : DEALLOCATE (bond_list(I)%array1)
1053 : END DO
1054 8885 : DEALLOCATE (bond_list)
1055 8885 : CALL timestop(handle)
1056 8885 : 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 29942 : 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 14971 : INTEGER, DIMENSION(:), POINTER :: atlist, Ilist1, Ilist2, Ilist3, Ilist4
1086 : LOGICAL :: explicit, ip3, ip4
1087 :
1088 14971 : natom = topology%natoms
1089 : ! Preliminary sort of arrays
1090 14971 : ip3 = PRESENT(Iarray3)
1091 14971 : ip4 = PRESENT(Iarray4)
1092 14971 : nsize = 2
1093 5772 : IF (ip3) nsize = nsize + 1
1094 14971 : IF (ip3 .AND. ip4) nsize = nsize + 1
1095 : ! Put the lists always in the canonical order
1096 14971 : CALL reorder_list_array(Iarray1, Iarray2, Iarray3, Iarray4, nsize, nvar)
1097 : ! Go on with external control
1098 14971 : CALL section_vals_get(section, explicit=explicit, n_repetition=n_rep)
1099 14971 : IF (explicit) THEN
1100 30 : NULLIFY (Ilist1, Ilist2, Ilist3, Ilist4, atlist)
1101 88 : ALLOCATE (Ilist1(nvar))
1102 58 : ALLOCATE (Ilist2(nvar))
1103 2702 : Ilist1 = Iarray1(1:nvar)
1104 2702 : Ilist2 = Iarray2(1:nvar)
1105 10 : SELECT CASE (nsize)
1106 : CASE (2) !do nothing
1107 : CASE (3)
1108 20 : ALLOCATE (Ilist3(nvar))
1109 706 : Ilist3 = Iarray3(1:nvar)
1110 : CASE (4)
1111 24 : ALLOCATE (Ilist3(nvar))
1112 24 : ALLOCATE (Ilist4(nvar))
1113 828 : Ilist3 = Iarray3(1:nvar)
1114 828 : Ilist4 = Iarray4(1:nvar)
1115 : CASE DEFAULT
1116 : ! Should never reach this point
1117 30 : CPABORT("Only 2, 3, 4 are supported as the value of nsize")
1118 : END SELECT
1119 30 : CALL list_canonical_order(Ilist1, Ilist2, Ilist3, Ilist4, nsize, is_impr)
1120 : !
1121 98 : DO i = 1, n_rep
1122 68 : 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 68 : i_val=do_action)
1125 : !
1126 180 : DO j = 1, n_rep_val
1127 : CALL section_vals_val_get(section, "ATOMS", i_rep_section=i, i_rep_val=j, &
1128 82 : i_vals=atlist)
1129 82 : CPASSERT(SIZE(atlist) == nsize)
1130 82 : CALL integer_to_string(nsize - 1, fmt)
1131 : CALL check_element_list(do_it, do_action, atlist, Ilist1, Ilist2, Ilist3, Ilist4, &
1132 82 : is_impr)
1133 150 : IF (do_action == do_add) THEN
1134 : ! Add to the element to the list
1135 42 : IF (do_it > 0) THEN
1136 26 : nvar = nvar + 1
1137 26 : IF (output_unit > 0) THEN
1138 : WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T64,A,I6)') &
1139 13 : "element (", &
1140 48 : atlist(1), (",", atlist(k), k=2, nsize), ") added.", " NEW size::", nvar
1141 : END IF
1142 26 : IF (nvar > SIZE(Iarray1)) THEN
1143 2 : new_size = INT(5 + 1.2*nvar)
1144 2 : CALL reallocate(Iarray1, 1, new_size)
1145 2 : CALL reallocate(Iarray2, 1, new_size)
1146 0 : SELECT CASE (nsize)
1147 : CASE (3)
1148 0 : CALL reallocate(Iarray3, 1, new_size)
1149 : CASE (4)
1150 0 : CALL reallocate(Iarray3, 1, new_size)
1151 2 : CALL reallocate(Iarray4, 1, new_size)
1152 : END SELECT
1153 : END IF
1154 : ! Using Ilist instead of atlist the canonical order is preserved..
1155 428 : Iarray1(do_it + 1:nvar) = Iarray1(do_it:nvar - 1)
1156 428 : Iarray2(do_it + 1:nvar) = Iarray2(do_it:nvar - 1)
1157 26 : Iarray1(do_it) = Ilist1(do_it)
1158 26 : Iarray2(do_it) = Ilist2(do_it)
1159 2 : SELECT CASE (nsize)
1160 : CASE (3)
1161 86 : Iarray3(do_it + 1:nvar) = Iarray3(do_it:nvar - 1)
1162 2 : Iarray3(do_it) = Ilist3(do_it)
1163 : CASE (4)
1164 230 : Iarray3(do_it + 1:nvar) = Iarray3(do_it:nvar - 1)
1165 230 : Iarray4(do_it + 1:nvar) = Iarray4(do_it:nvar - 1)
1166 8 : Iarray3(do_it) = Ilist3(do_it)
1167 34 : Iarray4(do_it) = Ilist4(do_it)
1168 : END SELECT
1169 : ELSE
1170 16 : IF (output_unit > 0) THEN
1171 : WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T80,A)') &
1172 8 : "element (", &
1173 30 : 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 40 : IF (do_it > 0) THEN
1179 34 : nvar = nvar - 1
1180 34 : IF (output_unit > 0) THEN
1181 : WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T64,A,I6)') &
1182 17 : "element (", &
1183 73 : atlist(1), (",", atlist(k), k=2, nsize), ") removed.", " NEW size::", nvar
1184 : END IF
1185 506 : Iarray1(do_it:nvar) = Iarray1(do_it + 1:nvar + 1)
1186 506 : Iarray2(do_it:nvar) = Iarray2(do_it + 1:nvar + 1)
1187 34 : Iarray1(nvar + 1) = -HUGE(0)
1188 34 : Iarray2(nvar + 1) = -HUGE(0)
1189 16 : SELECT CASE (nsize)
1190 : CASE (3)
1191 260 : Iarray3(do_it:nvar) = Iarray3(do_it + 1:nvar + 1)
1192 16 : Iarray3(nvar + 1) = -HUGE(0)
1193 : CASE (4)
1194 146 : Iarray3(do_it:nvar) = Iarray3(do_it + 1:nvar + 1)
1195 146 : Iarray4(do_it:nvar) = Iarray4(do_it + 1:nvar + 1)
1196 14 : Iarray3(nvar + 1) = -HUGE(0)
1197 48 : Iarray4(nvar + 1) = -HUGE(0)
1198 : END SELECT
1199 : ELSE
1200 6 : IF (output_unit > 0) THEN
1201 : WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T80,A)') &
1202 3 : "element (", &
1203 10 : 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 30 : DEALLOCATE (Ilist1)
1211 30 : DEALLOCATE (Ilist2)
1212 10 : SELECT CASE (nsize)
1213 : CASE (2) ! do nothing
1214 : CASE (3)
1215 10 : DEALLOCATE (Ilist3)
1216 : CASE (4)
1217 12 : DEALLOCATE (Ilist3)
1218 12 : DEALLOCATE (Ilist4)
1219 : CASE DEFAULT
1220 : ! Should never reach this point
1221 30 : CPABORT("Only 2, 3, 4 are supported as the value of nsize")
1222 : END SELECT
1223 : END IF
1224 14971 : 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 30 : 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 30 : do_impr = .FALSE.
1247 30 : IF (PRESENT(is_impr)) do_impr = is_impr
1248 38 : SELECT CASE (nsize)
1249 : CASE (2)
1250 588 : DO i = 1, SIZE(Ilist1)
1251 580 : tmp1 = Ilist1(i)
1252 580 : tmp2 = Ilist2(i)
1253 580 : Ilist1(i) = MIN(tmp1, tmp2)
1254 588 : Ilist2(i) = MAX(tmp1, tmp2)
1255 : END DO
1256 : CASE (3)
1257 358 : DO i = 1, SIZE(Ilist1)
1258 348 : tmp1 = Ilist1(i)
1259 348 : tmp2 = Ilist3(i)
1260 348 : Ilist1(i) = MIN(tmp1, tmp2)
1261 358 : Ilist3(i) = MAX(tmp1, tmp2)
1262 : END DO
1263 : CASE (4)
1264 438 : DO i = 1, SIZE(Ilist1)
1265 420 : IF (.NOT. do_impr) THEN
1266 372 : tmp1 = Ilist1(i)
1267 372 : tmp2 = Ilist4(i)
1268 372 : Ilist1(i) = MIN(tmp1, tmp2)
1269 372 : IF (Ilist1(i) == tmp2) THEN
1270 0 : tmp3 = Ilist3(i)
1271 0 : Ilist3(i) = Ilist2(i)
1272 0 : Ilist2(i) = tmp3
1273 : END IF
1274 372 : Ilist4(i) = MAX(tmp1, tmp2)
1275 : ELSE
1276 36 : tt(1) = Ilist2(i)
1277 36 : tt(2) = Ilist3(i)
1278 36 : tt(3) = Ilist4(i)
1279 36 : CALL sort(tt, 3, ss)
1280 36 : Ilist2(i) = tt(1)
1281 36 : Ilist3(i) = tt(2)
1282 36 : Ilist4(i) = tt(3)
1283 : END IF
1284 : END DO
1285 : END SELECT
1286 :
1287 30 : 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 82 : 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 82 : do_impr = .FALSE.
1315 82 : IF (PRESENT(is_impr)) do_impr = is_impr
1316 82 : found = .FALSE.
1317 82 : nsize = SIZE(atlist)
1318 82 : ndim = SIZE(Ilist1)
1319 322 : DO i = 1, nsize
1320 322 : tmp(i) = atlist(i)
1321 : END DO
1322 28 : SELECT CASE (nsize)
1323 : CASE (2)
1324 28 : tmp1 = tmp(1)
1325 28 : tmp2 = tmp(2)
1326 28 : tmp(1) = MIN(tmp1, tmp2)
1327 28 : tmp(2) = MAX(tmp1, tmp2)
1328 : CASE (3)
1329 32 : tmp1 = tmp(1)
1330 32 : tmp2 = tmp(3)
1331 32 : tmp(1) = MIN(tmp1, tmp2)
1332 32 : tmp(3) = MAX(tmp1, tmp2)
1333 : CASE (4)
1334 82 : IF (.NOT. do_impr) THEN
1335 10 : tmp1 = tmp(1)
1336 10 : tmp2 = tmp(4)
1337 10 : tmp(1) = MIN(tmp1, tmp2)
1338 10 : IF (tmp(1) == tmp2) THEN
1339 6 : tmp3 = tmp(3)
1340 6 : tmp(3) = tmp(2)
1341 6 : tmp(2) = tmp3
1342 : END IF
1343 10 : tmp(4) = MAX(tmp1, tmp2)
1344 : ELSE
1345 12 : tt(1) = tmp(2)
1346 12 : tt(2) = tmp(3)
1347 12 : tt(3) = tmp(4)
1348 12 : CALL sort(tt, 3, ss)
1349 12 : tmp(2) = tt(1)
1350 12 : tmp(3) = tt(2)
1351 12 : tmp(4) = tt(3)
1352 : END IF
1353 : END SELECT
1354 : ! boundary to search
1355 1788 : DO istart = 1, ndim
1356 1788 : IF (Ilist1(istart) >= tmp(1)) EXIT
1357 : END DO
1358 : ! if nothing there stay within bounds
1359 82 : IF (istart <= ndim) THEN
1360 76 : IF (Ilist1(istart) > tmp(1) .AND. (istart /= 1)) istart = istart - 1
1361 : END IF
1362 222 : DO iend = istart, ndim
1363 222 : IF (Ilist1(iend) /= tmp(1)) EXIT
1364 : END DO
1365 82 : IF (iend == ndim + 1) iend = ndim
1366 : ! Final search in array
1367 : SELECT CASE (nsize)
1368 : CASE (2)
1369 40 : DO i = istart, iend
1370 28 : IF ((Ilist1(i) > tmp(1)) .OR. (Ilist2(i) > tmp(2))) EXIT
1371 40 : 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 40 : DO i = istart, iend
1378 40 : IF ((Ilist1(i) > tmp(1)) .OR. (Ilist2(i) > tmp(2)) .OR. (Ilist3(i) > tmp(3))) EXIT
1379 40 : 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 106 : DO i = istart, iend
1386 22 : 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 24 : .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 124 : SELECT CASE (do_action)
1395 : CASE (do_add)
1396 42 : IF (found) THEN
1397 16 : 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 26 : do_it = i
1408 26 : new_size = ndim + 1
1409 26 : CALL reallocate(Ilist1, 1, new_size)
1410 26 : CALL reallocate(Ilist2, 1, new_size)
1411 428 : Ilist1(i + 1:new_size) = Ilist1(i:ndim)
1412 428 : Ilist2(i + 1:new_size) = Ilist2(i:ndim)
1413 26 : Ilist1(i) = tmp(1)
1414 26 : Ilist2(i) = tmp(2)
1415 2 : SELECT CASE (nsize)
1416 : CASE (3)
1417 2 : CALL reallocate(Ilist3, 1, new_size)
1418 86 : Ilist3(i + 1:new_size) = Ilist3(i:ndim)
1419 2 : Ilist3(i) = tmp(3)
1420 : CASE (4)
1421 8 : CALL reallocate(Ilist3, 1, new_size)
1422 8 : CALL reallocate(Ilist4, 1, new_size)
1423 230 : Ilist3(i + 1:new_size) = Ilist3(i:ndim)
1424 230 : Ilist4(i + 1:new_size) = Ilist4(i:ndim)
1425 8 : Ilist3(i) = tmp(3)
1426 34 : Ilist4(i) = tmp(4)
1427 : END SELECT
1428 : END IF
1429 : CASE (do_remove)
1430 82 : IF (found) THEN
1431 34 : do_it = i
1432 : ! Let's delete the element in position do_it
1433 34 : new_size = ndim - 1
1434 506 : Ilist1(i:new_size) = Ilist1(i + 1:ndim)
1435 506 : Ilist2(i:new_size) = Ilist2(i + 1:ndim)
1436 34 : CALL reallocate(Ilist1, 1, new_size)
1437 34 : CALL reallocate(Ilist2, 1, new_size)
1438 16 : SELECT CASE (nsize)
1439 : CASE (3)
1440 260 : Ilist3(i:new_size) = Ilist3(i + 1:ndim)
1441 16 : CALL reallocate(Ilist3, 1, new_size)
1442 : CASE (4)
1443 146 : Ilist3(i:new_size) = Ilist3(i + 1:ndim)
1444 146 : Ilist4(i:new_size) = Ilist4(i + 1:ndim)
1445 14 : CALL reallocate(Ilist3, 1, new_size)
1446 48 : CALL reallocate(Ilist4, 1, new_size)
1447 : END SELECT
1448 : ELSE
1449 6 : 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 82 : 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 117086 : 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 117086 : old_size = SIZE(conn_info%bond_a)
1472 117086 : IF (n_bonds > old_size) THEN
1473 6032 : new_size = INT(5 + 1.2*old_size)
1474 6032 : CALL reallocate(conn_info%bond_a, 1, new_size)
1475 6032 : CALL reallocate(conn_info%bond_b, 1, new_size)
1476 : END IF
1477 117086 : conn_info%bond_a(n_bonds) = atm1
1478 117086 : conn_info%bond_b(n_bonds) = atm2
1479 117086 : 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 21214 : 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 10607 : 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 10607 : NULLIFY (logger)
1501 21214 : logger => cp_get_default_logger()
1502 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1503 10607 : extension=".subsysLog")
1504 10607 : CALL timeset(routineN, handle)
1505 10607 : output_unit = cp_logger_get_default_io_unit(logger)
1506 10607 : conn_info => topology%conn_info
1507 10607 : nbond = 0
1508 10607 : ntheta = 0
1509 10607 : natom = topology%natoms
1510 : ! This call is for connectivity off
1511 10607 : IF (ASSOCIATED(conn_info%bond_a)) THEN
1512 8881 : nbond = SIZE(conn_info%bond_a)
1513 : ELSE
1514 1726 : CALL reallocate(conn_info%bond_a, 1, nbond)
1515 1726 : CALL reallocate(conn_info%bond_b, 1, nbond)
1516 : END IF
1517 10607 : IF (nbond /= 0) THEN
1518 1924 : nsize = INT(5 + 1.2*ntheta)
1519 1924 : CALL reallocate(conn_info%theta_a, 1, nsize)
1520 1924 : CALL reallocate(conn_info%theta_b, 1, nsize)
1521 1924 : CALL reallocate(conn_info%theta_c, 1, nsize)
1522 : ! Get list of bonds to pre-process theta
1523 156946 : ALLOCATE (bond_list(natom))
1524 153098 : DO I = 1, natom
1525 153098 : ALLOCATE (bond_list(I)%array1(0))
1526 : END DO
1527 1924 : 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 1924 : 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 1924 : Oarray3=conn_info%theta_c)
1537 1924 : CALL timestop(handle2)
1538 153098 : DO I = 1, natom
1539 153098 : DEALLOCATE (bond_list(I)%array1)
1540 : END DO
1541 1924 : DEALLOCATE (bond_list)
1542 1924 : IF (output_unit > 0) THEN
1543 1024 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Bends generated:", &
1544 2048 : ntheta
1545 : END IF
1546 : ! External defined bends (useful for complex connectivity)
1547 1924 : 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 3848 : output_unit=output_unit)
1555 : END IF
1556 : ! Resize arrays to their proper size..
1557 10607 : CALL reallocate(conn_info%theta_a, 1, ntheta)
1558 10607 : CALL reallocate(conn_info%theta_b, 1, ntheta)
1559 10607 : CALL reallocate(conn_info%theta_c, 1, ntheta)
1560 10607 : IF (output_unit > 0 .AND. ntheta > 0) THEN
1561 970 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bends generated:", &
1562 1940 : ntheta
1563 : END IF
1564 10607 : CALL timestop(handle)
1565 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1566 10607 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1567 10607 : 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 903320 : RECURSIVE SUBROUTINE match_iterative_path(Iarray1, Iarray2, Iarray3, &
1587 903320 : 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 903320 : INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list
1600 : LOGICAL :: check
1601 903320 : TYPE(array1_list_type), DIMENSION(:), POINTER :: wrk
1602 :
1603 903320 : check = max_levl >= 2 .AND. max_levl <= 4
1604 0 : CPASSERT(check)
1605 903320 : IF (.NOT. PRESENT(Ilist)) THEN
1606 0 : SELECT CASE (max_levl)
1607 : CASE (2)
1608 0 : CPASSERT(.NOT. PRESENT(Iarray2))
1609 0 : CPASSERT(.NOT. PRESENT(Iarray3))
1610 0 : CPASSERT(.NOT. PRESENT(Oarray3))
1611 0 : CPASSERT(.NOT. PRESENT(Oarray4))
1612 : CASE (3)
1613 1924 : CPASSERT(PRESENT(Iarray2))
1614 1924 : CPASSERT(.NOT. PRESENT(Iarray3))
1615 1924 : CPASSERT(PRESENT(Oarray3))
1616 1924 : CPASSERT(.NOT. PRESENT(Oarray4))
1617 : CASE (4)
1618 1924 : CPASSERT(PRESENT(Iarray2))
1619 1924 : CPASSERT(PRESENT(Iarray3))
1620 1924 : CPASSERT(PRESENT(Oarray3))
1621 5772 : CPASSERT(PRESENT(Oarray4))
1622 : END SELECT
1623 : END IF
1624 903320 : natom = SIZE(Iarray1)
1625 903320 : IF (.NOT. PRESENT(Ilist)) THEN
1626 : ! Start a new loop.. Only the first time the routine is called
1627 11544 : ALLOCATE (my_list(max_levl))
1628 306196 : DO i = 1, natom
1629 302348 : my_levl = 1
1630 1360566 : my_list = -1
1631 302348 : 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 306196 : Ilist=my_list)
1643 : END DO
1644 3848 : DEALLOCATE (my_list)
1645 : ELSE
1646 1201820 : SELECT CASE (it_levl)
1647 : CASE (2)
1648 302348 : wrk => Iarray1
1649 : CASE (3)
1650 429032 : wrk => Iarray2
1651 : CASE (4)
1652 899472 : wrk => Iarray3
1653 : END SELECT
1654 899472 : i = Ilist(it_levl - 1)
1655 2352140 : DO j = 1, SIZE(Iarray1(i)%array1)
1656 1452668 : ind = wrk(i)%array1(j)
1657 4658888 : IF (ANY(Ilist == ind)) CYCLE
1658 1754632 : IF (it_levl < max_levl) THEN
1659 597124 : 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 597124 : Ilist=Ilist)
1671 597124 : Ilist(it_levl) = -1
1672 258036 : ELSEIF (it_levl == max_levl) THEN
1673 258036 : IF (Ilist(1) > ind) CYCLE
1674 129018 : Ilist(it_levl) = ind
1675 129018 : nvar = nvar + 1
1676 0 : SELECT CASE (it_levl)
1677 : CASE (2)
1678 0 : IF (nvar > SIZE(Oarray1)) THEN
1679 0 : CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar))
1680 0 : CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar))
1681 : END IF
1682 0 : Oarray1(nvar) = Ilist(1)
1683 0 : Oarray2(nvar) = Ilist(2)
1684 : CASE (3)
1685 84046 : IF (nvar > SIZE(Oarray1)) THEN
1686 3204 : CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar))
1687 3204 : CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar))
1688 3204 : CALL reallocate(Oarray3, 1, INT(5 + 1.2*nvar))
1689 : END IF
1690 84046 : Oarray1(nvar) = Ilist(1)
1691 84046 : Oarray2(nvar) = Ilist(2)
1692 84046 : Oarray3(nvar) = Ilist(3)
1693 : CASE (4)
1694 44972 : IF (nvar > SIZE(Oarray1)) THEN
1695 1436 : CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar))
1696 1436 : CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar))
1697 1436 : CALL reallocate(Oarray3, 1, INT(5 + 1.2*nvar))
1698 1436 : CALL reallocate(Oarray4, 1, INT(5 + 1.2*nvar))
1699 : END IF
1700 44972 : Oarray1(nvar) = Ilist(1)
1701 44972 : Oarray2(nvar) = Ilist(2)
1702 44972 : Oarray3(nvar) = Ilist(3)
1703 44972 : Oarray4(nvar) = Ilist(4)
1704 : CASE DEFAULT
1705 : !should never reach this point
1706 129018 : CPABORT("Only 2, 3, 4 are supported as the value of it_levl")
1707 : END SELECT
1708 129018 : Ilist(it_levl) = -1
1709 : ELSE
1710 : !should never reach this point
1711 0 : CPABORT("it_levl exceeds max_levl in match_iterative_path")
1712 : END IF
1713 : END DO
1714 : END IF
1715 1806640 : 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 21214 : 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 10607 : NULLIFY (logger)
1735 10607 : logger => cp_get_default_logger()
1736 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1737 10607 : extension=".subsysLog")
1738 10607 : output_unit = cp_logger_get_default_io_unit(logger)
1739 10607 : CALL timeset(routineN, handle)
1740 10607 : conn_info => topology%conn_info
1741 10607 : ntheta = SIZE(conn_info%theta_a)
1742 10607 : CALL reallocate(conn_info%ub_a, 1, ntheta)
1743 10607 : CALL reallocate(conn_info%ub_b, 1, ntheta)
1744 10607 : CALL reallocate(conn_info%ub_c, 1, ntheta)
1745 :
1746 94639 : DO itheta = 1, ntheta
1747 84032 : conn_info%ub_a(itheta) = conn_info%theta_a(itheta)
1748 84032 : conn_info%ub_b(itheta) = conn_info%theta_b(itheta)
1749 94639 : conn_info%ub_c(itheta) = conn_info%theta_c(itheta)
1750 : END DO
1751 10607 : IF (output_unit > 0 .AND. ntheta > 0) THEN
1752 970 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of UB generated:", &
1753 1940 : ntheta
1754 : END IF
1755 10607 : CALL timestop(handle)
1756 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1757 10607 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1758 :
1759 10607 : 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 21214 : 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 10607 : 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 10607 : NULLIFY (logger)
1781 21214 : logger => cp_get_default_logger()
1782 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1783 10607 : extension=".subsysLog")
1784 10607 : output_unit = cp_logger_get_default_io_unit(logger)
1785 10607 : CALL timeset(routineN, handle)
1786 10607 : conn_info => topology%conn_info
1787 10607 : nphi = 0
1788 10607 : nbond = SIZE(conn_info%bond_a)
1789 10607 : IF (nbond /= 0) THEN
1790 1924 : nsize = INT(5 + 1.2*nphi)
1791 1924 : CALL reallocate(conn_info%phi_a, 1, nsize)
1792 1924 : CALL reallocate(conn_info%phi_b, 1, nsize)
1793 1924 : CALL reallocate(conn_info%phi_c, 1, nsize)
1794 1924 : CALL reallocate(conn_info%phi_d, 1, nsize)
1795 : ! Get list of bonds to pre-process phi
1796 1924 : natom = topology%natoms
1797 156946 : ALLOCATE (bond_list(natom))
1798 153098 : DO I = 1, natom
1799 153098 : ALLOCATE (bond_list(I)%array1(0))
1800 : END DO
1801 1924 : 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 1924 : Oarray4=conn_info%phi_d)
1812 153098 : DO I = 1, natom
1813 153098 : DEALLOCATE (bond_list(I)%array1)
1814 : END DO
1815 1924 : DEALLOCATE (bond_list)
1816 1924 : IF (output_unit > 0) THEN
1817 1024 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Torsions generated:", &
1818 2048 : nphi
1819 : END IF
1820 : ! External defined torsions (useful for complex connectivity)
1821 1924 : 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 1924 : output_unit=output_unit)
1830 : END IF
1831 : ! Resize arrays to their proper size..
1832 10607 : CALL reallocate(conn_info%phi_a, 1, nphi)
1833 10607 : CALL reallocate(conn_info%phi_b, 1, nphi)
1834 10607 : CALL reallocate(conn_info%phi_c, 1, nphi)
1835 10607 : CALL reallocate(conn_info%phi_d, 1, nphi)
1836 10607 : IF (output_unit > 0 .AND. nphi > 0) THEN
1837 222 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Torsions generated:", &
1838 444 : nphi
1839 : END IF
1840 10607 : CALL timestop(handle)
1841 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1842 10607 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1843 :
1844 10607 : 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 21214 : 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 10607 : 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 10607 : NULLIFY (logger)
1869 21214 : logger => cp_get_default_logger()
1870 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1871 10607 : extension=".subsysLog")
1872 10607 : output_unit = cp_logger_get_default_io_unit(logger)
1873 10607 : CALL timeset(routineN, handle)
1874 10607 : atom_info => topology%atom_info
1875 10607 : conn_info => topology%conn_info
1876 10607 : natom = topology%natoms
1877 10607 : nimpr = 0
1878 10607 : nbond = SIZE(conn_info%bond_a)
1879 10607 : IF (nbond /= 0) THEN
1880 1924 : nsize = INT(5 + 1.2*nimpr)
1881 1924 : CALL reallocate(conn_info%impr_a, 1, nsize)
1882 1924 : CALL reallocate(conn_info%impr_b, 1, nsize)
1883 1924 : CALL reallocate(conn_info%impr_c, 1, nsize)
1884 1924 : CALL reallocate(conn_info%impr_d, 1, nsize)
1885 : ! Get list of bonds to pre-process phi
1886 156946 : ALLOCATE (bond_list(natom))
1887 153098 : DO I = 1, natom
1888 153098 : ALLOCATE (bond_list(I)%array1(0))
1889 : END DO
1890 1924 : CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond)
1891 153098 : DO I = 1, natom
1892 : ! Count all atoms with three bonds
1893 153098 : IF (SIZE(bond_list(I)%array1) == 3) THEN
1894 : ! Problematic cases::
1895 : ! Nitrogen
1896 3342 : accept_impr = .TRUE.
1897 3342 : atm_symbol = TRIM(id2str(atom_info%id_element(i)))
1898 3342 : CALL uppercase(atm_symbol)
1899 3342 : 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 8712 : DO j = 1, 3
1904 6534 : ind = bond_list(I)%array1(j)
1905 8712 : IF (SIZE(bond_list(ind)%array1) == 3) accept_impr = .TRUE.
1906 : END DO
1907 : END IF
1908 2178 : IF (.NOT. accept_impr) CYCLE
1909 1914 : nimpr = nimpr + 1
1910 1914 : IF (nimpr > SIZE(conn_info%impr_a)) THEN
1911 136 : nsize = INT(5 + 1.2*nimpr)
1912 136 : CALL reallocate(conn_info%impr_a, 1, nsize)
1913 136 : CALL reallocate(conn_info%impr_b, 1, nsize)
1914 136 : CALL reallocate(conn_info%impr_c, 1, nsize)
1915 136 : CALL reallocate(conn_info%impr_d, 1, nsize)
1916 : END IF
1917 1914 : conn_info%impr_a(nimpr) = i
1918 1914 : conn_info%impr_b(nimpr) = bond_list(I)%array1(1)
1919 1914 : conn_info%impr_c(nimpr) = bond_list(I)%array1(2)
1920 1914 : conn_info%impr_d(nimpr) = bond_list(I)%array1(3)
1921 : END IF
1922 : END DO
1923 153098 : DO I = 1, natom
1924 153098 : DEALLOCATE (bond_list(I)%array1)
1925 : END DO
1926 1924 : DEALLOCATE (bond_list)
1927 : ! External defined impropers (useful for complex connectivity)
1928 1924 : 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 1924 : is_impr=.TRUE.)
1938 : END IF
1939 : ! Resize arrays to their proper size..
1940 10607 : CALL reallocate(conn_info%impr_a, 1, nimpr)
1941 10607 : CALL reallocate(conn_info%impr_b, 1, nimpr)
1942 10607 : CALL reallocate(conn_info%impr_c, 1, nimpr)
1943 10607 : CALL reallocate(conn_info%impr_d, 1, nimpr)
1944 10607 : IF (output_unit > 0 .AND. nimpr > 0) THEN
1945 44 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Impropers generated:", &
1946 88 : nimpr
1947 : END IF
1948 10607 : CALL timestop(handle)
1949 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1950 10607 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
1951 :
1952 10607 : 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 10607 : 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 10607 : 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 10607 : NULLIFY (logger)
1972 21214 : logger => cp_get_default_logger()
1973 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", &
1974 10607 : extension=".subsysLog")
1975 10607 : output_unit = cp_logger_get_default_io_unit(logger)
1976 10607 : CALL timeset(routineN, handle)
1977 :
1978 10607 : conn_info => topology%conn_info
1979 10607 : natom = topology%natoms
1980 :
1981 : ! Get list of bonds (sic). Get a list of bonded neighbors for every atom.
1982 313289 : ALLOCATE (bond_list(natom))
1983 292075 : DO i = 1, natom
1984 292075 : ALLOCATE (bond_list(i)%array1(0))
1985 : END DO
1986 10607 : nbond = SIZE(conn_info%bond_a)
1987 10607 : 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 302682 : ALLOCATE (theta_list(natom))
1991 292075 : DO i = 1, natom
1992 292075 : ALLOCATE (theta_list(i)%array1(0))
1993 : END DO
1994 10607 : ntheta = SIZE(conn_info%theta_a)
1995 10607 : 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 302682 : ALLOCATE (phi_list(natom))
1999 292075 : DO i = 1, natom
2000 292075 : ALLOCATE (phi_list(i)%array1(0))
2001 : END DO
2002 10607 : nphi = SIZE(conn_info%phi_a)
2003 10607 : CALL reorder_structure(phi_list, conn_info%phi_a, conn_info%phi_d, nphi)
2004 :
2005 : ! Allocate enough (possible too much)
2006 10607 : CALL reallocate(conn_info%onfo_a, 1, nphi)
2007 10607 : CALL reallocate(conn_info%onfo_b, 1, nphi)
2008 :
2009 10607 : ionfo = 0
2010 292075 : DO atom_a = 1, natom
2011 382015 : DO i = 1, SIZE(phi_list(atom_a)%array1)
2012 89940 : atom_b = phi_list(atom_a)%array1(i)
2013 : ! Avoid trivial duplicates.
2014 89940 : IF (atom_a > atom_b) CYCLE
2015 : ! Avoid onfo's in 4-rings.
2016 157018 : IF (ANY(atom_b == bond_list(atom_a)%array1)) CYCLE
2017 : ! Avoid onfo's in 5-rings.
2018 206964 : IF (ANY(atom_b == theta_list(atom_a)%array1)) CYCLE
2019 : ! Avoid onfo's in 6-rings.
2020 214908 : IF (ANY(atom_b == phi_list(atom_a)%array1(:i - 1))) CYCLE
2021 44578 : ionfo = ionfo + 1
2022 44578 : conn_info%onfo_a(ionfo) = atom_a
2023 371408 : conn_info%onfo_b(ionfo) = atom_b
2024 : END DO
2025 : END DO
2026 :
2027 : ! Reallocate such that just enough memory is used.
2028 10607 : CALL reallocate(conn_info%onfo_a, 1, ionfo)
2029 10607 : CALL reallocate(conn_info%onfo_b, 1, ionfo)
2030 :
2031 : ! Deallocate bond_list
2032 292075 : DO i = 1, natom
2033 292075 : DEALLOCATE (bond_list(i)%array1)
2034 : END DO
2035 10607 : DEALLOCATE (bond_list)
2036 : ! Deallocate theta_list
2037 292075 : DO i = 1, natom
2038 292075 : DEALLOCATE (theta_list(i)%array1)
2039 : END DO
2040 10607 : DEALLOCATE (theta_list)
2041 : ! Deallocate phi_list
2042 292075 : DO i = 1, natom
2043 292075 : DEALLOCATE (phi_list(i)%array1)
2044 : END DO
2045 10607 : DEALLOCATE (phi_list)
2046 :
2047 : ! Final output
2048 10607 : IF (output_unit > 0 .AND. ionfo > 0) THEN
2049 222 : WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of 1-4 interactions generated:", &
2050 444 : ionfo
2051 : END IF
2052 10607 : CALL timestop(handle)
2053 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
2054 10607 : "PRINT%TOPOLOGY_INFO/GENERATE_INFO")
2055 :
2056 21214 : END SUBROUTINE topology_generate_onfo
2057 :
2058 : END MODULE topology_generate_util
|
