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