(git:374b731)
Loading...
Searching...
No Matches
topology_connectivity_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Collection of subroutine needed for topology related things
10!> \par History
11!> jgh (23-05-2004) Last atom of molecule information added
12! **************************************************************************************************
13
14MODULE 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
52CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: c_var_a, c_var_b, c_var_c, c_var_d, c_var_type, &
79 first_list, last_list, map_atom_mol, map_atom_type, map_cvar_mol, map_cvars, map_var_mol, &
80 map_vars, molecule_list
81 INTEGER, DIMENSION(:, :), POINTER :: bnd_ctype, bnd_type
82 LOGICAL :: found, found_last
83 TYPE(atom_info_type), POINTER :: atom_info
84 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
85 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list
86 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
87 TYPE(connectivity_info_type), POINTER :: conn_info
88 TYPE(cp_logger_type), POINTER :: logger
89 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list
90 TYPE(local_states_type), DIMENSION(:), POINTER :: lmi
91 TYPE(molecule_kind_type), POINTER :: molecule_kind
92 TYPE(molecule_type), POINTER :: molecule
93 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list
94 TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list
95 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list
96
97 NULLIFY (logger)
98 logger => cp_get_default_logger()
99 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 extension=".subsysLog")
102 CALL timeset(routinen, handle)
103
104 atom_info => topology%atom_info
105 conn_info => topology%conn_info
106 ALLOCATE (map_atom_mol(topology%natoms))
107 ALLOCATE (map_atom_type(topology%natoms))
108 !-----------------------------------------------------------------------------
109 !-----------------------------------------------------------------------------
110 ! 1. Set the topology%[nmol_type,nmol,nmol_conn]
111 !-----------------------------------------------------------------------------
112 CALL timeset(routinen//"_1", handle2)
113 natom = topology%natoms
114 topology%nmol = 1
115 topology%nmol_type = 1
116 topology%nmol_conn = 0
117 map_atom_mol = -1
118 map_atom_type = -1
119 map_atom_mol(1) = 1
120 map_atom_type(1) = 1
121 DO i = 1, natom - 1
122 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 topology%nmol_type = topology%nmol_type + 1
126 END IF
127 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 (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 topology%nmol = topology%nmol + 1
132 END IF
133 map_atom_mol(i + 1) = topology%nmol
134 IF ((atom_info%map_mol_typ(i + 1) == atom_info%map_mol_typ(i)) .AND. &
135 (atom_info%map_mol_num(i + 1) == atom_info%map_mol_num(i)) .AND. &
136 (atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i))) THEN
137 topology%nmol_conn = topology%nmol_conn + 1
138 END IF
139 END DO
140 IF (iw > 0) WRITE (iw, *) "topology%nmol ::", topology%nmol
141 IF (iw > 0) WRITE (iw, *) "topology%nmol_type ::", topology%nmol_type
142 IF (iw > 0) WRITE (iw, *) "topology%nmol_conn ::", topology%nmol_conn
143 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 CALL timeset(routinen//"_1.1", handle2)
150 istart_mol = map_atom_mol(1)
151 istart_typ = map_atom_type(1)
152 DO i = 2, natom
153 IF ((map_atom_mol(i) /= istart_mol) .AND. (map_atom_type(i) == istart_typ)) THEN
154 map_atom_mol(i) = -map_atom_mol(i)
155 ELSE IF (map_atom_type(i) /= istart_typ) THEN
156 istart_mol = map_atom_mol(i)
157 istart_typ = map_atom_type(i)
158 END IF
159 END DO
160 CALL timestop(handle2)
161 !-----------------------------------------------------------------------------
162 !-----------------------------------------------------------------------------
163 ! 2. Allocate the molecule_kind_set
164 !-----------------------------------------------------------------------------
165 CALL timeset(routinen//"_2", handle2)
166 IF (topology%nmol_type <= 0) THEN
167 cpabort("No molecule kind defined")
168 ELSE
169 NULLIFY (molecule_kind_set)
170 i = topology%nmol_type
171 CALL allocate_molecule_kind_set(molecule_kind_set, i)
172 IF (iw > 0) WRITE (iw, *) " Allocated molecule_kind_set, Dimenstion of ", &
173 SIZE(molecule_kind_set)
174 END IF
175 CALL timestop(handle2)
176 !-----------------------------------------------------------------------------
177 !-----------------------------------------------------------------------------
178 ! 3. Allocate the molecule_set
179 !-----------------------------------------------------------------------------
180 CALL timeset(routinen//"_3", handle2)
181 IF (topology%nmol <= 0) THEN
182 cpabort("No molecule defined")
183 ELSE
184 NULLIFY (molecule_set)
185 i = topology%nmol
186 CALL allocate_molecule_set(molecule_set, i)
187 IF (iw > 0) WRITE (iw, *) " Allocated molecule_set, dimension of ", &
188 topology%nmol
189 END IF
190 CALL timestop(handle2)
191 !-----------------------------------------------------------------------------
192 !-----------------------------------------------------------------------------
193 ! 4. Set the molecule_kind_set%[kind_number,name,nsgf,nelectron]
194 !-----------------------------------------------------------------------------
195 CALL timeset(routinen//"_4", handle2)
196 natom = topology%natoms
197 ikind = -1
198 DO i = 1, natom
199 IF (ikind /= map_atom_type(i)) THEN
200 ikind = map_atom_type(i)
201 molecule_kind => molecule_kind_set(ikind)
202 nsgf = 0
203 nelectron = 0
204 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 nelectron=nelectron)
211 END IF
212 END DO
213 CALL timestop(handle2)
214 !-----------------------------------------------------------------------------
215 !-----------------------------------------------------------------------------
216 ! 5. Set the molecule_list for molecule_kind in molecule_kind_set
217 !-----------------------------------------------------------------------------
218 CALL timeset(routinen//"_5", handle2)
219 natom = topology%natoms
220 ikind = map_atom_type(1)
221 imol = abs(map_atom_mol(1))
222 counter = 0
223 DO i = 1, natom - 1
224 IF (ikind /= map_atom_type(i + 1)) THEN
225 found = .true.
226 found_last = .false.
227 imol = abs(map_atom_mol(i))
228 ELSEIF (ikind == topology%nmol_type) THEN
229 found = .true.
230 found_last = .true.
231 imol = abs(map_atom_mol(natom))
232 ELSE
233 found = .false.
234 found_last = .false.
235 END IF
236
237 IF (found) THEN
238 ALLOCATE (molecule_list(imol - counter))
239 DO j = 1, SIZE(molecule_list)
240 molecule_list(j) = j + counter
241 END DO
242 molecule_kind => molecule_kind_set(ikind)
243 CALL set_molecule_kind(molecule_kind=molecule_kind, &
244 molecule_list=molecule_list)
245 IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:)
246 IF (found_last) EXIT
247 counter = imol
248 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 IF (i == natom) THEN
253 imol = abs(map_atom_mol(natom))
254 ! Last atom is also a molecule by itself
255 ALLOCATE (molecule_list(imol - counter))
256 DO j = 1, SIZE(molecule_list)
257 molecule_list(j) = j + counter
258 END DO
259 molecule_kind => molecule_kind_set(ikind)
260 CALL set_molecule_kind(molecule_kind=molecule_kind, &
261 molecule_list=molecule_list)
262 IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:)
263 END IF
264 CALL timestop(handle2)
265 !-----------------------------------------------------------------------------
266 !-----------------------------------------------------------------------------
267 ! 6. Set the molecule_set(imol)%molecule_kind via set_molecule
268 !-----------------------------------------------------------------------------
269 CALL timeset(routinen//"_6", handle2)
270 DO ikind = 1, SIZE(molecule_kind_set)
271 molecule_kind => molecule_kind_set(ikind)
272 CALL get_molecule_kind(molecule_kind=molecule_kind, &
273 molecule_list=molecule_list)
274 DO i = 1, SIZE(molecule_list)
275 molecule => molecule_set(molecule_list(i))
276 CALL set_molecule(molecule, molecule_kind=molecule_kind)
277 END DO
278 END DO
279 CALL timestop(handle2)
280 !-----------------------------------------------------------------------------
281 !-----------------------------------------------------------------------------
282 ! 7. Set the molecule_set(imol)%[first_atom,last_atom] via set_molecule_set
283 !-----------------------------------------------------------------------------
284 ALLOCATE (first_list(SIZE(molecule_set)))
285 ALLOCATE (last_list(SIZE(molecule_set)))
286 CALL timeset(routinen//"_7", handle2)
287 first_list(:) = 0
288 last_list(:) = 0
289 ityp = atom_info%map_mol_typ(1)
290 inum = atom_info%map_mol_num(1)
291 ires = atom_info%map_mol_res(1)
292 imol = 1
293 first_list(1) = 1
294 DO j = 2, natom
295 IF ((atom_info%map_mol_typ(j) /= ityp) .OR. &
296 (atom_info%map_mol_num(j) /= inum) .OR. &
297 (atom_info%map_mol_res(j) /= ires)) THEN
298 ityp = atom_info%map_mol_typ(j)
299 inum = atom_info%map_mol_num(j)
300 ires = atom_info%map_mol_res(j)
301 imol = imol + 1
302 first_list(imol) = j
303 END IF
304 END DO
305 cpassert(imol == topology%nmol)
306 DO ikind = 1, topology%nmol - 1
307 last_list(ikind) = first_list(ikind + 1) - 1
308 END DO
309 last_list(topology%nmol) = topology%natoms
310 CALL set_molecule_set(molecule_set, first_list, last_list)
311 CALL timestop(handle2)
312 !-----------------------------------------------------------------------------
313 !-----------------------------------------------------------------------------
314 ! 8. Set and NULLIFY the molecule_set(imol)%lmi via set_molecule_set
315 !-----------------------------------------------------------------------------
316 CALL timeset(routinen//"_8", handle2)
317 DO imol = 1, SIZE(molecule_set)
318 molecule => molecule_set(imol)
319 NULLIFY (lmi)
320 CALL set_molecule(molecule, lmi=lmi)
321 END DO
322 CALL timestop(handle2)
323 !-----------------------------------------------------------------------------
324 !-----------------------------------------------------------------------------
325 ! 9. Set the atom_list for molecule_kind in molecule_kind_set (PART 1)
326 !-----------------------------------------------------------------------------
327 CALL timeset(routinen//"_9", handle2)
328 counter = 0
329 DO imol = 1, SIZE(molecule_set)
330 molecule => molecule_set(imol)
331 molecule_kind => molecule_set(imol)%molecule_kind
332 CALL get_molecule_kind(molecule_kind=molecule_kind, &
333 kind_number=i)
334 IF (counter /= i) THEN
335 counter = i
336 CALL get_molecule(molecule=molecule, &
337 first_atom=first, last_atom=last)
338 natom = 0
339 IF (first /= 0 .AND. last /= 0) natom = last - first + 1
340 ALLOCATE (atom_list(natom))
341 DO i = 1, natom
342 !Atomic kind information will be filled in (PART 2)
343 NULLIFY (atom_list(i)%atomic_kind)
344 atom_list(i)%id_name = atom_info%id_atmname(i + first - 1)
345 IF (iw > 0) WRITE (iw, '(5X,A,3I5,1X,A5)') "atom_list ", &
346 imol, counter, i, trim(id2str(atom_list(i)%id_name))
347 END DO
348 CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
349 END IF
350 END DO
351 CALL timestop(handle2)
352 !-----------------------------------------------------------------------------
353 !-----------------------------------------------------------------------------
354 ! 10. Set the molecule_kind%[nbond,bond_list] via set_molecule_kind
355 !-----------------------------------------------------------------------------
356 CALL timeset(routinen//"_10", handle2)
357 ! First map bonds on molecules
358 nvar1 = 0
359 nvar2 = 0
360 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
361 IF (ASSOCIATED(conn_info%bond_a)) nvar1 = SIZE(conn_info%bond_a)
362 IF (ASSOCIATED(conn_info%c_bond_a)) nvar2 = SIZE(conn_info%c_bond_a)
363 nval_tot1 = nvar1
364 nval_tot2 = 0
365 ALLOCATE (map_var_mol(nvar1))
366 ALLOCATE (map_cvar_mol(nvar2))
367 map_var_mol = -1
368 map_cvar_mol = -1
369 DO i = 1, nvar1
370 j1 = map_atom_mol(conn_info%bond_a(i))
371 j2 = map_atom_mol(conn_info%bond_b(i))
372 IF (j1 == j2) THEN
373 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%bond_a(i))
374 END IF
375 END DO
376 DO i = 1, nvar2
377 min_index = min(conn_info%c_bond_a(i), conn_info%c_bond_b(i))
378 j1 = map_atom_mol(min_index)
379 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
380 END DO
381 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
382 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
383 DO i = 1, topology%nmol_type
384 intra_bonds = 0
385 inter_bonds = 0
386 IF (all(bnd_type(:, i) > 0)) THEN
387 intra_bonds = bnd_type(2, i) - bnd_type(1, i) + 1
388 END IF
389 IF (all(bnd_ctype(:, i) > 0)) THEN
390 inter_bonds = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
391 END IF
392 ibond = intra_bonds + inter_bonds
393 IF (iw > 0) THEN
394 WRITE (iw, *) " Total number bonds for molecule type ", i, " :", ibond
395 WRITE (iw, *) " intra (bonds inside molecules) :: ", intra_bonds
396 WRITE (iw, *) " inter (bonds between molecules) :: ", inter_bonds
397 END IF
398 molecule_kind => molecule_kind_set(i)
399 nval_tot2 = nval_tot2 + ibond*SIZE(molecule_kind%molecule_list)
400
401 ALLOCATE (bond_list(ibond))
402 ibond = 0
403 DO j = bnd_type(1, i), bnd_type(2, i)
404 IF (j == 0) cycle
405 ibond = ibond + 1
406 jind = map_vars(j)
407 first = first_list(map_atom_mol(conn_info%bond_a(jind)))
408 bond_list(ibond)%a = conn_info%bond_a(jind) - first + 1
409 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 bond_list(ibond)%id_type = do_ff_charmm
412 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
413 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 NULLIFY (bond_list(ibond)%bond_kind)
417 IF (iw > 0) THEN
418 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", &
419 i, " intra bond", &
420 conn_info%bond_a(jind), &
421 conn_info%bond_b(jind), &
422 "offset number at", &
423 conn_info%bond_a(jind) - first + 1, &
424 conn_info%bond_b(jind) - first + 1
425 END IF
426 END DO
427 DO j = bnd_ctype(1, i), bnd_ctype(2, i)
428 IF (j == 0) cycle
429 ibond = ibond + 1
430 jind = map_cvars(j)
431 min_index = min(conn_info%c_bond_a(jind), conn_info%c_bond_b(jind))
432 first = first_list(map_atom_mol(min_index))
433 bond_list(ibond)%a = conn_info%c_bond_a(jind) - first + 1
434 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 bond_list(ibond)%id_type = do_ff_charmm
437 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
438 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 NULLIFY (bond_list(ibond)%bond_kind)
442 IF (iw > 0) THEN
443 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", &
444 i, " inter bond", &
445 conn_info%c_bond_a(jind), &
446 conn_info%c_bond_b(jind), &
447 "offset number at", &
448 conn_info%c_bond_a(jind) - first + 1, &
449 conn_info%c_bond_b(jind) - first + 1
450 END IF
451 END DO
452 CALL set_molecule_kind(molecule_kind=molecule_kind, &
453 nbond=SIZE(bond_list), bond_list=bond_list)
454 END DO
455 IF ((nval_tot1 /= nval_tot2) .AND. (output_unit > 0)) THEN
456 WRITE (output_unit, '(/)')
457 WRITE (output_unit, '(T5,A)') "ERROR| Mismatching found between the total number of atoms"
458 WRITE (output_unit, '(T5,A)') "ERROR| and the number of atoms computed multiplying the Nr."
459 WRITE (output_unit, '(T5,A)') "ERROR| of molecules by the number of atoms building that"
460 WRITE (output_unit, '(T5,A)') "ERROR| kind of molecule."
461 WRITE (output_unit, '(T5,A)') "ERROR| This happens when the connectivity is wrongly built"
462 WRITE (output_unit, '(T5,A)') "ERROR| One example could be two same kind of molecules have"
463 WRITE (output_unit, '(T5,A)') "ERROR| a different number of atoms. Check the connectivity!"
464 END IF
465 cpassert(nval_tot1 == nval_tot2)
466 DEALLOCATE (map_var_mol)
467 DEALLOCATE (map_cvar_mol)
468 DEALLOCATE (map_vars)
469 DEALLOCATE (map_cvars)
470 DEALLOCATE (bnd_type)
471 DEALLOCATE (bnd_ctype)
472 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 CALL timeset(routinen//"_11_pre", handle2)
479 idim = 0
480 ALLOCATE (c_var_a(idim))
481 ALLOCATE (c_var_b(idim))
482 ALLOCATE (c_var_c(idim))
483 found = ASSOCIATED(conn_info%theta_type)
484 IF (found) THEN
485 ALLOCATE (c_var_type(idim))
486 END IF
487 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%theta_a)) THEN
488 DO j = 1, SIZE(conn_info%theta_a)
489 j1 = map_atom_mol(conn_info%theta_a(j))
490 j2 = map_atom_mol(conn_info%theta_b(j))
491 j3 = map_atom_mol(conn_info%theta_c(j))
492 IF (j1 /= j2 .OR. j2 /= j3) THEN
493 idim = idim + 1
494 END IF
495 END DO
496 CALL reallocate(c_var_a, 1, idim)
497 CALL reallocate(c_var_b, 1, idim)
498 CALL reallocate(c_var_c, 1, idim)
499 IF (found) THEN
500 CALL reallocate(c_var_type, 1, idim)
501 END IF
502 idim = 0
503 DO j = 1, SIZE(conn_info%theta_a)
504 j1 = map_atom_mol(conn_info%theta_a(j))
505 j2 = map_atom_mol(conn_info%theta_b(j))
506 j3 = map_atom_mol(conn_info%theta_c(j))
507 IF (j1 /= j2 .OR. j2 /= j3) THEN
508 idim = idim + 1
509 c_var_a(idim) = conn_info%theta_a(j)
510 c_var_b(idim) = conn_info%theta_b(j)
511 c_var_c(idim) = conn_info%theta_c(j)
512 IF (found) THEN
513 c_var_type(idim) = conn_info%theta_type(j)
514 END IF
515 END IF
516 END DO
517 END IF
518 CALL timestop(handle2)
519 CALL timeset(routinen//"_11", handle2)
520 ! map bends on molecules
521 nvar1 = 0
522 nvar2 = 0
523 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
524 IF (ASSOCIATED(conn_info%theta_a)) nvar1 = SIZE(conn_info%theta_a)
525 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
526 nval_tot1 = nvar1
527 nval_tot2 = 0
528 ALLOCATE (map_var_mol(nvar1))
529 ALLOCATE (map_cvar_mol(nvar2))
530 map_var_mol = -1
531 map_cvar_mol = -1
532 DO i = 1, nvar1
533 j1 = map_atom_mol(conn_info%theta_a(i))
534 j2 = map_atom_mol(conn_info%theta_b(i))
535 j3 = map_atom_mol(conn_info%theta_c(i))
536 IF (j1 == j2 .AND. j2 == j3) THEN
537 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%theta_a(i))
538 END IF
539 END DO
540 DO i = 1, nvar2
541 min_index = min(c_var_a(i), c_var_b(i), c_var_c(i))
542 j1 = map_atom_mol(min_index)
543 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
544 END DO
545 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
546 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
547 DO i = 1, topology%nmol_type
548 intra_bends = 0
549 inter_bends = 0
550 IF (all(bnd_type(:, i) > 0)) THEN
551 intra_bends = bnd_type(2, i) - bnd_type(1, i) + 1
552 END IF
553 IF (all(bnd_ctype(:, i) > 0)) THEN
554 inter_bends = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
555 END IF
556 ibend = intra_bends + inter_bends
557 IF (iw > 0) THEN
558 WRITE (iw, *) " Total number of angles for molecule type ", i, " :", ibend
559 WRITE (iw, *) " intra (angles inside molecules) :: ", intra_bends
560 WRITE (iw, *) " inter (angles between molecules) :: ", inter_bends
561 END IF
562 molecule_kind => molecule_kind_set(i)
563 nval_tot2 = nval_tot2 + ibend*SIZE(molecule_kind%molecule_list)
564 ALLOCATE (bend_list(ibend))
565 ibend = 0
566 DO j = bnd_type(1, i), bnd_type(2, i)
567 IF (j == 0) cycle
568 ibend = ibend + 1
569 jind = map_vars(j)
570 first = first_list(map_atom_mol(conn_info%theta_a(jind)))
571 bend_list(ibend)%a = conn_info%theta_a(jind) - first + 1
572 bend_list(ibend)%b = conn_info%theta_b(jind) - first + 1
573 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 bend_list(ibend)%id_type = do_ff_charmm
576 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
577 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 NULLIFY (bend_list(ibend)%bend_kind)
581 IF (iw > 0) THEN
582 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
583 "molecule_kind", ikind, "intra bend", &
584 conn_info%theta_a(jind), &
585 conn_info%theta_b(jind), &
586 conn_info%theta_c(jind), &
587 "offset number at", &
588 conn_info%theta_a(jind) - first + 1, &
589 conn_info%theta_b(jind) - first + 1, &
590 conn_info%theta_c(jind) - first + 1
591 END IF
592 END DO
593 DO j = bnd_ctype(1, i), bnd_ctype(2, i)
594 IF (j == 0) cycle
595 ibend = ibend + 1
596 jind = map_cvars(j)
597 min_index = min(c_var_a(jind), c_var_b(jind), c_var_c(jind))
598 first = first_list(map_atom_mol(min_index))
599 bend_list(ibend)%a = c_var_a(jind) - first + 1
600 bend_list(ibend)%b = c_var_b(jind) - first + 1
601 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 bend_list(ibend)%id_type = do_ff_charmm
604 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
605 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 NULLIFY (bend_list(ibend)%bend_kind)
609 IF (iw > 0) THEN
610 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
611 "molecule_kind", ikind, "inter bend", &
612 c_var_a(jind), &
613 c_var_b(jind), &
614 c_var_c(jind), &
615 "offset number at", &
616 c_var_a(jind) - first + 1, &
617 c_var_b(jind) - first + 1, &
618 c_var_c(jind) - first + 1
619 END IF
620 END DO
621 CALL set_molecule_kind(molecule_kind=molecule_kind, &
622 nbend=SIZE(bend_list), bend_list=bend_list)
623 END DO
624 cpassert(nval_tot1 == nval_tot2)
625 DEALLOCATE (map_var_mol)
626 DEALLOCATE (map_cvar_mol)
627 DEALLOCATE (map_vars)
628 DEALLOCATE (map_cvars)
629 DEALLOCATE (bnd_type)
630 DEALLOCATE (bnd_ctype)
631 DEALLOCATE (c_var_a)
632 DEALLOCATE (c_var_b)
633 DEALLOCATE (c_var_c)
634 IF (found) THEN
635 DEALLOCATE (c_var_type)
636 END IF
637 CALL timestop(handle2)
638 !-----------------------------------------------------------------------------
639 !-----------------------------------------------------------------------------
640 ! 12. Set the molecule_kind%[nub,ub_list] via set_molecule_kind
641 !-----------------------------------------------------------------------------
642 CALL timeset(routinen//"_12_pre", handle2)
643 idim = 0
644 ALLOCATE (c_var_a(idim))
645 ALLOCATE (c_var_b(idim))
646 ALLOCATE (c_var_c(idim))
647 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%ub_a)) THEN
648 DO j = 1, SIZE(conn_info%ub_a)
649 j1 = map_atom_mol(conn_info%ub_a(j))
650 j2 = map_atom_mol(conn_info%ub_b(j))
651 j3 = map_atom_mol(conn_info%ub_c(j))
652 IF (j1 /= j2 .OR. j2 /= j3) THEN
653 idim = idim + 1
654 END IF
655 END DO
656 CALL reallocate(c_var_a, 1, idim)
657 CALL reallocate(c_var_b, 1, idim)
658 CALL reallocate(c_var_c, 1, idim)
659 idim = 0
660 DO j = 1, SIZE(conn_info%ub_a)
661 j1 = map_atom_mol(conn_info%ub_a(j))
662 j2 = map_atom_mol(conn_info%ub_b(j))
663 j3 = map_atom_mol(conn_info%ub_c(j))
664 IF (j1 /= j2 .OR. j2 /= j3) THEN
665 idim = idim + 1
666 c_var_a(idim) = conn_info%ub_a(j)
667 c_var_b(idim) = conn_info%ub_b(j)
668 c_var_c(idim) = conn_info%ub_c(j)
669 END IF
670 END DO
671 END IF
672 CALL timestop(handle2)
673 CALL timeset(routinen//"_12", handle2)
674 ! map UBs on molecules
675 nvar1 = 0
676 nvar2 = 0
677 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
678 IF (ASSOCIATED(conn_info%ub_a)) nvar1 = SIZE(conn_info%ub_a)
679 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
680 nval_tot1 = nvar1
681 nval_tot2 = 0
682 ALLOCATE (map_var_mol(nvar1))
683 ALLOCATE (map_cvar_mol(nvar2))
684 map_var_mol = -1
685 map_cvar_mol = -1
686 DO i = 1, nvar1
687 j1 = map_atom_mol(conn_info%ub_a(i))
688 j2 = map_atom_mol(conn_info%ub_b(i))
689 j3 = map_atom_mol(conn_info%ub_c(i))
690 IF (j1 == j2 .AND. j2 == j3) THEN
691 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%ub_a(i))
692 END IF
693 END DO
694 DO i = 1, nvar2
695 min_index = min(c_var_a(i), c_var_b(i), c_var_c(i))
696 j1 = map_atom_mol(min_index)
697 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
698 END DO
699 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
700 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
701 DO i = 1, topology%nmol_type
702 intra_ubs = 0
703 inter_ubs = 0
704 IF (all(bnd_type(:, i) > 0)) THEN
705 intra_ubs = bnd_type(2, i) - bnd_type(1, i) + 1
706 END IF
707 IF (all(bnd_ctype(:, i) > 0)) THEN
708 inter_ubs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
709 END IF
710 iub = intra_ubs + inter_ubs
711 IF (iw > 0) THEN
712 WRITE (iw, *) " Total number of Urey-Bradley for molecule type ", i, " :", iub
713 WRITE (iw, *) " intra (UB inside molecules) :: ", intra_ubs
714 WRITE (iw, *) " inter (UB between molecules) :: ", inter_ubs
715 END IF
716 molecule_kind => molecule_kind_set(i)
717 nval_tot2 = nval_tot2 + iub*SIZE(molecule_kind%molecule_list)
718 ALLOCATE (ub_list(iub))
719 iub = 0
720 DO j = bnd_type(1, i), bnd_type(2, i)
721 IF (j == 0) cycle
722 iub = iub + 1
723 jind = map_vars(j)
724 first = first_list(map_atom_mol(conn_info%ub_a(jind)))
725 ub_list(iub)%a = conn_info%ub_a(jind) - first + 1
726 ub_list(iub)%b = conn_info%ub_b(jind) - first + 1
727 ub_list(iub)%c = conn_info%ub_c(jind) - first + 1
728 ub_list(iub)%id_type = do_ff_charmm
729 !point this to the right ub_kind_type if using force field
730 NULLIFY (ub_list(iub)%ub_kind)
731 IF (iw > 0) THEN
732 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
733 "molecule_kind", i, "intra UB", &
734 conn_info%ub_a(jind), &
735 conn_info%ub_b(jind), &
736 conn_info%ub_c(jind), &
737 "offset number at", &
738 conn_info%ub_a(jind) - first + 1, &
739 conn_info%ub_b(jind) - first + 1, &
740 conn_info%ub_c(jind) - first + 1
741 END IF
742 END DO
743 DO j = bnd_ctype(1, i), bnd_ctype(2, i)
744 IF (j == 0) cycle
745 iub = iub + 1
746 jind = map_cvars(j)
747 min_index = min(c_var_a(jind), c_var_b(jind), c_var_c(jind))
748 first = first_list(map_atom_mol(min_index))
749 ub_list(iub)%a = c_var_a(jind) - first + 1
750 ub_list(iub)%b = c_var_b(jind) - first + 1
751 ub_list(iub)%c = c_var_c(jind) - first + 1
752 ub_list(iub)%id_type = do_ff_charmm
753 !point this to the right ub_kind_type if using force field
754 NULLIFY (ub_list(iub)%ub_kind)
755 IF (iw > 0) THEN
756 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') &
757 "molecule_kind", i, "inter UB", &
758 c_var_a(jind), &
759 c_var_b(jind), &
760 c_var_c(jind), &
761 "offset number at", &
762 c_var_a(jind) - first + 1, &
763 c_var_b(jind) - first + 1, &
764 c_var_c(jind) - first + 1
765 END IF
766 END DO
767 CALL set_molecule_kind(molecule_kind=molecule_kind, &
768 nub=SIZE(ub_list), ub_list=ub_list)
769 END DO
770 cpassert(nval_tot1 == nval_tot2)
771 DEALLOCATE (map_var_mol)
772 DEALLOCATE (map_cvar_mol)
773 DEALLOCATE (map_vars)
774 DEALLOCATE (map_cvars)
775 DEALLOCATE (bnd_type)
776 DEALLOCATE (bnd_ctype)
777 DEALLOCATE (c_var_a)
778 DEALLOCATE (c_var_b)
779 DEALLOCATE (c_var_c)
780 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 CALL timeset(routinen//"_13_pre", handle2)
787 idim = 0
788 ALLOCATE (c_var_a(idim))
789 ALLOCATE (c_var_b(idim))
790 ALLOCATE (c_var_c(idim))
791 ALLOCATE (c_var_d(idim))
792 found = ASSOCIATED(conn_info%phi_type)
793 IF (found) THEN
794 ALLOCATE (c_var_type(idim))
795 END IF
796 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%phi_a)) THEN
797 DO j = 1, SIZE(conn_info%phi_a)
798 j1 = map_atom_mol(conn_info%phi_a(j))
799 j2 = map_atom_mol(conn_info%phi_b(j))
800 j3 = map_atom_mol(conn_info%phi_c(j))
801 j4 = map_atom_mol(conn_info%phi_d(j))
802 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
803 idim = idim + 1
804 END IF
805 END DO
806 CALL reallocate(c_var_a, 1, idim)
807 CALL reallocate(c_var_b, 1, idim)
808 CALL reallocate(c_var_c, 1, idim)
809 CALL reallocate(c_var_d, 1, idim)
810 IF (found) THEN
811 CALL reallocate(c_var_type, 1, idim)
812 END IF
813 idim = 0
814 DO j = 1, SIZE(conn_info%phi_a)
815 j1 = map_atom_mol(conn_info%phi_a(j))
816 j2 = map_atom_mol(conn_info%phi_b(j))
817 j3 = map_atom_mol(conn_info%phi_c(j))
818 j4 = map_atom_mol(conn_info%phi_d(j))
819 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
820 idim = idim + 1
821 c_var_a(idim) = conn_info%phi_a(j)
822 c_var_b(idim) = conn_info%phi_b(j)
823 c_var_c(idim) = conn_info%phi_c(j)
824 c_var_d(idim) = conn_info%phi_d(j)
825 IF (found) THEN
826 c_var_type(idim) = conn_info%phi_type(j)
827 END IF
828 END IF
829 END DO
830 END IF
831 CALL timestop(handle2)
832 CALL timeset(routinen//"_13", handle2)
833 ! map torsions on molecules
834 nvar1 = 0
835 nvar2 = 0
836 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
837 IF (ASSOCIATED(conn_info%phi_a)) nvar1 = SIZE(conn_info%phi_a)
838 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
839 nval_tot1 = nvar1
840 nval_tot2 = 0
841 ALLOCATE (map_var_mol(nvar1))
842 ALLOCATE (map_cvar_mol(nvar2))
843 map_var_mol = -1
844 map_cvar_mol = -1
845 DO i = 1, nvar1
846 j1 = map_atom_mol(conn_info%phi_a(i))
847 j2 = map_atom_mol(conn_info%phi_b(i))
848 j3 = map_atom_mol(conn_info%phi_c(i))
849 j4 = map_atom_mol(conn_info%phi_d(i))
850 IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN
851 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%phi_a(i))
852 END IF
853 END DO
854 DO i = 1, nvar2
855 min_index = min(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i))
856 j1 = map_atom_mol(min_index)
857 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
858 END DO
859 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
860 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
861 DO i = 1, topology%nmol_type
862 intra_torsions = 0
863 inter_torsions = 0
864 IF (all(bnd_type(:, i) > 0)) THEN
865 intra_torsions = bnd_type(2, i) - bnd_type(1, i) + 1
866 END IF
867 IF (all(bnd_ctype(:, i) > 0)) THEN
868 inter_torsions = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
869 END IF
870 itorsion = intra_torsions + inter_torsions
871 IF (iw > 0) THEN
872 WRITE (iw, *) " Total number of torsions for molecule type ", i, " :", itorsion
873 WRITE (iw, *) " intra (torsions inside molecules) :: ", intra_torsions
874 WRITE (iw, *) " inter (torsions between molecules) :: ", inter_torsions
875 END IF
876 molecule_kind => molecule_kind_set(i)
877 nval_tot2 = nval_tot2 + itorsion*SIZE(molecule_kind%molecule_list)
878 ALLOCATE (torsion_list(itorsion))
879 itorsion = 0
880 DO j = bnd_type(1, i), bnd_type(2, i)
881 IF (j == 0) cycle
882 itorsion = itorsion + 1
883 jind = map_vars(j)
884 first = first_list(map_atom_mol(conn_info%phi_a(jind)))
885 torsion_list(itorsion)%a = conn_info%phi_a(jind) - first + 1
886 torsion_list(itorsion)%b = conn_info%phi_b(jind) - first + 1
887 torsion_list(itorsion)%c = conn_info%phi_c(jind) - first + 1
888 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 torsion_list(itorsion)%id_type = do_ff_charmm
891 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
892 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 NULLIFY (torsion_list(itorsion)%torsion_kind)
896 IF (iw > 0) THEN
897 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
898 "molecule_kind", i, "intra TOR", &
899 conn_info%phi_a(jind), &
900 conn_info%phi_b(jind), &
901 conn_info%phi_c(jind), &
902 conn_info%phi_d(jind), &
903 "offset number at", &
904 conn_info%phi_a(jind) - first + 1, &
905 conn_info%phi_b(jind) - first + 1, &
906 conn_info%phi_c(jind) - first + 1, &
907 conn_info%phi_d(jind) - first + 1
908 END IF
909 END DO
910 DO j = bnd_ctype(1, i), bnd_ctype(2, i)
911 IF (j == 0) cycle
912 itorsion = itorsion + 1
913 jind = map_cvars(j)
914 min_index = min(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind))
915 first = first_list(map_atom_mol(min_index))
916 torsion_list(itorsion)%a = c_var_a(jind) - first + 1
917 torsion_list(itorsion)%b = c_var_b(jind) - first + 1
918 torsion_list(itorsion)%c = c_var_c(jind) - first + 1
919 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 torsion_list(itorsion)%id_type = do_ff_charmm
922 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
923 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 NULLIFY (torsion_list(itorsion)%torsion_kind)
927 IF (iw > 0) THEN
928 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
929 "molecule_kind", i, "inter TOR", &
930 c_var_a(jind), &
931 c_var_b(jind), &
932 c_var_c(jind), &
933 c_var_d(jind), &
934 "offset number at", &
935 c_var_a(jind) - first + 1, &
936 c_var_b(jind) - first + 1, &
937 c_var_c(jind) - first + 1, &
938 c_var_d(jind) - first + 1
939 END IF
940 END DO
941 CALL set_molecule_kind(molecule_kind=molecule_kind, &
942 ntorsion=SIZE(torsion_list), torsion_list=torsion_list)
943 END DO
944 cpassert(nval_tot1 == nval_tot2)
945 DEALLOCATE (map_var_mol)
946 DEALLOCATE (map_cvar_mol)
947 DEALLOCATE (map_vars)
948 DEALLOCATE (map_cvars)
949 DEALLOCATE (bnd_type)
950 DEALLOCATE (bnd_ctype)
951 DEALLOCATE (c_var_a)
952 DEALLOCATE (c_var_b)
953 DEALLOCATE (c_var_c)
954 DEALLOCATE (c_var_d)
955 IF (found) THEN
956 DEALLOCATE (c_var_type)
957 END IF
958 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 CALL timeset(routinen//"_14_pre", handle2)
966 idim = 0
967 ALLOCATE (c_var_a(idim))
968 ALLOCATE (c_var_b(idim))
969 ALLOCATE (c_var_c(idim))
970 ALLOCATE (c_var_d(idim))
971 found = ASSOCIATED(conn_info%impr_type)
972 IF (found) THEN
973 ALLOCATE (c_var_type(idim))
974 END IF
975 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%impr_a)) THEN
976 DO j = 1, SIZE(conn_info%impr_a)
977 j1 = map_atom_mol(conn_info%impr_a(j))
978 j2 = map_atom_mol(conn_info%impr_b(j))
979 j3 = map_atom_mol(conn_info%impr_c(j))
980 j4 = map_atom_mol(conn_info%impr_d(j))
981 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
982 idim = idim + 1
983 END IF
984 END DO
985 CALL reallocate(c_var_a, 1, idim)
986 CALL reallocate(c_var_b, 1, idim)
987 CALL reallocate(c_var_c, 1, idim)
988 CALL reallocate(c_var_d, 1, idim)
989 IF (found) THEN
990 CALL reallocate(c_var_type, 1, idim)
991 END IF
992 idim = 0
993 DO j = 1, SIZE(conn_info%impr_a)
994 j1 = map_atom_mol(conn_info%impr_a(j))
995 j2 = map_atom_mol(conn_info%impr_b(j))
996 j3 = map_atom_mol(conn_info%impr_c(j))
997 j4 = map_atom_mol(conn_info%impr_d(j))
998 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN
999 idim = idim + 1
1000 c_var_a(idim) = conn_info%impr_a(j)
1001 c_var_b(idim) = conn_info%impr_b(j)
1002 c_var_c(idim) = conn_info%impr_c(j)
1003 c_var_d(idim) = conn_info%impr_d(j)
1004 IF (found) THEN
1005 c_var_type(idim) = conn_info%impr_type(j)
1006 END IF
1007 END IF
1008 END DO
1009 END IF
1010 CALL timestop(handle2)
1011 CALL timeset(routinen//"_14", handle2)
1012 ! map imprs on molecules
1013 nvar1 = 0
1014 nvar2 = 0
1015 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype)
1016 IF (ASSOCIATED(conn_info%impr_a)) nvar1 = SIZE(conn_info%impr_a)
1017 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a)
1018 nval_tot1 = nvar1
1019 nval_tot2 = 0
1020 ALLOCATE (map_var_mol(nvar1))
1021 ALLOCATE (map_cvar_mol(nvar2))
1022 map_var_mol = -1
1023 map_cvar_mol = -1
1024 DO i = 1, nvar1
1025 j1 = map_atom_mol(conn_info%impr_a(i))
1026 j2 = map_atom_mol(conn_info%impr_b(i))
1027 j3 = map_atom_mol(conn_info%impr_c(i))
1028 j4 = map_atom_mol(conn_info%impr_d(i))
1029 IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN
1030 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%impr_a(i))
1031 END IF
1032 END DO
1033 DO i = 1, nvar2
1034 min_index = min(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i))
1035 j1 = map_atom_mol(min_index)
1036 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index)
1037 END DO
1038 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1)
1039 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2)
1040 DO i = 1, topology%nmol_type
1041 intra_imprs = 0
1042 inter_imprs = 0
1043 IF (all(bnd_type(:, i) > 0)) THEN
1044 intra_imprs = bnd_type(2, i) - bnd_type(1, i) + 1
1045 END IF
1046 IF (all(bnd_ctype(:, i) > 0)) THEN
1047 inter_imprs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1
1048 END IF
1049 iimpr = intra_imprs + inter_imprs
1050 IF (iw > 0) THEN
1051 WRITE (iw, *) " Total number of imprs for molecule type ", i, " :", iimpr
1052 WRITE (iw, *) " intra (imprs inside molecules) :: ", intra_imprs
1053 WRITE (iw, *) " inter (imprs between molecules) :: ", inter_imprs
1054 WRITE (iw, *) " Total number of opbends for molecule type ", i, " :", iimpr
1055 WRITE (iw, *) " intra (opbends inside molecules) :: ", intra_imprs
1056 WRITE (iw, *) " inter (opbends between molecules) :: ", inter_imprs
1057 END IF
1058 molecule_kind => molecule_kind_set(i)
1059 nval_tot2 = nval_tot2 + iimpr*SIZE(molecule_kind%molecule_list)
1060 ALLOCATE (impr_list(iimpr), stat=stat)
1061 ALLOCATE (opbend_list(iimpr), stat=stat)
1062 cpassert(stat == 0)
1063 iimpr = 0
1064 DO j = bnd_type(1, i), bnd_type(2, i)
1065 IF (j == 0) cycle
1066 iimpr = iimpr + 1
1067 jind = map_vars(j)
1068 first = first_list(map_atom_mol(conn_info%impr_a(jind)))
1069 impr_list(iimpr)%a = conn_info%impr_a(jind) - first + 1
1070 impr_list(iimpr)%b = conn_info%impr_b(jind) - first + 1
1071 impr_list(iimpr)%c = conn_info%impr_c(jind) - first + 1
1072 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 opbend_list(iimpr)%a = conn_info%impr_b(jind) - first + 1
1078 opbend_list(iimpr)%b = conn_info%impr_d(jind) - first + 1
1079 opbend_list(iimpr)%c = conn_info%impr_c(jind) - first + 1
1080 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 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 opbend_list(iimpr)%id_type = do_ff_harmonic
1085 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
1086 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 NULLIFY (impr_list(iimpr)%impr_kind)
1090 NULLIFY (opbend_list(iimpr)%opbend_kind)
1091 IF (iw > 0) THEN
1092 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1093 "molecule_kind", i, "intra IMPR", &
1094 conn_info%impr_a(jind), &
1095 conn_info%impr_b(jind), &
1096 conn_info%impr_c(jind), &
1097 conn_info%impr_d(jind), &
1098 "offset number at", &
1099 conn_info%impr_a(jind) - first + 1, &
1100 conn_info%impr_b(jind) - first + 1, &
1101 conn_info%impr_c(jind) - first + 1, &
1102 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 "molecule_kind", i, "intra OPBEND", &
1105 conn_info%impr_b(jind), &
1106 conn_info%impr_d(jind), &
1107 conn_info%impr_c(jind), &
1108 conn_info%impr_a(jind), &
1109 "offset number at", &
1110 conn_info%impr_b(jind) - first + 1, &
1111 conn_info%impr_d(jind) - first + 1, &
1112 conn_info%impr_c(jind) - first + 1, &
1113 conn_info%impr_a(jind) - first + 1
1114 END IF
1115 END DO
1116 DO j = bnd_ctype(1, i), bnd_ctype(2, i)
1117 IF (j == 0) cycle
1118 iimpr = iimpr + 1
1119 jind = map_cvars(j)
1120 min_index = min(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind))
1121 first = first_list(map_atom_mol(min_index))
1122 impr_list(iimpr)%a = c_var_a(jind) - first + 1
1123 impr_list(iimpr)%b = c_var_b(jind) - first + 1
1124 impr_list(iimpr)%c = c_var_c(jind) - first + 1
1125 impr_list(iimpr)%d = c_var_d(jind) - first + 1
1126 opbend_list(iimpr)%a = c_var_b(jind) - first + 1
1127 opbend_list(iimpr)%b = c_var_d(jind) - first + 1
1128 opbend_list(iimpr)%c = c_var_c(jind) - first + 1
1129 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 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 opbend_list(iimpr)%id_type = do_ff_harmonic
1134 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN
1135 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 NULLIFY (impr_list(iimpr)%impr_kind)
1139 NULLIFY (opbend_list(iimpr)%opbend_kind)
1140 IF (iw > 0) THEN
1141 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') &
1142 "molecule_kind", i, "inter IMPR", &
1143 c_var_a(jind), &
1144 c_var_b(jind), &
1145 c_var_c(jind), &
1146 c_var_d(jind), &
1147 "offset number at", &
1148 c_var_a(jind) - first + 1, &
1149 c_var_b(jind) - first + 1, &
1150 c_var_c(jind) - first + 1, &
1151 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 "molecule_kind", i, "inter OPBEND", &
1154 c_var_b(jind), &
1155 c_var_d(jind), &
1156 c_var_c(jind), &
1157 c_var_a(jind), &
1158 "offset number at", &
1159 c_var_b(jind) - first + 1, &
1160 c_var_d(jind) - first + 1, &
1161 c_var_c(jind) - first + 1, &
1162 c_var_a(jind) - first + 1
1163 END IF
1164 END DO
1165 CALL set_molecule_kind(molecule_kind=molecule_kind, &
1166 nimpr=SIZE(impr_list), impr_list=impr_list)
1167 CALL set_molecule_kind(molecule_kind=molecule_kind, &
1168 nopbend=SIZE(opbend_list), opbend_list=opbend_list)
1169 END DO
1170 cpassert(nval_tot1 == nval_tot2)
1171 DEALLOCATE (map_var_mol)
1172 DEALLOCATE (map_cvar_mol)
1173 DEALLOCATE (map_vars)
1174 DEALLOCATE (map_cvars)
1175 DEALLOCATE (bnd_type)
1176 DEALLOCATE (bnd_ctype)
1177 DEALLOCATE (c_var_a)
1178 DEALLOCATE (c_var_b)
1179 DEALLOCATE (c_var_c)
1180 DEALLOCATE (c_var_d)
1181 IF (found) THEN
1182 DEALLOCATE (c_var_type)
1183 END IF
1184 CALL timestop(handle2)
1185 !-----------------------------------------------------------------------------
1186 !-----------------------------------------------------------------------------
1187 ! Finally deallocate some stuff.
1188 !-----------------------------------------------------------------------------
1189 DEALLOCATE (first_list)
1190 DEALLOCATE (last_list)
1191 DEALLOCATE (map_atom_mol)
1192 DEALLOCATE (map_atom_type)
1193 CALL timestop(handle)
1194 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1195 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
1196 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 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 ALLOCATE (map_vars(nvar1))
1217 CALL sort(map_var_mol, nvar1, map_vars)
1218 ALLOCATE (bnd_type(2, nmol_type))
1219 bnd_type = 0
1220 IF (nvar1 == 0) RETURN
1221 DO j = 1, nvar1
1222 IF (map_var_mol(j) /= -1) EXIT
1223 END DO
1224 IF (j == nvar1 + 1) RETURN
1225 i = map_var_mol(j)
1226 bnd_type(1, i) = j
1227 DO ibond = j, nvar1
1228 IF (map_var_mol(ibond) /= i) THEN
1229 bnd_type(2, i) = ibond - 1
1230 i = map_var_mol(ibond)
1231 bnd_type(1, i) = ibond
1232 END IF
1233 END DO
1234 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 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 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
1252 TYPE(connectivity_info_type), POINTER :: conn_info
1253
1254 NULLIFY (multiple_unit_cell)
1255 CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", &
1256 i_vals=multiple_unit_cell)
1257 IF (any(multiple_unit_cell /= 1)) THEN
1258 fac = product(multiple_unit_cell)
1259 conn_info => topology%conn_info
1260
1261 nbond = 0
1262 IF (ASSOCIATED(conn_info%bond_a)) THEN
1263 nbond = SIZE(conn_info%bond_a)
1264 CALL reallocate(conn_info%bond_a, 1, nbond*fac)
1265 CALL reallocate(conn_info%bond_b, 1, nbond*fac)
1266 END IF
1267
1268 ntheta = 0
1269 IF (ASSOCIATED(conn_info%theta_a)) THEN
1270 ntheta = SIZE(conn_info%theta_a)
1271 CALL reallocate(conn_info%theta_a, 1, ntheta*fac)
1272 CALL reallocate(conn_info%theta_b, 1, ntheta*fac)
1273 CALL reallocate(conn_info%theta_c, 1, ntheta*fac)
1274 END IF
1275
1276 nphi = 0
1277 IF (ASSOCIATED(conn_info%phi_a)) THEN
1278 nphi = SIZE(conn_info%phi_a)
1279 CALL reallocate(conn_info%phi_a, 1, nphi*fac)
1280 CALL reallocate(conn_info%phi_b, 1, nphi*fac)
1281 CALL reallocate(conn_info%phi_c, 1, nphi*fac)
1282 CALL reallocate(conn_info%phi_d, 1, nphi*fac)
1283 END IF
1284
1285 nimpr = 0
1286 IF (ASSOCIATED(conn_info%impr_a)) THEN
1287 nimpr = SIZE(conn_info%impr_a)
1288 CALL reallocate(conn_info%impr_a, 1, nimpr*fac)
1289 CALL reallocate(conn_info%impr_b, 1, nimpr*fac)
1290 CALL reallocate(conn_info%impr_c, 1, nimpr*fac)
1291 CALL reallocate(conn_info%impr_d, 1, nimpr*fac)
1292 END IF
1293
1294 nbond_c = 0
1295 IF (ASSOCIATED(conn_info%c_bond_a)) THEN
1296 nbond_c = SIZE(conn_info%c_bond_a)
1297 CALL reallocate(conn_info%c_bond_a, 1, nbond_c*fac)
1298 CALL reallocate(conn_info%c_bond_b, 1, nbond_c*fac)
1299 END IF
1300
1301 nub = 0
1302 IF (ASSOCIATED(conn_info%ub_a)) THEN
1303 nub = SIZE(conn_info%ub_a)
1304 CALL reallocate(conn_info%ub_a, 1, nub*fac)
1305 CALL reallocate(conn_info%ub_b, 1, nub*fac)
1306 CALL reallocate(conn_info%ub_c, 1, nub*fac)
1307 END IF
1308
1309 nonfo = 0
1310 IF (ASSOCIATED(conn_info%onfo_a)) THEN
1311 nonfo = SIZE(conn_info%onfo_a)
1312 CALL reallocate(conn_info%onfo_a, 1, nonfo*fac)
1313 CALL reallocate(conn_info%onfo_b, 1, nonfo*fac)
1314 END IF
1315
1316 natoms_orig = topology%natoms/fac
1317 ind = 0
1318 DO k = 1, multiple_unit_cell(3)
1319 DO j = 1, multiple_unit_cell(2)
1320 DO i = 1, multiple_unit_cell(1)
1321 ind = ind + 1
1322 IF (ind == 1) cycle
1323 a = (ind - 1)*natoms_orig
1324
1325 ! Bonds
1326 IF (nbond > 0) THEN
1327 m = (ind - 1)*nbond
1328 conn_info%bond_a(m + 1:m + nbond) = conn_info%bond_a(1:nbond) + a
1329 conn_info%bond_b(m + 1:m + nbond) = conn_info%bond_b(1:nbond) + a
1330 END IF
1331 ! Theta
1332 IF (ntheta > 0) THEN
1333 m = (ind - 1)*ntheta
1334 conn_info%theta_a(m + 1:m + ntheta) = conn_info%theta_a(1:ntheta) + a
1335 conn_info%theta_b(m + 1:m + ntheta) = conn_info%theta_b(1:ntheta) + a
1336 conn_info%theta_c(m + 1:m + ntheta) = conn_info%theta_c(1:ntheta) + a
1337 END IF
1338 ! Phi
1339 IF (nphi > 0) THEN
1340 m = (ind - 1)*nphi
1341 conn_info%phi_a(m + 1:m + nphi) = conn_info%phi_a(1:nphi) + a
1342 conn_info%phi_b(m + 1:m + nphi) = conn_info%phi_b(1:nphi) + a
1343 conn_info%phi_c(m + 1:m + nphi) = conn_info%phi_c(1:nphi) + a
1344 conn_info%phi_d(m + 1:m + nphi) = conn_info%phi_d(1:nphi) + a
1345 END IF
1346 ! Impropers
1347 IF (nimpr > 0) THEN
1348 m = (ind - 1)*nimpr
1349 conn_info%impr_a(m + 1:m + nimpr) = conn_info%impr_a(1:nimpr) + a
1350 conn_info%impr_b(m + 1:m + nimpr) = conn_info%impr_b(1:nimpr) + a
1351 conn_info%impr_c(m + 1:m + nimpr) = conn_info%impr_c(1:nimpr) + a
1352 conn_info%impr_d(m + 1:m + nimpr) = conn_info%impr_d(1:nimpr) + a
1353 END IF
1354 ! Para_res
1355 IF (nbond_c > 0) THEN
1356 m = (ind - 1)*nbond_c
1357 conn_info%c_bond_a(m + 1:m + nbond_c) = conn_info%c_bond_a(1:nbond_c) + a
1358 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 IF (nub > 0) THEN
1362 m = (ind - 1)*nub
1363 conn_info%ub_a(m + 1:m + nub) = conn_info%ub_a(1:nub) + a
1364 conn_info%ub_b(m + 1:m + nub) = conn_info%ub_b(1:nub) + a
1365 conn_info%ub_c(m + 1:m + nub) = conn_info%ub_c(1:nub) + a
1366 END IF
1367 ! ONFO
1368 IF (nonfo > 0) THEN
1369 m = (ind - 1)*nonfo
1370 conn_info%onfo_a(m + 1:m + nonfo) = conn_info%onfo_a(1:nonfo) + a
1371 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 END SUBROUTINE topology_conn_multiple
1379
1380END MODULE topology_connectivity_util
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
memory_utilities::reallocate
Definition memory_utilities.F:27
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
force_field_kind_types
Define all structure types related to force field kinds.
Definition force_field_kind_types.F:14
force_field_kind_types::do_ff_charmm
integer, parameter, public do_ff_charmm
Definition force_field_kind_types.F:25
force_field_kind_types::do_ff_harmonic
integer, parameter, public do_ff_harmonic
Definition force_field_kind_types.F:25
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_conn_g87
integer, parameter, public do_conn_g87
Definition input_constants.F:412
input_constants::do_conn_user
integer, parameter, public do_conn_user
Definition input_constants.F:412
input_constants::do_conn_g96
integer, parameter, public do_conn_g96
Definition input_constants.F:412
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
molecule_kind_types::set_molecule_kind
subroutine, public set_molecule_kind(molecule_kind, name, mass, charge, kind_number, molecule_list, atom_list, nbond, bond_list, nbend, bend_list, nub, ub_list, nimpr, impr_list, nopbend, opbend_list, ntorsion, torsion_list, fixd_list, ncolv, colv_list, ng3x3, g3x3_list, ng4x6, nfixd, g4x6_list, nvsite, vsite_list, ng3x3_restraint, ng4x6_restraint, nfixd_restraint, nshell, shell_list, nvsite_restraint, bond_kind_set, bend_kind_set, ub_kind_set, torsion_kind_set, impr_kind_set, opbend_kind_set, nelectron, nsgf, molname_generated)
Set the components of a molecule kind.
Definition molecule_kind_types.F:779
molecule_kind_types::allocate_molecule_kind_set
subroutine, public allocate_molecule_kind_set(molecule_kind_set, nmolecule_kind)
Allocate and initialize a molecule kind set.
Definition molecule_kind_types.F:323
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
molecule_types::allocate_molecule_set
subroutine, public allocate_molecule_set(molecule_set, nmolecule)
Allocate a molecule set.
Definition molecule_types.F:202
molecule_types::set_molecule_set
subroutine, public set_molecule_set(molecule_set, first_atom, last_atom)
Set a molecule data set.
Definition molecule_types.F:386
molecule_types::set_molecule
subroutine, public set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
Set a molecule data set.
Definition molecule_types.F:337
string_table
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition string_table.F:22
string_table::id2str
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition string_table.F:115
topology_connectivity_util
Collection of subroutine needed for topology related things.
Definition topology_connectivity_util.F:14
topology_connectivity_util::topology_conn_multiple
subroutine, public topology_conn_multiple(topology, subsys_section)
Handles the multiple unit cell option for the connectivity.
Definition topology_connectivity_util.F:1245
topology_connectivity_util::topology_connectivity_pack
subroutine, public topology_connectivity_pack(molecule_kind_set, molecule_set, topology, subsys_section)
topology connectivity pack
Definition topology_connectivity_util.F:65
topology_types
Definition topology_types.F:18
topology
Control for reading in different topologies and coordinates.
Definition topology.F:13
util
All kind of helpful little routines.
Definition util.F:14
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
molecule_kind_types::atom_type
Definition molecule_kind_types.F:56
molecule_kind_types::bend_type
Definition molecule_kind_types.F:76
molecule_kind_types::bond_type
Definition molecule_kind_types.F:69
molecule_kind_types::impr_type
Definition molecule_kind_types.F:97
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_kind_types::opbend_type
Definition molecule_kind_types.F:104
molecule_kind_types::torsion_type
Definition molecule_kind_types.F:90
molecule_kind_types::ub_type
Definition molecule_kind_types.F:83
molecule_types::local_states_type
Definition molecule_types.F:80
molecule_types::molecule_type
Definition molecule_types.F:110
topology_types::atom_info_type
Definition topology_types.F:35
topology_types::connectivity_info_type
Definition topology_types.F:53
topology_types::topology_parameters_type
Definition topology_types.F:142