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_4x6
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE constraint_fxd, ONLY: check_fixed_atom_cns_g4x6
17 : USE kinds, ONLY: dp
18 : USE linear_systems, ONLY: solve_system
19 : USE molecule_kind_types, ONLY: fixd_constraint_type,&
20 : g4x6_constraint_type,&
21 : get_molecule_kind,&
22 : molecule_kind_type
23 : USE molecule_types, ONLY: get_molecule,&
24 : global_constraint_type,&
25 : local_g4x6_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_4x6_int, &
34 : rattle_4x6_int, &
35 : shake_roll_4x6_int, &
36 : rattle_roll_4x6_int, &
37 : shake_4x6_ext, &
38 : rattle_4x6_ext, &
39 : shake_roll_4x6_ext, &
40 : rattle_roll_4x6_ext
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_4x6'
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Intramolecular shake_4x6
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 2466 : SUBROUTINE shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake, &
59 : max_sigma)
60 : TYPE(molecule_type), POINTER :: molecule
61 : TYPE(particle_type), POINTER :: particle_set(:)
62 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
63 : REAL(kind=dp), INTENT(in) :: dt
64 : INTEGER, INTENT(IN) :: ishake
65 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
66 :
67 : INTEGER :: first_atom, ng4x6
68 2466 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
69 2466 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
70 2466 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
71 : TYPE(molecule_kind_type), POINTER :: molecule_kind
72 :
73 2466 : NULLIFY (fixd_list)
74 2466 : molecule_kind => molecule%molecule_kind
75 : CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
76 2466 : fixd_list=fixd_list, g4x6_list=g4x6_list)
77 2466 : CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
78 : ! Real Shake
79 : CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
80 2466 : particle_set, pos, vel, dt, ishake, max_sigma)
81 :
82 2466 : END SUBROUTINE shake_4x6_int
83 :
84 : ! **************************************************************************************************
85 : !> \brief Intramolecular shake_4x6_roll
86 : !> \param molecule ...
87 : !> \param particle_set ...
88 : !> \param pos ...
89 : !> \param vel ...
90 : !> \param r_shake ...
91 : !> \param dt ...
92 : !> \param ishake ...
93 : !> \param max_sigma ...
94 : !> \par History
95 : !> none
96 : ! **************************************************************************************************
97 2225 : SUBROUTINE shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, &
98 : dt, ishake, max_sigma)
99 : TYPE(molecule_type), POINTER :: molecule
100 : TYPE(particle_type), POINTER :: particle_set(:)
101 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
102 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake
103 : REAL(kind=dp), INTENT(in) :: dt
104 : INTEGER, INTENT(IN) :: ishake
105 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
106 :
107 : INTEGER :: first_atom, ng4x6
108 2225 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
109 2225 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
110 2225 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
111 : TYPE(molecule_kind_type), POINTER :: molecule_kind
112 :
113 2225 : NULLIFY (fixd_list)
114 2225 : molecule_kind => molecule%molecule_kind
115 : CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
116 2225 : fixd_list=fixd_list, g4x6_list=g4x6_list)
117 2225 : CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
118 : ! Real Shake
119 : CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
120 2225 : particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
121 :
122 2225 : END SUBROUTINE shake_roll_4x6_int
123 :
124 : ! **************************************************************************************************
125 : !> \brief Intramolecular rattle_4x6
126 : !> \param molecule ...
127 : !> \param particle_set ...
128 : !> \param vel ...
129 : !> \param dt ...
130 : !> \par History
131 : !> none
132 : ! **************************************************************************************************
133 682 : SUBROUTINE rattle_4x6_int(molecule, particle_set, vel, dt)
134 : TYPE(molecule_type), POINTER :: molecule
135 : TYPE(particle_type), POINTER :: particle_set(:)
136 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
137 : REAL(kind=dp), INTENT(in) :: dt
138 :
139 : INTEGER :: first_atom, ng4x6
140 682 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
141 682 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
142 682 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
143 : TYPE(molecule_kind_type), POINTER :: molecule_kind
144 :
145 682 : NULLIFY (fixd_list)
146 682 : molecule_kind => molecule%molecule_kind
147 : CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
148 682 : fixd_list=fixd_list, g4x6_list=g4x6_list)
149 682 : CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
150 : ! Real Rattle
151 : CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
152 682 : particle_set, vel, dt)
153 :
154 682 : END SUBROUTINE rattle_4x6_int
155 :
156 : ! **************************************************************************************************
157 : !> \brief Intramolecular rattle_4x6_roll
158 : !> \param molecule ...
159 : !> \param particle_set ...
160 : !> \param vel ...
161 : !> \param r_rattle ...
162 : !> \param dt ...
163 : !> \param veps ...
164 : !> \par History
165 : !> none
166 : ! **************************************************************************************************
167 1024 : SUBROUTINE rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, veps)
168 : TYPE(molecule_type), POINTER :: molecule
169 : TYPE(particle_type), POINTER :: particle_set(:)
170 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
171 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
172 : REAL(kind=dp), INTENT(in) :: dt
173 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
174 :
175 : INTEGER :: first_atom, ng4x6
176 1024 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
177 1024 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
178 1024 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
179 : TYPE(molecule_kind_type), POINTER :: molecule_kind
180 :
181 1024 : NULLIFY (fixd_list)
182 1024 : molecule_kind => molecule%molecule_kind
183 : CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
184 1024 : fixd_list=fixd_list, g4x6_list=g4x6_list)
185 1024 : CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
186 : ! Real Rattle
187 : CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
188 1024 : particle_set, vel, r_rattle, dt, veps)
189 :
190 1024 : END SUBROUTINE rattle_roll_4x6_int
191 :
192 : ! **************************************************************************************************
193 : !> \brief Intramolecular shake_4x6
194 : !> \param gci ...
195 : !> \param particle_set ...
196 : !> \param pos ...
197 : !> \param vel ...
198 : !> \param dt ...
199 : !> \param ishake ...
200 : !> \param max_sigma ...
201 : !> \par History
202 : !> none
203 : ! **************************************************************************************************
204 48 : SUBROUTINE shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake, &
205 : max_sigma)
206 :
207 : TYPE(global_constraint_type), POINTER :: gci
208 : TYPE(particle_type), POINTER :: particle_set(:)
209 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
210 : REAL(kind=dp), INTENT(in) :: dt
211 : INTEGER, INTENT(IN) :: ishake
212 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
213 :
214 : INTEGER :: first_atom, ng4x6
215 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
216 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
217 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
218 :
219 48 : first_atom = 1
220 48 : ng4x6 = gci%ng4x6
221 48 : fixd_list => gci%fixd_list
222 48 : g4x6_list => gci%g4x6_list
223 48 : lg4x6 => gci%lg4x6
224 : ! Real Shake
225 : CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
226 48 : particle_set, pos, vel, dt, ishake, max_sigma)
227 :
228 48 : END SUBROUTINE shake_4x6_ext
229 :
230 : ! **************************************************************************************************
231 : !> \brief Intramolecular shake_4x6_roll
232 : !> \param gci ...
233 : !> \param particle_set ...
234 : !> \param pos ...
235 : !> \param vel ...
236 : !> \param r_shake ...
237 : !> \param dt ...
238 : !> \param ishake ...
239 : !> \param max_sigma ...
240 : !> \par History
241 : !> none
242 : ! **************************************************************************************************
243 0 : SUBROUTINE shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, &
244 : dt, ishake, max_sigma)
245 :
246 : TYPE(global_constraint_type), POINTER :: gci
247 : TYPE(particle_type), POINTER :: particle_set(:)
248 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
249 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake
250 : REAL(kind=dp), INTENT(in) :: dt
251 : INTEGER, INTENT(IN) :: ishake
252 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
253 :
254 : INTEGER :: first_atom, ng4x6
255 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
256 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
257 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
258 :
259 0 : first_atom = 1
260 0 : ng4x6 = gci%ng4x6
261 0 : fixd_list => gci%fixd_list
262 0 : g4x6_list => gci%g4x6_list
263 0 : lg4x6 => gci%lg4x6
264 : ! Real Shake
265 : CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
266 0 : particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
267 :
268 0 : END SUBROUTINE shake_roll_4x6_ext
269 :
270 : ! **************************************************************************************************
271 : !> \brief Intramolecular rattle_4x6
272 : !> \param gci ...
273 : !> \param particle_set ...
274 : !> \param vel ...
275 : !> \param dt ...
276 : !> \par History
277 : !> none
278 : ! **************************************************************************************************
279 20 : SUBROUTINE rattle_4x6_ext(gci, particle_set, vel, dt)
280 :
281 : TYPE(global_constraint_type), POINTER :: gci
282 : TYPE(particle_type), POINTER :: particle_set(:)
283 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
284 : REAL(kind=dp), INTENT(in) :: dt
285 :
286 : INTEGER :: first_atom, ng4x6
287 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
288 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
289 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
290 :
291 20 : first_atom = 1
292 20 : ng4x6 = gci%ng4x6
293 20 : fixd_list => gci%fixd_list
294 20 : g4x6_list => gci%g4x6_list
295 20 : lg4x6 => gci%lg4x6
296 : ! Real Rattle
297 : CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
298 20 : particle_set, vel, dt)
299 :
300 20 : END SUBROUTINE rattle_4x6_ext
301 :
302 : ! **************************************************************************************************
303 : !> \brief Intramolecular rattle_4x6_roll
304 : !> \param gci ...
305 : !> \param particle_set ...
306 : !> \param vel ...
307 : !> \param r_rattle ...
308 : !> \param dt ...
309 : !> \param veps ...
310 : !> \par History
311 : !> none
312 : ! **************************************************************************************************
313 0 : SUBROUTINE rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, veps)
314 :
315 : TYPE(global_constraint_type), POINTER :: gci
316 : TYPE(particle_type), POINTER :: particle_set(:)
317 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
318 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
319 : REAL(kind=dp), INTENT(in) :: dt
320 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
321 :
322 : INTEGER :: first_atom, ng4x6
323 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
324 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
325 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
326 :
327 0 : first_atom = 1
328 0 : ng4x6 = gci%ng4x6
329 0 : fixd_list => gci%fixd_list
330 0 : g4x6_list => gci%g4x6_list
331 0 : lg4x6 => gci%lg4x6
332 : ! Real Rattle
333 : CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
334 0 : particle_set, vel, r_rattle, dt, veps)
335 :
336 0 : END SUBROUTINE rattle_roll_4x6_ext
337 :
338 : ! **************************************************************************************************
339 : !> \brief ...
340 : !> \param fixd_list ...
341 : !> \param g4x6_list ...
342 : !> \param lg4x6 ...
343 : !> \param ng4x6 ...
344 : !> \param first_atom ...
345 : !> \param particle_set ...
346 : !> \param pos ...
347 : !> \param vel ...
348 : !> \param dt ...
349 : !> \param ishake ...
350 : !> \param max_sigma ...
351 : !> \par History
352 : !> none
353 : ! **************************************************************************************************
354 2514 : SUBROUTINE shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
355 2514 : particle_set, pos, vel, dt, ishake, max_sigma)
356 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
357 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
358 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
359 : INTEGER, INTENT(IN) :: ng4x6, first_atom
360 : TYPE(particle_type), POINTER :: particle_set(:)
361 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
362 : REAL(kind=dp), INTENT(in) :: dt
363 : INTEGER, INTENT(IN) :: ishake
364 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
365 :
366 : INTEGER :: iconst, index_a, index_b, index_c, &
367 : index_d
368 : REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, &
369 : imass3, imass4, sigma
370 : REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r0_12, r0_13, r0_14, &
371 : r0_23, r0_24, r0_34, r1, r12, r13, &
372 : r14, r2, r23, r24, r3, r34, r4, v1, &
373 : v2, v3, v4, vec
374 : REAL(KIND=dp), DIMENSION(6, 1) :: bvec
375 : REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp
376 : TYPE(atomic_kind_type), POINTER :: atomic_kind
377 :
378 2514 : dtsqby2 = dt*dt*.5_dp
379 2514 : idtsq = 1.0_dp/dt/dt
380 2514 : dtby2 = dt*.5_dp
381 5028 : DO iconst = 1, ng4x6
382 2514 : IF (g4x6_list(iconst)%restraint%active) CYCLE
383 2504 : index_a = g4x6_list(iconst)%a + first_atom - 1
384 2504 : index_b = g4x6_list(iconst)%b + first_atom - 1
385 2504 : index_c = g4x6_list(iconst)%c + first_atom - 1
386 2504 : index_d = g4x6_list(iconst)%d + first_atom - 1
387 2504 : IF (ishake == 1) THEN
388 2768 : r0_12(:) = pos(:, index_a) - pos(:, index_b)
389 2768 : r0_13(:) = pos(:, index_a) - pos(:, index_c)
390 2768 : r0_14(:) = pos(:, index_a) - pos(:, index_d)
391 2768 : r0_23(:) = pos(:, index_b) - pos(:, index_c)
392 2768 : r0_24(:) = pos(:, index_b) - pos(:, index_d)
393 2768 : r0_34(:) = pos(:, index_c) - pos(:, index_d)
394 692 : atomic_kind => particle_set(index_a)%atomic_kind
395 692 : imass1 = 1.0_dp/atomic_kind%mass
396 692 : atomic_kind => particle_set(index_b)%atomic_kind
397 692 : imass2 = 1.0_dp/atomic_kind%mass
398 692 : atomic_kind => particle_set(index_c)%atomic_kind
399 692 : imass3 = 1.0_dp/atomic_kind%mass
400 692 : atomic_kind => particle_set(index_d)%atomic_kind
401 692 : imass4 = 1.0_dp/atomic_kind%mass
402 : lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
403 2768 : lg4x6(iconst)%rb_old)
404 : lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
405 2768 : lg4x6(iconst)%rc_old)
406 : lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
407 2768 : lg4x6(iconst)%rd_old)
408 : lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
409 2768 : lg4x6(iconst)%rc_old)
410 : lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
411 2768 : lg4x6(iconst)%rd_old)
412 : lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
413 2768 : lg4x6(iconst)%rd_old)
414 : ! Check for fixed atom constraints
415 : CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
416 692 : index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
417 : ! construct matrix
418 2768 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, lg4x6(iconst)%fa)
419 2768 : amat(1, 2) = imass1*DOT_PRODUCT(r0_12, lg4x6(iconst)%fb)
420 2768 : amat(1, 3) = imass1*DOT_PRODUCT(r0_12, lg4x6(iconst)%fc)
421 2768 : amat(1, 4) = -imass2*DOT_PRODUCT(r0_12, lg4x6(iconst)%fd)
422 2768 : amat(1, 5) = -imass2*DOT_PRODUCT(r0_12, lg4x6(iconst)%fe)
423 692 : amat(1, 6) = 0.0_dp
424 :
425 2768 : amat(2, 1) = imass1*DOT_PRODUCT(r0_13, lg4x6(iconst)%fa)
426 2768 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, lg4x6(iconst)%fb)
427 2768 : amat(2, 3) = imass1*DOT_PRODUCT(r0_13, lg4x6(iconst)%fc)
428 2768 : amat(2, 4) = imass3*DOT_PRODUCT(r0_13, lg4x6(iconst)%fd)
429 692 : amat(2, 5) = 0.0_dp
430 2768 : amat(2, 6) = -imass3*DOT_PRODUCT(r0_13, lg4x6(iconst)%ff)
431 :
432 2768 : amat(3, 1) = imass1*DOT_PRODUCT(r0_14, lg4x6(iconst)%fa)
433 2768 : amat(3, 2) = imass1*DOT_PRODUCT(r0_14, lg4x6(iconst)%fb)
434 2768 : amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r0_14, lg4x6(iconst)%fc)
435 692 : amat(3, 4) = 0.0_dp
436 2768 : amat(3, 5) = imass4*DOT_PRODUCT(r0_14, lg4x6(iconst)%fe)
437 2768 : amat(3, 6) = imass4*DOT_PRODUCT(r0_14, lg4x6(iconst)%ff)
438 :
439 2768 : amat(4, 1) = -imass2*DOT_PRODUCT(r0_23, lg4x6(iconst)%fa)
440 2768 : amat(4, 2) = imass3*DOT_PRODUCT(r0_23, lg4x6(iconst)%fb)
441 692 : amat(4, 3) = 0.0_dp
442 2768 : amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r0_23, lg4x6(iconst)%fd)
443 2768 : amat(4, 5) = imass2*DOT_PRODUCT(r0_23, lg4x6(iconst)%fe)
444 2768 : amat(4, 6) = -imass3*DOT_PRODUCT(r0_23, lg4x6(iconst)%ff)
445 :
446 2768 : amat(5, 1) = -imass2*DOT_PRODUCT(r0_24, lg4x6(iconst)%fa)
447 692 : amat(5, 2) = 0.0_dp
448 2768 : amat(5, 3) = imass4*DOT_PRODUCT(r0_24, lg4x6(iconst)%fc)
449 2768 : amat(5, 4) = imass2*DOT_PRODUCT(r0_24, lg4x6(iconst)%fd)
450 2768 : amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r0_24, lg4x6(iconst)%fe)
451 2768 : amat(5, 6) = imass4*DOT_PRODUCT(r0_24, lg4x6(iconst)%ff)
452 :
453 692 : amat(6, 1) = 0.0_dp
454 2768 : amat(6, 2) = -imass3*DOT_PRODUCT(r0_34, lg4x6(iconst)%fb)
455 2768 : amat(6, 3) = imass4*DOT_PRODUCT(r0_34, lg4x6(iconst)%fc)
456 2768 : amat(6, 4) = -imass3*DOT_PRODUCT(r0_34, lg4x6(iconst)%fd)
457 2768 : amat(6, 5) = imass4*DOT_PRODUCT(r0_34, lg4x6(iconst)%fe)
458 2768 : amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r0_34, lg4x6(iconst)%ff)
459 : ! Store values
460 2768 : lg4x6(iconst)%r0_12 = r0_12
461 2768 : lg4x6(iconst)%r0_13 = r0_13
462 2768 : lg4x6(iconst)%r0_14 = r0_14
463 2768 : lg4x6(iconst)%r0_23 = r0_23
464 2768 : lg4x6(iconst)%r0_24 = r0_24
465 2768 : lg4x6(iconst)%r0_34 = r0_34
466 29756 : lg4x6(iconst)%amat = amat
467 692 : lg4x6(iconst)%imass1 = imass1
468 692 : lg4x6(iconst)%imass2 = imass2
469 692 : lg4x6(iconst)%imass3 = imass3
470 692 : lg4x6(iconst)%imass4 = imass4
471 4844 : lg4x6(iconst)%lambda_old = 0.0_dp
472 4844 : lg4x6(iconst)%del_lambda = 0.0_dp
473 : ELSE
474 : ! Retrieve values
475 7248 : r0_12 = lg4x6(iconst)%r0_12
476 7248 : r0_13 = lg4x6(iconst)%r0_13
477 7248 : r0_14 = lg4x6(iconst)%r0_14
478 7248 : r0_23 = lg4x6(iconst)%r0_23
479 7248 : r0_24 = lg4x6(iconst)%r0_24
480 7248 : r0_34 = lg4x6(iconst)%r0_34
481 77916 : amat = lg4x6(iconst)%amat
482 1812 : imass1 = lg4x6(iconst)%imass1
483 1812 : imass2 = lg4x6(iconst)%imass2
484 1812 : imass3 = lg4x6(iconst)%imass3
485 1812 : imass4 = lg4x6(iconst)%imass4
486 : END IF
487 :
488 : ! Iterate until convergence:
489 : vec = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa*(imass1 + imass2) + &
490 : lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + &
491 : lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc - &
492 : lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd - &
493 10016 : lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe
494 : bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
495 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
496 : vec = lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb*(imass1 + imass3) + &
497 : lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + &
498 : lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc + &
499 : lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd - &
500 10016 : lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
501 : bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
502 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
503 : vec = lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc*(imass1 + imass4) + &
504 : lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + &
505 : lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + &
506 : lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe + &
507 10016 : lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
508 : bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
509 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_14, r0_14)
510 : vec = lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd*(imass2 + imass3) - &
511 : lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + &
512 : lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + &
513 : lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe - &
514 10016 : lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
515 : bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
516 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
517 : vec = lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe*(imass2 + imass4) - &
518 : lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + &
519 : lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc + &
520 : lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd + &
521 10016 : lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
522 : bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
523 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_24, r0_24)
524 : vec = lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff*(imass3 + imass4) - &
525 : lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + &
526 : lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc - &
527 : lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd + &
528 10016 : lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe
529 : bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
530 20032 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_34, r0_34)
531 20032 : bvec = bvec*idtsq
532 :
533 : ! get lambda
534 2504 : atemp = amat
535 2504 : CALL solve_system(atemp, 6, bvec)
536 17528 : lg4x6(iconst)%lambda(:) = bvec(:, 1)
537 : lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
538 17528 : lg4x6(iconst)%lambda_old(:)
539 17528 : lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:)
540 :
541 : fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
542 : lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + &
543 10016 : lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc
544 : fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
545 : lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
546 10016 : lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe
547 : fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - &
548 : lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
549 10016 : lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
550 : fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - &
551 : lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - &
552 10016 : lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
553 10016 : r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
554 10016 : r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
555 10016 : r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
556 10016 : r4(:) = pos(:, index_d) + imass4*dtsqby2*fc4(:)
557 10016 : v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
558 10016 : v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
559 10016 : v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
560 10016 : v4(:) = vel(:, index_d) + imass4*dtby2*fc4(:)
561 10016 : r12 = r1 - r2
562 10016 : r13 = r1 - r3
563 10016 : r14 = r1 - r4
564 10016 : r23 = r2 - r3
565 10016 : r24 = r2 - r4
566 10016 : r34 = r3 - r4
567 : ! compute the tolerance:
568 : sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* &
569 10016 : g4x6_list(iconst)%dab
570 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
571 : sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* &
572 10016 : g4x6_list(iconst)%dac
573 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
574 : sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* &
575 10016 : g4x6_list(iconst)%dad
576 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
577 : sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* &
578 10016 : g4x6_list(iconst)%dbc
579 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
580 : sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* &
581 10016 : g4x6_list(iconst)%dbd
582 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
583 : sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* &
584 10016 : g4x6_list(iconst)%dcd
585 2504 : max_sigma = MAX(max_sigma, ABS(sigma))
586 :
587 : ! update positions with full multiplier
588 10016 : pos(:, index_a) = r1(:)
589 10016 : pos(:, index_b) = r2(:)
590 10016 : pos(:, index_c) = r3(:)
591 10016 : pos(:, index_d) = r4(:)
592 :
593 : ! update velocites with full multiplier
594 10016 : vel(:, index_a) = v1(:)
595 10016 : vel(:, index_b) = v2(:)
596 10016 : vel(:, index_c) = v3(:)
597 12540 : vel(:, index_d) = v4(:)
598 : END DO
599 2514 : END SUBROUTINE shake_4x6_low
600 :
601 : ! **************************************************************************************************
602 : !> \brief ...
603 : !> \param fixd_list ...
604 : !> \param g4x6_list ...
605 : !> \param lg4x6 ...
606 : !> \param ng4x6 ...
607 : !> \param first_atom ...
608 : !> \param particle_set ...
609 : !> \param pos ...
610 : !> \param vel ...
611 : !> \param r_shake ...
612 : !> \param dt ...
613 : !> \param ishake ...
614 : !> \param max_sigma ...
615 : !> \par History
616 : !> none
617 : ! **************************************************************************************************
618 2225 : SUBROUTINE shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
619 2225 : particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
620 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
621 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
622 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
623 : INTEGER, INTENT(IN) :: ng4x6, first_atom
624 : TYPE(particle_type), POINTER :: particle_set(:)
625 : REAL(KIND=dp), INTENT(INOUT) :: pos(:, :), vel(:, :)
626 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_shake
627 : REAL(kind=dp), INTENT(in) :: dt
628 : INTEGER, INTENT(IN) :: ishake
629 : REAL(KIND=dp), INTENT(INOUT) :: max_sigma
630 :
631 : INTEGER :: iconst, index_a, index_b, index_c, &
632 : index_d
633 : REAL(KIND=dp) :: dtby2, dtsqby2, idtsq, imass1, imass2, &
634 : imass3, imass4, sigma
635 : REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, f_roll5, f_roll6, fc1, &
636 : fc2, fc3, fc4, r0_12, r0_13, r0_14, r0_23, r0_24, r0_34, r1, r12, r13, r14, r2, r23, r24, &
637 : r3, r34, r4, v1, v2, v3, v4, vec
638 : REAL(KIND=dp), DIMENSION(6, 1) :: bvec
639 : REAL(KIND=dp), DIMENSION(6, 6) :: amat, atemp
640 : TYPE(atomic_kind_type), POINTER :: atomic_kind
641 :
642 2225 : dtsqby2 = dt*dt*.5_dp
643 2225 : idtsq = 1.0_dp/dt/dt
644 2225 : dtby2 = dt*.5_dp
645 4450 : DO iconst = 1, ng4x6
646 2225 : IF (g4x6_list(iconst)%restraint%active) CYCLE
647 2225 : index_a = g4x6_list(iconst)%a + first_atom - 1
648 2225 : index_b = g4x6_list(iconst)%b + first_atom - 1
649 2225 : index_c = g4x6_list(iconst)%c + first_atom - 1
650 2225 : index_d = g4x6_list(iconst)%d + first_atom - 1
651 2225 : IF (ishake == 1) THEN
652 4352 : r0_12(:) = pos(:, index_a) - pos(:, index_b)
653 4352 : r0_13(:) = pos(:, index_a) - pos(:, index_c)
654 4352 : r0_23(:) = pos(:, index_b) - pos(:, index_c)
655 4352 : r0_14(:) = pos(:, index_a) - pos(:, index_d)
656 4352 : r0_24(:) = pos(:, index_b) - pos(:, index_d)
657 4352 : r0_34(:) = pos(:, index_c) - pos(:, index_d)
658 1088 : atomic_kind => particle_set(index_a)%atomic_kind
659 1088 : imass1 = 1.0_dp/atomic_kind%mass
660 1088 : atomic_kind => particle_set(index_b)%atomic_kind
661 1088 : imass2 = 1.0_dp/atomic_kind%mass
662 1088 : atomic_kind => particle_set(index_c)%atomic_kind
663 1088 : imass3 = 1.0_dp/atomic_kind%mass
664 1088 : atomic_kind => particle_set(index_d)%atomic_kind
665 1088 : imass4 = 1.0_dp/atomic_kind%mass
666 : lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
667 4352 : lg4x6(iconst)%rb_old)
668 : lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
669 4352 : lg4x6(iconst)%rc_old)
670 : lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
671 4352 : lg4x6(iconst)%rd_old)
672 : lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
673 4352 : lg4x6(iconst)%rc_old)
674 : lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
675 4352 : lg4x6(iconst)%rd_old)
676 : lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
677 4352 : lg4x6(iconst)%rd_old)
678 : ! Check for fixed atom constraints
679 : CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
680 1088 : index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
681 : ! rotate fconst:
682 14144 : f_roll1 = MATMUL(r_shake, lg4x6(iconst)%fa)
683 14144 : f_roll2 = MATMUL(r_shake, lg4x6(iconst)%fb)
684 14144 : f_roll3 = MATMUL(r_shake, lg4x6(iconst)%fc)
685 14144 : f_roll4 = MATMUL(r_shake, lg4x6(iconst)%fd)
686 14144 : f_roll5 = MATMUL(r_shake, lg4x6(iconst)%fe)
687 14144 : f_roll6 = MATMUL(r_shake, lg4x6(iconst)%ff)
688 :
689 : ! construct matrix
690 4352 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, f_roll1)
691 4352 : amat(1, 2) = imass1*DOT_PRODUCT(r0_12, f_roll2)
692 4352 : amat(1, 3) = imass1*DOT_PRODUCT(r0_12, f_roll3)
693 4352 : amat(1, 4) = -imass2*DOT_PRODUCT(r0_12, f_roll4)
694 4352 : amat(1, 5) = -imass2*DOT_PRODUCT(r0_12, f_roll5)
695 1088 : amat(1, 6) = 0.0_dp
696 :
697 4352 : amat(2, 1) = imass1*DOT_PRODUCT(r0_13, f_roll1)
698 4352 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, f_roll2)
699 4352 : amat(2, 3) = imass1*DOT_PRODUCT(r0_13, f_roll3)
700 4352 : amat(2, 4) = imass3*DOT_PRODUCT(r0_13, f_roll4)
701 1088 : amat(2, 5) = 0.0_dp
702 4352 : amat(2, 6) = -imass3*DOT_PRODUCT(r0_13, f_roll6)
703 :
704 4352 : amat(3, 1) = imass1*DOT_PRODUCT(r0_14, f_roll1)
705 4352 : amat(3, 2) = imass1*DOT_PRODUCT(r0_14, f_roll2)
706 4352 : amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r0_14, f_roll3)
707 1088 : amat(3, 4) = 0.0_dp
708 4352 : amat(3, 5) = imass4*DOT_PRODUCT(r0_14, f_roll5)
709 4352 : amat(3, 6) = imass4*DOT_PRODUCT(r0_14, f_roll6)
710 :
711 4352 : amat(4, 1) = -imass2*DOT_PRODUCT(r0_23, f_roll1)
712 4352 : amat(4, 2) = imass3*DOT_PRODUCT(r0_23, f_roll2)
713 1088 : amat(4, 3) = 0.0_dp
714 4352 : amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r0_23, f_roll4)
715 4352 : amat(4, 5) = imass2*DOT_PRODUCT(r0_23, f_roll5)
716 4352 : amat(4, 6) = -imass3*DOT_PRODUCT(r0_23, f_roll6)
717 :
718 4352 : amat(5, 1) = -imass2*DOT_PRODUCT(r0_24, f_roll1)
719 1088 : amat(5, 2) = 0.0_dp
720 4352 : amat(5, 3) = imass4*DOT_PRODUCT(r0_24, f_roll3)
721 4352 : amat(5, 4) = imass2*DOT_PRODUCT(r0_24, f_roll4)
722 4352 : amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r0_24, f_roll5)
723 4352 : amat(5, 6) = imass4*DOT_PRODUCT(r0_24, f_roll6)
724 :
725 1088 : amat(6, 1) = 0.0_dp
726 4352 : amat(6, 2) = -imass3*DOT_PRODUCT(r0_34, f_roll2)
727 4352 : amat(6, 3) = imass4*DOT_PRODUCT(r0_34, f_roll3)
728 4352 : amat(6, 4) = -imass3*DOT_PRODUCT(r0_34, f_roll4)
729 4352 : amat(6, 5) = imass4*DOT_PRODUCT(r0_34, f_roll5)
730 4352 : amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r0_34, f_roll6)
731 : ! Store values
732 4352 : lg4x6(iconst)%r0_12 = r0_12
733 4352 : lg4x6(iconst)%r0_13 = r0_13
734 4352 : lg4x6(iconst)%r0_14 = r0_14
735 4352 : lg4x6(iconst)%r0_23 = r0_23
736 4352 : lg4x6(iconst)%r0_24 = r0_24
737 4352 : lg4x6(iconst)%r0_34 = r0_34
738 4352 : lg4x6(iconst)%f_roll1 = f_roll1
739 4352 : lg4x6(iconst)%f_roll2 = f_roll2
740 4352 : lg4x6(iconst)%f_roll3 = f_roll3
741 4352 : lg4x6(iconst)%f_roll4 = f_roll4
742 4352 : lg4x6(iconst)%f_roll5 = f_roll5
743 4352 : lg4x6(iconst)%f_roll6 = f_roll6
744 46784 : lg4x6(iconst)%amat = amat
745 1088 : lg4x6(iconst)%imass1 = imass1
746 1088 : lg4x6(iconst)%imass2 = imass2
747 1088 : lg4x6(iconst)%imass3 = imass3
748 1088 : lg4x6(iconst)%imass4 = imass4
749 7616 : lg4x6(iconst)%lambda_old = 0.0_dp
750 7616 : lg4x6(iconst)%del_lambda = 0.0_dp
751 : ELSE
752 : ! Retrieve values
753 4548 : r0_12 = lg4x6(iconst)%r0_12
754 4548 : r0_13 = lg4x6(iconst)%r0_13
755 4548 : r0_14 = lg4x6(iconst)%r0_14
756 4548 : r0_23 = lg4x6(iconst)%r0_23
757 4548 : r0_24 = lg4x6(iconst)%r0_24
758 4548 : r0_34 = lg4x6(iconst)%r0_34
759 4548 : f_roll1 = lg4x6(iconst)%f_roll1
760 4548 : f_roll2 = lg4x6(iconst)%f_roll2
761 4548 : f_roll3 = lg4x6(iconst)%f_roll3
762 4548 : f_roll4 = lg4x6(iconst)%f_roll4
763 4548 : f_roll5 = lg4x6(iconst)%f_roll5
764 4548 : f_roll6 = lg4x6(iconst)%f_roll6
765 48891 : amat = lg4x6(iconst)%amat
766 1137 : imass1 = lg4x6(iconst)%imass1
767 1137 : imass2 = lg4x6(iconst)%imass2
768 1137 : imass3 = lg4x6(iconst)%imass3
769 1137 : imass4 = lg4x6(iconst)%imass4
770 : END IF
771 :
772 : ! Iterate until convergence:
773 : vec = lg4x6(iconst)%lambda(1)*f_roll1*(imass1 + imass2) + &
774 : lg4x6(iconst)%lambda(2)*imass1*f_roll2 + &
775 : lg4x6(iconst)%lambda(3)*imass1*f_roll3 - &
776 : lg4x6(iconst)%lambda(4)*imass2*f_roll4 - &
777 8900 : lg4x6(iconst)%lambda(5)*imass2*f_roll5
778 : bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
779 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
780 : vec = lg4x6(iconst)%lambda(2)*f_roll2*(imass1 + imass3) + &
781 : lg4x6(iconst)%lambda(1)*imass1*f_roll1 + &
782 : lg4x6(iconst)%lambda(3)*imass1*f_roll3 + &
783 : lg4x6(iconst)%lambda(4)*imass3*f_roll4 - &
784 8900 : lg4x6(iconst)%lambda(6)*imass3*f_roll6
785 : bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
786 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
787 : vec = lg4x6(iconst)%lambda(3)*f_roll3*(imass1 + imass4) + &
788 : lg4x6(iconst)%lambda(1)*imass1*f_roll1 + &
789 : lg4x6(iconst)%lambda(2)*imass1*f_roll2 + &
790 : lg4x6(iconst)%lambda(5)*imass4*f_roll5 + &
791 8900 : lg4x6(iconst)%lambda(6)*imass4*f_roll6
792 : bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
793 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_14, r0_14)
794 : vec = lg4x6(iconst)%lambda(4)*f_roll4*(imass2 + imass3) - &
795 : lg4x6(iconst)%lambda(1)*imass2*f_roll1 + &
796 : lg4x6(iconst)%lambda(2)*imass3*f_roll2 + &
797 : lg4x6(iconst)%lambda(5)*imass2*f_roll5 - &
798 8900 : lg4x6(iconst)%lambda(6)*imass3*f_roll6
799 : bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
800 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
801 : vec = lg4x6(iconst)%lambda(5)*f_roll5*(imass2 + imass4) - &
802 : lg4x6(iconst)%lambda(1)*imass2*f_roll1 + &
803 : lg4x6(iconst)%lambda(3)*imass4*f_roll3 + &
804 : lg4x6(iconst)%lambda(4)*imass2*f_roll4 + &
805 8900 : lg4x6(iconst)%lambda(6)*imass4*f_roll6
806 : bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
807 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_24, r0_24)
808 : vec = lg4x6(iconst)%lambda(6)*f_roll6*(imass3 + imass4) - &
809 : lg4x6(iconst)%lambda(2)*imass3*f_roll2 + &
810 : lg4x6(iconst)%lambda(3)*imass4*f_roll3 - &
811 : lg4x6(iconst)%lambda(4)*imass3*f_roll4 + &
812 8900 : lg4x6(iconst)%lambda(5)*imass4*f_roll5
813 : bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
814 17800 : - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_34, r0_34)
815 17800 : bvec = bvec*idtsq
816 :
817 : ! get lambda
818 2225 : atemp = amat
819 2225 : CALL solve_system(atemp, 6, bvec)
820 15575 : lg4x6(iconst)%lambda(:) = bvec(:, 1)
821 : lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
822 15575 : lg4x6(iconst)%lambda_old(:)
823 15575 : lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:)
824 :
825 : fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
826 : lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + &
827 8900 : lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc
828 : fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
829 : lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
830 8900 : lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe
831 : fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - &
832 : lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
833 8900 : lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
834 : fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - &
835 : lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - &
836 8900 : lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
837 35600 : r1(:) = pos(:, index_a) + imass1*dtsqby2*MATMUL(r_shake, fc1)
838 35600 : r2(:) = pos(:, index_b) + imass2*dtsqby2*MATMUL(r_shake, fc2)
839 35600 : r3(:) = pos(:, index_c) + imass3*dtsqby2*MATMUL(r_shake, fc3)
840 35600 : r4(:) = pos(:, index_d) + imass4*dtsqby2*MATMUL(r_shake, fc4)
841 35600 : v1(:) = vel(:, index_a) + imass1*dtby2*MATMUL(r_shake, fc1)
842 35600 : v2(:) = vel(:, index_b) + imass2*dtby2*MATMUL(r_shake, fc2)
843 35600 : v3(:) = vel(:, index_c) + imass3*dtby2*MATMUL(r_shake, fc3)
844 35600 : v4(:) = vel(:, index_d) + imass4*dtby2*MATMUL(r_shake, fc4)
845 :
846 8900 : r12 = r1 - r2
847 8900 : r13 = r1 - r3
848 8900 : r23 = r2 - r3
849 8900 : r14 = r1 - r4
850 8900 : r24 = r2 - r4
851 8900 : r34 = r3 - r4
852 : ! compute the tolerance:
853 : sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* &
854 8900 : g4x6_list(iconst)%dab
855 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
856 : sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* &
857 8900 : g4x6_list(iconst)%dac
858 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
859 : sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* &
860 8900 : g4x6_list(iconst)%dad
861 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
862 : sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* &
863 8900 : g4x6_list(iconst)%dbc
864 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
865 : sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* &
866 8900 : g4x6_list(iconst)%dbd
867 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
868 : sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* &
869 8900 : g4x6_list(iconst)%dcd
870 2225 : max_sigma = MAX(max_sigma, ABS(sigma))
871 :
872 : ! update positions with full multiplier
873 8900 : pos(:, index_a) = r1(:)
874 8900 : pos(:, index_b) = r2(:)
875 8900 : pos(:, index_c) = r3(:)
876 8900 : pos(:, index_d) = r4(:)
877 :
878 : ! update velocites with full multiplier
879 8900 : vel(:, index_a) = v1(:)
880 8900 : vel(:, index_b) = v2(:)
881 8900 : vel(:, index_c) = v3(:)
882 11125 : vel(:, index_d) = v4(:)
883 : END DO
884 2225 : END SUBROUTINE shake_roll_4x6_low
885 :
886 : ! **************************************************************************************************
887 : !> \brief ...
888 : !> \param fixd_list ...
889 : !> \param g4x6_list ...
890 : !> \param lg4x6 ...
891 : !> \param first_atom ...
892 : !> \param particle_set ...
893 : !> \param vel ...
894 : !> \param dt ...
895 : !> \par History
896 : !> none
897 : ! **************************************************************************************************
898 702 : SUBROUTINE rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
899 702 : particle_set, vel, dt)
900 :
901 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
902 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
903 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
904 : INTEGER, INTENT(IN) :: first_atom
905 : TYPE(particle_type), POINTER :: particle_set(:)
906 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
907 : REAL(kind=dp), INTENT(in) :: dt
908 :
909 : INTEGER :: iconst, index_a, index_b, index_c, &
910 : index_d
911 : REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, &
912 : imass4, mass
913 : REAL(KIND=dp), DIMENSION(3) :: fc1, fc2, fc3, fc4, r12, r13, r14, r23, &
914 : r24, r34, v12, v13, v14, v23, v24, v34
915 : REAL(KIND=dp), DIMENSION(6, 1) :: bvec
916 : REAL(KIND=dp), DIMENSION(6, 6) :: amat
917 : TYPE(atomic_kind_type), POINTER :: atomic_kind
918 :
919 702 : idt = 1.0_dp/dt
920 702 : dtby2 = dt*.5_dp
921 1404 : DO iconst = 1, SIZE(g4x6_list)
922 702 : IF (g4x6_list(iconst)%restraint%active) CYCLE
923 692 : index_a = g4x6_list(iconst)%a + first_atom - 1
924 692 : index_b = g4x6_list(iconst)%b + first_atom - 1
925 692 : index_c = g4x6_list(iconst)%c + first_atom - 1
926 692 : index_d = g4x6_list(iconst)%d + first_atom - 1
927 2768 : v12(:) = vel(:, index_a) - vel(:, index_b)
928 2768 : v13(:) = vel(:, index_a) - vel(:, index_c)
929 2768 : v14(:) = vel(:, index_a) - vel(:, index_d)
930 2768 : v23(:) = vel(:, index_b) - vel(:, index_c)
931 2768 : v24(:) = vel(:, index_b) - vel(:, index_d)
932 2768 : v34(:) = vel(:, index_c) - vel(:, index_d)
933 :
934 2768 : r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
935 2768 : r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
936 2768 : r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
937 2768 : r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
938 2768 : r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
939 2768 : r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
940 692 : atomic_kind => particle_set(index_a)%atomic_kind
941 692 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
942 692 : imass1 = 1.0_dp/mass
943 692 : atomic_kind => particle_set(index_b)%atomic_kind
944 692 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
945 692 : imass2 = 1.0_dp/mass
946 692 : atomic_kind => particle_set(index_c)%atomic_kind
947 692 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
948 692 : imass3 = 1.0_dp/mass
949 692 : atomic_kind => particle_set(index_d)%atomic_kind
950 692 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
951 692 : imass4 = 1.0_dp/mass
952 2768 : lg4x6(iconst)%fa = -2.0_dp*r12
953 2768 : lg4x6(iconst)%fb = -2.0_dp*r13
954 2768 : lg4x6(iconst)%fc = -2.0_dp*r14
955 2768 : lg4x6(iconst)%fd = -2.0_dp*r23
956 2768 : lg4x6(iconst)%fe = -2.0_dp*r24
957 2768 : lg4x6(iconst)%ff = -2.0_dp*r34
958 : ! Check for fixed atom constraints
959 : CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
960 692 : index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
961 : ! construct matrix
962 2768 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, lg4x6(iconst)%fa)
963 2768 : amat(1, 2) = imass1*DOT_PRODUCT(r12, lg4x6(iconst)%fb)
964 2768 : amat(1, 3) = imass1*DOT_PRODUCT(r12, lg4x6(iconst)%fc)
965 2768 : amat(1, 4) = -imass2*DOT_PRODUCT(r12, lg4x6(iconst)%fd)
966 2768 : amat(1, 5) = -imass2*DOT_PRODUCT(r12, lg4x6(iconst)%fe)
967 692 : amat(1, 6) = 0.0_dp
968 :
969 2768 : amat(2, 1) = imass1*DOT_PRODUCT(r13, lg4x6(iconst)%fa)
970 2768 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, lg4x6(iconst)%fb)
971 2768 : amat(2, 3) = imass1*DOT_PRODUCT(r13, lg4x6(iconst)%fc)
972 2768 : amat(2, 4) = imass3*DOT_PRODUCT(r13, lg4x6(iconst)%fd)
973 692 : amat(2, 5) = 0.0_dp
974 2768 : amat(2, 6) = -imass3*DOT_PRODUCT(r13, lg4x6(iconst)%ff)
975 :
976 2768 : amat(3, 1) = imass1*DOT_PRODUCT(r14, lg4x6(iconst)%fa)
977 2768 : amat(3, 2) = imass1*DOT_PRODUCT(r14, lg4x6(iconst)%fb)
978 2768 : amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r14, lg4x6(iconst)%fc)
979 692 : amat(3, 4) = 0.0_dp
980 2768 : amat(3, 5) = imass4*DOT_PRODUCT(r14, lg4x6(iconst)%fe)
981 2768 : amat(3, 6) = imass4*DOT_PRODUCT(r14, lg4x6(iconst)%ff)
982 :
983 2768 : amat(4, 1) = -imass2*DOT_PRODUCT(r23, lg4x6(iconst)%fa)
984 2768 : amat(4, 2) = imass3*DOT_PRODUCT(r23, lg4x6(iconst)%fb)
985 692 : amat(4, 3) = 0.0_dp
986 2768 : amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r23, lg4x6(iconst)%fd)
987 2768 : amat(4, 5) = imass2*DOT_PRODUCT(r23, lg4x6(iconst)%fe)
988 2768 : amat(4, 6) = -imass3*DOT_PRODUCT(r23, lg4x6(iconst)%ff)
989 :
990 2768 : amat(5, 1) = -imass2*DOT_PRODUCT(r24, lg4x6(iconst)%fa)
991 692 : amat(5, 2) = 0.0_dp
992 2768 : amat(5, 3) = imass4*DOT_PRODUCT(r24, lg4x6(iconst)%fc)
993 2768 : amat(5, 4) = imass2*DOT_PRODUCT(r24, lg4x6(iconst)%fd)
994 2768 : amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r24, lg4x6(iconst)%fe)
995 2768 : amat(5, 6) = imass4*DOT_PRODUCT(r24, lg4x6(iconst)%ff)
996 :
997 692 : amat(6, 1) = 0.0_dp
998 2768 : amat(6, 2) = -imass3*DOT_PRODUCT(r34, lg4x6(iconst)%fb)
999 2768 : amat(6, 3) = imass4*DOT_PRODUCT(r34, lg4x6(iconst)%fc)
1000 2768 : amat(6, 4) = -imass3*DOT_PRODUCT(r34, lg4x6(iconst)%fd)
1001 2768 : amat(6, 5) = imass4*DOT_PRODUCT(r34, lg4x6(iconst)%fe)
1002 2768 : amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r34, lg4x6(iconst)%ff)
1003 :
1004 : ! construct solution vector
1005 2768 : bvec(1, 1) = DOT_PRODUCT(r12, v12)
1006 2768 : bvec(2, 1) = DOT_PRODUCT(r13, v13)
1007 2768 : bvec(3, 1) = DOT_PRODUCT(r14, v14)
1008 2768 : bvec(4, 1) = DOT_PRODUCT(r23, v23)
1009 2768 : bvec(5, 1) = DOT_PRODUCT(r24, v24)
1010 2768 : bvec(6, 1) = DOT_PRODUCT(r34, v34)
1011 5536 : bvec = -bvec*2.0_dp*idt
1012 :
1013 : ! get lambda
1014 692 : CALL solve_system(amat, 6, bvec)
1015 4844 : lg4x6(iconst)%lambda(:) = bvec(:, 1)
1016 :
1017 : fc1 = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + &
1018 : lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb + &
1019 2768 : lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc
1020 : fc2 = -lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + &
1021 : lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + &
1022 2768 : lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe
1023 : fc3 = -lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb - &
1024 : lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + &
1025 2768 : lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff
1026 : fc4 = -lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc - &
1027 : lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe - &
1028 2768 : lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff
1029 2768 : vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
1030 2768 : vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
1031 2768 : vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
1032 4172 : vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
1033 : END DO
1034 702 : END SUBROUTINE rattle_4x6_low
1035 :
1036 : ! **************************************************************************************************
1037 : !> \brief ...
1038 : !> \param fixd_list ...
1039 : !> \param g4x6_list ...
1040 : !> \param lg4x6 ...
1041 : !> \param first_atom ...
1042 : !> \param particle_set ...
1043 : !> \param vel ...
1044 : !> \param r_rattle ...
1045 : !> \param dt ...
1046 : !> \param veps ...
1047 : !> \par History
1048 : !> none
1049 : ! **************************************************************************************************
1050 1024 : SUBROUTINE rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
1051 1024 : particle_set, vel, r_rattle, dt, veps)
1052 :
1053 : TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
1054 : TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
1055 : TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
1056 : INTEGER, INTENT(IN) :: first_atom
1057 : TYPE(particle_type), POINTER :: particle_set(:)
1058 : REAL(KIND=dp), INTENT(INOUT) :: vel(:, :)
1059 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r_rattle
1060 : REAL(kind=dp), INTENT(in) :: dt
1061 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: veps
1062 :
1063 : INTEGER :: iconst, index_a, index_b, index_c, &
1064 : index_d
1065 : REAL(KIND=dp) :: dtby2, idt, imass1, imass2, imass3, &
1066 : imass4, mass
1067 : REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, &
1068 : f_roll5, f_roll6, fc1, fc2, fc3, fc4, &
1069 : r12, r13, r14, r23, r24, r34, v12, &
1070 : v13, v14, v23, v24, v34
1071 : REAL(KIND=dp), DIMENSION(6) :: lambda
1072 : REAL(KIND=dp), DIMENSION(6, 1) :: bvec
1073 : REAL(KIND=dp), DIMENSION(6, 6) :: amat
1074 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1075 :
1076 1024 : idt = 1.0_dp/dt
1077 1024 : dtby2 = dt*.5_dp
1078 2048 : DO iconst = 1, SIZE(g4x6_list)
1079 1024 : IF (g4x6_list(iconst)%restraint%active) CYCLE
1080 1024 : index_a = g4x6_list(iconst)%a + first_atom - 1
1081 1024 : index_b = g4x6_list(iconst)%b + first_atom - 1
1082 1024 : index_c = g4x6_list(iconst)%c + first_atom - 1
1083 1024 : index_d = g4x6_list(iconst)%d + first_atom - 1
1084 4096 : v12(:) = vel(:, index_a) - vel(:, index_b)
1085 4096 : v13(:) = vel(:, index_a) - vel(:, index_c)
1086 4096 : v14(:) = vel(:, index_a) - vel(:, index_d)
1087 4096 : v23(:) = vel(:, index_b) - vel(:, index_c)
1088 4096 : v24(:) = vel(:, index_b) - vel(:, index_d)
1089 4096 : v34(:) = vel(:, index_c) - vel(:, index_d)
1090 :
1091 4096 : r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
1092 4096 : r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
1093 4096 : r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
1094 4096 : r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
1095 4096 : r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
1096 4096 : r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
1097 1024 : atomic_kind => particle_set(index_a)%atomic_kind
1098 1024 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1099 1024 : imass1 = 1.0_dp/mass
1100 1024 : atomic_kind => particle_set(index_b)%atomic_kind
1101 1024 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1102 1024 : imass2 = 1.0_dp/mass
1103 1024 : atomic_kind => particle_set(index_c)%atomic_kind
1104 1024 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1105 1024 : imass3 = 1.0_dp/mass
1106 1024 : atomic_kind => particle_set(index_d)%atomic_kind
1107 1024 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1108 1024 : imass4 = 1.0_dp/mass
1109 4096 : lg4x6(iconst)%fa = -2.0_dp*r12
1110 4096 : lg4x6(iconst)%fb = -2.0_dp*r13
1111 4096 : lg4x6(iconst)%fc = -2.0_dp*r14
1112 4096 : lg4x6(iconst)%fd = -2.0_dp*r23
1113 4096 : lg4x6(iconst)%fe = -2.0_dp*r24
1114 4096 : lg4x6(iconst)%ff = -2.0_dp*r34
1115 : ! Check for fixed atom constraints
1116 : CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
1117 1024 : index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
1118 : ! roll the fc
1119 13312 : f_roll1 = MATMUL(r_rattle, lg4x6(iconst)%fa)
1120 13312 : f_roll2 = MATMUL(r_rattle, lg4x6(iconst)%fb)
1121 13312 : f_roll3 = MATMUL(r_rattle, lg4x6(iconst)%fc)
1122 13312 : f_roll4 = MATMUL(r_rattle, lg4x6(iconst)%fd)
1123 13312 : f_roll5 = MATMUL(r_rattle, lg4x6(iconst)%fe)
1124 13312 : f_roll6 = MATMUL(r_rattle, lg4x6(iconst)%ff)
1125 : ! construct matrix
1126 4096 : amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, f_roll1)
1127 4096 : amat(1, 2) = imass1*DOT_PRODUCT(r12, f_roll2)
1128 4096 : amat(1, 3) = imass1*DOT_PRODUCT(r12, f_roll3)
1129 4096 : amat(1, 4) = -imass2*DOT_PRODUCT(r12, f_roll4)
1130 4096 : amat(1, 5) = -imass2*DOT_PRODUCT(r12, f_roll5)
1131 1024 : amat(1, 6) = 0.0_dp
1132 :
1133 4096 : amat(2, 1) = imass1*DOT_PRODUCT(r13, f_roll1)
1134 4096 : amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, f_roll2)
1135 4096 : amat(2, 3) = imass1*DOT_PRODUCT(r13, f_roll3)
1136 4096 : amat(2, 4) = imass3*DOT_PRODUCT(r13, f_roll4)
1137 1024 : amat(2, 5) = 0.0_dp
1138 4096 : amat(2, 6) = -imass3*DOT_PRODUCT(r13, f_roll6)
1139 :
1140 4096 : amat(3, 1) = imass1*DOT_PRODUCT(r14, f_roll1)
1141 4096 : amat(3, 2) = imass1*DOT_PRODUCT(r14, f_roll2)
1142 4096 : amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r14, f_roll3)
1143 1024 : amat(3, 4) = 0.0_dp
1144 4096 : amat(3, 5) = imass4*DOT_PRODUCT(r14, f_roll5)
1145 4096 : amat(3, 6) = imass4*DOT_PRODUCT(r14, f_roll6)
1146 :
1147 4096 : amat(4, 1) = -imass2*DOT_PRODUCT(r23, f_roll1)
1148 4096 : amat(4, 2) = imass3*DOT_PRODUCT(r23, f_roll2)
1149 1024 : amat(4, 3) = 0.0_dp
1150 4096 : amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r23, f_roll4)
1151 4096 : amat(4, 5) = imass2*DOT_PRODUCT(r23, f_roll5)
1152 4096 : amat(4, 6) = -imass3*DOT_PRODUCT(r23, f_roll6)
1153 :
1154 4096 : amat(5, 1) = -imass2*DOT_PRODUCT(r24, f_roll1)
1155 1024 : amat(5, 2) = 0.0_dp
1156 4096 : amat(5, 3) = imass4*DOT_PRODUCT(r24, f_roll3)
1157 4096 : amat(5, 4) = imass2*DOT_PRODUCT(r24, f_roll4)
1158 4096 : amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r24, f_roll5)
1159 4096 : amat(5, 6) = imass4*DOT_PRODUCT(r24, f_roll6)
1160 :
1161 1024 : amat(6, 1) = 0.0_dp
1162 4096 : amat(6, 2) = -imass3*DOT_PRODUCT(r34, f_roll2)
1163 4096 : amat(6, 3) = imass4*DOT_PRODUCT(r34, f_roll3)
1164 4096 : amat(6, 4) = -imass3*DOT_PRODUCT(r34, f_roll4)
1165 4096 : amat(6, 5) = imass4*DOT_PRODUCT(r34, f_roll5)
1166 4096 : amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r34, f_roll6)
1167 :
1168 : ! construct solution vector
1169 22528 : bvec(1, 1) = DOT_PRODUCT(r12, v12 + MATMUL(veps, r12))
1170 22528 : bvec(2, 1) = DOT_PRODUCT(r13, v13 + MATMUL(veps, r13))
1171 22528 : bvec(3, 1) = DOT_PRODUCT(r14, v14 + MATMUL(veps, r14))
1172 22528 : bvec(4, 1) = DOT_PRODUCT(r23, v23 + MATMUL(veps, r23))
1173 22528 : bvec(5, 1) = DOT_PRODUCT(r24, v24 + MATMUL(veps, r24))
1174 22528 : bvec(6, 1) = DOT_PRODUCT(r34, v34 + MATMUL(veps, r34))
1175 8192 : bvec = -bvec*2.0_dp*idt
1176 :
1177 : ! get lambda
1178 1024 : CALL solve_system(amat, 6, bvec)
1179 1024 : lambda(:) = bvec(:, 1)
1180 7168 : lg4x6(iconst)%lambda(:) = lambda
1181 :
1182 : fc1 = lambda(1)*f_roll1 + &
1183 : lambda(2)*f_roll2 + &
1184 4096 : lambda(3)*f_roll3
1185 : fc2 = -lambda(1)*f_roll1 + &
1186 : lambda(4)*f_roll4 + &
1187 4096 : lambda(5)*f_roll5
1188 : fc3 = -lambda(2)*f_roll2 - &
1189 : lambda(4)*f_roll4 + &
1190 4096 : lambda(6)*f_roll6
1191 : fc4 = -lambda(3)*f_roll3 - &
1192 : lambda(5)*f_roll5 - &
1193 4096 : lambda(6)*f_roll6
1194 4096 : vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
1195 4096 : vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
1196 4096 : vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
1197 6144 : vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
1198 : END DO
1199 1024 : END SUBROUTINE rattle_roll_4x6_low
1200 :
1201 6144 : END MODULE constraint_4x6
|
