(git:e7e05ae)
topology_constraint_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 MODULE topology_constraint_util
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
15  get_atomic_kind,&
16  is_hydrogen
17  USE cell_types, ONLY: use_perd_x,&
18  use_perd_xy,&
19  use_perd_xyz,&
20  use_perd_xz,&
21  use_perd_y,&
22  use_perd_yz,&
23  use_perd_z
24  USE colvar_methods, ONLY: colvar_eval_mol_f
25  USE colvar_types, ONLY: &
26  colvar_clone, colvar_counters, colvar_create, colvar_p_reallocate, colvar_release, &
27  colvar_setup, colvar_type, dist_colvar_id, torsion_colvar_id, xyz_diag_colvar_id, &
28  xyz_outerdiag_colvar_id
29  USE colvar_utils, ONLY: post_process_colvar
30  USE cp_log_handling, ONLY: cp_get_default_logger,&
31  cp_logger_type,&
32  cp_to_string
33  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34  cp_print_key_unit_nr
35  USE input_constants, ONLY: do_constr_atomic,&
36  do_constr_molec
37  USE input_section_types, ONLY: section_vals_get,&
38  section_vals_get_subs_vals,&
39  section_vals_type,&
40  section_vals_val_get,&
41  section_vals_val_set
42  USE kinds, ONLY: default_string_length,&
43  dp
44  USE memory_utilities, ONLY: reallocate
45  USE molecule_kind_types, ONLY: &
46  atom_type, bond_type, colvar_constraint_type, fixd_constraint_type, g3x3_constraint_type, &
47  g4x6_constraint_type, get_molecule_kind, molecule_kind_type, set_molecule_kind, &
48  setup_colvar_counters, vsite_constraint_type
49  USE molecule_types, ONLY: get_molecule,&
50  global_constraint_type,&
51  local_colvar_constraint_type,&
52  local_constraint_type,&
53  local_g3x3_constraint_type,&
54  local_g4x6_constraint_type,&
55  molecule_type,&
56  set_molecule
57  USE particle_types, ONLY: particle_type
58  USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
59  USE qmmm_types_low, ONLY: qmmm_env_mm_type
60  USE topology_types, ONLY: constr_list_type,&
61  constraint_info_type,&
62  topology_parameters_type
63 #include "./base/base_uses.f90"
64 
65  IMPLICIT NONE
66 
67  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_constraint_util'
68 
69  PRIVATE
70  PUBLIC :: topology_constraint_pack
71 
72 CONTAINS
73 
74 ! **************************************************************************************************
75 !> \brief Pack in all the information needed for the constraints
76 !> \param molecule_kind_set ...
77 !> \param molecule_set ...
78 !> \param topology ...
79 !> \param qmmm_env ...
80 !> \param particle_set ...
81 !> \param input_file ...
82 !> \param subsys_section ...
83 !> \param gci ...
84 ! **************************************************************************************************
85  SUBROUTINE topology_constraint_pack(molecule_kind_set, molecule_set, &
86  topology, qmmm_env, particle_set, input_file, subsys_section, gci)
87  TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
88  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
89  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
90  TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
91  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
92  TYPE(section_vals_type), POINTER :: input_file, subsys_section
93  TYPE(global_constraint_type), POINTER :: gci
94 
95  CHARACTER(len=*), PARAMETER :: routinen = 'topology_constraint_pack'
96 
97  CHARACTER(LEN=2) :: element_symbol
98  CHARACTER(LEN=default_string_length) :: molname, name
99  CHARACTER(LEN=default_string_length), &
100  DIMENSION(:), POINTER :: atom_typeh, cnds
101  INTEGER :: cind, first, first_atom, gind, handle, handle2, i, ii, itype, iw, j, k, k1loc, &
102  k2loc, kk, last, last_atom, m, n_start_colv, natom, nbond, ncolv_glob, ncolv_mol, &
103  nfixd_list_gci, nfixd_restart, nfixd_restraint, nfixed_atoms, ng3x3, ng3x3_restraint, &
104  ng4x6, ng4x6_restraint, nhdist, nmolecule, nrep, nvsite, nvsite_restraint, offset
105  INTEGER, DIMENSION(:), POINTER :: constr_x_glob, inds, molecule_list
106  LOGICAL :: exclude_mm, exclude_qm, fix_atom_mm, fix_atom_molname, fix_atom_qm, &
107  fix_atom_qmmm, fix_fixed_atom, found_molname, is_qm, ishbond, ldummy, &
108  restart_restraint_clv, restart_restraint_pos, use_clv_info
109  LOGICAL, ALLOCATABLE, DIMENSION(:) :: missed_molname
110  REAL(kind=dp) :: rmod, rvec(3)
111  REAL(kind=dp), DIMENSION(:), POINTER :: hdist, r
112  TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
113  TYPE(atomic_kind_type), POINTER :: atomic_kind
114  TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
115  TYPE(colvar_constraint_type), DIMENSION(:), &
116  POINTER :: colv_list
117  TYPE(colvar_counters) :: ncolv
118  TYPE(constr_list_type), DIMENSION(:), POINTER :: constr_x_mol
119  TYPE(constraint_info_type), POINTER :: cons_info
120  TYPE(cp_logger_type), POINTER :: logger
121  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list, fixd_list_gci
122  TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
123  TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
124  TYPE(local_colvar_constraint_type), DIMENSION(:), &
125  POINTER :: lcolv
126  TYPE(local_constraint_type), POINTER :: lci
127  TYPE(local_g3x3_constraint_type), DIMENSION(:), &
128  POINTER :: lg3x3
129  TYPE(local_g4x6_constraint_type), DIMENSION(:), &
130  POINTER :: lg4x6
131  TYPE(molecule_kind_type), POINTER :: molecule_kind
132  TYPE(molecule_type), POINTER :: molecule
133  TYPE(section_vals_type), POINTER :: colvar_func_info, colvar_rest, &
134  fixd_restr_rest, hbonds_section
135  TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
136 
137  NULLIFY (logger, constr_x_mol, constr_x_glob)
138  logger => cp_get_default_logger()
139  iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
140  extension=".subsysLog")
141  CALL timeset(routinen, handle)
142  CALL timeset(routinen//"_1", handle2)
143 
144  cons_info => topology%cons_info
145  hbonds_section => section_vals_get_subs_vals(input_file, &
146  "MOTION%CONSTRAINT%HBONDS")
147  fixd_restr_rest => section_vals_get_subs_vals(input_file, &
148  "MOTION%CONSTRAINT%FIX_ATOM_RESTART")
149  CALL section_vals_get(fixd_restr_rest, explicit=restart_restraint_pos)
150  colvar_rest => section_vals_get_subs_vals(input_file, &
151  "MOTION%CONSTRAINT%COLVAR_RESTART")
152  CALL section_vals_get(colvar_rest, explicit=restart_restraint_clv)
153  colvar_func_info => section_vals_get_subs_vals(subsys_section, &
154  "COLVAR%COLVAR_FUNC_INFO")
155  CALL section_vals_get(colvar_func_info, explicit=use_clv_info)
156  !-----------------------------------------------------------------------------
157  !-----------------------------------------------------------------------------
158  ! 1. NULLIFY the molecule_set(imol)%lci via set_molecule_set
159  !-----------------------------------------------------------------------------
160  DO i = 1, topology%nmol
161  molecule => molecule_set(i)
162  NULLIFY (lci)
163  ! only allocate the lci if constraints are active. Can this stuff be distributed ?
164  IF (topology%const_atom .OR. topology%const_hydr .OR. &
165  topology%const_33 .OR. topology%const_46 .OR. &
166  topology%const_colv .OR. topology%const_vsite) THEN
167  ALLOCATE (lci)
168  NULLIFY (lci%lcolv)
169  NULLIFY (lci%lg3x3)
170  NULLIFY (lci%lg4x6)
171  END IF
172  CALL set_molecule(molecule, lci=lci)
173  END DO
174  ALLOCATE (gci)
175  NULLIFY (gci%lcolv, &
176  gci%lg3x3, &
177  gci%lg4x6, &
178  gci%fixd_list, &
179  gci%colv_list, &
180  gci%g3x3_list, &
181  gci%g4x6_list, &
182  gci%vsite_list)
183  gci%ntot = 0
184  gci%ng3x3 = 0
185  gci%ng4x6 = 0
186  gci%nvsite = 0
187  gci%ng3x3_restraint = 0
188  gci%ng4x6_restraint = 0
189  gci%nvsite_restraint = 0
190  CALL setup_colvar_counters(gci%colv_list, gci%ncolv)
191  gci%nrestraint = gci%ng3x3_restraint + &
192  gci%ng4x6_restraint + &
193  gci%nvsite_restraint + &
194  gci%ncolv%nrestraint
195  CALL timestop(handle2)
196  CALL timeset(routinen//"_2", handle2)
197  !-----------------------------------------------------------------------------
198  !-----------------------------------------------------------------------------
199  ! 2. Add more stuff to COLVAR constraint if constraint hydrogen is on
200  !-----------------------------------------------------------------------------
201  IF (topology%const_hydr) THEN
202  topology%const_colv = .true.
203  NULLIFY (atom_typeh, hdist)
204  ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
205  DO i = 1, SIZE(molecule_kind_set)
206  ALLOCATE (constr_x_mol(i)%constr(1))
207  constr_x_mol(i)%constr(1) = 1
208  END DO
209  CALL section_vals_val_get(hbonds_section, "MOLECULE", n_rep_val=nrep)
210  IF (nrep /= 0) THEN
211  NULLIFY (inds)
212  DO i = 1, SIZE(molecule_kind_set)
213  constr_x_mol(i)%constr(1) = 0
214  END DO
215  CALL section_vals_val_get(hbonds_section, "MOLECULE", i_vals=inds)
216  DO i = 1, SIZE(inds)
217  constr_x_mol(inds(i))%constr(1) = 1
218  END DO
219  ELSE
220  CALL section_vals_val_get(hbonds_section, "MOLNAME", n_rep_val=nrep)
221  IF (nrep /= 0) THEN
222  NULLIFY (cnds)
223  DO i = 1, SIZE(molecule_kind_set)
224  constr_x_mol(i)%constr(1) = 0
225  END DO
226  CALL section_vals_val_get(hbonds_section, "MOLNAME", c_vals=cnds)
227  DO i = 1, SIZE(cnds)
228  found_molname = .false.
229  DO k = 1, SIZE(molecule_kind_set)
230  molecule_kind => molecule_kind_set(k)
231  name = molecule_kind%name
232  ldummy = qmmm_ff_precond_only_qm(id1=name)
233  IF (cnds(i) == name) THEN
234  constr_x_mol(k)%constr(1) = 1
235  found_molname = .true.
236  END IF
237  END DO
238  CALL print_warning_molname(found_molname, cnds(i))
239  END DO
240  END IF
241  END IF
242  CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", n_rep_val=nrep)
243  IF (nrep /= 0) &
244  CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", c_vals=atom_typeh)
245  CALL section_vals_val_get(hbonds_section, "TARGETS", n_rep_val=nrep)
246  IF (nrep /= 0) &
247  CALL section_vals_val_get(hbonds_section, "TARGETS", r_vals=hdist)
248  IF (ASSOCIATED(hdist)) THEN
249  cpassert(SIZE(hdist) == SIZE(atom_typeh))
250  END IF
251  CALL section_vals_val_get(hbonds_section, "exclude_qm", l_val=exclude_qm)
252  CALL section_vals_val_get(hbonds_section, "exclude_mm", l_val=exclude_mm)
253  nhdist = 0
254  DO i = 1, SIZE(molecule_kind_set)
255  molecule_kind => molecule_kind_set(i)
256  IF (constr_x_mol(i)%constr(1) == 0) cycle
257  CALL get_molecule_kind(molecule_kind=molecule_kind, &
258  bond_list=bond_list, nbond=nbond, atom_list=atom_list, &
259  molecule_list=molecule_list)
260  ! Let's tag all requested atoms involving Hydrogen
261  ! on the first molecule of this kind
262  molecule => molecule_set(molecule_list(1))
263  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
264  natom = last_atom - first_atom + 1
265  DO k = 1, nbond
266  ishbond = .false.
267  j = bond_list(k)%a
268  IF (j < 1 .OR. j > natom) cycle
269  atomic_kind => atom_list(j)%atomic_kind
270  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
271  is_qm = qmmm_ff_precond_only_qm(id1=name)
272  IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
273  IF (is_qm .AND. exclude_qm) ishbond = .false.
274  IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
275  IF (.NOT. ishbond) THEN
276  j = bond_list(k)%b
277  IF (j < 1 .OR. j > natom) cycle
278  atomic_kind => atom_list(j)%atomic_kind
279  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
280  is_qm = qmmm_ff_precond_only_qm(id1=name)
281  IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
282  IF (is_qm .AND. exclude_qm) ishbond = .false.
283  IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
284  END IF
285  IF (ishbond) THEN
286  nhdist = nhdist + 1
287  END IF
288  END DO
289  END DO
290  n_start_colv = cons_info%nconst_colv
291  cons_info%nconst_colv = nhdist + n_start_colv
292  CALL reallocate(cons_info%const_colv_mol, 1, cons_info%nconst_colv)
293  CALL reallocate(cons_info%const_colv_molname, 1, cons_info%nconst_colv)
294  CALL reallocate(cons_info%const_colv_target, 1, cons_info%nconst_colv)
295  CALL reallocate(cons_info%const_colv_target_growth, 1, cons_info%nconst_colv)
296  CALL colvar_p_reallocate(cons_info%colvar_set, 1, cons_info%nconst_colv)
297  ! Fill in Restraints info
298  CALL reallocate(cons_info%colv_intermolecular, 1, cons_info%nconst_colv)
299  CALL reallocate(cons_info%colv_restraint, 1, cons_info%nconst_colv)
300  CALL reallocate(cons_info%colv_k0, 1, cons_info%nconst_colv)
301  CALL reallocate(cons_info%colv_exclude_qm, 1, cons_info%nconst_colv)
302  CALL reallocate(cons_info%colv_exclude_mm, 1, cons_info%nconst_colv)
303  ! Bonds involving hydrogens are by their nature only intramolecular
304  cons_info%colv_intermolecular(n_start_colv + 1:cons_info%nconst_colv) = .false.
305  cons_info%colv_exclude_qm(n_start_colv + 1:cons_info%nconst_colv) = .false.
306  cons_info%colv_exclude_mm(n_start_colv + 1:cons_info%nconst_colv) = .false.
307  cons_info%colv_restraint(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_restraint
308  cons_info%colv_k0(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_k0
309  !
310  nhdist = 0
311  DO i = 1, SIZE(molecule_kind_set)
312  IF (constr_x_mol(i)%constr(1) == 0) cycle
313  molecule_kind => molecule_kind_set(i)
314  CALL get_molecule_kind(molecule_kind=molecule_kind, &
315  bond_list=bond_list, nbond=nbond, atom_list=atom_list, &
316  molecule_list=molecule_list)
317  molecule => molecule_set(molecule_list(1))
318  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
319  natom = last_atom - first_atom + 1
320  offset = first_atom - 1
321  DO k = 1, nbond
322  ishbond = .false.
323  j = bond_list(k)%a
324  IF (j < 1 .OR. j > natom) cycle
325  atomic_kind => atom_list(j)%atomic_kind
326  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
327  is_qm = qmmm_ff_precond_only_qm(id1=name)
328  IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
329  IF (is_qm .AND. exclude_qm) ishbond = .false.
330  IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
331  IF (.NOT. ishbond) THEN
332  j = bond_list(k)%b
333  IF (j < 1 .OR. j > natom) cycle
334  atomic_kind => atom_list(j)%atomic_kind
335  CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
336  is_qm = qmmm_ff_precond_only_qm(id1=name)
337  IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
338  IF (is_qm .AND. exclude_qm) ishbond = .false.
339  IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
340  END IF
341  IF (ishbond) THEN
342  nhdist = nhdist + 1
343  rvec = particle_set(offset + bond_list(k)%a)%r - particle_set(offset + bond_list(k)%b)%r
344  rmod = sqrt(dot_product(rvec, rvec))
345  IF (ASSOCIATED(hdist)) THEN
346  IF (SIZE(hdist) > 0) THEN
347  IF (bond_list(k)%a == j) atomic_kind => atom_list(bond_list(k)%b)%atomic_kind
348  IF (bond_list(k)%b == j) atomic_kind => atom_list(bond_list(k)%a)%atomic_kind
349  CALL get_atomic_kind(atomic_kind=atomic_kind, &
350  name=name, element_symbol=element_symbol)
351  ldummy = qmmm_ff_precond_only_qm(id1=name)
352  DO m = 1, SIZE(hdist)
353  IF (trim(name) == trim(atom_typeh(m))) EXIT
354  IF (trim(element_symbol) == trim(atom_typeh(m))) EXIT
355  END DO
356  IF (m <= SIZE(hdist)) THEN
357  rmod = hdist(m)
358  END IF
359  END IF
360  END IF
361  cons_info%const_colv_mol(nhdist + n_start_colv) = i
362  cons_info%const_colv_molname(nhdist + n_start_colv) = "UNDEF"
363  cons_info%const_colv_target(nhdist + n_start_colv) = rmod
364  cons_info%const_colv_target_growth(nhdist + n_start_colv) = 0.0_dp
365  CALL colvar_create(cons_info%colvar_set(nhdist + n_start_colv)%colvar, &
366  dist_colvar_id)
367  cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%i_at = bond_list(k)%a
368  cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%j_at = bond_list(k)%b
369  CALL colvar_setup(cons_info%colvar_set(nhdist + n_start_colv)%colvar)
370  END IF
371  END DO
372  END DO
373  DO j = 1, SIZE(constr_x_mol)
374  DEALLOCATE (constr_x_mol(j)%constr)
375  END DO
376  DEALLOCATE (constr_x_mol)
377  END IF
378 
379  CALL timestop(handle2)
380  CALL timeset(routinen//"_3", handle2)
381  !-----------------------------------------------------------------------------
382  !-----------------------------------------------------------------------------
383  ! 3. Set the COLVAR constraint molecule_kind_set(ikind)%colv_list
384  !-----------------------------------------------------------------------------
385  IF (topology%const_colv) THEN
386  ! Post Process of COLVARS..
387  DO ii = 1, SIZE(cons_info%colvar_set)
388  CALL post_process_colvar(cons_info%colvar_set(ii)%colvar, particle_set)
389  END DO
390  ! Real constraint/restraint part..
391  CALL give_constraint_array(cons_info%const_colv_mol, &
392  cons_info%const_colv_molname, &
393  cons_info%colv_intermolecular, &
394  constr_x_mol, &
395  constr_x_glob, &
396  molecule_kind_set, &
397  cons_info%colv_exclude_qm, &
398  cons_info%colv_exclude_mm)
399  ! Intramolecular constraints
400  gind = 0
401  cind = 0
402  DO ii = 1, SIZE(molecule_kind_set)
403  molecule_kind => molecule_kind_set(ii)
404  CALL get_molecule_kind(molecule_kind=molecule_kind, &
405  nmolecule=nmolecule, molecule_list=molecule_list)
406  ncolv_mol = SIZE(constr_x_mol(ii)%constr)
407  ALLOCATE (colv_list(ncolv_mol))
408  ! Starting index of the first molecule of this kind.
409  ! We need the index if no target is provided in the input file
410  ! for the collective variable.. The target will be computed on the
411  ! first molecule of the kind...
412  molecule => molecule_set(molecule_list(1))
413  CALL get_molecule(molecule, first_atom=first_atom)
414  CALL setup_colv_list(colv_list, constr_x_mol(ii)%constr, gind, &
415  cons_info, topology, particle_set, restart_restraint_clv, &
416  colvar_rest, first_atom)
417  CALL setup_colvar_counters(colv_list, ncolv)
418  CALL set_molecule_kind(molecule_kind, colv_list=colv_list, ncolv=ncolv)
419  DO j = 1, nmolecule
420  molecule => molecule_set(molecule_list(j))
421  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
422  ALLOCATE (lcolv(ncolv_mol))
423  CALL setup_lcolv(lcolv, constr_x_mol(ii)%constr, first_atom, last_atom, &
424  cons_info, particle_set, colvar_func_info, use_clv_info, cind)
425  CALL set_molecule(molecule=molecule, lcolv=lcolv)
426  END DO
427  END DO
428  DO j = 1, SIZE(constr_x_mol)
429  DEALLOCATE (constr_x_mol(j)%constr)
430  END DO
431  DEALLOCATE (constr_x_mol)
432  ! Intermolecular constraints
433  ncolv_glob = 0
434  IF (ASSOCIATED(constr_x_glob)) THEN
435  ncolv_glob = SIZE(constr_x_glob)
436  ALLOCATE (colv_list(ncolv_glob))
437  CALL setup_colv_list(colv_list, constr_x_glob, gind, cons_info, &
438  topology, particle_set, restart_restraint_clv, colvar_rest, &
439  first_atom=1)
440  CALL setup_colvar_counters(colv_list, ncolv)
441  ALLOCATE (lcolv(ncolv_glob))
442  CALL setup_lcolv(lcolv, constr_x_glob, 1, SIZE(particle_set), cons_info, &
443  particle_set, colvar_func_info, use_clv_info, cind)
444  gci%colv_list => colv_list
445  gci%lcolv => lcolv
446  gci%ncolv = ncolv
447  ! Total number of Intermolecular constraints
448  gci%ntot = gci%ncolv%ntot + gci%ntot
449  DEALLOCATE (constr_x_glob)
450  END IF
451  END IF
452 
453  CALL timestop(handle2)
454  CALL timeset(routinen//"_4", handle2)
455  !-----------------------------------------------------------------------------
456  !-----------------------------------------------------------------------------
457  ! 4. Set the group 3x3 constraint g3x3_list
458  !-----------------------------------------------------------------------------
459  IF (topology%const_33) THEN
460  CALL give_constraint_array(cons_info%const_g33_mol, &
461  cons_info%const_g33_molname, &
462  cons_info%g33_intermolecular, &
463  constr_x_mol, &
464  constr_x_glob, &
465  molecule_kind_set, &
466  cons_info%g33_exclude_qm, &
467  cons_info%g33_exclude_mm)
468  ! Intramolecular constraints
469  DO ii = 1, SIZE(molecule_kind_set)
470  molecule_kind => molecule_kind_set(ii)
471  CALL get_molecule_kind(molecule_kind=molecule_kind, &
472  nmolecule=nmolecule, &
473  molecule_list=molecule_list)
474  ng3x3 = SIZE(constr_x_mol(ii)%constr)
475  ALLOCATE (g3x3_list(ng3x3))
476  CALL setup_g3x3_list(g3x3_list, constr_x_mol(ii)%constr, cons_info, ng3x3_restraint)
477  CALL set_molecule_kind(molecule_kind, ng3x3=ng3x3, ng3x3_restraint=ng3x3_restraint, g3x3_list=g3x3_list)
478  DO j = 1, nmolecule
479  molecule => molecule_set(molecule_list(j))
480  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
481  ALLOCATE (lg3x3(ng3x3))
482  CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
483  CALL set_molecule(molecule=molecule, lg3x3=lg3x3)
484  END DO
485  END DO
486  DO j = 1, SIZE(constr_x_mol)
487  DEALLOCATE (constr_x_mol(j)%constr)
488  END DO
489  DEALLOCATE (constr_x_mol)
490  ! Intermolecular constraints
491  IF (ASSOCIATED(constr_x_glob)) THEN
492  ng3x3 = SIZE(constr_x_glob)
493  ALLOCATE (g3x3_list(ng3x3))
494  CALL setup_g3x3_list(g3x3_list, constr_x_glob, cons_info, ng3x3_restraint)
495  ALLOCATE (lg3x3(ng3x3))
496  CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
497  gci%g3x3_list => g3x3_list
498  gci%lg3x3 => lg3x3
499  gci%ng3x3 = ng3x3
500  gci%ng3x3_restraint = ng3x3_restraint
501  ! Total number of Intermolecular constraints
502  gci%ntot = 3*gci%ng3x3 + gci%ntot
503  DEALLOCATE (constr_x_glob)
504  END IF
505  END IF
506 
507  CALL timestop(handle2)
508  CALL timeset(routinen//"_5", handle2)
509  !-----------------------------------------------------------------------------
510  !-----------------------------------------------------------------------------
511  ! 5. Set the group 4x6 constraint g4x6_list
512  !-----------------------------------------------------------------------------
513  IF (topology%const_46) THEN
514  CALL give_constraint_array(cons_info%const_g46_mol, &
515  cons_info%const_g46_molname, &
516  cons_info%g46_intermolecular, &
517  constr_x_mol, &
518  constr_x_glob, &
519  molecule_kind_set, &
520  cons_info%g46_exclude_qm, &
521  cons_info%g46_exclude_mm)
522  ! Intramolecular constraints
523  DO ii = 1, SIZE(molecule_kind_set)
524  molecule_kind => molecule_kind_set(ii)
525  CALL get_molecule_kind(molecule_kind=molecule_kind, &
526  nmolecule=nmolecule, molecule_list=molecule_list)
527  ng4x6 = SIZE(constr_x_mol(ii)%constr)
528  ALLOCATE (g4x6_list(ng4x6))
529  CALL setup_g4x6_list(g4x6_list, constr_x_mol(ii)%constr, cons_info, ng4x6_restraint)
530  CALL set_molecule_kind(molecule_kind, ng4x6=ng4x6, ng4x6_restraint=ng4x6_restraint, g4x6_list=g4x6_list)
531  DO j = 1, nmolecule
532  molecule => molecule_set(molecule_list(j))
533  CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
534  ALLOCATE (lg4x6(ng4x6))
535  CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
536  CALL set_molecule(molecule=molecule, lg4x6=lg4x6)
537  END DO
538  END DO
539  DO j = 1, SIZE(constr_x_mol)
540  DEALLOCATE (constr_x_mol(j)%constr)
541  END DO
542  DEALLOCATE (constr_x_mol)
543  ! Intermolecular constraints
544  IF (ASSOCIATED(constr_x_glob)) THEN
545  ng4x6 = SIZE(constr_x_glob)
546  ALLOCATE (g4x6_list(ng4x6))
547  CALL setup_g4x6_list(g4x6_list, constr_x_glob, cons_info, ng4x6_restraint)
548  ALLOCATE (lg4x6(ng4x6))
549  CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
550  gci%g4x6_list => g4x6_list
551  gci%lg4x6 => lg4x6
552  gci%ng4x6 = ng4x6
553  gci%ng4x6_restraint = ng4x6_restraint
554  ! Total number of Intermolecular constraints
555  gci%ntot = 6*gci%ng4x6 + gci%ntot
556  DEALLOCATE (constr_x_glob)
557  END IF
558  END IF
559 
560  CALL timestop(handle2)
561  CALL timeset(routinen//"_6", handle2)
562  !-----------------------------------------------------------------------------
563  !-----------------------------------------------------------------------------
564  ! 6. Set the group vsite constraint vsite_list
565  !-----------------------------------------------------------------------------
566  IF (topology%const_vsite) THEN
567  CALL give_constraint_array(cons_info%const_vsite_mol, &
568  cons_info%const_vsite_molname, &
569  cons_info%vsite_intermolecular, &
570  constr_x_mol, &
571  constr_x_glob, &
572  molecule_kind_set, &
573  cons_info%vsite_exclude_qm, &
574  cons_info%vsite_exclude_mm)
575  ! Intramolecular constraints
576  DO ii = 1, SIZE(molecule_kind_set)
577  molecule_kind => molecule_kind_set(ii)
578  CALL get_molecule_kind(molecule_kind=molecule_kind, &
579  nmolecule=nmolecule, molecule_list=molecule_list)
580  nvsite = SIZE(constr_x_mol(ii)%constr)
581  ALLOCATE (vsite_list(nvsite))
582  CALL setup_vsite_list(vsite_list, constr_x_mol(ii)%constr, cons_info, nvsite_restraint)
583  CALL set_molecule_kind(molecule_kind, nvsite=nvsite, nvsite_restraint=nvsite_restraint, &
584  vsite_list=vsite_list)
585  END DO
586  DO j = 1, SIZE(constr_x_mol)
587  DEALLOCATE (constr_x_mol(j)%constr)
588  END DO
589  DEALLOCATE (constr_x_mol)
590  ! Intermolecular constraints
591  IF (ASSOCIATED(constr_x_glob)) THEN
592  nvsite = SIZE(constr_x_glob)
593  ALLOCATE (vsite_list(nvsite))
594  CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint)
595  gci%vsite_list => vsite_list
596  gci%nvsite = nvsite
597  gci%nvsite_restraint = nvsite_restraint
598  ! Total number of Intermolecular constraints
599  gci%ntot = gci%nvsite + gci%ntot
600  DEALLOCATE (constr_x_glob)
601  END IF
602  END IF
603  CALL timestop(handle2)
604  CALL timeset(routinen//"_7", handle2)
605  !-----------------------------------------------------------------------------
606  !-----------------------------------------------------------------------------
607  ! 7. Set the group fixed_atom constraint fixd_list
608  !-----------------------------------------------------------------------------
609  IF (topology%const_atom) THEN
610  ALLOCATE (fixd_list_gci(SIZE(particle_set)))
611  nfixd_list_gci = 0
612  ALLOCATE (missed_molname(SIZE(cons_info%fixed_molnames, 1)))
613  missed_molname = .true.
614  nfixd_restart = 0
615  DO i = 1, SIZE(molecule_kind_set)
616  molecule_kind => molecule_kind_set(i)
617  CALL get_molecule_kind(molecule_kind=molecule_kind, &
618  nmolecule=nmolecule, molecule_list=molecule_list, name=molname)
619  is_qm = qmmm_ff_precond_only_qm(id1=molname)
620  WHERE (molname .EQ. cons_info%fixed_molnames)
621  missed_molname = .false.
622  END WHERE
623  ! Try to figure out how many atoms of the list belong to this molecule_kind
624  nfixed_atoms = 0
625  DO j = 1, nmolecule
626  molecule => molecule_set(molecule_list(j))
627  CALL get_molecule(molecule, first_atom=first, last_atom=last)
628  fix_atom_molname = .false.
629  IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
630  DO k = 1, SIZE(cons_info%fixed_molnames)
631  IF (cons_info%fixed_molnames(k) .EQ. molname) THEN
632  fix_atom_molname = .true.
633  IF (is_qm .AND. cons_info%fixed_exclude_qm(k)) fix_atom_molname = .false.
634  IF ((.NOT. is_qm) .AND. cons_info%fixed_exclude_mm(k)) fix_atom_molname = .false.
635  END IF
636  END DO
637  END IF
638  DO k = first, last
639  fix_atom_qmmm = .false.
640  IF (PRESENT(qmmm_env)) THEN
641  SELECT CASE (cons_info%freeze_qm)
642  CASE (do_constr_atomic)
643  IF (any(qmmm_env%qm_atom_index == k)) fix_atom_qmmm = .true.
644  CASE (do_constr_molec)
645  IF (any(qmmm_env%qm_molecule_index == molecule_list(j))) fix_atom_qmmm = .true.
646  END SELECT
647  SELECT CASE (cons_info%freeze_mm)
648  CASE (do_constr_atomic)
649  IF (all(qmmm_env%qm_atom_index /= k)) fix_atom_qmmm = .true.
650  CASE (do_constr_molec)
651  IF (all(qmmm_env%qm_molecule_index /= molecule_list(j))) fix_atom_qmmm = .true.
652  END SELECT
653  END IF
654  IF (any(cons_info%fixed_atoms == k) .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
655  nfixed_atoms = nfixed_atoms + 1
656  END IF
657  END DO
658  END DO
659  ALLOCATE (fixd_list(nfixed_atoms))
660  kk = 0
661  nfixd_restraint = 0
662  IF (nfixed_atoms /= 0) THEN
663  DO j = 1, nmolecule
664  molecule => molecule_set(molecule_list(j))
665  CALL get_molecule(molecule, first_atom=first, last_atom=last)
666  fix_atom_molname = .false.
667  IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
668  DO k1loc = 1, SIZE(cons_info%fixed_molnames)
669  IF (cons_info%fixed_molnames(k1loc) .EQ. molname) THEN
670  fix_atom_molname = .true.
671  itype = cons_info%fixed_mol_type(k1loc)
672  EXIT
673  END IF
674  END DO
675  END IF
676  DO k = first, last
677  ! FIXED LIST ATOMS
678  fix_fixed_atom = .false.
679  DO k2loc = 1, SIZE(cons_info%fixed_atoms)
680  IF (cons_info%fixed_atoms(k2loc) == k) THEN
681  fix_fixed_atom = .true.
682  itype = cons_info%fixed_type(k2loc)
683  EXIT
684  END IF
685  END DO
686  ! QMMM FIXED ATOMS (QM OR MM)
687  fix_atom_qmmm = .false.
688  fix_atom_mm = .false.
689  fix_atom_qm = .false.
690  IF (PRESENT(qmmm_env)) THEN
691  SELECT CASE (cons_info%freeze_qm)
692  CASE (do_constr_atomic)
693  IF (any(qmmm_env%qm_atom_index == k)) THEN
694  fix_atom_qmmm = .true.
695  fix_atom_qm = .true.
696  itype = cons_info%freeze_qm_type
697  END IF
698  CASE (do_constr_molec)
699  IF (any(qmmm_env%qm_molecule_index == molecule_list(j))) THEN
700  fix_atom_qmmm = .true.
701  fix_atom_qm = .true.
702  itype = cons_info%freeze_qm_type
703  END IF
704  END SELECT
705  SELECT CASE (cons_info%freeze_mm)
706  CASE (do_constr_atomic)
707  IF (all(qmmm_env%qm_atom_index /= k)) THEN
708  fix_atom_qmmm = .true.
709  fix_atom_mm = .true.
710  itype = cons_info%freeze_mm_type
711  END IF
712  CASE (do_constr_molec)
713  IF (all(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN
714  fix_atom_qmmm = .true.
715  fix_atom_mm = .true.
716  itype = cons_info%freeze_mm_type
717  END IF
718  END SELECT
719  ! We should never reach this point but let's check it anyway
720  IF (fix_atom_qm .AND. fix_atom_mm) THEN
721  CALL cp_abort(__location__, &
722  "Atom number: "//cp_to_string(k)// &
723  " has been defined both QM and MM. General Error!")
724  END IF
725  END IF
726  ! Check that the fixed atom constraint/restraint is unique
727  IF ((fix_fixed_atom .AND. fix_atom_qmmm) .OR. (fix_fixed_atom .AND. fix_atom_molname) &
728  .OR. (fix_atom_qmmm .AND. fix_atom_molname)) THEN
729  CALL cp_abort(__location__, &
730  "Atom number: "//cp_to_string(k)// &
731  " has been constrained/restrained to be fixed in more than one"// &
732  " input section. Check and correct your input file!")
733  END IF
734  ! Let's store the atom index
735  IF (fix_fixed_atom .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
736  kk = kk + 1
737  fixd_list(kk)%fixd = k
738  fixd_list(kk)%coord = particle_set(k)%r
739  fixd_list(kk)%itype = itype
740  ! Possibly Restraint
741  IF (fix_fixed_atom) THEN
742  fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc)
743  fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc)
744  ELSEIF (fix_atom_qm) THEN
745  fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint
746  fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0
747  ELSEIF (fix_atom_mm) THEN
748  fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint
749  fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0
750  ELSEIF (fix_atom_molname) THEN
751  fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc)
752  fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc)
753  ELSE
754  ! Should never reach this point
755  cpabort("")
756  END IF
757  IF (fixd_list(kk)%restraint%active) THEN
758  nfixd_restraint = nfixd_restraint + 1
759  nfixd_restart = nfixd_restart + 1
760  ! Check that we use the components that we really want..
761  SELECT CASE (itype)
762  CASE (use_perd_x)
763  fixd_list(kk)%coord(2) = huge(0.0_dp)
764  fixd_list(kk)%coord(3) = huge(0.0_dp)
765  CASE (use_perd_y)
766  fixd_list(kk)%coord(1) = huge(0.0_dp)
767  fixd_list(kk)%coord(3) = huge(0.0_dp)
768  CASE (use_perd_z)
769  fixd_list(kk)%coord(1) = huge(0.0_dp)
770  fixd_list(kk)%coord(2) = huge(0.0_dp)
771  CASE (use_perd_xy)
772  fixd_list(kk)%coord(3) = huge(0.0_dp)
773  CASE (use_perd_xz)
774  fixd_list(kk)%coord(2) = huge(0.0_dp)
775  CASE (use_perd_yz)
776  fixd_list(kk)%coord(1) = huge(0.0_dp)
777  END SELECT
778  IF (restart_restraint_pos) THEN
779  ! Read coord0 value for restraint
780  CALL section_vals_val_get(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
781  i_rep_val=nfixd_restart, r_vals=r)
782  SELECT CASE (itype)
783  CASE (use_perd_x)
784  cpassert(SIZE(r) == 1)
785  fixd_list(kk)%coord(1) = r(1)
786  CASE (use_perd_y)
787  cpassert(SIZE(r) == 1)
788  fixd_list(kk)%coord(2) = r(1)
789  CASE (use_perd_z)
790  cpassert(SIZE(r) == 1)
791  fixd_list(kk)%coord(3) = r(1)
792  CASE (use_perd_xy)
793  cpassert(SIZE(r) == 2)
794  fixd_list(kk)%coord(1) = r(1)
795  fixd_list(kk)%coord(2) = r(2)
796  CASE (use_perd_xz)
797  cpassert(SIZE(r) == 2)
798  fixd_list(kk)%coord(1) = r(1)
799  fixd_list(kk)%coord(3) = r(2)
800  CASE (use_perd_yz)
801  cpassert(SIZE(r) == 2)
802  fixd_list(kk)%coord(2) = r(1)
803  fixd_list(kk)%coord(3) = r(2)
804  CASE (use_perd_xyz)
805  cpassert(SIZE(r) == 3)
806  fixd_list(kk)%coord(1:3) = r(1:3)
807  END SELECT
808  ELSE
809  ! Write coord0 value for restraint
810  SELECT CASE (itype)
811  CASE (use_perd_x)
812  ALLOCATE (r(1))
813  r(1) = fixd_list(kk)%coord(1)
814  CASE (use_perd_y)
815  ALLOCATE (r(1))
816  r(1) = fixd_list(kk)%coord(2)
817  CASE (use_perd_z)
818  ALLOCATE (r(1))
819  r(1) = fixd_list(kk)%coord(3)
820  CASE (use_perd_xy)
821  ALLOCATE (r(2))
822  r(1) = fixd_list(kk)%coord(1)
823  r(2) = fixd_list(kk)%coord(2)
824  CASE (use_perd_xz)
825  ALLOCATE (r(2))
826  r(1) = fixd_list(kk)%coord(1)
827  r(2) = fixd_list(kk)%coord(3)
828  CASE (use_perd_yz)
829  ALLOCATE (r(2))
830  r(1) = fixd_list(kk)%coord(1)
831  r(2) = fixd_list(kk)%coord(3)
832  CASE (use_perd_xyz)
833  ALLOCATE (r(3))
834  r(1:3) = fixd_list(kk)%coord(1:3)
835  END SELECT
836  CALL section_vals_val_set(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
837  i_rep_val=nfixd_restart, r_vals_ptr=r)
838  END IF
839  END IF
840  END IF
841  END DO
842  END DO
843  END IF
844  IF (iw > 0) THEN
845  WRITE (iw, *) "MOLECULE KIND:", i, " NR. FIXED ATOMS:", SIZE(fixd_list(:)%fixd), " LIST::", fixd_list(:)%fixd
846  END IF
847  CALL set_molecule_kind(molecule_kind, nfixd=nfixed_atoms, nfixd_restraint=nfixd_restraint, &
848  fixd_list=fixd_list)
849  fixd_list_gci(nfixd_list_gci + 1:nfixd_list_gci + nfixed_atoms) = fixd_list
850  nfixd_list_gci = nfixd_list_gci + nfixed_atoms
851  END DO
852  IF (iw > 0) THEN
853  WRITE (iw, *) "TOTAL NUMBER OF FIXED ATOMS:", nfixd_list_gci
854  END IF
855  cpassert(count(missed_molname) == 0)
856  DEALLOCATE (missed_molname)
857  ! Intermolecular constraints
858  IF (gci%ntot /= 0) THEN
859  ALLOCATE (fixd_list(nfixd_list_gci))
860  fixd_list(1:nfixd_list_gci) = fixd_list_gci(1:nfixd_list_gci)
861  gci%fixd_list => fixd_list
862  END IF
863  DEALLOCATE (fixd_list_gci)
864  END IF
865  ! Final setup of the number of possible restraints
866  gci%nrestraint = gci%ng3x3_restraint + &
867  gci%ng4x6_restraint + &
868  gci%nvsite_restraint + &
869  gci%ncolv%nrestraint
870  CALL cp_print_key_finished_output(iw, logger, subsys_section, &
871  "PRINT%TOPOLOGY_INFO/UTIL_INFO")
872  CALL timestop(handle2)
873  CALL timestop(handle)
874  END SUBROUTINE topology_constraint_pack
875 
876 ! **************************************************************************************************
877 !> \brief Setup the colv_list for the packing of constraints
878 !> \param colv_list ...
879 !> \param ilist ...
880 !> \param gind ...
881 !> \param cons_info ...
882 !> \param topology ...
883 !> \param particle_set ...
884 !> \param restart_restraint_clv ...
885 !> \param colvar_rest ...
886 !> \param first_atom ...
887 !> \par History
888 !> Updated 2007 for intermolecular constraints
889 !> \author Teodoro Laino [2007]
890 ! **************************************************************************************************
891  SUBROUTINE setup_colv_list(colv_list, ilist, gind, cons_info, topology, &
892  particle_set, restart_restraint_clv, colvar_rest, first_atom)
893 
894  TYPE(colvar_constraint_type), DIMENSION(:), &
895  POINTER :: colv_list
896  INTEGER, DIMENSION(:), POINTER :: ilist
897  INTEGER, INTENT(INOUT) :: gind
898  TYPE(constraint_info_type), POINTER :: cons_info
899  TYPE(topology_parameters_type), INTENT(INOUT) :: topology
900  TYPE(particle_type), DIMENSION(:), OPTIONAL, &
901  POINTER :: particle_set
902  LOGICAL, INTENT(IN) :: restart_restraint_clv
903  TYPE(section_vals_type), POINTER :: colvar_rest
904  INTEGER, INTENT(IN) :: first_atom
905 
906  INTEGER :: j, kdim, kk, ncolv_mol
907  REAL(kind=dp) :: rmod
908  TYPE(colvar_type), POINTER :: local_colvar
909 
910  ncolv_mol = 0
911  DO kk = 1, SIZE(ilist)
912  j = ilist(kk)
913  ncolv_mol = ncolv_mol + 1
914  kdim = SIZE(cons_info%colvar_set(j)%colvar%i_atom)
915  ALLOCATE (colv_list(ncolv_mol)%i_atoms(kdim))
916  colv_list(ncolv_mol)%inp_seq_num = j
917  colv_list(ncolv_mol)%type_id = cons_info%colvar_set(j)%colvar%type_id
918  colv_list(ncolv_mol)%i_atoms = cons_info%colvar_set(j)%colvar%i_atom
919  colv_list(ncolv_mol)%use_points = cons_info%colvar_set(j)%colvar%use_points
920  ! Restraint
921  colv_list(ncolv_mol)%restraint%active = cons_info%colv_restraint(j)
922  colv_list(ncolv_mol)%restraint%k0 = cons_info%colv_k0(j)
923  IF (cons_info%const_colv_target(j) == -huge(0.0_dp)) THEN
924  ! Let's compute the value..
925  NULLIFY (local_colvar)
926  CALL colvar_clone(local_colvar, cons_info%colvar_set(j)%colvar, &
927  i_atom_offset=first_atom - 1)
928  CALL colvar_eval_mol_f(local_colvar, topology%cell, particle_set)
929  colv_list(ncolv_mol)%expected_value = local_colvar%ss
930  CALL colvar_release(local_colvar)
931  ELSE
932  colv_list(ncolv_mol)%expected_value = cons_info%const_colv_target(j)
933  END IF
934  colv_list(ncolv_mol)%expected_value_growth_speed = cons_info%const_colv_target_growth(j)
935  ! In case of Restraint let's check for possible restart values
936  IF (colv_list(ncolv_mol)%restraint%active .AND. &
937  (colv_list(ncolv_mol)%expected_value_growth_speed == 0.0_dp)) THEN
938  gind = gind + 1
939  IF (restart_restraint_clv) THEN
940  CALL section_vals_val_get(colvar_rest, "_DEFAULT_KEYWORD_", &
941  i_rep_val=gind, r_val=rmod)
942  colv_list(ncolv_mol)%expected_value = rmod
943  ELSE
944  rmod = colv_list(ncolv_mol)%expected_value
945  CALL section_vals_val_set(colvar_rest, "_DEFAULT_KEYWORD_", &
946  i_rep_val=gind, r_val=rmod)
947  END IF
948  END IF
949  ! Only if torsion let's take into account the singularity in the definition
950  ! of the dihedral
951  IF (cons_info%colvar_set(j)%colvar%type_id == torsion_colvar_id) THEN
952  cons_info%colvar_set(j)%colvar%torsion_param%o0 = colv_list(ncolv_mol)%expected_value
953  END IF
954  END DO
955  END SUBROUTINE setup_colv_list
956 
957 ! **************************************************************************************************
958 !> \brief Setup the g3x3_list for the packing of constraints
959 !> \param g3x3_list ...
960 !> \param ilist ...
961 !> \param cons_info ...
962 !> \param ng3x3_restraint ...
963 !> \par History
964 !> Updated 2007 for intermolecular constraints
965 !> \author Teodoro Laino [2007]
966 ! **************************************************************************************************
967  SUBROUTINE setup_g3x3_list(g3x3_list, ilist, cons_info, ng3x3_restraint)
968  TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
969  INTEGER, DIMENSION(:), POINTER :: ilist
970  TYPE(constraint_info_type), POINTER :: cons_info
971  INTEGER, INTENT(OUT) :: ng3x3_restraint
972 
973  INTEGER :: j, ng3x3
974 
975  ng3x3_restraint = 0
976  DO ng3x3 = 1, SIZE(ilist)
977  j = ilist(ng3x3)
978  g3x3_list(ng3x3)%a = cons_info%const_g33_a(j)
979  g3x3_list(ng3x3)%b = cons_info%const_g33_b(j)
980  g3x3_list(ng3x3)%c = cons_info%const_g33_c(j)
981  g3x3_list(ng3x3)%dab = cons_info%const_g33_dab(j)
982  g3x3_list(ng3x3)%dac = cons_info%const_g33_dac(j)
983  g3x3_list(ng3x3)%dbc = cons_info%const_g33_dbc(j)
984  ! Restraint
985  g3x3_list(ng3x3)%restraint%active = cons_info%g33_restraint(j)
986  g3x3_list(ng3x3)%restraint%k0 = cons_info%g33_k0(j)
987  IF (g3x3_list(ng3x3)%restraint%active) ng3x3_restraint = ng3x3_restraint + 1
988  END DO
989 
990  END SUBROUTINE setup_g3x3_list
991 
992 ! **************************************************************************************************
993 !> \brief Setup the g4x6_list for the packing of constraints
994 !> \param g4x6_list ...
995 !> \param ilist ...
996 !> \param cons_info ...
997 !> \param ng4x6_restraint ...
998 !> \par History
999 !> Updated 2007 for intermolecular constraints
1000 !> \author Teodoro Laino [2007]
1001 ! **************************************************************************************************
1002  SUBROUTINE setup_g4x6_list(g4x6_list, ilist, cons_info, ng4x6_restraint)
1003  TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
1004  INTEGER, DIMENSION(:), POINTER :: ilist
1005  TYPE(constraint_info_type), POINTER :: cons_info
1006  INTEGER, INTENT(OUT) :: ng4x6_restraint
1007 
1008  INTEGER :: j, ng4x6
1009 
1010  ng4x6 = 0
1011  ng4x6_restraint = 0
1012  DO ng4x6 = 1, SIZE(ilist)
1013  j = ilist(ng4x6)
1014  g4x6_list(ng4x6)%a = cons_info%const_g46_a(j)
1015  g4x6_list(ng4x6)%b = cons_info%const_g46_b(j)
1016  g4x6_list(ng4x6)%c = cons_info%const_g46_c(j)
1017  g4x6_list(ng4x6)%d = cons_info%const_g46_d(j)
1018  g4x6_list(ng4x6)%dab = cons_info%const_g46_dab(j)
1019  g4x6_list(ng4x6)%dac = cons_info%const_g46_dac(j)
1020  g4x6_list(ng4x6)%dbc = cons_info%const_g46_dbc(j)
1021  g4x6_list(ng4x6)%dad = cons_info%const_g46_dad(j)
1022  g4x6_list(ng4x6)%dbd = cons_info%const_g46_dbd(j)
1023  g4x6_list(ng4x6)%dcd = cons_info%const_g46_dcd(j)
1024  ! Restraint
1025  g4x6_list(ng4x6)%restraint%active = cons_info%g46_restraint(j)
1026  g4x6_list(ng4x6)%restraint%k0 = cons_info%g46_k0(j)
1027  IF (g4x6_list(ng4x6)%restraint%active) ng4x6_restraint = ng4x6_restraint + 1
1028  END DO
1029 
1030  END SUBROUTINE setup_g4x6_list
1031 
1032 ! **************************************************************************************************
1033 !> \brief Setup the vsite_list for the packing of constraints
1034 !> \param vsite_list ...
1035 !> \param ilist ...
1036 !> \param cons_info ...
1037 !> \param nvsite_restraint ...
1038 !> \par History
1039 !> \author Marcel Baer [2008]
1040 ! **************************************************************************************************
1041  SUBROUTINE setup_vsite_list(vsite_list, ilist, cons_info, nvsite_restraint)
1042  TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
1043  INTEGER, DIMENSION(:), POINTER :: ilist
1044  TYPE(constraint_info_type), POINTER :: cons_info
1045  INTEGER, INTENT(OUT) :: nvsite_restraint
1046 
1047  INTEGER :: j, nvsite
1048 
1049  nvsite = 0
1050  nvsite_restraint = 0
1051  DO nvsite = 1, SIZE(ilist)
1052  j = ilist(nvsite)
1053  vsite_list(nvsite)%a = cons_info%const_vsite_a(j)
1054  vsite_list(nvsite)%b = cons_info%const_vsite_b(j)
1055  vsite_list(nvsite)%c = cons_info%const_vsite_c(j)
1056  vsite_list(nvsite)%d = cons_info%const_vsite_d(j)
1057  vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j)
1058  vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j)
1059  ! Restraint
1060  vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j)
1061  vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j)
1062  IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint + 1
1063  END DO
1064 
1065  END SUBROUTINE setup_vsite_list
1066 ! **************************************************************************************************
1067 !> \brief Setup the lcolv for the packing of constraints
1068 !> \param lcolv ...
1069 !> \param ilist ...
1070 !> \param first_atom ...
1071 !> \param last_atom ...
1072 !> \param cons_info ...
1073 !> \param particle_set ...
1074 !> \param colvar_func_info ...
1075 !> \param use_clv_info ...
1076 !> \param cind ...
1077 !> \par History
1078 !> Updated 2007 for intermolecular constraints
1079 !> \author Teodoro Laino [2007]
1080 ! **************************************************************************************************
1081  SUBROUTINE setup_lcolv(lcolv, ilist, first_atom, last_atom, cons_info, &
1082  particle_set, colvar_func_info, use_clv_info, &
1083  cind)
1084  TYPE(local_colvar_constraint_type), DIMENSION(:), &
1085  POINTER :: lcolv
1086  INTEGER, DIMENSION(:), POINTER :: ilist
1087  INTEGER, INTENT(IN) :: first_atom, last_atom
1088  TYPE(constraint_info_type), POINTER :: cons_info
1089  TYPE(particle_type), DIMENSION(:), OPTIONAL, &
1090  POINTER :: particle_set
1091  TYPE(section_vals_type), POINTER :: colvar_func_info
1092  LOGICAL, INTENT(IN) :: use_clv_info
1093  INTEGER, INTENT(INOUT) :: cind
1094 
1095  INTEGER :: ind, k, kk
1096  REAL(kind=dp), DIMENSION(:), POINTER :: r_vals
1097 
1098  DO kk = 1, SIZE(ilist)
1099  k = ilist(kk)
1100  lcolv(kk)%init = .false.
1101  lcolv(kk)%lambda = 0.0_dp
1102  lcolv(kk)%sigma = 0.0_dp
1103 
1104  ! Set Up colvar variable
1105  NULLIFY (lcolv(kk)%colvar, lcolv(kk)%colvar_old)
1106  ! Colvar
1107  CALL colvar_clone(lcolv(kk)%colvar, cons_info%colvar_set(k)%colvar, &
1108  i_atom_offset=first_atom - 1)
1109 
1110  ! Some COLVARS may need additional information for evaluating the
1111  ! functional form: this is the case for COLVARS which depend on the
1112  ! initial position of the atoms: This information is stored in a proper
1113  ! container in the COLVAR_RESTART section..
1114  IF ((lcolv(kk)%colvar%type_id == xyz_diag_colvar_id) .OR. &
1115  (lcolv(kk)%colvar%type_id == xyz_outerdiag_colvar_id)) THEN
1116  cind = cind + 1
1117  IF (use_clv_info) THEN
1118  CALL section_vals_val_get(colvar_func_info, "_DEFAULT_KEYWORD_", &
1119  i_rep_val=cind, r_vals=r_vals)
1120  SELECT CASE (lcolv(kk)%colvar%type_id)
1121  CASE (xyz_diag_colvar_id)
1122  cpassert(SIZE(r_vals) == 3)
1123  lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1124  CASE (xyz_outerdiag_colvar_id)
1125  cpassert(SIZE(r_vals) == 6)
1126  lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1127  lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6)
1128  END SELECT
1129  ELSE
1130  SELECT CASE (lcolv(kk)%colvar%type_id)
1131  CASE (xyz_diag_colvar_id)
1132  ALLOCATE (r_vals(3))
1133  ind = first_atom - 1 + lcolv(kk)%colvar%xyz_diag_param%i_atom
1134  r_vals = particle_set(ind)%r
1135  lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1136  CASE (xyz_outerdiag_colvar_id)
1137  ALLOCATE (r_vals(6))
1138  ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(1)
1139  r_vals(1:3) = particle_set(ind)%r
1140  ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(2)
1141  r_vals(4:6) = particle_set(ind)%r
1142  lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1143  lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6)
1144  END SELECT
1145  CALL section_vals_val_set(colvar_func_info, "_DEFAULT_KEYWORD_", &
1146  i_rep_val=cind, r_vals_ptr=r_vals)
1147  END IF
1148  END IF
1149 
1150  ! Setup Colvar_old
1151  CALL colvar_clone(lcolv(kk)%colvar_old, lcolv(kk)%colvar)
1152 
1153  ! Check for consistency in the constraint definition
1154  IF (any(lcolv(kk)%colvar%i_atom > last_atom) .OR. &
1155  any(lcolv(kk)%colvar%i_atom < first_atom)) THEN
1156  WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1157  WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", &
1158  " but the atoms specified in the constraint and the atoms defined for", &
1159  " the molecule DO NOT match!", &
1160  "This could be very probable due to a wrong connectivity, or an error", &
1161  " in the constraint specification in the input file.", &
1162  " Please check it carefully!"
1163  cpabort("")
1164  END IF
1165  END DO
1166  END SUBROUTINE setup_lcolv
1167 
1168 ! **************************************************************************************************
1169 !> \brief Setup the lg3x3 for the packing of constraints
1170 !> \param lg3x3 ...
1171 !> \param g3x3_list ...
1172 !> \param first_atom ...
1173 !> \param last_atom ...
1174 !> \par History
1175 !> Updated 2007 for intermolecular constraints
1176 !> \author Teodoro Laino [2007]
1177 ! **************************************************************************************************
1178  SUBROUTINE setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
1179  TYPE(local_g3x3_constraint_type), DIMENSION(:), &
1180  POINTER :: lg3x3
1181  TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
1182  INTEGER, INTENT(IN) :: first_atom, last_atom
1183 
1184  INTEGER :: kk
1185 
1186  DO kk = 1, SIZE(lg3x3)
1187  lg3x3(kk)%init = .false.
1188  lg3x3(kk)%scale = 0.0_dp
1189  lg3x3(kk)%scale_old = 0.0_dp
1190  lg3x3(kk)%fa = 0.0_dp
1191  lg3x3(kk)%fb = 0.0_dp
1192  lg3x3(kk)%fc = 0.0_dp
1193  lg3x3(kk)%ra_old = 0.0_dp
1194  lg3x3(kk)%rb_old = 0.0_dp
1195  lg3x3(kk)%rc_old = 0.0_dp
1196  lg3x3(kk)%va = 0.0_dp
1197  lg3x3(kk)%vb = 0.0_dp
1198  lg3x3(kk)%vc = 0.0_dp
1199  lg3x3(kk)%lambda = 0.0_dp
1200  IF ((g3x3_list(kk)%a + first_atom - 1 < first_atom) .OR. &
1201  (g3x3_list(kk)%b + first_atom - 1 < first_atom) .OR. &
1202  (g3x3_list(kk)%c + first_atom - 1 < first_atom) .OR. &
1203  (g3x3_list(kk)%a + first_atom - 1 > last_atom) .OR. &
1204  (g3x3_list(kk)%b + first_atom - 1 > last_atom) .OR. &
1205  (g3x3_list(kk)%c + first_atom - 1 > last_atom)) THEN
1206  WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1207  WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", &
1208  " but the atoms specified in the constraint and the atoms defined for", &
1209  " the molecule DO NOT match!", &
1210  "This could be very probable due to a wrong connectivity, or an error", &
1211  " in the constraint specification in the input file.", &
1212  " Please check it carefully!"
1213  cpabort("")
1214  END IF
1215  END DO
1216 
1217  END SUBROUTINE setup_lg3x3
1218 
1219 ! **************************************************************************************************
1220 !> \brief Setup the lg4x6 for the packing of constraints
1221 !> \param lg4x6 ...
1222 !> \param g4x6_list ...
1223 !> \param first_atom ...
1224 !> \param last_atom ...
1225 !> \par History
1226 !> Updated 2007 for intermolecular constraints
1227 !> \author Teodoro Laino [2007]
1228 ! **************************************************************************************************
1229  SUBROUTINE setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
1230  TYPE(local_g4x6_constraint_type), DIMENSION(:), &
1231  POINTER :: lg4x6
1232  TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
1233  INTEGER, INTENT(IN) :: first_atom, last_atom
1234 
1235  INTEGER :: kk
1236 
1237  DO kk = 1, SIZE(lg4x6)
1238  lg4x6(kk)%init = .false.
1239  lg4x6(kk)%scale = 0.0_dp
1240  lg4x6(kk)%scale_old = 0.0_dp
1241  lg4x6(kk)%fa = 0.0_dp
1242  lg4x6(kk)%fb = 0.0_dp
1243  lg4x6(kk)%fc = 0.0_dp
1244  lg4x6(kk)%fd = 0.0_dp
1245  lg4x6(kk)%fe = 0.0_dp
1246  lg4x6(kk)%ff = 0.0_dp
1247  lg4x6(kk)%ra_old = 0.0_dp
1248  lg4x6(kk)%rb_old = 0.0_dp
1249  lg4x6(kk)%rc_old = 0.0_dp
1250  lg4x6(kk)%rd_old = 0.0_dp
1251  lg4x6(kk)%re_old = 0.0_dp
1252  lg4x6(kk)%rf_old = 0.0_dp
1253  lg4x6(kk)%va = 0.0_dp
1254  lg4x6(kk)%vb = 0.0_dp
1255  lg4x6(kk)%vc = 0.0_dp
1256  lg4x6(kk)%vd = 0.0_dp
1257  lg4x6(kk)%ve = 0.0_dp
1258  lg4x6(kk)%vf = 0.0_dp
1259  lg4x6(kk)%lambda = 0.0_dp
1260  IF ((g4x6_list(kk)%a + first_atom - 1 < first_atom) .OR. &
1261  (g4x6_list(kk)%b + first_atom - 1 < first_atom) .OR. &
1262  (g4x6_list(kk)%c + first_atom - 1 < first_atom) .OR. &
1263  (g4x6_list(kk)%d + first_atom - 1 < first_atom) .OR. &
1264  (g4x6_list(kk)%a + first_atom - 1 > last_atom) .OR. &
1265  (g4x6_list(kk)%b + first_atom - 1 > last_atom) .OR. &
1266  (g4x6_list(kk)%c + first_atom - 1 > last_atom) .OR. &
1267  (g4x6_list(kk)%d + first_atom - 1 > last_atom)) THEN
1268  WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1269  WRITE (*, '(T5,"|",T8,A)') "A constrained has been defined for a molecule type", &
1270  " but the atoms specified in the constraint and the atoms defined for", &
1271  " the molecule DO NOT match!", &
1272  "This could be very probable due to a wrong connectivity, or an error", &
1273  " in the constraint specification in the input file.", &
1274  " Please check it carefully!"
1275  cpabort("")
1276  END IF
1277  END DO
1278 
1279  END SUBROUTINE setup_lg4x6
1280 
1281 ! **************************************************************************************************
1282 !> \brief Gives back a list of molecule to which apply the constraint
1283 !> \param const_mol ...
1284 !> \param const_molname ...
1285 !> \param const_intermolecular ...
1286 !> \param constr_x_mol ...
1287 !> \param constr_x_glob ...
1288 !> \param molecule_kind_set ...
1289 !> \param exclude_qm ...
1290 !> \param exclude_mm ...
1291 !> \par History
1292 !> Updated 2007 for intermolecular constraints
1293 !> \author Teodoro Laino [2006]
1294 ! **************************************************************************************************
1295  SUBROUTINE give_constraint_array(const_mol, const_molname, const_intermolecular, &
1296  constr_x_mol, constr_x_glob, molecule_kind_set, exclude_qm, exclude_mm)
1297 
1298  INTEGER, DIMENSION(:), POINTER :: const_mol
1299  CHARACTER(LEN=default_string_length), &
1300  DIMENSION(:), POINTER :: const_molname
1301  LOGICAL, DIMENSION(:), POINTER :: const_intermolecular
1302  TYPE(constr_list_type), DIMENSION(:), POINTER :: constr_x_mol
1303  INTEGER, DIMENSION(:), POINTER :: constr_x_glob
1304  TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
1305  LOGICAL, DIMENSION(:), POINTER :: exclude_qm, exclude_mm
1306 
1307  CHARACTER(len=*), PARAMETER :: routinen = 'give_constraint_array'
1308 
1309  CHARACTER(LEN=default_string_length) :: myname, name
1310  INTEGER :: handle, i, iglob, isize, k
1311  LOGICAL :: found_molname, is_qm
1312  TYPE(molecule_kind_type), POINTER :: molecule_kind
1313 
1314  CALL timeset(routinen, handle)
1315  NULLIFY (molecule_kind)
1316  ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
1317  DO i = 1, SIZE(constr_x_mol)
1318  NULLIFY (constr_x_mol(i)%constr)
1319  ALLOCATE (constr_x_mol(i)%constr(0))
1320  END DO
1321  cpassert(SIZE(const_mol) == SIZE(const_molname))
1322  iglob = 0
1323  DO i = 1, SIZE(const_mol)
1324  IF (const_intermolecular(i)) THEN
1325  ! Intermolecular constraint
1326  iglob = iglob + 1
1327  CALL reallocate(constr_x_glob, 1, iglob)
1328  constr_x_glob(iglob) = i
1329  ELSE
1330  ! Intramolecular constraint
1331  IF (const_mol(i) /= 0) THEN
1332  k = const_mol(i)
1333  IF (k > SIZE(molecule_kind_set)) &
1334  CALL cp_abort(__location__, &
1335  "A constraint has been specified providing the molecule index. But the"// &
1336  " molecule index ("//cp_to_string(k)//") is out of range of the possible"// &
1337  " molecule kinds ("//cp_to_string(SIZE(molecule_kind_set))//").")
1338  isize = SIZE(constr_x_mol(k)%constr)
1339  CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1340  constr_x_mol(k)%constr(isize + 1) = i
1341  ELSE
1342  myname = const_molname(i)
1343  found_molname = .false.
1344  DO k = 1, SIZE(molecule_kind_set)
1345  molecule_kind => molecule_kind_set(k)
1346  name = molecule_kind%name
1347  is_qm = qmmm_ff_precond_only_qm(id1=name)
1348  IF (is_qm .AND. exclude_qm(i)) cycle
1349  IF (.NOT. is_qm .AND. exclude_mm(i)) cycle
1350  IF (name == myname) THEN
1351  isize = SIZE(constr_x_mol(k)%constr)
1352  CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1353  constr_x_mol(k)%constr(isize + 1) = i
1354  found_molname = .true.
1355  END IF
1356  END DO
1357  CALL print_warning_molname(found_molname, myname)
1358  END IF
1359  END IF
1360  END DO
1361  CALL timestop(handle)
1362  END SUBROUTINE give_constraint_array
1363 
1364 ! **************************************************************************************************
1365 !> \brief Prints a warning message if undefined molnames are used to define constraints
1366 !> \param found ...
1367 !> \param name ...
1368 !> \author Teodoro Laino [2007] - Zurich University
1369 ! **************************************************************************************************
1370  SUBROUTINE print_warning_molname(found, name)
1371  LOGICAL, INTENT(IN) :: found
1372  CHARACTER(LEN=*), INTENT(IN) :: name
1373 
1374  IF (.NOT. found) &
1375  CALL cp_warn(__location__, &
1376  " MOLNAME ("//trim(name)//") was defined for constraints, but this molecule name "// &
1377  "is not defined. Please check carefully your PDB, PSF (has priority over PDB) or "// &
1378  "input driven CP2K coordinates. In case you may not find the reason for this warning "// &
1379  "it may be a good idea to print all molecule information (including kind name) activating "// &
1380  "the print_key MOLECULES specific of the SUBSYS%PRINT section. ")
1381 
1382  END SUBROUTINE print_warning_molname
1383 
1384 END MODULE topology_constraint_util
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
atomic_kind_types::is_hydrogen
elemental logical function, public is_hydrogen(atomic_kind)
Determines if the atomic_kind is HYDROGEN.
Definition: atomic_kind_types.F:380
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition: cell_types.F:42
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition: cell_types.F:42
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition: cell_types.F:42
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition: cell_types.F:42
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition: cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition: cell_types.F:42
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition: cell_types.F:42
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition: colvar_methods.F:16
colvar_methods::colvar_eval_mol_f
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Definition: colvar_methods.F:1631
colvar_types
Initialize the collective variables types.
Definition: colvar_types.F:15
colvar_types::colvar_release
recursive subroutine, public colvar_release(colvar)
releases the memory that might have been allocated by the colvar
Definition: colvar_types.F:1176
colvar_types::xyz_outerdiag_colvar_id
integer, parameter, public xyz_outerdiag_colvar_id
Definition: colvar_types.F:56
colvar_types::colvar_create
subroutine, public colvar_create(colvar, colvar_id)
initializes a colvar_param type
Definition: colvar_types.F:421
colvar_types::colvar_clone
recursive subroutine, public colvar_clone(colvar_out, colvar_in, i_atom_offset)
Clone a colvar type.
Definition: colvar_types.F:1353
colvar_types::colvar_p_reallocate
subroutine, public colvar_p_reallocate(colvar_set, lb1_new, ub1_new)
Change the dimension of a colvar_p_type.
Definition: colvar_types.F:1791
colvar_types::dist_colvar_id
integer, parameter, public dist_colvar_id
Definition: colvar_types.F:56
colvar_types::colvar_setup
subroutine, public colvar_setup(colvar)
Finalize the setup of the collective variable.
Definition: colvar_types.F:502
colvar_types::torsion_colvar_id
integer, parameter, public torsion_colvar_id
Definition: colvar_types.F:56
colvar_types::xyz_diag_colvar_id
integer, parameter, public xyz_diag_colvar_id
Definition: colvar_types.F:56
colvar_utils
evaluations of colvar for internal coordinates schemes
Definition: colvar_utils.F:14
colvar_utils::post_process_colvar
subroutine, public post_process_colvar(colvar, particles)
Complete the description of the COORDINATION colvar when defined using KINDS.
Definition: colvar_utils.F:566
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_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
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_constr_atomic
integer, parameter, public do_constr_atomic
Definition: input_constants.F:79
input_constants::do_constr_molec
integer, parameter, public do_constr_molec
Definition: input_constants.F:79
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_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition: input_section_types.F:1216
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
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::dp
integer, parameter, public dp
Definition: kinds.F:34
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::setup_colvar_counters
subroutine, public setup_colvar_counters(colv_list, ncolv)
...
Definition: molecule_kind_types.F:239
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_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::set_molecule
subroutine, public set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
Set a molecule data set.
Definition: molecule_types.F:337
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qmmm_ff_fist
Definition: qmmm_ff_fist.F:11
qmmm_ff_fist::qmmm_ff_precond_only_qm
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
Definition: qmmm_ff_fist.F:39
qmmm_types_low
Definition: qmmm_types_low.F:13
topology_constraint_util
Collection of subroutine needed for topology related things.
Definition: topology_constraint_util.F:13
topology_constraint_util::topology_constraint_pack
subroutine, public topology_constraint_pack(molecule_kind_set, molecule_set, topology, qmmm_env, particle_set, input_file, subsys_section, gci)
Pack in all the information needed for the constraints.
Definition: topology_constraint_util.F:87
topology_types
Definition: topology_types.F:18
topology
Control for reading in different topologies and coordinates.
Definition: topology.F:13
topology_types::constr_list_type
Definition: topology_types.F:179
topology_types::constraint_info_type
Definition: topology_types.F:64
topology_types::topology_parameters_type
Definition: topology_types.F:142