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 : !> \par History
10 : !> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints
11 : ! **************************************************************************************************
12 : MODULE constraint
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE cell_types, ONLY: cell_type
16 : USE colvar_types, ONLY: colvar_counters
17 : USE constraint_3x3, ONLY: rattle_3x3_ext,&
18 : rattle_3x3_int,&
19 : rattle_roll_3x3_ext,&
20 : rattle_roll_3x3_int,&
21 : shake_3x3_ext,&
22 : shake_3x3_int,&
23 : shake_roll_3x3_ext,&
24 : shake_roll_3x3_int
25 : USE constraint_4x6, ONLY: rattle_4x6_ext,&
26 : rattle_4x6_int,&
27 : rattle_roll_4x6_ext,&
28 : rattle_roll_4x6_int,&
29 : shake_4x6_ext,&
30 : shake_4x6_int,&
31 : shake_roll_4x6_ext,&
32 : shake_roll_4x6_int
33 : USE constraint_clv, ONLY: &
34 : check_colv_int, rattle_colv_ext, rattle_colv_int, rattle_roll_colv_ext, &
35 : rattle_roll_colv_int, shake_colv_ext, shake_colv_int, shake_roll_colv_ext, &
36 : shake_roll_colv_int, shake_update_colv_ext, shake_update_colv_low, &
37 : shake_update_colv_targets
38 : USE constraint_util, ONLY: check_tol,&
39 : get_roll_matrix,&
40 : restore_temporary_set,&
41 : update_temporary_set
42 : USE constraint_vsite, ONLY: shake_vsite_ext,&
43 : shake_vsite_int
44 : USE cp_log_handling, ONLY: cp_to_string
45 : USE distribution_1d_types, ONLY: distribution_1d_type
46 : USE input_constants, ONLY: npt_f_ensemble,&
47 : npt_i_ensemble,&
48 : npt_ia_ensemble
49 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
50 : section_vals_type
51 : USE kinds, ONLY: default_string_length,&
52 : dp
53 : USE memory_utilities, ONLY: reallocate
54 : USE message_passing, ONLY: mp_comm_type,&
55 : mp_para_env_type
56 : USE molecule_kind_types, ONLY: get_molecule_kind,&
57 : get_molecule_kind_set,&
58 : molecule_kind_type
59 : USE molecule_types, ONLY: global_constraint_type,&
60 : molecule_type
61 : USE particle_types, ONLY: particle_type
62 : USE simpar_types, ONLY: simpar_type
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 :
67 : PRIVATE
68 : PUBLIC :: shake_control, &
69 : rattle_control, &
70 : shake_roll_control, &
71 : rattle_roll_control, &
72 : shake_update_targets
73 :
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint'
75 : INTEGER, PARAMETER, PRIVATE :: Max_Shake_Iter = 1000
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief ...
81 : !> \param gci ...
82 : !> \param local_molecules ...
83 : !> \param molecule_set ...
84 : !> \param molecule_kind_set ...
85 : !> \param particle_set ...
86 : !> \param pos ...
87 : !> \param vel ...
88 : !> \param dt ...
89 : !> \param shake_tol ...
90 : !> \param log_unit ...
91 : !> \param lagrange_mult ...
92 : !> \param dump_lm ...
93 : !> \param cell ...
94 : !> \param group ...
95 : !> \param local_particles ...
96 : !> \par History
97 : !> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints
98 : ! **************************************************************************************************
99 16394 : SUBROUTINE shake_control(gci, local_molecules, molecule_set, molecule_kind_set, &
100 16394 : particle_set, pos, vel, dt, shake_tol, log_unit, lagrange_mult, dump_lm, &
101 : cell, group, local_particles)
102 :
103 : TYPE(global_constraint_type), POINTER :: gci
104 : TYPE(distribution_1d_type), POINTER :: local_molecules
105 : TYPE(molecule_type), POINTER :: molecule_set(:)
106 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
107 : TYPE(particle_type), POINTER :: particle_set(:)
108 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
109 : REAL(kind=dp), INTENT(in) :: dt, shake_tol
110 : INTEGER, INTENT(in) :: log_unit, lagrange_mult
111 : LOGICAL, INTENT(IN) :: dump_lm
112 : TYPE(cell_type), POINTER :: cell
113 :
114 : CLASS(mp_comm_type), INTENT(in) :: group
115 : TYPE(distribution_1d_type), POINTER :: local_particles
116 :
117 : CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_control'
118 :
119 : INTEGER :: handle, i, ikind, imol, ishake_ext, &
120 : ishake_int, k, n3x3con, n4x6con, &
121 : nconstraint, nkind, nmol_per_kind, &
122 : nvsitecon
123 : LOGICAL :: do_ext_constraint
124 : REAL(KIND=dp) :: int_max_sigma, mass, max_sigma
125 32788 : REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass
126 : TYPE(atomic_kind_type), POINTER :: atomic_kind
127 : TYPE(colvar_counters) :: ncolv
128 : TYPE(molecule_kind_type), POINTER :: molecule_kind
129 : TYPE(molecule_type), POINTER :: molecule
130 :
131 16394 : CALL timeset(routineN, handle)
132 16394 : nkind = SIZE(molecule_kind_set)
133 1595302 : DO k = 1, SIZE(pos, 2)
134 1578908 : atomic_kind => particle_set(k)%atomic_kind
135 1578908 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
136 1595302 : imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp))
137 : END DO
138 16394 : do_ext_constraint = (gci%ntot /= 0)
139 16394 : ishake_ext = 0
140 16394 : max_sigma = -1.0E+10_dp
141 33491 : Shake_Inter_Loop: DO WHILE ((ABS(max_sigma) >= shake_tol) .AND. (ishake_ext <= Max_Shake_Iter))
142 17097 : max_sigma = 0.0_dp
143 17097 : ishake_ext = ishake_ext + 1
144 : ! Intramolecular Constraints
145 84962 : MOL: DO ikind = 1, nkind
146 67865 : nmol_per_kind = local_molecules%n_el(ikind)
147 358381 : DO imol = 1, nmol_per_kind
148 273419 : i = local_molecules%list(ikind)%array(imol)
149 273419 : molecule => molecule_set(i)
150 273419 : molecule_kind => molecule%molecule_kind
151 : CALL get_molecule_kind(molecule_kind, ncolv=ncolv, &
152 : ng3x3=n3x3con, ng4x6=n4x6con, &
153 273419 : nconstraint=nconstraint, nvsite=nvsitecon)
154 273419 : IF (nconstraint == 0) CYCLE
155 146095 : ishake_int = MERGE(0, 1, ishake_ext == 1 .OR. .NOT. gci%joint_colvar_checks)
156 146095 : int_max_sigma = -1.0E+10_dp
157 447479 : Shake_Intra_Loop: DO WHILE ((ABS(int_max_sigma) >= shake_tol) .AND. (ishake_int <= Max_Shake_Iter))
158 301384 : int_max_sigma = 0.0_dp
159 301384 : ishake_int = ishake_int + 1
160 : ! 3x3
161 301384 : IF (n3x3con /= 0) THEN
162 : CALL shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake_int, &
163 280301 : int_max_sigma)
164 : END IF
165 : ! 4x6
166 301384 : IF (n4x6con /= 0) THEN
167 : CALL shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake_int, &
168 2466 : int_max_sigma)
169 : END IF
170 : ! Collective Variables
171 447479 : IF (ncolv%ntot /= 0) THEN
172 : CALL shake_colv_int(molecule, particle_set, pos, vel, dt, ishake_int, &
173 18617 : cell, imass, int_max_sigma)
174 : END IF
175 : END DO Shake_Intra_Loop
176 146095 : max_sigma = MAX(max_sigma, int_max_sigma)
177 146095 : CALL shake_int_info(log_unit, i, ishake_int, max_sigma)
178 : ! Virtual Site
179 487379 : IF (nvsitecon /= 0) THEN
180 838 : CALL shake_vsite_int(molecule, pos)
181 : END IF
182 : END DO
183 : END DO MOL
184 : ! Intermolecular constraints
185 17097 : IF (do_ext_constraint) THEN
186 2011 : CALL update_temporary_set(group, pos=pos, vel=vel)
187 : ! 3x3
188 2011 : IF (gci%ng3x3 /= 0) THEN
189 : CALL shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
190 76 : max_sigma)
191 : END IF
192 : ! 4x6
193 2011 : IF (gci%ng4x6 /= 0) THEN
194 : CALL shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
195 48 : max_sigma)
196 : END IF
197 : ! Collective Variables
198 2011 : IF (gci%ncolv%ntot /= 0) THEN
199 : CALL shake_colv_ext(gci, particle_set, pos, vel, dt, ishake_ext, &
200 1887 : cell, imass, max_sigma)
201 : END IF
202 : ! Virtual Site
203 2011 : IF (gci%nvsite /= 0) THEN
204 0 : CALL shake_vsite_ext(gci, pos)
205 : END IF
206 2011 : IF (gci%joint_colvar_checks) THEN
207 : CALL check_local_colvars(local_molecules, molecule_set, particle_set, cell, group, &
208 126 : max_sigma, pos=pos)
209 : END IF
210 2011 : CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel)
211 : END IF
212 17097 : CALL shake_ext_info(log_unit, ishake_ext, max_sigma)
213 : END DO Shake_Inter_Loop
214 : CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, &
215 16394 : molecule_kind_set, group, "S")
216 :
217 16394 : CALL timestop(handle)
218 16394 : END SUBROUTINE shake_control
219 :
220 : ! **************************************************************************************************
221 : !> \brief ...
222 : !> \param gci ...
223 : !> \param local_molecules ...
224 : !> \param molecule_set ...
225 : !> \param molecule_kind_set ...
226 : !> \param particle_set ...
227 : !> \param vel ...
228 : !> \param dt ...
229 : !> \param rattle_tol ...
230 : !> \param log_unit ...
231 : !> \param lagrange_mult ...
232 : !> \param dump_lm ...
233 : !> \param cell ...
234 : !> \param group ...
235 : !> \param local_particles ...
236 : !> \par History
237 : !> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints
238 : ! **************************************************************************************************
239 16410 : SUBROUTINE rattle_control(gci, local_molecules, molecule_set, molecule_kind_set, &
240 16410 : particle_set, vel, dt, rattle_tol, log_unit, lagrange_mult, dump_lm, cell, group, &
241 : local_particles)
242 :
243 : TYPE(global_constraint_type), POINTER :: gci
244 : TYPE(distribution_1d_type), POINTER :: local_molecules
245 : TYPE(molecule_type), POINTER :: molecule_set(:)
246 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
247 : TYPE(particle_type), POINTER :: particle_set(:)
248 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
249 : REAL(kind=dp), INTENT(in) :: dt, rattle_tol
250 : INTEGER, INTENT(in) :: log_unit, lagrange_mult
251 : LOGICAL, INTENT(IN) :: dump_lm
252 : TYPE(cell_type), POINTER :: cell
253 :
254 : CLASS(mp_comm_type), INTENT(in) :: group
255 : TYPE(distribution_1d_type), POINTER :: local_particles
256 :
257 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_control'
258 :
259 : INTEGER :: handle, i, ikind, imol, irattle_ext, &
260 : irattle_int, k, n3x3con, n4x6con, &
261 : nconstraint, nkind, nmol_per_kind
262 : LOGICAL :: do_ext_constraint
263 : REAL(KIND=dp) :: int_max_sigma, mass, max_sigma
264 32820 : REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass
265 : TYPE(atomic_kind_type), POINTER :: atomic_kind
266 : TYPE(colvar_counters) :: ncolv
267 : TYPE(molecule_kind_type), POINTER :: molecule_kind
268 : TYPE(molecule_type), POINTER :: molecule
269 :
270 16410 : CALL timeset(routineN, handle)
271 16410 : nkind = SIZE(molecule_kind_set)
272 1595372 : DO k = 1, SIZE(vel, 2)
273 1578962 : atomic_kind => particle_set(k)%atomic_kind
274 1578962 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
275 1595372 : imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp))
276 : END DO
277 16410 : do_ext_constraint = (gci%ntot /= 0)
278 16410 : irattle_ext = 0
279 16410 : max_sigma = -1.0E+10_dp
280 33294 : Rattle_Inter_Loop: DO WHILE (ABS(max_sigma) >= rattle_tol)
281 16884 : max_sigma = 0.0_dp
282 16884 : irattle_ext = irattle_ext + 1
283 : ! Intramolecular Constraints
284 84660 : MOL: DO ikind = 1, nkind
285 67776 : nmol_per_kind = local_molecules%n_el(ikind)
286 357006 : DO imol = 1, nmol_per_kind
287 272346 : i = local_molecules%list(ikind)%array(imol)
288 272346 : molecule => molecule_set(i)
289 272346 : molecule_kind => molecule%molecule_kind
290 : CALL get_molecule_kind(molecule_kind, ncolv=ncolv, ng3x3=n3x3con, &
291 272346 : ng4x6=n4x6con, nconstraint=nconstraint)
292 272346 : IF (nconstraint == 0) CYCLE
293 146097 : irattle_int = MERGE(0, 1, irattle_ext == 1 .OR. .NOT. gci%joint_colvar_checks)
294 146097 : int_max_sigma = -1.0E+10_dp
295 298103 : Rattle_Intra_Loop: DO WHILE (ABS(int_max_sigma) >= rattle_tol)
296 152006 : int_max_sigma = 0.0_dp
297 152006 : irattle_int = irattle_int + 1
298 : ! 3x3
299 152006 : IF (n3x3con /= 0) THEN
300 143476 : CALL rattle_3x3_int(molecule, particle_set, vel, dt)
301 : END IF
302 : ! 4x6
303 152006 : IF (n4x6con /= 0) THEN
304 682 : CALL rattle_4x6_int(molecule, particle_set, vel, dt)
305 : END IF
306 : ! Collective Variables
307 298103 : IF (ncolv%ntot /= 0) THEN
308 : CALL rattle_colv_int(molecule, particle_set, vel, dt, &
309 7848 : irattle_int, cell, imass, int_max_sigma)
310 : END IF
311 : END DO Rattle_Intra_Loop
312 146097 : max_sigma = MAX(max_sigma, int_max_sigma)
313 486219 : CALL rattle_int_info(log_unit, i, irattle_int, max_sigma)
314 : END DO
315 : END DO MOL
316 : ! Intermolecular Constraints
317 16884 : IF (do_ext_constraint) THEN
318 1794 : CALL update_temporary_set(group, vel=vel)
319 : ! 3x3
320 1794 : IF (gci%ng3x3 /= 0) THEN
321 40 : CALL rattle_3x3_ext(gci, particle_set, vel, dt)
322 : END IF
323 : ! 4x6
324 1794 : IF (gci%ng4x6 /= 0) THEN
325 20 : CALL rattle_4x6_ext(gci, particle_set, vel, dt)
326 : END IF
327 : ! Collective Variables
328 1794 : IF (gci%ncolv%ntot /= 0) THEN
329 : CALL rattle_colv_ext(gci, particle_set, vel, dt, &
330 1734 : irattle_ext, cell, imass, max_sigma)
331 : END IF
332 1794 : IF (gci%joint_colvar_checks) THEN
333 : CALL check_local_colvars(local_molecules, molecule_set, particle_set, cell, group, &
334 252 : max_sigma, vel=vel)
335 : END IF
336 1794 : CALL restore_temporary_set(particle_set, local_particles, vel=vel)
337 : END IF
338 16884 : CALL rattle_ext_info(log_unit, irattle_ext, max_sigma)
339 : END DO Rattle_Inter_Loop
340 : CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, &
341 16410 : molecule_kind_set, group, "R")
342 16410 : CALL timestop(handle)
343 :
344 16410 : END SUBROUTINE rattle_control
345 :
346 : ! **************************************************************************************************
347 : !> \brief ...
348 : !> \param gci ...
349 : !> \param local_molecules ...
350 : !> \param molecule_set ...
351 : !> \param molecule_kind_set ...
352 : !> \param particle_set ...
353 : !> \param pos ...
354 : !> \param vel ...
355 : !> \param dt ...
356 : !> \param simpar ...
357 : !> \param roll_tol ...
358 : !> \param iroll ...
359 : !> \param vector_r ...
360 : !> \param vector_v ...
361 : !> \param group ...
362 : !> \param u ...
363 : !> \param cell ...
364 : !> \param local_particles ...
365 : !> \par History
366 : !> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints
367 : ! **************************************************************************************************
368 1826 : SUBROUTINE shake_roll_control(gci, local_molecules, molecule_set, &
369 3652 : molecule_kind_set, particle_set, pos, vel, dt, simpar, roll_tol, iroll, &
370 1826 : vector_r, vector_v, group, u, cell, local_particles)
371 :
372 : TYPE(global_constraint_type), POINTER :: gci
373 : TYPE(distribution_1d_type), POINTER :: local_molecules
374 : TYPE(molecule_type), POINTER :: molecule_set(:)
375 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
376 : TYPE(particle_type), POINTER :: particle_set(:)
377 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
378 : REAL(KIND=dp), INTENT(IN) :: dt
379 : TYPE(simpar_type), INTENT(IN) :: simpar
380 : REAL(KIND=dp), INTENT(OUT) :: roll_tol
381 : INTEGER, INTENT(INOUT) :: iroll
382 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector_r, vector_v
383 :
384 : CLASS(mp_comm_type), INTENT(IN) :: group
385 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
386 : OPTIONAL :: u
387 : TYPE(cell_type), POINTER :: cell
388 : TYPE(distribution_1d_type), POINTER :: local_particles
389 :
390 : CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_roll_control'
391 :
392 : INTEGER :: handle, i, ikind, imol, ishake_ext, ishake_int, k, lagrange_mult, log_unit, &
393 : n3x3con, n4x6con, nconstraint, nkind, nmol_per_kind, nvsitecon
394 : LOGICAL :: do_ext_constraint, dump_lm
395 : REAL(KIND=dp) :: int_max_sigma, mass, max_sigma, shake_tol
396 : REAL(KIND=dp), DIMENSION(3, 3) :: r_shake, v_shake
397 3652 : REAL(KIND=dp), DIMENSION(SIZE(pos, 2)) :: imass
398 : TYPE(atomic_kind_type), POINTER :: atomic_kind
399 : TYPE(colvar_counters) :: ncolv
400 : TYPE(molecule_kind_type), POINTER :: molecule_kind
401 : TYPE(molecule_type), POINTER :: molecule
402 :
403 1826 : CALL timeset(routineN, handle)
404 1826 : nkind = SIZE(molecule_kind_set)
405 1826 : shake_tol = simpar%shake_tol
406 1826 : log_unit = simpar%info_constraint
407 1826 : lagrange_mult = simpar%lagrange_multipliers
408 1826 : dump_lm = simpar%dump_lm
409 : ! setting up for roll
410 1826 : IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npt_ia_ensemble) THEN
411 1806 : CALL get_roll_matrix('SHAKE', r_shake, v_shake, vector_r, vector_v)
412 20 : ELSE IF (simpar%ensemble == npt_f_ensemble) THEN
413 20 : CALL get_roll_matrix('SHAKE', r_shake, v_shake, vector_r, vector_v, u)
414 : END IF
415 713974 : DO k = 1, SIZE(pos, 2)
416 712148 : atomic_kind => particle_set(k)%atomic_kind
417 712148 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
418 713974 : imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp))
419 : END DO
420 1826 : do_ext_constraint = (gci%ntot /= 0)
421 1826 : ishake_ext = 0
422 1826 : max_sigma = -1.0E+10_dp
423 3652 : Shake_Inter_Loop: DO WHILE (ABS(max_sigma) >= shake_tol)
424 1826 : max_sigma = 0.0_dp
425 1826 : ishake_ext = ishake_ext + 1
426 : ! Intramolecular Constraints
427 3672 : MOL: DO ikind = 1, nkind
428 1846 : nmol_per_kind = local_molecules%n_el(ikind)
429 119630 : DO imol = 1, nmol_per_kind
430 115958 : i = local_molecules%list(ikind)%array(imol)
431 115958 : molecule => molecule_set(i)
432 115958 : molecule_kind => molecule%molecule_kind
433 : CALL get_molecule_kind(molecule_kind, ncolv=ncolv, &
434 : ng3x3=n3x3con, ng4x6=n4x6con, &
435 115958 : nconstraint=nconstraint, nvsite=nvsitecon)
436 115958 : IF (nconstraint == 0) CYCLE
437 115118 : ishake_int = MERGE(0, 1, ishake_ext == 1 .OR. .NOT. gci%joint_colvar_checks)
438 115118 : int_max_sigma = -1.0E+10_dp
439 261471 : Shake_Roll_Intra_Loop: DO WHILE (ABS(int_max_sigma) >= shake_tol)
440 146353 : int_max_sigma = 0.0_dp
441 146353 : ishake_int = ishake_int + 1
442 : ! 3x3
443 146353 : IF (n3x3con /= 0) THEN
444 : CALL shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, &
445 128355 : v_shake, dt, ishake_int, int_max_sigma)
446 : END IF
447 : ! 4x6
448 146353 : IF (n4x6con /= 0) THEN
449 : CALL shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, &
450 2225 : dt, ishake_int, int_max_sigma)
451 : END IF
452 : ! Collective Variables
453 261471 : IF (ncolv%ntot /= 0) THEN
454 : CALL shake_roll_colv_int(molecule, particle_set, pos, vel, r_shake, &
455 15773 : v_shake, dt, ishake_int, cell, imass, int_max_sigma)
456 : END IF
457 : END DO Shake_Roll_Intra_Loop
458 115118 : max_sigma = MAX(max_sigma, int_max_sigma)
459 115118 : CALL shake_int_info(log_unit, i, ishake_int, max_sigma)
460 : ! Virtual Site
461 232922 : IF (nvsitecon /= 0) THEN
462 0 : CPABORT("Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!")
463 : END IF
464 : END DO
465 : END DO MOL
466 : ! Intermolecular constraints
467 1826 : IF (do_ext_constraint) THEN
468 0 : CALL update_temporary_set(group, pos=pos, vel=vel)
469 : ! 3x3
470 0 : IF (gci%ng3x3 /= 0) THEN
471 : CALL shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, &
472 0 : v_shake, dt, ishake_ext, max_sigma)
473 : END IF
474 : ! 4x6
475 0 : IF (gci%ng4x6 /= 0) THEN
476 : CALL shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, &
477 0 : dt, ishake_ext, max_sigma)
478 : END IF
479 : ! Collective Variables
480 0 : IF (gci%ncolv%ntot /= 0) THEN
481 : CALL shake_roll_colv_ext(gci, particle_set, pos, vel, r_shake, &
482 0 : v_shake, dt, ishake_ext, cell, imass, max_sigma)
483 : END IF
484 : ! Virtual Site
485 0 : IF (gci%nvsite /= 0) THEN
486 0 : CPABORT("Virtual Site Constraint/Restraint not implemented for SHAKE_ROLL!")
487 : END IF
488 0 : IF (gci%joint_colvar_checks) THEN
489 : CALL check_local_colvars(local_molecules, molecule_set, particle_set, cell, group, &
490 0 : max_sigma, pos=pos)
491 : END IF
492 0 : CALL restore_temporary_set(particle_set, local_particles, pos=pos, vel=vel)
493 : END IF
494 1826 : CALL shake_ext_info(log_unit, ishake_ext, max_sigma)
495 : END DO Shake_Inter_Loop
496 : CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, &
497 1826 : molecule_kind_set, group, "S")
498 1826 : CALL check_tol(roll_tol, iroll, 'SHAKE', r_shake)
499 1826 : CALL timestop(handle)
500 :
501 1826 : END SUBROUTINE shake_roll_control
502 :
503 : ! **************************************************************************************************
504 : !> \brief ...
505 : !> \param gci ...
506 : !> \param local_molecules ...
507 : !> \param molecule_set ...
508 : !> \param molecule_kind_set ...
509 : !> \param particle_set ...
510 : !> \param vel ...
511 : !> \param dt ...
512 : !> \param simpar ...
513 : !> \param vector ...
514 : !> \param veps ...
515 : !> \param roll_tol ...
516 : !> \param iroll ...
517 : !> \param para_env ...
518 : !> \param u ...
519 : !> \param cell ...
520 : !> \param local_particles ...
521 : !> \par History
522 : !> Teodoro Laino [tlaino] 2007 - Extension to Intermolecular constraints
523 : ! **************************************************************************************************
524 1794 : SUBROUTINE rattle_roll_control(gci, local_molecules, molecule_set, &
525 1794 : molecule_kind_set, particle_set, vel, dt, simpar, vector, &
526 1794 : veps, roll_tol, iroll, para_env, u, cell, local_particles)
527 :
528 : TYPE(global_constraint_type), POINTER :: gci
529 : TYPE(distribution_1d_type), POINTER :: local_molecules
530 : TYPE(molecule_type), POINTER :: molecule_set(:)
531 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
532 : TYPE(particle_type), POINTER :: particle_set(:)
533 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
534 : REAL(KIND=dp), INTENT(IN) :: dt
535 : TYPE(simpar_type), INTENT(IN) :: simpar
536 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vector
537 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: veps
538 : REAL(KIND=dp), INTENT(OUT) :: roll_tol
539 : INTEGER, INTENT(INOUT) :: iroll
540 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
541 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
542 : OPTIONAL :: u
543 : TYPE(cell_type), POINTER :: cell
544 : TYPE(distribution_1d_type), POINTER :: local_particles
545 :
546 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rattle_roll_control'
547 :
548 : INTEGER :: handle, i, ikind, imol, irattle_ext, irattle_int, k, lagrange_mult, log_unit, &
549 : n3x3con, n4x6con, nconstraint, nkind, nmol_per_kind
550 : LOGICAL :: do_ext_constraint, dump_lm
551 : REAL(KIND=dp) :: int_max_sigma, mass, max_sigma, &
552 : rattle_tol
553 : REAL(KIND=dp), DIMENSION(3, 3) :: r_rattle
554 3588 : REAL(KIND=dp), DIMENSION(SIZE(vel, 2)) :: imass
555 : TYPE(atomic_kind_type), POINTER :: atomic_kind
556 : TYPE(colvar_counters) :: ncolv
557 : TYPE(molecule_kind_type), POINTER :: molecule_kind
558 : TYPE(molecule_type), POINTER :: molecule
559 :
560 1794 : CALL timeset(routineN, handle)
561 : ! initialize locals
562 1794 : nkind = SIZE(molecule_kind_set)
563 1794 : rattle_tol = simpar%shake_tol
564 1794 : log_unit = simpar%info_constraint
565 1794 : lagrange_mult = simpar%lagrange_multipliers
566 1794 : dump_lm = simpar%dump_lm
567 : ! setting up for roll
568 1794 : IF (simpar%ensemble == npt_i_ensemble .OR. simpar%ensemble == npt_ia_ensemble) THEN
569 1774 : CALL get_roll_matrix('RATTLE', v_shake=r_rattle, vector_v=vector)
570 20 : ELSE IF (simpar%ensemble == npt_f_ensemble) THEN
571 20 : CALL get_roll_matrix('RATTLE', v_shake=r_rattle, vector_v=vector, u=u)
572 : END IF
573 625048 : DO k = 1, SIZE(vel, 2)
574 623254 : atomic_kind => particle_set(k)%atomic_kind
575 623254 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
576 625048 : imass(k) = 1.0_dp/MAX(mass, EPSILON(0.0_dp))
577 : END DO
578 1794 : do_ext_constraint = (gci%ntot /= 0)
579 1794 : irattle_ext = 0
580 1794 : max_sigma = -1.0E+10_dp
581 3588 : Rattle_Inter_Loop: DO WHILE (ABS(max_sigma) >= rattle_tol)
582 1794 : max_sigma = 0.0_dp
583 1794 : irattle_ext = irattle_ext + 1
584 : ! Intramolecular Constraints
585 3606 : MOL: DO ikind = 1, nkind
586 1812 : nmol_per_kind = local_molecules%n_el(ikind)
587 104736 : DO imol = 1, nmol_per_kind
588 101130 : i = local_molecules%list(ikind)%array(imol)
589 101130 : molecule => molecule_set(i)
590 101130 : molecule_kind => molecule%molecule_kind
591 : CALL get_molecule_kind(molecule_kind, ncolv=ncolv, &
592 : ng3x3=n3x3con, ng4x6=n4x6con, &
593 101130 : nconstraint=nconstraint)
594 101130 : IF (nconstraint == 0) CYCLE
595 100310 : int_max_sigma = -1.0E+10_dp
596 100310 : irattle_int = MERGE(0, 1, irattle_ext == 1 .OR. .NOT. gci%joint_colvar_checks)
597 204854 : Rattle_Roll_Intramolecular: DO WHILE (ABS(int_max_sigma) >= rattle_tol)
598 104544 : int_max_sigma = 0.0_dp
599 104544 : irattle_int = irattle_int + 1
600 : ! 3x3
601 104544 : IF (n3x3con /= 0) THEN
602 : CALL rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, &
603 97999 : veps)
604 : END IF
605 : ! 4x6
606 104544 : IF (n4x6con /= 0) THEN
607 : CALL rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, &
608 1024 : veps)
609 : END IF
610 : ! Collective Variables
611 204854 : IF (ncolv%ntot /= 0) THEN
612 : CALL rattle_roll_colv_int(molecule, particle_set, vel, r_rattle, dt, &
613 5521 : irattle_int, veps, cell, imass, int_max_sigma)
614 : END IF
615 : END DO Rattle_Roll_Intramolecular
616 100310 : max_sigma = MAX(max_sigma, int_max_sigma)
617 203252 : CALL rattle_int_info(log_unit, i, irattle_int, max_sigma)
618 : END DO
619 : END DO MOL
620 : ! Intermolecular Constraints
621 1794 : IF (do_ext_constraint) THEN
622 0 : CALL update_temporary_set(para_env, vel=vel)
623 : ! 3x3
624 0 : IF (gci%ng3x3 /= 0) THEN
625 : CALL rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, &
626 0 : veps)
627 : END IF
628 : ! 4x6
629 0 : IF (gci%ng4x6 /= 0) THEN
630 : CALL rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, &
631 0 : veps)
632 : END IF
633 : ! Collective Variables
634 0 : IF (gci%ncolv%ntot /= 0) THEN
635 : CALL rattle_roll_colv_ext(gci, particle_set, vel, r_rattle, dt, &
636 0 : irattle_ext, veps, cell, imass, max_sigma)
637 : END IF
638 0 : IF (gci%joint_colvar_checks) THEN
639 : CALL check_local_colvars(local_molecules, molecule_set, particle_set, cell, para_env, &
640 0 : max_sigma, vel=vel, rmat=r_rattle, veps=veps)
641 : END IF
642 0 : CALL restore_temporary_set(particle_set, local_particles, vel=vel)
643 : END IF
644 1794 : CALL rattle_ext_info(log_unit, irattle_ext, max_sigma)
645 : END DO Rattle_Inter_Loop
646 : CALL dump_lagrange_mult(dump_lm, lagrange_mult, local_molecules, molecule_set, gci, &
647 1794 : molecule_kind_set, para_env, "R")
648 1794 : CALL check_tol(roll_tol, iroll, 'RATTLE', veps=veps)
649 1794 : CALL timestop(handle)
650 1794 : END SUBROUTINE rattle_roll_control
651 :
652 : ! **************************************************************************************************
653 : !> \brief Recheck local collective variables after global constraint corrections.
654 : !> \param local_molecules ...
655 : !> \param molecule_set ...
656 : !> \param particle_set ...
657 : !> \param cell ...
658 : !> \param group ...
659 : !> \param max_sigma ...
660 : !> \param pos ...
661 : !> \param vel ...
662 : !> \param rmat ...
663 : !> \param veps ...
664 : ! **************************************************************************************************
665 378 : SUBROUTINE check_local_colvars(local_molecules, molecule_set, particle_set, cell, group, &
666 378 : max_sigma, pos, vel, rmat, veps)
667 : TYPE(distribution_1d_type), POINTER :: local_molecules
668 : TYPE(molecule_type), POINTER :: molecule_set(:)
669 : TYPE(particle_type), POINTER :: particle_set(:)
670 : TYPE(cell_type), POINTER :: cell
671 :
672 : CLASS(mp_comm_type), INTENT(IN) :: group
673 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
674 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: pos(:, :), vel(:, :), rmat(:, :), veps(:, :)
675 :
676 : INTEGER :: ikind, imol, i
677 : TYPE(molecule_type), POINTER :: molecule
678 :
679 : ! Global corrections can invalidate already-converged local collective variables.
680 1134 : DO ikind = 1, SIZE(local_molecules%n_el)
681 1512 : DO imol = 1, local_molecules%n_el(ikind)
682 378 : i = local_molecules%list(ikind)%array(imol)
683 378 : molecule => molecule_set(i)
684 : CALL check_colv_int(molecule, particle_set, cell, max_sigma, pos=pos, vel=vel, &
685 2268 : rmat=rmat, veps=veps)
686 : END DO
687 : END DO
688 378 : CALL group%max(max_sigma)
689 378 : END SUBROUTINE check_local_colvars
690 :
691 : ! **************************************************************************************************
692 : !> \brief ...
693 : !> \param dump_lm ...
694 : !> \param log_unit ...
695 : !> \param local_molecules ...
696 : !> \param molecule_set ...
697 : !> \param gci ...
698 : !> \param molecule_kind_set ...
699 : !> \param group ...
700 : !> \param id_type ...
701 : !> \par History
702 : !> Teodoro Laino [tlaino] 2007 - Dumps lagrange multipliers
703 : ! **************************************************************************************************
704 72848 : SUBROUTINE dump_lagrange_mult(dump_lm, log_unit, local_molecules, molecule_set, gci, &
705 : molecule_kind_set, group, id_type)
706 : LOGICAL, INTENT(IN) :: dump_lm
707 : INTEGER, INTENT(IN) :: log_unit
708 : TYPE(distribution_1d_type), POINTER :: local_molecules
709 : TYPE(molecule_type), POINTER :: molecule_set(:)
710 : TYPE(global_constraint_type), POINTER :: gci
711 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
712 :
713 : CLASS(mp_comm_type), INTENT(IN) :: group
714 : CHARACTER(LEN=1), INTENT(IN) :: id_type
715 :
716 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_lagrange_mult'
717 :
718 : CHARACTER(LEN=default_string_length) :: label
719 : INTEGER :: handle, i, ikind, imol, j, my_index, &
720 : n3x3con, n4x6con, nconstraint, nkind
721 : LOGICAL :: do_ext_constraint, do_int_constraint
722 36424 : REAL(KIND=dp), DIMENSION(:), POINTER :: lagr
723 : TYPE(colvar_counters) :: ncolv
724 : TYPE(molecule_kind_type), POINTER :: molecule_kind
725 : TYPE(molecule_type), POINTER :: molecule
726 :
727 36424 : CALL timeset(routineN, handle)
728 : ! Total number of intramolecular constraints (distributed)
729 : CALL get_molecule_kind_set(molecule_kind_set=molecule_kind_set, &
730 36424 : nconstraint=nconstraint)
731 36424 : do_int_constraint = (nconstraint > 0)
732 36424 : do_ext_constraint = (gci%ntot > 0)
733 36424 : IF (dump_lm .AND. (do_int_constraint .OR. do_ext_constraint)) THEN
734 402 : nkind = SIZE(molecule_kind_set)
735 830 : ALLOCATE (lagr(nconstraint))
736 2898 : lagr = 0.0_dp
737 : ! Dump lagrange multipliers for Intramolecular Constraints
738 402 : my_index = 0
739 402 : IF (do_int_constraint) THEN
740 52 : MOL: DO ikind = 1, nkind
741 26 : molecule_kind => molecule_kind_set(ikind)
742 : CALL get_molecule_kind(molecule_kind, &
743 : ncolv=ncolv, &
744 : ng3x3=n3x3con, &
745 26 : ng4x6=n4x6con)
746 884 : DO imol = 1, molecule_kind%nmolecule
747 832 : i = molecule_kind%molecule_list(imol)
748 10634 : IF (ANY(local_molecules%list(ikind)%array == i)) THEN
749 416 : molecule => molecule_set(i)
750 : ! Collective Variables
751 416 : DO j = 1, ncolv%ntot
752 0 : lagr(my_index + 1) = molecule%lci%lcolv(j)%lambda
753 416 : my_index = my_index + 1
754 : END DO
755 : ! 3x3
756 832 : DO j = 1, n3x3con
757 1664 : lagr(my_index + 1:my_index + 3) = molecule%lci%lg3x3(j)%lambda(:)
758 832 : my_index = my_index + 3
759 : END DO
760 : ! 4x6
761 416 : DO j = 1, n4x6con
762 0 : lagr(my_index + 1:my_index + 6) = molecule%lci%lg4x6(j)%lambda(:)
763 416 : my_index = my_index + 6
764 : END DO
765 : ELSE
766 416 : my_index = my_index + ncolv%ntot + 3*n3x3con + 6*n4x6con
767 : END IF
768 : END DO
769 : END DO MOL
770 5018 : CALL group%sum(lagr)
771 : END IF
772 : ! Intermolecular constraints
773 402 : IF (do_ext_constraint) THEN
774 376 : CALL reallocate(lagr, 1, SIZE(lagr) + gci%ntot)
775 : ! Collective Variables
776 752 : DO j = 1, gci%ncolv%ntot
777 376 : lagr(my_index + 1) = gci%lcolv(j)%lambda
778 752 : my_index = my_index + 1
779 : END DO
780 : ! 3x3
781 376 : DO j = 1, gci%ng3x3
782 0 : lagr(my_index + 1:my_index + 3) = gci%lg3x3(j)%lambda(:)
783 376 : my_index = my_index + 3
784 : END DO
785 : ! 4x6
786 376 : DO j = 1, gci%ng4x6
787 0 : lagr(my_index + 1:my_index + 6) = gci%lg4x6(j)%lambda(:)
788 376 : my_index = my_index + 6
789 : END DO
790 : END IF
791 402 : IF (log_unit > 0) THEN
792 229 : IF (id_type == "S") THEN
793 115 : label = "Shake Lagrangian Multipliers:"
794 114 : ELSE IF (id_type == "R") THEN
795 114 : label = "Rattle Lagrangian Multipliers:"
796 : ELSE
797 0 : CPABORT("Only S for Shake or R for Rattle are supported for Lagrangian Multipliers")
798 : END IF
799 497 : WRITE (log_unit, FMT='(A,T40,4F15.9)') TRIM(label), lagr(1:MIN(4, SIZE(lagr)))
800 528 : DO j = 5, SIZE(lagr), 4
801 1724 : WRITE (log_unit, FMT='(T40,4F15.9)') lagr(j:MIN(j + 3, SIZE(lagr)))
802 : END DO
803 : END IF
804 402 : DEALLOCATE (lagr)
805 : END IF
806 36424 : CALL timestop(handle)
807 :
808 36424 : END SUBROUTINE dump_lagrange_mult
809 :
810 : ! **************************************************************************************************
811 : !> \brief Dumps convergence info about shake - intramolecular constraint loop
812 : !> \param log_unit ...
813 : !> \param i ...
814 : !> \param ishake_int ...
815 : !> \param max_sigma ...
816 : !> \par History
817 : !> Teodoro Laino [tlaino] 2007 - University of Zurich
818 : ! **************************************************************************************************
819 261213 : SUBROUTINE shake_int_info(log_unit, i, ishake_int, max_sigma)
820 : INTEGER, INTENT(IN) :: log_unit, i, ishake_int
821 : REAL(KIND=dp), INTENT(IN) :: max_sigma
822 :
823 261213 : IF (log_unit > 0) THEN
824 : ! Dump info if requested
825 : WRITE (log_unit, '("SHAKE_INFO|",2X,2(A,I6),A,F15.9)') &
826 117 : "Molecule Nr.:", i, " Nr. Iterations:", ishake_int, " Max. Err.:", max_sigma
827 : END IF
828 : ! Notify a not converged SHAKE
829 261213 : IF (ishake_int > Max_Shake_Iter) THEN
830 : CALL cp_warn(__LOCATION__, &
831 : "Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
832 : "intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// &
833 0 : ". CP2K continues but results could be meaningless. ")
834 : END IF
835 261213 : END SUBROUTINE shake_int_info
836 :
837 : ! **************************************************************************************************
838 : !> \brief Dumps convergence info about shake - intermolecular constraint loop
839 : !> \param log_unit ...
840 : !> \param ishake_ext ...
841 : !> \param max_sigma ...
842 : !> \par History
843 : !> Teodoro Laino [tlaino] 2007 - University of Zurich
844 : ! **************************************************************************************************
845 18923 : SUBROUTINE shake_ext_info(log_unit, ishake_ext, max_sigma)
846 : INTEGER, INTENT(IN) :: log_unit, ishake_ext
847 : REAL(KIND=dp), INTENT(IN) :: max_sigma
848 :
849 18923 : IF (log_unit > 0) THEN
850 : ! Dump info if requested
851 : WRITE (log_unit, '("SHAKE_INFO|",2X,A,I6,A,F15.9)') &
852 12 : "External Shake Nr. Iterations:", ishake_ext, &
853 24 : " Max. Err.:", max_sigma
854 : END IF
855 : ! Notify a not converged SHAKE
856 18923 : IF (ishake_ext > Max_Shake_Iter) THEN
857 : CALL cp_warn(__LOCATION__, &
858 : "Shake NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
859 0 : "intermolecular constraint. CP2K continues but results could be meaningless.")
860 : END IF
861 18923 : END SUBROUTINE shake_ext_info
862 :
863 : ! **************************************************************************************************
864 : !> \brief Dumps convergence info about rattle - intramolecular constraint loop
865 : !> \param log_unit ...
866 : !> \param i ...
867 : !> \param irattle_int ...
868 : !> \param max_sigma ...
869 : !> \par History
870 : !> Teodoro Laino [tlaino] 2007 - University of Zurich
871 : ! **************************************************************************************************
872 246407 : SUBROUTINE rattle_int_info(log_unit, i, irattle_int, max_sigma)
873 : INTEGER, INTENT(IN) :: log_unit, i, irattle_int
874 : REAL(KIND=dp), INTENT(IN) :: max_sigma
875 :
876 246407 : IF (log_unit > 0) THEN
877 : ! Dump info if requested
878 : WRITE (log_unit, '("RATTLE_INFO|",1X,2(A,I6),A,F15.9)') &
879 101 : "Molecule Nr.:", i, " Nr. Iterations:", irattle_int, " Max. Err.:", max_sigma
880 : END IF
881 : ! Notify a not converged RATTLE
882 246407 : IF (irattle_int > Max_shake_Iter) THEN
883 : CALL cp_warn(__LOCATION__, &
884 : "Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
885 : "intramolecular constraint loop for Molecule nr. "//cp_to_string(i)// &
886 0 : ". CP2K continues but results could be meaningless.")
887 : END IF
888 246407 : END SUBROUTINE rattle_int_info
889 :
890 : ! **************************************************************************************************
891 : !> \brief Dumps convergence info about rattle - intermolecular constraint loop
892 : !> \param log_unit ...
893 : !> \param irattle_ext ...
894 : !> \param max_sigma ...
895 : !> \par History
896 : !> Teodoro Laino [tlaino] 2007 - University of Zurich
897 : ! **************************************************************************************************
898 18678 : SUBROUTINE rattle_ext_info(log_unit, irattle_ext, max_sigma)
899 : INTEGER, INTENT(IN) :: log_unit, irattle_ext
900 : REAL(KIND=dp), INTENT(IN) :: max_sigma
901 :
902 18678 : IF (log_unit > 0) THEN
903 : ! Dump info if requested
904 : WRITE (log_unit, '("RATTLE_INFO|",1X,A,I6,A,F15.9)') &
905 11 : "External Rattle Nr. Iterations:", irattle_ext, &
906 22 : " Max. Err.:", max_sigma
907 : END IF
908 : ! Notify a not converged RATTLE
909 18678 : IF (irattle_ext > Max_shake_Iter) THEN
910 : CALL cp_warn(__LOCATION__, &
911 : "Rattle NOT converged in "//cp_to_string(Max_Shake_Iter)//" iterations in the "// &
912 0 : "intermolecular constraint. CP2K continues but results could be meaningless.")
913 : END IF
914 18678 : END SUBROUTINE rattle_ext_info
915 :
916 : ! **************************************************************************************************
917 : !> \brief Updates the TARGET of the COLLECTIVE constraints if the growth speed
918 : !> is different from zero.
919 : !> \param gci ...
920 : !> \param local_molecules ...
921 : !> \param molecule_set ...
922 : !> \param molecule_kind_set ...
923 : !> \param dt ...
924 : !> \param root_section ...
925 : !> \date 02.2008
926 : !> \author Teodoro Laino [tlaino] - University of Zurich
927 : ! **************************************************************************************************
928 17054 : SUBROUTINE shake_update_targets(gci, local_molecules, molecule_set, &
929 : molecule_kind_set, dt, root_section)
930 :
931 : TYPE(global_constraint_type), POINTER :: gci
932 : TYPE(distribution_1d_type), POINTER :: local_molecules
933 : TYPE(molecule_type), POINTER :: molecule_set(:)
934 : TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
935 : REAL(kind=dp), INTENT(in) :: dt
936 : TYPE(section_vals_type), POINTER :: root_section
937 :
938 : CHARACTER(LEN=*), PARAMETER :: routineN = 'shake_update_targets'
939 :
940 : INTEGER :: handle, ikind, nkind
941 : LOGICAL :: do_ext_constraint
942 : TYPE(section_vals_type), POINTER :: motion_section
943 :
944 17054 : CALL timeset(routineN, handle)
945 17054 : CPASSERT(ASSOCIATED(local_molecules))
946 17054 : CPASSERT(ASSOCIATED(molecule_set))
947 17054 : motion_section => section_vals_get_subs_vals(root_section, "MOTION")
948 17054 : CALL shake_update_colv_targets(dt, motion_section)
949 17054 : nkind = SIZE(molecule_kind_set)
950 17054 : do_ext_constraint = (gci%ntot /= 0)
951 : ! Intramolecular Constraints
952 84856 : MOL: DO ikind = 1, nkind
953 84856 : IF (ASSOCIATED(molecule_kind_set(ikind)%colv_list)) THEN
954 51676 : CALL shake_update_colv_low(molecule_kind_set(ikind)%colv_list, dt, motion_section)
955 : END IF
956 : END DO MOL
957 : ! Intermolecular constraints
958 17054 : IF (do_ext_constraint) THEN
959 : ! Collective Variables
960 1314 : IF (gci%ncolv%ntot /= 0) CALL shake_update_colv_ext(gci, dt, motion_section)
961 : END IF
962 17054 : CALL timestop(handle)
963 17054 : END SUBROUTINE shake_update_targets
964 :
965 : END MODULE constraint
|
