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