(git:ed6f26b)
Loading...
Searching...
No Matches
restraint.F
Go to the documentation of this file.
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 Handles all possible kinds of restraints in CP2K
10!> \author Teodoro Laino 08.2006 [tlaino] - University of Zurich
11!> \par History
12!> Teodoro Laino [tlaino] - 11.2008 : Improved the fixd_list restraints
13! **************************************************************************************************
14MODULE restraint
15 USE cell_types, ONLY: cell_type,&
16 use_perd_x,&
17 use_perd_xy,&
18 use_perd_xyz,&
19 use_perd_xz,&
20 use_perd_y,&
21 use_perd_yz,&
22 use_perd_z
23 USE colvar_methods, ONLY: colvar_eval_mol_f
24 USE colvar_types, ONLY: colvar_counters,&
25 diff_colvar
26 USE constraint_fxd, ONLY: create_local_fixd_list,&
27 release_local_fixd_list
28 USE cp_subsys_types, ONLY: cp_subsys_get,&
29 cp_subsys_type
30 USE distribution_1d_types, ONLY: distribution_1d_type
31 USE force_env_types, ONLY: force_env_get,&
32 force_env_set,&
33 force_env_type
34 USE kinds, ONLY: dp
35 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
36 USE molecule_kind_types, ONLY: colvar_constraint_type,&
37 fixd_constraint_type,&
38 g3x3_constraint_type,&
39 g4x6_constraint_type,&
40 get_molecule_kind,&
41 local_fixd_constraint_type,&
42 molecule_kind_type
43 USE molecule_list_types, ONLY: molecule_list_type
44 USE molecule_types, ONLY: get_molecule,&
45 global_constraint_type,&
46 local_colvar_constraint_type,&
47 molecule_type
48 USE particle_list_types, ONLY: particle_list_type
49 USE particle_types, ONLY: particle_type,&
50 update_particle_set
51#include "./base/base_uses.f90"
52
53 IMPLICIT NONE
54
55 PRIVATE
56 PUBLIC :: restraint_control
57 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'restraint'
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief Computes restraints
63!> \param force_env ...
64!> \author Teodoro Laino 08.2006 [tlaino]
65! **************************************************************************************************
66 SUBROUTINE restraint_control(force_env)
67
68 TYPE(force_env_type), POINTER :: force_env
69
70 CHARACTER(LEN=*), PARAMETER :: routinen = 'restraint_control'
71
72 INTEGER :: handle, i, ifixd, ii, ikind, imol, &
73 iparticle, n3x3con_restraint, &
74 n4x6con_restraint, n_restraint, nkind, &
75 nmol_per_kind
76 REAL(kind=dp) :: energy_3x3, energy_4x6, energy_colv, &
77 energy_fixd, extended_energies, k0, &
78 rab(3), rab2, targ(3)
79 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: force
80 TYPE(cell_type), POINTER :: cell
81 TYPE(colvar_counters) :: ncolv
82 TYPE(cp_subsys_type), POINTER :: subsys
83 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
84 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
85 TYPE(global_constraint_type), POINTER :: gci
86 TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:)
87 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
88 TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:)
89 TYPE(molecule_list_type), POINTER :: molecules
90 TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
91 TYPE(particle_list_type), POINTER :: particles
92 TYPE(particle_type), POINTER :: particle_set(:)
93
94 NULLIFY (cell, subsys, local_molecules, local_particles, fixd_list, molecule_kinds, &
95 molecules, molecule_kind, molecule_kind_set, molecule, molecule_set, particles, &
96 particle_set, gci, lfixd_list)
97 CALL timeset(routinen, handle)
98 CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
99 energy_4x6 = 0.0_dp
100 energy_3x3 = 0.0_dp
101 energy_colv = 0.0_dp
102 energy_fixd = 0.0_dp
103 n_restraint = 0
104 CALL cp_subsys_get(subsys=subsys, particles=particles, molecules=molecules, &
105 local_particles=local_particles, local_molecules=local_molecules, &
106 gci=gci, molecule_kinds=molecule_kinds)
107
108 nkind = molecule_kinds%n_els
109 particle_set => particles%els
110 molecule_set => molecules%els
111 molecule_kind_set => molecule_kinds%els
112
113 ! Intramolecular Restraints
114 ALLOCATE (force(3, SIZE(particle_set)))
115 force = 0.0_dp
116
117 ! Create the list of locally fixed atoms
118 CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
119
120 DO ifixd = 1, SIZE(lfixd_list)
121 ikind = lfixd_list(ifixd)%ikind
122 ii = lfixd_list(ifixd)%ifixd_index
123 molecule_kind => molecule_kind_set(ikind)
124 CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
125 IF (fixd_list(ii)%restraint%active) THEN
126 n_restraint = n_restraint + 1
127 iparticle = fixd_list(ii)%fixd
128 k0 = fixd_list(ii)%restraint%k0
129 targ = fixd_list(ii)%coord
130 rab = 0.0_dp
131 SELECT CASE (fixd_list(ii)%itype)
132 CASE (use_perd_x)
133 rab(1) = particle_set(iparticle)%r(1) - targ(1)
134 CASE (use_perd_y)
135 rab(2) = particle_set(iparticle)%r(2) - targ(2)
136 CASE (use_perd_z)
137 rab(3) = particle_set(iparticle)%r(3) - targ(3)
138 CASE (use_perd_xy)
139 rab(1) = particle_set(iparticle)%r(1) - targ(1)
140 rab(2) = particle_set(iparticle)%r(2) - targ(2)
141 CASE (use_perd_xz)
142 rab(1) = particle_set(iparticle)%r(1) - targ(1)
143 rab(3) = particle_set(iparticle)%r(3) - targ(3)
144 CASE (use_perd_yz)
145 rab(2) = particle_set(iparticle)%r(2) - targ(2)
146 rab(3) = particle_set(iparticle)%r(3) - targ(3)
147 CASE (use_perd_xyz)
148 rab = particle_set(iparticle)%r - targ
149 END SELECT
150 rab2 = dot_product(rab, rab)
151 ! Energy
152 energy_fixd = energy_fixd + k0*rab2
153 ! Forces
154 force(:, iparticle) = force(:, iparticle) - 2.0_dp*k0*rab
155 END IF
156 END DO
157 CALL release_local_fixd_list(lfixd_list)
158
159 ! Loop over other kind of Restraints
160 mol: DO ikind = 1, nkind
161 molecule_kind => molecule_kind_set(ikind)
162 nmol_per_kind = local_molecules%n_el(ikind)
163 DO imol = 1, nmol_per_kind
164 i = local_molecules%list(ikind)%array(imol)
165 molecule => molecule_set(i)
166 molecule_kind => molecule%molecule_kind
167
168 CALL get_molecule_kind(molecule_kind, &
169 ncolv=ncolv, &
170 ng3x3_restraint=n3x3con_restraint, &
171 ng4x6_restraint=n4x6con_restraint)
172 ! 3x3
173 IF (n3x3con_restraint /= 0) THEN
174 n_restraint = n_restraint + n3x3con_restraint
175 CALL restraint_3x3_int(molecule, particle_set, energy_3x3, force)
176 END IF
177 ! 4x6
178 IF (n4x6con_restraint /= 0) THEN
179 n_restraint = n_restraint + n4x6con_restraint
180 CALL restraint_4x6_int(molecule, particle_set, energy_4x6, force)
181 END IF
182 ! collective variables
183 IF (ncolv%nrestraint /= 0) THEN
184 n_restraint = n_restraint + ncolv%nrestraint
185 CALL restraint_colv_int(molecule, particle_set, cell, energy_colv, force)
186 END IF
187 END DO
188 END DO mol
189 CALL force_env%para_env%sum(n_restraint)
190 IF (n_restraint > 0) THEN
191 CALL force_env%para_env%sum(energy_fixd)
192 CALL force_env%para_env%sum(energy_3x3)
193 CALL force_env%para_env%sum(energy_4x6)
194 CALL force_env%para_env%sum(energy_colv)
195 CALL update_particle_set(particle_set, force_env%para_env, for=force, add=.true.)
196 force = 0.0_dp
197 n_restraint = 0
198 END IF
199 ! Intermolecular Restraints
200 IF (ASSOCIATED(gci)) THEN
201 IF (gci%nrestraint > 0) THEN
202 ! 3x3
203 IF (gci%ng3x3_restraint /= 0) THEN
204 n_restraint = n_restraint + gci%ng3x3_restraint
205 CALL restraint_3x3_ext(gci, particle_set, energy_3x3, force)
206 END IF
207 ! 4x6
208 IF (gci%ng4x6_restraint /= 0) THEN
209 n_restraint = n_restraint + gci%ng4x6_restraint
210 CALL restraint_4x6_ext(gci, particle_set, energy_4x6, force)
211 END IF
212 ! collective variables
213 IF (gci%ncolv%nrestraint /= 0) THEN
214 n_restraint = n_restraint + gci%ncolv%nrestraint
215 CALL restraint_colv_ext(gci, particle_set, cell, energy_colv, force)
216 END IF
217 DO iparticle = 1, SIZE(particle_set)
218 particle_set(iparticle)%f = particle_set(iparticle)%f + force(:, iparticle)
219 END DO
220 END IF
221 END IF
222 DEALLOCATE (force)
223
224 ! Store restraint energies
225 CALL force_env_get(force_env=force_env, additional_potential=extended_energies)
226 extended_energies = extended_energies + energy_3x3 + &
227 energy_fixd + &
228 energy_4x6 + &
229 energy_colv
230 CALL force_env_set(force_env=force_env, additional_potential=extended_energies)
231 CALL timestop(handle)
232
233 END SUBROUTINE restraint_control
234
235! **************************************************************************************************
236!> \brief Computes restraints 3x3 - Intramolecular
237!> \param molecule ...
238!> \param particle_set ...
239!> \param energy ...
240!> \param force ...
241!> \author Teodoro Laino 08.2006 [tlaino]
242! **************************************************************************************************
243 SUBROUTINE restraint_3x3_int(molecule, particle_set, energy, force)
244
245 TYPE(molecule_type), POINTER :: molecule
246 TYPE(particle_type), POINTER :: particle_set(:)
247 REAL(kind=dp), INTENT(INOUT) :: energy
248 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
249
250 INTEGER :: first_atom, ng3x3
251 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
252 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
253 TYPE(molecule_kind_type), POINTER :: molecule_kind
254
255 molecule_kind => molecule%molecule_kind
256 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, g3x3_list=g3x3_list, &
257 fixd_list=fixd_list)
258 CALL get_molecule(molecule, first_atom=first_atom)
259 CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, &
260 energy, force)
261
262 END SUBROUTINE restraint_3x3_int
263
264! **************************************************************************************************
265!> \brief Computes restraints 4x6 - Intramolecular
266!> \param molecule ...
267!> \param particle_set ...
268!> \param energy ...
269!> \param force ...
270!> \author Teodoro Laino 08.2006 [tlaino]
271! **************************************************************************************************
272 SUBROUTINE restraint_4x6_int(molecule, particle_set, energy, force)
273
274 TYPE(molecule_type), POINTER :: molecule
275 TYPE(particle_type), POINTER :: particle_set(:)
276 REAL(kind=dp), INTENT(INOUT) :: energy
277 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
278
279 INTEGER :: first_atom, ng4x6
280 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
281 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
282 TYPE(molecule_kind_type), POINTER :: molecule_kind
283
284 molecule_kind => molecule%molecule_kind
285 CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, g4x6_list=g4x6_list, &
286 fixd_list=fixd_list)
287 CALL get_molecule(molecule, first_atom=first_atom)
288 CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, &
289 energy, force)
290
291 END SUBROUTINE restraint_4x6_int
292
293! **************************************************************************************************
294!> \brief Computes restraints colv - Intramolecular
295!> \param molecule ...
296!> \param particle_set ...
297!> \param cell ...
298!> \param energy ...
299!> \param force ...
300!> \author Teodoro Laino 08.2006 [tlaino]
301! **************************************************************************************************
302 SUBROUTINE restraint_colv_int(molecule, particle_set, cell, energy, force)
303
304 TYPE(molecule_type), POINTER :: molecule
305 TYPE(particle_type), POINTER :: particle_set(:)
306 TYPE(cell_type), POINTER :: cell
307 REAL(kind=dp), INTENT(INOUT) :: energy
308 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
309
310 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
311 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
312 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
313 TYPE(molecule_kind_type), POINTER :: molecule_kind
314
315 NULLIFY (fixd_list)
316
317 molecule_kind => molecule%molecule_kind
318 CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
319 CALL get_molecule(molecule, lcolv=lcolv)
320 CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, &
321 cell, energy, force)
322
323 END SUBROUTINE restraint_colv_int
324
325! **************************************************************************************************
326!> \brief Computes restraints 3x3 - Intermolecular
327!> \param gci ...
328!> \param particle_set ...
329!> \param energy ...
330!> \param force ...
331!> \author Teodoro Laino 08.2006 [tlaino]
332! **************************************************************************************************
333 SUBROUTINE restraint_3x3_ext(gci, particle_set, energy, force)
334
335 TYPE(global_constraint_type), POINTER :: gci
336 TYPE(particle_type), POINTER :: particle_set(:)
337 REAL(kind=dp), INTENT(INOUT) :: energy
338 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
339
340 INTEGER :: first_atom, ng3x3
341 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
342 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
343
344 first_atom = 1
345 ng3x3 = gci%ng3x3
346 g3x3_list => gci%g3x3_list
347 fixd_list => gci%fixd_list
348 CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, &
349 energy, force)
350
351 END SUBROUTINE restraint_3x3_ext
352
353! **************************************************************************************************
354!> \brief Computes restraints 4x6 - Intermolecular
355!> \param gci ...
356!> \param particle_set ...
357!> \param energy ...
358!> \param force ...
359!> \author Teodoro Laino 08.2006 [tlaino]
360! **************************************************************************************************
361 SUBROUTINE restraint_4x6_ext(gci, particle_set, energy, force)
362
363 TYPE(global_constraint_type), POINTER :: gci
364 TYPE(particle_type), POINTER :: particle_set(:)
365 REAL(kind=dp), INTENT(INOUT) :: energy
366 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
367
368 INTEGER :: first_atom, ng4x6
369 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
370 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
371
372 first_atom = 1
373 ng4x6 = gci%ng4x6
374 g4x6_list => gci%g4x6_list
375 fixd_list => gci%fixd_list
376 CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, &
377 energy, force)
378
379 END SUBROUTINE restraint_4x6_ext
380
381! **************************************************************************************************
382!> \brief Computes restraints colv - Intermolecular
383!> \param gci ...
384!> \param particle_set ...
385!> \param cell ...
386!> \param energy ...
387!> \param force ...
388!> \author Teodoro Laino 08.2006 [tlaino]
389! **************************************************************************************************
390 SUBROUTINE restraint_colv_ext(gci, particle_set, cell, energy, force)
391
392 TYPE(global_constraint_type), POINTER :: gci
393 TYPE(particle_type), POINTER :: particle_set(:)
394 TYPE(cell_type), POINTER :: cell
395 REAL(kind=dp), INTENT(INOUT) :: energy
396 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
397
398 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
399 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
400 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
401
402 colv_list => gci%colv_list
403 fixd_list => gci%fixd_list
404 lcolv => gci%lcolv
405 CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, &
406 cell, energy, force)
407
408 END SUBROUTINE restraint_colv_ext
409
410! **************************************************************************************************
411!> \brief Computes restraints 3x3 - Real 3x3 restraints
412!> \param ng3x3 ...
413!> \param g3x3_list ...
414!> \param fixd_list ...
415!> \param first_atom ...
416!> \param particle_set ...
417!> \param energy ...
418!> \param force ...
419!> \author Teodoro Laino 08.2006 [tlaino]
420! **************************************************************************************************
421 SUBROUTINE restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, &
422 particle_set, energy, force)
423
424 INTEGER :: ng3x3
425 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
426 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
427 INTEGER, INTENT(IN) :: first_atom
428 TYPE(particle_type), POINTER :: particle_set(:)
429 REAL(kind=dp), INTENT(INOUT) :: energy
430 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
431
432 INTEGER :: iconst, index_a, index_b, index_c
433 REAL(kind=dp) :: k, rab, rac, rbc, tab, tac, tbc
434 REAL(kind=dp), DIMENSION(3) :: r0_12, r0_13, r0_23
435
436 DO iconst = 1, ng3x3
437 IF (.NOT. g3x3_list(iconst)%restraint%active) cycle
438 index_a = g3x3_list(iconst)%a + first_atom - 1
439 index_b = g3x3_list(iconst)%b + first_atom - 1
440 index_c = g3x3_list(iconst)%c + first_atom - 1
441 r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r
442 r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r
443 r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r
444
445 rab = sqrt(dot_product(r0_12, r0_12))
446 rac = sqrt(dot_product(r0_13, r0_13))
447 rbc = sqrt(dot_product(r0_23, r0_23))
448 tab = rab - g3x3_list(ng3x3)%dab
449 tac = rac - g3x3_list(ng3x3)%dac
450 tbc = rbc - g3x3_list(ng3x3)%dbc
451 k = g3x3_list(iconst)%restraint%k0
452 ! Update Energy
453 energy = energy + k*(tab**2 + tac**2 + tbc**2)
454 ! Update Forces
455 force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac)
456 force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc)
457 force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc)
458 ! Fixed atoms
459 IF (ASSOCIATED(fixd_list)) THEN
460 IF (SIZE(fixd_list) > 0) THEN
461 IF (any(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp
462 IF (any(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp
463 IF (any(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp
464 END IF
465 END IF
466 END DO
467
468 END SUBROUTINE restraint_3x3_low
469
470! **************************************************************************************************
471!> \brief Computes restraints 4x6 - Real 4x6 restraints
472!> \param ng4x6 ...
473!> \param g4x6_list ...
474!> \param fixd_list ...
475!> \param first_atom ...
476!> \param particle_set ...
477!> \param energy ...
478!> \param force ...
479!> \author Teodoro Laino 08.2006 [tlaino]
480! **************************************************************************************************
481 SUBROUTINE restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, &
482 particle_set, energy, force)
483
484 INTEGER :: ng4x6
485 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
486 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
487 INTEGER, INTENT(IN) :: first_atom
488 TYPE(particle_type), POINTER :: particle_set(:)
489 REAL(kind=dp), INTENT(INOUT) :: energy
490 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
491
492 INTEGER :: iconst, index_a, index_b, index_c, &
493 index_d
494 REAL(kind=dp) :: k, rab, rac, rad, rbc, rbd, rcd, tab, &
495 tac, tad, tbc, tbd, tcd
496 REAL(kind=dp), DIMENSION(3) :: r0_12, r0_13, r0_14, r0_23, r0_24, r0_34
497
498 DO iconst = 1, ng4x6
499 IF (.NOT. g4x6_list(iconst)%restraint%active) cycle
500 index_a = g4x6_list(iconst)%a + first_atom - 1
501 index_b = g4x6_list(iconst)%b + first_atom - 1
502 index_c = g4x6_list(iconst)%c + first_atom - 1
503 index_d = g4x6_list(iconst)%d + first_atom - 1
504 r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r
505 r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r
506 r0_14(:) = particle_set(index_a)%r - particle_set(index_d)%r
507 r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r
508 r0_24(:) = particle_set(index_b)%r - particle_set(index_d)%r
509 r0_34(:) = particle_set(index_c)%r - particle_set(index_d)%r
510
511 rab = sqrt(dot_product(r0_12, r0_12))
512 rac = sqrt(dot_product(r0_13, r0_13))
513 rad = sqrt(dot_product(r0_14, r0_14))
514 rbc = sqrt(dot_product(r0_23, r0_23))
515 rbd = sqrt(dot_product(r0_24, r0_24))
516 rcd = sqrt(dot_product(r0_34, r0_34))
517
518 tab = rab - g4x6_list(ng4x6)%dab
519 tac = rac - g4x6_list(ng4x6)%dac
520 tad = rad - g4x6_list(ng4x6)%dad
521 tbc = rbc - g4x6_list(ng4x6)%dbc
522 tbd = rbd - g4x6_list(ng4x6)%dbd
523 tcd = rcd - g4x6_list(ng4x6)%dcd
524
525 k = g4x6_list(iconst)%restraint%k0
526 ! Update Energy
527 energy = energy + k*(tab**2 + tac**2 + tad**2 + tbc**2 + tbd**2 + tcd**2)
528 ! Update Forces
529 force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac + r0_14/rad*tad)
530 force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc + r0_24/rbd*tbd)
531 force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc + r0_34/rcd*tcd)
532 force(:, index_d) = force(:, index_d) - 2.0_dp*k*(-r0_14/rad*tad - r0_24/rbd*tbd - r0_34/rcd*tcd)
533 ! Fixed atoms
534 IF (ASSOCIATED(fixd_list)) THEN
535 IF (SIZE(fixd_list) > 0) THEN
536 IF (any(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp
537 IF (any(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp
538 IF (any(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp
539 IF (any(fixd_list(:)%fixd == index_d)) force(:, index_d) = 0.0_dp
540 END IF
541 END IF
542 END DO
543
544 END SUBROUTINE restraint_4x6_low
545
546! **************************************************************************************************
547!> \brief Computes restraints colv - Real COLVAR restraints
548!> \param colv_list ...
549!> \param fixd_list ...
550!> \param lcolv ...
551!> \param particle_set ...
552!> \param cell ...
553!> \param energy ...
554!> \param force ...
555!> \author Teodoro Laino 08.2006 [tlaino]
556! **************************************************************************************************
557 SUBROUTINE restraint_colv_low(colv_list, fixd_list, lcolv, &
558 particle_set, cell, energy, force)
559
560 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
561 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
562 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
563 TYPE(particle_type), POINTER :: particle_set(:)
564 TYPE(cell_type), POINTER :: cell
565 REAL(kind=dp), INTENT(INOUT) :: energy
566 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
567
568 INTEGER :: iatm, iconst, ind
569 REAL(kind=dp) :: k, tab, targ
570
571 DO iconst = 1, SIZE(colv_list)
572 IF (.NOT. colv_list(iconst)%restraint%active) cycle
573 ! Update colvar
574 CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, &
575 particles=particle_set, fixd_list=fixd_list)
576
577 k = colv_list(iconst)%restraint%k0
578 targ = colv_list(iconst)%expected_value
579 tab = diff_colvar(lcolv(iconst)%colvar, targ)
580 ! Update Energy
581 energy = energy + k*tab**2
582 ! Update Forces
583 DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom)
584 ind = lcolv(iconst)%colvar%i_atom(iatm)
585 force(:, ind) = force(:, ind) - 2.0_dp*k*tab*lcolv(iconst)%colvar%dsdr(:, iatm)
586 END DO
587 END DO
588
589 END SUBROUTINE restraint_colv_low
590
591END MODULE restraint
for
for(int lxp=0;lxp<=lp;lxp++)
Definition grid_cpu_collint.h:66
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition cell_types.F:42
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition cell_types.F:42
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition cell_types.F:42
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition cell_types.F:42
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition cell_types.F:42
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition colvar_methods.F:16
colvar_methods::colvar_eval_mol_f
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Definition colvar_methods.F:1630
colvar_types
Initialize the collective variables types.
Definition colvar_types.F:15
colvar_types::diff_colvar
real(kind=dp) function, public diff_colvar(colvar, b)
subtract b from the ss value of a colvar: general function for handling periodic/non-periodic colvar
Definition colvar_types.F:1947
constraint_fxd
Definition constraint_fxd.F:12
constraint_fxd::release_local_fixd_list
subroutine, public release_local_fixd_list(lfixd_list)
destroy the list of local atoms on which to apply constraints/restraints Teodoro Laino [tlaino] - 11....
Definition constraint_fxd.F:542
constraint_fxd::create_local_fixd_list
subroutine, public create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
setup a list of local atoms on which to apply constraints/restraints
Definition constraint_fxd.F:423
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_types::force_env_set
subroutine, public force_env_set(force_env, meta_env, fp_env, force_env_section, method_name_id, additional_potential)
changes some attributes of the force_env
Definition force_env_types.F:666
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
molecule_list_types
represent a simple array based list of the given type
Definition molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::update_particle_set
subroutine, public update_particle_set(particle_set, int_group, pos, vel, for, add)
...
Definition particle_types.F:106
restraint
Handles all possible kinds of restraints in CP2K.
Definition restraint.F:14
restraint::restraint_control
subroutine, public restraint_control(force_env)
Computes restraints.
Definition restraint.F:67
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
colvar_types::colvar_counters
Definition colvar_types.F:388
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::colvar_constraint_type
Definition molecule_kind_types.F:117
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:149
molecule_kind_types::g3x3_constraint_type
Definition molecule_kind_types.F:128
molecule_kind_types::g4x6_constraint_type
Definition molecule_kind_types.F:135
molecule_kind_types::local_fixd_constraint_type
Definition molecule_kind_types.F:156
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::local_colvar_constraint_type
Definition molecule_types.F:41
molecule_types::molecule_type
Definition molecule_types.F:110
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35