(git:374b731)
Loading...
Searching...
No Matches
constraint_3x3.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> none
11! **************************************************************************************************
13
17 USE kinds, ONLY: dp
23 USE molecule_types, ONLY: get_molecule,&
28#include "./base/base_uses.f90"
29
30 IMPLICIT NONE
31
32 PRIVATE
33 PUBLIC :: shake_3x3_int, &
41
42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_3x3'
43
44CONTAINS
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 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
70 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
71 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
72 TYPE(molecule_kind_type), POINTER :: molecule_kind
73
74 NULLIFY (fixd_list)
75 molecule_kind => molecule%molecule_kind
76 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
77 g3x3_list=g3x3_list, fixd_list=fixd_list)
78 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 particle_set, pos, vel, dt, ishake, max_sigma)
82
83 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 SUBROUTINE shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, &
100 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
112 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
113 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
114 TYPE(molecule_kind_type), POINTER :: molecule_kind
115
116 NULLIFY (fixd_list)
117 molecule_kind => molecule%molecule_kind
118 CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
119 g3x3_list=g3x3_list, fixd_list=fixd_list)
120 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 particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
124
125 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 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
149 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
150 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
151 TYPE(molecule_kind_type), POINTER :: molecule_kind
152
153 NULLIFY (fixd_list)
154 molecule_kind => molecule%molecule_kind
155 CALL get_molecule_kind(molecule_kind, &
156 g3x3_list=g3x3_list, fixd_list=fixd_list)
157 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 particle_set, vel, r_rattle, dt, veps)
161
162 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 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 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
182 TYPE(g3x3_constraint_type), POINTER :: g3x3_list(:)
183 TYPE(local_g3x3_constraint_type), POINTER :: lg3x3(:)
184 TYPE(molecule_kind_type), POINTER :: molecule_kind
185
186 NULLIFY (fixd_list)
187 molecule_kind => molecule%molecule_kind
188 CALL get_molecule_kind(molecule_kind, &
189 g3x3_list=g3x3_list, fixd_list=fixd_list)
190 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 particle_set, vel, dt)
194
195 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 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 first_atom = 1
225 ng3x3 = gci%ng3x3
226 g3x3_list => gci%g3x3_list
227 fixd_list => gci%fixd_list
228 lg3x3 => gci%lg3x3
229 ! Real Shake
230 CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
231 particle_set, pos, vel, dt, ishake, max_sigma)
232
233 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 SUBROUTINE shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, &
250 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 first_atom = 1
266 ng3x3 = gci%ng3x3
267 g3x3_list => gci%g3x3_list
268 fixd_list => gci%fixd_list
269 lg3x3 => gci%lg3x3
270 ! Real Shake
271 CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
272 particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
273
274 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 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 first_atom = 1
302 g3x3_list => gci%g3x3_list
303 fixd_list => gci%fixd_list
304 lg3x3 => gci%lg3x3
305 ! Real Rattle
306 CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
307 particle_set, vel, r_rattle, dt, veps)
308
309 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 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 first_atom = 1
333 g3x3_list => gci%g3x3_list
334 fixd_list => gci%fixd_list
335 lg3x3 => gci%lg3x3
336 ! Real Rattle
337 CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
338 particle_set, vel, dt)
339
340 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 SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
359 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 dtsqby2 = dt*dt*.5_dp
381 idtsq = 1.0_dp/dt/dt
382 dtby2 = dt*.5_dp
383 DO iconst = 1, ng3x3
384 IF (g3x3_list(iconst)%restraint%active) cycle
385 index_a = g3x3_list(iconst)%a + first_atom - 1
386 index_b = g3x3_list(iconst)%b + first_atom - 1
387 index_c = g3x3_list(iconst)%c + first_atom - 1
388 IF (ishake == 1) THEN
389 r0_12(:) = pos(:, index_a) - pos(:, index_b)
390 r0_13(:) = pos(:, index_a) - pos(:, index_c)
391 r0_23(:) = pos(:, index_b) - pos(:, index_c)
392 atomic_kind => particle_set(index_a)%atomic_kind
393 imass1 = 1.0_dp/atomic_kind%mass
394 atomic_kind => particle_set(index_b)%atomic_kind
395 imass2 = 1.0_dp/atomic_kind%mass
396 atomic_kind => particle_set(index_c)%atomic_kind
397 imass3 = 1.0_dp/atomic_kind%mass
398 lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
399 lg3x3(iconst)%rb_old)
400 lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
401 lg3x3(iconst)%rc_old)
402 lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
403 lg3x3(iconst)%rc_old)
404 ! Check for fixed atom constraints
405 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
406 index_a, index_b, index_c, fixd_list, lg3x3(iconst))
407 ! construct matrix
408 amat(1, 1) = (imass1 + imass2)*dot_product(r0_12, lg3x3(iconst)%fa)
409 amat(1, 2) = imass1*dot_product(r0_12, lg3x3(iconst)%fb)
410 amat(1, 3) = -imass2*dot_product(r0_12, lg3x3(iconst)%fc)
411 amat(2, 1) = imass1*dot_product(r0_13, lg3x3(iconst)%fa)
412 amat(2, 2) = (imass1 + imass3)*dot_product(r0_13, lg3x3(iconst)%fb)
413 amat(2, 3) = imass3*dot_product(r0_13, lg3x3(iconst)%fc)
414 amat(3, 1) = -imass2*dot_product(r0_23, lg3x3(iconst)%fa)
415 amat(3, 2) = imass3*dot_product(r0_23, lg3x3(iconst)%fb)
416 amat(3, 3) = (imass3 + imass2)*dot_product(r0_23, lg3x3(iconst)%fc)
417 ! Store values
418 lg3x3(iconst)%r0_12 = r0_12
419 lg3x3(iconst)%r0_13 = r0_13
420 lg3x3(iconst)%r0_23 = r0_23
421 lg3x3(iconst)%amat = amat
422 lg3x3(iconst)%lambda_old = 0.0_dp
423 lg3x3(iconst)%del_lambda = 0.0_dp
424 ELSE
425 ! Retrieve values
426 r0_12 = lg3x3(iconst)%r0_12
427 r0_13 = lg3x3(iconst)%r0_13
428 r0_23 = lg3x3(iconst)%r0_23
429 amat = lg3x3(iconst)%amat
430 imass1 = lg3x3(iconst)%imass1
431 imass2 = lg3x3(iconst)%imass2
432 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 lg3x3(iconst)%lambda(3)*imass2*lg3x3(iconst)%fc
438 bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
439 - 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 lg3x3(iconst)%lambda(3)*imass3*lg3x3(iconst)%fc
443 bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
444 - 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 lg3x3(iconst)%lambda(3)*(imass2 + imass3)*lg3x3(iconst)%fc
448 bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
449 - dtsqby2*dtsqby2*dot_product(vec, vec) - dot_product(r0_23, r0_23)
450 bvec = bvec*idtsq
451 ! get lambda
452 atemp = amat
453 CALL solve_system(atemp, 3, bvec)
454 lg3x3(iconst)%lambda(:) = bvec(:, 1)
455 lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
456 lg3x3(iconst)%lambda_old(:)
457 lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
458 fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
459 lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
460 fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
461 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
462 fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
463 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
464 r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
465 r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
466 r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
467 v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
468 v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
469 v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
470 r12 = r1 - r2
471 r13 = r1 - r3
472 r23 = r2 - r3
473
474 ! compute the tolerance:
475 sigma = dot_product(r12, r12) - g3x3_list(iconst)%dab* &
476 g3x3_list(iconst)%dab
477 max_sigma = max(max_sigma, abs(sigma))
478 sigma = dot_product(r13, r13) - g3x3_list(iconst)%dac* &
479 g3x3_list(iconst)%dac
480 max_sigma = max(max_sigma, abs(sigma))
481 sigma = dot_product(r23, r23) - g3x3_list(iconst)%dbc* &
482 g3x3_list(iconst)%dbc
483 max_sigma = max(max_sigma, abs(sigma))
484 ! update positions with full multiplier
485 pos(:, index_a) = r1(:)
486 pos(:, index_b) = r2(:)
487 pos(:, index_c) = r3(:)
488
489 ! update velocites with full multiplier
490 vel(:, index_a) = v1(:)
491 vel(:, index_b) = v2(:)
492 vel(:, index_c) = v3(:)
493 END DO
494
495 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 SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
516 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 dtsqby2 = dt*dt*.5_dp
540 idtsq = 1.0_dp/dt/dt
541 dtby2 = dt*.5_dp
542 DO iconst = 1, ng3x3
543 IF (g3x3_list(iconst)%restraint%active) cycle
544 index_a = g3x3_list(iconst)%a + first_atom - 1
545 index_b = g3x3_list(iconst)%b + first_atom - 1
546 index_c = g3x3_list(iconst)%c + first_atom - 1
547 IF (ishake == 1) THEN
548 r0_12(:) = pos(:, index_a) - pos(:, index_b)
549 r0_13(:) = pos(:, index_a) - pos(:, index_c)
550 r0_23(:) = pos(:, index_b) - pos(:, index_c)
551 atomic_kind => particle_set(index_a)%atomic_kind
552 imass1 = 1.0_dp/atomic_kind%mass
553 atomic_kind => particle_set(index_b)%atomic_kind
554 imass2 = 1.0_dp/atomic_kind%mass
555 atomic_kind => particle_set(index_c)%atomic_kind
556 imass3 = 1.0_dp/atomic_kind%mass
557 lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
558 lg3x3(iconst)%rb_old)
559 lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
560 lg3x3(iconst)%rc_old)
561 lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
562 lg3x3(iconst)%rc_old)
563 ! Check for fixed atom constraints
564 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
565 index_a, index_b, index_c, fixd_list, lg3x3(iconst))
566 ! rotate fconst:
567 f_roll1 = matmul(r_shake, lg3x3(iconst)%fa)
568 f_roll2 = matmul(r_shake, lg3x3(iconst)%fb)
569 f_roll3 = matmul(r_shake, lg3x3(iconst)%fc)
570 ! construct matrix
571 amat(1, 1) = (imass1 + imass2)*dot_product(r0_12, f_roll1)
572 amat(1, 2) = imass1*dot_product(r0_12, f_roll2)
573 amat(1, 3) = -imass2*dot_product(r0_12, f_roll3)
574 amat(2, 1) = imass1*dot_product(r0_13, f_roll1)
575 amat(2, 2) = (imass1 + imass3)*dot_product(r0_13, f_roll2)
576 amat(2, 3) = imass3*dot_product(r0_13, f_roll3)
577 amat(3, 1) = -imass2*dot_product(r0_23, f_roll1)
578 amat(3, 2) = imass3*dot_product(r0_23, f_roll2)
579 amat(3, 3) = (imass3 + imass2)*dot_product(r0_23, f_roll3)
580 ! Store values
581 lg3x3(iconst)%r0_12 = r0_12
582 lg3x3(iconst)%r0_13 = r0_13
583 lg3x3(iconst)%r0_23 = r0_23
584 lg3x3(iconst)%f_roll1 = f_roll1
585 lg3x3(iconst)%f_roll2 = f_roll2
586 lg3x3(iconst)%f_roll3 = f_roll3
587 lg3x3(iconst)%amat = amat
588 lg3x3(iconst)%lambda_old = 0.0_dp
589 lg3x3(iconst)%del_lambda = 0.0_dp
590 ELSE
591 ! Retrieve values
592 r0_12 = lg3x3(iconst)%r0_12
593 r0_13 = lg3x3(iconst)%r0_13
594 r0_23 = lg3x3(iconst)%r0_23
595 f_roll1 = lg3x3(iconst)%f_roll1
596 f_roll2 = lg3x3(iconst)%f_roll2
597 f_roll3 = lg3x3(iconst)%f_roll3
598 amat = lg3x3(iconst)%amat
599 imass1 = lg3x3(iconst)%imass1
600 imass2 = lg3x3(iconst)%imass2
601 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 lg3x3(iconst)%lambda(3)*imass2*f_roll3
607 bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
608 - 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 lg3x3(iconst)%lambda(3)*imass3*f_roll3
612 bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
613 - 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 lg3x3(iconst)%lambda(3)*(imass2 + imass3)*f_roll3
617 bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
618 - dtsqby2*dtsqby2*dot_product(vec, vec) - dot_product(r0_23, r0_23)
619 bvec = bvec*idtsq
620
621 ! get lambda
622 atemp = amat
623 CALL solve_system(atemp, 3, bvec)
624 lg3x3(iconst)%lambda(:) = bvec(:, 1)
625 lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
626 lg3x3(iconst)%lambda_old(:)
627 lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
628
629 fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
630 lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
631 fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
632 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
633 fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
634 lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
635 r1(:) = pos(:, index_a) + imass1*dtsqby2*matmul(r_shake, fc1)
636 r2(:) = pos(:, index_b) + imass2*dtsqby2*matmul(r_shake, fc2)
637 r3(:) = pos(:, index_c) + imass3*dtsqby2*matmul(r_shake, fc3)
638 v1(:) = vel(:, index_a) + imass1*dtby2*matmul(v_shake, fc1)
639 v2(:) = vel(:, index_b) + imass2*dtby2*matmul(v_shake, fc2)
640 v3(:) = vel(:, index_c) + imass3*dtby2*matmul(v_shake, fc3)
641 r12 = r1 - r2
642 r13 = r1 - r3
643 r23 = r2 - r3
644
645 ! compute the tolerance:
646 sigma = dot_product(r12, r12) - g3x3_list(iconst)%dab* &
647 g3x3_list(iconst)%dab
648 max_sigma = max(max_sigma, abs(sigma))
649 sigma = dot_product(r13, r13) - g3x3_list(iconst)%dac* &
650 g3x3_list(iconst)%dac
651 max_sigma = max(max_sigma, abs(sigma))
652 sigma = dot_product(r23, r23) - g3x3_list(iconst)%dbc* &
653 g3x3_list(iconst)%dbc
654 max_sigma = max(max_sigma, abs(sigma))
655
656 ! update positions with full multiplier
657 pos(:, index_a) = r1(:)
658 pos(:, index_b) = r2(:)
659 pos(:, index_c) = r3(:)
660
661 ! update velocites with full multiplier
662 vel(:, index_a) = v1(:)
663 vel(:, index_b) = v2(:)
664 vel(:, index_c) = v3(:)
665 END DO
666 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 SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
683 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 idt = 1.0_dp/dt
704 dtby2 = dt*.5_dp
705 DO iconst = 1, SIZE(g3x3_list)
706 IF (g3x3_list(iconst)%restraint%active) cycle
707 index_a = g3x3_list(iconst)%a + first_atom - 1
708 index_b = g3x3_list(iconst)%b + first_atom - 1
709 index_c = g3x3_list(iconst)%c + first_atom - 1
710 v12(:) = vel(:, index_a) - vel(:, index_b)
711 v13(:) = vel(:, index_a) - vel(:, index_c)
712 v23(:) = vel(:, index_b) - vel(:, index_c)
713 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
714 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
715 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
716 atomic_kind => particle_set(index_a)%atomic_kind
717 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
718 imass1 = 1.0_dp/mass
719 atomic_kind => particle_set(index_b)%atomic_kind
720 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
721 imass2 = 1.0_dp/mass
722 atomic_kind => particle_set(index_c)%atomic_kind
723 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
724 imass3 = 1.0_dp/mass
725 lg3x3(iconst)%fa = -2.0_dp*r12
726 lg3x3(iconst)%fb = -2.0_dp*r13
727 lg3x3(iconst)%fc = -2.0_dp*r23
728 ! Check for fixed atom constraints
729 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
730 index_a, index_b, index_c, fixd_list, lg3x3(iconst))
731 ! roll the fc
732 f_roll1 = matmul(r_rattle, lg3x3(iconst)%fa)
733 f_roll2 = matmul(r_rattle, lg3x3(iconst)%fb)
734 f_roll3 = matmul(r_rattle, lg3x3(iconst)%fc)
735
736 ! construct matrix
737 amat(1, 1) = (imass1 + imass2)*dot_product(r12, f_roll1)
738 amat(1, 2) = imass1*dot_product(r12, f_roll2)
739 amat(1, 3) = -imass2*dot_product(r12, f_roll3)
740 amat(2, 1) = imass1*dot_product(r13, f_roll1)
741 amat(2, 2) = (imass1 + imass3)*dot_product(r13, f_roll2)
742 amat(2, 3) = imass3*dot_product(r13, f_roll3)
743 amat(3, 1) = -imass2*dot_product(r23, f_roll1)
744 amat(3, 2) = imass3*dot_product(r23, f_roll2)
745 amat(3, 3) = (imass2 + imass3)*dot_product(r23, f_roll3)
746
747 ! construct solution vector
748 bvec(1, 1) = dot_product(r12, v12 + matmul(veps, r12))
749 bvec(2, 1) = dot_product(r13, v13 + matmul(veps, r13))
750 bvec(3, 1) = dot_product(r23, v23 + matmul(veps, r23))
751 bvec = -bvec*2.0_dp*idt
752
753 ! get lambda
754 CALL solve_system(amat, 3, bvec)
755 lambda(:) = bvec(:, 1)
756 lg3x3(iconst)%lambda(:) = lambda
757
758 fc1 = lambda(1)*f_roll1 + &
759 lambda(2)*f_roll2
760 fc2 = -lambda(1)*f_roll1 + &
761 lambda(3)*f_roll3
762 fc3 = -lambda(2)*f_roll2 - &
763 lambda(3)*f_roll3
764 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
765 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
766 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
767 END DO
768 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 SUBROUTINE rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
783 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 idt = 1.0_dp/dt
802 dtby2 = dt*.5_dp
803 DO iconst = 1, SIZE(g3x3_list)
804 IF (g3x3_list(iconst)%restraint%active) cycle
805 index_a = g3x3_list(iconst)%a + first_atom - 1
806 index_b = g3x3_list(iconst)%b + first_atom - 1
807 index_c = g3x3_list(iconst)%c + first_atom - 1
808 v12(:) = vel(:, index_a) - vel(:, index_b)
809 v13(:) = vel(:, index_a) - vel(:, index_c)
810 v23(:) = vel(:, index_b) - vel(:, index_c)
811 r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
812 r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
813 r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
814 atomic_kind => particle_set(index_a)%atomic_kind
815 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
816 imass1 = 1.0_dp/mass
817 atomic_kind => particle_set(index_b)%atomic_kind
818 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
819 imass2 = 1.0_dp/mass
820 atomic_kind => particle_set(index_c)%atomic_kind
821 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
822 imass3 = 1.0_dp/mass
823 lg3x3(iconst)%fa = -2.0_dp*r12
824 lg3x3(iconst)%fb = -2.0_dp*r13
825 lg3x3(iconst)%fc = -2.0_dp*r23
826 ! Check for fixed atom constraints
827 CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
828 index_a, index_b, index_c, fixd_list, lg3x3(iconst))
829 ! construct matrix
830 amat(1, 1) = (imass1 + imass2)*dot_product(r12, lg3x3(iconst)%fa)
831 amat(1, 2) = imass1*dot_product(r12, lg3x3(iconst)%fb)
832 amat(1, 3) = -imass2*dot_product(r12, lg3x3(iconst)%fc)
833 amat(2, 1) = imass1*dot_product(r13, lg3x3(iconst)%fa)
834 amat(2, 2) = (imass1 + imass3)*dot_product(r13, lg3x3(iconst)%fb)
835 amat(2, 3) = imass3*dot_product(r13, lg3x3(iconst)%fc)
836 amat(3, 1) = -imass2*dot_product(r23, lg3x3(iconst)%fa)
837 amat(3, 2) = imass3*dot_product(r23, lg3x3(iconst)%fb)
838 amat(3, 3) = (imass2 + imass3)*dot_product(r23, lg3x3(iconst)%fc)
839
840 ! construct solution vector
841 bvec(1, 1) = dot_product(r12, v12)
842 bvec(2, 1) = dot_product(r13, v13)
843 bvec(3, 1) = dot_product(r23, v23)
844 bvec = -bvec*2.0_dp*idt
845
846 ! get lambda
847 CALL solve_system(amat, 3, bvec)
848 lg3x3(iconst)%lambda(:) = bvec(:, 1)
849
850 fc1 = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
851 lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb
852 fc2 = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
853 lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
854 fc3 = -lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb - &
855 lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
856 vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
857 vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
858 vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
859 END DO
860 END SUBROUTINE rattle_3x3_low
861
862END MODULE constraint_3x3
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public shake_3x3_ext(gci, particle_set, pos, vel, dt, ishake, max_sigma)
...
subroutine, public rattle_roll_3x3_ext(gci, particle_set, vel, r_rattle, dt, veps)
...
subroutine, public shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
...
subroutine, public rattle_3x3_int(molecule, particle_set, vel, dt)
...
subroutine, public shake_3x3_int(molecule, particle_set, pos, vel, dt, ishake, max_sigma)
...
subroutine, public rattle_3x3_ext(gci, particle_set, vel, dt)
...
subroutine, public shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
...
subroutine, public rattle_roll_3x3_int(molecule, particle_set, vel, r_rattle, dt, veps)
...
subroutine, public check_fixed_atom_cns_g3x3(imass1, imass2, imass3, index_a, index_b, index_c, fixd_list, lg3x3)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Provides interfaces to LAPACK routines for factorisation and linear system solving.
subroutine, public solve_system(matrix, mysize, eigenvectors)
...
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Define the data structure for the molecule information.
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Define the data structure for the particle information.
Provides all information about an atomic kind.