Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Module that handles the COLLECTIVE constraints
10 : !> \par History
11 : !> none
12 : !> \author Teodoro Laino [tlaino]
13 : ! **************************************************************************************************
14 : MODULE constraint_clv
15 : USE cell_types, ONLY: cell_type
16 : USE colvar_methods, ONLY: colvar_eval_mol_f
17 : USE colvar_types, ONLY: colvar_type,&
18 : diff_colvar
19 : USE input_section_types, ONLY: section_vals_get,&
20 : section_vals_get_subs_vals,&
21 : section_vals_type,&
22 : section_vals_val_get,&
23 : section_vals_val_set
24 : USE kinds, ONLY: dp
25 : USE molecule_kind_types, ONLY: colvar_constraint_type,&
26 : fixd_constraint_type,&
27 : get_molecule_kind,&
28 : molecule_kind_type
29 : USE molecule_types, ONLY: get_molecule,&
30 : global_constraint_type,&
31 : local_colvar_constraint_type,&
32 : molecule_type
33 : USE particle_types, ONLY: particle_type
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 : PUBLIC :: shake_roll_colv_int, &
40 : rattle_roll_colv_int, &
41 : shake_colv_int, &
42 : rattle_colv_int, &
43 : shake_roll_colv_ext, &
44 : rattle_roll_colv_ext, &
45 : shake_colv_ext, &
46 : rattle_colv_ext, &
47 : shake_update_colv_int, &
48 : shake_update_colv_ext
49 :
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_clv'
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Intramolecular subroutine
56 : !> shake_colv algorithm for collective variables constraints
57 : !> updates the multiplier one molecule type at a time
58 : !> \param molecule ...
59 : !> \param particle_set ...
60 : !> \param pos ...
61 : !> \param vel ...
62 : !> \param dt ...
63 : !> \param ishake ...
64 : !> \param cell ...
65 : !> \param imass ...
66 : !> \param max_sigma ...
67 : !> \par History
68 : !> none
69 : !> \author Teodoro Laino [tlaino]
70 : ! **************************************************************************************************
71 18698 : SUBROUTINE shake_colv_int(molecule, particle_set, pos, vel, dt, ishake, &
72 18698 : cell, imass, max_sigma)
73 :
74 : TYPE(molecule_type), POINTER :: molecule
75 : TYPE(particle_type), POINTER :: particle_set(:)
76 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
77 : REAL(kind=dp), INTENT(in) :: dt
78 : INTEGER, INTENT(IN) :: ishake
79 : TYPE(cell_type), POINTER :: cell
80 : REAL(KIND=dp), DIMENSION(:) :: imass
81 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
82 :
83 18698 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
84 18698 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
85 18698 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
86 : TYPE(molecule_kind_type), POINTER :: molecule_kind
87 :
88 18698 : NULLIFY (fixd_list)
89 18698 : molecule_kind => molecule%molecule_kind
90 18698 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
91 18698 : CALL get_molecule(molecule, lcolv=lcolv)
92 : ! Real Shake
93 : CALL shake_colv_low(fixd_list, colv_list, lcolv, &
94 18698 : particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
95 :
96 18698 : END SUBROUTINE shake_colv_int
97 :
98 : ! **************************************************************************************************
99 : !> \brief Intramolecular subroutine for updating the TARGET value of collective
100 : !> constraints
101 : !> \param molecule ...
102 : !> \param dt ...
103 : !> \param motion_section ...
104 : !> \author Teodoro Laino [tlaino] - University of Zurich
105 : ! **************************************************************************************************
106 2591 : SUBROUTINE shake_update_colv_int(molecule, dt, motion_section)
107 :
108 : TYPE(molecule_type), POINTER :: molecule
109 : REAL(kind=dp), INTENT(in) :: dt
110 : TYPE(section_vals_type), POINTER :: motion_section
111 :
112 2591 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
113 : TYPE(molecule_kind_type), POINTER :: molecule_kind
114 :
115 2591 : molecule_kind => molecule%molecule_kind
116 2591 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list)
117 : ! Real update of the Shake target
118 2591 : CALL shake_update_colv_low(colv_list, dt, motion_section)
119 :
120 2591 : END SUBROUTINE shake_update_colv_int
121 :
122 : ! **************************************************************************************************
123 : !> \brief Intramolecular subroutine
124 : !> rattle algorithm for collective variables constraints
125 : !> updates the multiplier one molecule type at a time
126 : !> \param molecule ...
127 : !> \param particle_set ...
128 : !> \param vel ...
129 : !> \param dt ...
130 : !> \param irattle ...
131 : !> \param cell ...
132 : !> \param imass ...
133 : !> \param max_sigma ...
134 : !> \par History
135 : !> none
136 : !> \author Teodoro Laino [tlaino]
137 : ! **************************************************************************************************
138 7833 : SUBROUTINE rattle_colv_int(molecule, particle_set, vel, dt, irattle, &
139 7833 : cell, imass, max_sigma)
140 :
141 : TYPE(molecule_type), POINTER :: molecule
142 : TYPE(particle_type), POINTER :: particle_set(:)
143 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
144 : REAL(kind=dp), INTENT(in) :: dt
145 : INTEGER, INTENT(IN) :: irattle
146 : TYPE(cell_type), POINTER :: cell
147 : REAL(KIND=dp), DIMENSION(:) :: imass
148 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
149 :
150 7833 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
151 7833 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
152 7833 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
153 : TYPE(molecule_kind_type), POINTER :: molecule_kind
154 :
155 7833 : NULLIFY (fixd_list)
156 7833 : molecule_kind => molecule%molecule_kind
157 7833 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
158 7833 : CALL get_molecule(molecule, lcolv=lcolv)
159 : ! Real Rattle
160 : CALL rattle_colv_low(fixd_list, colv_list, lcolv, &
161 7833 : particle_set, vel, dt, irattle, cell, imass, max_sigma)
162 :
163 7833 : END SUBROUTINE rattle_colv_int
164 :
165 : ! **************************************************************************************************
166 : !> \brief Intramolecular subroutine
167 : !> shake algorithm (box allowed to change) for collective variables constraints
168 : !> updates the multiplier one molecule type at a time
169 : !> \param molecule ...
170 : !> \param particle_set ...
171 : !> \param pos ...
172 : !> \param vel ...
173 : !> \param r_shake ...
174 : !> \param v_shake ...
175 : !> \param dt ...
176 : !> \param ishake ...
177 : !> \param cell ...
178 : !> \param imass ...
179 : !> \param max_sigma ...
180 : !> \par History
181 : !> none
182 : !> \author Teodoro Laino [tlaino]
183 : ! **************************************************************************************************
184 15773 : SUBROUTINE shake_roll_colv_int(molecule, particle_set, pos, vel, r_shake, v_shake, &
185 15773 : dt, ishake, cell, imass, max_sigma)
186 :
187 : TYPE(molecule_type), POINTER :: molecule
188 : TYPE(particle_type), POINTER :: particle_set(:)
189 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
190 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
191 : REAL(kind=dp), INTENT(in) :: dt
192 : INTEGER, INTENT(in) :: ishake
193 : TYPE(cell_type), POINTER :: cell
194 : REAL(KIND=dp), DIMENSION(:) :: imass
195 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
196 :
197 15773 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
198 15773 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
199 15773 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
200 : TYPE(molecule_kind_type), POINTER :: molecule_kind
201 :
202 15773 : NULLIFY (fixd_list)
203 15773 : molecule_kind => molecule%molecule_kind
204 15773 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
205 15773 : CALL get_molecule(molecule, lcolv=lcolv)
206 : ! Real Shake
207 : CALL shake_roll_colv_low(fixd_list, colv_list, lcolv, &
208 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
209 15773 : imass, max_sigma)
210 :
211 15773 : END SUBROUTINE shake_roll_colv_int
212 :
213 : ! **************************************************************************************************
214 : !> \brief Intramolecular subroutine
215 : !> rattle algorithm (box allowed to change) for collective variables constraints
216 : !> updates the multiplier one molecule type at a time
217 : !> \param molecule ...
218 : !> \param particle_set ...
219 : !> \param vel ...
220 : !> \param r_rattle ...
221 : !> \param dt ...
222 : !> \param irattle ...
223 : !> \param veps ...
224 : !> \param cell ...
225 : !> \param imass ...
226 : !> \param max_sigma ...
227 : !> \par History
228 : !> none
229 : !> \author Teodoro Laino [tlaino]
230 : ! **************************************************************************************************
231 5521 : SUBROUTINE rattle_roll_colv_int(molecule, particle_set, vel, r_rattle, &
232 5521 : dt, irattle, veps, cell, imass, max_sigma)
233 :
234 : TYPE(molecule_type), POINTER :: molecule
235 : TYPE(particle_type), POINTER :: particle_set(:)
236 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
237 : REAL(KIND=dp), INTENT(IN) :: r_rattle(:, :), dt
238 : INTEGER, INTENT(in) :: irattle
239 : REAL(KIND=dp), INTENT(IN) :: veps(:, :)
240 : TYPE(cell_type), POINTER :: cell
241 : REAL(KIND=dp), DIMENSION(:) :: imass
242 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
243 :
244 5521 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
245 5521 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
246 5521 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
247 : TYPE(molecule_kind_type), POINTER :: molecule_kind
248 :
249 5521 : NULLIFY (fixd_list)
250 5521 : molecule_kind => molecule%molecule_kind
251 5521 : CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
252 5521 : CALL get_molecule(molecule, lcolv=lcolv)
253 : ! Real Rattle
254 : CALL rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
255 : particle_set, vel, r_rattle, dt, irattle, veps, cell, &
256 5521 : imass, max_sigma)
257 :
258 5521 : END SUBROUTINE rattle_roll_colv_int
259 :
260 : ! **************************************************************************************************
261 : !> \brief Intermolecular subroutine
262 : !> shake_colv algorithm for collective variables constraints
263 : !> updates the multiplier one molecule type at a time
264 : !> \param gci ...
265 : !> \param particle_set ...
266 : !> \param pos ...
267 : !> \param vel ...
268 : !> \param dt ...
269 : !> \param ishake ...
270 : !> \param cell ...
271 : !> \param imass ...
272 : !> \param max_sigma ...
273 : !> \par History
274 : !> none
275 : !> \author Teodoro Laino [tlaino]
276 : ! **************************************************************************************************
277 1719 : SUBROUTINE shake_colv_ext(gci, particle_set, pos, vel, dt, ishake, &
278 1719 : cell, imass, max_sigma)
279 :
280 : TYPE(global_constraint_type), POINTER :: gci
281 : TYPE(particle_type), POINTER :: particle_set(:)
282 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
283 : REAL(kind=dp), INTENT(in) :: dt
284 : INTEGER, INTENT(IN) :: ishake
285 : TYPE(cell_type), POINTER :: cell
286 : REAL(KIND=dp), DIMENSION(:) :: imass
287 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
288 :
289 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
290 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
291 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
292 :
293 1719 : colv_list => gci%colv_list
294 1719 : fixd_list => gci%fixd_list
295 1719 : lcolv => gci%lcolv
296 : ! Real Shake
297 : CALL shake_colv_low(fixd_list, colv_list, lcolv, &
298 1719 : particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
299 :
300 1719 : END SUBROUTINE shake_colv_ext
301 :
302 : ! **************************************************************************************************
303 : !> \brief Intermolecular subroutine for updating the TARGET value for collective
304 : !> constraints
305 : !> \param gci ...
306 : !> \param dt ...
307 : !> \param motion_section ...
308 : !> \author Teodoro Laino [tlaino]
309 : ! **************************************************************************************************
310 1094 : SUBROUTINE shake_update_colv_ext(gci, dt, motion_section)
311 :
312 : TYPE(global_constraint_type), POINTER :: gci
313 : REAL(kind=dp), INTENT(in) :: dt
314 : TYPE(section_vals_type), POINTER :: motion_section
315 :
316 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
317 :
318 1094 : colv_list => gci%colv_list
319 : ! Real update of the Shake target
320 1094 : CALL shake_update_colv_low(colv_list, dt, motion_section)
321 :
322 1094 : END SUBROUTINE shake_update_colv_ext
323 :
324 : ! **************************************************************************************************
325 : !> \brief Intermolecular subroutine
326 : !> rattle algorithm for collective variables constraints
327 : !> updates the multiplier one molecule type at a time
328 : !> \param gci ...
329 : !> \param particle_set ...
330 : !> \param vel ...
331 : !> \param dt ...
332 : !> \param irattle ...
333 : !> \param cell ...
334 : !> \param imass ...
335 : !> \param max_sigma ...
336 : !> \par History
337 : !> none
338 : !> \author Teodoro Laino [tlaino]
339 : ! **************************************************************************************************
340 1442 : SUBROUTINE rattle_colv_ext(gci, particle_set, vel, dt, irattle, &
341 1442 : cell, imass, max_sigma)
342 :
343 : TYPE(global_constraint_type), POINTER :: gci
344 : TYPE(particle_type), POINTER :: particle_set(:)
345 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
346 : REAL(kind=dp), INTENT(in) :: dt
347 : INTEGER, INTENT(IN) :: irattle
348 : TYPE(cell_type), POINTER :: cell
349 : REAL(KIND=dp), DIMENSION(:) :: imass
350 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
351 :
352 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
353 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
354 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
355 :
356 1442 : colv_list => gci%colv_list
357 1442 : fixd_list => gci%fixd_list
358 1442 : lcolv => gci%lcolv
359 : ! Real Rattle
360 : CALL rattle_colv_low(fixd_list, colv_list, lcolv, &
361 1442 : particle_set, vel, dt, irattle, cell, imass, max_sigma)
362 :
363 1442 : END SUBROUTINE rattle_colv_ext
364 :
365 : ! **************************************************************************************************
366 : !> \brief Intermolecular subroutine
367 : !> shake algorithm (box allowed to change) for collective variables constraints
368 : !> updates the multiplier one molecule type at a time
369 : !> \param gci ...
370 : !> \param particle_set ...
371 : !> \param pos ...
372 : !> \param vel ...
373 : !> \param r_shake ...
374 : !> \param v_shake ...
375 : !> \param dt ...
376 : !> \param ishake ...
377 : !> \param cell ...
378 : !> \param imass ...
379 : !> \param max_sigma ...
380 : !> \par History
381 : !> none
382 : !> \author Teodoro Laino [tlaino]
383 : ! **************************************************************************************************
384 0 : SUBROUTINE shake_roll_colv_ext(gci, particle_set, pos, vel, r_shake, v_shake, &
385 0 : dt, ishake, cell, imass, max_sigma)
386 :
387 : TYPE(global_constraint_type), POINTER :: gci
388 : TYPE(particle_type), POINTER :: particle_set(:)
389 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
390 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
391 : REAL(kind=dp), INTENT(in) :: dt
392 : INTEGER, INTENT(in) :: ishake
393 : TYPE(cell_type), POINTER :: cell
394 : REAL(KIND=dp), DIMENSION(:) :: imass
395 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
396 :
397 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
398 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
399 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
400 :
401 0 : colv_list => gci%colv_list
402 0 : fixd_list => gci%fixd_list
403 0 : lcolv => gci%lcolv
404 : ! Real Shake
405 : CALL shake_roll_colv_low(fixd_list, colv_list, lcolv, &
406 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
407 0 : imass, max_sigma)
408 :
409 0 : END SUBROUTINE shake_roll_colv_ext
410 :
411 : ! **************************************************************************************************
412 : !> \brief Intermolecular subroutine
413 : !> rattle algorithm (box allowed to change) for collective variables constraints
414 : !> updates the multiplier one molecule type at a time
415 : !> \param gci ...
416 : !> \param particle_set ...
417 : !> \param vel ...
418 : !> \param r_rattle ...
419 : !> \param dt ...
420 : !> \param irattle ...
421 : !> \param veps ...
422 : !> \param cell ...
423 : !> \param imass ...
424 : !> \param max_sigma ...
425 : !> \par History
426 : !> none
427 : !> \author Teodoro Laino [tlaino]
428 : ! **************************************************************************************************
429 0 : SUBROUTINE rattle_roll_colv_ext(gci, particle_set, vel, r_rattle, &
430 0 : dt, irattle, veps, cell, imass, max_sigma)
431 :
432 : TYPE(global_constraint_type), POINTER :: gci
433 : TYPE(particle_type), POINTER :: particle_set(:)
434 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
435 : REAL(KIND=dp), INTENT(IN) :: r_rattle(:, :), dt
436 : INTEGER, INTENT(in) :: irattle
437 : REAL(KIND=dp), INTENT(IN) :: veps(:, :)
438 : TYPE(cell_type), POINTER :: cell
439 : REAL(KIND=dp), DIMENSION(:) :: imass
440 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
441 :
442 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
443 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
444 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
445 :
446 0 : colv_list => gci%colv_list
447 0 : fixd_list => gci%fixd_list
448 0 : lcolv => gci%lcolv
449 : ! Real Rattle
450 : CALL rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
451 : particle_set, vel, r_rattle, dt, irattle, veps, cell, &
452 0 : imass, max_sigma)
453 :
454 0 : END SUBROUTINE rattle_roll_colv_ext
455 :
456 : ! **************************************************************************************************
457 : !> \brief Real Shake subroutine - Low Level
458 : !> shake_colv algorithm for collective variables constraints
459 : !> updates the multiplier one molecule type at a time
460 : !> \param fixd_list ...
461 : !> \param colv_list ...
462 : !> \param lcolv ...
463 : !> \param particle_set ...
464 : !> \param pos ...
465 : !> \param vel ...
466 : !> \param dt ...
467 : !> \param ishake ...
468 : !> \param cell ...
469 : !> \param imass ...
470 : !> \param max_sigma ...
471 : !> \par History
472 : !> none
473 : !> \author Teodoro Laino [tlaino]
474 : ! **************************************************************************************************
475 20417 : SUBROUTINE shake_colv_low(fixd_list, colv_list, lcolv, &
476 20417 : particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
477 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
478 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
479 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
480 : TYPE(particle_type), POINTER :: particle_set(:)
481 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
482 : REAL(kind=dp), INTENT(in) :: dt
483 : INTEGER, INTENT(IN) :: ishake
484 : TYPE(cell_type), POINTER :: cell
485 : REAL(KIND=dp), DIMENSION(:) :: imass
486 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
487 :
488 : INTEGER :: iconst
489 : REAL(KIND=dp) :: del_lam, dtby2, dtsqby2, fdotf_sum
490 :
491 20417 : dtsqby2 = dt*dt*.5_dp
492 20417 : dtby2 = dt*.5_dp
493 20417 : IF (ishake == 1) THEN
494 9800 : DO iconst = 1, SIZE(colv_list)
495 6775 : IF (colv_list(iconst)%restraint%active) CYCLE
496 : ! Update positions
497 : CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
498 : lambda=lcolv(iconst)%lambda, &
499 5855 : imass=imass)
500 : ! Update velocities
501 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
502 : lambda=lcolv(iconst)%lambda, &
503 9800 : imass=imass)
504 : END DO
505 : ELSE
506 107106 : DO iconst = 1, SIZE(colv_list)
507 89714 : IF (colv_list(iconst)%restraint%active) CYCLE
508 : ! Update colvar
509 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
510 89680 : pos=pos, fixd_list=fixd_list)
511 : lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
512 89680 : colv_list(iconst)%expected_value)
513 : fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar, &
514 89680 : lcolv(iconst)%colvar_old, imass=imass)
515 89680 : del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*dt*fdotf_sum)
516 89680 : lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
517 :
518 : ! Update positions
519 : CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
520 : lambda=del_lam, &
521 89680 : imass=imass)
522 : ! Update velocities
523 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
524 : lambda=del_lam, &
525 107106 : imass=imass)
526 : END DO
527 : END IF
528 : ! computing the constraint and value of tolerance
529 116906 : DO iconst = 1, SIZE(colv_list)
530 96489 : IF (colv_list(iconst)%restraint%active) CYCLE
531 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
532 95535 : pos=pos, fixd_list=fixd_list)
533 : lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
534 95535 : colv_list(iconst)%expected_value)
535 116906 : max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
536 : END DO
537 20417 : END SUBROUTINE shake_colv_low
538 :
539 : ! **************************************************************************************************
540 : !> \brief Real Shake subroutine - Low Level - for updating the TARGET value
541 : !> \param colv_list ...
542 : !> \param dt ...
543 : !> \param motion_section ...
544 : !> \date 02.2008
545 : !> \author Teodoro Laino [tlaino] - University of Zurich
546 : ! **************************************************************************************************
547 7370 : SUBROUTINE shake_update_colv_low(colv_list, dt, motion_section)
548 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
549 : REAL(kind=dp), INTENT(in) :: dt
550 : TYPE(section_vals_type), POINTER :: motion_section
551 :
552 : INTEGER :: iconst, irep, n_rep
553 : LOGICAL :: do_update_colvar, explicit
554 : REAL(KIND=dp) :: clv_target, limit, new_clv_target, value
555 : TYPE(section_vals_type), POINTER :: collective_sections
556 :
557 : ! Update globally for restart
558 :
559 3685 : collective_sections => section_vals_get_subs_vals(motion_section, "CONSTRAINT%COLLECTIVE")
560 3685 : CALL section_vals_get(collective_sections, n_repetition=n_rep)
561 3685 : IF (n_rep /= 0) THEN
562 10798 : DO irep = 1, n_rep
563 : CALL section_vals_val_get(collective_sections, "TARGET_GROWTH", r_val=value, &
564 7813 : i_rep_section=irep)
565 10798 : IF (value /= 0.0_dp) THEN
566 : CALL section_vals_val_get(collective_sections, "TARGET", r_val=clv_target, &
567 200 : i_rep_section=irep)
568 200 : new_clv_target = clv_target + value*dt
569 : ! Check limits..
570 : CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", explicit=explicit, &
571 200 : i_rep_section=irep)
572 200 : do_update_colvar = .TRUE.
573 200 : IF (explicit) THEN
574 : CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", r_val=limit, &
575 100 : i_rep_section=irep)
576 100 : IF (value > 0.0_dp) THEN
577 50 : IF (clv_target == limit) THEN
578 : do_update_colvar = .FALSE.
579 22 : ELSE IF (new_clv_target >= limit) THEN
580 1 : new_clv_target = limit
581 : END IF
582 : ELSE
583 50 : IF (clv_target == limit) THEN
584 : do_update_colvar = .FALSE.
585 39 : ELSE IF (new_clv_target <= limit) THEN
586 1 : new_clv_target = limit
587 : END IF
588 : END IF
589 : END IF
590 : IF (do_update_colvar) THEN
591 : CALL section_vals_val_set(collective_sections, "TARGET", r_val=new_clv_target, &
592 161 : i_rep_section=irep)
593 : END IF
594 : END IF
595 : END DO
596 : END IF
597 :
598 : ! Update locally the value to each processor
599 13318 : DO iconst = 1, SIZE(colv_list)
600 : ! Update local to each processor
601 9633 : IF (colv_list(iconst)%expected_value_growth_speed == 0.0_dp) CYCLE
602 : CALL section_vals_val_get(collective_sections, "TARGET", &
603 : r_val=colv_list(iconst)%expected_value, &
604 13318 : i_rep_section=colv_list(iconst)%inp_seq_num)
605 : END DO
606 :
607 3685 : END SUBROUTINE shake_update_colv_low
608 :
609 : ! **************************************************************************************************
610 : !> \brief Real Rattle - Low Level
611 : !> rattle algorithm for collective variables constraints
612 : !> updates the multiplier one molecule type at a time
613 : !> \param fixd_list ...
614 : !> \param colv_list ...
615 : !> \param lcolv ...
616 : !> \param particle_set ...
617 : !> \param vel ...
618 : !> \param dt ...
619 : !> \param irattle ...
620 : !> \param cell ...
621 : !> \param imass ...
622 : !> \param max_sigma ...
623 : !> \par History
624 : !> none
625 : !> \author Teodoro Laino [tlaino]
626 : ! **************************************************************************************************
627 9275 : SUBROUTINE rattle_colv_low(fixd_list, colv_list, lcolv, &
628 9275 : particle_set, vel, dt, irattle, cell, imass, max_sigma)
629 :
630 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
631 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
632 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
633 : TYPE(particle_type), POINTER :: particle_set(:)
634 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
635 : REAL(kind=dp), INTENT(in) :: dt
636 : INTEGER, INTENT(IN) :: irattle
637 : TYPE(cell_type), POINTER :: cell
638 : REAL(KIND=dp), DIMENSION(:) :: imass
639 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
640 :
641 : INTEGER :: iconst
642 : REAL(KIND=dp) :: del_lam, dtby2, fdotf_sum
643 :
644 9275 : dtby2 = dt*.5_dp
645 9275 : IF (irattle == 1) THEN
646 9812 : DO iconst = 1, SIZE(colv_list)
647 6781 : IF (colv_list(iconst)%restraint%active) CYCLE
648 : ! Update colvar_old
649 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
650 5861 : particles=particle_set, fixd_list=fixd_list)
651 : ! Update velocities
652 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
653 : lambda=lcolv(iconst)%lambda, &
654 9812 : imass=imass)
655 : END DO
656 : ELSE
657 35597 : DO iconst = 1, SIZE(colv_list)
658 29353 : IF (colv_list(iconst)%restraint%active) CYCLE
659 29333 : lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel)
660 : fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar_old, &
661 29333 : lcolv(iconst)%colvar_old, imass=imass)
662 29333 : del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*fdotf_sum)
663 29333 : lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
664 :
665 : ! Update velocities
666 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
667 : lambda=del_lam, &
668 35597 : imass=imass)
669 : END DO
670 : END IF
671 :
672 45409 : DO iconst = 1, SIZE(colv_list)
673 36134 : IF (colv_list(iconst)%restraint%active) CYCLE
674 35194 : lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel)
675 45409 : max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
676 : END DO
677 :
678 9275 : END SUBROUTINE rattle_colv_low
679 :
680 : ! **************************************************************************************************
681 : !> \brief Real shake_roll - Low Level
682 : !> shake algorithm (box allowed to change) for collective variables constraints
683 : !> updates the multiplier one molecule type at a time
684 : !> \param fixd_list ...
685 : !> \param colv_list ...
686 : !> \param lcolv ...
687 : !> \param particle_set ...
688 : !> \param pos ...
689 : !> \param vel ...
690 : !> \param r_shake ...
691 : !> \param v_shake ...
692 : !> \param dt ...
693 : !> \param ishake ...
694 : !> \param cell ...
695 : !> \param imass ...
696 : !> \param max_sigma ...
697 : !> \par History
698 : !> none
699 : !> \author Teodoro Laino [tlaino]
700 : ! **************************************************************************************************
701 15773 : SUBROUTINE shake_roll_colv_low(fixd_list, colv_list, lcolv, &
702 31546 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
703 15773 : imass, max_sigma)
704 :
705 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
706 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
707 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
708 : TYPE(particle_type), POINTER :: particle_set(:)
709 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
710 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
711 : REAL(kind=dp), INTENT(in) :: dt
712 : INTEGER, INTENT(in) :: ishake
713 : TYPE(cell_type), POINTER :: cell
714 : REAL(KIND=dp), DIMENSION(:) :: imass
715 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
716 :
717 : INTEGER :: iconst
718 : REAL(KIND=dp) :: del_lam, dtby2, dtsqby2, fdotf_sum
719 :
720 15773 : dtsqby2 = dt*dt*.5_dp
721 15773 : dtby2 = dt*.5_dp
722 15773 : IF (ishake == 1) THEN
723 8664 : DO iconst = 1, SIZE(colv_list)
724 7312 : IF (colv_list(iconst)%restraint%active) CYCLE
725 : ! Update positions
726 : CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
727 : lambda=lcolv(iconst)%lambda, &
728 7312 : roll=.TRUE., rmat=r_shake, imass=imass)
729 : ! Update velocities
730 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
731 : lambda=lcolv(iconst)%lambda, &
732 8664 : roll=.TRUE., rmat=v_shake, imass=imass)
733 : END DO
734 : ELSE
735 100275 : DO iconst = 1, SIZE(colv_list)
736 85854 : IF (colv_list(iconst)%restraint%active) CYCLE
737 : ! Update colvar
738 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
739 85854 : pos=pos, fixd_list=fixd_list)
740 : lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
741 85854 : colv_list(iconst)%expected_value)
742 : fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar, &
743 : lcolv(iconst)%colvar_old, roll=.TRUE., rmat=r_shake, &
744 85854 : imass=imass)
745 85854 : del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*dt*fdotf_sum)
746 85854 : lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
747 :
748 : ! Update positions
749 : CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
750 : lambda=del_lam, &
751 85854 : roll=.TRUE., rmat=r_shake, imass=imass)
752 : ! Update velocities
753 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
754 : lambda=del_lam, &
755 100275 : roll=.TRUE., rmat=v_shake, imass=imass)
756 : END DO
757 : END IF
758 : ! computing the constraint and value of tolerance
759 108939 : DO iconst = 1, SIZE(colv_list)
760 93166 : IF (colv_list(iconst)%restraint%active) CYCLE
761 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
762 93166 : pos=pos, fixd_list=fixd_list)
763 : lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
764 93166 : colv_list(iconst)%expected_value)
765 108939 : max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
766 : END DO
767 :
768 15773 : END SUBROUTINE shake_roll_colv_low
769 :
770 : ! **************************************************************************************************
771 : !> \brief Real Rattle_roll - Low Level
772 : !> rattle algorithm (box allowed to change) for collective variables constraints
773 : !> updates the multiplier one molecule type at a time
774 : !> \param fixd_list ...
775 : !> \param colv_list ...
776 : !> \param lcolv ...
777 : !> \param particle_set ...
778 : !> \param vel ...
779 : !> \param r_rattle ...
780 : !> \param dt ...
781 : !> \param irattle ...
782 : !> \param veps ...
783 : !> \param cell ...
784 : !> \param imass ...
785 : !> \param max_sigma ...
786 : !> \par History
787 : !> none
788 : !> \author Teodoro Laino [tlaino]
789 : ! **************************************************************************************************
790 5521 : SUBROUTINE rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
791 11042 : particle_set, vel, r_rattle, dt, irattle, veps, cell, &
792 5521 : imass, max_sigma)
793 :
794 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
795 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
796 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
797 : TYPE(particle_type), POINTER :: particle_set(:)
798 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
799 : REAL(KIND=dp), INTENT(IN) :: r_rattle(:, :), dt
800 : INTEGER, INTENT(in) :: irattle
801 : REAL(KIND=dp), INTENT(IN) :: veps(:, :)
802 : TYPE(cell_type), POINTER :: cell
803 : REAL(KIND=dp), DIMENSION(:) :: imass
804 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
805 :
806 : INTEGER :: iconst
807 : REAL(KIND=dp) :: del_lam, dtby2, fdotf_sum
808 :
809 5521 : dtby2 = dt*.5_dp
810 5521 : IF (irattle == 1) THEN
811 7957 : DO iconst = 1, SIZE(colv_list)
812 6670 : IF (colv_list(iconst)%restraint%active) CYCLE
813 : ! Update colvar_old
814 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
815 6670 : particles=particle_set, fixd_list=fixd_list)
816 : ! Update velocities
817 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
818 : lambda=lcolv(iconst)%lambda, &
819 7957 : imass=imass)
820 : END DO
821 : ELSE
822 29110 : DO iconst = 1, SIZE(colv_list)
823 24876 : IF (colv_list(iconst)%restraint%active) CYCLE
824 : lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel, &
825 24876 : roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set)
826 : fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar_old, &
827 : lcolv(iconst)%colvar_old, roll=.TRUE., &
828 24876 : rmat=r_rattle, imass=imass)
829 24876 : del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*fdotf_sum)
830 24876 : lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
831 : ! Update velocities
832 : CALL update_con_colv(vel, dtby2, lcolv(iconst), &
833 : lambda=del_lam, &
834 29110 : roll=.TRUE., rmat=r_rattle, imass=imass)
835 : END DO
836 : END IF
837 : ! computing the constraint and value of the tolerance
838 37067 : DO iconst = 1, SIZE(colv_list)
839 31546 : IF (colv_list(iconst)%restraint%active) CYCLE
840 : lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel, &
841 31546 : roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set)
842 37067 : max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
843 : END DO
844 :
845 5521 : END SUBROUTINE rattle_roll_colv_low
846 :
847 : ! **************************************************************************************************
848 : !> \brief Update position/velocities
849 : !> \param wrk ...
850 : !> \param fac ...
851 : !> \param lcolv ...
852 : !> \param lambda ...
853 : !> \param roll ...
854 : !> \param rmat ...
855 : !> \param imass ...
856 : !> \par History
857 : !> Teodoro Laino [teo] created 04.2006
858 : !> \author Teodoro Laino [tlaino]
859 : ! **************************************************************************************************
860 444142 : SUBROUTINE update_con_colv(wrk, fac, lcolv, lambda, roll, rmat, imass)
861 : REAL(KIND=dp), INTENT(INOUT) :: wrk(:, :)
862 : REAL(KIND=dp), INTENT(IN) :: fac
863 : TYPE(local_colvar_constraint_type), INTENT(IN) :: lcolv
864 : REAL(KIND=dp), INTENT(IN) :: lambda
865 : LOGICAL, INTENT(in), OPTIONAL :: roll
866 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
867 : OPTIONAL :: rmat
868 : REAL(KIND=dp), DIMENSION(:) :: imass
869 :
870 : INTEGER :: iatm, ind
871 : LOGICAL :: my_roll
872 : REAL(KIND=dp), DIMENSION(3) :: f_roll
873 :
874 444142 : my_roll = .FALSE.
875 444142 : IF (PRESENT(roll)) THEN
876 211208 : my_roll = roll
877 211208 : IF (my_roll) THEN
878 211208 : CPASSERT(PRESENT(rmat))
879 : END IF
880 : END IF
881 1344082 : DO iatm = 1, SIZE(lcolv%colvar_old%i_atom)
882 899940 : ind = lcolv%colvar_old%i_atom(iatm)
883 : !
884 899940 : IF (my_roll) THEN
885 : ! If ROLL rotate forces
886 5520684 : f_roll = MATMUL(rmat, lcolv%colvar_old%dsdr(:, iatm))
887 : ELSE
888 1901088 : f_roll = lcolv%colvar_old%dsdr(:, iatm)
889 : END IF
890 4043902 : wrk(:, ind) = wrk(:, ind) - imass(ind)*fac*lambda*f_roll
891 : END DO
892 444142 : END SUBROUTINE update_con_colv
893 :
894 : ! **************************************************************************************************
895 : !> \brief Evaluates the Jacobian of the collective variables constraints
896 : !> \param colvar ...
897 : !> \param colvar_old ...
898 : !> \param roll ...
899 : !> \param rmat ...
900 : !> \param imass ...
901 : !> \return ...
902 : !> \par History
903 : !> Teodoro Laino [teo] created 04.2006
904 : !> \author Teodoro Laino [tlaino]
905 : ! **************************************************************************************************
906 229743 : FUNCTION eval_Jac_colvar(colvar, colvar_old, roll, rmat, imass) RESULT(res)
907 : TYPE(colvar_type), POINTER :: colvar, colvar_old
908 : LOGICAL, INTENT(IN), OPTIONAL :: roll
909 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
910 : OPTIONAL :: rmat
911 : REAL(KIND=dp), DIMENSION(:) :: imass
912 : REAL(KIND=dp) :: res
913 :
914 : INTEGER :: i, iatom
915 : LOGICAL :: my_roll
916 : REAL(KIND=dp), DIMENSION(3) :: tmp1, tmp2, tmp3
917 :
918 229743 : my_roll = .FALSE.
919 229743 : IF (PRESENT(roll)) THEN
920 110730 : my_roll = roll
921 110730 : IF (my_roll) THEN
922 110730 : CPASSERT(PRESENT(rmat))
923 : END IF
924 : END IF
925 :
926 229743 : res = 0.0_dp
927 : IF (my_roll) THEN
928 332989 : DO i = 1, SIZE(colvar%i_atom)
929 222259 : iatom = colvar%i_atom(i)
930 889036 : tmp1 = colvar%dsdr(1:3, i)
931 889036 : tmp3 = colvar_old%dsdr(1:3, i)
932 2889367 : tmp2 = MATMUL(rmat, tmp3)
933 999766 : res = res + DOT_PRODUCT(tmp1, tmp2)*imass(iatom)
934 : END DO
935 : ELSE
936 360221 : DO i = 1, SIZE(colvar%i_atom)
937 241208 : iatom = colvar%i_atom(i)
938 964832 : tmp1 = colvar%dsdr(1:3, i)
939 964832 : tmp2 = colvar_old%dsdr(1:3, i)
940 1083845 : res = res + DOT_PRODUCT(tmp1, tmp2)*imass(iatom)
941 : END DO
942 : END IF
943 :
944 229743 : END FUNCTION eval_Jac_colvar
945 :
946 : ! **************************************************************************************************
947 : !> \brief Evaluates the constraint for the rattle scheme
948 : !> \param colvar ...
949 : !> \param vel ...
950 : !> \param roll ...
951 : !> \param veps ...
952 : !> \param rmat ...
953 : !> \param particles ...
954 : !> \return ...
955 : !> \par History
956 : !> Teodoro Laino [teo] created 04.2006
957 : !> \author Teodoro Laino [tlaino]
958 : ! **************************************************************************************************
959 120949 : FUNCTION rattle_con_eval(colvar, vel, roll, veps, rmat, particles) RESULT(res)
960 : TYPE(colvar_type), POINTER :: colvar
961 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
962 : LOGICAL, INTENT(IN), OPTIONAL :: roll
963 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
964 : OPTIONAL :: veps, rmat
965 : TYPE(particle_type), DIMENSION(:), OPTIONAL, &
966 : POINTER :: particles
967 : REAL(KIND=dp) :: res
968 :
969 : INTEGER :: iatm, ind
970 : LOGICAL :: my_roll
971 : REAL(KIND=dp), DIMENSION(3) :: f_roll, pos, v_roll
972 :
973 120949 : my_roll = .FALSE.
974 120949 : IF (PRESENT(roll)) THEN
975 56422 : my_roll = roll
976 56422 : IF (my_roll) THEN
977 56422 : CPASSERT(PRESENT(rmat))
978 56422 : CPASSERT(PRESENT(veps))
979 56422 : CPASSERT(PRESENT(particles))
980 : END IF
981 : END IF
982 120949 : res = 0.0_dp
983 368027 : DO iatm = 1, SIZE(colvar%i_atom)
984 247078 : ind = colvar%i_atom(iatm)
985 247078 : IF (my_roll) THEN
986 456848 : pos = particles(ind)%r
987 1484756 : f_roll = MATMUL(rmat, colvar%dsdr(:, iatm))
988 1827392 : v_roll = vel(:, ind) + MATMUL(veps, pos)
989 : ELSE
990 531464 : f_roll = colvar%dsdr(:, iatm)
991 531464 : v_roll = vel(:, ind)
992 : END IF
993 1109261 : res = res + DOT_PRODUCT(f_roll, v_roll)
994 : END DO
995 :
996 120949 : END FUNCTION rattle_con_eval
997 :
998 : END MODULE constraint_clv
|
