Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Collection of subroutine needed for topology related things
10 : !> \par History
11 : !> jgh (23-05-2004) Last atom of molecule information added
12 : ! **************************************************************************************************
13 :
14 : MODULE topology_connectivity_util
15 : USE cp_log_handling, ONLY: cp_get_default_logger,&
16 : cp_logger_get_default_io_unit,&
17 : cp_logger_type
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE force_field_kind_types, ONLY: do_ff_charmm,&
21 : do_ff_harmonic
22 : USE input_constants, ONLY: do_conn_g87,&
23 : do_conn_g96,&
24 : do_conn_user
25 : USE input_section_types, ONLY: section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: default_string_length
28 : USE memory_utilities, ONLY: reallocate
29 : USE molecule_kind_types, ONLY: &
30 : allocate_molecule_kind_set, atom_type, bend_type, bond_type, get_molecule_kind, impr_type, &
31 : molecule_kind_type, opbend_type, set_molecule_kind, torsion_type, ub_type
32 : USE molecule_types, ONLY: allocate_molecule_set,&
33 : get_molecule,&
34 : local_states_type,&
35 : molecule_type,&
36 : set_molecule,&
37 : set_molecule_set
38 : USE string_table, ONLY: id2str
39 : USE topology_types, ONLY: atom_info_type,&
40 : connectivity_info_type,&
41 : topology_parameters_type
42 : USE util, ONLY: sort
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_connectivity_util'
48 :
49 : PRIVATE
50 : PUBLIC :: topology_connectivity_pack, topology_conn_multiple
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief topology connectivity pack
56 : !> \param molecule_kind_set ...
57 : !> \param molecule_set ...
58 : !> \param topology ...
59 : !> \param subsys_section ...
60 : !> \par History 11/2009 (Louis Vanduyhuys): added Out of Plane bends based on
61 : !> impropers in topology
62 : ! **************************************************************************************************
63 10274 : SUBROUTINE topology_connectivity_pack(molecule_kind_set, molecule_set, &
64 : topology, subsys_section)
65 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
66 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
67 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
68 : TYPE(section_vals_type), POINTER :: subsys_section
69 :
70 : CHARACTER(len=*), PARAMETER :: routineN = 'topology_connectivity_pack'
71 :
72 : CHARACTER(LEN=default_string_length) :: name
73 : INTEGER :: counter, first, handle, handle2, i, ibend, ibond, idim, iimpr, ikind, imol, &
74 : inter_bends, inter_bonds, inter_imprs, inter_torsions, inter_ubs, intra_bends, &
75 : intra_bonds, intra_imprs, intra_torsions, intra_ubs, inum, ires, istart_mol, istart_typ, &
76 : itorsion, ityp, iub, iw, j, j1, j2, j3, j4, jind, last, min_index, natom, nelectron, &
77 : nsgf, nval_tot1, nval_tot2, nvar1, nvar2, output_unit, stat
78 10274 : INTEGER, DIMENSION(:), POINTER :: c_var_a, c_var_b, c_var_c, c_var_d, c_var_type, &
79 10274 : first_list, last_list, map_atom_mol, map_atom_type, map_cvar_mol, map_cvars, map_var_mol, &
80 10274 : map_vars, molecule_list
81 10274 : INTEGER, DIMENSION(:, :), POINTER :: bnd_ctype, bnd_type
82 : LOGICAL :: found, found_last
83 : TYPE(atom_info_type), POINTER :: atom_info
84 10274 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
85 10274 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
86 10274 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
87 : TYPE(connectivity_info_type), POINTER :: conn_info
88 : TYPE(cp_logger_type), POINTER :: logger
89 10274 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
90 10274 : TYPE(local_states_type), DIMENSION(:), POINTER :: lmi
91 : TYPE(molecule_kind_type), POINTER :: molecule_kind
92 : TYPE(molecule_type), POINTER :: molecule
93 10274 : TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
94 10274 : TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
95 10274 : TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
96 :
97 10274 : NULLIFY (logger)
98 20548 : logger => cp_get_default_logger()
99 10274 : output_unit = cp_logger_get_default_io_unit(logger)
100 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
101 10274 : extension=".subsysLog")
102 10274 : CALL timeset(routineN, handle)
103 :
104 10274 : atom_info => topology%atom_info
105 10274 : conn_info => topology%conn_info
106 30822 : ALLOCATE (map_atom_mol(topology%natoms))
107 20548 : ALLOCATE (map_atom_type(topology%natoms))
108 : !-----------------------------------------------------------------------------
109 : !-----------------------------------------------------------------------------
110 : ! 1. Set the topology%[nmol_type,nmol,nmol_conn]
111 : !-----------------------------------------------------------------------------
112 10274 : CALL timeset(routineN//"_1", handle2)
113 10274 : natom = topology%natoms
114 10274 : topology%nmol = 1
115 10274 : topology%nmol_type = 1
116 10274 : topology%nmol_conn = 0
117 764843 : map_atom_mol = -1
118 764843 : map_atom_type = -1
119 10274 : map_atom_mol(1) = 1
120 10274 : map_atom_type(1) = 1
121 754569 : DO i = 1, natom - 1
122 744295 : IF ((atom_info%map_mol_typ(i + 1) /= atom_info%map_mol_typ(i)) .OR. &
123 : ((atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i)) .AND. &
124 : (.NOT. (topology%conn_type == do_conn_user)))) THEN
125 126155 : topology%nmol_type = topology%nmol_type + 1
126 : END IF
127 744295 : map_atom_type(i + 1) = topology%nmol_type
128 : IF ((atom_info%map_mol_typ(i + 1) /= atom_info%map_mol_typ(i)) .OR. &
129 744295 : (atom_info%map_mol_num(i + 1) /= atom_info%map_mol_num(i)) .OR. &
130 : (atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i))) THEN
131 294998 : topology%nmol = topology%nmol + 1
132 : END IF
133 744295 : map_atom_mol(i + 1) = topology%nmol
134 : IF ((atom_info%map_mol_typ(i + 1) == atom_info%map_mol_typ(i)) .AND. &
135 744295 : (atom_info%map_mol_num(i + 1) == atom_info%map_mol_num(i)) .AND. &
136 10274 : (atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i))) THEN
137 2832 : topology%nmol_conn = topology%nmol_conn + 1
138 : END IF
139 : END DO
140 10274 : IF (iw > 0) WRITE (iw, *) "topology%nmol ::", topology%nmol
141 10274 : IF (iw > 0) WRITE (iw, *) "topology%nmol_type ::", topology%nmol_type
142 10274 : IF (iw > 0) WRITE (iw, *) "topology%nmol_conn ::", topology%nmol_conn
143 10274 : CALL timestop(handle2)
144 : !-----------------------------------------------------------------------------
145 : !-----------------------------------------------------------------------------
146 : ! 1.1 Clean the temporary arrays to avoid quadratic loops around
147 : ! after this fix all topology_pack will be linear scaling
148 : !-----------------------------------------------------------------------------
149 10274 : CALL timeset(routineN//"_1.1", handle2)
150 10274 : istart_mol = map_atom_mol(1)
151 10274 : istart_typ = map_atom_type(1)
152 754569 : DO i = 2, natom
153 754569 : IF ((map_atom_mol(i) /= istart_mol) .AND. (map_atom_type(i) == istart_typ)) THEN
154 505249 : map_atom_mol(i) = -map_atom_mol(i)
155 239046 : ELSE IF (map_atom_type(i) /= istart_typ) THEN
156 126155 : istart_mol = map_atom_mol(i)
157 126155 : istart_typ = map_atom_type(i)
158 : END IF
159 : END DO
160 10274 : CALL timestop(handle2)
161 : !-----------------------------------------------------------------------------
162 : !-----------------------------------------------------------------------------
163 : ! 2. Allocate the molecule_kind_set
164 : !-----------------------------------------------------------------------------
165 10274 : CALL timeset(routineN//"_2", handle2)
166 10274 : IF (topology%nmol_type <= 0) THEN
167 0 : CPABORT("No molecule kind defined")
168 : ELSE
169 10274 : NULLIFY (molecule_kind_set)
170 10274 : i = topology%nmol_type
171 10274 : CALL allocate_molecule_kind_set(molecule_kind_set, i)
172 10300 : IF (iw > 0) WRITE (iw, *) " Allocated molecule_kind_set, Dimenstion of ", &
173 52 : SIZE(molecule_kind_set)
174 : END IF
175 10274 : CALL timestop(handle2)
176 : !-----------------------------------------------------------------------------
177 : !-----------------------------------------------------------------------------
178 : ! 3. Allocate the molecule_set
179 : !-----------------------------------------------------------------------------
180 10274 : CALL timeset(routineN//"_3", handle2)
181 10274 : IF (topology%nmol <= 0) THEN
182 0 : CPABORT("No molecule defined")
183 : ELSE
184 10274 : NULLIFY (molecule_set)
185 10274 : i = topology%nmol
186 10274 : CALL allocate_molecule_set(molecule_set, i)
187 10300 : IF (iw > 0) WRITE (iw, *) " Allocated molecule_set, dimension of ", &
188 52 : topology%nmol
189 : END IF
190 10274 : CALL timestop(handle2)
191 : !-----------------------------------------------------------------------------
192 : !-----------------------------------------------------------------------------
193 : ! 4. Set the molecule_kind_set%[kind_number,name,nsgf,nelectron]
194 : !-----------------------------------------------------------------------------
195 10274 : CALL timeset(routineN//"_4", handle2)
196 10274 : natom = topology%natoms
197 10274 : ikind = -1
198 764843 : DO i = 1, natom
199 764843 : IF (ikind /= map_atom_type(i)) THEN
200 136429 : ikind = map_atom_type(i)
201 136429 : molecule_kind => molecule_kind_set(ikind)
202 136429 : nsgf = 0
203 136429 : nelectron = 0
204 136429 : name = TRIM(id2str(atom_info%id_molname(i)))
205 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
206 : kind_number=ikind, &
207 : molname_generated=topology%molname_generated, &
208 : name=TRIM(name), &
209 : nsgf=nsgf, &
210 136429 : nelectron=nelectron)
211 : END IF
212 : END DO
213 10274 : CALL timestop(handle2)
214 : !-----------------------------------------------------------------------------
215 : !-----------------------------------------------------------------------------
216 : ! 5. Set the molecule_list for molecule_kind in molecule_kind_set
217 : !-----------------------------------------------------------------------------
218 10274 : CALL timeset(routineN//"_5", handle2)
219 10274 : natom = topology%natoms
220 10274 : ikind = map_atom_type(1)
221 10274 : imol = ABS(map_atom_mol(1))
222 10274 : counter = 0
223 290002 : DO i = 1, natom - 1
224 282473 : IF (ikind /= map_atom_type(i + 1)) THEN
225 126155 : found = .TRUE.
226 126155 : found_last = .FALSE.
227 126155 : imol = ABS(map_atom_mol(i))
228 156318 : ELSEIF (ikind == topology%nmol_type) THEN
229 2745 : found = .TRUE.
230 2745 : found_last = .TRUE.
231 2745 : imol = ABS(map_atom_mol(natom))
232 : ELSE
233 : found = .FALSE.
234 : found_last = .FALSE.
235 : END IF
236 :
237 10274 : IF (found) THEN
238 386700 : ALLOCATE (molecule_list(imol - counter))
239 426643 : DO j = 1, SIZE(molecule_list)
240 426643 : molecule_list(j) = j + counter
241 : END DO
242 128900 : molecule_kind => molecule_kind_set(ikind)
243 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
244 128900 : molecule_list=molecule_list)
245 128900 : IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:)
246 128900 : IF (found_last) EXIT
247 126155 : counter = imol
248 126155 : ikind = map_atom_type(i + 1)
249 : END IF
250 : END DO
251 : ! Treat separately the case in which the last atom is also a molecule
252 10274 : IF (i == natom) THEN
253 7529 : imol = ABS(map_atom_mol(natom))
254 : ! Last atom is also a molecule by itself
255 22587 : ALLOCATE (molecule_list(imol - counter))
256 15058 : DO j = 1, SIZE(molecule_list)
257 15058 : molecule_list(j) = j + counter
258 : END DO
259 7529 : molecule_kind => molecule_kind_set(ikind)
260 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
261 7529 : molecule_list=molecule_list)
262 7529 : IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:)
263 : END IF
264 10274 : CALL timestop(handle2)
265 : !-----------------------------------------------------------------------------
266 : !-----------------------------------------------------------------------------
267 : ! 6. Set the molecule_set(imol)%molecule_kind via set_molecule
268 : !-----------------------------------------------------------------------------
269 10274 : CALL timeset(routineN//"_6", handle2)
270 146703 : DO ikind = 1, SIZE(molecule_kind_set)
271 136429 : molecule_kind => molecule_kind_set(ikind)
272 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
273 136429 : molecule_list=molecule_list)
274 451975 : DO i = 1, SIZE(molecule_list)
275 305272 : molecule => molecule_set(molecule_list(i))
276 441701 : CALL set_molecule(molecule, molecule_kind=molecule_kind)
277 : END DO
278 : END DO
279 10274 : CALL timestop(handle2)
280 : !-----------------------------------------------------------------------------
281 : !-----------------------------------------------------------------------------
282 : ! 7. Set the molecule_set(imol)%[first_atom,last_atom] via set_molecule_set
283 : !-----------------------------------------------------------------------------
284 30822 : ALLOCATE (first_list(SIZE(molecule_set)))
285 20548 : ALLOCATE (last_list(SIZE(molecule_set)))
286 10274 : CALL timeset(routineN//"_7", handle2)
287 315546 : first_list(:) = 0
288 315546 : last_list(:) = 0
289 10274 : ityp = atom_info%map_mol_typ(1)
290 10274 : inum = atom_info%map_mol_num(1)
291 10274 : ires = atom_info%map_mol_res(1)
292 10274 : imol = 1
293 10274 : first_list(1) = 1
294 754569 : DO j = 2, natom
295 : IF ((atom_info%map_mol_typ(j) /= ityp) .OR. &
296 744295 : (atom_info%map_mol_num(j) /= inum) .OR. &
297 10274 : (atom_info%map_mol_res(j) /= ires)) THEN
298 294998 : ityp = atom_info%map_mol_typ(j)
299 294998 : inum = atom_info%map_mol_num(j)
300 294998 : ires = atom_info%map_mol_res(j)
301 294998 : imol = imol + 1
302 294998 : first_list(imol) = j
303 : END IF
304 : END DO
305 10274 : CPASSERT(imol == topology%nmol)
306 305272 : DO ikind = 1, topology%nmol - 1
307 305272 : last_list(ikind) = first_list(ikind + 1) - 1
308 : END DO
309 10274 : last_list(topology%nmol) = topology%natoms
310 10274 : CALL set_molecule_set(molecule_set, first_list, last_list)
311 10274 : CALL timestop(handle2)
312 : !-----------------------------------------------------------------------------
313 : !-----------------------------------------------------------------------------
314 : ! 8. Set and NULLIFY the molecule_set(imol)%lmi via set_molecule_set
315 : !-----------------------------------------------------------------------------
316 10274 : CALL timeset(routineN//"_8", handle2)
317 315546 : DO imol = 1, SIZE(molecule_set)
318 305272 : molecule => molecule_set(imol)
319 305272 : NULLIFY (lmi)
320 315546 : CALL set_molecule(molecule, lmi=lmi)
321 : END DO
322 10274 : CALL timestop(handle2)
323 : !-----------------------------------------------------------------------------
324 : !-----------------------------------------------------------------------------
325 : ! 9. Set the atom_list for molecule_kind in molecule_kind_set (PART 1)
326 : !-----------------------------------------------------------------------------
327 10274 : CALL timeset(routineN//"_9", handle2)
328 10274 : counter = 0
329 315546 : DO imol = 1, SIZE(molecule_set)
330 305272 : molecule => molecule_set(imol)
331 305272 : molecule_kind => molecule_set(imol)%molecule_kind
332 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
333 305272 : kind_number=i)
334 315546 : IF (counter /= i) THEN
335 136429 : counter = i
336 : CALL get_molecule(molecule=molecule, &
337 136429 : first_atom=first, last_atom=last)
338 136429 : natom = 0
339 136429 : IF (first /= 0 .AND. last /= 0) natom = last - first + 1
340 658607 : ALLOCATE (atom_list(natom))
341 385749 : DO i = 1, natom
342 : !Atomic kind information will be filled in (PART 2)
343 249320 : NULLIFY (atom_list(i)%atomic_kind)
344 249320 : atom_list(i)%id_name = atom_info%id_atmname(i + first - 1)
345 250110 : IF (iw > 0) WRITE (iw, '(5X,A,3I5,1X,A5)') "atom_list ", &
346 138009 : imol, counter, i, TRIM(id2str(atom_list(i)%id_name))
347 : END DO
348 136429 : CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
349 : END IF
350 : END DO
351 10274 : CALL timestop(handle2)
352 : !-----------------------------------------------------------------------------
353 : !-----------------------------------------------------------------------------
354 : ! 10. Set the molecule_kind%[nbond,bond_list] via set_molecule_kind
355 : !-----------------------------------------------------------------------------
356 10274 : CALL timeset(routineN//"_10", handle2)
357 : ! First map bonds on molecules
358 10274 : nvar1 = 0
359 10274 : nvar2 = 0
360 10274 : NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
361 10274 : IF (ASSOCIATED(conn_info%bond_a)) nvar1 = SIZE(conn_info%bond_a)
362 10274 : IF (ASSOCIATED(conn_info%c_bond_a)) nvar2 = SIZE(conn_info%c_bond_a)
363 10274 : nval_tot1 = nvar1
364 10274 : nval_tot2 = 0
365 23095 : ALLOCATE (map_var_mol(nvar1))
366 20702 : ALLOCATE (map_cvar_mol(nvar2))
367 518581 : map_var_mol = -1
368 13140 : map_cvar_mol = -1
369 518581 : DO i = 1, nvar1
370 508307 : j1 = map_atom_mol(conn_info%bond_a(i))
371 508307 : j2 = map_atom_mol(conn_info%bond_b(i))
372 518581 : IF (j1 == j2) THEN
373 505441 : IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%bond_a(i))
374 : END IF
375 : END DO
376 13140 : DO i = 1, nvar2
377 2866 : min_index = MIN(conn_info%c_bond_a(i), conn_info%c_bond_b(i))
378 2866 : j1 = map_atom_mol(min_index)
379 13140 : IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
380 : END DO
381 10274 : CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
382 10274 : CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
383 146703 : DO i = 1, topology%nmol_type
384 136429 : intra_bonds = 0
385 136429 : inter_bonds = 0
386 197055 : IF (ALL(bnd_type(:, i) > 0)) THEN
387 30313 : intra_bonds = bnd_type(2, i) - bnd_type(1, i) + 1
388 : END IF
389 142097 : IF (ALL(bnd_ctype(:, i) > 0)) THEN
390 2834 : inter_bonds = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
391 : END IF
392 136429 : ibond = intra_bonds + inter_bonds
393 136429 : IF (iw > 0) THEN
394 434 : WRITE (iw, *) " Total number bonds for molecule type ", i, " :", ibond
395 434 : WRITE (iw, *) " intra (bonds inside molecules) :: ", intra_bonds
396 434 : WRITE (iw, *) " inter (bonds between molecules) :: ", inter_bonds
397 : END IF
398 136429 : molecule_kind => molecule_kind_set(i)
399 136429 : nval_tot2 = nval_tot2 + ibond*SIZE(molecule_kind%molecule_list)
400 :
401 419468 : ALLOCATE (bond_list(ibond))
402 136429 : ibond = 0
403 355974 : DO j = bnd_type(1, i), bnd_type(2, i)
404 219545 : IF (j == 0) CYCLE
405 113429 : ibond = ibond + 1
406 113429 : jind = map_vars(j)
407 113429 : first = first_list(map_atom_mol(conn_info%bond_a(jind)))
408 113429 : bond_list(ibond)%a = conn_info%bond_a(jind) - first + 1
409 113429 : bond_list(ibond)%b = conn_info%bond_b(jind) - first + 1
410 : ! Set by default id_type to charmm and modify when handling the forcefield
411 113429 : bond_list(ibond)%id_type = do_ff_charmm
412 113429 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
413 132 : bond_list(ibond)%itype = conn_info%bond_type(jind)
414 : END IF
415 : !point this to the right bond_kind_type if using force field
416 113429 : NULLIFY (bond_list(ibond)%bond_kind)
417 249858 : IF (iw > 0) THEN
418 135 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", &
419 135 : i, " intra bond", &
420 135 : conn_info%bond_a(jind), &
421 135 : conn_info%bond_b(jind), &
422 135 : "offset number at", &
423 135 : conn_info%bond_a(jind) - first + 1, &
424 270 : conn_info%bond_b(jind) - first + 1
425 : END IF
426 : END DO
427 272890 : DO j = bnd_ctype(1, i), bnd_ctype(2, i)
428 136461 : IF (j == 0) CYCLE
429 2866 : ibond = ibond + 1
430 2866 : jind = map_cvars(j)
431 2866 : min_index = MIN(conn_info%c_bond_a(jind), conn_info%c_bond_b(jind))
432 2866 : first = first_list(map_atom_mol(min_index))
433 2866 : bond_list(ibond)%a = conn_info%c_bond_a(jind) - first + 1
434 2866 : bond_list(ibond)%b = conn_info%c_bond_b(jind) - first + 1
435 : ! Set by default id_type to charmm and modify when handling the forcefield
436 2866 : bond_list(ibond)%id_type = do_ff_charmm
437 2866 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
438 0 : bond_list(ibond)%itype = conn_info%c_bond_type(jind)
439 : END IF
440 : !point this to the right bond_kind_type if using force field
441 2866 : NULLIFY (bond_list(ibond)%bond_kind)
442 139295 : IF (iw > 0) THEN
443 2 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", &
444 2 : i, " inter bond", &
445 2 : conn_info%c_bond_a(jind), &
446 2 : conn_info%c_bond_b(jind), &
447 2 : "offset number at", &
448 2 : conn_info%c_bond_a(jind) - first + 1, &
449 4 : conn_info%c_bond_b(jind) - first + 1
450 : END IF
451 : END DO
452 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
453 146703 : nbond=SIZE(bond_list), bond_list=bond_list)
454 : END DO
455 10274 : IF ((nval_tot1 /= nval_tot2) .AND. (output_unit > 0)) THEN
456 0 : WRITE (output_unit, '(/)')
457 0 : WRITE (output_unit, '(T5,A)') "ERROR| Mismatching found between the total number of atoms"
458 0 : WRITE (output_unit, '(T5,A)') "ERROR| and the number of atoms computed multiplying the Nr."
459 0 : WRITE (output_unit, '(T5,A)') "ERROR| of molecules by the number of atoms building that"
460 0 : WRITE (output_unit, '(T5,A)') "ERROR| kind of molecule."
461 0 : WRITE (output_unit, '(T5,A)') "ERROR| This happens when the connectivity is wrongly built"
462 0 : WRITE (output_unit, '(T5,A)') "ERROR| One example could be two same kind of molecules have"
463 0 : WRITE (output_unit, '(T5,A)') "ERROR| a different number of atoms. Check the connectivity!"
464 : END IF
465 10274 : CPASSERT(nval_tot1 == nval_tot2)
466 10274 : DEALLOCATE (map_var_mol)
467 10274 : DEALLOCATE (map_cvar_mol)
468 10274 : DEALLOCATE (map_vars)
469 10274 : DEALLOCATE (map_cvars)
470 10274 : DEALLOCATE (bnd_type)
471 10274 : DEALLOCATE (bnd_ctype)
472 10274 : CALL timestop(handle2)
473 : !-----------------------------------------------------------------------------
474 : !-----------------------------------------------------------------------------
475 : ! 11. Set the molecule_kind%[nbend,bend_list] via set_molecule_kind
476 : !-----------------------------------------------------------------------------
477 : ! Allocate c_var_a, c_var_b, c_var_c, c_var_type
478 10274 : CALL timeset(routineN//"_11_pre", handle2)
479 10274 : idim = 0
480 10274 : ALLOCATE (c_var_a(idim))
481 10274 : ALLOCATE (c_var_b(idim))
482 10274 : ALLOCATE (c_var_c(idim))
483 10274 : found = ASSOCIATED(conn_info%theta_type)
484 10274 : IF (found) THEN
485 14 : ALLOCATE (c_var_type(idim))
486 : END IF
487 10274 : IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%theta_a)) THEN
488 230290 : DO j = 1, SIZE(conn_info%theta_a)
489 222343 : j1 = map_atom_mol(conn_info%theta_a(j))
490 222343 : j2 = map_atom_mol(conn_info%theta_b(j))
491 222343 : j3 = map_atom_mol(conn_info%theta_c(j))
492 230290 : IF (j1 /= j2 .OR. j2 /= j3) THEN
493 11392 : idim = idim + 1
494 : END IF
495 : END DO
496 7947 : CALL reallocate(c_var_a, 1, idim)
497 7947 : CALL reallocate(c_var_b, 1, idim)
498 7947 : CALL reallocate(c_var_c, 1, idim)
499 7947 : IF (found) THEN
500 0 : CALL reallocate(c_var_type, 1, idim)
501 : END IF
502 7947 : idim = 0
503 230290 : DO j = 1, SIZE(conn_info%theta_a)
504 222343 : j1 = map_atom_mol(conn_info%theta_a(j))
505 222343 : j2 = map_atom_mol(conn_info%theta_b(j))
506 222343 : j3 = map_atom_mol(conn_info%theta_c(j))
507 230290 : IF (j1 /= j2 .OR. j2 /= j3) THEN
508 11392 : idim = idim + 1
509 11392 : c_var_a(idim) = conn_info%theta_a(j)
510 11392 : c_var_b(idim) = conn_info%theta_b(j)
511 11392 : c_var_c(idim) = conn_info%theta_c(j)
512 11392 : IF (found) THEN
513 0 : c_var_type(idim) = conn_info%theta_type(j)
514 : END IF
515 : END IF
516 : END DO
517 : END IF
518 10274 : CALL timestop(handle2)
519 10274 : CALL timeset(routineN//"_11", handle2)
520 : ! map bends on molecules
521 10274 : nvar1 = 0
522 10274 : nvar2 = 0
523 : NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
524 10274 : IF (ASSOCIATED(conn_info%theta_a)) nvar1 = SIZE(conn_info%theta_a)
525 10274 : IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
526 10274 : nval_tot1 = nvar1
527 10274 : nval_tot2 = 0
528 22925 : ALLOCATE (map_var_mol(nvar1))
529 20702 : ALLOCATE (map_cvar_mol(nvar2))
530 281849 : map_var_mol = -1
531 21666 : map_cvar_mol = -1
532 281849 : DO i = 1, nvar1
533 271575 : j1 = map_atom_mol(conn_info%theta_a(i))
534 271575 : j2 = map_atom_mol(conn_info%theta_b(i))
535 271575 : j3 = map_atom_mol(conn_info%theta_c(i))
536 281849 : IF (j1 == j2 .AND. j2 == j3) THEN
537 260183 : IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%theta_a(i))
538 : END IF
539 : END DO
540 21666 : DO i = 1, nvar2
541 11392 : min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i))
542 11392 : j1 = map_atom_mol(min_index)
543 21666 : IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
544 : END DO
545 10274 : CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
546 10274 : CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
547 146703 : DO i = 1, topology%nmol_type
548 136429 : intra_bends = 0
549 136429 : inter_bends = 0
550 196043 : IF (ALL(bnd_type(:, i) > 0)) THEN
551 29807 : intra_bends = bnd_type(2, i) - bnd_type(1, i) + 1
552 : END IF
553 142097 : IF (ALL(bnd_ctype(:, i) > 0)) THEN
554 2834 : inter_bends = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
555 : END IF
556 136429 : ibend = intra_bends + inter_bends
557 136429 : IF (iw > 0) THEN
558 434 : WRITE (iw, *) " Total number of angles for molecule type ", i, " :", ibend
559 434 : WRITE (iw, *) " intra (angles inside molecules) :: ", intra_bends
560 434 : WRITE (iw, *) " inter (angles between molecules) :: ", inter_bends
561 : END IF
562 136429 : molecule_kind => molecule_kind_set(i)
563 136429 : nval_tot2 = nval_tot2 + ibend*SIZE(molecule_kind%molecule_list)
564 443459 : ALLOCATE (bend_list(ibend))
565 136429 : ibend = 0
566 372451 : DO j = bnd_type(1, i), bnd_type(2, i)
567 236022 : IF (j == 0) CYCLE
568 129400 : ibend = ibend + 1
569 129400 : jind = map_vars(j)
570 129400 : first = first_list(map_atom_mol(conn_info%theta_a(jind)))
571 129400 : bend_list(ibend)%a = conn_info%theta_a(jind) - first + 1
572 129400 : bend_list(ibend)%b = conn_info%theta_b(jind) - first + 1
573 129400 : bend_list(ibend)%c = conn_info%theta_c(jind) - first + 1
574 : ! Set by default id_type to charmm and modify when handling the forcefield
575 129400 : bend_list(ibend)%id_type = do_ff_charmm
576 129400 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
577 158 : bend_list(ibend)%itype = conn_info%theta_type(jind)
578 : END IF
579 : !point this to the right bend_kind_type if using force field
580 129400 : NULLIFY (bend_list(ibend)%bend_kind)
581 265829 : IF (iw > 0) THEN
582 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
583 229 : "molecule_kind", ikind, "intra bend", &
584 229 : conn_info%theta_a(jind), &
585 229 : conn_info%theta_b(jind), &
586 229 : conn_info%theta_c(jind), &
587 229 : "offset number at", &
588 229 : conn_info%theta_a(jind) - first + 1, &
589 229 : conn_info%theta_b(jind) - first + 1, &
590 458 : conn_info%theta_c(jind) - first + 1
591 : END IF
592 : END DO
593 281416 : DO j = bnd_ctype(1, i), bnd_ctype(2, i)
594 144987 : IF (j == 0) CYCLE
595 11392 : ibend = ibend + 1
596 11392 : jind = map_cvars(j)
597 11392 : min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind))
598 11392 : first = first_list(map_atom_mol(min_index))
599 11392 : bend_list(ibend)%a = c_var_a(jind) - first + 1
600 11392 : bend_list(ibend)%b = c_var_b(jind) - first + 1
601 11392 : bend_list(ibend)%c = c_var_c(jind) - first + 1
602 : ! Set by default id_type to charmm and modify when handling the forcefield
603 11392 : bend_list(ibend)%id_type = do_ff_charmm
604 11392 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
605 0 : bend_list(ibend)%itype = c_var_type(jind)
606 : END IF
607 : !point this to the right bend_kind_type if using force field
608 11392 : NULLIFY (bend_list(ibend)%bend_kind)
609 147821 : IF (iw > 0) THEN
610 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
611 8 : "molecule_kind", ikind, "inter bend", &
612 8 : c_var_a(jind), &
613 8 : c_var_b(jind), &
614 8 : c_var_c(jind), &
615 8 : "offset number at", &
616 8 : c_var_a(jind) - first + 1, &
617 8 : c_var_b(jind) - first + 1, &
618 16 : c_var_c(jind) - first + 1
619 : END IF
620 : END DO
621 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
622 146703 : nbend=SIZE(bend_list), bend_list=bend_list)
623 : END DO
624 10274 : CPASSERT(nval_tot1 == nval_tot2)
625 10274 : DEALLOCATE (map_var_mol)
626 10274 : DEALLOCATE (map_cvar_mol)
627 10274 : DEALLOCATE (map_vars)
628 10274 : DEALLOCATE (map_cvars)
629 10274 : DEALLOCATE (bnd_type)
630 10274 : DEALLOCATE (bnd_ctype)
631 10274 : DEALLOCATE (c_var_a)
632 10274 : DEALLOCATE (c_var_b)
633 10274 : DEALLOCATE (c_var_c)
634 10274 : IF (found) THEN
635 14 : DEALLOCATE (c_var_type)
636 : END IF
637 10274 : CALL timestop(handle2)
638 : !-----------------------------------------------------------------------------
639 : !-----------------------------------------------------------------------------
640 : ! 12. Set the molecule_kind%[nub,ub_list] via set_molecule_kind
641 : !-----------------------------------------------------------------------------
642 10274 : CALL timeset(routineN//"_12_pre", handle2)
643 10274 : idim = 0
644 10274 : ALLOCATE (c_var_a(idim))
645 10274 : ALLOCATE (c_var_b(idim))
646 10274 : ALLOCATE (c_var_c(idim))
647 10274 : IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%ub_a)) THEN
648 230290 : DO j = 1, SIZE(conn_info%ub_a)
649 222343 : j1 = map_atom_mol(conn_info%ub_a(j))
650 222343 : j2 = map_atom_mol(conn_info%ub_b(j))
651 222343 : j3 = map_atom_mol(conn_info%ub_c(j))
652 230290 : IF (j1 /= j2 .OR. j2 /= j3) THEN
653 11392 : idim = idim + 1
654 : END IF
655 : END DO
656 7947 : CALL reallocate(c_var_a, 1, idim)
657 7947 : CALL reallocate(c_var_b, 1, idim)
658 7947 : CALL reallocate(c_var_c, 1, idim)
659 7947 : idim = 0
660 230290 : DO j = 1, SIZE(conn_info%ub_a)
661 222343 : j1 = map_atom_mol(conn_info%ub_a(j))
662 222343 : j2 = map_atom_mol(conn_info%ub_b(j))
663 222343 : j3 = map_atom_mol(conn_info%ub_c(j))
664 230290 : IF (j1 /= j2 .OR. j2 /= j3) THEN
665 11392 : idim = idim + 1
666 11392 : c_var_a(idim) = conn_info%ub_a(j)
667 11392 : c_var_b(idim) = conn_info%ub_b(j)
668 11392 : c_var_c(idim) = conn_info%ub_c(j)
669 : END IF
670 : END DO
671 : END IF
672 10274 : CALL timestop(handle2)
673 10274 : CALL timeset(routineN//"_12", handle2)
674 : ! map UBs on molecules
675 10274 : nvar1 = 0
676 10274 : nvar2 = 0
677 : NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
678 10274 : IF (ASSOCIATED(conn_info%ub_a)) nvar1 = SIZE(conn_info%ub_a)
679 10274 : IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
680 10274 : nval_tot1 = nvar1
681 10274 : nval_tot2 = 0
682 22911 : ALLOCATE (map_var_mol(nvar1))
683 20702 : ALLOCATE (map_cvar_mol(nvar2))
684 281691 : map_var_mol = -1
685 21666 : map_cvar_mol = -1
686 281691 : DO i = 1, nvar1
687 271417 : j1 = map_atom_mol(conn_info%ub_a(i))
688 271417 : j2 = map_atom_mol(conn_info%ub_b(i))
689 271417 : j3 = map_atom_mol(conn_info%ub_c(i))
690 281691 : IF (j1 == j2 .AND. j2 == j3) THEN
691 260025 : IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%ub_a(i))
692 : END IF
693 : END DO
694 21666 : DO i = 1, nvar2
695 11392 : min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i))
696 11392 : j1 = map_atom_mol(min_index)
697 21666 : IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
698 : END DO
699 10274 : CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
700 10274 : CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
701 146703 : DO i = 1, topology%nmol_type
702 136429 : intra_ubs = 0
703 136429 : inter_ubs = 0
704 196015 : IF (ALL(bnd_type(:, i) > 0)) THEN
705 29793 : intra_ubs = bnd_type(2, i) - bnd_type(1, i) + 1
706 : END IF
707 142097 : IF (ALL(bnd_ctype(:, i) > 0)) THEN
708 2834 : inter_ubs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
709 : END IF
710 136429 : iub = intra_ubs + inter_ubs
711 136429 : IF (iw > 0) THEN
712 434 : WRITE (iw, *) " Total number of Urey-Bradley for molecule type ", i, " :", iub
713 434 : WRITE (iw, *) " intra (UB inside molecules) :: ", intra_ubs
714 434 : WRITE (iw, *) " inter (UB between molecules) :: ", inter_ubs
715 : END IF
716 136429 : molecule_kind => molecule_kind_set(i)
717 136429 : nval_tot2 = nval_tot2 + iub*SIZE(molecule_kind%molecule_list)
718 443287 : ALLOCATE (ub_list(iub))
719 136429 : iub = 0
720 372307 : DO j = bnd_type(1, i), bnd_type(2, i)
721 235878 : IF (j == 0) CYCLE
722 129242 : iub = iub + 1
723 129242 : jind = map_vars(j)
724 129242 : first = first_list(map_atom_mol(conn_info%ub_a(jind)))
725 129242 : ub_list(iub)%a = conn_info%ub_a(jind) - first + 1
726 129242 : ub_list(iub)%b = conn_info%ub_b(jind) - first + 1
727 129242 : ub_list(iub)%c = conn_info%ub_c(jind) - first + 1
728 129242 : ub_list(iub)%id_type = do_ff_charmm
729 : !point this to the right ub_kind_type if using force field
730 129242 : NULLIFY (ub_list(iub)%ub_kind)
731 265671 : IF (iw > 0) THEN
732 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
733 229 : "molecule_kind", i, "intra UB", &
734 229 : conn_info%ub_a(jind), &
735 229 : conn_info%ub_b(jind), &
736 229 : conn_info%ub_c(jind), &
737 229 : "offset number at", &
738 229 : conn_info%ub_a(jind) - first + 1, &
739 229 : conn_info%ub_b(jind) - first + 1, &
740 458 : conn_info%ub_c(jind) - first + 1
741 : END IF
742 : END DO
743 281416 : DO j = bnd_ctype(1, i), bnd_ctype(2, i)
744 144987 : IF (j == 0) CYCLE
745 11392 : iub = iub + 1
746 11392 : jind = map_cvars(j)
747 11392 : min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind))
748 11392 : first = first_list(map_atom_mol(min_index))
749 11392 : ub_list(iub)%a = c_var_a(jind) - first + 1
750 11392 : ub_list(iub)%b = c_var_b(jind) - first + 1
751 11392 : ub_list(iub)%c = c_var_c(jind) - first + 1
752 11392 : ub_list(iub)%id_type = do_ff_charmm
753 : !point this to the right ub_kind_type if using force field
754 11392 : NULLIFY (ub_list(iub)%ub_kind)
755 147821 : IF (iw > 0) THEN
756 : WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
757 8 : "molecule_kind", i, "inter UB", &
758 8 : c_var_a(jind), &
759 8 : c_var_b(jind), &
760 8 : c_var_c(jind), &
761 8 : "offset number at", &
762 8 : c_var_a(jind) - first + 1, &
763 8 : c_var_b(jind) - first + 1, &
764 16 : c_var_c(jind) - first + 1
765 : END IF
766 : END DO
767 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
768 146703 : nub=SIZE(ub_list), ub_list=ub_list)
769 : END DO
770 10274 : CPASSERT(nval_tot1 == nval_tot2)
771 10274 : DEALLOCATE (map_var_mol)
772 10274 : DEALLOCATE (map_cvar_mol)
773 10274 : DEALLOCATE (map_vars)
774 10274 : DEALLOCATE (map_cvars)
775 10274 : DEALLOCATE (bnd_type)
776 10274 : DEALLOCATE (bnd_ctype)
777 10274 : DEALLOCATE (c_var_a)
778 10274 : DEALLOCATE (c_var_b)
779 10274 : DEALLOCATE (c_var_c)
780 10274 : CALL timestop(handle2)
781 : !-----------------------------------------------------------------------------
782 : !-----------------------------------------------------------------------------
783 : ! 13. Set the molecule_kind%[ntorsion,torsion_list] via set_molecule_kind
784 : !-----------------------------------------------------------------------------
785 : ! Allocate c_var_a, c_var_b, c_var_c, c_var_d, c_var_type
786 10274 : CALL timeset(routineN//"_13_pre", handle2)
787 10274 : idim = 0
788 10274 : ALLOCATE (c_var_a(idim))
789 10274 : ALLOCATE (c_var_b(idim))
790 10274 : ALLOCATE (c_var_c(idim))
791 10274 : ALLOCATE (c_var_d(idim))
792 10274 : found = ASSOCIATED(conn_info%phi_type)
793 10274 : IF (found) THEN
794 14 : ALLOCATE (c_var_type(idim))
795 : END IF
796 10274 : IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%phi_a)) THEN
797 142338 : DO j = 1, SIZE(conn_info%phi_a)
798 134391 : j1 = map_atom_mol(conn_info%phi_a(j))
799 134391 : j2 = map_atom_mol(conn_info%phi_b(j))
800 134391 : j3 = map_atom_mol(conn_info%phi_c(j))
801 134391 : j4 = map_atom_mol(conn_info%phi_d(j))
802 142338 : IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
803 31630 : idim = idim + 1
804 : END IF
805 : END DO
806 7947 : CALL reallocate(c_var_a, 1, idim)
807 7947 : CALL reallocate(c_var_b, 1, idim)
808 7947 : CALL reallocate(c_var_c, 1, idim)
809 7947 : CALL reallocate(c_var_d, 1, idim)
810 7947 : IF (found) THEN
811 0 : CALL reallocate(c_var_type, 1, idim)
812 : END IF
813 7947 : idim = 0
814 142338 : DO j = 1, SIZE(conn_info%phi_a)
815 134391 : j1 = map_atom_mol(conn_info%phi_a(j))
816 134391 : j2 = map_atom_mol(conn_info%phi_b(j))
817 134391 : j3 = map_atom_mol(conn_info%phi_c(j))
818 134391 : j4 = map_atom_mol(conn_info%phi_d(j))
819 142338 : IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
820 31630 : idim = idim + 1
821 31630 : c_var_a(idim) = conn_info%phi_a(j)
822 31630 : c_var_b(idim) = conn_info%phi_b(j)
823 31630 : c_var_c(idim) = conn_info%phi_c(j)
824 31630 : c_var_d(idim) = conn_info%phi_d(j)
825 31630 : IF (found) THEN
826 0 : c_var_type(idim) = conn_info%phi_type(j)
827 : END IF
828 : END IF
829 : END DO
830 : END IF
831 10274 : CALL timestop(handle2)
832 10274 : CALL timeset(routineN//"_13", handle2)
833 : ! map torsions on molecules
834 10274 : nvar1 = 0
835 10274 : nvar2 = 0
836 : NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
837 10274 : IF (ASSOCIATED(conn_info%phi_a)) nvar1 = SIZE(conn_info%phi_a)
838 10274 : IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
839 10274 : nval_tot1 = nvar1
840 10274 : nval_tot2 = 0
841 21202 : ALLOCATE (map_var_mol(nvar1))
842 20700 : ALLOCATE (map_cvar_mol(nvar2))
843 175743 : map_var_mol = -1
844 41904 : map_cvar_mol = -1
845 175743 : DO i = 1, nvar1
846 165469 : j1 = map_atom_mol(conn_info%phi_a(i))
847 165469 : j2 = map_atom_mol(conn_info%phi_b(i))
848 165469 : j3 = map_atom_mol(conn_info%phi_c(i))
849 165469 : j4 = map_atom_mol(conn_info%phi_d(i))
850 175743 : IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN
851 133839 : IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%phi_a(i))
852 : END IF
853 : END DO
854 41904 : DO i = 1, nvar2
855 31630 : min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i))
856 31630 : j1 = map_atom_mol(min_index)
857 41904 : IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
858 : END DO
859 10274 : CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
860 10274 : CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
861 146703 : DO i = 1, topology%nmol_type
862 136429 : intra_torsions = 0
863 136429 : inter_torsions = 0
864 147669 : IF (ALL(bnd_type(:, i) > 0)) THEN
865 5620 : intra_torsions = bnd_type(2, i) - bnd_type(1, i) + 1
866 : END IF
867 142093 : IF (ALL(bnd_ctype(:, i) > 0)) THEN
868 2832 : inter_torsions = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
869 : END IF
870 136429 : itorsion = intra_torsions + inter_torsions
871 136429 : IF (iw > 0) THEN
872 434 : WRITE (iw, *) " Total number of torsions for molecule type ", i, " :", itorsion
873 434 : WRITE (iw, *) " intra (torsions inside molecules) :: ", intra_torsions
874 434 : WRITE (iw, *) " inter (torsions between molecules) :: ", inter_torsions
875 : END IF
876 136429 : molecule_kind => molecule_kind_set(i)
877 136429 : nval_tot2 = nval_tot2 + itorsion*SIZE(molecule_kind%molecule_list)
878 435099 : ALLOCATE (torsion_list(itorsion))
879 136429 : itorsion = 0
880 392229 : DO j = bnd_type(1, i), bnd_type(2, i)
881 255800 : IF (j == 0) CYCLE
882 124991 : itorsion = itorsion + 1
883 124991 : jind = map_vars(j)
884 124991 : first = first_list(map_atom_mol(conn_info%phi_a(jind)))
885 124991 : torsion_list(itorsion)%a = conn_info%phi_a(jind) - first + 1
886 124991 : torsion_list(itorsion)%b = conn_info%phi_b(jind) - first + 1
887 124991 : torsion_list(itorsion)%c = conn_info%phi_c(jind) - first + 1
888 124991 : torsion_list(itorsion)%d = conn_info%phi_d(jind) - first + 1
889 : ! Set by default id_type to charmm and modify when handling the forcefield
890 124991 : torsion_list(itorsion)%id_type = do_ff_charmm
891 124991 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
892 312 : torsion_list(itorsion)%itype = conn_info%phi_type(jind)
893 : END IF
894 : !point this to the right torsion_kind_type if using force field
895 124991 : NULLIFY (torsion_list(itorsion)%torsion_kind)
896 261420 : IF (iw > 0) THEN
897 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
898 366 : "molecule_kind", i, "intra TOR", &
899 366 : conn_info%phi_a(jind), &
900 366 : conn_info%phi_b(jind), &
901 366 : conn_info%phi_c(jind), &
902 366 : conn_info%phi_d(jind), &
903 366 : "offset number at", &
904 366 : conn_info%phi_a(jind) - first + 1, &
905 366 : conn_info%phi_b(jind) - first + 1, &
906 366 : conn_info%phi_c(jind) - first + 1, &
907 732 : conn_info%phi_d(jind) - first + 1
908 : END IF
909 : END DO
910 301656 : DO j = bnd_ctype(1, i), bnd_ctype(2, i)
911 165227 : IF (j == 0) CYCLE
912 31630 : itorsion = itorsion + 1
913 31630 : jind = map_cvars(j)
914 31630 : min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind))
915 31630 : first = first_list(map_atom_mol(min_index))
916 31630 : torsion_list(itorsion)%a = c_var_a(jind) - first + 1
917 31630 : torsion_list(itorsion)%b = c_var_b(jind) - first + 1
918 31630 : torsion_list(itorsion)%c = c_var_c(jind) - first + 1
919 31630 : torsion_list(itorsion)%d = c_var_d(jind) - first + 1
920 : ! Set by default id_type to charmm and modify when handling the forcefield
921 31630 : torsion_list(itorsion)%id_type = do_ff_charmm
922 31630 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
923 0 : torsion_list(itorsion)%itype = c_var_type(jind)
924 : END IF
925 : !point this to the right torsion_kind_type if using force field
926 31630 : NULLIFY (torsion_list(itorsion)%torsion_kind)
927 168059 : IF (iw > 0) THEN
928 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
929 34 : "molecule_kind", i, "inter TOR", &
930 34 : c_var_a(jind), &
931 34 : c_var_b(jind), &
932 34 : c_var_c(jind), &
933 34 : c_var_d(jind), &
934 34 : "offset number at", &
935 34 : c_var_a(jind) - first + 1, &
936 34 : c_var_b(jind) - first + 1, &
937 34 : c_var_c(jind) - first + 1, &
938 68 : c_var_d(jind) - first + 1
939 : END IF
940 : END DO
941 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
942 146703 : ntorsion=SIZE(torsion_list), torsion_list=torsion_list)
943 : END DO
944 10274 : CPASSERT(nval_tot1 == nval_tot2)
945 10274 : DEALLOCATE (map_var_mol)
946 10274 : DEALLOCATE (map_cvar_mol)
947 10274 : DEALLOCATE (map_vars)
948 10274 : DEALLOCATE (map_cvars)
949 10274 : DEALLOCATE (bnd_type)
950 10274 : DEALLOCATE (bnd_ctype)
951 10274 : DEALLOCATE (c_var_a)
952 10274 : DEALLOCATE (c_var_b)
953 10274 : DEALLOCATE (c_var_c)
954 10274 : DEALLOCATE (c_var_d)
955 10274 : IF (found) THEN
956 14 : DEALLOCATE (c_var_type)
957 : END IF
958 10274 : CALL timestop(handle2)
959 : !-----------------------------------------------------------------------------
960 : !-----------------------------------------------------------------------------
961 : ! 14. Set the molecule_kind%[nimpr,impr_list] via set_molecule_kind
962 : ! Also set the molecule_kind%[nopbend,opbend_list]
963 : !-----------------------------------------------------------------------------
964 : ! Allocate c_var_a, c_var_b, c_var_c, c_var_d, c_var_type
965 10274 : CALL timeset(routineN//"_14_pre", handle2)
966 10274 : idim = 0
967 10274 : ALLOCATE (c_var_a(idim))
968 10274 : ALLOCATE (c_var_b(idim))
969 10274 : ALLOCATE (c_var_c(idim))
970 10274 : ALLOCATE (c_var_d(idim))
971 10274 : found = ASSOCIATED(conn_info%impr_type)
972 10274 : IF (found) THEN
973 14 : ALLOCATE (c_var_type(idim))
974 : END IF
975 10274 : IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%impr_a)) THEN
976 12361 : DO j = 1, SIZE(conn_info%impr_a)
977 4414 : j1 = map_atom_mol(conn_info%impr_a(j))
978 4414 : j2 = map_atom_mol(conn_info%impr_b(j))
979 4414 : j3 = map_atom_mol(conn_info%impr_c(j))
980 4414 : j4 = map_atom_mol(conn_info%impr_d(j))
981 12361 : IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
982 3124 : idim = idim + 1
983 : END IF
984 : END DO
985 7947 : CALL reallocate(c_var_a, 1, idim)
986 7947 : CALL reallocate(c_var_b, 1, idim)
987 7947 : CALL reallocate(c_var_c, 1, idim)
988 7947 : CALL reallocate(c_var_d, 1, idim)
989 7947 : IF (found) THEN
990 0 : CALL reallocate(c_var_type, 1, idim)
991 : END IF
992 7947 : idim = 0
993 12361 : DO j = 1, SIZE(conn_info%impr_a)
994 4414 : j1 = map_atom_mol(conn_info%impr_a(j))
995 4414 : j2 = map_atom_mol(conn_info%impr_b(j))
996 4414 : j3 = map_atom_mol(conn_info%impr_c(j))
997 4414 : j4 = map_atom_mol(conn_info%impr_d(j))
998 12361 : IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
999 3124 : idim = idim + 1
1000 3124 : c_var_a(idim) = conn_info%impr_a(j)
1001 3124 : c_var_b(idim) = conn_info%impr_b(j)
1002 3124 : c_var_c(idim) = conn_info%impr_c(j)
1003 3124 : c_var_d(idim) = conn_info%impr_d(j)
1004 3124 : IF (found) THEN
1005 0 : c_var_type(idim) = conn_info%impr_type(j)
1006 : END IF
1007 : END IF
1008 : END DO
1009 : END IF
1010 10274 : CALL timestop(handle2)
1011 10274 : CALL timeset(routineN//"_14", handle2)
1012 : ! map imprs on molecules
1013 10274 : nvar1 = 0
1014 10274 : nvar2 = 0
1015 : NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
1016 10274 : IF (ASSOCIATED(conn_info%impr_a)) nvar1 = SIZE(conn_info%impr_a)
1017 10274 : IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
1018 10274 : nval_tot1 = nvar1
1019 10274 : nval_tot2 = 0
1020 20694 : ALLOCATE (map_var_mol(nvar1))
1021 20588 : ALLOCATE (map_cvar_mol(nvar2))
1022 15684 : map_var_mol = -1
1023 13398 : map_cvar_mol = -1
1024 15684 : DO i = 1, nvar1
1025 5410 : j1 = map_atom_mol(conn_info%impr_a(i))
1026 5410 : j2 = map_atom_mol(conn_info%impr_b(i))
1027 5410 : j3 = map_atom_mol(conn_info%impr_c(i))
1028 5410 : j4 = map_atom_mol(conn_info%impr_d(i))
1029 15684 : IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN
1030 2286 : IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%impr_a(i))
1031 : END IF
1032 : END DO
1033 13398 : DO i = 1, nvar2
1034 3124 : min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i))
1035 3124 : j1 = map_atom_mol(min_index)
1036 13398 : IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
1037 : END DO
1038 10274 : CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
1039 10274 : CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
1040 146703 : DO i = 1, topology%nmol_type
1041 136429 : intra_imprs = 0
1042 136429 : inter_imprs = 0
1043 137517 : IF (ALL(bnd_type(:, i) > 0)) THEN
1044 544 : intra_imprs = bnd_type(2, i) - bnd_type(1, i) + 1
1045 : END IF
1046 139553 : IF (ALL(bnd_ctype(:, i) > 0)) THEN
1047 1562 : inter_imprs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
1048 : END IF
1049 136429 : iimpr = intra_imprs + inter_imprs
1050 136429 : IF (iw > 0) THEN
1051 434 : WRITE (iw, *) " Total number of imprs for molecule type ", i, " :", iimpr
1052 434 : WRITE (iw, *) " intra (imprs inside molecules) :: ", intra_imprs
1053 434 : WRITE (iw, *) " inter (imprs between molecules) :: ", inter_imprs
1054 434 : WRITE (iw, *) " Total number of opbends for molecule type ", i, " :", iimpr
1055 434 : WRITE (iw, *) " intra (opbends inside molecules) :: ", intra_imprs
1056 434 : WRITE (iw, *) " inter (opbends between molecules) :: ", inter_imprs
1057 : END IF
1058 136429 : molecule_kind => molecule_kind_set(i)
1059 136429 : nval_tot2 = nval_tot2 + iimpr*SIZE(molecule_kind%molecule_list)
1060 279930 : ALLOCATE (impr_list(iimpr), STAT=stat)
1061 143501 : ALLOCATE (opbend_list(iimpr), STAT=stat)
1062 136429 : CPASSERT(stat == 0)
1063 136429 : iimpr = 0
1064 274576 : DO j = bnd_type(1, i), bnd_type(2, i)
1065 138147 : IF (j == 0) CYCLE
1066 2262 : iimpr = iimpr + 1
1067 2262 : jind = map_vars(j)
1068 2262 : first = first_list(map_atom_mol(conn_info%impr_a(jind)))
1069 2262 : impr_list(iimpr)%a = conn_info%impr_a(jind) - first + 1
1070 2262 : impr_list(iimpr)%b = conn_info%impr_b(jind) - first + 1
1071 2262 : impr_list(iimpr)%c = conn_info%impr_c(jind) - first + 1
1072 2262 : impr_list(iimpr)%d = conn_info%impr_d(jind) - first + 1
1073 : ! Atom sequence for improper is A B C D in which A is central atom,
1074 : ! B is deviating atom and C & D are secondairy atoms. Atom sequence for
1075 : ! opbend is B D C A in which A is central atom, B is deviating. Hence
1076 : ! to create an opbend out of an improper, B and D need to be interchanged.
1077 2262 : opbend_list(iimpr)%a = conn_info%impr_b(jind) - first + 1
1078 2262 : opbend_list(iimpr)%b = conn_info%impr_d(jind) - first + 1
1079 2262 : opbend_list(iimpr)%c = conn_info%impr_c(jind) - first + 1
1080 2262 : opbend_list(iimpr)%d = conn_info%impr_a(jind) - first + 1
1081 : ! Set by default id_type of improper to charmm and modify when handling the forcefield
1082 2262 : impr_list(iimpr)%id_type = do_ff_charmm
1083 : ! Set by default id_type of opbend to harmonic and modify when handling the forcefield
1084 2262 : opbend_list(iimpr)%id_type = do_ff_harmonic
1085 2262 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
1086 0 : impr_list(iimpr)%itype = conn_info%impr_type(jind)
1087 : END IF
1088 : !point this to the right impr_kind_type if using force field
1089 2262 : NULLIFY (impr_list(iimpr)%impr_kind)
1090 2262 : NULLIFY (opbend_list(iimpr)%opbend_kind)
1091 138691 : IF (iw > 0) THEN
1092 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1093 0 : "molecule_kind", i, "intra IMPR", &
1094 0 : conn_info%impr_a(jind), &
1095 0 : conn_info%impr_b(jind), &
1096 0 : conn_info%impr_c(jind), &
1097 0 : conn_info%impr_d(jind), &
1098 0 : "offset number at", &
1099 0 : conn_info%impr_a(jind) - first + 1, &
1100 0 : conn_info%impr_b(jind) - first + 1, &
1101 0 : conn_info%impr_c(jind) - first + 1, &
1102 0 : conn_info%impr_d(jind) - first + 1
1103 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1104 0 : "molecule_kind", i, "intra OPBEND", &
1105 0 : conn_info%impr_b(jind), &
1106 0 : conn_info%impr_d(jind), &
1107 0 : conn_info%impr_c(jind), &
1108 0 : conn_info%impr_a(jind), &
1109 0 : "offset number at", &
1110 0 : conn_info%impr_b(jind) - first + 1, &
1111 0 : conn_info%impr_d(jind) - first + 1, &
1112 0 : conn_info%impr_c(jind) - first + 1, &
1113 0 : conn_info%impr_a(jind) - first + 1
1114 : END IF
1115 : END DO
1116 274420 : DO j = bnd_ctype(1, i), bnd_ctype(2, i)
1117 137991 : IF (j == 0) CYCLE
1118 3124 : iimpr = iimpr + 1
1119 3124 : jind = map_cvars(j)
1120 3124 : min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind))
1121 3124 : first = first_list(map_atom_mol(min_index))
1122 3124 : impr_list(iimpr)%a = c_var_a(jind) - first + 1
1123 3124 : impr_list(iimpr)%b = c_var_b(jind) - first + 1
1124 3124 : impr_list(iimpr)%c = c_var_c(jind) - first + 1
1125 3124 : impr_list(iimpr)%d = c_var_d(jind) - first + 1
1126 3124 : opbend_list(iimpr)%a = c_var_b(jind) - first + 1
1127 3124 : opbend_list(iimpr)%b = c_var_d(jind) - first + 1
1128 3124 : opbend_list(iimpr)%c = c_var_c(jind) - first + 1
1129 3124 : opbend_list(iimpr)%d = c_var_a(jind) - first + 1
1130 : ! Set by default id_type of improper to charmm and modify when handling the forcefield
1131 3124 : impr_list(iimpr)%id_type = do_ff_charmm
1132 : ! Set by default id_type of opbend to harmonic and modify when handling the forcefield
1133 3124 : opbend_list(iimpr)%id_type = do_ff_harmonic
1134 3124 : IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
1135 0 : impr_list(iimpr)%itype = c_var_type(jind)
1136 : END IF
1137 : !point this to the right impr_kind_type and opbend_kind_type if using force field
1138 3124 : NULLIFY (impr_list(iimpr)%impr_kind)
1139 3124 : NULLIFY (opbend_list(iimpr)%opbend_kind)
1140 139553 : IF (iw > 0) THEN
1141 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1142 0 : "molecule_kind", i, "inter IMPR", &
1143 0 : c_var_a(jind), &
1144 0 : c_var_b(jind), &
1145 0 : c_var_c(jind), &
1146 0 : c_var_d(jind), &
1147 0 : "offset number at", &
1148 0 : c_var_a(jind) - first + 1, &
1149 0 : c_var_b(jind) - first + 1, &
1150 0 : c_var_c(jind) - first + 1, &
1151 0 : c_var_d(jind) - first + 1
1152 : WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1153 0 : "molecule_kind", i, "inter OPBEND", &
1154 0 : c_var_b(jind), &
1155 0 : c_var_d(jind), &
1156 0 : c_var_c(jind), &
1157 0 : c_var_a(jind), &
1158 0 : "offset number at", &
1159 0 : c_var_b(jind) - first + 1, &
1160 0 : c_var_d(jind) - first + 1, &
1161 0 : c_var_c(jind) - first + 1, &
1162 0 : c_var_a(jind) - first + 1
1163 : END IF
1164 : END DO
1165 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1166 136429 : nimpr=SIZE(impr_list), impr_list=impr_list)
1167 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1168 146703 : nopbend=SIZE(opbend_list), opbend_list=opbend_list)
1169 : END DO
1170 10274 : CPASSERT(nval_tot1 == nval_tot2)
1171 10274 : DEALLOCATE (map_var_mol)
1172 10274 : DEALLOCATE (map_cvar_mol)
1173 10274 : DEALLOCATE (map_vars)
1174 10274 : DEALLOCATE (map_cvars)
1175 10274 : DEALLOCATE (bnd_type)
1176 10274 : DEALLOCATE (bnd_ctype)
1177 10274 : DEALLOCATE (c_var_a)
1178 10274 : DEALLOCATE (c_var_b)
1179 10274 : DEALLOCATE (c_var_c)
1180 10274 : DEALLOCATE (c_var_d)
1181 10274 : IF (found) THEN
1182 14 : DEALLOCATE (c_var_type)
1183 : END IF
1184 10274 : CALL timestop(handle2)
1185 : !-----------------------------------------------------------------------------
1186 : !-----------------------------------------------------------------------------
1187 : ! Finally deallocate some stuff.
1188 : !-----------------------------------------------------------------------------
1189 10274 : DEALLOCATE (first_list)
1190 10274 : DEALLOCATE (last_list)
1191 10274 : DEALLOCATE (map_atom_mol)
1192 10274 : DEALLOCATE (map_atom_type)
1193 10274 : CALL timestop(handle)
1194 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1195 10274 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
1196 215754 : END SUBROUTINE topology_connectivity_pack
1197 :
1198 : ! **************************************************************************************************
1199 : !> \brief used to achieve linear scaling in the connectivity_pack
1200 : !> \param nmol_type ...
1201 : !> \param map_vars ...
1202 : !> \param map_var_mol ...
1203 : !> \param bnd_type ...
1204 : !> \param nvar1 ...
1205 : !> \par History
1206 : !> Teodoro Laino
1207 : ! **************************************************************************************************
1208 102740 : SUBROUTINE find_bnd_typ(nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
1209 : INTEGER, INTENT(IN) :: nmol_type
1210 : INTEGER, DIMENSION(:), POINTER :: map_vars, map_var_mol
1211 : INTEGER, DIMENSION(:, :), POINTER :: bnd_type
1212 : INTEGER, INTENT(IN) :: nvar1
1213 :
1214 : INTEGER :: i, ibond, j
1215 :
1216 214221 : ALLOCATE (map_vars(nvar1))
1217 102740 : CALL sort(map_var_mol, nvar1, map_vars)
1218 308220 : ALLOCATE (bnd_type(2, nmol_type))
1219 4195610 : bnd_type = 0
1220 102740 : IF (nvar1 == 0) RETURN
1221 731595 : DO j = 1, nvar1
1222 731595 : IF (map_var_mol(j) /= -1) EXIT
1223 : END DO
1224 8741 : IF (j == nvar1 + 1) RETURN
1225 8709 : i = map_var_mol(j)
1226 8709 : bnd_type(1, i) = j
1227 568437 : DO ibond = j, nvar1
1228 568437 : IF (map_var_mol(ibond) /= i) THEN
1229 100264 : bnd_type(2, i) = ibond - 1
1230 100264 : i = map_var_mol(ibond)
1231 100264 : bnd_type(1, i) = ibond
1232 : END IF
1233 : END DO
1234 8709 : bnd_type(2, i) = nvar1
1235 :
1236 : END SUBROUTINE find_bnd_typ
1237 :
1238 : ! **************************************************************************************************
1239 : !> \brief Handles the multiple unit cell option for the connectivity
1240 : !> \param topology ...
1241 : !> \param subsys_section ...
1242 : !> \author Teodoro Laino [tlaino] - 06.2009
1243 : ! **************************************************************************************************
1244 10274 : SUBROUTINE topology_conn_multiple(topology, subsys_section)
1245 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
1246 : TYPE(section_vals_type), POINTER :: subsys_section
1247 :
1248 : INTEGER :: a, fac, i, ind, j, k, m, natoms_orig, &
1249 : nbond, nbond_c, nimpr, nonfo, nphi, &
1250 : ntheta, nub
1251 10274 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
1252 : TYPE(connectivity_info_type), POINTER :: conn_info
1253 :
1254 10274 : NULLIFY (multiple_unit_cell)
1255 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
1256 10274 : i_vals=multiple_unit_cell)
1257 40662 : IF (ANY(multiple_unit_cell /= 1)) THEN
1258 592 : fac = PRODUCT(multiple_unit_cell)
1259 148 : conn_info => topology%conn_info
1260 :
1261 148 : nbond = 0
1262 148 : IF (ASSOCIATED(conn_info%bond_a)) THEN
1263 148 : nbond = SIZE(conn_info%bond_a)
1264 148 : CALL reallocate(conn_info%bond_a, 1, nbond*fac)
1265 148 : CALL reallocate(conn_info%bond_b, 1, nbond*fac)
1266 : END IF
1267 :
1268 148 : ntheta = 0
1269 148 : IF (ASSOCIATED(conn_info%theta_a)) THEN
1270 148 : ntheta = SIZE(conn_info%theta_a)
1271 148 : CALL reallocate(conn_info%theta_a, 1, ntheta*fac)
1272 148 : CALL reallocate(conn_info%theta_b, 1, ntheta*fac)
1273 148 : CALL reallocate(conn_info%theta_c, 1, ntheta*fac)
1274 : END IF
1275 :
1276 148 : nphi = 0
1277 148 : IF (ASSOCIATED(conn_info%phi_a)) THEN
1278 148 : nphi = SIZE(conn_info%phi_a)
1279 148 : CALL reallocate(conn_info%phi_a, 1, nphi*fac)
1280 148 : CALL reallocate(conn_info%phi_b, 1, nphi*fac)
1281 148 : CALL reallocate(conn_info%phi_c, 1, nphi*fac)
1282 148 : CALL reallocate(conn_info%phi_d, 1, nphi*fac)
1283 : END IF
1284 :
1285 148 : nimpr = 0
1286 148 : IF (ASSOCIATED(conn_info%impr_a)) THEN
1287 148 : nimpr = SIZE(conn_info%impr_a)
1288 148 : CALL reallocate(conn_info%impr_a, 1, nimpr*fac)
1289 148 : CALL reallocate(conn_info%impr_b, 1, nimpr*fac)
1290 148 : CALL reallocate(conn_info%impr_c, 1, nimpr*fac)
1291 148 : CALL reallocate(conn_info%impr_d, 1, nimpr*fac)
1292 : END IF
1293 :
1294 148 : nbond_c = 0
1295 148 : IF (ASSOCIATED(conn_info%c_bond_a)) THEN
1296 116 : nbond_c = SIZE(conn_info%c_bond_a)
1297 116 : CALL reallocate(conn_info%c_bond_a, 1, nbond_c*fac)
1298 116 : CALL reallocate(conn_info%c_bond_b, 1, nbond_c*fac)
1299 : END IF
1300 :
1301 148 : nub = 0
1302 148 : IF (ASSOCIATED(conn_info%ub_a)) THEN
1303 148 : nub = SIZE(conn_info%ub_a)
1304 148 : CALL reallocate(conn_info%ub_a, 1, nub*fac)
1305 148 : CALL reallocate(conn_info%ub_b, 1, nub*fac)
1306 148 : CALL reallocate(conn_info%ub_c, 1, nub*fac)
1307 : END IF
1308 :
1309 148 : nonfo = 0
1310 148 : IF (ASSOCIATED(conn_info%onfo_a)) THEN
1311 148 : nonfo = SIZE(conn_info%onfo_a)
1312 148 : CALL reallocate(conn_info%onfo_a, 1, nonfo*fac)
1313 148 : CALL reallocate(conn_info%onfo_b, 1, nonfo*fac)
1314 : END IF
1315 :
1316 148 : natoms_orig = topology%natoms/fac
1317 148 : ind = 0
1318 556 : DO k = 1, multiple_unit_cell(3)
1319 1870 : DO j = 1, multiple_unit_cell(2)
1320 6698 : DO i = 1, multiple_unit_cell(1)
1321 4976 : ind = ind + 1
1322 4976 : IF (ind == 1) CYCLE
1323 4828 : a = (ind - 1)*natoms_orig
1324 :
1325 : ! Bonds
1326 4828 : IF (nbond > 0) THEN
1327 3208 : m = (ind - 1)*nbond
1328 34024 : conn_info%bond_a(m + 1:m + nbond) = conn_info%bond_a(1:nbond) + a
1329 37232 : conn_info%bond_b(m + 1:m + nbond) = conn_info%bond_b(1:nbond) + a
1330 : END IF
1331 : ! Theta
1332 4828 : IF (ntheta > 0) THEN
1333 2168 : m = (ind - 1)*ntheta
1334 37428 : conn_info%theta_a(m + 1:m + ntheta) = conn_info%theta_a(1:ntheta) + a
1335 39596 : conn_info%theta_b(m + 1:m + ntheta) = conn_info%theta_b(1:ntheta) + a
1336 39596 : conn_info%theta_c(m + 1:m + ntheta) = conn_info%theta_c(1:ntheta) + a
1337 : END IF
1338 : ! Phi
1339 4828 : IF (nphi > 0) THEN
1340 14 : m = (ind - 1)*nphi
1341 35756 : conn_info%phi_a(m + 1:m + nphi) = conn_info%phi_a(1:nphi) + a
1342 35770 : conn_info%phi_b(m + 1:m + nphi) = conn_info%phi_b(1:nphi) + a
1343 35770 : conn_info%phi_c(m + 1:m + nphi) = conn_info%phi_c(1:nphi) + a
1344 35770 : conn_info%phi_d(m + 1:m + nphi) = conn_info%phi_d(1:nphi) + a
1345 : END IF
1346 : ! Impropers
1347 4828 : IF (nimpr > 0) THEN
1348 0 : m = (ind - 1)*nimpr
1349 0 : conn_info%impr_a(m + 1:m + nimpr) = conn_info%impr_a(1:nimpr) + a
1350 0 : conn_info%impr_b(m + 1:m + nimpr) = conn_info%impr_b(1:nimpr) + a
1351 0 : conn_info%impr_c(m + 1:m + nimpr) = conn_info%impr_c(1:nimpr) + a
1352 0 : conn_info%impr_d(m + 1:m + nimpr) = conn_info%impr_d(1:nimpr) + a
1353 : END IF
1354 : ! Para_res
1355 4828 : IF (nbond_c > 0) THEN
1356 0 : m = (ind - 1)*nbond_c
1357 0 : conn_info%c_bond_a(m + 1:m + nbond_c) = conn_info%c_bond_a(1:nbond_c) + a
1358 0 : conn_info%c_bond_b(m + 1:m + nbond_c) = conn_info%c_bond_b(1:nbond_c) + a
1359 : END IF
1360 : ! Urey-Bradley
1361 4828 : IF (nub > 0) THEN
1362 2168 : m = (ind - 1)*nub
1363 37428 : conn_info%ub_a(m + 1:m + nub) = conn_info%ub_a(1:nub) + a
1364 39596 : conn_info%ub_b(m + 1:m + nub) = conn_info%ub_b(1:nub) + a
1365 39596 : conn_info%ub_c(m + 1:m + nub) = conn_info%ub_c(1:nub) + a
1366 : END IF
1367 : ! ONFO
1368 6142 : IF (nonfo > 0) THEN
1369 14 : m = (ind - 1)*nonfo
1370 28364 : conn_info%onfo_a(m + 1:m + nonfo) = conn_info%onfo_a(1:nonfo) + a
1371 28378 : conn_info%onfo_b(m + 1:m + nonfo) = conn_info%onfo_b(1:nonfo) + a
1372 : END IF
1373 : END DO
1374 : END DO
1375 : END DO
1376 : END IF
1377 :
1378 10274 : END SUBROUTINE topology_conn_multiple
1379 :
1380 : END MODULE topology_connectivity_util
|
