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