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