(git:374b731)
Loading...
Searching...
No Matches
constraint_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Contains routines useful for the application of constraints during MD
10!> \par History
11!> none
12! **************************************************************************************************
13MODULE constraint_util
14 USE cell_types, ONLY: cell_type
15 USE colvar_methods, ONLY: colvar_eval_mol_f
16 USE distribution_1d_types, ONLY: distribution_1d_type
17 USE kinds, ONLY: dp
18 USE message_passing, ONLY: mp_comm_type
19 USE molecule_kind_types, ONLY: colvar_constraint_type,&
20 fixd_constraint_type,&
21 g3x3_constraint_type,&
22 g4x6_constraint_type,&
23 get_molecule_kind,&
24 molecule_kind_type
25 USE molecule_types, ONLY: get_molecule,&
26 global_constraint_type,&
27 local_colvar_constraint_type,&
28 local_constraint_type,&
29 local_g3x3_constraint_type,&
30 local_g4x6_constraint_type,&
31 molecule_type
32 USE particle_types, ONLY: particle_type
33 USE virial_types, ONLY: virial_type
34#include "./base/base_uses.f90"
35
36 IMPLICIT NONE
37
38 PRIVATE
39 PUBLIC :: getold, &
40 pv_constraint, &
41 check_tol, &
42 get_roll_matrix, &
43 restore_temporary_set, &
44 update_temporary_set
45
46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_util'
47
48CONTAINS
49
50! **************************************************************************************************
51!> \brief saves all of the old variables
52!> \param gci ...
53!> \param local_molecules ...
54!> \param molecule_set ...
55!> \param molecule_kind_set ...
56!> \param particle_set ...
57!> \param cell ...
58!> \par History
59!> none
60! **************************************************************************************************
61 SUBROUTINE getold(gci, local_molecules, molecule_set, molecule_kind_set, &
62 particle_set, cell)
63
64 TYPE(global_constraint_type), POINTER :: gci
65 TYPE(distribution_1d_type), POINTER :: local_molecules
66 TYPE(molecule_type), POINTER :: molecule_set(:)
67 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
68 TYPE(particle_type), POINTER :: particle_set(:)
69 TYPE(cell_type), POINTER :: cell
70
71 INTEGER :: first_atom, i, ikind, imol, n3x3con, &
72 n4x6con, nkind, nmol_per_kind
73 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
74 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
75 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
76 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
77 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
78 TYPE(local_constraint_type), POINTER :: lci
79 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
80 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
81 TYPE(molecule_kind_type), POINTER :: molecule_kind
82 TYPE(molecule_type), POINTER :: molecule
83
84 NULLIFY (fixd_list)
85 nkind = SIZE(molecule_kind_set)
86 ! Intramolecular constraints
87 mol: DO ikind = 1, nkind
88 nmol_per_kind = local_molecules%n_el(ikind)
89 DO imol = 1, nmol_per_kind
90 i = local_molecules%list(ikind)%array(imol)
91 molecule => molecule_set(i)
92 molecule_kind => molecule%molecule_kind
93 CALL get_molecule_kind(molecule_kind, ng3x3=n3x3con, ng4x6=n4x6con, &
94 colv_list=colv_list, g3x3_list=g3x3_list, g4x6_list=g4x6_list, &
95 fixd_list=fixd_list)
96 CALL get_molecule(molecule, lci=lci)
97 IF (.NOT. ASSOCIATED(lci)) cycle
98 CALL get_molecule(molecule, first_atom=first_atom, &
99 lcolv=lcolv, lg3x3=lg3x3, lg4x6=lg4x6)
100 CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
101 lcolv, lg3x3, lg4x6, first_atom, particle_set, cell)
102 END DO
103 END DO mol
104 ! Intermolecular constraints
105 IF (gci%ntot > 0) THEN
106 n3x3con = gci%ng3x3
107 n4x6con = gci%ng4x6
108 colv_list => gci%colv_list
109 g3x3_list => gci%g3x3_list
110 g4x6_list => gci%g4x6_list
111 fixd_list => gci%fixd_list
112 lcolv => gci%lcolv
113 lg3x3 => gci%lg3x3
114 lg4x6 => gci%lg4x6
115 CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
116 lcolv, lg3x3, lg4x6, 1, particle_set, cell)
117 END IF
118 END SUBROUTINE getold
119
120! **************************************************************************************************
121!> \brief saves all of the old variables - Low Level
122!> \param n3x3con ...
123!> \param n4x6con ...
124!> \param colv_list ...
125!> \param g3x3_list ...
126!> \param g4x6_list ...
127!> \param fixd_list ...
128!> \param lcolv ...
129!> \param lg3x3 ...
130!> \param lg4x6 ...
131!> \param first_atom ...
132!> \param particle_set ...
133!> \param cell ...
134!> \par History
135!> none
136! **************************************************************************************************
137 SUBROUTINE getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
138 lcolv, lg3x3, lg4x6, first_atom, particle_set, cell)
139
140 INTEGER, INTENT(IN) :: n3x3con, n4x6con
141 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
142 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
143 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
144 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
145 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
146 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
147 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
148 INTEGER, INTENT(IN) :: first_atom
149 TYPE(particle_type), POINTER :: particle_set(:)
150 TYPE(cell_type), POINTER :: cell
151
152 INTEGER :: iconst, index
153
154 IF (ASSOCIATED(colv_list)) THEN
155 ! Collective constraints
156 DO iconst = 1, SIZE(colv_list)
157 CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
158 particles=particle_set, fixd_list=fixd_list)
159 END DO
160 END IF
161 ! 3x3 constraints
162 DO iconst = 1, n3x3con
163 index = g3x3_list(iconst)%a + first_atom - 1
164 lg3x3(iconst)%ra_old = particle_set(index)%r
165 index = g3x3_list(iconst)%b + first_atom - 1
166 lg3x3(iconst)%rb_old = particle_set(index)%r
167 index = g3x3_list(iconst)%c + first_atom - 1
168 lg3x3(iconst)%rc_old = particle_set(index)%r
169 END DO
170 ! 4x6 constraints
171 DO iconst = 1, n4x6con
172 index = g4x6_list(iconst)%a + first_atom - 1
173 lg4x6(iconst)%ra_old = particle_set(index)%r
174 index = g4x6_list(iconst)%b + first_atom - 1
175 lg4x6(iconst)%rb_old = particle_set(index)%r
176 index = g4x6_list(iconst)%c + first_atom - 1
177 lg4x6(iconst)%rc_old = particle_set(index)%r
178 index = g4x6_list(iconst)%d + first_atom - 1
179 lg4x6(iconst)%rd_old = particle_set(index)%r
180 END DO
181
182 END SUBROUTINE getold_low
183
184! **************************************************************************************************
185!> \brief ...
186!> \param gci ...
187!> \param local_molecules ...
188!> \param molecule_set ...
189!> \param molecule_kind_set ...
190!> \param particle_set ...
191!> \param virial ...
192!> \param group ...
193!> \par History
194!> none
195! **************************************************************************************************
196 SUBROUTINE pv_constraint(gci, local_molecules, molecule_set, molecule_kind_set, &
197 particle_set, virial, group)
198
199 TYPE(global_constraint_type), POINTER :: gci
200 TYPE(distribution_1d_type), POINTER :: local_molecules
201 TYPE(molecule_type), POINTER :: molecule_set(:)
202 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:)
203 TYPE(particle_type), POINTER :: particle_set(:)
204 TYPE(virial_type), INTENT(INOUT) :: virial
205
206 CLASS(mp_comm_type), INTENT(IN) :: group
207
208 INTEGER :: first_atom, i, ikind, imol, ng3x3, &
209 ng4x6, nkind, nmol_per_kind
210 REAL(kind=dp) :: pv(3, 3)
211 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
212 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
213 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
214 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
215 TYPE(local_constraint_type), POINTER :: lci
216 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
217 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
218 TYPE(molecule_kind_type), POINTER :: molecule_kind
219 TYPE(molecule_type), POINTER :: molecule
220
221 pv = 0.0_dp
222 nkind = SIZE(molecule_kind_set)
223 ! Intramolecular Constraints
224 mol: DO ikind = 1, nkind
225 nmol_per_kind = local_molecules%n_el(ikind)
226 DO imol = 1, nmol_per_kind
227 i = local_molecules%list(ikind)%array(imol)
228 molecule => molecule_set(i)
229 molecule_kind => molecule%molecule_kind
230 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
231 ng4x6=ng4x6, g3x3_list=g3x3_list, g4x6_list=g4x6_list, &
232 colv_list=colv_list)
233 CALL get_molecule(molecule, lci=lci)
234 IF (.NOT. ASSOCIATED(lci)) cycle
235 CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3, &
236 lg4x6=lg4x6, lcolv=lcolv)
237 CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
238 first_atom, lg3x3, lg4x6, lcolv, particle_set, pv)
239 END DO
240 END DO mol
241 ! Intermolecular constraints
242 IF (gci%ntot > 0) THEN
243 ng3x3 = gci%ng3x3
244 ng4x6 = gci%ng4x6
245 colv_list => gci%colv_list
246 g3x3_list => gci%g3x3_list
247 g4x6_list => gci%g4x6_list
248 lcolv => gci%lcolv
249 lg3x3 => gci%lg3x3
250 lg4x6 => gci%lg4x6
251 CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
252 1, lg3x3, lg4x6, lcolv, particle_set, pv)
253 END IF
254 CALL group%sum(pv)
255 ! Symmetrize PV
256 pv(1, 2) = 0.5_dp*(pv(1, 2) + pv(2, 1))
257 pv(2, 1) = pv(1, 2)
258 pv(1, 3) = 0.5_dp*(pv(1, 3) + pv(3, 1))
259 pv(3, 1) = pv(1, 3)
260 pv(3, 2) = 0.5_dp*(pv(3, 2) + pv(2, 3))
261 pv(2, 3) = pv(3, 2)
262 ! Store in virial type
263 virial%pv_constraint = pv
264
265 END SUBROUTINE pv_constraint
266
267! **************************************************************************************************
268!> \brief ...
269!> \param ng3x3 ...
270!> \param ng4x6 ...
271!> \param g3x3_list ...
272!> \param g4x6_list ...
273!> \param colv_list ...
274!> \param first_atom ...
275!> \param lg3x3 ...
276!> \param lg4x6 ...
277!> \param lcolv ...
278!> \param particle_set ...
279!> \param pv ...
280!> \par History
281!> none
282! **************************************************************************************************
283 SUBROUTINE pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
284 first_atom, lg3x3, lg4x6, lcolv, particle_set, pv)
285
286 INTEGER, INTENT(IN) :: ng3x3, ng4x6
287 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
288 TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
289 TYPE(colvar_constraint_type), POINTER :: colv_list(:)
290 INTEGER, INTENT(IN) :: first_atom
291 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
292 TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
293 TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
294 TYPE(particle_type), POINTER :: particle_set(:)
295 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv
296
297 INTEGER :: iconst, index_a, index_b, index_c, &
298 index_d
299 REAL(kind=dp) :: fc1(3), fc2(3), fc3(3), fc4(3), &
300 lambda_3x3(3), lambda_4x6(6)
301
302 IF (ASSOCIATED(colv_list)) THEN
303 ! Colvar Constraints
304 DO iconst = 1, SIZE(colv_list)
305 CALL pv_colv_eval(pv, lcolv(iconst), particle_set)
306 END DO
307 END IF
308 ! 3x3
309 DO iconst = 1, ng3x3
310 ! pv gets updated with FULL multiplier
311 lambda_3x3 = lg3x3(iconst)%lambda
312
313 fc1 = lambda_3x3(1)*lg3x3(iconst)%fa + &
314 lambda_3x3(2)*lg3x3(iconst)%fb
315 fc2 = -lambda_3x3(1)*lg3x3(iconst)%fa + &
316 lambda_3x3(3)*lg3x3(iconst)%fc
317 fc3 = -lambda_3x3(2)*lg3x3(iconst)%fb - &
318 lambda_3x3(3)*lg3x3(iconst)%fc
319 index_a = g3x3_list(iconst)%a + first_atom - 1
320 index_b = g3x3_list(iconst)%b + first_atom - 1
321 index_c = g3x3_list(iconst)%c + first_atom - 1
322
323 !pv(1,1)
324 pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
325 pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
326 pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
327 !pv(1,2)
328 pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
329 pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
330 pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
331 !pv(1,3)
332 pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
333 pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
334 pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
335 !pv(2,1)
336 pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
337 pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
338 pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
339 !pv(2,2)
340 pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
341 pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
342 pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
343 !pv(2,3)
344 pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
345 pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
346 pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
347 !pv(3,1)
348 pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
349 pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
350 pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
351 !pv(3,2)
352 pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
353 pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
354 pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
355 !pv(3,3)
356 pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
357 pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
358 pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
359 END DO
360
361 ! 4x6
362 DO iconst = 1, ng4x6
363 ! pv gets updated with FULL multiplier
364 lambda_4x6 = lg4x6(iconst)%lambda
365
366 fc1 = lambda_4x6(1)*lg4x6(iconst)%fa + &
367 lambda_4x6(2)*lg4x6(iconst)%fb + &
368 lambda_4x6(3)*lg4x6(iconst)%fc
369 fc2 = -lambda_4x6(1)*lg4x6(iconst)%fa + &
370 lambda_4x6(4)*lg4x6(iconst)%fd + &
371 lambda_4x6(5)*lg4x6(iconst)%fe
372 fc3 = -lambda_4x6(2)*lg4x6(iconst)%fb - &
373 lambda_4x6(4)*lg4x6(iconst)%fd + &
374 lambda_4x6(6)*lg4x6(iconst)%ff
375 fc4 = -lambda_4x6(3)*lg4x6(iconst)%fc - &
376 lambda_4x6(5)*lg4x6(iconst)%fe - &
377 lambda_4x6(6)*lg4x6(iconst)%ff
378 index_a = g4x6_list(iconst)%a + first_atom - 1
379 index_b = g4x6_list(iconst)%b + first_atom - 1
380 index_c = g4x6_list(iconst)%c + first_atom - 1
381 index_d = g4x6_list(iconst)%d + first_atom - 1
382
383 !pv(1,1)
384 pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
385 pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
386 pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
387 pv(1, 1) = pv(1, 1) + fc4(1)*particle_set(index_d)%r(1)
388 !pv(1,2)
389 pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
390 pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
391 pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
392 pv(1, 2) = pv(1, 2) + fc4(1)*particle_set(index_d)%r(2)
393 !pv(1,3)
394 pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
395 pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
396 pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
397 pv(1, 3) = pv(1, 3) + fc4(1)*particle_set(index_d)%r(3)
398 !pv(2,1)
399 pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
400 pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
401 pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
402 pv(2, 1) = pv(2, 1) + fc4(2)*particle_set(index_d)%r(1)
403 !pv(2,2)
404 pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
405 pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
406 pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
407 pv(2, 2) = pv(2, 2) + fc4(2)*particle_set(index_d)%r(2)
408 !pv(2,3)
409 pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
410 pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
411 pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
412 pv(2, 3) = pv(2, 3) + fc4(2)*particle_set(index_d)%r(3)
413 !pv(3,1)
414 pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
415 pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
416 pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
417 pv(3, 1) = pv(3, 1) + fc4(3)*particle_set(index_d)%r(1)
418 !pv(3,2)
419 pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
420 pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
421 pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
422 pv(3, 2) = pv(3, 2) + fc4(3)*particle_set(index_d)%r(2)
423 !pv(3,3)
424 pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
425 pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
426 pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
427 pv(3, 3) = pv(3, 3) + fc4(3)*particle_set(index_d)%r(3)
428 END DO
429
430 END SUBROUTINE pv_constraint_low
431
432! **************************************************************************************************
433!> \brief ...
434!> \param pv ...
435!> \param lcolv ...
436!> \param particle_set ...
437!> \par History
438!> none
439! **************************************************************************************************
440 SUBROUTINE pv_colv_eval(pv, lcolv, particle_set)
441 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv
442 TYPE(local_colvar_constraint_type), INTENT(IN) :: lcolv
443 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
444
445 INTEGER :: i, iatm, ind, j
446 REAL(kind=dp) :: lambda, tmp
447 REAL(kind=dp), DIMENSION(3) :: f
448
449 DO iatm = 1, SIZE(lcolv%colvar_old%i_atom)
450 ind = lcolv%colvar_old%i_atom(iatm)
451 f = -lcolv%colvar_old%dsdr(:, iatm)
452 ! pv gets updated with FULL multiplier
453 lambda = lcolv%lambda
454 DO i = 1, 3
455 tmp = lambda*particle_set(ind)%r(i)
456 DO j = 1, 3
457 pv(j, i) = pv(j, i) + f(j)*tmp
458 END DO
459 END DO
460 END DO
461
462 END SUBROUTINE pv_colv_eval
463
464! **************************************************************************************************
465!> \brief ...
466!> \param roll_tol ...
467!> \param iroll ...
468!> \param char ...
469!> \param matrix ...
470!> \param veps ...
471!> \par History
472!> none
473! **************************************************************************************************
474 SUBROUTINE check_tol(roll_tol, iroll, char, matrix, veps)
475
476 REAL(kind=dp), INTENT(OUT) :: roll_tol
477 INTEGER, INTENT(INOUT) :: iroll
478 CHARACTER(LEN=*), INTENT(IN) :: char
479 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
480 OPTIONAL :: matrix, veps
481
482 REAL(kind=dp) :: local_tol
483 REAL(kind=dp), DIMENSION(3, 3) :: diff_rattle, diff_shake
484 REAL(kind=dp), DIMENSION(3, 3), SAVE :: matrix_old, veps_old
485
486 SELECT CASE (char)
487 CASE ('SHAKE')
488 IF (iroll == 1) THEN
489 matrix_old = matrix
490 roll_tol = -1.e10_dp
491 ELSE
492 roll_tol = 0.0_dp
493 diff_shake = abs(matrix_old - matrix)
494 local_tol = maxval(diff_shake)
495 roll_tol = max(roll_tol, local_tol)
496 matrix_old = matrix
497 END IF
498 iroll = iroll + 1
499 CASE ('RATTLE')
500 IF (iroll == 1) THEN
501 veps_old = veps
502 roll_tol = -1.e+10_dp
503 ELSE
504 roll_tol = 0.0_dp
505 ! compute tolerance on veps
506 diff_rattle = abs(veps - veps_old)
507 local_tol = maxval(diff_rattle)
508 roll_tol = max(roll_tol, local_tol)
509 veps_old = veps
510 END IF
511 iroll = iroll + 1
512 END SELECT
513
514 END SUBROUTINE check_tol
515
516! **************************************************************************************************
517!> \brief ...
518!> \param char ...
519!> \param r_shake ...
520!> \param v_shake ...
521!> \param vector_r ...
522!> \param vector_v ...
523!> \param u ...
524!> \par History
525!> none
526! **************************************************************************************************
527 SUBROUTINE get_roll_matrix(char, r_shake, v_shake, vector_r, vector_v, u)
528
529 CHARACTER(len=*), INTENT(IN) :: char
530 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
531 OPTIONAL :: r_shake, v_shake
532 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: vector_r, vector_v
533 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
534 OPTIONAL :: u
535
536 INTEGER :: i
537 REAL(kind=dp), DIMENSION(3, 3) :: diag
538
539 IF (PRESENT(r_shake)) r_shake = 0.0_dp
540 IF (PRESENT(v_shake)) v_shake = 0.0_dp
541 diag = 0.0_dp
542
543 SELECT CASE (char)
544 CASE ('SHAKE')
545 IF (PRESENT(u) .AND. PRESENT(vector_v) .AND. &
546 PRESENT(vector_r)) THEN
547 diag(1, 1) = vector_r(1)
548 diag(2, 2) = vector_r(2)
549 diag(3, 3) = vector_r(3)
550 r_shake = matmul(matmul(u, diag), transpose(u))
551 diag(1, 1) = vector_v(1)
552 diag(2, 2) = vector_v(2)
553 diag(3, 3) = vector_v(3)
554 v_shake = matmul(matmul(u, diag), transpose(u))
555 diag = matmul(r_shake, v_shake)
556 r_shake = diag
557 ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. &
558 PRESENT(vector_r)) THEN
559 DO i = 1, 3
560 r_shake(i, i) = vector_r(i)*vector_v(i)
561 v_shake(i, i) = vector_v(i)
562 END DO
563 ELSE
564 cpabort("Not sufficient parameters")
565 END IF
566 CASE ('RATTLE')
567 IF (PRESENT(u) .AND. PRESENT(vector_v)) THEN
568 diag(1, 1) = vector_v(1)
569 diag(2, 2) = vector_v(2)
570 diag(3, 3) = vector_v(3)
571 v_shake = matmul(matmul(u, diag), transpose(u))
572 ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN
573 DO i = 1, 3
574 v_shake(i, i) = vector_v(i)
575 END DO
576 ELSE
577 cpabort("Not sufficient parameters")
578 END IF
579 END SELECT
580
581 END SUBROUTINE get_roll_matrix
582
583! **************************************************************************************************
584!> \brief ...
585!> \param particle_set ...
586!> \param local_particles ...
587!> \param pos ...
588!> \param vel ...
589!> \par History
590!> Teodoro Laino [tlaino] 2007
591! **************************************************************************************************
592 SUBROUTINE restore_temporary_set(particle_set, local_particles, pos, vel)
593 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
594 TYPE(distribution_1d_type), POINTER :: local_particles
595 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
596 OPTIONAL :: pos, vel
597
598 INTEGER :: iparticle, iparticle_kind, &
599 iparticle_local, nparticle_local
600 LOGICAL, ALLOCATABLE, DIMENSION(:) :: wrk
601
602 ALLOCATE (wrk(SIZE(particle_set)))
603 wrk = .true.
604 DO iparticle_kind = 1, SIZE(local_particles%n_el)
605 nparticle_local = local_particles%n_el(iparticle_kind)
606 DO iparticle_local = 1, nparticle_local
607 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
608 wrk(iparticle) = .false.
609 END DO
610 END DO
611 IF (PRESENT(vel)) THEN
612 DO iparticle = 1, SIZE(particle_set)
613 IF (wrk(iparticle)) THEN
614 vel(:, iparticle) = 0.0_dp
615 END IF
616 END DO
617 END IF
618 IF (PRESENT(pos)) THEN
619 DO iparticle = 1, SIZE(particle_set)
620 IF (wrk(iparticle)) THEN
621 pos(:, iparticle) = 0.0_dp
622 END IF
623 END DO
624 END IF
625 DEALLOCATE (wrk)
626 END SUBROUTINE restore_temporary_set
627
628! **************************************************************************************************
629!> \brief ...
630!> \param group ...
631!> \param pos ...
632!> \param vel ...
633!> \par History
634!> Teodoro Laino [tlaino] 2007
635! **************************************************************************************************
636 SUBROUTINE update_temporary_set(group, pos, vel)
637 CLASS(mp_comm_type), INTENT(IN) :: group
638 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
639 OPTIONAL :: pos, vel
640
641 IF (PRESENT(pos)) THEN
642 CALL group%sum(pos)
643 END IF
644 IF (PRESENT(vel)) THEN
645 CALL group%sum(vel)
646 END IF
647 END SUBROUTINE update_temporary_set
648
649END MODULE constraint_util
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
colvar_methods
defines collective variables s({R}) and the derivative of this variable wrt R these can then be used ...
Definition colvar_methods.F:16
colvar_methods::colvar_eval_mol_f
subroutine, public colvar_eval_mol_f(colvar, cell, particles, pos, fixd_list)
evaluates the derivatives (dsdr) given and due to the given colvar variables in a molecular environme...
Definition colvar_methods.F:1631
constraint_util
Contains routines useful for the application of constraints during MD.
Definition constraint_util.F:13
constraint_util::pv_constraint
subroutine, public pv_constraint(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, virial, group)
...
Definition constraint_util.F:198
constraint_util::get_roll_matrix
subroutine, public get_roll_matrix(char, r_shake, v_shake, vector_r, vector_v, u)
...
Definition constraint_util.F:528
constraint_util::restore_temporary_set
subroutine, public restore_temporary_set(particle_set, local_particles, pos, vel)
...
Definition constraint_util.F:593
constraint_util::update_temporary_set
subroutine, public update_temporary_set(group, pos, vel)
...
Definition constraint_util.F:637
constraint_util::getold
subroutine, public getold(gci, local_molecules, molecule_set, molecule_kind_set, particle_set, cell)
saves all of the old variables
Definition constraint_util.F:63
constraint_util::check_tol
subroutine, public check_tol(roll_tol, iroll, char, matrix, veps)
...
Definition constraint_util.F:475
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
virial_types
Definition virial_types.F:13
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
message_passing::mp_comm_type
Definition message_passing.F:189
molecule_kind_types::colvar_constraint_type
Definition molecule_kind_types.F:117
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:149
molecule_kind_types::g3x3_constraint_type
Definition molecule_kind_types.F:128
molecule_kind_types::g4x6_constraint_type
Definition molecule_kind_types.F:135
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::local_colvar_constraint_type
Definition molecule_types.F:41
molecule_types::local_constraint_type
Definition molecule_types.F:86
molecule_types::local_g3x3_constraint_type
Definition molecule_types.F:49
molecule_types::local_g4x6_constraint_type
Definition molecule_types.F:63
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
virial_types::virial_type
Definition virial_types.F:26