Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Collection of subroutine needed for topology related things
10 : !> \par History
11 : !> jgh (23-05-2004) Last atom of molecule information added
12 : ! **************************************************************************************************
13 : 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 9736 : 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 9736 : 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 9736 : 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 9736 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: missed_molname
110 : REAL(KIND=dp) :: rmod, rvec(3)
111 9736 : REAL(KIND=dp), DIMENSION(:), POINTER :: hdist, r
112 9736 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
113 : TYPE(atomic_kind_type), POINTER :: atomic_kind
114 9736 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
115 : TYPE(colvar_constraint_type), DIMENSION(:), &
116 9736 : POINTER :: colv_list
117 : TYPE(colvar_counters) :: ncolv
118 9736 : 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 9736 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list, fixd_list_gci
122 9736 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
123 9736 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
124 : TYPE(local_colvar_constraint_type), DIMENSION(:), &
125 9736 : POINTER :: lcolv
126 : TYPE(local_constraint_type), POINTER :: lci
127 : TYPE(local_g3x3_constraint_type), DIMENSION(:), &
128 9736 : POINTER :: lg3x3
129 : TYPE(local_g4x6_constraint_type), DIMENSION(:), &
130 9736 : 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 9736 : TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
136 :
137 9736 : NULLIFY (logger, constr_x_mol, constr_x_glob)
138 19472 : logger => cp_get_default_logger()
139 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
140 9736 : extension=".subsysLog")
141 9736 : CALL timeset(routineN, handle)
142 9736 : CALL timeset(routineN//"_1", handle2)
143 :
144 9736 : cons_info => topology%cons_info
145 : hbonds_section => section_vals_get_subs_vals(input_file, &
146 9736 : "MOTION%CONSTRAINT%HBONDS")
147 : fixd_restr_rest => section_vals_get_subs_vals(input_file, &
148 9736 : "MOTION%CONSTRAINT%FIX_ATOM_RESTART")
149 9736 : CALL section_vals_get(fixd_restr_rest, explicit=restart_restraint_pos)
150 : colvar_rest => section_vals_get_subs_vals(input_file, &
151 9736 : "MOTION%CONSTRAINT%COLVAR_RESTART")
152 9736 : CALL section_vals_get(colvar_rest, explicit=restart_restraint_clv)
153 : colvar_func_info => section_vals_get_subs_vals(subsys_section, &
154 9736 : "COLVAR%COLVAR_FUNC_INFO")
155 9736 : 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 311296 : DO i = 1, topology%nmol
161 301560 : molecule => molecule_set(i)
162 301560 : 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 301560 : topology%const_colv .OR. topology%const_vsite) THEN
167 43692 : ALLOCATE (lci)
168 : NULLIFY (lci%lcolv)
169 : NULLIFY (lci%lg3x3)
170 : NULLIFY (lci%lg4x6)
171 : END IF
172 311296 : CALL set_molecule(molecule, lci=lci)
173 : END DO
174 9736 : 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 9736 : CALL setup_colvar_counters(gci%colv_list, gci%ncolv)
191 : gci%nrestraint = gci%ng3x3_restraint + &
192 : gci%ng4x6_restraint + &
193 : gci%nvsite_restraint + &
194 9736 : gci%ncolv%nrestraint
195 9736 : CALL timestop(handle2)
196 9736 : CALL timeset(routineN//"_2", handle2)
197 : !-----------------------------------------------------------------------------
198 : !-----------------------------------------------------------------------------
199 : ! 2. Add more stuff to COLVAR constraint if constraint hydrogen is on
200 : !-----------------------------------------------------------------------------
201 9736 : IF (topology%const_hydr) THEN
202 16 : topology%const_colv = .TRUE.
203 16 : NULLIFY (atom_typeh, hdist)
204 98 : ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
205 66 : DO i = 1, SIZE(molecule_kind_set)
206 50 : ALLOCATE (constr_x_mol(i)%constr(1))
207 66 : constr_x_mol(i)%constr(1) = 1
208 : END DO
209 16 : CALL section_vals_val_get(hbonds_section, "MOLECULE", n_rep_val=nrep)
210 16 : IF (nrep /= 0) THEN
211 4 : NULLIFY (inds)
212 36 : DO i = 1, SIZE(molecule_kind_set)
213 36 : constr_x_mol(i)%constr(1) = 0
214 : END DO
215 4 : CALL section_vals_val_get(hbonds_section, "MOLECULE", i_vals=inds)
216 32 : DO i = 1, SIZE(inds)
217 32 : constr_x_mol(inds(i))%constr(1) = 1
218 : END DO
219 : ELSE
220 12 : CALL section_vals_val_get(hbonds_section, "MOLNAME", n_rep_val=nrep)
221 12 : IF (nrep /= 0) THEN
222 2 : NULLIFY (cnds)
223 10 : DO i = 1, SIZE(molecule_kind_set)
224 10 : constr_x_mol(i)%constr(1) = 0
225 : END DO
226 2 : CALL section_vals_val_get(hbonds_section, "MOLNAME", c_vals=cnds)
227 4 : DO i = 1, SIZE(cnds)
228 2 : found_molname = .FALSE.
229 10 : DO k = 1, SIZE(molecule_kind_set)
230 8 : molecule_kind => molecule_kind_set(k)
231 8 : name = molecule_kind%name
232 8 : ldummy = qmmm_ff_precond_only_qm(id1=name)
233 10 : IF (cnds(i) == name) THEN
234 4 : constr_x_mol(k)%constr(1) = 1
235 4 : found_molname = .TRUE.
236 : END IF
237 : END DO
238 4 : CALL print_warning_molname(found_molname, cnds(i))
239 : END DO
240 : END IF
241 : END IF
242 16 : CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", n_rep_val=nrep)
243 16 : IF (nrep /= 0) &
244 8 : CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", c_vals=atom_typeh)
245 16 : CALL section_vals_val_get(hbonds_section, "TARGETS", n_rep_val=nrep)
246 16 : IF (nrep /= 0) &
247 8 : CALL section_vals_val_get(hbonds_section, "TARGETS", r_vals=hdist)
248 16 : IF (ASSOCIATED(hdist)) THEN
249 8 : CPASSERT(SIZE(hdist) == SIZE(atom_typeh))
250 : END IF
251 16 : CALL section_vals_val_get(hbonds_section, "exclude_qm", l_val=exclude_qm)
252 16 : CALL section_vals_val_get(hbonds_section, "exclude_mm", l_val=exclude_mm)
253 16 : nhdist = 0
254 66 : DO i = 1, SIZE(molecule_kind_set)
255 50 : molecule_kind => molecule_kind_set(i)
256 50 : 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 42 : molecule_list=molecule_list)
260 : ! Let's tag all requested atoms involving Hydrogen
261 : ! on the first molecule of this kind
262 42 : molecule => molecule_set(molecule_list(1))
263 42 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
264 42 : natom = last_atom - first_atom + 1
265 464 : DO k = 1, nbond
266 364 : ishbond = .FALSE.
267 364 : j = bond_list(k)%a
268 364 : IF (j < 1 .OR. j > natom) CYCLE
269 364 : atomic_kind => atom_list(j)%atomic_kind
270 364 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
271 364 : is_qm = qmmm_ff_precond_only_qm(id1=name)
272 364 : IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE.
273 364 : IF (is_qm .AND. exclude_qm) ishbond = .FALSE.
274 364 : IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE.
275 332 : IF (.NOT. ishbond) THEN
276 364 : j = bond_list(k)%b
277 364 : IF (j < 1 .OR. j > natom) CYCLE
278 344 : atomic_kind => atom_list(j)%atomic_kind
279 344 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
280 344 : is_qm = qmmm_ff_precond_only_qm(id1=name)
281 344 : IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE.
282 344 : IF (is_qm .AND. exclude_qm) ishbond = .FALSE.
283 344 : IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE.
284 : END IF
285 354 : IF (ishbond) THEN
286 180 : nhdist = nhdist + 1
287 : END IF
288 : END DO
289 : END DO
290 16 : n_start_colv = cons_info%nconst_colv
291 16 : cons_info%nconst_colv = nhdist + n_start_colv
292 16 : CALL reallocate(cons_info%const_colv_mol, 1, cons_info%nconst_colv)
293 16 : CALL reallocate(cons_info%const_colv_molname, 1, cons_info%nconst_colv)
294 16 : CALL reallocate(cons_info%const_colv_target, 1, cons_info%nconst_colv)
295 16 : CALL reallocate(cons_info%const_colv_target_growth, 1, cons_info%nconst_colv)
296 16 : CALL colvar_p_reallocate(cons_info%colvar_set, 1, cons_info%nconst_colv)
297 : ! Fill in Restraints info
298 16 : CALL reallocate(cons_info%colv_intermolecular, 1, cons_info%nconst_colv)
299 16 : CALL reallocate(cons_info%colv_restraint, 1, cons_info%nconst_colv)
300 16 : CALL reallocate(cons_info%colv_k0, 1, cons_info%nconst_colv)
301 16 : CALL reallocate(cons_info%colv_exclude_qm, 1, cons_info%nconst_colv)
302 16 : CALL reallocate(cons_info%colv_exclude_mm, 1, cons_info%nconst_colv)
303 : ! Bonds involving hydrogens are by their nature only intramolecular
304 196 : cons_info%colv_intermolecular(n_start_colv + 1:cons_info%nconst_colv) = .FALSE.
305 196 : cons_info%colv_exclude_qm(n_start_colv + 1:cons_info%nconst_colv) = .FALSE.
306 196 : cons_info%colv_exclude_mm(n_start_colv + 1:cons_info%nconst_colv) = .FALSE.
307 196 : cons_info%colv_restraint(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_restraint
308 196 : cons_info%colv_k0(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_k0
309 : !
310 16 : nhdist = 0
311 66 : DO i = 1, SIZE(molecule_kind_set)
312 50 : IF (constr_x_mol(i)%constr(1) == 0) CYCLE
313 42 : 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 42 : molecule_list=molecule_list)
317 42 : molecule => molecule_set(molecule_list(1))
318 42 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
319 42 : natom = last_atom - first_atom + 1
320 42 : offset = first_atom - 1
321 464 : DO k = 1, nbond
322 364 : ishbond = .FALSE.
323 364 : j = bond_list(k)%a
324 364 : IF (j < 1 .OR. j > natom) CYCLE
325 364 : atomic_kind => atom_list(j)%atomic_kind
326 364 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
327 364 : is_qm = qmmm_ff_precond_only_qm(id1=name)
328 364 : IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE.
329 364 : IF (is_qm .AND. exclude_qm) ishbond = .FALSE.
330 364 : IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE.
331 332 : IF (.NOT. ishbond) THEN
332 364 : j = bond_list(k)%b
333 364 : IF (j < 1 .OR. j > natom) CYCLE
334 344 : atomic_kind => atom_list(j)%atomic_kind
335 344 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
336 344 : is_qm = qmmm_ff_precond_only_qm(id1=name)
337 344 : IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE.
338 344 : IF (is_qm .AND. exclude_qm) ishbond = .FALSE.
339 344 : IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE.
340 : END IF
341 394 : IF (ishbond) THEN
342 180 : nhdist = nhdist + 1
343 720 : rvec = particle_set(offset + bond_list(k)%a)%r - particle_set(offset + bond_list(k)%b)%r
344 720 : rmod = SQRT(DOT_PRODUCT(rvec, rvec))
345 180 : IF (ASSOCIATED(hdist)) THEN
346 32 : IF (SIZE(hdist) > 0) THEN
347 32 : IF (bond_list(k)%a == j) atomic_kind => atom_list(bond_list(k)%b)%atomic_kind
348 32 : 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 32 : name=name, element_symbol=element_symbol)
351 32 : ldummy = qmmm_ff_precond_only_qm(id1=name)
352 32 : DO m = 1, SIZE(hdist)
353 32 : IF (TRIM(name) == TRIM(atom_typeh(m))) EXIT
354 32 : IF (TRIM(element_symbol) == TRIM(atom_typeh(m))) EXIT
355 : END DO
356 32 : IF (m <= SIZE(hdist)) THEN
357 32 : rmod = hdist(m)
358 : END IF
359 : END IF
360 : END IF
361 180 : cons_info%const_colv_mol(nhdist + n_start_colv) = i
362 180 : cons_info%const_colv_molname(nhdist + n_start_colv) = "UNDEF"
363 180 : cons_info%const_colv_target(nhdist + n_start_colv) = rmod
364 180 : 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 180 : dist_colvar_id)
367 180 : cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%i_at = bond_list(k)%a
368 180 : cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%j_at = bond_list(k)%b
369 180 : CALL colvar_setup(cons_info%colvar_set(nhdist + n_start_colv)%colvar)
370 : END IF
371 : END DO
372 : END DO
373 66 : DO j = 1, SIZE(constr_x_mol)
374 66 : DEALLOCATE (constr_x_mol(j)%constr)
375 : END DO
376 80 : DEALLOCATE (constr_x_mol)
377 : END IF
378 :
379 9736 : CALL timestop(handle2)
380 9736 : CALL timeset(routineN//"_3", handle2)
381 : !-----------------------------------------------------------------------------
382 : !-----------------------------------------------------------------------------
383 : ! 3. Set the COLVAR constraint molecule_kind_set(ikind)%colv_list
384 : !-----------------------------------------------------------------------------
385 9736 : IF (topology%const_colv) THEN
386 : ! Post Process of COLVARS..
387 586 : DO ii = 1, SIZE(cons_info%colvar_set)
388 586 : 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 136 : cons_info%colv_exclude_mm)
399 : ! Intramolecular constraints
400 136 : gind = 0
401 136 : cind = 0
402 714 : DO ii = 1, SIZE(molecule_kind_set)
403 578 : molecule_kind => molecule_kind_set(ii)
404 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
405 578 : nmolecule=nmolecule, molecule_list=molecule_list)
406 578 : ncolv_mol = SIZE(constr_x_mol(ii)%constr)
407 1660 : 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 578 : molecule => molecule_set(molecule_list(1))
413 578 : 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 578 : colvar_rest, first_atom)
417 578 : CALL setup_colvar_counters(colv_list, ncolv)
418 578 : CALL set_molecule_kind(molecule_kind, colv_list=colv_list, ncolv=ncolv)
419 3978 : DO j = 1, nmolecule
420 2108 : molecule => molecule_set(molecule_list(j))
421 2108 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
422 7014 : ALLOCATE (lcolv(ncolv_mol))
423 : CALL setup_lcolv(lcolv, constr_x_mol(ii)%constr, first_atom, last_atom, &
424 2108 : cons_info, particle_set, colvar_func_info, use_clv_info, cind)
425 2686 : CALL set_molecule(molecule=molecule, lcolv=lcolv)
426 : END DO
427 : END DO
428 714 : DO j = 1, SIZE(constr_x_mol)
429 714 : DEALLOCATE (constr_x_mol(j)%constr)
430 : END DO
431 136 : DEALLOCATE (constr_x_mol)
432 : ! Intermolecular constraints
433 136 : ncolv_glob = 0
434 136 : IF (ASSOCIATED(constr_x_glob)) THEN
435 44 : ncolv_glob = SIZE(constr_x_glob)
436 198 : 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 44 : first_atom=1)
440 44 : CALL setup_colvar_counters(colv_list, ncolv)
441 198 : ALLOCATE (lcolv(ncolv_glob))
442 : CALL setup_lcolv(lcolv, constr_x_glob, 1, SIZE(particle_set), cons_info, &
443 44 : particle_set, colvar_func_info, use_clv_info, cind)
444 44 : gci%colv_list => colv_list
445 44 : gci%lcolv => lcolv
446 44 : gci%ncolv = ncolv
447 : ! Total number of Intermolecular constraints
448 44 : gci%ntot = gci%ncolv%ntot + gci%ntot
449 88 : DEALLOCATE (constr_x_glob)
450 : END IF
451 : END IF
452 :
453 9736 : CALL timestop(handle2)
454 9736 : CALL timeset(routineN//"_4", handle2)
455 : !-----------------------------------------------------------------------------
456 : !-----------------------------------------------------------------------------
457 : ! 4. Set the group 3x3 constraint g3x3_list
458 : !-----------------------------------------------------------------------------
459 9736 : 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 156 : cons_info%g33_exclude_mm)
468 : ! Intramolecular constraints
469 426 : DO ii = 1, SIZE(molecule_kind_set)
470 270 : molecule_kind => molecule_kind_set(ii)
471 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
472 : nmolecule=nmolecule, &
473 270 : molecule_list=molecule_list)
474 270 : ng3x3 = SIZE(constr_x_mol(ii)%constr)
475 852 : ALLOCATE (g3x3_list(ng3x3))
476 270 : CALL setup_g3x3_list(g3x3_list, constr_x_mol(ii)%constr, cons_info, ng3x3_restraint)
477 270 : CALL set_molecule_kind(molecule_kind, ng3x3=ng3x3, ng3x3_restraint=ng3x3_restraint, g3x3_list=g3x3_list)
478 37320 : DO j = 1, nmolecule
479 36354 : molecule => molecule_set(molecule_list(j))
480 36354 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
481 2524548 : ALLOCATE (lg3x3(ng3x3))
482 36354 : CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
483 36624 : CALL set_molecule(molecule=molecule, lg3x3=lg3x3)
484 : END DO
485 : END DO
486 426 : DO j = 1, SIZE(constr_x_mol)
487 426 : DEALLOCATE (constr_x_mol(j)%constr)
488 : END DO
489 156 : DEALLOCATE (constr_x_mol)
490 : ! Intermolecular constraints
491 156 : IF (ASSOCIATED(constr_x_glob)) THEN
492 4 : ng3x3 = SIZE(constr_x_glob)
493 16 : ALLOCATE (g3x3_list(ng3x3))
494 4 : CALL setup_g3x3_list(g3x3_list, constr_x_glob, cons_info, ng3x3_restraint)
495 280 : ALLOCATE (lg3x3(ng3x3))
496 4 : CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
497 4 : gci%g3x3_list => g3x3_list
498 4 : gci%lg3x3 => lg3x3
499 4 : gci%ng3x3 = ng3x3
500 4 : gci%ng3x3_restraint = ng3x3_restraint
501 : ! Total number of Intermolecular constraints
502 4 : gci%ntot = 3*gci%ng3x3 + gci%ntot
503 8 : DEALLOCATE (constr_x_glob)
504 : END IF
505 : END IF
506 :
507 9736 : CALL timestop(handle2)
508 9736 : CALL timeset(routineN//"_5", handle2)
509 : !-----------------------------------------------------------------------------
510 : !-----------------------------------------------------------------------------
511 : ! 5. Set the group 4x6 constraint g4x6_list
512 : !-----------------------------------------------------------------------------
513 9736 : 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 16 : cons_info%g46_exclude_mm)
522 : ! Intramolecular constraints
523 36 : DO ii = 1, SIZE(molecule_kind_set)
524 20 : molecule_kind => molecule_kind_set(ii)
525 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
526 20 : nmolecule=nmolecule, molecule_list=molecule_list)
527 20 : ng4x6 = SIZE(constr_x_mol(ii)%constr)
528 64 : ALLOCATE (g4x6_list(ng4x6))
529 20 : CALL setup_g4x6_list(g4x6_list, constr_x_mol(ii)%constr, cons_info, ng4x6_restraint)
530 20 : CALL set_molecule_kind(molecule_kind, ng4x6=ng4x6, ng4x6_restraint=ng4x6_restraint, g4x6_list=g4x6_list)
531 726 : DO j = 1, nmolecule
532 650 : molecule => molecule_set(molecule_list(j))
533 650 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
534 99580 : ALLOCATE (lg4x6(ng4x6))
535 650 : CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
536 670 : CALL set_molecule(molecule=molecule, lg4x6=lg4x6)
537 : END DO
538 : END DO
539 36 : DO j = 1, SIZE(constr_x_mol)
540 36 : DEALLOCATE (constr_x_mol(j)%constr)
541 : END DO
542 16 : DEALLOCATE (constr_x_mol)
543 : ! Intermolecular constraints
544 16 : IF (ASSOCIATED(constr_x_glob)) THEN
545 4 : ng4x6 = SIZE(constr_x_glob)
546 16 : ALLOCATE (g4x6_list(ng4x6))
547 4 : CALL setup_g4x6_list(g4x6_list, constr_x_glob, cons_info, ng4x6_restraint)
548 616 : ALLOCATE (lg4x6(ng4x6))
549 4 : CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
550 4 : gci%g4x6_list => g4x6_list
551 4 : gci%lg4x6 => lg4x6
552 4 : gci%ng4x6 = ng4x6
553 4 : gci%ng4x6_restraint = ng4x6_restraint
554 : ! Total number of Intermolecular constraints
555 4 : gci%ntot = 6*gci%ng4x6 + gci%ntot
556 8 : DEALLOCATE (constr_x_glob)
557 : END IF
558 : END IF
559 :
560 9736 : CALL timestop(handle2)
561 9736 : CALL timeset(routineN//"_6", handle2)
562 : !-----------------------------------------------------------------------------
563 : !-----------------------------------------------------------------------------
564 : ! 6. Set the group vsite constraint vsite_list
565 : !-----------------------------------------------------------------------------
566 9736 : 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 8 : cons_info%vsite_exclude_mm)
575 : ! Intramolecular constraints
576 18 : DO ii = 1, SIZE(molecule_kind_set)
577 10 : molecule_kind => molecule_kind_set(ii)
578 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
579 10 : nmolecule=nmolecule, molecule_list=molecule_list)
580 10 : nvsite = SIZE(constr_x_mol(ii)%constr)
581 36 : ALLOCATE (vsite_list(nvsite))
582 10 : 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 28 : vsite_list=vsite_list)
585 : END DO
586 18 : DO j = 1, SIZE(constr_x_mol)
587 18 : DEALLOCATE (constr_x_mol(j)%constr)
588 : END DO
589 8 : DEALLOCATE (constr_x_mol)
590 : ! Intermolecular constraints
591 8 : IF (ASSOCIATED(constr_x_glob)) THEN
592 0 : nvsite = SIZE(constr_x_glob)
593 0 : ALLOCATE (vsite_list(nvsite))
594 0 : CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint)
595 0 : gci%vsite_list => vsite_list
596 0 : gci%nvsite = nvsite
597 0 : gci%nvsite_restraint = nvsite_restraint
598 : ! Total number of Intermolecular constraints
599 0 : gci%ntot = gci%nvsite + gci%ntot
600 0 : DEALLOCATE (constr_x_glob)
601 : END IF
602 : END IF
603 9736 : CALL timestop(handle2)
604 9736 : CALL timeset(routineN//"_7", handle2)
605 : !-----------------------------------------------------------------------------
606 : !-----------------------------------------------------------------------------
607 : ! 7. Set the group fixed_atom constraint fixd_list
608 : !-----------------------------------------------------------------------------
609 9736 : IF (topology%const_atom) THEN
610 30574 : ALLOCATE (fixd_list_gci(SIZE(particle_set)))
611 110 : nfixd_list_gci = 0
612 226 : ALLOCATE (missed_molname(SIZE(cons_info%fixed_molnames, 1)))
613 116 : missed_molname = .TRUE.
614 110 : nfixd_restart = 0
615 5036 : DO i = 1, SIZE(molecule_kind_set)
616 4926 : molecule_kind => molecule_kind_set(i)
617 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
618 4926 : nmolecule=nmolecule, molecule_list=molecule_list, name=molname)
619 4926 : is_qm = qmmm_ff_precond_only_qm(id1=molname)
620 4938 : WHERE (molname == 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 4926 : nfixed_atoms = 0
625 17634 : DO j = 1, nmolecule
626 12708 : molecule => molecule_set(molecule_list(j))
627 12708 : CALL get_molecule(molecule, first_atom=first, last_atom=last)
628 12708 : fix_atom_molname = .FALSE.
629 12708 : IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
630 14274 : DO k = 1, SIZE(cons_info%fixed_molnames)
631 14274 : IF (cons_info%fixed_molnames(k) == molname) THEN
632 48 : fix_atom_molname = .TRUE.
633 48 : IF (is_qm .AND. cons_info%fixed_exclude_qm(k)) fix_atom_molname = .FALSE.
634 48 : 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 47548 : DO k = first, last
639 29914 : fix_atom_qmmm = .FALSE.
640 29914 : IF (PRESENT(qmmm_env)) THEN
641 324 : SELECT CASE (cons_info%freeze_qm)
642 : CASE (do_constr_atomic)
643 0 : IF (ANY(qmmm_env%qm_atom_index == k)) fix_atom_qmmm = .TRUE.
644 : CASE (do_constr_molec)
645 336 : IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) fix_atom_qmmm = .TRUE.
646 : END SELECT
647 394 : SELECT CASE (cons_info%freeze_mm)
648 : CASE (do_constr_atomic)
649 840 : IF (ALL(qmmm_env%qm_atom_index /= k)) fix_atom_qmmm = .TRUE.
650 : CASE (do_constr_molec)
651 408 : IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) fix_atom_qmmm = .TRUE.
652 : END SELECT
653 : END IF
654 3861838 : IF (ANY(cons_info%fixed_atoms == k) .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
655 10196 : nfixed_atoms = nfixed_atoms + 1
656 : END IF
657 : END DO
658 : END DO
659 39006 : ALLOCATE (fixd_list(nfixed_atoms))
660 4926 : kk = 0
661 4926 : nfixd_restraint = 0
662 4926 : IF (nfixed_atoms /= 0) THEN
663 10122 : DO j = 1, nmolecule
664 5942 : molecule => molecule_set(molecule_list(j))
665 5942 : CALL get_molecule(molecule, first_atom=first, last_atom=last)
666 5942 : fix_atom_molname = .FALSE.
667 5942 : IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
668 5942 : DO k1loc = 1, SIZE(cons_info%fixed_molnames)
669 5942 : IF (cons_info%fixed_molnames(k1loc) == molname) THEN
670 44 : fix_atom_molname = .TRUE.
671 44 : itype = cons_info%fixed_mol_type(k1loc)
672 44 : EXIT
673 : END IF
674 : END DO
675 : END IF
676 21162 : DO k = first, last
677 : ! FIXED LIST ATOMS
678 11040 : fix_fixed_atom = .FALSE.
679 2891634 : DO k2loc = 1, SIZE(cons_info%fixed_atoms)
680 2891634 : IF (cons_info%fixed_atoms(k2loc) == k) THEN
681 10012 : fix_fixed_atom = .TRUE.
682 10012 : itype = cons_info%fixed_type(k2loc)
683 10012 : EXIT
684 : END IF
685 : END DO
686 : ! QMMM FIXED ATOMS (QM OR MM)
687 11040 : fix_atom_qmmm = .FALSE.
688 11040 : fix_atom_mm = .FALSE.
689 11040 : fix_atom_qm = .FALSE.
690 11040 : IF (PRESENT(qmmm_env)) THEN
691 224 : SELECT CASE (cons_info%freeze_qm)
692 : CASE (do_constr_atomic)
693 0 : IF (ANY(qmmm_env%qm_atom_index == k)) THEN
694 0 : fix_atom_qmmm = .TRUE.
695 0 : fix_atom_qm = .TRUE.
696 0 : itype = cons_info%freeze_qm_type
697 : END IF
698 : CASE (do_constr_molec)
699 224 : IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) THEN
700 6 : fix_atom_qmmm = .TRUE.
701 6 : fix_atom_qm = .TRUE.
702 6 : itype = cons_info%freeze_qm_type
703 : END IF
704 : END SELECT
705 294 : SELECT CASE (cons_info%freeze_mm)
706 : CASE (do_constr_atomic)
707 840 : IF (ALL(qmmm_env%qm_atom_index /= k)) THEN
708 42 : fix_atom_qmmm = .TRUE.
709 42 : fix_atom_mm = .TRUE.
710 42 : itype = cons_info%freeze_mm_type
711 : END IF
712 : CASE (do_constr_molec)
713 308 : IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN
714 84 : fix_atom_qmmm = .TRUE.
715 84 : fix_atom_mm = .TRUE.
716 84 : 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 224 : IF (fix_atom_qm .AND. fix_atom_mm) THEN
721 : CALL cp_abort(__LOCATION__, &
722 : "Atom number: "//cp_to_string(k)// &
723 0 : " 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 11040 : .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 0 : " input section. Check and correct your input file!")
733 : END IF
734 : ! Let's store the atom index
735 16982 : IF (fix_fixed_atom .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
736 10196 : kk = kk + 1
737 10196 : fixd_list(kk)%fixd = k
738 71372 : fixd_list(kk)%coord = particle_set(k)%r
739 10196 : fixd_list(kk)%itype = itype
740 : ! Possibly Restraint
741 10196 : IF (fix_fixed_atom) THEN
742 10012 : fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc)
743 10012 : fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc)
744 184 : ELSEIF (fix_atom_qm) THEN
745 6 : fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint
746 6 : fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0
747 178 : ELSEIF (fix_atom_mm) THEN
748 126 : fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint
749 126 : fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0
750 52 : ELSEIF (fix_atom_molname) THEN
751 52 : fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc)
752 52 : fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc)
753 : ELSE
754 : ! Should never reach this point
755 0 : CPABORT("")
756 : END IF
757 10196 : IF (fixd_list(kk)%restraint%active) THEN
758 38 : nfixd_restraint = nfixd_restraint + 1
759 38 : nfixd_restart = nfixd_restart + 1
760 : ! Check that we use the components that we really want..
761 0 : SELECT CASE (itype)
762 : CASE (use_perd_x)
763 0 : fixd_list(kk)%coord(2) = HUGE(0.0_dp)
764 0 : fixd_list(kk)%coord(3) = HUGE(0.0_dp)
765 : CASE (use_perd_y)
766 0 : fixd_list(kk)%coord(1) = HUGE(0.0_dp)
767 0 : fixd_list(kk)%coord(3) = HUGE(0.0_dp)
768 : CASE (use_perd_z)
769 0 : fixd_list(kk)%coord(1) = HUGE(0.0_dp)
770 0 : fixd_list(kk)%coord(2) = HUGE(0.0_dp)
771 : CASE (use_perd_xy)
772 0 : fixd_list(kk)%coord(3) = HUGE(0.0_dp)
773 : CASE (use_perd_xz)
774 0 : fixd_list(kk)%coord(2) = HUGE(0.0_dp)
775 : CASE (use_perd_yz)
776 38 : fixd_list(kk)%coord(1) = HUGE(0.0_dp)
777 : END SELECT
778 38 : IF (restart_restraint_pos) THEN
779 : ! Read coord0 value for restraint
780 : CALL section_vals_val_get(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
781 14 : i_rep_val=nfixd_restart, r_vals=r)
782 0 : SELECT CASE (itype)
783 : CASE (use_perd_x)
784 0 : CPASSERT(SIZE(r) == 1)
785 0 : fixd_list(kk)%coord(1) = r(1)
786 : CASE (use_perd_y)
787 0 : CPASSERT(SIZE(r) == 1)
788 0 : fixd_list(kk)%coord(2) = r(1)
789 : CASE (use_perd_z)
790 0 : CPASSERT(SIZE(r) == 1)
791 0 : fixd_list(kk)%coord(3) = r(1)
792 : CASE (use_perd_xy)
793 0 : CPASSERT(SIZE(r) == 2)
794 0 : fixd_list(kk)%coord(1) = r(1)
795 0 : fixd_list(kk)%coord(2) = r(2)
796 : CASE (use_perd_xz)
797 0 : CPASSERT(SIZE(r) == 2)
798 0 : fixd_list(kk)%coord(1) = r(1)
799 0 : fixd_list(kk)%coord(3) = r(2)
800 : CASE (use_perd_yz)
801 0 : CPASSERT(SIZE(r) == 2)
802 0 : fixd_list(kk)%coord(2) = r(1)
803 0 : fixd_list(kk)%coord(3) = r(2)
804 : CASE (use_perd_xyz)
805 14 : CPASSERT(SIZE(r) == 3)
806 112 : fixd_list(kk)%coord(1:3) = r(1:3)
807 : END SELECT
808 : ELSE
809 : ! Write coord0 value for restraint
810 0 : SELECT CASE (itype)
811 : CASE (use_perd_x)
812 0 : ALLOCATE (r(1))
813 0 : r(1) = fixd_list(kk)%coord(1)
814 : CASE (use_perd_y)
815 0 : ALLOCATE (r(1))
816 0 : r(1) = fixd_list(kk)%coord(2)
817 : CASE (use_perd_z)
818 0 : ALLOCATE (r(1))
819 0 : r(1) = fixd_list(kk)%coord(3)
820 : CASE (use_perd_xy)
821 0 : ALLOCATE (r(2))
822 0 : r(1) = fixd_list(kk)%coord(1)
823 0 : r(2) = fixd_list(kk)%coord(2)
824 : CASE (use_perd_xz)
825 0 : ALLOCATE (r(2))
826 0 : r(1) = fixd_list(kk)%coord(1)
827 0 : r(2) = fixd_list(kk)%coord(3)
828 : CASE (use_perd_yz)
829 0 : ALLOCATE (r(2))
830 0 : r(1) = fixd_list(kk)%coord(1)
831 0 : r(2) = fixd_list(kk)%coord(3)
832 : CASE (use_perd_xyz)
833 24 : ALLOCATE (r(3))
834 120 : r(1:3) = fixd_list(kk)%coord(1:3)
835 : END SELECT
836 : CALL section_vals_val_set(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
837 24 : 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 4926 : IF (iw > 0) THEN
845 0 : 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 4926 : fixd_list=fixd_list)
849 25318 : fixd_list_gci(nfixd_list_gci + 1:nfixd_list_gci + nfixed_atoms) = fixd_list
850 9962 : nfixd_list_gci = nfixd_list_gci + nfixed_atoms
851 : END DO
852 110 : IF (iw > 0) THEN
853 0 : WRITE (iw, *) "TOTAL NUMBER OF FIXED ATOMS:", nfixd_list_gci
854 : END IF
855 116 : CPASSERT(COUNT(missed_molname) == 0)
856 110 : DEALLOCATE (missed_molname)
857 : ! Intermolecular constraints
858 110 : IF (gci%ntot /= 0) THEN
859 16 : ALLOCATE (fixd_list(nfixd_list_gci))
860 10 : fixd_list(1:nfixd_list_gci) = fixd_list_gci(1:nfixd_list_gci)
861 2 : gci%fixd_list => fixd_list
862 : END IF
863 110 : 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 9736 : gci%ncolv%nrestraint
870 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
871 9736 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
872 9736 : CALL timestop(handle2)
873 9736 : CALL timestop(handle)
874 9736 : 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 622 : 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 622 : ncolv_mol = 0
911 1070 : DO kk = 1, SIZE(ilist)
912 448 : j = ilist(kk)
913 448 : ncolv_mol = ncolv_mol + 1
914 448 : kdim = SIZE(cons_info%colvar_set(j)%colvar%i_atom)
915 1344 : ALLOCATE (colv_list(ncolv_mol)%i_atoms(kdim))
916 448 : colv_list(ncolv_mol)%inp_seq_num = j
917 448 : colv_list(ncolv_mol)%type_id = cons_info%colvar_set(j)%colvar%type_id
918 3020 : colv_list(ncolv_mol)%i_atoms = cons_info%colvar_set(j)%colvar%i_atom
919 448 : colv_list(ncolv_mol)%use_points = cons_info%colvar_set(j)%colvar%use_points
920 : ! Restraint
921 448 : colv_list(ncolv_mol)%restraint%active = cons_info%colv_restraint(j)
922 448 : colv_list(ncolv_mol)%restraint%k0 = cons_info%colv_k0(j)
923 448 : IF (cons_info%const_colv_target(j) == -HUGE(0.0_dp)) THEN
924 : ! Let's compute the value..
925 100 : NULLIFY (local_colvar)
926 : CALL colvar_clone(local_colvar, cons_info%colvar_set(j)%colvar, &
927 100 : i_atom_offset=first_atom - 1)
928 100 : CALL colvar_eval_mol_f(local_colvar, topology%cell, particle_set)
929 100 : colv_list(ncolv_mol)%expected_value = local_colvar%ss
930 100 : CALL colvar_release(local_colvar)
931 : ELSE
932 348 : colv_list(ncolv_mol)%expected_value = cons_info%const_colv_target(j)
933 : END IF
934 448 : 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 448 : IF (colv_list(ncolv_mol)%restraint%active .AND. &
937 : (colv_list(ncolv_mol)%expected_value_growth_speed == 0.0_dp)) THEN
938 96 : gind = gind + 1
939 96 : IF (restart_restraint_clv) THEN
940 : CALL section_vals_val_get(colvar_rest, "_DEFAULT_KEYWORD_", &
941 14 : i_rep_val=gind, r_val=rmod)
942 14 : colv_list(ncolv_mol)%expected_value = rmod
943 : ELSE
944 82 : rmod = colv_list(ncolv_mol)%expected_value
945 : CALL section_vals_val_set(colvar_rest, "_DEFAULT_KEYWORD_", &
946 82 : 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 1070 : IF (cons_info%colvar_set(j)%colvar%type_id == torsion_colvar_id) THEN
952 38 : cons_info%colvar_set(j)%colvar%torsion_param%o0 = colv_list(ncolv_mol)%expected_value
953 : END IF
954 : END DO
955 622 : 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 274 : 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 274 : ng3x3_restraint = 0
976 434 : DO ng3x3 = 1, SIZE(ilist)
977 160 : j = ilist(ng3x3)
978 160 : g3x3_list(ng3x3)%a = cons_info%const_g33_a(j)
979 160 : g3x3_list(ng3x3)%b = cons_info%const_g33_b(j)
980 160 : g3x3_list(ng3x3)%c = cons_info%const_g33_c(j)
981 160 : g3x3_list(ng3x3)%dab = cons_info%const_g33_dab(j)
982 160 : g3x3_list(ng3x3)%dac = cons_info%const_g33_dac(j)
983 160 : g3x3_list(ng3x3)%dbc = cons_info%const_g33_dbc(j)
984 : ! Restraint
985 160 : g3x3_list(ng3x3)%restraint%active = cons_info%g33_restraint(j)
986 160 : g3x3_list(ng3x3)%restraint%k0 = cons_info%g33_k0(j)
987 434 : IF (g3x3_list(ng3x3)%restraint%active) ng3x3_restraint = ng3x3_restraint + 1
988 : END DO
989 :
990 274 : 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 24 : 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 24 : ng4x6 = 0
1011 24 : ng4x6_restraint = 0
1012 40 : DO ng4x6 = 1, SIZE(ilist)
1013 16 : j = ilist(ng4x6)
1014 16 : g4x6_list(ng4x6)%a = cons_info%const_g46_a(j)
1015 16 : g4x6_list(ng4x6)%b = cons_info%const_g46_b(j)
1016 16 : g4x6_list(ng4x6)%c = cons_info%const_g46_c(j)
1017 16 : g4x6_list(ng4x6)%d = cons_info%const_g46_d(j)
1018 16 : g4x6_list(ng4x6)%dab = cons_info%const_g46_dab(j)
1019 16 : g4x6_list(ng4x6)%dac = cons_info%const_g46_dac(j)
1020 16 : g4x6_list(ng4x6)%dbc = cons_info%const_g46_dbc(j)
1021 16 : g4x6_list(ng4x6)%dad = cons_info%const_g46_dad(j)
1022 16 : g4x6_list(ng4x6)%dbd = cons_info%const_g46_dbd(j)
1023 16 : g4x6_list(ng4x6)%dcd = cons_info%const_g46_dcd(j)
1024 : ! Restraint
1025 16 : g4x6_list(ng4x6)%restraint%active = cons_info%g46_restraint(j)
1026 16 : g4x6_list(ng4x6)%restraint%k0 = cons_info%g46_k0(j)
1027 40 : IF (g4x6_list(ng4x6)%restraint%active) ng4x6_restraint = ng4x6_restraint + 1
1028 : END DO
1029 :
1030 24 : 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 10 : 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 10 : nvsite = 0
1050 10 : nvsite_restraint = 0
1051 18 : DO nvsite = 1, SIZE(ilist)
1052 8 : j = ilist(nvsite)
1053 8 : vsite_list(nvsite)%a = cons_info%const_vsite_a(j)
1054 8 : vsite_list(nvsite)%b = cons_info%const_vsite_b(j)
1055 8 : vsite_list(nvsite)%c = cons_info%const_vsite_c(j)
1056 8 : vsite_list(nvsite)%d = cons_info%const_vsite_d(j)
1057 8 : vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j)
1058 8 : vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j)
1059 : ! Restraint
1060 8 : vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j)
1061 8 : vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j)
1062 18 : IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint + 1
1063 : END DO
1064 :
1065 10 : 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 2152 : 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 2152 : REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals
1097 :
1098 4446 : DO kk = 1, SIZE(ilist)
1099 2294 : k = ilist(kk)
1100 2294 : lcolv(kk)%init = .FALSE.
1101 2294 : lcolv(kk)%lambda = 0.0_dp
1102 2294 : lcolv(kk)%sigma = 0.0_dp
1103 :
1104 : ! Set Up colvar variable
1105 2294 : NULLIFY (lcolv(kk)%colvar, lcolv(kk)%colvar_old)
1106 : ! Colvar
1107 : CALL colvar_clone(lcolv(kk)%colvar, cons_info%colvar_set(k)%colvar, &
1108 2294 : 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 2294 : IF ((lcolv(kk)%colvar%type_id == xyz_diag_colvar_id) .OR. &
1115 : (lcolv(kk)%colvar%type_id == xyz_outerdiag_colvar_id)) THEN
1116 12 : cind = cind + 1
1117 12 : IF (use_clv_info) THEN
1118 : CALL section_vals_val_get(colvar_func_info, "_DEFAULT_KEYWORD_", &
1119 0 : i_rep_val=cind, r_vals=r_vals)
1120 0 : SELECT CASE (lcolv(kk)%colvar%type_id)
1121 : CASE (xyz_diag_colvar_id)
1122 0 : CPASSERT(SIZE(r_vals) == 3)
1123 0 : lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1124 : CASE (xyz_outerdiag_colvar_id)
1125 0 : CPASSERT(SIZE(r_vals) == 6)
1126 0 : lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1127 0 : lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6)
1128 : END SELECT
1129 : ELSE
1130 6 : SELECT CASE (lcolv(kk)%colvar%type_id)
1131 : CASE (xyz_diag_colvar_id)
1132 6 : ALLOCATE (r_vals(3))
1133 6 : ind = first_atom - 1 + lcolv(kk)%colvar%xyz_diag_param%i_atom
1134 24 : r_vals = particle_set(ind)%r
1135 48 : lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1136 : CASE (xyz_outerdiag_colvar_id)
1137 6 : ALLOCATE (r_vals(6))
1138 6 : ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(1)
1139 24 : r_vals(1:3) = particle_set(ind)%r
1140 6 : ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(2)
1141 24 : r_vals(4:6) = particle_set(ind)%r
1142 48 : lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1143 60 : 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 12 : i_rep_val=cind, r_vals_ptr=r_vals)
1147 : END IF
1148 : END IF
1149 :
1150 : ! Setup Colvar_old
1151 2294 : CALL colvar_clone(lcolv(kk)%colvar_old, lcolv(kk)%colvar)
1152 :
1153 : ! Check for consistency in the constraint definition
1154 14100 : IF (ANY(lcolv(kk)%colvar%i_atom > last_atom) .OR. &
1155 2152 : ANY(lcolv(kk)%colvar%i_atom < first_atom)) THEN
1156 0 : WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1157 0 : WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", &
1158 0 : " but the atoms specified in the constraint and the atoms defined for", &
1159 0 : " the molecule DO NOT match!", &
1160 0 : "This could be very probable due to a wrong connectivity, or an error", &
1161 0 : " in the constraint specification in the input file.", &
1162 0 : " Please check it carefully!"
1163 0 : CPABORT("")
1164 : END IF
1165 : END DO
1166 2152 : 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 36358 : 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 62600 : DO kk = 1, SIZE(lg3x3)
1187 26242 : lg3x3(kk)%init = .FALSE.
1188 26242 : lg3x3(kk)%scale = 0.0_dp
1189 26242 : lg3x3(kk)%scale_old = 0.0_dp
1190 104968 : lg3x3(kk)%fa = 0.0_dp
1191 104968 : lg3x3(kk)%fb = 0.0_dp
1192 104968 : lg3x3(kk)%fc = 0.0_dp
1193 104968 : lg3x3(kk)%ra_old = 0.0_dp
1194 104968 : lg3x3(kk)%rb_old = 0.0_dp
1195 104968 : lg3x3(kk)%rc_old = 0.0_dp
1196 104968 : lg3x3(kk)%va = 0.0_dp
1197 104968 : lg3x3(kk)%vb = 0.0_dp
1198 104968 : lg3x3(kk)%vc = 0.0_dp
1199 104968 : 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 26242 : (g3x3_list(kk)%b + first_atom - 1 > last_atom) .OR. &
1205 36358 : (g3x3_list(kk)%c + first_atom - 1 > last_atom)) THEN
1206 0 : WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1207 0 : WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", &
1208 0 : " but the atoms specified in the constraint and the atoms defined for", &
1209 0 : " the molecule DO NOT match!", &
1210 0 : "This could be very probable due to a wrong connectivity, or an error", &
1211 0 : " in the constraint specification in the input file.", &
1212 0 : " Please check it carefully!"
1213 0 : CPABORT("")
1214 : END IF
1215 : END DO
1216 :
1217 36358 : 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 654 : 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 1048 : DO kk = 1, SIZE(lg4x6)
1238 394 : lg4x6(kk)%init = .FALSE.
1239 394 : lg4x6(kk)%scale = 0.0_dp
1240 394 : lg4x6(kk)%scale_old = 0.0_dp
1241 1576 : lg4x6(kk)%fa = 0.0_dp
1242 1576 : lg4x6(kk)%fb = 0.0_dp
1243 1576 : lg4x6(kk)%fc = 0.0_dp
1244 1576 : lg4x6(kk)%fd = 0.0_dp
1245 1576 : lg4x6(kk)%fe = 0.0_dp
1246 1576 : lg4x6(kk)%ff = 0.0_dp
1247 1576 : lg4x6(kk)%ra_old = 0.0_dp
1248 1576 : lg4x6(kk)%rb_old = 0.0_dp
1249 1576 : lg4x6(kk)%rc_old = 0.0_dp
1250 1576 : lg4x6(kk)%rd_old = 0.0_dp
1251 1576 : lg4x6(kk)%re_old = 0.0_dp
1252 1576 : lg4x6(kk)%rf_old = 0.0_dp
1253 1576 : lg4x6(kk)%va = 0.0_dp
1254 1576 : lg4x6(kk)%vb = 0.0_dp
1255 1576 : lg4x6(kk)%vc = 0.0_dp
1256 1576 : lg4x6(kk)%vd = 0.0_dp
1257 1576 : lg4x6(kk)%ve = 0.0_dp
1258 1576 : lg4x6(kk)%vf = 0.0_dp
1259 2758 : 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 394 : (g4x6_list(kk)%c + first_atom - 1 > last_atom) .OR. &
1267 654 : (g4x6_list(kk)%d + first_atom - 1 > last_atom)) THEN
1268 0 : WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!"
1269 0 : WRITE (*, '(T5,"|",T8,A)') "A constrained has been defined for a molecule type", &
1270 0 : " but the atoms specified in the constraint and the atoms defined for", &
1271 0 : " the molecule DO NOT match!", &
1272 0 : "This could be very probable due to a wrong connectivity, or an error", &
1273 0 : " in the constraint specification in the input file.", &
1274 0 : " Please check it carefully!"
1275 0 : CPABORT("")
1276 : END IF
1277 : END DO
1278 :
1279 654 : 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 316 : 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 316 : CALL timeset(routineN, handle)
1315 316 : NULLIFY (molecule_kind)
1316 1826 : ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
1317 1194 : DO i = 1, SIZE(constr_x_mol)
1318 878 : NULLIFY (constr_x_mol(i)%constr)
1319 1194 : ALLOCATE (constr_x_mol(i)%constr(0))
1320 : END DO
1321 316 : CPASSERT(SIZE(const_mol) == SIZE(const_molname))
1322 316 : iglob = 0
1323 950 : DO i = 1, SIZE(const_mol)
1324 950 : IF (const_intermolecular(i)) THEN
1325 : ! Intermolecular constraint
1326 74 : iglob = iglob + 1
1327 74 : CALL reallocate(constr_x_glob, 1, iglob)
1328 74 : constr_x_glob(iglob) = i
1329 : ELSE
1330 : ! Intramolecular constraint
1331 560 : IF (const_mol(i) /= 0) THEN
1332 476 : k = const_mol(i)
1333 476 : 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 0 : " molecule kinds ("//cp_to_string(SIZE(molecule_kind_set))//").")
1338 476 : isize = SIZE(constr_x_mol(k)%constr)
1339 476 : CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1340 476 : constr_x_mol(k)%constr(isize + 1) = i
1341 : ELSE
1342 84 : myname = const_molname(i)
1343 84 : found_molname = .FALSE.
1344 304 : DO k = 1, SIZE(molecule_kind_set)
1345 220 : molecule_kind => molecule_kind_set(k)
1346 220 : name = molecule_kind%name
1347 220 : is_qm = qmmm_ff_precond_only_qm(id1=name)
1348 220 : IF (is_qm .AND. exclude_qm(i)) CYCLE
1349 152 : IF (.NOT. is_qm .AND. exclude_mm(i)) CYCLE
1350 292 : IF (name == myname) THEN
1351 82 : isize = SIZE(constr_x_mol(k)%constr)
1352 82 : CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1353 82 : constr_x_mol(k)%constr(isize + 1) = i
1354 82 : found_molname = .TRUE.
1355 : END IF
1356 : END DO
1357 84 : CALL print_warning_molname(found_molname, myname)
1358 : END IF
1359 : END IF
1360 : END DO
1361 316 : CALL timestop(handle)
1362 316 : 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 86 : SUBROUTINE print_warning_molname(found, name)
1371 : LOGICAL, INTENT(IN) :: found
1372 : CHARACTER(LEN=*), INTENT(IN) :: name
1373 :
1374 86 : 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 6 : "the print_key MOLECULES specific of the SUBSYS%PRINT section. ")
1381 :
1382 86 : END SUBROUTINE print_warning_molname
1383 :
1384 : END MODULE topology_constraint_util
|
