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