(git:21ef868)
Loading...
Searching...
No Matches
topology_constraint_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-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! **************************************************************************************************
26 USE colvar_types, ONLY: &
43 USE kinds, ONLY: default_string_length,&
44 dp
46 USE molecule_kind_types, ONLY: &
50 USE molecule_types, ONLY: get_molecule,&
64#include "./base/base_uses.f90"
65
66 IMPLICIT NONE
67
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_constraint_util'
69
70 PRIVATE
72
73CONTAINS
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 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 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 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 LOGICAL, ALLOCATABLE, DIMENSION(:) :: missed_molname
111 REAL(kind=dp) :: rmod, rvec(3)
112 REAL(kind=dp), DIMENSION(:), POINTER :: hdist, r
113 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
114 TYPE(atomic_kind_type), POINTER :: atomic_kind
115 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list
116 TYPE(colvar_constraint_type), DIMENSION(:), &
117 POINTER :: colv_list
118 TYPE(colvar_counters) :: ncolv
119 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list, fixd_list_gci
123 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
124 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
125 TYPE(local_colvar_constraint_type), DIMENSION(:), &
126 POINTER :: lcolv
127 TYPE(local_constraint_type), POINTER :: lci
128 TYPE(local_g3x3_constraint_type), DIMENSION(:), &
129 POINTER :: lg3x3
130 TYPE(local_g4x6_constraint_type), DIMENSION(:), &
131 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 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list
137
138 NULLIFY (logger, constr_x_mol, constr_x_glob)
139 logger => cp_get_default_logger()
140 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
141 extension=".subsysLog")
142 CALL timeset(routinen, handle)
143 CALL timeset(routinen//"_1", handle2)
144
145 cons_info => topology%cons_info
146 hbonds_section => section_vals_get_subs_vals(input_file, &
147 "MOTION%CONSTRAINT%HBONDS")
148 fixd_restr_rest => section_vals_get_subs_vals(input_file, &
149 "MOTION%CONSTRAINT%FIX_ATOM_RESTART")
150 CALL section_vals_get(fixd_restr_rest, explicit=restart_restraint_pos)
151 colvar_rest => section_vals_get_subs_vals(input_file, &
152 "MOTION%CONSTRAINT%COLVAR_RESTART")
153 CALL section_vals_get(colvar_rest, explicit=restart_restraint_clv)
154 colvar_func_info => section_vals_get_subs_vals(subsys_section, &
155 "COLVAR%COLVAR_FUNC_INFO")
156 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 DO i = 1, topology%nmol
162 molecule => molecule_set(i)
163 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 topology%const_colv .OR. topology%const_vsite) THEN
168 ALLOCATE (lci)
169 NULLIFY (lci%lcolv)
170 NULLIFY (lci%lg3x3)
171 NULLIFY (lci%lg4x6)
172 END IF
173 CALL set_molecule(molecule, lci=lci)
174 END DO
175 ALLOCATE (gci)
176 NULLIFY (gci%lcolv, &
177 gci%lg3x3, &
178 gci%lg4x6, &
179 gci%fixd_list, &
180 gci%colv_list, &
181 gci%g3x3_list, &
182 gci%g4x6_list, &
183 gci%vsite_list)
184 gci%ntot = 0
185 gci%ng3x3 = 0
186 gci%ng4x6 = 0
187 gci%nvsite = 0
188 gci%ng3x3_restraint = 0
189 gci%ng4x6_restraint = 0
190 gci%nvsite_restraint = 0
191 CALL setup_colvar_counters(gci%colv_list, gci%ncolv)
192 gci%nrestraint = gci%ng3x3_restraint + &
193 gci%ng4x6_restraint + &
194 gci%nvsite_restraint + &
195 gci%ncolv%nrestraint
196 CALL timestop(handle2)
197 CALL timeset(routinen//"_2", handle2)
198 !-----------------------------------------------------------------------------
199 !-----------------------------------------------------------------------------
200 ! 2. Add more stuff to COLVAR constraint if constraint hydrogen is on
201 !-----------------------------------------------------------------------------
202 IF (topology%const_hydr) THEN
203 topology%const_colv = .true.
204 NULLIFY (atom_typeh, hdist)
205 ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
206 DO i = 1, SIZE(molecule_kind_set)
207 ALLOCATE (constr_x_mol(i)%constr(1))
208 constr_x_mol(i)%constr(1) = 1
209 END DO
210 CALL section_vals_val_get(hbonds_section, "MOLECULE", n_rep_val=nrep)
211 IF (nrep /= 0) THEN
212 NULLIFY (inds)
213 DO i = 1, SIZE(molecule_kind_set)
214 constr_x_mol(i)%constr(1) = 0
215 END DO
216 CALL section_vals_val_get(hbonds_section, "MOLECULE", i_vals=inds)
217 DO i = 1, SIZE(inds)
218 constr_x_mol(inds(i))%constr(1) = 1
219 END DO
220 ELSE
221 CALL section_vals_val_get(hbonds_section, "MOLNAME", n_rep_val=nrep)
222 IF (nrep /= 0) THEN
223 NULLIFY (cnds)
224 DO i = 1, SIZE(molecule_kind_set)
225 constr_x_mol(i)%constr(1) = 0
226 END DO
227 CALL section_vals_val_get(hbonds_section, "MOLNAME", c_vals=cnds)
228 DO i = 1, SIZE(cnds)
229 found_molname = .false.
230 DO k = 1, SIZE(molecule_kind_set)
231 molecule_kind => molecule_kind_set(k)
232 name = molecule_kind%name
233 ldummy = qmmm_ff_precond_only_qm(id1=name)
234 IF (cnds(i) == name) THEN
235 constr_x_mol(k)%constr(1) = 1
236 found_molname = .true.
237 END IF
238 END DO
239 CALL print_warning_molname(found_molname, cnds(i))
240 END DO
241 END IF
242 END IF
243 CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", n_rep_val=nrep)
244 IF (nrep /= 0) THEN
245 CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", c_vals=atom_typeh)
246 END IF
247 CALL section_vals_val_get(hbonds_section, "TARGETS", n_rep_val=nrep)
248 IF (nrep /= 0) THEN
249 CALL section_vals_val_get(hbonds_section, "TARGETS", r_vals=hdist)
250 END IF
251 IF (ASSOCIATED(hdist)) THEN
252 cpassert(SIZE(hdist) == SIZE(atom_typeh))
253 END IF
254 CALL section_vals_val_get(hbonds_section, "exclude_qm", l_val=exclude_qm)
255 CALL section_vals_val_get(hbonds_section, "exclude_mm", l_val=exclude_mm)
256 nhdist = 0
257 DO i = 1, SIZE(molecule_kind_set)
258 molecule_kind => molecule_kind_set(i)
259 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 molecule_list=molecule_list)
263 ! Let's tag all requested atoms involving Hydrogen
264 ! on the first molecule of this kind
265 molecule => molecule_set(molecule_list(1))
266 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
267 natom = last_atom - first_atom + 1
268 DO k = 1, nbond
269 ishbond = .false.
270 j = bond_list(k)%a
271 IF (j < 1 .OR. j > natom) cycle
272 atomic_kind => atom_list(j)%atomic_kind
273 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
274 is_qm = qmmm_ff_precond_only_qm(id1=name)
275 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
276 IF (is_qm .AND. exclude_qm) ishbond = .false.
277 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
278 IF (.NOT. ishbond) THEN
279 j = bond_list(k)%b
280 IF (j < 1 .OR. j > natom) cycle
281 atomic_kind => atom_list(j)%atomic_kind
282 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
283 is_qm = qmmm_ff_precond_only_qm(id1=name)
284 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
285 IF (is_qm .AND. exclude_qm) ishbond = .false.
286 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
287 END IF
288 IF (ishbond) THEN
289 nhdist = nhdist + 1
290 END IF
291 END DO
292 END DO
293 n_start_colv = cons_info%nconst_colv
294 cons_info%nconst_colv = nhdist + n_start_colv
295 CALL reallocate(cons_info%const_colv_mol, 1, cons_info%nconst_colv)
296 CALL reallocate(cons_info%const_colv_molname, 1, cons_info%nconst_colv)
297 CALL reallocate(cons_info%const_colv_target, 1, cons_info%nconst_colv)
298 CALL reallocate(cons_info%const_colv_target_growth, 1, cons_info%nconst_colv)
299 CALL colvar_p_reallocate(cons_info%colvar_set, 1, cons_info%nconst_colv)
300 ! Fill in Restraints info
301 CALL reallocate(cons_info%colv_intermolecular, 1, cons_info%nconst_colv)
302 CALL reallocate(cons_info%colv_restraint, 1, cons_info%nconst_colv)
303 CALL reallocate(cons_info%colv_k0, 1, cons_info%nconst_colv)
304 CALL reallocate(cons_info%colv_exclude_qm, 1, cons_info%nconst_colv)
305 CALL reallocate(cons_info%colv_exclude_mm, 1, cons_info%nconst_colv)
306 ! Bonds involving hydrogens are by their nature only intramolecular
307 cons_info%colv_intermolecular(n_start_colv + 1:cons_info%nconst_colv) = .false.
308 cons_info%colv_exclude_qm(n_start_colv + 1:cons_info%nconst_colv) = .false.
309 cons_info%colv_exclude_mm(n_start_colv + 1:cons_info%nconst_colv) = .false.
310 cons_info%colv_restraint(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_restraint
311 cons_info%colv_k0(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_k0
312 !
313 nhdist = 0
314 DO i = 1, SIZE(molecule_kind_set)
315 IF (constr_x_mol(i)%constr(1) == 0) cycle
316 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 molecule_list=molecule_list)
320 molecule => molecule_set(molecule_list(1))
321 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
322 natom = last_atom - first_atom + 1
323 offset = first_atom - 1
324 DO k = 1, nbond
325 ishbond = .false.
326 j = bond_list(k)%a
327 IF (j < 1 .OR. j > natom) cycle
328 atomic_kind => atom_list(j)%atomic_kind
329 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
330 is_qm = qmmm_ff_precond_only_qm(id1=name)
331 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
332 IF (is_qm .AND. exclude_qm) ishbond = .false.
333 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
334 IF (.NOT. ishbond) THEN
335 j = bond_list(k)%b
336 IF (j < 1 .OR. j > natom) cycle
337 atomic_kind => atom_list(j)%atomic_kind
338 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name)
339 is_qm = qmmm_ff_precond_only_qm(id1=name)
340 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .true.
341 IF (is_qm .AND. exclude_qm) ishbond = .false.
342 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .false.
343 END IF
344 IF (ishbond) THEN
345 nhdist = nhdist + 1
346 rvec = particle_set(offset + bond_list(k)%a)%r - particle_set(offset + bond_list(k)%b)%r
347 rmod = norm2(rvec)
348 IF (ASSOCIATED(hdist)) THEN
349 IF (SIZE(hdist) > 0) THEN
350 IF (bond_list(k)%a == j) atomic_kind => atom_list(bond_list(k)%b)%atomic_kind
351 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 name=name, element_symbol=element_symbol)
354 ldummy = qmmm_ff_precond_only_qm(id1=name)
355 DO m = 1, SIZE(hdist)
356 IF (trim(name) == trim(atom_typeh(m))) EXIT
357 IF (trim(element_symbol) == trim(atom_typeh(m))) EXIT
358 END DO
359 IF (m <= SIZE(hdist)) THEN
360 rmod = hdist(m)
361 END IF
362 END IF
363 END IF
364 cons_info%const_colv_mol(nhdist + n_start_colv) = i
365 cons_info%const_colv_molname(nhdist + n_start_colv) = "UNDEF"
366 cons_info%const_colv_target(nhdist + n_start_colv) = rmod
367 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, &
370 cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%i_at = bond_list(k)%a
371 cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%j_at = bond_list(k)%b
372 CALL colvar_setup(cons_info%colvar_set(nhdist + n_start_colv)%colvar)
373 END IF
374 END DO
375 END DO
376 DO j = 1, SIZE(constr_x_mol)
377 DEALLOCATE (constr_x_mol(j)%constr)
378 END DO
379 DEALLOCATE (constr_x_mol)
380 END IF
381
382 CALL timestop(handle2)
383 CALL timeset(routinen//"_3", handle2)
384 !-----------------------------------------------------------------------------
385 !-----------------------------------------------------------------------------
386 ! 3. Set the COLVAR constraint molecule_kind_set(ikind)%colv_list
387 !-----------------------------------------------------------------------------
388 IF (topology%const_colv) THEN
389 ! Post Process of COLVARS..
390 DO ii = 1, SIZE(cons_info%colvar_set)
391 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 cons_info%colv_exclude_mm)
402 ! Intramolecular constraints
403 gind = 0
404 cind = 0
405 DO ii = 1, SIZE(molecule_kind_set)
406 molecule_kind => molecule_kind_set(ii)
407 CALL get_molecule_kind(molecule_kind=molecule_kind, &
408 nmolecule=nmolecule, molecule_list=molecule_list)
409 ncolv_mol = SIZE(constr_x_mol(ii)%constr)
410 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 molecule => molecule_set(molecule_list(1))
416 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 colvar_rest, first_atom)
420 CALL setup_colvar_counters(colv_list, ncolv)
421 CALL set_molecule_kind(molecule_kind, colv_list=colv_list, ncolv=ncolv)
422 DO j = 1, nmolecule
423 molecule => molecule_set(molecule_list(j))
424 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
425 ALLOCATE (lcolv(ncolv_mol))
426 CALL setup_lcolv(lcolv, constr_x_mol(ii)%constr, first_atom, last_atom, &
427 cons_info, particle_set, colvar_func_info, use_clv_info, cind)
428 CALL set_molecule(molecule=molecule, lcolv=lcolv)
429 END DO
430 END DO
431 DO j = 1, SIZE(constr_x_mol)
432 DEALLOCATE (constr_x_mol(j)%constr)
433 END DO
434 DEALLOCATE (constr_x_mol)
435 ! Intermolecular constraints
436 ncolv_glob = 0
437 IF (ASSOCIATED(constr_x_glob)) THEN
438 ncolv_glob = SIZE(constr_x_glob)
439 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 first_atom=1)
443 CALL setup_colvar_counters(colv_list, ncolv)
444 ALLOCATE (lcolv(ncolv_glob))
445 CALL setup_lcolv(lcolv, constr_x_glob, 1, SIZE(particle_set), cons_info, &
446 particle_set, colvar_func_info, use_clv_info, cind)
447 gci%colv_list => colv_list
448 gci%lcolv => lcolv
449 gci%ncolv = ncolv
450 ! Total number of Intermolecular constraints
451 gci%ntot = gci%ncolv%ntot + gci%ntot
452 DEALLOCATE (constr_x_glob)
453 END IF
454 END IF
455
456 CALL timestop(handle2)
457 CALL timeset(routinen//"_4", handle2)
458 !-----------------------------------------------------------------------------
459 !-----------------------------------------------------------------------------
460 ! 4. Set the group 3x3 constraint g3x3_list
461 !-----------------------------------------------------------------------------
462 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 cons_info%g33_exclude_mm)
471 ! Intramolecular constraints
472 DO ii = 1, SIZE(molecule_kind_set)
473 molecule_kind => molecule_kind_set(ii)
474 CALL get_molecule_kind(molecule_kind=molecule_kind, &
475 nmolecule=nmolecule, &
476 molecule_list=molecule_list)
477 ng3x3 = SIZE(constr_x_mol(ii)%constr)
478 ALLOCATE (g3x3_list(ng3x3))
479 CALL setup_g3x3_list(g3x3_list, constr_x_mol(ii)%constr, cons_info, ng3x3_restraint)
480 CALL set_molecule_kind(molecule_kind, ng3x3=ng3x3, ng3x3_restraint=ng3x3_restraint, g3x3_list=g3x3_list)
481 DO j = 1, nmolecule
482 molecule => molecule_set(molecule_list(j))
483 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
484 ALLOCATE (lg3x3(ng3x3))
485 CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
486 CALL set_molecule(molecule=molecule, lg3x3=lg3x3)
487 END DO
488 END DO
489 DO j = 1, SIZE(constr_x_mol)
490 DEALLOCATE (constr_x_mol(j)%constr)
491 END DO
492 DEALLOCATE (constr_x_mol)
493 ! Intermolecular constraints
494 IF (ASSOCIATED(constr_x_glob)) THEN
495 ng3x3 = SIZE(constr_x_glob)
496 ALLOCATE (g3x3_list(ng3x3))
497 CALL setup_g3x3_list(g3x3_list, constr_x_glob, cons_info, ng3x3_restraint)
498 ALLOCATE (lg3x3(ng3x3))
499 CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom)
500 gci%g3x3_list => g3x3_list
501 gci%lg3x3 => lg3x3
502 gci%ng3x3 = ng3x3
503 gci%ng3x3_restraint = ng3x3_restraint
504 ! Total number of Intermolecular constraints
505 gci%ntot = 3*gci%ng3x3 + gci%ntot
506 DEALLOCATE (constr_x_glob)
507 END IF
508 END IF
509
510 CALL timestop(handle2)
511 CALL timeset(routinen//"_5", handle2)
512 !-----------------------------------------------------------------------------
513 !-----------------------------------------------------------------------------
514 ! 5. Set the group 4x6 constraint g4x6_list
515 !-----------------------------------------------------------------------------
516 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 cons_info%g46_exclude_mm)
525 ! Intramolecular constraints
526 DO ii = 1, SIZE(molecule_kind_set)
527 molecule_kind => molecule_kind_set(ii)
528 CALL get_molecule_kind(molecule_kind=molecule_kind, &
529 nmolecule=nmolecule, molecule_list=molecule_list)
530 ng4x6 = SIZE(constr_x_mol(ii)%constr)
531 ALLOCATE (g4x6_list(ng4x6))
532 CALL setup_g4x6_list(g4x6_list, constr_x_mol(ii)%constr, cons_info, ng4x6_restraint)
533 CALL set_molecule_kind(molecule_kind, ng4x6=ng4x6, ng4x6_restraint=ng4x6_restraint, g4x6_list=g4x6_list)
534 DO j = 1, nmolecule
535 molecule => molecule_set(molecule_list(j))
536 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
537 ALLOCATE (lg4x6(ng4x6))
538 CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
539 CALL set_molecule(molecule=molecule, lg4x6=lg4x6)
540 END DO
541 END DO
542 DO j = 1, SIZE(constr_x_mol)
543 DEALLOCATE (constr_x_mol(j)%constr)
544 END DO
545 DEALLOCATE (constr_x_mol)
546 ! Intermolecular constraints
547 IF (ASSOCIATED(constr_x_glob)) THEN
548 ng4x6 = SIZE(constr_x_glob)
549 ALLOCATE (g4x6_list(ng4x6))
550 CALL setup_g4x6_list(g4x6_list, constr_x_glob, cons_info, ng4x6_restraint)
551 ALLOCATE (lg4x6(ng4x6))
552 CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom)
553 gci%g4x6_list => g4x6_list
554 gci%lg4x6 => lg4x6
555 gci%ng4x6 = ng4x6
556 gci%ng4x6_restraint = ng4x6_restraint
557 ! Total number of Intermolecular constraints
558 gci%ntot = 6*gci%ng4x6 + gci%ntot
559 DEALLOCATE (constr_x_glob)
560 END IF
561 END IF
562
563 CALL timestop(handle2)
564 CALL timeset(routinen//"_6", handle2)
565 !-----------------------------------------------------------------------------
566 !-----------------------------------------------------------------------------
567 ! 6. Set the group vsite constraint vsite_list
568 !-----------------------------------------------------------------------------
569 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 cons_info%vsite_exclude_mm)
578 ! Intramolecular constraints
579 DO ii = 1, SIZE(molecule_kind_set)
580 molecule_kind => molecule_kind_set(ii)
581 CALL get_molecule_kind(molecule_kind=molecule_kind, &
582 nmolecule=nmolecule, molecule_list=molecule_list)
583 nvsite = SIZE(constr_x_mol(ii)%constr)
584 ALLOCATE (vsite_list(nvsite))
585 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 vsite_list=vsite_list)
588 END DO
589 DO j = 1, SIZE(constr_x_mol)
590 DEALLOCATE (constr_x_mol(j)%constr)
591 END DO
592 DEALLOCATE (constr_x_mol)
593 ! Intermolecular constraints
594 IF (ASSOCIATED(constr_x_glob)) THEN
595 nvsite = SIZE(constr_x_glob)
596 ALLOCATE (vsite_list(nvsite))
597 CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint)
598 gci%vsite_list => vsite_list
599 gci%nvsite = nvsite
600 gci%nvsite_restraint = nvsite_restraint
601 ! Total number of Intermolecular constraints
602 gci%ntot = gci%nvsite + gci%ntot
603 DEALLOCATE (constr_x_glob)
604 END IF
605 END IF
606 CALL timestop(handle2)
607 CALL timeset(routinen//"_7", handle2)
608 !-----------------------------------------------------------------------------
609 !-----------------------------------------------------------------------------
610 ! 7. Set the group fixed_atom constraint fixd_list
611 !-----------------------------------------------------------------------------
612 IF (topology%const_atom) THEN
613 ALLOCATE (fixd_list_gci(SIZE(particle_set)))
614 nfixd_list_gci = 0
615 ALLOCATE (missed_molname(SIZE(cons_info%fixed_molnames, 1)))
616 missed_molname = .true.
617 nfixd_restart = 0
618 DO i = 1, SIZE(molecule_kind_set)
619 molecule_kind => molecule_kind_set(i)
620 CALL get_molecule_kind(molecule_kind=molecule_kind, &
621 nmolecule=nmolecule, molecule_list=molecule_list, name=molname)
622 is_qm = qmmm_ff_precond_only_qm(id1=molname)
623 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 nfixed_atoms = 0
628 DO j = 1, nmolecule
629 molecule => molecule_set(molecule_list(j))
630 CALL get_molecule(molecule, first_atom=first, last_atom=last)
631 fix_atom_molname = .false.
632 IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
633 DO k = 1, SIZE(cons_info%fixed_molnames)
634 IF (cons_info%fixed_molnames(k) == molname) THEN
635 fix_atom_molname = .true.
636 IF (is_qm .AND. cons_info%fixed_exclude_qm(k)) fix_atom_molname = .false.
637 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 DO k = first, last
642 fix_atom_qmmm = .false.
643 IF (PRESENT(qmmm_env)) THEN
644 SELECT CASE (cons_info%freeze_qm)
645 CASE (do_constr_atomic)
646 IF (any(qmmm_env%qm_atom_index == k)) fix_atom_qmmm = .true.
647 CASE (do_constr_molec)
648 IF (any(qmmm_env%qm_molecule_index == molecule_list(j))) fix_atom_qmmm = .true.
649 END SELECT
650 SELECT CASE (cons_info%freeze_mm)
651 CASE (do_constr_atomic)
652 IF (all(qmmm_env%qm_atom_index /= k)) fix_atom_qmmm = .true.
653 CASE (do_constr_molec)
654 IF (all(qmmm_env%qm_molecule_index /= molecule_list(j))) fix_atom_qmmm = .true.
655 END SELECT
656 END IF
657 IF (any(cons_info%fixed_atoms == k) .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
658 nfixed_atoms = nfixed_atoms + 1
659 END IF
660 END DO
661 END DO
662 ALLOCATE (fixd_list(nfixed_atoms))
663 kk = 0
664 nfixd_restraint = 0
665 IF (nfixed_atoms /= 0) THEN
666 DO j = 1, nmolecule
667 molecule => molecule_set(molecule_list(j))
668 CALL get_molecule(molecule, first_atom=first, last_atom=last)
669 fix_atom_molname = .false.
670 IF (ASSOCIATED(cons_info%fixed_molnames)) THEN
671 DO k1loc = 1, SIZE(cons_info%fixed_molnames)
672 IF (cons_info%fixed_molnames(k1loc) == molname) THEN
673 fix_atom_molname = .true.
674 itype = cons_info%fixed_mol_type(k1loc)
675 EXIT
676 END IF
677 END DO
678 END IF
679 DO k = first, last
680 ! FIXED LIST ATOMS
681 fix_fixed_atom = .false.
682 DO k2loc = 1, SIZE(cons_info%fixed_atoms)
683 IF (cons_info%fixed_atoms(k2loc) == k) THEN
684 fix_fixed_atom = .true.
685 itype = cons_info%fixed_type(k2loc)
686 EXIT
687 END IF
688 END DO
689 ! QMMM FIXED ATOMS (QM OR MM)
690 fix_atom_qmmm = .false.
691 fix_atom_mm = .false.
692 fix_atom_qm = .false.
693 IF (PRESENT(qmmm_env)) THEN
694 SELECT CASE (cons_info%freeze_qm)
695 CASE (do_constr_atomic)
696 IF (any(qmmm_env%qm_atom_index == k)) THEN
697 fix_atom_qmmm = .true.
698 fix_atom_qm = .true.
699 itype = cons_info%freeze_qm_type
700 END IF
701 CASE (do_constr_molec)
702 IF (any(qmmm_env%qm_molecule_index == molecule_list(j))) THEN
703 fix_atom_qmmm = .true.
704 fix_atom_qm = .true.
705 itype = cons_info%freeze_qm_type
706 END IF
707 END SELECT
708 SELECT CASE (cons_info%freeze_mm)
709 CASE (do_constr_atomic)
710 IF (all(qmmm_env%qm_atom_index /= k)) THEN
711 fix_atom_qmmm = .true.
712 fix_atom_mm = .true.
713 itype = cons_info%freeze_mm_type
714 END IF
715 CASE (do_constr_molec)
716 IF (all(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN
717 fix_atom_qmmm = .true.
718 fix_atom_mm = .true.
719 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 IF (fix_atom_qm .AND. fix_atom_mm) THEN
724 CALL cp_abort(__location__, &
725 "Atom number: "//cp_to_string(k)// &
726 " 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 .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 " input section. Check and correct your input file!")
736 END IF
737 ! Let's store the atom index
738 IF (fix_fixed_atom .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN
739 IF (ASSOCIATED(topology%cell_muc)) THEN
740 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 "disable CELL%CANONICALIZE for this input.")
745 END IF
746 END IF
747 kk = kk + 1
748 fixd_list(kk)%fixd = k
749 fixd_list(kk)%coord = particle_set(k)%r
750 fixd_list(kk)%itype = itype
751 ! Possibly Restraint
752 IF (fix_fixed_atom) THEN
753 fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc)
754 fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc)
755 ELSE IF (fix_atom_qm) THEN
756 fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint
757 fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0
758 ELSE IF (fix_atom_mm) THEN
759 fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint
760 fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0
761 ELSE IF (fix_atom_molname) THEN
762 fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc)
763 fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc)
764 ELSE
765 ! Should never reach this point
766 cpabort("Unknown fix atom specification")
767 END IF
768 IF (fixd_list(kk)%restraint%active) THEN
769 nfixd_restraint = nfixd_restraint + 1
770 nfixd_restart = nfixd_restart + 1
771 ! Check that we use the components that we really want..
772 SELECT CASE (itype)
773 CASE (use_perd_x)
774 fixd_list(kk)%coord(2) = huge(0.0_dp)
775 fixd_list(kk)%coord(3) = huge(0.0_dp)
776 CASE (use_perd_y)
777 fixd_list(kk)%coord(1) = huge(0.0_dp)
778 fixd_list(kk)%coord(3) = huge(0.0_dp)
779 CASE (use_perd_z)
780 fixd_list(kk)%coord(1) = huge(0.0_dp)
781 fixd_list(kk)%coord(2) = huge(0.0_dp)
782 CASE (use_perd_xy)
783 fixd_list(kk)%coord(3) = huge(0.0_dp)
784 CASE (use_perd_xz)
785 fixd_list(kk)%coord(2) = huge(0.0_dp)
786 CASE (use_perd_yz)
787 fixd_list(kk)%coord(1) = huge(0.0_dp)
788 END SELECT
789 IF (restart_restraint_pos) THEN
790 ! Read coord0 value for restraint
791 CALL section_vals_val_get(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
792 i_rep_val=nfixd_restart, r_vals=r)
793 SELECT CASE (itype)
794 CASE (use_perd_x)
795 cpassert(SIZE(r) == 1)
796 fixd_list(kk)%coord(1) = r(1)
797 CASE (use_perd_y)
798 cpassert(SIZE(r) == 1)
799 fixd_list(kk)%coord(2) = r(1)
800 CASE (use_perd_z)
801 cpassert(SIZE(r) == 1)
802 fixd_list(kk)%coord(3) = r(1)
803 CASE (use_perd_xy)
804 cpassert(SIZE(r) == 2)
805 fixd_list(kk)%coord(1) = r(1)
806 fixd_list(kk)%coord(2) = r(2)
807 CASE (use_perd_xz)
808 cpassert(SIZE(r) == 2)
809 fixd_list(kk)%coord(1) = r(1)
810 fixd_list(kk)%coord(3) = r(2)
811 CASE (use_perd_yz)
812 cpassert(SIZE(r) == 2)
813 fixd_list(kk)%coord(2) = r(1)
814 fixd_list(kk)%coord(3) = r(2)
815 CASE (use_perd_xyz)
816 cpassert(SIZE(r) == 3)
817 fixd_list(kk)%coord(1:3) = r(1:3)
818 IF (ASSOCIATED(topology%cell_muc)) THEN
819 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 SELECT CASE (itype)
825 CASE (use_perd_x)
826 ALLOCATE (r(1))
827 r(1) = fixd_list(kk)%coord(1)
828 CASE (use_perd_y)
829 ALLOCATE (r(1))
830 r(1) = fixd_list(kk)%coord(2)
831 CASE (use_perd_z)
832 ALLOCATE (r(1))
833 r(1) = fixd_list(kk)%coord(3)
834 CASE (use_perd_xy)
835 ALLOCATE (r(2))
836 r(1) = fixd_list(kk)%coord(1)
837 r(2) = fixd_list(kk)%coord(2)
838 CASE (use_perd_xz)
839 ALLOCATE (r(2))
840 r(1) = fixd_list(kk)%coord(1)
841 r(2) = fixd_list(kk)%coord(3)
842 CASE (use_perd_yz)
843 ALLOCATE (r(2))
844 r(1) = fixd_list(kk)%coord(1)
845 r(2) = fixd_list(kk)%coord(3)
846 CASE (use_perd_xyz)
847 ALLOCATE (r(3))
848 r(1:3) = fixd_list(kk)%coord(1:3)
849 END SELECT
850 CALL section_vals_val_set(fixd_restr_rest, "_DEFAULT_KEYWORD_", &
851 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 IF (iw > 0) THEN
859 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 fixd_list=fixd_list)
863 fixd_list_gci(nfixd_list_gci + 1:nfixd_list_gci + nfixed_atoms) = fixd_list
864 nfixd_list_gci = nfixd_list_gci + nfixed_atoms
865 END DO
866 IF (iw > 0) THEN
867 WRITE (iw, *) "TOTAL NUMBER OF FIXED ATOMS:", nfixd_list_gci
868 END IF
869 cpassert(count(missed_molname) == 0)
870 DEALLOCATE (missed_molname)
871 ! Intermolecular constraints
872 IF (gci%ntot /= 0) THEN
873 ALLOCATE (fixd_list(nfixd_list_gci))
874 fixd_list(1:nfixd_list_gci) = fixd_list_gci(1:nfixd_list_gci)
875 gci%fixd_list => fixd_list
876 END IF
877 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 gci%ncolv%nrestraint
884 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
885 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
886 CALL timestop(handle2)
887 CALL timestop(handle)
888 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 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 ncolv_mol = 0
925 DO kk = 1, SIZE(ilist)
926 j = ilist(kk)
927 ncolv_mol = ncolv_mol + 1
928 kdim = SIZE(cons_info%colvar_set(j)%colvar%i_atom)
929 ALLOCATE (colv_list(ncolv_mol)%i_atoms(kdim))
930 colv_list(ncolv_mol)%inp_seq_num = j
931 colv_list(ncolv_mol)%type_id = cons_info%colvar_set(j)%colvar%type_id
932 colv_list(ncolv_mol)%i_atoms = cons_info%colvar_set(j)%colvar%i_atom
933 colv_list(ncolv_mol)%use_points = cons_info%colvar_set(j)%colvar%use_points
934 ! Restraint
935 colv_list(ncolv_mol)%restraint%active = cons_info%colv_restraint(j)
936 colv_list(ncolv_mol)%restraint%k0 = cons_info%colv_k0(j)
937 IF (cons_info%const_colv_target(j) == -huge(0.0_dp)) THEN
938 ! Let's compute the value..
939 NULLIFY (local_colvar)
940 CALL colvar_clone(local_colvar, cons_info%colvar_set(j)%colvar, &
941 i_atom_offset=first_atom - 1)
942 CALL colvar_eval_mol_f(local_colvar, topology%cell, particle_set)
943 colv_list(ncolv_mol)%expected_value = local_colvar%ss
944 CALL colvar_release(local_colvar)
945 ELSE
946 colv_list(ncolv_mol)%expected_value = cons_info%const_colv_target(j)
947 END IF
948 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 IF (colv_list(ncolv_mol)%restraint%active .AND. &
951 (colv_list(ncolv_mol)%expected_value_growth_speed == 0.0_dp)) THEN
952 gind = gind + 1
953 IF (restart_restraint_clv) THEN
954 CALL section_vals_val_get(colvar_rest, "_DEFAULT_KEYWORD_", &
955 i_rep_val=gind, r_val=rmod)
956 colv_list(ncolv_mol)%expected_value = rmod
957 ELSE
958 rmod = colv_list(ncolv_mol)%expected_value
959 CALL section_vals_val_set(colvar_rest, "_DEFAULT_KEYWORD_", &
960 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 IF (cons_info%colvar_set(j)%colvar%type_id == torsion_colvar_id) THEN
966 cons_info%colvar_set(j)%colvar%torsion_param%o0 = colv_list(ncolv_mol)%expected_value
967 END IF
968 END DO
969 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 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 ng3x3_restraint = 0
990 DO ng3x3 = 1, SIZE(ilist)
991 j = ilist(ng3x3)
992 g3x3_list(ng3x3)%a = cons_info%const_g33_a(j)
993 g3x3_list(ng3x3)%b = cons_info%const_g33_b(j)
994 g3x3_list(ng3x3)%c = cons_info%const_g33_c(j)
995 g3x3_list(ng3x3)%dab = cons_info%const_g33_dab(j)
996 g3x3_list(ng3x3)%dac = cons_info%const_g33_dac(j)
997 g3x3_list(ng3x3)%dbc = cons_info%const_g33_dbc(j)
998 ! Restraint
999 g3x3_list(ng3x3)%restraint%active = cons_info%g33_restraint(j)
1000 g3x3_list(ng3x3)%restraint%k0 = cons_info%g33_k0(j)
1001 IF (g3x3_list(ng3x3)%restraint%active) ng3x3_restraint = ng3x3_restraint + 1
1002 END DO
1003
1004 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 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 ng4x6 = 0
1025 ng4x6_restraint = 0
1026 DO ng4x6 = 1, SIZE(ilist)
1027 j = ilist(ng4x6)
1028 g4x6_list(ng4x6)%a = cons_info%const_g46_a(j)
1029 g4x6_list(ng4x6)%b = cons_info%const_g46_b(j)
1030 g4x6_list(ng4x6)%c = cons_info%const_g46_c(j)
1031 g4x6_list(ng4x6)%d = cons_info%const_g46_d(j)
1032 g4x6_list(ng4x6)%dab = cons_info%const_g46_dab(j)
1033 g4x6_list(ng4x6)%dac = cons_info%const_g46_dac(j)
1034 g4x6_list(ng4x6)%dbc = cons_info%const_g46_dbc(j)
1035 g4x6_list(ng4x6)%dad = cons_info%const_g46_dad(j)
1036 g4x6_list(ng4x6)%dbd = cons_info%const_g46_dbd(j)
1037 g4x6_list(ng4x6)%dcd = cons_info%const_g46_dcd(j)
1038 ! Restraint
1039 g4x6_list(ng4x6)%restraint%active = cons_info%g46_restraint(j)
1040 g4x6_list(ng4x6)%restraint%k0 = cons_info%g46_k0(j)
1041 IF (g4x6_list(ng4x6)%restraint%active) ng4x6_restraint = ng4x6_restraint + 1
1042 END DO
1043
1044 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 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 nvsite = 0
1064 nvsite_restraint = 0
1065 DO nvsite = 1, SIZE(ilist)
1066 j = ilist(nvsite)
1067 vsite_list(nvsite)%a = cons_info%const_vsite_a(j)
1068 vsite_list(nvsite)%b = cons_info%const_vsite_b(j)
1069 vsite_list(nvsite)%c = cons_info%const_vsite_c(j)
1070 vsite_list(nvsite)%d = cons_info%const_vsite_d(j)
1071 vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j)
1072 vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j)
1073 ! Restraint
1074 vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j)
1075 vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j)
1076 IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint + 1
1077 END DO
1078
1079 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: r_vals
1111
1112 DO kk = 1, SIZE(ilist)
1113 k = ilist(kk)
1114 lcolv(kk)%init = .false.
1115 lcolv(kk)%lambda = 0.0_dp
1116 lcolv(kk)%sigma = 0.0_dp
1117
1118 ! Set Up colvar variable
1119 NULLIFY (lcolv(kk)%colvar, lcolv(kk)%colvar_old)
1120 ! Colvar
1121 CALL colvar_clone(lcolv(kk)%colvar, cons_info%colvar_set(k)%colvar, &
1122 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 IF ((lcolv(kk)%colvar%type_id == xyz_diag_colvar_id) .OR. &
1129 (lcolv(kk)%colvar%type_id == xyz_outerdiag_colvar_id)) THEN
1130 cind = cind + 1
1131 IF (use_clv_info) THEN
1132 CALL section_vals_val_get(colvar_func_info, "_DEFAULT_KEYWORD_", &
1133 i_rep_val=cind, r_vals=r_vals)
1134 SELECT CASE (lcolv(kk)%colvar%type_id)
1135 CASE (xyz_diag_colvar_id)
1136 cpassert(SIZE(r_vals) == 3)
1137 lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1139 cpassert(SIZE(r_vals) == 6)
1140 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1141 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6)
1142 END SELECT
1143 ELSE
1144 SELECT CASE (lcolv(kk)%colvar%type_id)
1145 CASE (xyz_diag_colvar_id)
1146 ALLOCATE (r_vals(3))
1147 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_diag_param%i_atom
1148 r_vals = particle_set(ind)%r
1149 lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals
1151 ALLOCATE (r_vals(6))
1152 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(1)
1153 r_vals(1:3) = particle_set(ind)%r
1154 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(2)
1155 r_vals(4:6) = particle_set(ind)%r
1156 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3)
1157 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 i_rep_val=cind, r_vals_ptr=r_vals)
1161 END IF
1162 END IF
1163
1164 ! Setup Colvar_old
1165 CALL colvar_clone(lcolv(kk)%colvar_old, lcolv(kk)%colvar)
1166
1167 ! Check for consistency in the constraint definition
1168 IF (any(lcolv(kk)%colvar%i_atom > last_atom) .OR. &
1169 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 "in the input file. Please check it carefully!")
1177 END IF
1178 END DO
1179 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 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 DO kk = 1, SIZE(lg3x3)
1200 lg3x3(kk)%init = .false.
1201 lg3x3(kk)%scale = 0.0_dp
1202 lg3x3(kk)%scale_old = 0.0_dp
1203 lg3x3(kk)%fa = 0.0_dp
1204 lg3x3(kk)%fb = 0.0_dp
1205 lg3x3(kk)%fc = 0.0_dp
1206 lg3x3(kk)%ra_old = 0.0_dp
1207 lg3x3(kk)%rb_old = 0.0_dp
1208 lg3x3(kk)%rc_old = 0.0_dp
1209 lg3x3(kk)%va = 0.0_dp
1210 lg3x3(kk)%vb = 0.0_dp
1211 lg3x3(kk)%vc = 0.0_dp
1212 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 (g3x3_list(kk)%b + first_atom - 1 > last_atom) .OR. &
1218 (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 "in the input file. Please check it carefully!")
1226 END IF
1227 END DO
1228
1229 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 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 DO kk = 1, SIZE(lg4x6)
1250 lg4x6(kk)%init = .false.
1251 lg4x6(kk)%scale = 0.0_dp
1252 lg4x6(kk)%scale_old = 0.0_dp
1253 lg4x6(kk)%fa = 0.0_dp
1254 lg4x6(kk)%fb = 0.0_dp
1255 lg4x6(kk)%fc = 0.0_dp
1256 lg4x6(kk)%fd = 0.0_dp
1257 lg4x6(kk)%fe = 0.0_dp
1258 lg4x6(kk)%ff = 0.0_dp
1259 lg4x6(kk)%ra_old = 0.0_dp
1260 lg4x6(kk)%rb_old = 0.0_dp
1261 lg4x6(kk)%rc_old = 0.0_dp
1262 lg4x6(kk)%rd_old = 0.0_dp
1263 lg4x6(kk)%re_old = 0.0_dp
1264 lg4x6(kk)%rf_old = 0.0_dp
1265 lg4x6(kk)%va = 0.0_dp
1266 lg4x6(kk)%vb = 0.0_dp
1267 lg4x6(kk)%vc = 0.0_dp
1268 lg4x6(kk)%vd = 0.0_dp
1269 lg4x6(kk)%ve = 0.0_dp
1270 lg4x6(kk)%vf = 0.0_dp
1271 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 (g4x6_list(kk)%c + first_atom - 1 > last_atom) .OR. &
1279 (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 "in the input file. Please check it carefully!")
1287 END IF
1288 END DO
1289
1290 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 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 CALL timeset(routinen, handle)
1326 NULLIFY (molecule_kind)
1327 ALLOCATE (constr_x_mol(SIZE(molecule_kind_set)))
1328 DO i = 1, SIZE(constr_x_mol)
1329 NULLIFY (constr_x_mol(i)%constr)
1330 ALLOCATE (constr_x_mol(i)%constr(0))
1331 END DO
1332 cpassert(SIZE(const_mol) == SIZE(const_molname))
1333 iglob = 0
1334 DO i = 1, SIZE(const_mol)
1335 IF (const_intermolecular(i)) THEN
1336 ! Intermolecular constraint
1337 iglob = iglob + 1
1338 CALL reallocate(constr_x_glob, 1, iglob)
1339 constr_x_glob(iglob) = i
1340 ELSE
1341 ! Intramolecular constraint
1342 IF (const_mol(i) /= 0) THEN
1343 k = const_mol(i)
1344 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 " molecule kinds ("//cp_to_string(SIZE(molecule_kind_set))//").")
1349 END IF
1350 isize = SIZE(constr_x_mol(k)%constr)
1351 CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1352 constr_x_mol(k)%constr(isize + 1) = i
1353 ELSE
1354 myname = const_molname(i)
1355 found_molname = .false.
1356 DO k = 1, SIZE(molecule_kind_set)
1357 molecule_kind => molecule_kind_set(k)
1358 name = molecule_kind%name
1359 is_qm = qmmm_ff_precond_only_qm(id1=name)
1360 IF (is_qm .AND. exclude_qm(i)) cycle
1361 IF (.NOT. is_qm .AND. exclude_mm(i)) cycle
1362 IF (name == myname) THEN
1363 isize = SIZE(constr_x_mol(k)%constr)
1364 CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1)
1365 constr_x_mol(k)%constr(isize + 1) = i
1366 found_molname = .true.
1367 END IF
1368 END DO
1369 CALL print_warning_molname(found_molname, myname)
1370 END IF
1371 END IF
1372 END DO
1373 CALL timestop(handle)
1374 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 SUBROUTINE print_warning_molname(found, name)
1383 LOGICAL, INTENT(IN) :: found
1384 CHARACTER(LEN=*), INTENT(IN) :: name
1385
1386 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 "the print_key MOLECULES specific of the SUBSYS%PRINT section. ")
1393 END IF
1394
1395 END SUBROUTINE print_warning_molname
1396
1397END MODULE topology_constraint_util
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
elemental logical function, public is_hydrogen(atomic_kind)
Determines if the atomic_kind is HYDROGEN.
Handles all functions related to the CELL.
Definition cell_types.F:15
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
integer, parameter, public use_perd_y
Definition cell_types.F:42
integer, parameter, public use_perd_xz
Definition cell_types.F:42
integer, parameter, public use_perd_x
Definition cell_types.F:42
integer, parameter, public use_perd_z
Definition cell_types.F:42
integer, parameter, public use_perd_yz
Definition cell_types.F:42
integer, parameter, public use_perd_xy
Definition cell_types.F:42
subroutine, public cell_transform_input_cartesian(cell, vector)
Transform a Cartesian real-space vector from the user input cell frame into CP2K's canonical internal...
Definition cell_types.F:284
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Initialize the collective variables types.
recursive subroutine, public colvar_release(colvar)
releases the memory that might have been allocated by the colvar
integer, parameter, public xyz_outerdiag_colvar_id
subroutine, public colvar_create(colvar, colvar_id)
initializes a colvar_param type
recursive subroutine, public colvar_clone(colvar_out, colvar_in, i_atom_offset)
Clone a colvar type.
subroutine, public colvar_p_reallocate(colvar_set, lb1_new, ub1_new)
Change the dimension of a colvar_p_type.
integer, parameter, public dist_colvar_id
subroutine, public colvar_setup(colvar)
Finalize the setup of the collective variable.
integer, parameter, public torsion_colvar_id
integer, parameter, public xyz_diag_colvar_id
evaluations of colvar for internal coordinates schemes
recursive subroutine, public post_process_colvar(colvar, particles)
Complete the description of the COORDINATION colvar when defined using KINDS.
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_constr_atomic
integer, parameter, public do_constr_molec
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Utility routines for the memory handling.
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
subroutine, public setup_colvar_counters(colv_list, ncolv)
...
subroutine, public set_molecule_kind(molecule_kind, name, mass, charge, kind_number, molecule_list, atom_list, nbond, bond_list, nbend, bend_list, nub, ub_list, nimpr, impr_list, nopbend, opbend_list, ntorsion, torsion_list, fixd_list, ncolv, colv_list, ng3x3, g3x3_list, ng4x6, nfixd, g4x6_list, nvsite, vsite_list, ng3x3_restraint, ng4x6_restraint, nfixd_restraint, nshell, shell_list, nvsite_restraint, bond_kind_set, bend_kind_set, ub_kind_set, torsion_kind_set, impr_kind_set, opbend_kind_set, nelectron, nsgf, molname_generated)
Set the components of a molecule kind.
Define the data structure for the molecule information.
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
subroutine, public set_molecule(molecule, molecule_kind, lmi, lci, lcolv, lg3x3, lg4x6)
Set a molecule data set.
Define the data structure for the particle information.
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
Collection of subroutine needed for topology related things.
subroutine, public topology_constraint_pack(molecule_kind_set, molecule_set, topology, qmmm_env, particle_set, input_file, subsys_section, gci)
Pack in all the information needed for the constraints.
Control for reading in different topologies and coordinates.
Definition topology.F:13
Provides all information about an atomic kind.
parameters for a collective variable
type of a logger, at the moment it contains just a print level starting at which level it should be l...