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 : !> \brief Contains routines useful for the application of constraints during MD
10 : !> \par History
11 : !> none
12 : ! **************************************************************************************************
13 : MODULE 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 :
48 : CONTAINS
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 16934 : 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 16934 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
74 16934 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
75 16934 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
76 16934 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
77 16934 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
78 : TYPE(local_constraint_type), POINTER :: lci
79 16934 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
80 16934 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
81 : TYPE(molecule_kind_type), POINTER :: molecule_kind
82 : TYPE(molecule_type), POINTER :: molecule
83 :
84 16934 : NULLIFY (fixd_list)
85 16934 : nkind = SIZE(molecule_kind_set)
86 : ! Intramolecular constraints
87 84474 : MOL: DO ikind = 1, nkind
88 67540 : nmol_per_kind = local_molecules%n_el(ikind)
89 396005 : DO imol = 1, nmol_per_kind
90 311531 : i = local_molecules%list(ikind)%array(imol)
91 311531 : molecule => molecule_set(i)
92 311531 : 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 311531 : fixd_list=fixd_list)
96 311531 : CALL get_molecule(molecule, lci=lci)
97 311531 : IF (.NOT. ASSOCIATED(lci)) CYCLE
98 : CALL get_molecule(molecule, first_atom=first_atom, &
99 311531 : lcolv=lcolv, lg3x3=lg3x3, lg4x6=lg4x6)
100 : CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
101 690602 : lcolv, lg3x3, lg4x6, first_atom, particle_set, cell)
102 : END DO
103 : END DO MOL
104 : ! Intermolecular constraints
105 16934 : IF (gci%ntot > 0) THEN
106 1150 : n3x3con = gci%ng3x3
107 1150 : n4x6con = gci%ng4x6
108 1150 : colv_list => gci%colv_list
109 1150 : g3x3_list => gci%g3x3_list
110 1150 : g4x6_list => gci%g4x6_list
111 1150 : fixd_list => gci%fixd_list
112 1150 : lcolv => gci%lcolv
113 1150 : lg3x3 => gci%lg3x3
114 1150 : lg4x6 => gci%lg4x6
115 : CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
116 1150 : lcolv, lg3x3, lg4x6, 1, particle_set, cell)
117 : END IF
118 16934 : 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 312681 : 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 312681 : IF (ASSOCIATED(colv_list)) THEN
155 : ! Collective constraints
156 43227 : DO iconst = 1, SIZE(colv_list)
157 : CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
158 43227 : particles=particle_set, fixd_list=fixd_list)
159 : END DO
160 : END IF
161 : ! 3x3 constraints
162 494908 : DO iconst = 1, n3x3con
163 182227 : index = g3x3_list(iconst)%a + first_atom - 1
164 1457816 : lg3x3(iconst)%ra_old = particle_set(index)%r
165 182227 : index = g3x3_list(iconst)%b + first_atom - 1
166 1457816 : lg3x3(iconst)%rb_old = particle_set(index)%r
167 182227 : index = g3x3_list(iconst)%c + first_atom - 1
168 1770497 : lg3x3(iconst)%rc_old = particle_set(index)%r
169 : END DO
170 : ! 4x6 constraints
171 313708 : DO iconst = 1, n4x6con
172 1027 : index = g4x6_list(iconst)%a + first_atom - 1
173 8216 : lg4x6(iconst)%ra_old = particle_set(index)%r
174 1027 : index = g4x6_list(iconst)%b + first_atom - 1
175 8216 : lg4x6(iconst)%rb_old = particle_set(index)%r
176 1027 : index = g4x6_list(iconst)%c + first_atom - 1
177 8216 : lg4x6(iconst)%rc_old = particle_set(index)%r
178 1027 : index = g4x6_list(iconst)%d + first_atom - 1
179 320897 : lg4x6(iconst)%rd_old = particle_set(index)%r
180 : END DO
181 :
182 312681 : 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 21178 : 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 21178 : TYPE(colvar_constraint_type), POINTER :: colv_list(:)
212 21178 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
213 21178 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
214 21178 : TYPE(local_colvar_constraint_type), POINTER :: lcolv(:)
215 : TYPE(local_constraint_type), POINTER :: lci
216 21178 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
217 21178 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
218 : TYPE(molecule_kind_type), POINTER :: molecule_kind
219 : TYPE(molecule_type), POINTER :: molecule
220 :
221 21178 : pv = 0.0_dp
222 21178 : nkind = SIZE(molecule_kind_set)
223 : ! Intramolecular Constraints
224 92988 : MOL: DO ikind = 1, nkind
225 71810 : nmol_per_kind = local_molecules%n_el(ikind)
226 660848 : DO imol = 1, nmol_per_kind
227 567860 : i = local_molecules%list(ikind)%array(imol)
228 567860 : molecule => molecule_set(i)
229 567860 : 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 567860 : colv_list=colv_list)
233 567860 : CALL get_molecule(molecule, lci=lci)
234 567860 : IF (.NOT. ASSOCIATED(lci)) CYCLE
235 : CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3, &
236 567860 : lg4x6=lg4x6, lcolv=lcolv)
237 : CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
238 1207530 : first_atom, lg3x3, lg4x6, lcolv, particle_set, pv)
239 : END DO
240 : END DO MOL
241 : ! Intermolecular constraints
242 21178 : IF (gci%ntot > 0) THEN
243 1120 : ng3x3 = gci%ng3x3
244 1120 : ng4x6 = gci%ng4x6
245 1120 : colv_list => gci%colv_list
246 1120 : g3x3_list => gci%g3x3_list
247 1120 : g4x6_list => gci%g4x6_list
248 1120 : lcolv => gci%lcolv
249 1120 : lg3x3 => gci%lg3x3
250 1120 : lg4x6 => gci%lg4x6
251 : CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
252 1120 : 1, lg3x3, lg4x6, lcolv, particle_set, pv)
253 : END IF
254 21178 : CALL group%sum(pv)
255 : ! Symmetrize PV
256 21178 : pv(1, 2) = 0.5_dp*(pv(1, 2) + pv(2, 1))
257 21178 : pv(2, 1) = pv(1, 2)
258 21178 : pv(1, 3) = 0.5_dp*(pv(1, 3) + pv(3, 1))
259 21178 : pv(3, 1) = pv(1, 3)
260 21178 : pv(3, 2) = 0.5_dp*(pv(3, 2) + pv(2, 3))
261 21178 : pv(2, 3) = pv(3, 2)
262 : ! Store in virial type
263 275314 : virial%pv_constraint = pv
264 :
265 21178 : 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 568980 : 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 568980 : IF (ASSOCIATED(colv_list)) THEN
303 : ! Colvar Constraints
304 62284 : DO iconst = 1, SIZE(colv_list)
305 62284 : CALL pv_colv_eval(pv, lcolv(iconst), particle_set)
306 : END DO
307 : END IF
308 : ! 3x3
309 1000052 : DO iconst = 1, ng3x3
310 : ! pv gets updated with FULL multiplier
311 1724288 : lambda_3x3 = lg3x3(iconst)%lambda
312 :
313 : fc1 = lambda_3x3(1)*lg3x3(iconst)%fa + &
314 1724288 : lambda_3x3(2)*lg3x3(iconst)%fb
315 : fc2 = -lambda_3x3(1)*lg3x3(iconst)%fa + &
316 1724288 : lambda_3x3(3)*lg3x3(iconst)%fc
317 : fc3 = -lambda_3x3(2)*lg3x3(iconst)%fb - &
318 1724288 : lambda_3x3(3)*lg3x3(iconst)%fc
319 431072 : index_a = g3x3_list(iconst)%a + first_atom - 1
320 431072 : index_b = g3x3_list(iconst)%b + first_atom - 1
321 431072 : index_c = g3x3_list(iconst)%c + first_atom - 1
322 :
323 : !pv(1,1)
324 431072 : pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
325 431072 : pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
326 431072 : pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
327 : !pv(1,2)
328 431072 : pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
329 431072 : pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
330 431072 : pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
331 : !pv(1,3)
332 431072 : pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
333 431072 : pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
334 431072 : pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
335 : !pv(2,1)
336 431072 : pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
337 431072 : pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
338 431072 : pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
339 : !pv(2,2)
340 431072 : pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
341 431072 : pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
342 431072 : pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
343 : !pv(2,3)
344 431072 : pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
345 431072 : pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
346 431072 : pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
347 : !pv(3,1)
348 431072 : pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
349 431072 : pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
350 431072 : pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
351 : !pv(3,2)
352 431072 : pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
353 431072 : pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
354 431072 : pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
355 : !pv(3,3)
356 431072 : pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
357 431072 : pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
358 1000052 : pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
359 : END DO
360 :
361 : ! 4x6
362 572437 : DO iconst = 1, ng4x6
363 : ! pv gets updated with FULL multiplier
364 24199 : lambda_4x6 = lg4x6(iconst)%lambda
365 :
366 : fc1 = lambda_4x6(1)*lg4x6(iconst)%fa + &
367 : lambda_4x6(2)*lg4x6(iconst)%fb + &
368 13828 : lambda_4x6(3)*lg4x6(iconst)%fc
369 : fc2 = -lambda_4x6(1)*lg4x6(iconst)%fa + &
370 : lambda_4x6(4)*lg4x6(iconst)%fd + &
371 13828 : lambda_4x6(5)*lg4x6(iconst)%fe
372 : fc3 = -lambda_4x6(2)*lg4x6(iconst)%fb - &
373 : lambda_4x6(4)*lg4x6(iconst)%fd + &
374 13828 : lambda_4x6(6)*lg4x6(iconst)%ff
375 : fc4 = -lambda_4x6(3)*lg4x6(iconst)%fc - &
376 : lambda_4x6(5)*lg4x6(iconst)%fe - &
377 13828 : lambda_4x6(6)*lg4x6(iconst)%ff
378 3457 : index_a = g4x6_list(iconst)%a + first_atom - 1
379 3457 : index_b = g4x6_list(iconst)%b + first_atom - 1
380 3457 : index_c = g4x6_list(iconst)%c + first_atom - 1
381 3457 : index_d = g4x6_list(iconst)%d + first_atom - 1
382 :
383 : !pv(1,1)
384 3457 : pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
385 3457 : pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
386 3457 : pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
387 3457 : pv(1, 1) = pv(1, 1) + fc4(1)*particle_set(index_d)%r(1)
388 : !pv(1,2)
389 3457 : pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
390 3457 : pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
391 3457 : pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
392 3457 : pv(1, 2) = pv(1, 2) + fc4(1)*particle_set(index_d)%r(2)
393 : !pv(1,3)
394 3457 : pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
395 3457 : pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
396 3457 : pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
397 3457 : pv(1, 3) = pv(1, 3) + fc4(1)*particle_set(index_d)%r(3)
398 : !pv(2,1)
399 3457 : pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
400 3457 : pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
401 3457 : pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
402 3457 : pv(2, 1) = pv(2, 1) + fc4(2)*particle_set(index_d)%r(1)
403 : !pv(2,2)
404 3457 : pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
405 3457 : pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
406 3457 : pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
407 3457 : pv(2, 2) = pv(2, 2) + fc4(2)*particle_set(index_d)%r(2)
408 : !pv(2,3)
409 3457 : pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
410 3457 : pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
411 3457 : pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
412 3457 : pv(2, 3) = pv(2, 3) + fc4(2)*particle_set(index_d)%r(3)
413 : !pv(3,1)
414 3457 : pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
415 3457 : pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
416 3457 : pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
417 3457 : pv(3, 1) = pv(3, 1) + fc4(3)*particle_set(index_d)%r(1)
418 : !pv(3,2)
419 3457 : pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
420 3457 : pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
421 3457 : pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
422 3457 : pv(3, 2) = pv(3, 2) + fc4(3)*particle_set(index_d)%r(2)
423 : !pv(3,3)
424 3457 : pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
425 3457 : pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
426 3457 : pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
427 572437 : pv(3, 3) = pv(3, 3) + fc4(3)*particle_set(index_d)%r(3)
428 : END DO
429 :
430 568980 : 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 25701 : 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 81829 : DO iatm = 1, SIZE(lcolv%colvar_old%i_atom)
450 56128 : ind = lcolv%colvar_old%i_atom(iatm)
451 224512 : f = -lcolv%colvar_old%dsdr(:, iatm)
452 : ! pv gets updated with FULL multiplier
453 56128 : lambda = lcolv%lambda
454 250213 : DO i = 1, 3
455 168384 : tmp = lambda*particle_set(ind)%r(i)
456 729664 : DO j = 1, 3
457 673536 : pv(j, i) = pv(j, i) + f(j)*tmp
458 : END DO
459 : END DO
460 : END DO
461 :
462 25701 : 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 3620 : 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 1826 : SELECT CASE (char)
487 : CASE ('SHAKE')
488 1826 : IF (iroll == 1) THEN
489 8814 : matrix_old = matrix
490 678 : roll_tol = -1.E10_dp
491 : ELSE
492 : roll_tol = 0.0_dp
493 14924 : diff_shake = ABS(matrix_old - matrix)
494 14924 : local_tol = MAXVAL(diff_shake)
495 1148 : roll_tol = MAX(roll_tol, local_tol)
496 14924 : matrix_old = matrix
497 : END IF
498 1826 : iroll = iroll + 1
499 : CASE ('RATTLE')
500 1794 : IF (iroll == 1) THEN
501 8814 : veps_old = veps
502 678 : roll_tol = -1.E+10_dp
503 : ELSE
504 : roll_tol = 0.0_dp
505 : ! compute tolerance on veps
506 14508 : diff_rattle = ABS(veps - veps_old)
507 14508 : local_tol = MAXVAL(diff_rattle)
508 1116 : roll_tol = MAX(roll_tol, local_tol)
509 14508 : veps_old = veps
510 : END IF
511 5414 : iroll = iroll + 1
512 : END SELECT
513 :
514 3620 : 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 3620 : 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 25532 : IF (PRESENT(r_shake)) r_shake = 0.0_dp
540 47060 : IF (PRESENT(v_shake)) v_shake = 0.0_dp
541 3620 : diag = 0.0_dp
542 :
543 1826 : SELECT CASE (char)
544 : CASE ('SHAKE')
545 1826 : IF (PRESENT(u) .AND. PRESENT(vector_v) .AND. &
546 1794 : PRESENT(vector_r)) THEN
547 20 : diag(1, 1) = vector_r(1)
548 20 : diag(2, 2) = vector_r(2)
549 20 : diag(3, 3) = vector_r(3)
550 2660 : r_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
551 20 : diag(1, 1) = vector_v(1)
552 20 : diag(2, 2) = vector_v(2)
553 20 : diag(3, 3) = vector_v(3)
554 2660 : v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
555 800 : diag = MATMUL(r_shake, v_shake)
556 260 : r_shake = diag
557 1806 : ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. &
558 : PRESENT(vector_r)) THEN
559 7224 : DO i = 1, 3
560 5418 : r_shake(i, i) = vector_r(i)*vector_v(i)
561 7224 : v_shake(i, i) = vector_v(i)
562 : END DO
563 : ELSE
564 0 : CPABORT("Not sufficient parameters")
565 : END IF
566 : CASE ('RATTLE')
567 3620 : IF (PRESENT(u) .AND. PRESENT(vector_v)) THEN
568 20 : diag(1, 1) = vector_v(1)
569 20 : diag(2, 2) = vector_v(2)
570 20 : diag(3, 3) = vector_v(3)
571 2660 : v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
572 1774 : ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN
573 7096 : DO i = 1, 3
574 7096 : v_shake(i, i) = vector_v(i)
575 : END DO
576 : ELSE
577 0 : CPABORT("Not sufficient parameters")
578 : END IF
579 : END SELECT
580 :
581 3620 : 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 3345 : 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 3345 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: wrk
601 :
602 10035 : ALLOCATE (wrk(SIZE(particle_set)))
603 187160 : wrk = .TRUE.
604 8555 : DO iparticle_kind = 1, SIZE(local_particles%n_el)
605 5210 : nparticle_local = local_particles%n_el(iparticle_kind)
606 100608 : DO iparticle_local = 1, nparticle_local
607 92053 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
608 97263 : wrk(iparticle) = .FALSE.
609 : END DO
610 : END DO
611 3345 : IF (PRESENT(vel)) THEN
612 187160 : DO iparticle = 1, SIZE(particle_set)
613 187160 : IF (wrk(iparticle)) THEN
614 367048 : vel(:, iparticle) = 0.0_dp
615 : END IF
616 : END DO
617 : END IF
618 3345 : IF (PRESENT(pos)) THEN
619 98242 : DO iparticle = 1, SIZE(particle_set)
620 98242 : IF (wrk(iparticle)) THEN
621 192504 : pos(:, iparticle) = 0.0_dp
622 : END IF
623 : END DO
624 : END IF
625 3345 : DEALLOCATE (wrk)
626 3345 : 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 3345 : 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 3345 : IF (PRESENT(pos)) THEN
642 773035 : CALL group%sum(pos)
643 : END IF
644 3345 : IF (PRESENT(vel)) THEN
645 1473865 : CALL group%sum(vel)
646 : END IF
647 3345 : END SUBROUTINE update_temporary_set
648 :
649 60 : END MODULE constraint_util
|
