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