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