(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
12MODULE 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
44CONTAINS
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
1201END 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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
molecule_kind_types::fixd_constraint_type
Definition molecule_kind_types.F:149
molecule_kind_types::g4x6_constraint_type
Definition molecule_kind_types.F:135
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_types::global_constraint_type
Definition molecule_types.F:93
molecule_types::local_g4x6_constraint_type
Definition molecule_types.F:63
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35