Line data Source code
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 : !> \par History
10 : !> none
11 : ! **************************************************************************************************
12 : MODULE constraint_3x3
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE constraint_fxd, ONLY: check_fixed_atom_cns_g3x3
17 : USE kinds, ONLY: dp
18 : USE linear_systems, ONLY: solve_system
19 : USE molecule_kind_types, ONLY: fixd_constraint_type,&
20 : g3x3_constraint_type,&
21 : get_molecule_kind,&
22 : molecule_kind_type
23 : USE molecule_types, ONLY: get_molecule,&
24 : global_constraint_type,&
25 : local_g3x3_constraint_type,&
26 : molecule_type
27 : USE particle_types, ONLY: particle_type
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 : PUBLIC :: shake_3x3_int, &
34 : rattle_3x3_int, &
35 : shake_roll_3x3_int, &
36 : rattle_roll_3x3_int, &
37 : shake_3x3_ext, &
38 : rattle_3x3_ext, &
39 : shake_roll_3x3_ext, &
40 : rattle_roll_3x3_ext
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_3x3'
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param molecule ...
49 : !> \param particle_set ...
50 : !> \param pos ...
51 : !> \param vel ...
52 : !> \param dt ...
53 : !> \param ishake ...
54 : !> \param max_sigma ...
55 : !> \par History
56 : !> none
57 : ! **************************************************************************************************
58 281363 : SUBROUTINE shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake, &
59 : max_sigma)
60 :
61 : TYPE(molecule_type), POINTER :: molecule
62 : TYPE(particle_type), POINTER :: particle_set(:)
63 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
64 : REAL(kind=dp), INTENT(in) :: dt
65 : INTEGER, INTENT(IN) :: ishake
66 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
67 :
68 : INTEGER :: first_atom, ng3x3
69 281363 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
70 281363 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
71 281363 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
72 : TYPE(molecule_kind_type), POINTER :: molecule_kind
73 :
74 281363 : NULLIFY (fixd_list)
75 281363 : molecule_kind => molecule%molecule_kind
76 : CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
77 281363 : g3x3_list=g3x3_list, fixd_list=fixd_list)
78 281363 : CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3)
79 : ! Real Shake
80 : CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
81 281363 : particle_set, pos, vel, dt, ishake, max_sigma)
82 :
83 281363 : END SUBROUTINE shake_3x3_int
84 :
85 : ! **************************************************************************************************
86 : !> \brief ...
87 : !> \param molecule ...
88 : !> \param particle_set ...
89 : !> \param pos ...
90 : !> \param vel ...
91 : !> \param r_shake ...
92 : !> \param v_shake ...
93 : !> \param dt ...
94 : !> \param ishake ...
95 : !> \param max_sigma ...
96 : !> \par History
97 : !> none
98 : ! **************************************************************************************************
99 128355 : SUBROUTINE shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, &
100 128355 : v_shake, dt, ishake, max_sigma)
101 :
102 : TYPE(molecule_type), POINTER :: molecule
103 : TYPE(particle_type), POINTER :: particle_set(:)
104 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
105 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
106 : REAL(kind=dp), INTENT(in) :: dt
107 : INTEGER, INTENT(IN) :: ishake
108 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
109 :
110 : INTEGER :: first_atom, ng3x3
111 128355 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
112 128355 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
113 128355 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
114 : TYPE(molecule_kind_type), POINTER :: molecule_kind
115 :
116 128355 : NULLIFY (fixd_list)
117 128355 : molecule_kind => molecule%molecule_kind
118 : CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
119 128355 : g3x3_list=g3x3_list, fixd_list=fixd_list)
120 128355 : CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3)
121 : ! Real Shake
122 : CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
123 128355 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
124 :
125 128355 : END SUBROUTINE shake_roll_3x3_int
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param molecule ...
130 : !> \param particle_set ...
131 : !> \param vel ...
132 : !> \param r_rattle ...
133 : !> \param dt ...
134 : !> \param veps ...
135 : !> \par History
136 : !> none
137 : ! **************************************************************************************************
138 97999 : SUBROUTINE rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, veps)
139 :
140 : TYPE(molecule_type), POINTER :: molecule
141 : TYPE(particle_type), POINTER :: particle_set(:)
142 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
143 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
144 : REAL(kind=dp), INTENT(in) :: dt
145 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
146 :
147 : INTEGER :: first_atom
148 97999 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
149 97999 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
150 97999 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
151 : TYPE(molecule_kind_type), POINTER :: molecule_kind
152 :
153 97999 : NULLIFY (fixd_list)
154 97999 : molecule_kind => molecule%molecule_kind
155 : CALL get_molecule_kind(molecule_kind, &
156 97999 : g3x3_list=g3x3_list, fixd_list=fixd_list)
157 97999 : CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3)
158 : ! Real Rattle
159 : CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
160 97999 : particle_set, vel, r_rattle, dt, veps)
161 :
162 97999 : END SUBROUTINE rattle_roll_3x3_int
163 :
164 : ! **************************************************************************************************
165 : !> \brief ...
166 : !> \param molecule ...
167 : !> \param particle_set ...
168 : !> \param vel ...
169 : !> \param dt ...
170 : !> \par History
171 : !> none
172 : ! **************************************************************************************************
173 143796 : SUBROUTINE rattle_3x3_int(molecule, particle_set, vel, dt)
174 :
175 : TYPE(molecule_type), POINTER :: molecule
176 : TYPE(particle_type), POINTER :: particle_set(:)
177 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
178 : REAL(kind=dp), INTENT(in) :: dt
179 :
180 : INTEGER :: first_atom
181 143796 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
182 143796 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
183 143796 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
184 : TYPE(molecule_kind_type), POINTER :: molecule_kind
185 :
186 143796 : NULLIFY (fixd_list)
187 143796 : molecule_kind => molecule%molecule_kind
188 : CALL get_molecule_kind(molecule_kind, &
189 143796 : g3x3_list=g3x3_list, fixd_list=fixd_list)
190 143796 : CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3)
191 : ! Real Rattle
192 : CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
193 143796 : particle_set, vel, dt)
194 :
195 143796 : END SUBROUTINE rattle_3x3_int
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param gci ...
200 : !> \param particle_set ...
201 : !> \param pos ...
202 : !> \param vel ...
203 : !> \param dt ...
204 : !> \param ishake ...
205 : !> \param max_sigma ...
206 : !> \par History
207 : !> none
208 : ! **************************************************************************************************
209 76 : SUBROUTINE shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake, &
210 : max_sigma)
211 :
212 : TYPE(global_constraint_type), POINTER :: gci
213 : TYPE(particle_type), POINTER :: particle_set(:)
214 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
215 : REAL(kind=dp), INTENT(in) :: dt
216 : INTEGER, INTENT(IN) :: ishake
217 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
218 :
219 : INTEGER :: first_atom, ng3x3
220 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
221 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
222 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
223 :
224 76 : first_atom = 1
225 76 : ng3x3 = gci%ng3x3
226 76 : g3x3_list => gci%g3x3_list
227 76 : fixd_list => gci%fixd_list
228 76 : lg3x3 => gci%lg3x3
229 : ! Real Shake
230 : CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
231 76 : particle_set, pos, vel, dt, ishake, max_sigma)
232 :
233 76 : END SUBROUTINE shake_3x3_ext
234 :
235 : ! **************************************************************************************************
236 : !> \brief ...
237 : !> \param gci ...
238 : !> \param particle_set ...
239 : !> \param pos ...
240 : !> \param vel ...
241 : !> \param r_shake ...
242 : !> \param v_shake ...
243 : !> \param dt ...
244 : !> \param ishake ...
245 : !> \param max_sigma ...
246 : !> \par History
247 : !> none
248 : ! **************************************************************************************************
249 0 : SUBROUTINE shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, &
250 0 : v_shake, dt, ishake, max_sigma)
251 :
252 : TYPE(global_constraint_type), POINTER :: gci
253 : TYPE(particle_type), POINTER :: particle_set(:)
254 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
255 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
256 : REAL(kind=dp), INTENT(in) :: dt
257 : INTEGER, INTENT(IN) :: ishake
258 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
259 :
260 : INTEGER :: first_atom, ng3x3
261 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
262 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
263 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
264 :
265 0 : first_atom = 1
266 0 : ng3x3 = gci%ng3x3
267 0 : g3x3_list => gci%g3x3_list
268 0 : fixd_list => gci%fixd_list
269 0 : lg3x3 => gci%lg3x3
270 : ! Real Shake
271 : CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
272 0 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
273 :
274 0 : END SUBROUTINE shake_roll_3x3_ext
275 :
276 : ! **************************************************************************************************
277 : !> \brief ...
278 : !> \param gci ...
279 : !> \param particle_set ...
280 : !> \param vel ...
281 : !> \param r_rattle ...
282 : !> \param dt ...
283 : !> \param veps ...
284 : !> \par History
285 : !> none
286 : ! **************************************************************************************************
287 0 : SUBROUTINE rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, veps)
288 :
289 : TYPE(global_constraint_type), POINTER :: gci
290 : TYPE(particle_type), POINTER :: particle_set(:)
291 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
292 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
293 : REAL(kind=dp), INTENT(in) :: dt
294 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
295 :
296 : INTEGER :: first_atom
297 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
298 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
299 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
300 :
301 0 : first_atom = 1
302 0 : g3x3_list => gci%g3x3_list
303 0 : fixd_list => gci%fixd_list
304 0 : lg3x3 => gci%lg3x3
305 : ! Real Rattle
306 : CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
307 0 : particle_set, vel, r_rattle, dt, veps)
308 :
309 0 : END SUBROUTINE rattle_roll_3x3_ext
310 :
311 : ! **************************************************************************************************
312 : !> \brief ...
313 : !> \param gci ...
314 : !> \param particle_set ...
315 : !> \param vel ...
316 : !> \param dt ...
317 : !> \par History
318 : !> none
319 : ! **************************************************************************************************
320 40 : SUBROUTINE rattle_3x3_ext(gci, particle_set, vel, dt)
321 :
322 : TYPE(global_constraint_type), POINTER :: gci
323 : TYPE(particle_type), POINTER :: particle_set(:)
324 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
325 : REAL(kind=dp), INTENT(in) :: dt
326 :
327 : INTEGER :: first_atom
328 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
329 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
330 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
331 :
332 40 : first_atom = 1
333 40 : g3x3_list => gci%g3x3_list
334 40 : fixd_list => gci%fixd_list
335 40 : lg3x3 => gci%lg3x3
336 : ! Real Rattle
337 : CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
338 40 : particle_set, vel, dt)
339 :
340 40 : END SUBROUTINE rattle_3x3_ext
341 :
342 : ! **************************************************************************************************
343 : !> \brief ...
344 : !> \param fixd_list ...
345 : !> \param g3x3_list ...
346 : !> \param lg3x3 ...
347 : !> \param first_atom ...
348 : !> \param ng3x3 ...
349 : !> \param particle_set ...
350 : !> \param pos ...
351 : !> \param vel ...
352 : !> \param dt ...
353 : !> \param ishake ...
354 : !> \param max_sigma ...
355 : !> \par History
356 : !> none
357 : ! **************************************************************************************************
358 281439 : SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
359 281439 : particle_set, pos, vel, dt, ishake, max_sigma)
360 :
361 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
362 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
363 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
364 : INTEGER, INTENT(IN) :: first_atom, ng3x3
365 : TYPE(particle_type), POINTER :: particle_set(:)
366 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
367 : REAL(kind=dp), INTENT(in) :: dt
368 : INTEGER, INTENT(IN) :: ishake
369 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
370 :
371 : INTEGER :: iconst, index_a, index_b, index_c
372 : REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, &
373 : imass3, sigma
374 : REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, r0_12, r0_13, r0_23, r1, &
375 : r12, r13, r2, r23, r3, v1, v2, v3, vec
376 : REAL(KIND=dp), DIMENSION(3, 1) :: bvec
377 : REAL(KIND=dp), DIMENSION(3, 3) :: amat, atemp
378 : TYPE(atomic_kind_type), POINTER :: atomic_kind
379 :
380 281439 : dtsqby2 = dt*dt*.5_dp
381 281439 : idtsq = 1.0_dp/dt/dt
382 281439 : dtby2 = dt*.5_dp
383 562878 : DO iconst = 1, ng3x3
384 281439 : IF (g3x3_list(iconst)%restraint%active) CYCLE
385 281419 : index_a = g3x3_list(iconst)%a + first_atom - 1
386 281419 : index_b = g3x3_list(iconst)%b + first_atom - 1
387 281419 : index_c = g3x3_list(iconst)%c + first_atom - 1
388 281419 : IF (ishake == 1) THEN
389 575264 : r0_12(:) = pos(:, index_a) - pos(:, index_b)
390 575264 : r0_13(:) = pos(:, index_a) - pos(:, index_c)
391 575264 : r0_23(:) = pos(:, index_b) - pos(:, index_c)
392 143816 : atomic_kind => particle_set(index_a)%atomic_kind
393 143816 : imass1 = 1.0_dp/atomic_kind%mass
394 143816 : atomic_kind => particle_set(index_b)%atomic_kind
395 143816 : imass2 = 1.0_dp/atomic_kind%mass
396 143816 : atomic_kind => particle_set(index_c)%atomic_kind
397 143816 : imass3 = 1.0_dp/atomic_kind%mass
398 : lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
399 575264 : lg3x3(iconst)%rb_old)
400 : lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
401 575264 : lg3x3(iconst)%rc_old)
402 : lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
403 575264 : lg3x3(iconst)%rc_old)
404 : ! Check for fixed atom constraints
405 : CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
406 143816 : index_a, index_b, index_c, fixd_list, lg3x3(iconst))
407 : ! construct matrix
408 575264 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, lg3x3(iconst)%fa)
409 575264 : amat(1, 2) = imass1*DOT_PRODUCT(r0_12, lg3x3(iconst)%fb)
410 575264 : amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, lg3x3(iconst)%fc)
411 575264 : amat(2, 1) = imass1*DOT_PRODUCT(r0_13, lg3x3(iconst)%fa)
412 575264 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, lg3x3(iconst)%fb)
413 575264 : amat(2, 3) = imass3*DOT_PRODUCT(r0_13, lg3x3(iconst)%fc)
414 575264 : amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, lg3x3(iconst)%fa)
415 575264 : amat(3, 2) = imass3*DOT_PRODUCT(r0_23, lg3x3(iconst)%fb)
416 575264 : amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, lg3x3(iconst)%fc)
417 : ! Store values
418 575264 : lg3x3(iconst)%r0_12 = r0_12
419 575264 : lg3x3(iconst)%r0_13 = r0_13
420 575264 : lg3x3(iconst)%r0_23 = r0_23
421 1869608 : lg3x3(iconst)%amat = amat
422 575264 : lg3x3(iconst)%lambda_old = 0.0_dp
423 575264 : lg3x3(iconst)%del_lambda = 0.0_dp
424 : ELSE
425 : ! Retrieve values
426 550412 : r0_12 = lg3x3(iconst)%r0_12
427 550412 : r0_13 = lg3x3(iconst)%r0_13
428 550412 : r0_23 = lg3x3(iconst)%r0_23
429 1788839 : amat = lg3x3(iconst)%amat
430 137603 : imass1 = lg3x3(iconst)%imass1
431 137603 : imass2 = lg3x3(iconst)%imass2
432 137603 : imass3 = lg3x3(iconst)%imass3
433 : END IF
434 : ! Iterate until convergence
435 : vec = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*(imass1 + imass2) + &
436 : lg3x3(iconst)%lambda(2)*imass1*lg3x3(iconst)%fb - &
437 1125676 : lg3x3(iconst)%lambda(3)*imass2*lg3x3(iconst)%fc
438 : bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
439 2251352 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
440 : vec = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass1 + &
441 : lg3x3(iconst)%lambda(2)*(imass1 + imass3)*lg3x3(iconst)%fb + &
442 1125676 : lg3x3(iconst)%lambda(3)*imass3*lg3x3(iconst)%fc
443 : bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
444 2251352 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
445 : vec = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa*imass2 + &
446 : lg3x3(iconst)%lambda(2)*imass3*lg3x3(iconst)%fb + &
447 1125676 : lg3x3(iconst)%lambda(3)*(imass2 + imass3)*lg3x3(iconst)%fc
448 : bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
449 2251352 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
450 1407095 : bvec = bvec*idtsq
451 : ! get lambda
452 281419 : atemp = amat
453 281419 : CALL solve_system(atemp, 3, bvec)
454 1125676 : lg3x3(iconst)%lambda(:) = bvec(:, 1)
455 : lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
456 1125676 : lg3x3(iconst)%lambda_old(:)
457 1125676 : lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
458 : fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
459 1125676 : lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
460 : fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
461 1125676 : lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
462 : fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
463 1125676 : lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
464 1125676 : r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
465 1125676 : r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
466 1125676 : r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
467 1125676 : v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
468 1125676 : v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
469 1125676 : v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
470 1125676 : r12 = r1 - r2
471 1125676 : r13 = r1 - r3
472 1125676 : r23 = r2 - r3
473 :
474 : ! compute the tolerance:
475 : sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* &
476 1125676 : g3x3_list(iconst)%dab
477 281419 : max_sigma = MAX(max_sigma, ABS(sigma))
478 : sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* &
479 1125676 : g3x3_list(iconst)%dac
480 281419 : max_sigma = MAX(max_sigma, ABS(sigma))
481 : sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* &
482 1125676 : g3x3_list(iconst)%dbc
483 281419 : max_sigma = MAX(max_sigma, ABS(sigma))
484 : ! update positions with full multiplier
485 1125676 : pos(:, index_a) = r1(:)
486 1125676 : pos(:, index_b) = r2(:)
487 1125676 : pos(:, index_c) = r3(:)
488 :
489 : ! update velocites with full multiplier
490 1125676 : vel(:, index_a) = v1(:)
491 1125676 : vel(:, index_b) = v2(:)
492 1407115 : vel(:, index_c) = v3(:)
493 : END DO
494 :
495 281439 : END SUBROUTINE shake_3x3_low
496 :
497 : ! **************************************************************************************************
498 : !> \brief ...
499 : !> \param fixd_list ...
500 : !> \param g3x3_list ...
501 : !> \param lg3x3 ...
502 : !> \param first_atom ...
503 : !> \param ng3x3 ...
504 : !> \param particle_set ...
505 : !> \param pos ...
506 : !> \param vel ...
507 : !> \param r_shake ...
508 : !> \param v_shake ...
509 : !> \param dt ...
510 : !> \param ishake ...
511 : !> \param max_sigma ...
512 : !> \par History
513 : !> none
514 : ! **************************************************************************************************
515 128355 : SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
516 128355 : particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
517 :
518 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
519 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
520 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
521 : INTEGER, INTENT(IN) :: first_atom, ng3x3
522 : TYPE(particle_type), POINTER :: particle_set(:)
523 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
524 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake, v_shake
525 : REAL(kind=dp), INTENT(in) :: dt
526 : INTEGER, INTENT(IN) :: ishake
527 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
528 :
529 : INTEGER :: iconst, index_a, index_b, index_c
530 : REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, &
531 : imass3, sigma
532 : REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, fc1, fc2, &
533 : fc3, r0_12, r0_13, r0_23, r1, r12, &
534 : r13, r2, r23, r3, v1, v2, v3, vec
535 : REAL(KIND=dp), DIMENSION(3, 1) :: bvec
536 : REAL(KIND=dp), DIMENSION(3, 3) :: amat, atemp
537 : TYPE(atomic_kind_type), POINTER :: atomic_kind
538 :
539 128355 : dtsqby2 = dt*dt*.5_dp
540 128355 : idtsq = 1.0_dp/dt/dt
541 128355 : dtby2 = dt*.5_dp
542 256710 : DO iconst = 1, ng3x3
543 128355 : IF (g3x3_list(iconst)%restraint%active) CYCLE
544 128355 : index_a = g3x3_list(iconst)%a + first_atom - 1
545 128355 : index_b = g3x3_list(iconst)%b + first_atom - 1
546 128355 : index_c = g3x3_list(iconst)%c + first_atom - 1
547 128355 : IF (ishake == 1) THEN
548 450712 : r0_12(:) = pos(:, index_a) - pos(:, index_b)
549 450712 : r0_13(:) = pos(:, index_a) - pos(:, index_c)
550 450712 : r0_23(:) = pos(:, index_b) - pos(:, index_c)
551 112678 : atomic_kind => particle_set(index_a)%atomic_kind
552 112678 : imass1 = 1.0_dp/atomic_kind%mass
553 112678 : atomic_kind => particle_set(index_b)%atomic_kind
554 112678 : imass2 = 1.0_dp/atomic_kind%mass
555 112678 : atomic_kind => particle_set(index_c)%atomic_kind
556 112678 : imass3 = 1.0_dp/atomic_kind%mass
557 : lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
558 450712 : lg3x3(iconst)%rb_old)
559 : lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
560 450712 : lg3x3(iconst)%rc_old)
561 : lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
562 450712 : lg3x3(iconst)%rc_old)
563 : ! Check for fixed atom constraints
564 : CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
565 112678 : index_a, index_b, index_c, fixd_list, lg3x3(iconst))
566 : ! rotate fconst:
567 1464814 : f_roll1 = MATMUL(r_shake, lg3x3(iconst)%fa)
568 1464814 : f_roll2 = MATMUL(r_shake, lg3x3(iconst)%fb)
569 1464814 : f_roll3 = MATMUL(r_shake, lg3x3(iconst)%fc)
570 : ! construct matrix
571 450712 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, f_roll1)
572 450712 : amat(1, 2) = imass1*DOT_PRODUCT(r0_12, f_roll2)
573 450712 : amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, f_roll3)
574 450712 : amat(2, 1) = imass1*DOT_PRODUCT(r0_13, f_roll1)
575 450712 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, f_roll2)
576 450712 : amat(2, 3) = imass3*DOT_PRODUCT(r0_13, f_roll3)
577 450712 : amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, f_roll1)
578 450712 : amat(3, 2) = imass3*DOT_PRODUCT(r0_23, f_roll2)
579 450712 : amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, f_roll3)
580 : ! Store values
581 450712 : lg3x3(iconst)%r0_12 = r0_12
582 450712 : lg3x3(iconst)%r0_13 = r0_13
583 450712 : lg3x3(iconst)%r0_23 = r0_23
584 450712 : lg3x3(iconst)%f_roll1 = f_roll1
585 450712 : lg3x3(iconst)%f_roll2 = f_roll2
586 450712 : lg3x3(iconst)%f_roll3 = f_roll3
587 1464814 : lg3x3(iconst)%amat = amat
588 450712 : lg3x3(iconst)%lambda_old = 0.0_dp
589 450712 : lg3x3(iconst)%del_lambda = 0.0_dp
590 : ELSE
591 : ! Retrieve values
592 62708 : r0_12 = lg3x3(iconst)%r0_12
593 62708 : r0_13 = lg3x3(iconst)%r0_13
594 62708 : r0_23 = lg3x3(iconst)%r0_23
595 62708 : f_roll1 = lg3x3(iconst)%f_roll1
596 62708 : f_roll2 = lg3x3(iconst)%f_roll2
597 62708 : f_roll3 = lg3x3(iconst)%f_roll3
598 203801 : amat = lg3x3(iconst)%amat
599 15677 : imass1 = lg3x3(iconst)%imass1
600 15677 : imass2 = lg3x3(iconst)%imass2
601 15677 : imass3 = lg3x3(iconst)%imass3
602 : END IF
603 : ! Iterate until convergence
604 : vec = lg3x3(iconst)%lambda(1)*f_roll1*(imass1 + imass2) + &
605 : lg3x3(iconst)%lambda(2)*imass1*f_roll2 - &
606 513420 : lg3x3(iconst)%lambda(3)*imass2*f_roll3
607 : bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
608 1026840 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
609 : vec = lg3x3(iconst)%lambda(1)*f_roll1*imass1 + &
610 : lg3x3(iconst)%lambda(2)*(imass1 + imass3)*f_roll2 + &
611 513420 : lg3x3(iconst)%lambda(3)*imass3*f_roll3
612 : bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
613 1026840 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
614 : vec = -lg3x3(iconst)%lambda(1)*f_roll1*imass2 + &
615 : lg3x3(iconst)%lambda(2)*imass3*f_roll2 + &
616 513420 : lg3x3(iconst)%lambda(3)*(imass2 + imass3)*f_roll3
617 : bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
618 1026840 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
619 641775 : bvec = bvec*idtsq
620 :
621 : ! get lambda
622 128355 : atemp = amat
623 128355 : CALL solve_system(atemp, 3, bvec)
624 513420 : lg3x3(iconst)%lambda(:) = bvec(:, 1)
625 : lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
626 513420 : lg3x3(iconst)%lambda_old(:)
627 513420 : lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
628 :
629 : fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
630 513420 : lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
631 : fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
632 513420 : lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
633 : fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
634 513420 : lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
635 2053680 : r1(:) = pos(:, index_a) + imass1*dtsqby2*MATMUL(r_shake, fc1)
636 2053680 : r2(:) = pos(:, index_b) + imass2*dtsqby2*MATMUL(r_shake, fc2)
637 2053680 : r3(:) = pos(:, index_c) + imass3*dtsqby2*MATMUL(r_shake, fc3)
638 2053680 : v1(:) = vel(:, index_a) + imass1*dtby2*MATMUL(v_shake, fc1)
639 2053680 : v2(:) = vel(:, index_b) + imass2*dtby2*MATMUL(v_shake, fc2)
640 2053680 : v3(:) = vel(:, index_c) + imass3*dtby2*MATMUL(v_shake, fc3)
641 513420 : r12 = r1 - r2
642 513420 : r13 = r1 - r3
643 513420 : r23 = r2 - r3
644 :
645 : ! compute the tolerance:
646 : sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* &
647 513420 : g3x3_list(iconst)%dab
648 128355 : max_sigma = MAX(max_sigma, ABS(sigma))
649 : sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* &
650 513420 : g3x3_list(iconst)%dac
651 128355 : max_sigma = MAX(max_sigma, ABS(sigma))
652 : sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* &
653 513420 : g3x3_list(iconst)%dbc
654 128355 : max_sigma = MAX(max_sigma, ABS(sigma))
655 :
656 : ! update positions with full multiplier
657 513420 : pos(:, index_a) = r1(:)
658 513420 : pos(:, index_b) = r2(:)
659 513420 : pos(:, index_c) = r3(:)
660 :
661 : ! update velocites with full multiplier
662 513420 : vel(:, index_a) = v1(:)
663 513420 : vel(:, index_b) = v2(:)
664 641775 : vel(:, index_c) = v3(:)
665 : END DO
666 128355 : END SUBROUTINE shake_roll_3x3_low
667 :
668 : ! **************************************************************************************************
669 : !> \brief ...
670 : !> \param fixd_list ...
671 : !> \param g3x3_list ...
672 : !> \param lg3x3 ...
673 : !> \param first_atom ...
674 : !> \param particle_set ...
675 : !> \param vel ...
676 : !> \param r_rattle ...
677 : !> \param dt ...
678 : !> \param veps ...
679 : !> \par History
680 : !> none
681 : ! **************************************************************************************************
682 97999 : SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
683 97999 : particle_set, vel, r_rattle, dt, veps)
684 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
685 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
686 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
687 : INTEGER, INTENT(IN) :: first_atom
688 : TYPE(particle_type), POINTER :: particle_set(:)
689 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
690 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
691 : REAL(kind=dp), INTENT(in) :: dt
692 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
693 :
694 : INTEGER :: iconst, index_a, index_b, index_c
695 : REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, mass
696 : REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, fc1, fc2, &
697 : fc3, lambda, r12, r13, r23, v12, v13, &
698 : v23
699 : REAL(KIND=dp), DIMENSION(3, 1) :: bvec
700 : REAL(KIND=dp), DIMENSION(3, 3) :: amat
701 : TYPE(atomic_kind_type), POINTER :: atomic_kind
702 :
703 97999 : idt = 1.0_dp/dt
704 97999 : dtby2 = dt*.5_dp
705 195998 : DO iconst = 1, SIZE(g3x3_list)
706 97999 : IF (g3x3_list(iconst)%restraint%active) CYCLE
707 97999 : index_a = g3x3_list(iconst)%a + first_atom - 1
708 97999 : index_b = g3x3_list(iconst)%b + first_atom - 1
709 97999 : index_c = g3x3_list(iconst)%c + first_atom - 1
710 391996 : v12(:) = vel(:, index_a) - vel(:, index_b)
711 391996 : v13(:) = vel(:, index_a) - vel(:, index_c)
712 391996 : v23(:) = vel(:, index_b) - vel(:, index_c)
713 391996 : r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
714 391996 : r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
715 391996 : r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
716 97999 : atomic_kind => particle_set(index_a)%atomic_kind
717 97999 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
718 97999 : imass1 = 1.0_dp/mass
719 97999 : atomic_kind => particle_set(index_b)%atomic_kind
720 97999 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
721 97999 : imass2 = 1.0_dp/mass
722 97999 : atomic_kind => particle_set(index_c)%atomic_kind
723 97999 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
724 97999 : imass3 = 1.0_dp/mass
725 391996 : lg3x3(iconst)%fa = -2.0_dp*r12
726 391996 : lg3x3(iconst)%fb = -2.0_dp*r13
727 391996 : lg3x3(iconst)%fc = -2.0_dp*r23
728 : ! Check for fixed atom constraints
729 : CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
730 97999 : index_a, index_b, index_c, fixd_list, lg3x3(iconst))
731 : ! roll the fc
732 1273987 : f_roll1 = MATMUL(r_rattle, lg3x3(iconst)%fa)
733 1273987 : f_roll2 = MATMUL(r_rattle, lg3x3(iconst)%fb)
734 1273987 : f_roll3 = MATMUL(r_rattle, lg3x3(iconst)%fc)
735 :
736 : ! construct matrix
737 391996 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, f_roll1)
738 391996 : amat(1, 2) = imass1*DOT_PRODUCT(r12, f_roll2)
739 391996 : amat(1, 3) = -imass2*DOT_PRODUCT(r12, f_roll3)
740 391996 : amat(2, 1) = imass1*DOT_PRODUCT(r13, f_roll1)
741 391996 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, f_roll2)
742 391996 : amat(2, 3) = imass3*DOT_PRODUCT(r13, f_roll3)
743 391996 : amat(3, 1) = -imass2*DOT_PRODUCT(r23, f_roll1)
744 391996 : amat(3, 2) = imass3*DOT_PRODUCT(r23, f_roll2)
745 391996 : amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, f_roll3)
746 :
747 : ! construct solution vector
748 2155978 : bvec(1, 1) = DOT_PRODUCT(r12, v12 + MATMUL(veps, r12))
749 2155978 : bvec(2, 1) = DOT_PRODUCT(r13, v13 + MATMUL(veps, r13))
750 2155978 : bvec(3, 1) = DOT_PRODUCT(r23, v23 + MATMUL(veps, r23))
751 489995 : bvec = -bvec*2.0_dp*idt
752 :
753 : ! get lambda
754 97999 : CALL solve_system(amat, 3, bvec)
755 97999 : lambda(:) = bvec(:, 1)
756 391996 : lg3x3(iconst)%lambda(:) = lambda
757 :
758 : fc1 = lambda(1)*f_roll1 + &
759 391996 : lambda(2)*f_roll2
760 : fc2 = -lambda(1)*f_roll1 + &
761 391996 : lambda(3)*f_roll3
762 : fc3 = -lambda(2)*f_roll2 - &
763 391996 : lambda(3)*f_roll3
764 391996 : vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
765 391996 : vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
766 587994 : vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
767 : END DO
768 97999 : END SUBROUTINE rattle_roll_3x3_low
769 :
770 : ! **************************************************************************************************
771 : !> \brief ...
772 : !> \param fixd_list ...
773 : !> \param g3x3_list ...
774 : !> \param lg3x3 ...
775 : !> \param first_atom ...
776 : !> \param particle_set ...
777 : !> \param vel ...
778 : !> \param dt ...
779 : !> \par History
780 : !> none
781 : ! **************************************************************************************************
782 143836 : SUBROUTINE rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
783 143836 : particle_set, vel, dt)
784 :
785 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
786 : TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
787 : TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
788 : INTEGER, INTENT(IN) :: first_atom
789 : TYPE(particle_type), POINTER :: particle_set(:)
790 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
791 : REAL(kind=dp), INTENT(in) :: dt
792 :
793 : INTEGER :: iconst, index_a, index_b, index_c
794 : REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, mass
795 : REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, r12, r13, r23, v12, v13, &
796 : v23
797 : REAL(KIND=dp), DIMENSION(3, 1) :: bvec
798 : REAL(KIND=dp), DIMENSION(3, 3) :: amat
799 : TYPE(atomic_kind_type), POINTER :: atomic_kind
800 :
801 143836 : idt = 1.0_dp/dt
802 143836 : dtby2 = dt*.5_dp
803 287672 : DO iconst = 1, SIZE(g3x3_list)
804 143836 : IF (g3x3_list(iconst)%restraint%active) CYCLE
805 143816 : index_a = g3x3_list(iconst)%a + first_atom - 1
806 143816 : index_b = g3x3_list(iconst)%b + first_atom - 1
807 143816 : index_c = g3x3_list(iconst)%c + first_atom - 1
808 575264 : v12(:) = vel(:, index_a) - vel(:, index_b)
809 575264 : v13(:) = vel(:, index_a) - vel(:, index_c)
810 575264 : v23(:) = vel(:, index_b) - vel(:, index_c)
811 575264 : r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
812 575264 : r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
813 575264 : r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
814 143816 : atomic_kind => particle_set(index_a)%atomic_kind
815 143816 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
816 143816 : imass1 = 1.0_dp/mass
817 143816 : atomic_kind => particle_set(index_b)%atomic_kind
818 143816 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
819 143816 : imass2 = 1.0_dp/mass
820 143816 : atomic_kind => particle_set(index_c)%atomic_kind
821 143816 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
822 143816 : imass3 = 1.0_dp/mass
823 575264 : lg3x3(iconst)%fa = -2.0_dp*r12
824 575264 : lg3x3(iconst)%fb = -2.0_dp*r13
825 575264 : lg3x3(iconst)%fc = -2.0_dp*r23
826 : ! Check for fixed atom constraints
827 : CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
828 143816 : index_a, index_b, index_c, fixd_list, lg3x3(iconst))
829 : ! construct matrix
830 575264 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, lg3x3(iconst)%fa)
831 575264 : amat(1, 2) = imass1*DOT_PRODUCT(r12, lg3x3(iconst)%fb)
832 575264 : amat(1, 3) = -imass2*DOT_PRODUCT(r12, lg3x3(iconst)%fc)
833 575264 : amat(2, 1) = imass1*DOT_PRODUCT(r13, lg3x3(iconst)%fa)
834 575264 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, lg3x3(iconst)%fb)
835 575264 : amat(2, 3) = imass3*DOT_PRODUCT(r13, lg3x3(iconst)%fc)
836 575264 : amat(3, 1) = -imass2*DOT_PRODUCT(r23, lg3x3(iconst)%fa)
837 575264 : amat(3, 2) = imass3*DOT_PRODUCT(r23, lg3x3(iconst)%fb)
838 575264 : amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, lg3x3(iconst)%fc)
839 :
840 : ! construct solution vector
841 575264 : bvec(1, 1) = DOT_PRODUCT(r12, v12)
842 575264 : bvec(2, 1) = DOT_PRODUCT(r13, v13)
843 575264 : bvec(3, 1) = DOT_PRODUCT(r23, v23)
844 719080 : bvec = -bvec*2.0_dp*idt
845 :
846 : ! get lambda
847 143816 : CALL solve_system(amat, 3, bvec)
848 575264 : lg3x3(iconst)%lambda(:) = bvec(:, 1)
849 :
850 : fc1 = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
851 575264 : lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb
852 : fc2 = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
853 575264 : lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
854 : fc3 = -lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb - &
855 575264 : lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
856 575264 : vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
857 575264 : vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
858 862916 : vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
859 : END DO
860 143836 : END SUBROUTINE rattle_3x3_low
861 :
862 293997 : END MODULE constraint_3x3
|
