(git:b195825)
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 
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  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 
1380 END 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
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
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