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