Line data Source code
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,&
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 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Computes restraints
63 : !> \param force_env ...
64 : !> \author Teodoro Laino 08.2006 [tlaino]
65 : ! **************************************************************************************************
66 99502 : 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 99502 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
85 : TYPE(global_constraint_type), POINTER :: gci
86 99502 : TYPE(local_fixd_constraint_type), POINTER :: lfixd_list(:)
87 : TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
88 99502 : TYPE(molecule_kind_type), POINTER :: molecule_kind, molecule_kind_set(:)
89 : TYPE(molecule_list_type), POINTER :: molecules
90 99502 : TYPE(molecule_type), POINTER :: molecule, molecule_set(:)
91 : TYPE(particle_list_type), POINTER :: particles
92 99502 : TYPE(particle_type), POINTER :: particle_set(:)
93 :
94 99502 : NULLIFY (cell, subsys, local_molecules, local_particles, fixd_list, molecule_kinds, &
95 99502 : molecules, molecule_kind, molecule_kind_set, molecule, molecule_set, particles, &
96 99502 : particle_set, gci, lfixd_list)
97 99502 : CALL timeset(routineN, handle)
98 99502 : CALL force_env_get(force_env=force_env, subsys=subsys, cell=cell)
99 99502 : energy_4x6 = 0.0_dp
100 99502 : energy_3x3 = 0.0_dp
101 99502 : energy_colv = 0.0_dp
102 99502 : energy_fixd = 0.0_dp
103 99502 : n_restraint = 0
104 : CALL cp_subsys_get(subsys=subsys, particles=particles, molecules=molecules, &
105 : local_particles=local_particles, local_molecules=local_molecules, &
106 99502 : gci=gci, molecule_kinds=molecule_kinds)
107 :
108 99502 : nkind = molecule_kinds%n_els
109 99502 : particle_set => particles%els
110 99502 : molecule_set => molecules%els
111 99502 : molecule_kind_set => molecule_kinds%els
112 :
113 : ! Intramolecular Restraints
114 298506 : ALLOCATE (force(3, SIZE(particle_set)))
115 33489530 : force = 0.0_dp
116 :
117 : ! Create the list of locally fixed atoms
118 99502 : CALL create_local_fixd_list(lfixd_list, nkind, molecule_kind_set, local_particles)
119 :
120 156293 : DO ifixd = 1, SIZE(lfixd_list)
121 56791 : ikind = lfixd_list(ifixd)%ikind
122 56791 : ii = lfixd_list(ifixd)%ifixd_index
123 56791 : molecule_kind => molecule_kind_set(ikind)
124 56791 : CALL get_molecule_kind(molecule_kind, fixd_list=fixd_list)
125 156293 : IF (fixd_list(ii)%restraint%active) THEN
126 322 : n_restraint = n_restraint + 1
127 322 : iparticle = fixd_list(ii)%fixd
128 322 : k0 = fixd_list(ii)%restraint%k0
129 1288 : targ = fixd_list(ii)%coord
130 322 : rab = 0.0_dp
131 322 : SELECT CASE (fixd_list(ii)%itype)
132 : CASE (use_perd_x)
133 0 : rab(1) = particle_set(iparticle)%r(1) - targ(1)
134 : CASE (use_perd_y)
135 0 : rab(2) = particle_set(iparticle)%r(2) - targ(2)
136 : CASE (use_perd_z)
137 0 : rab(3) = particle_set(iparticle)%r(3) - targ(3)
138 : CASE (use_perd_xy)
139 0 : rab(1) = particle_set(iparticle)%r(1) - targ(1)
140 0 : rab(2) = particle_set(iparticle)%r(2) - targ(2)
141 : CASE (use_perd_xz)
142 0 : rab(1) = particle_set(iparticle)%r(1) - targ(1)
143 0 : rab(3) = particle_set(iparticle)%r(3) - targ(3)
144 : CASE (use_perd_yz)
145 0 : rab(2) = particle_set(iparticle)%r(2) - targ(2)
146 0 : rab(3) = particle_set(iparticle)%r(3) - targ(3)
147 : CASE (use_perd_xyz)
148 1288 : rab = particle_set(iparticle)%r - targ
149 : END SELECT
150 1288 : rab2 = DOT_PRODUCT(rab, rab)
151 : ! Energy
152 322 : energy_fixd = energy_fixd + k0*rab2
153 : ! Forces
154 1288 : force(:, iparticle) = force(:, iparticle) - 2.0_dp*k0*rab
155 : END IF
156 : END DO
157 99502 : CALL release_local_fixd_list(lfixd_list)
158 :
159 : ! Loop over other kind of Restraints
160 1849047 : MOL: DO ikind = 1, nkind
161 1749545 : molecule_kind => molecule_kind_set(ikind)
162 1749545 : nmol_per_kind = local_molecules%n_el(ikind)
163 3987711 : DO imol = 1, nmol_per_kind
164 2138664 : i = local_molecules%list(ikind)%array(imol)
165 2138664 : molecule => molecule_set(i)
166 2138664 : molecule_kind => molecule%molecule_kind
167 :
168 : CALL get_molecule_kind(molecule_kind, &
169 : ncolv=ncolv, &
170 : ng3x3_restraint=n3x3con_restraint, &
171 2138664 : ng4x6_restraint=n4x6con_restraint)
172 : ! 3x3
173 2138664 : IF (n3x3con_restraint /= 0) THEN
174 3 : n_restraint = n_restraint + n3x3con_restraint
175 3 : CALL restraint_3x3_int(molecule, particle_set, energy_3x3, force)
176 : END IF
177 : ! 4x6
178 2138664 : IF (n4x6con_restraint /= 0) THEN
179 6 : n_restraint = n_restraint + n4x6con_restraint
180 6 : CALL restraint_4x6_int(molecule, particle_set, energy_4x6, force)
181 : END IF
182 : ! collective variables
183 6026873 : IF (ncolv%nrestraint /= 0) THEN
184 1057 : n_restraint = n_restraint + ncolv%nrestraint
185 1057 : CALL restraint_colv_int(molecule, particle_set, cell, energy_colv, force)
186 : END IF
187 : END DO
188 : END DO MOL
189 99502 : CALL force_env%para_env%sum(n_restraint)
190 99502 : IF (n_restraint > 0) THEN
191 1530 : CALL force_env%para_env%sum(energy_fixd)
192 1530 : CALL force_env%para_env%sum(energy_3x3)
193 1530 : CALL force_env%para_env%sum(energy_4x6)
194 1530 : CALL force_env%para_env%sum(energy_colv)
195 1530 : CALL update_particle_set(particle_set, force_env%para_env, for=force, add=.TRUE.)
196 66498 : force = 0.0_dp
197 1530 : n_restraint = 0
198 : END IF
199 : ! Intermolecular Restraints
200 99502 : IF (ASSOCIATED(gci)) THEN
201 99502 : IF (gci%nrestraint > 0) THEN
202 : ! 3x3
203 912 : IF (gci%ng3x3_restraint /= 0) THEN
204 22 : n_restraint = n_restraint + gci%ng3x3_restraint
205 22 : CALL restraint_3x3_ext(gci, particle_set, energy_3x3, force)
206 : END IF
207 : ! 4x6
208 912 : IF (gci%ng4x6_restraint /= 0) THEN
209 12 : n_restraint = n_restraint + gci%ng4x6_restraint
210 12 : CALL restraint_4x6_ext(gci, particle_set, energy_4x6, force)
211 : END IF
212 : ! collective variables
213 912 : IF (gci%ncolv%nrestraint /= 0) THEN
214 878 : n_restraint = n_restraint + gci%ncolv%nrestraint
215 878 : CALL restraint_colv_ext(gci, particle_set, cell, energy_colv, force)
216 : END IF
217 86686 : DO iparticle = 1, SIZE(particle_set)
218 344008 : particle_set(iparticle)%f = particle_set(iparticle)%f + force(:, iparticle)
219 : END DO
220 : END IF
221 : END IF
222 99502 : DEALLOCATE (force)
223 :
224 : ! Store restraint energies
225 99502 : 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 99502 : energy_colv
230 99502 : CALL force_env_set(force_env=force_env, additional_potential=extended_energies)
231 99502 : CALL timestop(handle)
232 :
233 199004 : 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 3 : 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 3 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
252 3 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
253 : TYPE(molecule_kind_type), POINTER :: molecule_kind
254 :
255 3 : molecule_kind => molecule%molecule_kind
256 : CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, g3x3_list=g3x3_list, &
257 3 : fixd_list=fixd_list)
258 3 : CALL get_molecule(molecule, first_atom=first_atom)
259 : CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, &
260 3 : energy, force)
261 :
262 3 : 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 6 : 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 6 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
281 6 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
282 : TYPE(molecule_kind_type), POINTER :: molecule_kind
283 :
284 6 : molecule_kind => molecule%molecule_kind
285 : CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, g4x6_list=g4x6_list, &
286 6 : fixd_list=fixd_list)
287 6 : CALL get_molecule(molecule, first_atom=first_atom)
288 : CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, &
289 6 : energy, force)
290 :
291 6 : 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 1057 : 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 1057 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
311 1057 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
312 1057 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
313 : TYPE(molecule_kind_type), POINTER :: molecule_kind
314 :
315 1057 : NULLIFY (fixd_list)
316 :
317 1057 : molecule_kind => molecule%molecule_kind
318 1057 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
319 1057 : CALL get_molecule(molecule, lcolv=lcolv)
320 : CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, &
321 1057 : cell, energy, force)
322 :
323 1057 : 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 22 : 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 22 : first_atom = 1
345 22 : ng3x3 = gci%ng3x3
346 22 : g3x3_list => gci%g3x3_list
347 22 : fixd_list => gci%fixd_list
348 : CALL restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, particle_set, &
349 22 : energy, force)
350 :
351 22 : 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 12 : 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 12 : first_atom = 1
373 12 : ng4x6 = gci%ng4x6
374 12 : g4x6_list => gci%g4x6_list
375 12 : fixd_list => gci%fixd_list
376 : CALL restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, particle_set, &
377 12 : energy, force)
378 :
379 12 : 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 878 : 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 878 : colv_list => gci%colv_list
403 878 : fixd_list => gci%fixd_list
404 878 : lcolv => gci%lcolv
405 : CALL restraint_colv_low(colv_list, fixd_list, lcolv, particle_set, &
406 878 : cell, energy, force)
407 :
408 878 : 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 25 : SUBROUTINE restraint_3x3_low(ng3x3, g3x3_list, fixd_list, first_atom, &
422 25 : 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 50 : DO iconst = 1, ng3x3
437 25 : IF (.NOT. g3x3_list(iconst)%restraint%active) CYCLE
438 25 : index_a = g3x3_list(iconst)%a + first_atom - 1
439 25 : index_b = g3x3_list(iconst)%b + first_atom - 1
440 25 : index_c = g3x3_list(iconst)%c + first_atom - 1
441 100 : r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r
442 100 : r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r
443 100 : r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r
444 :
445 100 : rab = SQRT(DOT_PRODUCT(r0_12, r0_12))
446 100 : rac = SQRT(DOT_PRODUCT(r0_13, r0_13))
447 100 : rbc = SQRT(DOT_PRODUCT(r0_23, r0_23))
448 25 : tab = rab - g3x3_list(ng3x3)%dab
449 25 : tac = rac - g3x3_list(ng3x3)%dac
450 25 : tbc = rbc - g3x3_list(ng3x3)%dbc
451 25 : k = g3x3_list(iconst)%restraint%k0
452 : ! Update Energy
453 25 : energy = energy + k*(tab**2 + tac**2 + tbc**2)
454 : ! Update Forces
455 100 : force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac)
456 100 : force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc)
457 100 : force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc)
458 : ! Fixed atoms
459 50 : IF (ASSOCIATED(fixd_list)) THEN
460 0 : IF (SIZE(fixd_list) > 0) THEN
461 0 : IF (ANY(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp
462 0 : IF (ANY(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp
463 0 : IF (ANY(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp
464 : END IF
465 : END IF
466 : END DO
467 :
468 25 : 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 18 : SUBROUTINE restraint_4x6_low(ng4x6, g4x6_list, fixd_list, first_atom, &
482 18 : 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 36 : DO iconst = 1, ng4x6
499 18 : IF (.NOT. g4x6_list(iconst)%restraint%active) CYCLE
500 18 : index_a = g4x6_list(iconst)%a + first_atom - 1
501 18 : index_b = g4x6_list(iconst)%b + first_atom - 1
502 18 : index_c = g4x6_list(iconst)%c + first_atom - 1
503 18 : index_d = g4x6_list(iconst)%d + first_atom - 1
504 72 : r0_12(:) = particle_set(index_a)%r - particle_set(index_b)%r
505 72 : r0_13(:) = particle_set(index_a)%r - particle_set(index_c)%r
506 72 : r0_14(:) = particle_set(index_a)%r - particle_set(index_d)%r
507 72 : r0_23(:) = particle_set(index_b)%r - particle_set(index_c)%r
508 72 : r0_24(:) = particle_set(index_b)%r - particle_set(index_d)%r
509 72 : r0_34(:) = particle_set(index_c)%r - particle_set(index_d)%r
510 :
511 72 : rab = SQRT(DOT_PRODUCT(r0_12, r0_12))
512 72 : rac = SQRT(DOT_PRODUCT(r0_13, r0_13))
513 72 : rad = SQRT(DOT_PRODUCT(r0_14, r0_14))
514 72 : rbc = SQRT(DOT_PRODUCT(r0_23, r0_23))
515 72 : rbd = SQRT(DOT_PRODUCT(r0_24, r0_24))
516 72 : rcd = SQRT(DOT_PRODUCT(r0_34, r0_34))
517 :
518 18 : tab = rab - g4x6_list(ng4x6)%dab
519 18 : tac = rac - g4x6_list(ng4x6)%dac
520 18 : tad = rad - g4x6_list(ng4x6)%dad
521 18 : tbc = rbc - g4x6_list(ng4x6)%dbc
522 18 : tbd = rbd - g4x6_list(ng4x6)%dbd
523 18 : tcd = rcd - g4x6_list(ng4x6)%dcd
524 :
525 18 : k = g4x6_list(iconst)%restraint%k0
526 : ! Update Energy
527 18 : energy = energy + k*(tab**2 + tac**2 + tad**2 + tbc**2 + tbd**2 + tcd**2)
528 : ! Update Forces
529 72 : force(:, index_a) = force(:, index_a) - 2.0_dp*k*(r0_12/rab*tab + r0_13/rac*tac + r0_14/rad*tad)
530 72 : force(:, index_b) = force(:, index_b) - 2.0_dp*k*(-r0_12/rab*tab + r0_23/rbc*tbc + r0_24/rbd*tbd)
531 72 : force(:, index_c) = force(:, index_c) - 2.0_dp*k*(-r0_13/rac*tac - r0_23/rbc*tbc + r0_34/rcd*tcd)
532 72 : 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 36 : IF (ASSOCIATED(fixd_list)) THEN
535 0 : IF (SIZE(fixd_list) > 0) THEN
536 0 : IF (ANY(fixd_list(:)%fixd == index_a)) force(:, index_a) = 0.0_dp
537 0 : IF (ANY(fixd_list(:)%fixd == index_b)) force(:, index_b) = 0.0_dp
538 0 : IF (ANY(fixd_list(:)%fixd == index_c)) force(:, index_c) = 0.0_dp
539 0 : IF (ANY(fixd_list(:)%fixd == index_d)) force(:, index_d) = 0.0_dp
540 : END IF
541 : END IF
542 : END DO
543 :
544 18 : 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 1935 : SUBROUTINE restraint_colv_low(colv_list, fixd_list, lcolv, &
558 1935 : 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 4483 : DO iconst = 1, SIZE(colv_list)
572 2548 : IF (.NOT. colv_list(iconst)%restraint%active) CYCLE
573 : ! Update colvar
574 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, &
575 2504 : particles=particle_set, fixd_list=fixd_list)
576 :
577 2504 : k = colv_list(iconst)%restraint%k0
578 2504 : targ = colv_list(iconst)%expected_value
579 2504 : tab = diff_colvar(lcolv(iconst)%colvar, targ)
580 : ! Update Energy
581 2504 : energy = energy + k*tab**2
582 : ! Update Forces
583 12887 : DO iatm = 1, SIZE(lcolv(iconst)%colvar%i_atom)
584 8448 : ind = lcolv(iconst)%colvar%i_atom(iatm)
585 36340 : force(:, ind) = force(:, ind) - 2.0_dp*k*tab*lcolv(iconst)%colvar%dsdr(:, iatm)
586 : END DO
587 : END DO
588 :
589 1935 : END SUBROUTINE restraint_colv_low
590 :
591 : END MODULE restraint
|
