LCOV - code coverage report
Current view: top level - src - constraint_3x3.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 95.5 % 352 336
Test Date: 2025-07-25 12:55:17 Functions: 83.3 % 12 10

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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       281363 :    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       281363 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
      70       281363 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
      71       281363 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
      72              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      73              : 
      74       281363 :       NULLIFY (fixd_list)
      75       281363 :       molecule_kind => molecule%molecule_kind
      76              :       CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
      77       281363 :                              g3x3_list=g3x3_list, fixd_list=fixd_list)
      78       281363 :       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       281363 :                          particle_set, pos, vel, dt, ishake, max_sigma)
      82              : 
      83       281363 :    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       128355 :    SUBROUTINE shake_roll_3x3_int(molecule, particle_set, pos, vel, r_shake, &
     100       128355 :                                  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       128355 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     112       128355 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     113       128355 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     114              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     115              : 
     116       128355 :       NULLIFY (fixd_list)
     117       128355 :       molecule_kind => molecule%molecule_kind
     118              :       CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
     119       128355 :                              g3x3_list=g3x3_list, fixd_list=fixd_list)
     120       128355 :       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       128355 :                               particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
     124              : 
     125       128355 :    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        97999 :    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        97999 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     149        97999 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     150        97999 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     151              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     152              : 
     153        97999 :       NULLIFY (fixd_list)
     154        97999 :       molecule_kind => molecule%molecule_kind
     155              :       CALL get_molecule_kind(molecule_kind, &
     156        97999 :                              g3x3_list=g3x3_list, fixd_list=fixd_list)
     157        97999 :       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        97999 :                                particle_set, vel, r_rattle, dt, veps)
     161              : 
     162        97999 :    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       143796 :    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       143796 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     182       143796 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     183       143796 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     184              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     185              : 
     186       143796 :       NULLIFY (fixd_list)
     187       143796 :       molecule_kind => molecule%molecule_kind
     188              :       CALL get_molecule_kind(molecule_kind, &
     189       143796 :                              g3x3_list=g3x3_list, fixd_list=fixd_list)
     190       143796 :       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       143796 :                           particle_set, vel, dt)
     194              : 
     195       143796 :    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           76 :    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           76 :       first_atom = 1
     225           76 :       ng3x3 = gci%ng3x3
     226           76 :       g3x3_list => gci%g3x3_list
     227           76 :       fixd_list => gci%fixd_list
     228           76 :       lg3x3 => gci%lg3x3
     229              :       ! Real Shake
     230              :       CALL shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
     231           76 :                          particle_set, pos, vel, dt, ishake, max_sigma)
     232              : 
     233           76 :    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            0 :    SUBROUTINE shake_roll_3x3_ext(gci, particle_set, pos, vel, r_shake, &
     250            0 :                                  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            0 :       first_atom = 1
     266            0 :       ng3x3 = gci%ng3x3
     267            0 :       g3x3_list => gci%g3x3_list
     268            0 :       fixd_list => gci%fixd_list
     269            0 :       lg3x3 => gci%lg3x3
     270              :       ! Real Shake
     271              :       CALL shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
     272            0 :                               particle_set, pos, vel, r_shake, v_shake, dt, ishake, max_sigma)
     273              : 
     274            0 :    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            0 :    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            0 :       first_atom = 1
     302            0 :       g3x3_list => gci%g3x3_list
     303            0 :       fixd_list => gci%fixd_list
     304            0 :       lg3x3 => gci%lg3x3
     305              :       ! Real Rattle
     306              :       CALL rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
     307            0 :                                particle_set, vel, r_rattle, dt, veps)
     308              : 
     309            0 :    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           40 :    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           40 :       first_atom = 1
     333           40 :       g3x3_list => gci%g3x3_list
     334           40 :       fixd_list => gci%fixd_list
     335           40 :       lg3x3 => gci%lg3x3
     336              :       ! Real Rattle
     337              :       CALL rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
     338           40 :                           particle_set, vel, dt)
     339              : 
     340           40 :    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       281439 :    SUBROUTINE shake_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
     359       281439 :                             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       281439 :       dtsqby2 = dt*dt*.5_dp
     381       281439 :       idtsq = 1.0_dp/dt/dt
     382       281439 :       dtby2 = dt*.5_dp
     383       562878 :       DO iconst = 1, ng3x3
     384       281439 :          IF (g3x3_list(iconst)%restraint%active) CYCLE
     385       281419 :          index_a = g3x3_list(iconst)%a + first_atom - 1
     386       281419 :          index_b = g3x3_list(iconst)%b + first_atom - 1
     387       281419 :          index_c = g3x3_list(iconst)%c + first_atom - 1
     388       281419 :          IF (ishake == 1) THEN
     389       575264 :             r0_12(:) = pos(:, index_a) - pos(:, index_b)
     390       575264 :             r0_13(:) = pos(:, index_a) - pos(:, index_c)
     391       575264 :             r0_23(:) = pos(:, index_b) - pos(:, index_c)
     392       143816 :             atomic_kind => particle_set(index_a)%atomic_kind
     393       143816 :             imass1 = 1.0_dp/atomic_kind%mass
     394       143816 :             atomic_kind => particle_set(index_b)%atomic_kind
     395       143816 :             imass2 = 1.0_dp/atomic_kind%mass
     396       143816 :             atomic_kind => particle_set(index_c)%atomic_kind
     397       143816 :             imass3 = 1.0_dp/atomic_kind%mass
     398              :             lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
     399       575264 :                                         lg3x3(iconst)%rb_old)
     400              :             lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
     401       575264 :                                         lg3x3(iconst)%rc_old)
     402              :             lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
     403       575264 :                                         lg3x3(iconst)%rc_old)
     404              :             ! Check for fixed atom constraints
     405              :             CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
     406       143816 :                                            index_a, index_b, index_c, fixd_list, lg3x3(iconst))
     407              :             ! construct matrix
     408       575264 :             amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, lg3x3(iconst)%fa)
     409       575264 :             amat(1, 2) = imass1*DOT_PRODUCT(r0_12, lg3x3(iconst)%fb)
     410       575264 :             amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, lg3x3(iconst)%fc)
     411       575264 :             amat(2, 1) = imass1*DOT_PRODUCT(r0_13, lg3x3(iconst)%fa)
     412       575264 :             amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, lg3x3(iconst)%fb)
     413       575264 :             amat(2, 3) = imass3*DOT_PRODUCT(r0_13, lg3x3(iconst)%fc)
     414       575264 :             amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, lg3x3(iconst)%fa)
     415       575264 :             amat(3, 2) = imass3*DOT_PRODUCT(r0_23, lg3x3(iconst)%fb)
     416       575264 :             amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, lg3x3(iconst)%fc)
     417              :             ! Store values
     418       575264 :             lg3x3(iconst)%r0_12 = r0_12
     419       575264 :             lg3x3(iconst)%r0_13 = r0_13
     420       575264 :             lg3x3(iconst)%r0_23 = r0_23
     421      1869608 :             lg3x3(iconst)%amat = amat
     422       575264 :             lg3x3(iconst)%lambda_old = 0.0_dp
     423       575264 :             lg3x3(iconst)%del_lambda = 0.0_dp
     424              :          ELSE
     425              :             ! Retrieve values
     426       550412 :             r0_12 = lg3x3(iconst)%r0_12
     427       550412 :             r0_13 = lg3x3(iconst)%r0_13
     428       550412 :             r0_23 = lg3x3(iconst)%r0_23
     429      1788839 :             amat = lg3x3(iconst)%amat
     430       137603 :             imass1 = lg3x3(iconst)%imass1
     431       137603 :             imass2 = lg3x3(iconst)%imass2
     432       137603 :             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      1125676 :                lg3x3(iconst)%lambda(3)*imass2*lg3x3(iconst)%fc
     438              :          bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
     439      2251352 :                       - 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      1125676 :                lg3x3(iconst)%lambda(3)*imass3*lg3x3(iconst)%fc
     443              :          bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
     444      2251352 :                       - 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      1125676 :                lg3x3(iconst)%lambda(3)*(imass2 + imass3)*lg3x3(iconst)%fc
     448              :          bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
     449      2251352 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
     450      1407095 :          bvec = bvec*idtsq
     451              :          ! get lambda
     452       281419 :          atemp = amat
     453       281419 :          CALL solve_system(atemp, 3, bvec)
     454      1125676 :          lg3x3(iconst)%lambda(:) = bvec(:, 1)
     455              :          lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
     456      1125676 :                                        lg3x3(iconst)%lambda_old(:)
     457      1125676 :          lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
     458              :          fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
     459      1125676 :                lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
     460              :          fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
     461      1125676 :                lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
     462              :          fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
     463      1125676 :                lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
     464      1125676 :          r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
     465      1125676 :          r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
     466      1125676 :          r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
     467      1125676 :          v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
     468      1125676 :          v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
     469      1125676 :          v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
     470      1125676 :          r12 = r1 - r2
     471      1125676 :          r13 = r1 - r3
     472      1125676 :          r23 = r2 - r3
     473              : 
     474              :          ! compute the tolerance:
     475              :          sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* &
     476      1125676 :                  g3x3_list(iconst)%dab
     477       281419 :          max_sigma = MAX(max_sigma, ABS(sigma))
     478              :          sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* &
     479      1125676 :                  g3x3_list(iconst)%dac
     480       281419 :          max_sigma = MAX(max_sigma, ABS(sigma))
     481              :          sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* &
     482      1125676 :                  g3x3_list(iconst)%dbc
     483       281419 :          max_sigma = MAX(max_sigma, ABS(sigma))
     484              :          ! update positions with full multiplier
     485      1125676 :          pos(:, index_a) = r1(:)
     486      1125676 :          pos(:, index_b) = r2(:)
     487      1125676 :          pos(:, index_c) = r3(:)
     488              : 
     489              :          ! update velocites with full multiplier
     490      1125676 :          vel(:, index_a) = v1(:)
     491      1125676 :          vel(:, index_b) = v2(:)
     492      1407135 :          vel(:, index_c) = v3(:)
     493              :       END DO
     494              : 
     495       281439 :    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       128355 :    SUBROUTINE shake_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, ng3x3, &
     516       128355 :                                  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       128355 :       dtsqby2 = dt*dt*.5_dp
     540       128355 :       idtsq = 1.0_dp/dt/dt
     541       128355 :       dtby2 = dt*.5_dp
     542       256710 :       DO iconst = 1, ng3x3
     543       128355 :          IF (g3x3_list(iconst)%restraint%active) CYCLE
     544       128355 :          index_a = g3x3_list(iconst)%a + first_atom - 1
     545       128355 :          index_b = g3x3_list(iconst)%b + first_atom - 1
     546       128355 :          index_c = g3x3_list(iconst)%c + first_atom - 1
     547       128355 :          IF (ishake == 1) THEN
     548       450712 :             r0_12(:) = pos(:, index_a) - pos(:, index_b)
     549       450712 :             r0_13(:) = pos(:, index_a) - pos(:, index_c)
     550       450712 :             r0_23(:) = pos(:, index_b) - pos(:, index_c)
     551       112678 :             atomic_kind => particle_set(index_a)%atomic_kind
     552       112678 :             imass1 = 1.0_dp/atomic_kind%mass
     553       112678 :             atomic_kind => particle_set(index_b)%atomic_kind
     554       112678 :             imass2 = 1.0_dp/atomic_kind%mass
     555       112678 :             atomic_kind => particle_set(index_c)%atomic_kind
     556       112678 :             imass3 = 1.0_dp/atomic_kind%mass
     557              :             lg3x3(iconst)%fa = -2.0_dp*(lg3x3(iconst)%ra_old - &
     558       450712 :                                         lg3x3(iconst)%rb_old)
     559              :             lg3x3(iconst)%fb = -2.0_dp*(lg3x3(iconst)%ra_old - &
     560       450712 :                                         lg3x3(iconst)%rc_old)
     561              :             lg3x3(iconst)%fc = -2.0_dp*(lg3x3(iconst)%rb_old - &
     562       450712 :                                         lg3x3(iconst)%rc_old)
     563              :             ! Check for fixed atom constraints
     564              :             CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
     565       112678 :                                            index_a, index_b, index_c, fixd_list, lg3x3(iconst))
     566              :             ! rotate fconst:
     567      1464814 :             f_roll1 = MATMUL(r_shake, lg3x3(iconst)%fa)
     568      1464814 :             f_roll2 = MATMUL(r_shake, lg3x3(iconst)%fb)
     569      1464814 :             f_roll3 = MATMUL(r_shake, lg3x3(iconst)%fc)
     570              :             ! construct matrix
     571       450712 :             amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, f_roll1)
     572       450712 :             amat(1, 2) = imass1*DOT_PRODUCT(r0_12, f_roll2)
     573       450712 :             amat(1, 3) = -imass2*DOT_PRODUCT(r0_12, f_roll3)
     574       450712 :             amat(2, 1) = imass1*DOT_PRODUCT(r0_13, f_roll1)
     575       450712 :             amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, f_roll2)
     576       450712 :             amat(2, 3) = imass3*DOT_PRODUCT(r0_13, f_roll3)
     577       450712 :             amat(3, 1) = -imass2*DOT_PRODUCT(r0_23, f_roll1)
     578       450712 :             amat(3, 2) = imass3*DOT_PRODUCT(r0_23, f_roll2)
     579       450712 :             amat(3, 3) = (imass3 + imass2)*DOT_PRODUCT(r0_23, f_roll3)
     580              :             ! Store values
     581       450712 :             lg3x3(iconst)%r0_12 = r0_12
     582       450712 :             lg3x3(iconst)%r0_13 = r0_13
     583       450712 :             lg3x3(iconst)%r0_23 = r0_23
     584       450712 :             lg3x3(iconst)%f_roll1 = f_roll1
     585       450712 :             lg3x3(iconst)%f_roll2 = f_roll2
     586       450712 :             lg3x3(iconst)%f_roll3 = f_roll3
     587      1464814 :             lg3x3(iconst)%amat = amat
     588       450712 :             lg3x3(iconst)%lambda_old = 0.0_dp
     589       450712 :             lg3x3(iconst)%del_lambda = 0.0_dp
     590              :          ELSE
     591              :             ! Retrieve values
     592        62708 :             r0_12 = lg3x3(iconst)%r0_12
     593        62708 :             r0_13 = lg3x3(iconst)%r0_13
     594        62708 :             r0_23 = lg3x3(iconst)%r0_23
     595        62708 :             f_roll1 = lg3x3(iconst)%f_roll1
     596        62708 :             f_roll2 = lg3x3(iconst)%f_roll2
     597        62708 :             f_roll3 = lg3x3(iconst)%f_roll3
     598       203801 :             amat = lg3x3(iconst)%amat
     599        15677 :             imass1 = lg3x3(iconst)%imass1
     600        15677 :             imass2 = lg3x3(iconst)%imass2
     601        15677 :             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       513420 :                lg3x3(iconst)%lambda(3)*imass2*f_roll3
     607              :          bvec(1, 1) = g3x3_list(iconst)%dab*g3x3_list(iconst)%dab &
     608      1026840 :                       - 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       513420 :                lg3x3(iconst)%lambda(3)*imass3*f_roll3
     612              :          bvec(2, 1) = g3x3_list(iconst)%dac*g3x3_list(iconst)%dac &
     613      1026840 :                       - 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       513420 :                lg3x3(iconst)%lambda(3)*(imass2 + imass3)*f_roll3
     617              :          bvec(3, 1) = g3x3_list(iconst)%dbc*g3x3_list(iconst)%dbc &
     618      1026840 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
     619       641775 :          bvec = bvec*idtsq
     620              : 
     621              :          ! get lambda
     622       128355 :          atemp = amat
     623       128355 :          CALL solve_system(atemp, 3, bvec)
     624       513420 :          lg3x3(iconst)%lambda(:) = bvec(:, 1)
     625              :          lg3x3(iconst)%del_lambda(:) = lg3x3(iconst)%lambda(:) - &
     626       513420 :                                        lg3x3(iconst)%lambda_old(:)
     627       513420 :          lg3x3(iconst)%lambda_old(:) = lg3x3(iconst)%lambda(:)
     628              : 
     629              :          fc1 = lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
     630       513420 :                lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb
     631              :          fc2 = -lg3x3(iconst)%del_lambda(1)*lg3x3(iconst)%fa + &
     632       513420 :                lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
     633              :          fc3 = -lg3x3(iconst)%del_lambda(2)*lg3x3(iconst)%fb - &
     634       513420 :                lg3x3(iconst)%del_lambda(3)*lg3x3(iconst)%fc
     635      2053680 :          r1(:) = pos(:, index_a) + imass1*dtsqby2*MATMUL(r_shake, fc1)
     636      2053680 :          r2(:) = pos(:, index_b) + imass2*dtsqby2*MATMUL(r_shake, fc2)
     637      2053680 :          r3(:) = pos(:, index_c) + imass3*dtsqby2*MATMUL(r_shake, fc3)
     638      2053680 :          v1(:) = vel(:, index_a) + imass1*dtby2*MATMUL(v_shake, fc1)
     639      2053680 :          v2(:) = vel(:, index_b) + imass2*dtby2*MATMUL(v_shake, fc2)
     640      2053680 :          v3(:) = vel(:, index_c) + imass3*dtby2*MATMUL(v_shake, fc3)
     641       513420 :          r12 = r1 - r2
     642       513420 :          r13 = r1 - r3
     643       513420 :          r23 = r2 - r3
     644              : 
     645              :          ! compute the tolerance:
     646              :          sigma = DOT_PRODUCT(r12, r12) - g3x3_list(iconst)%dab* &
     647       513420 :                  g3x3_list(iconst)%dab
     648       128355 :          max_sigma = MAX(max_sigma, ABS(sigma))
     649              :          sigma = DOT_PRODUCT(r13, r13) - g3x3_list(iconst)%dac* &
     650       513420 :                  g3x3_list(iconst)%dac
     651       128355 :          max_sigma = MAX(max_sigma, ABS(sigma))
     652              :          sigma = DOT_PRODUCT(r23, r23) - g3x3_list(iconst)%dbc* &
     653       513420 :                  g3x3_list(iconst)%dbc
     654       128355 :          max_sigma = MAX(max_sigma, ABS(sigma))
     655              : 
     656              :          ! update positions with full multiplier
     657       513420 :          pos(:, index_a) = r1(:)
     658       513420 :          pos(:, index_b) = r2(:)
     659       513420 :          pos(:, index_c) = r3(:)
     660              : 
     661              :          ! update velocites with full multiplier
     662       513420 :          vel(:, index_a) = v1(:)
     663       513420 :          vel(:, index_b) = v2(:)
     664       641775 :          vel(:, index_c) = v3(:)
     665              :       END DO
     666       128355 :    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        97999 :    SUBROUTINE rattle_roll_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
     683        97999 :                                   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        97999 :       idt = 1.0_dp/dt
     704        97999 :       dtby2 = dt*.5_dp
     705       195998 :       DO iconst = 1, SIZE(g3x3_list)
     706        97999 :          IF (g3x3_list(iconst)%restraint%active) CYCLE
     707        97999 :          index_a = g3x3_list(iconst)%a + first_atom - 1
     708        97999 :          index_b = g3x3_list(iconst)%b + first_atom - 1
     709        97999 :          index_c = g3x3_list(iconst)%c + first_atom - 1
     710       391996 :          v12(:) = vel(:, index_a) - vel(:, index_b)
     711       391996 :          v13(:) = vel(:, index_a) - vel(:, index_c)
     712       391996 :          v23(:) = vel(:, index_b) - vel(:, index_c)
     713       391996 :          r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
     714       391996 :          r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
     715       391996 :          r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
     716        97999 :          atomic_kind => particle_set(index_a)%atomic_kind
     717        97999 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     718        97999 :          imass1 = 1.0_dp/mass
     719        97999 :          atomic_kind => particle_set(index_b)%atomic_kind
     720        97999 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     721        97999 :          imass2 = 1.0_dp/mass
     722        97999 :          atomic_kind => particle_set(index_c)%atomic_kind
     723        97999 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     724        97999 :          imass3 = 1.0_dp/mass
     725       391996 :          lg3x3(iconst)%fa = -2.0_dp*r12
     726       391996 :          lg3x3(iconst)%fb = -2.0_dp*r13
     727       391996 :          lg3x3(iconst)%fc = -2.0_dp*r23
     728              :          ! Check for fixed atom constraints
     729              :          CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
     730        97999 :                                         index_a, index_b, index_c, fixd_list, lg3x3(iconst))
     731              :          ! roll the fc
     732      1273987 :          f_roll1 = MATMUL(r_rattle, lg3x3(iconst)%fa)
     733      1273987 :          f_roll2 = MATMUL(r_rattle, lg3x3(iconst)%fb)
     734      1273987 :          f_roll3 = MATMUL(r_rattle, lg3x3(iconst)%fc)
     735              : 
     736              :          ! construct matrix
     737       391996 :          amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, f_roll1)
     738       391996 :          amat(1, 2) = imass1*DOT_PRODUCT(r12, f_roll2)
     739       391996 :          amat(1, 3) = -imass2*DOT_PRODUCT(r12, f_roll3)
     740       391996 :          amat(2, 1) = imass1*DOT_PRODUCT(r13, f_roll1)
     741       391996 :          amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, f_roll2)
     742       391996 :          amat(2, 3) = imass3*DOT_PRODUCT(r13, f_roll3)
     743       391996 :          amat(3, 1) = -imass2*DOT_PRODUCT(r23, f_roll1)
     744       391996 :          amat(3, 2) = imass3*DOT_PRODUCT(r23, f_roll2)
     745       391996 :          amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, f_roll3)
     746              : 
     747              :          ! construct solution vector
     748      2155978 :          bvec(1, 1) = DOT_PRODUCT(r12, v12 + MATMUL(veps, r12))
     749      2155978 :          bvec(2, 1) = DOT_PRODUCT(r13, v13 + MATMUL(veps, r13))
     750      2155978 :          bvec(3, 1) = DOT_PRODUCT(r23, v23 + MATMUL(veps, r23))
     751       489995 :          bvec = -bvec*2.0_dp*idt
     752              : 
     753              :          ! get lambda
     754        97999 :          CALL solve_system(amat, 3, bvec)
     755        97999 :          lambda(:) = bvec(:, 1)
     756       391996 :          lg3x3(iconst)%lambda(:) = lambda
     757              : 
     758              :          fc1 = lambda(1)*f_roll1 + &
     759       391996 :                lambda(2)*f_roll2
     760              :          fc2 = -lambda(1)*f_roll1 + &
     761       391996 :                lambda(3)*f_roll3
     762              :          fc3 = -lambda(2)*f_roll2 - &
     763       391996 :                lambda(3)*f_roll3
     764       391996 :          vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
     765       391996 :          vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
     766       587994 :          vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
     767              :       END DO
     768        97999 :    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       143836 :    SUBROUTINE rattle_3x3_low(fixd_list, g3x3_list, lg3x3, first_atom, &
     783       143836 :                              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       143836 :       idt = 1.0_dp/dt
     802       143836 :       dtby2 = dt*.5_dp
     803       287672 :       DO iconst = 1, SIZE(g3x3_list)
     804       143836 :          IF (g3x3_list(iconst)%restraint%active) CYCLE
     805       143816 :          index_a = g3x3_list(iconst)%a + first_atom - 1
     806       143816 :          index_b = g3x3_list(iconst)%b + first_atom - 1
     807       143816 :          index_c = g3x3_list(iconst)%c + first_atom - 1
     808       575264 :          v12(:) = vel(:, index_a) - vel(:, index_b)
     809       575264 :          v13(:) = vel(:, index_a) - vel(:, index_c)
     810       575264 :          v23(:) = vel(:, index_b) - vel(:, index_c)
     811       575264 :          r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
     812       575264 :          r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
     813       575264 :          r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
     814       143816 :          atomic_kind => particle_set(index_a)%atomic_kind
     815       143816 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     816       143816 :          imass1 = 1.0_dp/mass
     817       143816 :          atomic_kind => particle_set(index_b)%atomic_kind
     818       143816 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     819       143816 :          imass2 = 1.0_dp/mass
     820       143816 :          atomic_kind => particle_set(index_c)%atomic_kind
     821       143816 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     822       143816 :          imass3 = 1.0_dp/mass
     823       575264 :          lg3x3(iconst)%fa = -2.0_dp*r12
     824       575264 :          lg3x3(iconst)%fb = -2.0_dp*r13
     825       575264 :          lg3x3(iconst)%fc = -2.0_dp*r23
     826              :          ! Check for fixed atom constraints
     827              :          CALL check_fixed_atom_cns_g3x3(imass1, imass2, imass3, &
     828       143816 :                                         index_a, index_b, index_c, fixd_list, lg3x3(iconst))
     829              :          ! construct matrix
     830       575264 :          amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, lg3x3(iconst)%fa)
     831       575264 :          amat(1, 2) = imass1*DOT_PRODUCT(r12, lg3x3(iconst)%fb)
     832       575264 :          amat(1, 3) = -imass2*DOT_PRODUCT(r12, lg3x3(iconst)%fc)
     833       575264 :          amat(2, 1) = imass1*DOT_PRODUCT(r13, lg3x3(iconst)%fa)
     834       575264 :          amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, lg3x3(iconst)%fb)
     835       575264 :          amat(2, 3) = imass3*DOT_PRODUCT(r13, lg3x3(iconst)%fc)
     836       575264 :          amat(3, 1) = -imass2*DOT_PRODUCT(r23, lg3x3(iconst)%fa)
     837       575264 :          amat(3, 2) = imass3*DOT_PRODUCT(r23, lg3x3(iconst)%fb)
     838       575264 :          amat(3, 3) = (imass2 + imass3)*DOT_PRODUCT(r23, lg3x3(iconst)%fc)
     839              : 
     840              :          ! construct solution vector
     841       575264 :          bvec(1, 1) = DOT_PRODUCT(r12, v12)
     842       575264 :          bvec(2, 1) = DOT_PRODUCT(r13, v13)
     843       575264 :          bvec(3, 1) = DOT_PRODUCT(r23, v23)
     844       719080 :          bvec = -bvec*2.0_dp*idt
     845              : 
     846              :          ! get lambda
     847       143816 :          CALL solve_system(amat, 3, bvec)
     848       575264 :          lg3x3(iconst)%lambda(:) = bvec(:, 1)
     849              : 
     850              :          fc1 = lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
     851       575264 :                lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb
     852              :          fc2 = -lg3x3(iconst)%lambda(1)*lg3x3(iconst)%fa + &
     853       575264 :                lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
     854              :          fc3 = -lg3x3(iconst)%lambda(2)*lg3x3(iconst)%fb - &
     855       575264 :                lg3x3(iconst)%lambda(3)*lg3x3(iconst)%fc
     856       575264 :          vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
     857       575264 :          vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
     858       862936 :          vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
     859              :       END DO
     860       143836 :    END SUBROUTINE rattle_3x3_low
     861              : 
     862       293997 : END MODULE constraint_3x3
        

Generated by: LCOV version 2.0-1