(git:58e3e09)
restraint.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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 ! **************************************************************************************************
14 MODULE 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,&
24  USE colvar_types, ONLY: colvar_counters,&
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,&
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,&
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,&
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 
59 CONTAINS
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 
591 END MODULE restraint
for(int lxp=0;lxp<=lp;lxp++)
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
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.
Definition: colvar_types.F:15
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
subroutine, public release_local_fixd_list(lfixd_list)
destroy the list of local atoms on which to apply constraints/restraints Teodoro Laino [tlaino] - 11....
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
types that represent a subsys, i.e. a part of the system
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
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Interface for the force calculations.
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)
returns various attributes about the force environment
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
represent a simple array based list of the given type
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.
represent a simple array based list of the given type
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.
represent a simple array based list of the given type
Define the data structure for the particle information.
subroutine, public update_particle_set(particle_set, int_group, pos, vel, for, add)
...
Handles all possible kinds of restraints in CP2K.
Definition: restraint.F:14
subroutine, public restraint_control(force_env)
Computes restraints.
Definition: restraint.F:67