(git:1a29073)
constraint_4x6.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 ! **************************************************************************************************
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  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  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
69  TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
70  TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
71  TYPE(molecule_kind_type), POINTER :: molecule_kind
72 
73  NULLIFY (fixd_list)
74  molecule_kind => molecule%molecule_kind
75  CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
76  fixd_list=fixd_list, g4x6_list=g4x6_list)
77  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  particle_set, pos, vel, dt, ishake, max_sigma)
81 
82  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  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  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
109  TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
110  TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
111  TYPE(molecule_kind_type), POINTER :: molecule_kind
112 
113  NULLIFY (fixd_list)
114  molecule_kind => molecule%molecule_kind
115  CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
116  fixd_list=fixd_list, g4x6_list=g4x6_list)
117  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  particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
121 
122  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  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  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
141  TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
142  TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
143  TYPE(molecule_kind_type), POINTER :: molecule_kind
144 
145  NULLIFY (fixd_list)
146  molecule_kind => molecule%molecule_kind
147  CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
148  fixd_list=fixd_list, g4x6_list=g4x6_list)
149  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  particle_set, vel, dt)
153 
154  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  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  TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list
177  TYPE(g4x6_constraint_type), POINTER :: g4x6_list(:)
178  TYPE(local_g4x6_constraint_type), POINTER :: lg4x6(:)
179  TYPE(molecule_kind_type), POINTER :: molecule_kind
180 
181  NULLIFY (fixd_list)
182  molecule_kind => molecule%molecule_kind
183  CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
184  fixd_list=fixd_list, g4x6_list=g4x6_list)
185  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  particle_set, vel, r_rattle, dt, veps)
189 
190  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  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  first_atom = 1
220  ng4x6 = gci%ng4x6
221  fixd_list => gci%fixd_list
222  g4x6_list => gci%g4x6_list
223  lg4x6 => gci%lg4x6
224  ! Real Shake
225  CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
226  particle_set, pos, vel, dt, ishake, max_sigma)
227 
228  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  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  first_atom = 1
260  ng4x6 = gci%ng4x6
261  fixd_list => gci%fixd_list
262  g4x6_list => gci%g4x6_list
263  lg4x6 => gci%lg4x6
264  ! Real Shake
265  CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
266  particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
267 
268  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  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  first_atom = 1
292  ng4x6 = gci%ng4x6
293  fixd_list => gci%fixd_list
294  g4x6_list => gci%g4x6_list
295  lg4x6 => gci%lg4x6
296  ! Real Rattle
297  CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
298  particle_set, vel, dt)
299 
300  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  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  first_atom = 1
328  ng4x6 = gci%ng4x6
329  fixd_list => gci%fixd_list
330  g4x6_list => gci%g4x6_list
331  lg4x6 => gci%lg4x6
332  ! Real Rattle
333  CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
334  particle_set, vel, r_rattle, dt, veps)
335 
336  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  SUBROUTINE shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
355  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  dtsqby2 = dt*dt*.5_dp
379  idtsq = 1.0_dp/dt/dt
380  dtby2 = dt*.5_dp
381  DO iconst = 1, ng4x6
382  IF (g4x6_list(iconst)%restraint%active) cycle
383  index_a = g4x6_list(iconst)%a + first_atom - 1
384  index_b = g4x6_list(iconst)%b + first_atom - 1
385  index_c = g4x6_list(iconst)%c + first_atom - 1
386  index_d = g4x6_list(iconst)%d + first_atom - 1
387  IF (ishake == 1) THEN
388  r0_12(:) = pos(:, index_a) - pos(:, index_b)
389  r0_13(:) = pos(:, index_a) - pos(:, index_c)
390  r0_14(:) = pos(:, index_a) - pos(:, index_d)
391  r0_23(:) = pos(:, index_b) - pos(:, index_c)
392  r0_24(:) = pos(:, index_b) - pos(:, index_d)
393  r0_34(:) = pos(:, index_c) - pos(:, index_d)
394  atomic_kind => particle_set(index_a)%atomic_kind
395  imass1 = 1.0_dp/atomic_kind%mass
396  atomic_kind => particle_set(index_b)%atomic_kind
397  imass2 = 1.0_dp/atomic_kind%mass
398  atomic_kind => particle_set(index_c)%atomic_kind
399  imass3 = 1.0_dp/atomic_kind%mass
400  atomic_kind => particle_set(index_d)%atomic_kind
401  imass4 = 1.0_dp/atomic_kind%mass
402  lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
403  lg4x6(iconst)%rb_old)
404  lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
405  lg4x6(iconst)%rc_old)
406  lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
407  lg4x6(iconst)%rd_old)
408  lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
409  lg4x6(iconst)%rc_old)
410  lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
411  lg4x6(iconst)%rd_old)
412  lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
413  lg4x6(iconst)%rd_old)
414  ! Check for fixed atom constraints
415  CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
416  index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
417  ! construct matrix
418  amat(1, 1) = (imass1 + imass2)*dot_product(r0_12, lg4x6(iconst)%fa)
419  amat(1, 2) = imass1*dot_product(r0_12, lg4x6(iconst)%fb)
420  amat(1, 3) = imass1*dot_product(r0_12, lg4x6(iconst)%fc)
421  amat(1, 4) = -imass2*dot_product(r0_12, lg4x6(iconst)%fd)
422  amat(1, 5) = -imass2*dot_product(r0_12, lg4x6(iconst)%fe)
423  amat(1, 6) = 0.0_dp
424 
425  amat(2, 1) = imass1*dot_product(r0_13, lg4x6(iconst)%fa)
426  amat(2, 2) = (imass1 + imass3)*dot_product(r0_13, lg4x6(iconst)%fb)
427  amat(2, 3) = imass1*dot_product(r0_13, lg4x6(iconst)%fc)
428  amat(2, 4) = imass3*dot_product(r0_13, lg4x6(iconst)%fd)
429  amat(2, 5) = 0.0_dp
430  amat(2, 6) = -imass3*dot_product(r0_13, lg4x6(iconst)%ff)
431 
432  amat(3, 1) = imass1*dot_product(r0_14, lg4x6(iconst)%fa)
433  amat(3, 2) = imass1*dot_product(r0_14, lg4x6(iconst)%fb)
434  amat(3, 3) = (imass1 + imass4)*dot_product(r0_14, lg4x6(iconst)%fc)
435  amat(3, 4) = 0.0_dp
436  amat(3, 5) = imass4*dot_product(r0_14, lg4x6(iconst)%fe)
437  amat(3, 6) = imass4*dot_product(r0_14, lg4x6(iconst)%ff)
438 
439  amat(4, 1) = -imass2*dot_product(r0_23, lg4x6(iconst)%fa)
440  amat(4, 2) = imass3*dot_product(r0_23, lg4x6(iconst)%fb)
441  amat(4, 3) = 0.0_dp
442  amat(4, 4) = (imass3 + imass2)*dot_product(r0_23, lg4x6(iconst)%fd)
443  amat(4, 5) = imass2*dot_product(r0_23, lg4x6(iconst)%fe)
444  amat(4, 6) = -imass3*dot_product(r0_23, lg4x6(iconst)%ff)
445 
446  amat(5, 1) = -imass2*dot_product(r0_24, lg4x6(iconst)%fa)
447  amat(5, 2) = 0.0_dp
448  amat(5, 3) = imass4*dot_product(r0_24, lg4x6(iconst)%fc)
449  amat(5, 4) = imass2*dot_product(r0_24, lg4x6(iconst)%fd)
450  amat(5, 5) = (imass4 + imass2)*dot_product(r0_24, lg4x6(iconst)%fe)
451  amat(5, 6) = imass4*dot_product(r0_24, lg4x6(iconst)%ff)
452 
453  amat(6, 1) = 0.0_dp
454  amat(6, 2) = -imass3*dot_product(r0_34, lg4x6(iconst)%fb)
455  amat(6, 3) = imass4*dot_product(r0_34, lg4x6(iconst)%fc)
456  amat(6, 4) = -imass3*dot_product(r0_34, lg4x6(iconst)%fd)
457  amat(6, 5) = imass4*dot_product(r0_34, lg4x6(iconst)%fe)
458  amat(6, 6) = (imass3 + imass4)*dot_product(r0_34, lg4x6(iconst)%ff)
459  ! Store values
460  lg4x6(iconst)%r0_12 = r0_12
461  lg4x6(iconst)%r0_13 = r0_13
462  lg4x6(iconst)%r0_14 = r0_14
463  lg4x6(iconst)%r0_23 = r0_23
464  lg4x6(iconst)%r0_24 = r0_24
465  lg4x6(iconst)%r0_34 = r0_34
466  lg4x6(iconst)%amat = amat
467  lg4x6(iconst)%imass1 = imass1
468  lg4x6(iconst)%imass2 = imass2
469  lg4x6(iconst)%imass3 = imass3
470  lg4x6(iconst)%imass4 = imass4
471  lg4x6(iconst)%lambda_old = 0.0_dp
472  lg4x6(iconst)%del_lambda = 0.0_dp
473  ELSE
474  ! Retrieve values
475  r0_12 = lg4x6(iconst)%r0_12
476  r0_13 = lg4x6(iconst)%r0_13
477  r0_14 = lg4x6(iconst)%r0_14
478  r0_23 = lg4x6(iconst)%r0_23
479  r0_24 = lg4x6(iconst)%r0_24
480  r0_34 = lg4x6(iconst)%r0_34
481  amat = lg4x6(iconst)%amat
482  imass1 = lg4x6(iconst)%imass1
483  imass2 = lg4x6(iconst)%imass2
484  imass3 = lg4x6(iconst)%imass3
485  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  lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe
494  bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
495  - 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  lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
501  bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
502  - 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  lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
508  bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
509  - 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  lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
515  bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
516  - 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  lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
522  bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
523  - 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  lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe
529  bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
530  - dtsqby2*dtsqby2*dot_product(vec, vec) - dot_product(r0_34, r0_34)
531  bvec = bvec*idtsq
532 
533  ! get lambda
534  atemp = amat
535  CALL solve_system(atemp, 6, bvec)
536  lg4x6(iconst)%lambda(:) = bvec(:, 1)
537  lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
538  lg4x6(iconst)%lambda_old(:)
539  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  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  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  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  lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
553  r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
554  r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
555  r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
556  r4(:) = pos(:, index_d) + imass4*dtsqby2*fc4(:)
557  v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
558  v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
559  v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
560  v4(:) = vel(:, index_d) + imass4*dtby2*fc4(:)
561  r12 = r1 - r2
562  r13 = r1 - r3
563  r14 = r1 - r4
564  r23 = r2 - r3
565  r24 = r2 - r4
566  r34 = r3 - r4
567  ! compute the tolerance:
568  sigma = dot_product(r12, r12) - g4x6_list(iconst)%dab* &
569  g4x6_list(iconst)%dab
570  max_sigma = max(max_sigma, abs(sigma))
571  sigma = dot_product(r13, r13) - g4x6_list(iconst)%dac* &
572  g4x6_list(iconst)%dac
573  max_sigma = max(max_sigma, abs(sigma))
574  sigma = dot_product(r14, r14) - g4x6_list(iconst)%dad* &
575  g4x6_list(iconst)%dad
576  max_sigma = max(max_sigma, abs(sigma))
577  sigma = dot_product(r23, r23) - g4x6_list(iconst)%dbc* &
578  g4x6_list(iconst)%dbc
579  max_sigma = max(max_sigma, abs(sigma))
580  sigma = dot_product(r24, r24) - g4x6_list(iconst)%dbd* &
581  g4x6_list(iconst)%dbd
582  max_sigma = max(max_sigma, abs(sigma))
583  sigma = dot_product(r34, r34) - g4x6_list(iconst)%dcd* &
584  g4x6_list(iconst)%dcd
585  max_sigma = max(max_sigma, abs(sigma))
586 
587  ! update positions with full multiplier
588  pos(:, index_a) = r1(:)
589  pos(:, index_b) = r2(:)
590  pos(:, index_c) = r3(:)
591  pos(:, index_d) = r4(:)
592 
593  ! update velocites with full multiplier
594  vel(:, index_a) = v1(:)
595  vel(:, index_b) = v2(:)
596  vel(:, index_c) = v3(:)
597  vel(:, index_d) = v4(:)
598  END DO
599  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  SUBROUTINE shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
619  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  dtsqby2 = dt*dt*.5_dp
643  idtsq = 1.0_dp/dt/dt
644  dtby2 = dt*.5_dp
645  DO iconst = 1, ng4x6
646  IF (g4x6_list(iconst)%restraint%active) cycle
647  index_a = g4x6_list(iconst)%a + first_atom - 1
648  index_b = g4x6_list(iconst)%b + first_atom - 1
649  index_c = g4x6_list(iconst)%c + first_atom - 1
650  index_d = g4x6_list(iconst)%d + first_atom - 1
651  IF (ishake == 1) THEN
652  r0_12(:) = pos(:, index_a) - pos(:, index_b)
653  r0_13(:) = pos(:, index_a) - pos(:, index_c)
654  r0_23(:) = pos(:, index_b) - pos(:, index_c)
655  r0_14(:) = pos(:, index_a) - pos(:, index_d)
656  r0_24(:) = pos(:, index_b) - pos(:, index_d)
657  r0_34(:) = pos(:, index_c) - pos(:, index_d)
658  atomic_kind => particle_set(index_a)%atomic_kind
659  imass1 = 1.0_dp/atomic_kind%mass
660  atomic_kind => particle_set(index_b)%atomic_kind
661  imass2 = 1.0_dp/atomic_kind%mass
662  atomic_kind => particle_set(index_c)%atomic_kind
663  imass3 = 1.0_dp/atomic_kind%mass
664  atomic_kind => particle_set(index_d)%atomic_kind
665  imass4 = 1.0_dp/atomic_kind%mass
666  lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
667  lg4x6(iconst)%rb_old)
668  lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
669  lg4x6(iconst)%rc_old)
670  lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
671  lg4x6(iconst)%rd_old)
672  lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
673  lg4x6(iconst)%rc_old)
674  lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
675  lg4x6(iconst)%rd_old)
676  lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
677  lg4x6(iconst)%rd_old)
678  ! Check for fixed atom constraints
679  CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
680  index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
681  ! rotate fconst:
682  f_roll1 = matmul(r_shake, lg4x6(iconst)%fa)
683  f_roll2 = matmul(r_shake, lg4x6(iconst)%fb)
684  f_roll3 = matmul(r_shake, lg4x6(iconst)%fc)
685  f_roll4 = matmul(r_shake, lg4x6(iconst)%fd)
686  f_roll5 = matmul(r_shake, lg4x6(iconst)%fe)
687  f_roll6 = matmul(r_shake, lg4x6(iconst)%ff)
688 
689  ! construct matrix
690  amat(1, 1) = (imass1 + imass2)*dot_product(r0_12, f_roll1)
691  amat(1, 2) = imass1*dot_product(r0_12, f_roll2)
692  amat(1, 3) = imass1*dot_product(r0_12, f_roll3)
693  amat(1, 4) = -imass2*dot_product(r0_12, f_roll4)
694  amat(1, 5) = -imass2*dot_product(r0_12, f_roll5)
695  amat(1, 6) = 0.0_dp
696 
697  amat(2, 1) = imass1*dot_product(r0_13, f_roll1)
698  amat(2, 2) = (imass1 + imass3)*dot_product(r0_13, f_roll2)
699  amat(2, 3) = imass1*dot_product(r0_13, f_roll3)
700  amat(2, 4) = imass3*dot_product(r0_13, f_roll4)
701  amat(2, 5) = 0.0_dp
702  amat(2, 6) = -imass3*dot_product(r0_13, f_roll6)
703 
704  amat(3, 1) = imass1*dot_product(r0_14, f_roll1)
705  amat(3, 2) = imass1*dot_product(r0_14, f_roll2)
706  amat(3, 3) = (imass1 + imass4)*dot_product(r0_14, f_roll3)
707  amat(3, 4) = 0.0_dp
708  amat(3, 5) = imass4*dot_product(r0_14, f_roll5)
709  amat(3, 6) = imass4*dot_product(r0_14, f_roll6)
710 
711  amat(4, 1) = -imass2*dot_product(r0_23, f_roll1)
712  amat(4, 2) = imass3*dot_product(r0_23, f_roll2)
713  amat(4, 3) = 0.0_dp
714  amat(4, 4) = (imass3 + imass2)*dot_product(r0_23, f_roll4)
715  amat(4, 5) = imass2*dot_product(r0_23, f_roll5)
716  amat(4, 6) = -imass3*dot_product(r0_23, f_roll6)
717 
718  amat(5, 1) = -imass2*dot_product(r0_24, f_roll1)
719  amat(5, 2) = 0.0_dp
720  amat(5, 3) = imass4*dot_product(r0_24, f_roll3)
721  amat(5, 4) = imass2*dot_product(r0_24, f_roll4)
722  amat(5, 5) = (imass4 + imass2)*dot_product(r0_24, f_roll5)
723  amat(5, 6) = imass4*dot_product(r0_24, f_roll6)
724 
725  amat(6, 1) = 0.0_dp
726  amat(6, 2) = -imass3*dot_product(r0_34, f_roll2)
727  amat(6, 3) = imass4*dot_product(r0_34, f_roll3)
728  amat(6, 4) = -imass3*dot_product(r0_34, f_roll4)
729  amat(6, 5) = imass4*dot_product(r0_34, f_roll5)
730  amat(6, 6) = (imass3 + imass4)*dot_product(r0_34, f_roll6)
731  ! Store values
732  lg4x6(iconst)%r0_12 = r0_12
733  lg4x6(iconst)%r0_13 = r0_13
734  lg4x6(iconst)%r0_14 = r0_14
735  lg4x6(iconst)%r0_23 = r0_23
736  lg4x6(iconst)%r0_24 = r0_24
737  lg4x6(iconst)%r0_34 = r0_34
738  lg4x6(iconst)%f_roll1 = f_roll1
739  lg4x6(iconst)%f_roll2 = f_roll2
740  lg4x6(iconst)%f_roll3 = f_roll3
741  lg4x6(iconst)%f_roll4 = f_roll4
742  lg4x6(iconst)%f_roll5 = f_roll5
743  lg4x6(iconst)%f_roll6 = f_roll6
744  lg4x6(iconst)%amat = amat
745  lg4x6(iconst)%imass1 = imass1
746  lg4x6(iconst)%imass2 = imass2
747  lg4x6(iconst)%imass3 = imass3
748  lg4x6(iconst)%imass4 = imass4
749  lg4x6(iconst)%lambda_old = 0.0_dp
750  lg4x6(iconst)%del_lambda = 0.0_dp
751  ELSE
752  ! Retrieve values
753  r0_12 = lg4x6(iconst)%r0_12
754  r0_13 = lg4x6(iconst)%r0_13
755  r0_14 = lg4x6(iconst)%r0_14
756  r0_23 = lg4x6(iconst)%r0_23
757  r0_24 = lg4x6(iconst)%r0_24
758  r0_34 = lg4x6(iconst)%r0_34
759  f_roll1 = lg4x6(iconst)%f_roll1
760  f_roll2 = lg4x6(iconst)%f_roll2
761  f_roll3 = lg4x6(iconst)%f_roll3
762  f_roll4 = lg4x6(iconst)%f_roll4
763  f_roll5 = lg4x6(iconst)%f_roll5
764  f_roll6 = lg4x6(iconst)%f_roll6
765  amat = lg4x6(iconst)%amat
766  imass1 = lg4x6(iconst)%imass1
767  imass2 = lg4x6(iconst)%imass2
768  imass3 = lg4x6(iconst)%imass3
769  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  lg4x6(iconst)%lambda(5)*imass2*f_roll5
778  bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
779  - 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  lg4x6(iconst)%lambda(6)*imass3*f_roll6
785  bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
786  - 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  lg4x6(iconst)%lambda(6)*imass4*f_roll6
792  bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
793  - 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  lg4x6(iconst)%lambda(6)*imass3*f_roll6
799  bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
800  - 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  lg4x6(iconst)%lambda(6)*imass4*f_roll6
806  bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
807  - 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  lg4x6(iconst)%lambda(5)*imass4*f_roll5
813  bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
814  - dtsqby2*dtsqby2*dot_product(vec, vec) - dot_product(r0_34, r0_34)
815  bvec = bvec*idtsq
816 
817  ! get lambda
818  atemp = amat
819  CALL solve_system(atemp, 6, bvec)
820  lg4x6(iconst)%lambda(:) = bvec(:, 1)
821  lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
822  lg4x6(iconst)%lambda_old(:)
823  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  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  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  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  lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
837  r1(:) = pos(:, index_a) + imass1*dtsqby2*matmul(r_shake, fc1)
838  r2(:) = pos(:, index_b) + imass2*dtsqby2*matmul(r_shake, fc2)
839  r3(:) = pos(:, index_c) + imass3*dtsqby2*matmul(r_shake, fc3)
840  r4(:) = pos(:, index_d) + imass4*dtsqby2*matmul(r_shake, fc4)
841  v1(:) = vel(:, index_a) + imass1*dtby2*matmul(r_shake, fc1)
842  v2(:) = vel(:, index_b) + imass2*dtby2*matmul(r_shake, fc2)
843  v3(:) = vel(:, index_c) + imass3*dtby2*matmul(r_shake, fc3)
844  v4(:) = vel(:, index_d) + imass4*dtby2*matmul(r_shake, fc4)
845 
846  r12 = r1 - r2
847  r13 = r1 - r3
848  r23 = r2 - r3
849  r14 = r1 - r4
850  r24 = r2 - r4
851  r34 = r3 - r4
852  ! compute the tolerance:
853  sigma = dot_product(r12, r12) - g4x6_list(iconst)%dab* &
854  g4x6_list(iconst)%dab
855  max_sigma = max(max_sigma, abs(sigma))
856  sigma = dot_product(r13, r13) - g4x6_list(iconst)%dac* &
857  g4x6_list(iconst)%dac
858  max_sigma = max(max_sigma, abs(sigma))
859  sigma = dot_product(r14, r14) - g4x6_list(iconst)%dad* &
860  g4x6_list(iconst)%dad
861  max_sigma = max(max_sigma, abs(sigma))
862  sigma = dot_product(r23, r23) - g4x6_list(iconst)%dbc* &
863  g4x6_list(iconst)%dbc
864  max_sigma = max(max_sigma, abs(sigma))
865  sigma = dot_product(r24, r24) - g4x6_list(iconst)%dbd* &
866  g4x6_list(iconst)%dbd
867  max_sigma = max(max_sigma, abs(sigma))
868  sigma = dot_product(r34, r34) - g4x6_list(iconst)%dcd* &
869  g4x6_list(iconst)%dcd
870  max_sigma = max(max_sigma, abs(sigma))
871 
872  ! update positions with full multiplier
873  pos(:, index_a) = r1(:)
874  pos(:, index_b) = r2(:)
875  pos(:, index_c) = r3(:)
876  pos(:, index_d) = r4(:)
877 
878  ! update velocites with full multiplier
879  vel(:, index_a) = v1(:)
880  vel(:, index_b) = v2(:)
881  vel(:, index_c) = v3(:)
882  vel(:, index_d) = v4(:)
883  END DO
884  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  SUBROUTINE rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
899  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  idt = 1.0_dp/dt
920  dtby2 = dt*.5_dp
921  DO iconst = 1, SIZE(g4x6_list)
922  IF (g4x6_list(iconst)%restraint%active) cycle
923  index_a = g4x6_list(iconst)%a + first_atom - 1
924  index_b = g4x6_list(iconst)%b + first_atom - 1
925  index_c = g4x6_list(iconst)%c + first_atom - 1
926  index_d = g4x6_list(iconst)%d + first_atom - 1
927  v12(:) = vel(:, index_a) - vel(:, index_b)
928  v13(:) = vel(:, index_a) - vel(:, index_c)
929  v14(:) = vel(:, index_a) - vel(:, index_d)
930  v23(:) = vel(:, index_b) - vel(:, index_c)
931  v24(:) = vel(:, index_b) - vel(:, index_d)
932  v34(:) = vel(:, index_c) - vel(:, index_d)
933 
934  r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
935  r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
936  r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
937  r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
938  r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
939  r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
940  atomic_kind => particle_set(index_a)%atomic_kind
941  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
942  imass1 = 1.0_dp/mass
943  atomic_kind => particle_set(index_b)%atomic_kind
944  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
945  imass2 = 1.0_dp/mass
946  atomic_kind => particle_set(index_c)%atomic_kind
947  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
948  imass3 = 1.0_dp/mass
949  atomic_kind => particle_set(index_d)%atomic_kind
950  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
951  imass4 = 1.0_dp/mass
952  lg4x6(iconst)%fa = -2.0_dp*r12
953  lg4x6(iconst)%fb = -2.0_dp*r13
954  lg4x6(iconst)%fc = -2.0_dp*r14
955  lg4x6(iconst)%fd = -2.0_dp*r23
956  lg4x6(iconst)%fe = -2.0_dp*r24
957  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  index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
961  ! construct matrix
962  amat(1, 1) = (imass1 + imass2)*dot_product(r12, lg4x6(iconst)%fa)
963  amat(1, 2) = imass1*dot_product(r12, lg4x6(iconst)%fb)
964  amat(1, 3) = imass1*dot_product(r12, lg4x6(iconst)%fc)
965  amat(1, 4) = -imass2*dot_product(r12, lg4x6(iconst)%fd)
966  amat(1, 5) = -imass2*dot_product(r12, lg4x6(iconst)%fe)
967  amat(1, 6) = 0.0_dp
968 
969  amat(2, 1) = imass1*dot_product(r13, lg4x6(iconst)%fa)
970  amat(2, 2) = (imass1 + imass3)*dot_product(r13, lg4x6(iconst)%fb)
971  amat(2, 3) = imass1*dot_product(r13, lg4x6(iconst)%fc)
972  amat(2, 4) = imass3*dot_product(r13, lg4x6(iconst)%fd)
973  amat(2, 5) = 0.0_dp
974  amat(2, 6) = -imass3*dot_product(r13, lg4x6(iconst)%ff)
975 
976  amat(3, 1) = imass1*dot_product(r14, lg4x6(iconst)%fa)
977  amat(3, 2) = imass1*dot_product(r14, lg4x6(iconst)%fb)
978  amat(3, 3) = (imass1 + imass4)*dot_product(r14, lg4x6(iconst)%fc)
979  amat(3, 4) = 0.0_dp
980  amat(3, 5) = imass4*dot_product(r14, lg4x6(iconst)%fe)
981  amat(3, 6) = imass4*dot_product(r14, lg4x6(iconst)%ff)
982 
983  amat(4, 1) = -imass2*dot_product(r23, lg4x6(iconst)%fa)
984  amat(4, 2) = imass3*dot_product(r23, lg4x6(iconst)%fb)
985  amat(4, 3) = 0.0_dp
986  amat(4, 4) = (imass3 + imass2)*dot_product(r23, lg4x6(iconst)%fd)
987  amat(4, 5) = imass2*dot_product(r23, lg4x6(iconst)%fe)
988  amat(4, 6) = -imass3*dot_product(r23, lg4x6(iconst)%ff)
989 
990  amat(5, 1) = -imass2*dot_product(r24, lg4x6(iconst)%fa)
991  amat(5, 2) = 0.0_dp
992  amat(5, 3) = imass4*dot_product(r24, lg4x6(iconst)%fc)
993  amat(5, 4) = imass2*dot_product(r24, lg4x6(iconst)%fd)
994  amat(5, 5) = (imass4 + imass2)*dot_product(r24, lg4x6(iconst)%fe)
995  amat(5, 6) = imass4*dot_product(r24, lg4x6(iconst)%ff)
996 
997  amat(6, 1) = 0.0_dp
998  amat(6, 2) = -imass3*dot_product(r34, lg4x6(iconst)%fb)
999  amat(6, 3) = imass4*dot_product(r34, lg4x6(iconst)%fc)
1000  amat(6, 4) = -imass3*dot_product(r34, lg4x6(iconst)%fd)
1001  amat(6, 5) = imass4*dot_product(r34, lg4x6(iconst)%fe)
1002  amat(6, 6) = (imass3 + imass4)*dot_product(r34, lg4x6(iconst)%ff)
1003 
1004  ! construct solution vector
1005  bvec(1, 1) = dot_product(r12, v12)
1006  bvec(2, 1) = dot_product(r13, v13)
1007  bvec(3, 1) = dot_product(r14, v14)
1008  bvec(4, 1) = dot_product(r23, v23)
1009  bvec(5, 1) = dot_product(r24, v24)
1010  bvec(6, 1) = dot_product(r34, v34)
1011  bvec = -bvec*2.0_dp*idt
1012 
1013  ! get lambda
1014  CALL solve_system(amat, 6, bvec)
1015  lg4x6(iconst)%lambda(:) = bvec(:, 1)
1016 
1017  fc1 = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + &
1018  lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb + &
1019  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  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  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  lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff
1029  vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
1030  vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
1031  vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
1032  vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
1033  END DO
1034  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  SUBROUTINE rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
1051  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  idt = 1.0_dp/dt
1077  dtby2 = dt*.5_dp
1078  DO iconst = 1, SIZE(g4x6_list)
1079  IF (g4x6_list(iconst)%restraint%active) cycle
1080  index_a = g4x6_list(iconst)%a + first_atom - 1
1081  index_b = g4x6_list(iconst)%b + first_atom - 1
1082  index_c = g4x6_list(iconst)%c + first_atom - 1
1083  index_d = g4x6_list(iconst)%d + first_atom - 1
1084  v12(:) = vel(:, index_a) - vel(:, index_b)
1085  v13(:) = vel(:, index_a) - vel(:, index_c)
1086  v14(:) = vel(:, index_a) - vel(:, index_d)
1087  v23(:) = vel(:, index_b) - vel(:, index_c)
1088  v24(:) = vel(:, index_b) - vel(:, index_d)
1089  v34(:) = vel(:, index_c) - vel(:, index_d)
1090 
1091  r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
1092  r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
1093  r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
1094  r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
1095  r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
1096  r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
1097  atomic_kind => particle_set(index_a)%atomic_kind
1098  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1099  imass1 = 1.0_dp/mass
1100  atomic_kind => particle_set(index_b)%atomic_kind
1101  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1102  imass2 = 1.0_dp/mass
1103  atomic_kind => particle_set(index_c)%atomic_kind
1104  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1105  imass3 = 1.0_dp/mass
1106  atomic_kind => particle_set(index_d)%atomic_kind
1107  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
1108  imass4 = 1.0_dp/mass
1109  lg4x6(iconst)%fa = -2.0_dp*r12
1110  lg4x6(iconst)%fb = -2.0_dp*r13
1111  lg4x6(iconst)%fc = -2.0_dp*r14
1112  lg4x6(iconst)%fd = -2.0_dp*r23
1113  lg4x6(iconst)%fe = -2.0_dp*r24
1114  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  index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
1118  ! roll the fc
1119  f_roll1 = matmul(r_rattle, lg4x6(iconst)%fa)
1120  f_roll2 = matmul(r_rattle, lg4x6(iconst)%fb)
1121  f_roll3 = matmul(r_rattle, lg4x6(iconst)%fc)
1122  f_roll4 = matmul(r_rattle, lg4x6(iconst)%fd)
1123  f_roll5 = matmul(r_rattle, lg4x6(iconst)%fe)
1124  f_roll6 = matmul(r_rattle, lg4x6(iconst)%ff)
1125  ! construct matrix
1126  amat(1, 1) = (imass1 + imass2)*dot_product(r12, f_roll1)
1127  amat(1, 2) = imass1*dot_product(r12, f_roll2)
1128  amat(1, 3) = imass1*dot_product(r12, f_roll3)
1129  amat(1, 4) = -imass2*dot_product(r12, f_roll4)
1130  amat(1, 5) = -imass2*dot_product(r12, f_roll5)
1131  amat(1, 6) = 0.0_dp
1132 
1133  amat(2, 1) = imass1*dot_product(r13, f_roll1)
1134  amat(2, 2) = (imass1 + imass3)*dot_product(r13, f_roll2)
1135  amat(2, 3) = imass1*dot_product(r13, f_roll3)
1136  amat(2, 4) = imass3*dot_product(r13, f_roll4)
1137  amat(2, 5) = 0.0_dp
1138  amat(2, 6) = -imass3*dot_product(r13, f_roll6)
1139 
1140  amat(3, 1) = imass1*dot_product(r14, f_roll1)
1141  amat(3, 2) = imass1*dot_product(r14, f_roll2)
1142  amat(3, 3) = (imass1 + imass4)*dot_product(r14, f_roll3)
1143  amat(3, 4) = 0.0_dp
1144  amat(3, 5) = imass4*dot_product(r14, f_roll5)
1145  amat(3, 6) = imass4*dot_product(r14, f_roll6)
1146 
1147  amat(4, 1) = -imass2*dot_product(r23, f_roll1)
1148  amat(4, 2) = imass3*dot_product(r23, f_roll2)
1149  amat(4, 3) = 0.0_dp
1150  amat(4, 4) = (imass3 + imass2)*dot_product(r23, f_roll4)
1151  amat(4, 5) = imass2*dot_product(r23, f_roll5)
1152  amat(4, 6) = -imass3*dot_product(r23, f_roll6)
1153 
1154  amat(5, 1) = -imass2*dot_product(r24, f_roll1)
1155  amat(5, 2) = 0.0_dp
1156  amat(5, 3) = imass4*dot_product(r24, f_roll3)
1157  amat(5, 4) = imass2*dot_product(r24, f_roll4)
1158  amat(5, 5) = (imass4 + imass2)*dot_product(r24, f_roll5)
1159  amat(5, 6) = imass4*dot_product(r24, f_roll6)
1160 
1161  amat(6, 1) = 0.0_dp
1162  amat(6, 2) = -imass3*dot_product(r34, f_roll2)
1163  amat(6, 3) = imass4*dot_product(r34, f_roll3)
1164  amat(6, 4) = -imass3*dot_product(r34, f_roll4)
1165  amat(6, 5) = imass4*dot_product(r34, f_roll5)
1166  amat(6, 6) = (imass3 + imass4)*dot_product(r34, f_roll6)
1167 
1168  ! construct solution vector
1169  bvec(1, 1) = dot_product(r12, v12 + matmul(veps, r12))
1170  bvec(2, 1) = dot_product(r13, v13 + matmul(veps, r13))
1171  bvec(3, 1) = dot_product(r14, v14 + matmul(veps, r14))
1172  bvec(4, 1) = dot_product(r23, v23 + matmul(veps, r23))
1173  bvec(5, 1) = dot_product(r24, v24 + matmul(veps, r24))
1174  bvec(6, 1) = dot_product(r34, v34 + matmul(veps, r34))
1175  bvec = -bvec*2.0_dp*idt
1176 
1177  ! get lambda
1178  CALL solve_system(amat, 6, bvec)
1179  lambda(:) = bvec(:, 1)
1180  lg4x6(iconst)%lambda(:) = lambda
1181 
1182  fc1 = lambda(1)*f_roll1 + &
1183  lambda(2)*f_roll2 + &
1184  lambda(3)*f_roll3
1185  fc2 = -lambda(1)*f_roll1 + &
1186  lambda(4)*f_roll4 + &
1187  lambda(5)*f_roll5
1188  fc3 = -lambda(2)*f_roll2 - &
1189  lambda(4)*f_roll4 + &
1190  lambda(6)*f_roll6
1191  fc4 = -lambda(3)*f_roll3 - &
1192  lambda(5)*f_roll5 - &
1193  lambda(6)*f_roll6
1194  vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
1195  vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
1196  vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
1197  vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
1198  END DO
1199  END SUBROUTINE rattle_roll_4x6_low
1200 
1201 END MODULE constraint_4x6
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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.
Definition: atomic_kind_types.F:138
constraint_4x6
Definition: constraint_4x6.F:12
constraint_4x6::rattle_roll_4x6_int
subroutine, public rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, veps)
Intramolecular rattle_4x6_roll.
Definition: constraint_4x6.F:168
constraint_4x6::rattle_4x6_int
subroutine, public rattle_4x6_int(molecule, particle_set, vel, dt)
Intramolecular rattle_4x6.
Definition: constraint_4x6.F:134
constraint_4x6::shake_roll_4x6_ext
subroutine, public shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
Intramolecular shake_4x6_roll.
Definition: constraint_4x6.F:245
constraint_4x6::rattle_4x6_ext
subroutine, public rattle_4x6_ext(gci, particle_set, vel, dt)
Intramolecular rattle_4x6.
Definition: constraint_4x6.F:280
constraint_4x6::shake_4x6_ext
subroutine, public shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake, max_sigma)
Intramolecular shake_4x6.
Definition: constraint_4x6.F:206
constraint_4x6::shake_roll_4x6_int
subroutine, public shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
Intramolecular shake_4x6_roll.
Definition: constraint_4x6.F:99
constraint_4x6::rattle_roll_4x6_ext
subroutine, public rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, veps)
Intramolecular rattle_4x6_roll.
Definition: constraint_4x6.F:314
constraint_4x6::shake_4x6_int
subroutine, public shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake, max_sigma)
Intramolecular shake_4x6.
Definition: constraint_4x6.F:60
constraint_fxd
Definition: constraint_fxd.F:12
constraint_fxd::check_fixed_atom_cns_g4x6
subroutine, public check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, index_a, index_b, index_c, index_d, fixd_list, lg4x6)
...
Definition: constraint_fxd.F:327
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
linear_systems
Provides interfaces to LAPACK routines for factorisation and linear system solving.
Definition: linear_systems.F:18
linear_systems::solve_system
subroutine, public solve_system(matrix, mysize, eigenvectors)
...
Definition: linear_systems.F:41
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition: molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition: molecule_kind_types.F:510
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition: molecule_types.F:277
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19