LCOV - code coverage report
Current view: top level - src - constraint_4x6.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.3 % 588 572
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_4x6
      13              : 
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15              :                                               get_atomic_kind
      16              :    USE constraint_fxd,                  ONLY: check_fixed_atom_cns_g4x6
      17              :    USE kinds,                           ONLY: dp
      18              :    USE linear_systems,                  ONLY: solve_system
      19              :    USE molecule_kind_types,             ONLY: fixd_constraint_type,&
      20              :                                               g4x6_constraint_type,&
      21              :                                               get_molecule_kind,&
      22              :                                               molecule_kind_type
      23              :    USE molecule_types,                  ONLY: get_molecule,&
      24              :                                               global_constraint_type,&
      25              :                                               local_g4x6_constraint_type,&
      26              :                                               molecule_type
      27              :    USE particle_types,                  ONLY: particle_type
      28              : #include "./base/base_uses.f90"
      29              : 
      30              :    IMPLICIT NONE
      31              : 
      32              :    PRIVATE
      33              :    PUBLIC :: shake_4x6_int, &
      34              :              rattle_4x6_int, &
      35              :              shake_roll_4x6_int, &
      36              :              rattle_roll_4x6_int, &
      37              :              shake_4x6_ext, &
      38              :              rattle_4x6_ext, &
      39              :              shake_roll_4x6_ext, &
      40              :              rattle_roll_4x6_ext
      41              : 
      42              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_4x6'
      43              : 
      44              : CONTAINS
      45              : 
      46              : ! **************************************************************************************************
      47              : !> \brief Intramolecular shake_4x6
      48              : !> \param molecule ...
      49              : !> \param particle_set ...
      50              : !> \param pos ...
      51              : !> \param vel ...
      52              : !> \param dt ...
      53              : !> \param ishake ...
      54              : !> \param max_sigma ...
      55              : !> \par History
      56              : !>      none
      57              : ! **************************************************************************************************
      58         2466 :    SUBROUTINE shake_4x6_int(molecule, particle_set, pos, vel, dt, ishake, &
      59              :                             max_sigma)
      60              :       TYPE(molecule_type), POINTER                       :: molecule
      61              :       TYPE(particle_type), POINTER                       :: particle_set(:)
      62              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
      63              :       REAL(kind=dp), INTENT(in)                          :: dt
      64              :       INTEGER, INTENT(IN)                                :: ishake
      65              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
      66              : 
      67              :       INTEGER                                            :: first_atom, ng4x6
      68         2466 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
      69         2466 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
      70         2466 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
      71              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      72              : 
      73         2466 :       NULLIFY (fixd_list)
      74         2466 :       molecule_kind => molecule%molecule_kind
      75              :       CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
      76         2466 :                              fixd_list=fixd_list, g4x6_list=g4x6_list)
      77         2466 :       CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
      78              :       ! Real Shake
      79              :       CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
      80         2466 :                          particle_set, pos, vel, dt, ishake, max_sigma)
      81              : 
      82         2466 :    END SUBROUTINE shake_4x6_int
      83              : 
      84              : ! **************************************************************************************************
      85              : !> \brief Intramolecular shake_4x6_roll
      86              : !> \param molecule ...
      87              : !> \param particle_set ...
      88              : !> \param pos ...
      89              : !> \param vel ...
      90              : !> \param r_shake ...
      91              : !> \param dt ...
      92              : !> \param ishake ...
      93              : !> \param max_sigma ...
      94              : !> \par History
      95              : !>      none
      96              : ! **************************************************************************************************
      97         2225 :    SUBROUTINE shake_roll_4x6_int(molecule, particle_set, pos, vel, r_shake, &
      98              :                                  dt, ishake, max_sigma)
      99              :       TYPE(molecule_type), POINTER                       :: molecule
     100              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     101              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     102              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake
     103              :       REAL(kind=dp), INTENT(in)                          :: dt
     104              :       INTEGER, INTENT(IN)                                :: ishake
     105              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     106              : 
     107              :       INTEGER                                            :: first_atom, ng4x6
     108         2225 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     109         2225 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     110         2225 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     111              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     112              : 
     113         2225 :       NULLIFY (fixd_list)
     114         2225 :       molecule_kind => molecule%molecule_kind
     115              :       CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
     116         2225 :                              fixd_list=fixd_list, g4x6_list=g4x6_list)
     117         2225 :       CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
     118              :       ! Real Shake
     119              :       CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
     120         2225 :                               particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
     121              : 
     122         2225 :    END SUBROUTINE shake_roll_4x6_int
     123              : 
     124              : ! **************************************************************************************************
     125              : !> \brief Intramolecular rattle_4x6
     126              : !> \param molecule ...
     127              : !> \param particle_set ...
     128              : !> \param vel ...
     129              : !> \param dt ...
     130              : !> \par History
     131              : !>      none
     132              : ! **************************************************************************************************
     133          682 :    SUBROUTINE rattle_4x6_int(molecule, particle_set, vel, dt)
     134              :       TYPE(molecule_type), POINTER                       :: molecule
     135              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     136              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     137              :       REAL(kind=dp), INTENT(in)                          :: dt
     138              : 
     139              :       INTEGER                                            :: first_atom, ng4x6
     140          682 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     141          682 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     142          682 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     143              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     144              : 
     145          682 :       NULLIFY (fixd_list)
     146          682 :       molecule_kind => molecule%molecule_kind
     147              :       CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
     148          682 :                              fixd_list=fixd_list, g4x6_list=g4x6_list)
     149          682 :       CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
     150              :       ! Real Rattle
     151              :       CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
     152          682 :                           particle_set, vel, dt)
     153              : 
     154          682 :    END SUBROUTINE rattle_4x6_int
     155              : 
     156              : ! **************************************************************************************************
     157              : !> \brief Intramolecular rattle_4x6_roll
     158              : !> \param molecule ...
     159              : !> \param particle_set ...
     160              : !> \param vel ...
     161              : !> \param r_rattle ...
     162              : !> \param dt ...
     163              : !> \param veps ...
     164              : !> \par History
     165              : !>      none
     166              : ! **************************************************************************************************
     167         1024 :    SUBROUTINE rattle_roll_4x6_int(molecule, particle_set, vel, r_rattle, dt, veps)
     168              :       TYPE(molecule_type), POINTER                       :: molecule
     169              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     170              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     171              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_rattle
     172              :       REAL(kind=dp), INTENT(in)                          :: dt
     173              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: veps
     174              : 
     175              :       INTEGER                                            :: first_atom, ng4x6
     176         1024 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     177         1024 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     178         1024 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     179              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     180              : 
     181         1024 :       NULLIFY (fixd_list)
     182         1024 :       molecule_kind => molecule%molecule_kind
     183              :       CALL get_molecule_kind(molecule_kind, ng4x6=ng4x6, &
     184         1024 :                              fixd_list=fixd_list, g4x6_list=g4x6_list)
     185         1024 :       CALL get_molecule(molecule, first_atom=first_atom, lg4x6=lg4x6)
     186              :       ! Real Rattle
     187              :       CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
     188         1024 :                                particle_set, vel, r_rattle, dt, veps)
     189              : 
     190         1024 :    END SUBROUTINE rattle_roll_4x6_int
     191              : 
     192              : ! **************************************************************************************************
     193              : !> \brief Intramolecular shake_4x6
     194              : !> \param gci ...
     195              : !> \param particle_set ...
     196              : !> \param pos ...
     197              : !> \param vel ...
     198              : !> \param dt ...
     199              : !> \param ishake ...
     200              : !> \param max_sigma ...
     201              : !> \par History
     202              : !>      none
     203              : ! **************************************************************************************************
     204           48 :    SUBROUTINE shake_4x6_ext(gci, particle_set, pos, vel, dt, ishake, &
     205              :                             max_sigma)
     206              : 
     207              :       TYPE(global_constraint_type), POINTER              :: gci
     208              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     209              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     210              :       REAL(kind=dp), INTENT(in)                          :: dt
     211              :       INTEGER, INTENT(IN)                                :: ishake
     212              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     213              : 
     214              :       INTEGER                                            :: first_atom, ng4x6
     215              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     216              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     217              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     218              : 
     219           48 :       first_atom = 1
     220           48 :       ng4x6 = gci%ng4x6
     221           48 :       fixd_list => gci%fixd_list
     222           48 :       g4x6_list => gci%g4x6_list
     223           48 :       lg4x6 => gci%lg4x6
     224              :       ! Real Shake
     225              :       CALL shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
     226           48 :                          particle_set, pos, vel, dt, ishake, max_sigma)
     227              : 
     228           48 :    END SUBROUTINE shake_4x6_ext
     229              : 
     230              : ! **************************************************************************************************
     231              : !> \brief Intramolecular shake_4x6_roll
     232              : !> \param gci ...
     233              : !> \param particle_set ...
     234              : !> \param pos ...
     235              : !> \param vel ...
     236              : !> \param r_shake ...
     237              : !> \param dt ...
     238              : !> \param ishake ...
     239              : !> \param max_sigma ...
     240              : !> \par History
     241              : !>      none
     242              : ! **************************************************************************************************
     243            0 :    SUBROUTINE shake_roll_4x6_ext(gci, particle_set, pos, vel, r_shake, &
     244              :                                  dt, ishake, max_sigma)
     245              : 
     246              :       TYPE(global_constraint_type), POINTER              :: gci
     247              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     248              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     249              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake
     250              :       REAL(kind=dp), INTENT(in)                          :: dt
     251              :       INTEGER, INTENT(IN)                                :: ishake
     252              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     253              : 
     254              :       INTEGER                                            :: first_atom, ng4x6
     255              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     256              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     257              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     258              : 
     259            0 :       first_atom = 1
     260            0 :       ng4x6 = gci%ng4x6
     261            0 :       fixd_list => gci%fixd_list
     262            0 :       g4x6_list => gci%g4x6_list
     263            0 :       lg4x6 => gci%lg4x6
     264              :       ! Real Shake
     265              :       CALL shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
     266            0 :                               particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
     267              : 
     268            0 :    END SUBROUTINE shake_roll_4x6_ext
     269              : 
     270              : ! **************************************************************************************************
     271              : !> \brief Intramolecular rattle_4x6
     272              : !> \param gci ...
     273              : !> \param particle_set ...
     274              : !> \param vel ...
     275              : !> \param dt ...
     276              : !> \par History
     277              : !>      none
     278              : ! **************************************************************************************************
     279           20 :    SUBROUTINE rattle_4x6_ext(gci, particle_set, vel, dt)
     280              : 
     281              :       TYPE(global_constraint_type), POINTER              :: gci
     282              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     283              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     284              :       REAL(kind=dp), INTENT(in)                          :: dt
     285              : 
     286              :       INTEGER                                            :: first_atom, ng4x6
     287              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     288              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     289              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     290              : 
     291           20 :       first_atom = 1
     292           20 :       ng4x6 = gci%ng4x6
     293           20 :       fixd_list => gci%fixd_list
     294           20 :       g4x6_list => gci%g4x6_list
     295           20 :       lg4x6 => gci%lg4x6
     296              :       ! Real Rattle
     297              :       CALL rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
     298           20 :                           particle_set, vel, dt)
     299              : 
     300           20 :    END SUBROUTINE rattle_4x6_ext
     301              : 
     302              : ! **************************************************************************************************
     303              : !> \brief Intramolecular rattle_4x6_roll
     304              : !> \param gci ...
     305              : !> \param particle_set ...
     306              : !> \param vel ...
     307              : !> \param r_rattle ...
     308              : !> \param dt ...
     309              : !> \param veps ...
     310              : !> \par History
     311              : !>      none
     312              : ! **************************************************************************************************
     313            0 :    SUBROUTINE rattle_roll_4x6_ext(gci, particle_set, vel, r_rattle, dt, veps)
     314              : 
     315              :       TYPE(global_constraint_type), POINTER              :: gci
     316              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     317              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     318              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_rattle
     319              :       REAL(kind=dp), INTENT(in)                          :: dt
     320              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: veps
     321              : 
     322              :       INTEGER                                            :: first_atom, ng4x6
     323              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     324              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     325              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     326              : 
     327            0 :       first_atom = 1
     328            0 :       ng4x6 = gci%ng4x6
     329            0 :       fixd_list => gci%fixd_list
     330            0 :       g4x6_list => gci%g4x6_list
     331            0 :       lg4x6 => gci%lg4x6
     332              :       ! Real Rattle
     333              :       CALL rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
     334            0 :                                particle_set, vel, r_rattle, dt, veps)
     335              : 
     336            0 :    END SUBROUTINE rattle_roll_4x6_ext
     337              : 
     338              : ! **************************************************************************************************
     339              : !> \brief ...
     340              : !> \param fixd_list ...
     341              : !> \param g4x6_list ...
     342              : !> \param lg4x6 ...
     343              : !> \param ng4x6 ...
     344              : !> \param first_atom ...
     345              : !> \param particle_set ...
     346              : !> \param pos ...
     347              : !> \param vel ...
     348              : !> \param dt ...
     349              : !> \param ishake ...
     350              : !> \param max_sigma ...
     351              : !> \par History
     352              : !>      none
     353              : ! **************************************************************************************************
     354         2514 :    SUBROUTINE shake_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
     355         2514 :                             particle_set, pos, vel, dt, ishake, max_sigma)
     356              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     357              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     358              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     359              :       INTEGER, INTENT(IN)                                :: ng4x6, first_atom
     360              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     361              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     362              :       REAL(kind=dp), INTENT(in)                          :: dt
     363              :       INTEGER, INTENT(IN)                                :: ishake
     364              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     365              : 
     366              :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
     367              :                                                             index_d
     368              :       REAL(KIND=dp)                                      :: dtby2, dtsqby2, idtsq, imass1, imass2, &
     369              :                                                             imass3, imass4, sigma
     370              :       REAL(KIND=dp), DIMENSION(3)                        :: fc1, fc2, fc3, fc4, r0_12, r0_13, r0_14, &
     371              :                                                             r0_23, r0_24, r0_34, r1, r12, r13, &
     372              :                                                             r14, r2, r23, r24, r3, r34, r4, v1, &
     373              :                                                             v2, v3, v4, vec
     374              :       REAL(KIND=dp), DIMENSION(6, 1)                     :: bvec
     375              :       REAL(KIND=dp), DIMENSION(6, 6)                     :: amat, atemp
     376              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     377              : 
     378         2514 :       dtsqby2 = dt*dt*.5_dp
     379         2514 :       idtsq = 1.0_dp/dt/dt
     380         2514 :       dtby2 = dt*.5_dp
     381         5028 :       DO iconst = 1, ng4x6
     382         2514 :          IF (g4x6_list(iconst)%restraint%active) CYCLE
     383         2504 :          index_a = g4x6_list(iconst)%a + first_atom - 1
     384         2504 :          index_b = g4x6_list(iconst)%b + first_atom - 1
     385         2504 :          index_c = g4x6_list(iconst)%c + first_atom - 1
     386         2504 :          index_d = g4x6_list(iconst)%d + first_atom - 1
     387         2504 :          IF (ishake == 1) THEN
     388         2768 :             r0_12(:) = pos(:, index_a) - pos(:, index_b)
     389         2768 :             r0_13(:) = pos(:, index_a) - pos(:, index_c)
     390         2768 :             r0_14(:) = pos(:, index_a) - pos(:, index_d)
     391         2768 :             r0_23(:) = pos(:, index_b) - pos(:, index_c)
     392         2768 :             r0_24(:) = pos(:, index_b) - pos(:, index_d)
     393         2768 :             r0_34(:) = pos(:, index_c) - pos(:, index_d)
     394          692 :             atomic_kind => particle_set(index_a)%atomic_kind
     395          692 :             imass1 = 1.0_dp/atomic_kind%mass
     396          692 :             atomic_kind => particle_set(index_b)%atomic_kind
     397          692 :             imass2 = 1.0_dp/atomic_kind%mass
     398          692 :             atomic_kind => particle_set(index_c)%atomic_kind
     399          692 :             imass3 = 1.0_dp/atomic_kind%mass
     400          692 :             atomic_kind => particle_set(index_d)%atomic_kind
     401          692 :             imass4 = 1.0_dp/atomic_kind%mass
     402              :             lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
     403         2768 :                                         lg4x6(iconst)%rb_old)
     404              :             lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
     405         2768 :                                         lg4x6(iconst)%rc_old)
     406              :             lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
     407         2768 :                                         lg4x6(iconst)%rd_old)
     408              :             lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
     409         2768 :                                         lg4x6(iconst)%rc_old)
     410              :             lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
     411         2768 :                                         lg4x6(iconst)%rd_old)
     412              :             lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
     413         2768 :                                         lg4x6(iconst)%rd_old)
     414              :             ! Check for fixed atom constraints
     415              :             CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
     416          692 :                                            index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
     417              :             ! construct matrix
     418         2768 :             amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, lg4x6(iconst)%fa)
     419         2768 :             amat(1, 2) = imass1*DOT_PRODUCT(r0_12, lg4x6(iconst)%fb)
     420         2768 :             amat(1, 3) = imass1*DOT_PRODUCT(r0_12, lg4x6(iconst)%fc)
     421         2768 :             amat(1, 4) = -imass2*DOT_PRODUCT(r0_12, lg4x6(iconst)%fd)
     422         2768 :             amat(1, 5) = -imass2*DOT_PRODUCT(r0_12, lg4x6(iconst)%fe)
     423          692 :             amat(1, 6) = 0.0_dp
     424              : 
     425         2768 :             amat(2, 1) = imass1*DOT_PRODUCT(r0_13, lg4x6(iconst)%fa)
     426         2768 :             amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, lg4x6(iconst)%fb)
     427         2768 :             amat(2, 3) = imass1*DOT_PRODUCT(r0_13, lg4x6(iconst)%fc)
     428         2768 :             amat(2, 4) = imass3*DOT_PRODUCT(r0_13, lg4x6(iconst)%fd)
     429          692 :             amat(2, 5) = 0.0_dp
     430         2768 :             amat(2, 6) = -imass3*DOT_PRODUCT(r0_13, lg4x6(iconst)%ff)
     431              : 
     432         2768 :             amat(3, 1) = imass1*DOT_PRODUCT(r0_14, lg4x6(iconst)%fa)
     433         2768 :             amat(3, 2) = imass1*DOT_PRODUCT(r0_14, lg4x6(iconst)%fb)
     434         2768 :             amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r0_14, lg4x6(iconst)%fc)
     435          692 :             amat(3, 4) = 0.0_dp
     436         2768 :             amat(3, 5) = imass4*DOT_PRODUCT(r0_14, lg4x6(iconst)%fe)
     437         2768 :             amat(3, 6) = imass4*DOT_PRODUCT(r0_14, lg4x6(iconst)%ff)
     438              : 
     439         2768 :             amat(4, 1) = -imass2*DOT_PRODUCT(r0_23, lg4x6(iconst)%fa)
     440         2768 :             amat(4, 2) = imass3*DOT_PRODUCT(r0_23, lg4x6(iconst)%fb)
     441          692 :             amat(4, 3) = 0.0_dp
     442         2768 :             amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r0_23, lg4x6(iconst)%fd)
     443         2768 :             amat(4, 5) = imass2*DOT_PRODUCT(r0_23, lg4x6(iconst)%fe)
     444         2768 :             amat(4, 6) = -imass3*DOT_PRODUCT(r0_23, lg4x6(iconst)%ff)
     445              : 
     446         2768 :             amat(5, 1) = -imass2*DOT_PRODUCT(r0_24, lg4x6(iconst)%fa)
     447          692 :             amat(5, 2) = 0.0_dp
     448         2768 :             amat(5, 3) = imass4*DOT_PRODUCT(r0_24, lg4x6(iconst)%fc)
     449         2768 :             amat(5, 4) = imass2*DOT_PRODUCT(r0_24, lg4x6(iconst)%fd)
     450         2768 :             amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r0_24, lg4x6(iconst)%fe)
     451         2768 :             amat(5, 6) = imass4*DOT_PRODUCT(r0_24, lg4x6(iconst)%ff)
     452              : 
     453          692 :             amat(6, 1) = 0.0_dp
     454         2768 :             amat(6, 2) = -imass3*DOT_PRODUCT(r0_34, lg4x6(iconst)%fb)
     455         2768 :             amat(6, 3) = imass4*DOT_PRODUCT(r0_34, lg4x6(iconst)%fc)
     456         2768 :             amat(6, 4) = -imass3*DOT_PRODUCT(r0_34, lg4x6(iconst)%fd)
     457         2768 :             amat(6, 5) = imass4*DOT_PRODUCT(r0_34, lg4x6(iconst)%fe)
     458         2768 :             amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r0_34, lg4x6(iconst)%ff)
     459              :             ! Store values
     460         2768 :             lg4x6(iconst)%r0_12 = r0_12
     461         2768 :             lg4x6(iconst)%r0_13 = r0_13
     462         2768 :             lg4x6(iconst)%r0_14 = r0_14
     463         2768 :             lg4x6(iconst)%r0_23 = r0_23
     464         2768 :             lg4x6(iconst)%r0_24 = r0_24
     465         2768 :             lg4x6(iconst)%r0_34 = r0_34
     466        29756 :             lg4x6(iconst)%amat = amat
     467          692 :             lg4x6(iconst)%imass1 = imass1
     468          692 :             lg4x6(iconst)%imass2 = imass2
     469          692 :             lg4x6(iconst)%imass3 = imass3
     470          692 :             lg4x6(iconst)%imass4 = imass4
     471         4844 :             lg4x6(iconst)%lambda_old = 0.0_dp
     472         4844 :             lg4x6(iconst)%del_lambda = 0.0_dp
     473              :          ELSE
     474              :             ! Retrieve values
     475         7248 :             r0_12 = lg4x6(iconst)%r0_12
     476         7248 :             r0_13 = lg4x6(iconst)%r0_13
     477         7248 :             r0_14 = lg4x6(iconst)%r0_14
     478         7248 :             r0_23 = lg4x6(iconst)%r0_23
     479         7248 :             r0_24 = lg4x6(iconst)%r0_24
     480         7248 :             r0_34 = lg4x6(iconst)%r0_34
     481        77916 :             amat = lg4x6(iconst)%amat
     482         1812 :             imass1 = lg4x6(iconst)%imass1
     483         1812 :             imass2 = lg4x6(iconst)%imass2
     484         1812 :             imass3 = lg4x6(iconst)%imass3
     485         1812 :             imass4 = lg4x6(iconst)%imass4
     486              :          END IF
     487              : 
     488              :          ! Iterate until convergence:
     489              :          vec = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa*(imass1 + imass2) + &
     490              :                lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + &
     491              :                lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc - &
     492              :                lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd - &
     493        10016 :                lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe
     494              :          bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
     495        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
     496              :          vec = lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb*(imass1 + imass3) + &
     497              :                lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + &
     498              :                lg4x6(iconst)%lambda(3)*imass1*lg4x6(iconst)%fc + &
     499              :                lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd - &
     500        10016 :                lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
     501              :          bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
     502        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
     503              :          vec = lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc*(imass1 + imass4) + &
     504              :                lg4x6(iconst)%lambda(1)*imass1*lg4x6(iconst)%fa + &
     505              :                lg4x6(iconst)%lambda(2)*imass1*lg4x6(iconst)%fb + &
     506              :                lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe + &
     507        10016 :                lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
     508              :          bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
     509        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_14, r0_14)
     510              :          vec = lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd*(imass2 + imass3) - &
     511              :                lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + &
     512              :                lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + &
     513              :                lg4x6(iconst)%lambda(5)*imass2*lg4x6(iconst)%fe - &
     514        10016 :                lg4x6(iconst)%lambda(6)*imass3*lg4x6(iconst)%ff
     515              :          bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
     516        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
     517              :          vec = lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe*(imass2 + imass4) - &
     518              :                lg4x6(iconst)%lambda(1)*imass2*lg4x6(iconst)%fa + &
     519              :                lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc + &
     520              :                lg4x6(iconst)%lambda(4)*imass2*lg4x6(iconst)%fd + &
     521        10016 :                lg4x6(iconst)%lambda(6)*imass4*lg4x6(iconst)%ff
     522              :          bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
     523        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_24, r0_24)
     524              :          vec = lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff*(imass3 + imass4) - &
     525              :                lg4x6(iconst)%lambda(2)*imass3*lg4x6(iconst)%fb + &
     526              :                lg4x6(iconst)%lambda(3)*imass4*lg4x6(iconst)%fc - &
     527              :                lg4x6(iconst)%lambda(4)*imass3*lg4x6(iconst)%fd + &
     528        10016 :                lg4x6(iconst)%lambda(5)*imass4*lg4x6(iconst)%fe
     529              :          bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
     530        20032 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_34, r0_34)
     531        20032 :          bvec = bvec*idtsq
     532              : 
     533              :          ! get lambda
     534         2504 :          atemp = amat
     535         2504 :          CALL solve_system(atemp, 6, bvec)
     536        17528 :          lg4x6(iconst)%lambda(:) = bvec(:, 1)
     537              :          lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
     538        17528 :                                        lg4x6(iconst)%lambda_old(:)
     539        17528 :          lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:)
     540              : 
     541              :          fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
     542              :                lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + &
     543        10016 :                lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc
     544              :          fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
     545              :                lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
     546        10016 :                lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe
     547              :          fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - &
     548              :                lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
     549        10016 :                lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
     550              :          fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - &
     551              :                lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - &
     552        10016 :                lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
     553        10016 :          r1(:) = pos(:, index_a) + imass1*dtsqby2*fc1(:)
     554        10016 :          r2(:) = pos(:, index_b) + imass2*dtsqby2*fc2(:)
     555        10016 :          r3(:) = pos(:, index_c) + imass3*dtsqby2*fc3(:)
     556        10016 :          r4(:) = pos(:, index_d) + imass4*dtsqby2*fc4(:)
     557        10016 :          v1(:) = vel(:, index_a) + imass1*dtby2*fc1(:)
     558        10016 :          v2(:) = vel(:, index_b) + imass2*dtby2*fc2(:)
     559        10016 :          v3(:) = vel(:, index_c) + imass3*dtby2*fc3(:)
     560        10016 :          v4(:) = vel(:, index_d) + imass4*dtby2*fc4(:)
     561        10016 :          r12 = r1 - r2
     562        10016 :          r13 = r1 - r3
     563        10016 :          r14 = r1 - r4
     564        10016 :          r23 = r2 - r3
     565        10016 :          r24 = r2 - r4
     566        10016 :          r34 = r3 - r4
     567              :          ! compute the tolerance:
     568              :          sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* &
     569        10016 :                  g4x6_list(iconst)%dab
     570         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     571              :          sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* &
     572        10016 :                  g4x6_list(iconst)%dac
     573         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     574              :          sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* &
     575        10016 :                  g4x6_list(iconst)%dad
     576         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     577              :          sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* &
     578        10016 :                  g4x6_list(iconst)%dbc
     579         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     580              :          sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* &
     581        10016 :                  g4x6_list(iconst)%dbd
     582         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     583              :          sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* &
     584        10016 :                  g4x6_list(iconst)%dcd
     585         2504 :          max_sigma = MAX(max_sigma, ABS(sigma))
     586              : 
     587              :          ! update positions with full multiplier
     588        10016 :          pos(:, index_a) = r1(:)
     589        10016 :          pos(:, index_b) = r2(:)
     590        10016 :          pos(:, index_c) = r3(:)
     591        10016 :          pos(:, index_d) = r4(:)
     592              : 
     593              :          ! update velocites with full multiplier
     594        10016 :          vel(:, index_a) = v1(:)
     595        10016 :          vel(:, index_b) = v2(:)
     596        10016 :          vel(:, index_c) = v3(:)
     597        12540 :          vel(:, index_d) = v4(:)
     598              :       END DO
     599         2514 :    END SUBROUTINE shake_4x6_low
     600              : 
     601              : ! **************************************************************************************************
     602              : !> \brief ...
     603              : !> \param fixd_list ...
     604              : !> \param g4x6_list ...
     605              : !> \param lg4x6 ...
     606              : !> \param ng4x6 ...
     607              : !> \param first_atom ...
     608              : !> \param particle_set ...
     609              : !> \param pos ...
     610              : !> \param vel ...
     611              : !> \param r_shake ...
     612              : !> \param dt ...
     613              : !> \param ishake ...
     614              : !> \param max_sigma ...
     615              : !> \par History
     616              : !>      none
     617              : ! **************************************************************************************************
     618         2225 :    SUBROUTINE shake_roll_4x6_low(fixd_list, g4x6_list, lg4x6, ng4x6, first_atom, &
     619         2225 :                                  particle_set, pos, vel, r_shake, dt, ishake, max_sigma)
     620              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     621              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     622              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     623              :       INTEGER, INTENT(IN)                                :: ng4x6, first_atom
     624              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     625              :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     626              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake
     627              :       REAL(kind=dp), INTENT(in)                          :: dt
     628              :       INTEGER, INTENT(IN)                                :: ishake
     629              :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     630              : 
     631              :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
     632              :                                                             index_d
     633              :       REAL(KIND=dp)                                      :: dtby2, dtsqby2, idtsq, imass1, imass2, &
     634              :                                                             imass3, imass4, sigma
     635              :       REAL(KIND=dp), DIMENSION(3) :: f_roll1, f_roll2, f_roll3, f_roll4, f_roll5, f_roll6, fc1, &
     636              :          fc2, fc3, fc4, r0_12, r0_13, r0_14, r0_23, r0_24, r0_34, r1, r12, r13, r14, r2, r23, r24, &
     637              :          r3, r34, r4, v1, v2, v3, v4, vec
     638              :       REAL(KIND=dp), DIMENSION(6, 1)                     :: bvec
     639              :       REAL(KIND=dp), DIMENSION(6, 6)                     :: amat, atemp
     640              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     641              : 
     642         2225 :       dtsqby2 = dt*dt*.5_dp
     643         2225 :       idtsq = 1.0_dp/dt/dt
     644         2225 :       dtby2 = dt*.5_dp
     645         4450 :       DO iconst = 1, ng4x6
     646         2225 :          IF (g4x6_list(iconst)%restraint%active) CYCLE
     647         2225 :          index_a = g4x6_list(iconst)%a + first_atom - 1
     648         2225 :          index_b = g4x6_list(iconst)%b + first_atom - 1
     649         2225 :          index_c = g4x6_list(iconst)%c + first_atom - 1
     650         2225 :          index_d = g4x6_list(iconst)%d + first_atom - 1
     651         2225 :          IF (ishake == 1) THEN
     652         4352 :             r0_12(:) = pos(:, index_a) - pos(:, index_b)
     653         4352 :             r0_13(:) = pos(:, index_a) - pos(:, index_c)
     654         4352 :             r0_23(:) = pos(:, index_b) - pos(:, index_c)
     655         4352 :             r0_14(:) = pos(:, index_a) - pos(:, index_d)
     656         4352 :             r0_24(:) = pos(:, index_b) - pos(:, index_d)
     657         4352 :             r0_34(:) = pos(:, index_c) - pos(:, index_d)
     658         1088 :             atomic_kind => particle_set(index_a)%atomic_kind
     659         1088 :             imass1 = 1.0_dp/atomic_kind%mass
     660         1088 :             atomic_kind => particle_set(index_b)%atomic_kind
     661         1088 :             imass2 = 1.0_dp/atomic_kind%mass
     662         1088 :             atomic_kind => particle_set(index_c)%atomic_kind
     663         1088 :             imass3 = 1.0_dp/atomic_kind%mass
     664         1088 :             atomic_kind => particle_set(index_d)%atomic_kind
     665         1088 :             imass4 = 1.0_dp/atomic_kind%mass
     666              :             lg4x6(iconst)%fa = -2.0_dp*(lg4x6(iconst)%ra_old - &
     667         4352 :                                         lg4x6(iconst)%rb_old)
     668              :             lg4x6(iconst)%fb = -2.0_dp*(lg4x6(iconst)%ra_old - &
     669         4352 :                                         lg4x6(iconst)%rc_old)
     670              :             lg4x6(iconst)%fc = -2.0_dp*(lg4x6(iconst)%ra_old - &
     671         4352 :                                         lg4x6(iconst)%rd_old)
     672              :             lg4x6(iconst)%fd = -2.0_dp*(lg4x6(iconst)%rb_old - &
     673         4352 :                                         lg4x6(iconst)%rc_old)
     674              :             lg4x6(iconst)%fe = -2.0_dp*(lg4x6(iconst)%rb_old - &
     675         4352 :                                         lg4x6(iconst)%rd_old)
     676              :             lg4x6(iconst)%ff = -2.0_dp*(lg4x6(iconst)%rc_old - &
     677         4352 :                                         lg4x6(iconst)%rd_old)
     678              :             ! Check for fixed atom constraints
     679              :             CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
     680         1088 :                                            index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
     681              :             ! rotate fconst:
     682        14144 :             f_roll1 = MATMUL(r_shake, lg4x6(iconst)%fa)
     683        14144 :             f_roll2 = MATMUL(r_shake, lg4x6(iconst)%fb)
     684        14144 :             f_roll3 = MATMUL(r_shake, lg4x6(iconst)%fc)
     685        14144 :             f_roll4 = MATMUL(r_shake, lg4x6(iconst)%fd)
     686        14144 :             f_roll5 = MATMUL(r_shake, lg4x6(iconst)%fe)
     687        14144 :             f_roll6 = MATMUL(r_shake, lg4x6(iconst)%ff)
     688              : 
     689              :             ! construct matrix
     690         4352 :             amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r0_12, f_roll1)
     691         4352 :             amat(1, 2) = imass1*DOT_PRODUCT(r0_12, f_roll2)
     692         4352 :             amat(1, 3) = imass1*DOT_PRODUCT(r0_12, f_roll3)
     693         4352 :             amat(1, 4) = -imass2*DOT_PRODUCT(r0_12, f_roll4)
     694         4352 :             amat(1, 5) = -imass2*DOT_PRODUCT(r0_12, f_roll5)
     695         1088 :             amat(1, 6) = 0.0_dp
     696              : 
     697         4352 :             amat(2, 1) = imass1*DOT_PRODUCT(r0_13, f_roll1)
     698         4352 :             amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r0_13, f_roll2)
     699         4352 :             amat(2, 3) = imass1*DOT_PRODUCT(r0_13, f_roll3)
     700         4352 :             amat(2, 4) = imass3*DOT_PRODUCT(r0_13, f_roll4)
     701         1088 :             amat(2, 5) = 0.0_dp
     702         4352 :             amat(2, 6) = -imass3*DOT_PRODUCT(r0_13, f_roll6)
     703              : 
     704         4352 :             amat(3, 1) = imass1*DOT_PRODUCT(r0_14, f_roll1)
     705         4352 :             amat(3, 2) = imass1*DOT_PRODUCT(r0_14, f_roll2)
     706         4352 :             amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r0_14, f_roll3)
     707         1088 :             amat(3, 4) = 0.0_dp
     708         4352 :             amat(3, 5) = imass4*DOT_PRODUCT(r0_14, f_roll5)
     709         4352 :             amat(3, 6) = imass4*DOT_PRODUCT(r0_14, f_roll6)
     710              : 
     711         4352 :             amat(4, 1) = -imass2*DOT_PRODUCT(r0_23, f_roll1)
     712         4352 :             amat(4, 2) = imass3*DOT_PRODUCT(r0_23, f_roll2)
     713         1088 :             amat(4, 3) = 0.0_dp
     714         4352 :             amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r0_23, f_roll4)
     715         4352 :             amat(4, 5) = imass2*DOT_PRODUCT(r0_23, f_roll5)
     716         4352 :             amat(4, 6) = -imass3*DOT_PRODUCT(r0_23, f_roll6)
     717              : 
     718         4352 :             amat(5, 1) = -imass2*DOT_PRODUCT(r0_24, f_roll1)
     719         1088 :             amat(5, 2) = 0.0_dp
     720         4352 :             amat(5, 3) = imass4*DOT_PRODUCT(r0_24, f_roll3)
     721         4352 :             amat(5, 4) = imass2*DOT_PRODUCT(r0_24, f_roll4)
     722         4352 :             amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r0_24, f_roll5)
     723         4352 :             amat(5, 6) = imass4*DOT_PRODUCT(r0_24, f_roll6)
     724              : 
     725         1088 :             amat(6, 1) = 0.0_dp
     726         4352 :             amat(6, 2) = -imass3*DOT_PRODUCT(r0_34, f_roll2)
     727         4352 :             amat(6, 3) = imass4*DOT_PRODUCT(r0_34, f_roll3)
     728         4352 :             amat(6, 4) = -imass3*DOT_PRODUCT(r0_34, f_roll4)
     729         4352 :             amat(6, 5) = imass4*DOT_PRODUCT(r0_34, f_roll5)
     730         4352 :             amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r0_34, f_roll6)
     731              :             ! Store values
     732         4352 :             lg4x6(iconst)%r0_12 = r0_12
     733         4352 :             lg4x6(iconst)%r0_13 = r0_13
     734         4352 :             lg4x6(iconst)%r0_14 = r0_14
     735         4352 :             lg4x6(iconst)%r0_23 = r0_23
     736         4352 :             lg4x6(iconst)%r0_24 = r0_24
     737         4352 :             lg4x6(iconst)%r0_34 = r0_34
     738         4352 :             lg4x6(iconst)%f_roll1 = f_roll1
     739         4352 :             lg4x6(iconst)%f_roll2 = f_roll2
     740         4352 :             lg4x6(iconst)%f_roll3 = f_roll3
     741         4352 :             lg4x6(iconst)%f_roll4 = f_roll4
     742         4352 :             lg4x6(iconst)%f_roll5 = f_roll5
     743         4352 :             lg4x6(iconst)%f_roll6 = f_roll6
     744        46784 :             lg4x6(iconst)%amat = amat
     745         1088 :             lg4x6(iconst)%imass1 = imass1
     746         1088 :             lg4x6(iconst)%imass2 = imass2
     747         1088 :             lg4x6(iconst)%imass3 = imass3
     748         1088 :             lg4x6(iconst)%imass4 = imass4
     749         7616 :             lg4x6(iconst)%lambda_old = 0.0_dp
     750         7616 :             lg4x6(iconst)%del_lambda = 0.0_dp
     751              :          ELSE
     752              :             ! Retrieve values
     753         4548 :             r0_12 = lg4x6(iconst)%r0_12
     754         4548 :             r0_13 = lg4x6(iconst)%r0_13
     755         4548 :             r0_14 = lg4x6(iconst)%r0_14
     756         4548 :             r0_23 = lg4x6(iconst)%r0_23
     757         4548 :             r0_24 = lg4x6(iconst)%r0_24
     758         4548 :             r0_34 = lg4x6(iconst)%r0_34
     759         4548 :             f_roll1 = lg4x6(iconst)%f_roll1
     760         4548 :             f_roll2 = lg4x6(iconst)%f_roll2
     761         4548 :             f_roll3 = lg4x6(iconst)%f_roll3
     762         4548 :             f_roll4 = lg4x6(iconst)%f_roll4
     763         4548 :             f_roll5 = lg4x6(iconst)%f_roll5
     764         4548 :             f_roll6 = lg4x6(iconst)%f_roll6
     765        48891 :             amat = lg4x6(iconst)%amat
     766         1137 :             imass1 = lg4x6(iconst)%imass1
     767         1137 :             imass2 = lg4x6(iconst)%imass2
     768         1137 :             imass3 = lg4x6(iconst)%imass3
     769         1137 :             imass4 = lg4x6(iconst)%imass4
     770              :          END IF
     771              : 
     772              :          ! Iterate until convergence:
     773              :          vec = lg4x6(iconst)%lambda(1)*f_roll1*(imass1 + imass2) + &
     774              :                lg4x6(iconst)%lambda(2)*imass1*f_roll2 + &
     775              :                lg4x6(iconst)%lambda(3)*imass1*f_roll3 - &
     776              :                lg4x6(iconst)%lambda(4)*imass2*f_roll4 - &
     777         8900 :                lg4x6(iconst)%lambda(5)*imass2*f_roll5
     778              :          bvec(1, 1) = g4x6_list(iconst)%dab*g4x6_list(iconst)%dab &
     779        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_12, r0_12)
     780              :          vec = lg4x6(iconst)%lambda(2)*f_roll2*(imass1 + imass3) + &
     781              :                lg4x6(iconst)%lambda(1)*imass1*f_roll1 + &
     782              :                lg4x6(iconst)%lambda(3)*imass1*f_roll3 + &
     783              :                lg4x6(iconst)%lambda(4)*imass3*f_roll4 - &
     784         8900 :                lg4x6(iconst)%lambda(6)*imass3*f_roll6
     785              :          bvec(2, 1) = g4x6_list(iconst)%dac*g4x6_list(iconst)%dac &
     786        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_13, r0_13)
     787              :          vec = lg4x6(iconst)%lambda(3)*f_roll3*(imass1 + imass4) + &
     788              :                lg4x6(iconst)%lambda(1)*imass1*f_roll1 + &
     789              :                lg4x6(iconst)%lambda(2)*imass1*f_roll2 + &
     790              :                lg4x6(iconst)%lambda(5)*imass4*f_roll5 + &
     791         8900 :                lg4x6(iconst)%lambda(6)*imass4*f_roll6
     792              :          bvec(3, 1) = g4x6_list(iconst)%dad*g4x6_list(iconst)%dad &
     793        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_14, r0_14)
     794              :          vec = lg4x6(iconst)%lambda(4)*f_roll4*(imass2 + imass3) - &
     795              :                lg4x6(iconst)%lambda(1)*imass2*f_roll1 + &
     796              :                lg4x6(iconst)%lambda(2)*imass3*f_roll2 + &
     797              :                lg4x6(iconst)%lambda(5)*imass2*f_roll5 - &
     798         8900 :                lg4x6(iconst)%lambda(6)*imass3*f_roll6
     799              :          bvec(4, 1) = g4x6_list(iconst)%dbc*g4x6_list(iconst)%dbc &
     800        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_23, r0_23)
     801              :          vec = lg4x6(iconst)%lambda(5)*f_roll5*(imass2 + imass4) - &
     802              :                lg4x6(iconst)%lambda(1)*imass2*f_roll1 + &
     803              :                lg4x6(iconst)%lambda(3)*imass4*f_roll3 + &
     804              :                lg4x6(iconst)%lambda(4)*imass2*f_roll4 + &
     805         8900 :                lg4x6(iconst)%lambda(6)*imass4*f_roll6
     806              :          bvec(5, 1) = g4x6_list(iconst)%dbd*g4x6_list(iconst)%dbd &
     807        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_24, r0_24)
     808              :          vec = lg4x6(iconst)%lambda(6)*f_roll6*(imass3 + imass4) - &
     809              :                lg4x6(iconst)%lambda(2)*imass3*f_roll2 + &
     810              :                lg4x6(iconst)%lambda(3)*imass4*f_roll3 - &
     811              :                lg4x6(iconst)%lambda(4)*imass3*f_roll4 + &
     812         8900 :                lg4x6(iconst)%lambda(5)*imass4*f_roll5
     813              :          bvec(6, 1) = g4x6_list(iconst)%dcd*g4x6_list(iconst)%dcd &
     814        17800 :                       - dtsqby2*dtsqby2*DOT_PRODUCT(vec, vec) - DOT_PRODUCT(r0_34, r0_34)
     815        17800 :          bvec = bvec*idtsq
     816              : 
     817              :          ! get lambda
     818         2225 :          atemp = amat
     819         2225 :          CALL solve_system(atemp, 6, bvec)
     820        15575 :          lg4x6(iconst)%lambda(:) = bvec(:, 1)
     821              :          lg4x6(iconst)%del_lambda(:) = lg4x6(iconst)%lambda(:) - &
     822        15575 :                                        lg4x6(iconst)%lambda_old(:)
     823        15575 :          lg4x6(iconst)%lambda_old(:) = lg4x6(iconst)%lambda(:)
     824              : 
     825              :          fc1 = lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
     826              :                lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb + &
     827         8900 :                lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc
     828              :          fc2 = -lg4x6(iconst)%del_lambda(1)*lg4x6(iconst)%fa + &
     829              :                lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
     830         8900 :                lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe
     831              :          fc3 = -lg4x6(iconst)%del_lambda(2)*lg4x6(iconst)%fb - &
     832              :                lg4x6(iconst)%del_lambda(4)*lg4x6(iconst)%fd + &
     833         8900 :                lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
     834              :          fc4 = -lg4x6(iconst)%del_lambda(3)*lg4x6(iconst)%fc - &
     835              :                lg4x6(iconst)%del_lambda(5)*lg4x6(iconst)%fe - &
     836         8900 :                lg4x6(iconst)%del_lambda(6)*lg4x6(iconst)%ff
     837        35600 :          r1(:) = pos(:, index_a) + imass1*dtsqby2*MATMUL(r_shake, fc1)
     838        35600 :          r2(:) = pos(:, index_b) + imass2*dtsqby2*MATMUL(r_shake, fc2)
     839        35600 :          r3(:) = pos(:, index_c) + imass3*dtsqby2*MATMUL(r_shake, fc3)
     840        35600 :          r4(:) = pos(:, index_d) + imass4*dtsqby2*MATMUL(r_shake, fc4)
     841        35600 :          v1(:) = vel(:, index_a) + imass1*dtby2*MATMUL(r_shake, fc1)
     842        35600 :          v2(:) = vel(:, index_b) + imass2*dtby2*MATMUL(r_shake, fc2)
     843        35600 :          v3(:) = vel(:, index_c) + imass3*dtby2*MATMUL(r_shake, fc3)
     844        35600 :          v4(:) = vel(:, index_d) + imass4*dtby2*MATMUL(r_shake, fc4)
     845              : 
     846         8900 :          r12 = r1 - r2
     847         8900 :          r13 = r1 - r3
     848         8900 :          r23 = r2 - r3
     849         8900 :          r14 = r1 - r4
     850         8900 :          r24 = r2 - r4
     851         8900 :          r34 = r3 - r4
     852              :          ! compute the tolerance:
     853              :          sigma = DOT_PRODUCT(r12, r12) - g4x6_list(iconst)%dab* &
     854         8900 :                  g4x6_list(iconst)%dab
     855         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     856              :          sigma = DOT_PRODUCT(r13, r13) - g4x6_list(iconst)%dac* &
     857         8900 :                  g4x6_list(iconst)%dac
     858         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     859              :          sigma = DOT_PRODUCT(r14, r14) - g4x6_list(iconst)%dad* &
     860         8900 :                  g4x6_list(iconst)%dad
     861         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     862              :          sigma = DOT_PRODUCT(r23, r23) - g4x6_list(iconst)%dbc* &
     863         8900 :                  g4x6_list(iconst)%dbc
     864         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     865              :          sigma = DOT_PRODUCT(r24, r24) - g4x6_list(iconst)%dbd* &
     866         8900 :                  g4x6_list(iconst)%dbd
     867         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     868              :          sigma = DOT_PRODUCT(r34, r34) - g4x6_list(iconst)%dcd* &
     869         8900 :                  g4x6_list(iconst)%dcd
     870         2225 :          max_sigma = MAX(max_sigma, ABS(sigma))
     871              : 
     872              :          ! update positions with full multiplier
     873         8900 :          pos(:, index_a) = r1(:)
     874         8900 :          pos(:, index_b) = r2(:)
     875         8900 :          pos(:, index_c) = r3(:)
     876         8900 :          pos(:, index_d) = r4(:)
     877              : 
     878              :          ! update velocites with full multiplier
     879         8900 :          vel(:, index_a) = v1(:)
     880         8900 :          vel(:, index_b) = v2(:)
     881         8900 :          vel(:, index_c) = v3(:)
     882        11125 :          vel(:, index_d) = v4(:)
     883              :       END DO
     884         2225 :    END SUBROUTINE shake_roll_4x6_low
     885              : 
     886              : ! **************************************************************************************************
     887              : !> \brief ...
     888              : !> \param fixd_list ...
     889              : !> \param g4x6_list ...
     890              : !> \param lg4x6 ...
     891              : !> \param first_atom ...
     892              : !> \param particle_set ...
     893              : !> \param vel ...
     894              : !> \param dt ...
     895              : !> \par History
     896              : !>      none
     897              : ! **************************************************************************************************
     898          702 :    SUBROUTINE rattle_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
     899          702 :                              particle_set, vel, dt)
     900              : 
     901              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     902              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     903              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     904              :       INTEGER, INTENT(IN)                                :: first_atom
     905              :       TYPE(particle_type), POINTER                       :: particle_set(:)
     906              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     907              :       REAL(kind=dp), INTENT(in)                          :: dt
     908              : 
     909              :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
     910              :                                                             index_d
     911              :       REAL(KIND=dp)                                      :: dtby2, idt, imass1, imass2, imass3, &
     912              :                                                             imass4, mass
     913              :       REAL(KIND=dp), DIMENSION(3)                        :: fc1, fc2, fc3, fc4, r12, r13, r14, r23, &
     914              :                                                             r24, r34, v12, v13, v14, v23, v24, v34
     915              :       REAL(KIND=dp), DIMENSION(6, 1)                     :: bvec
     916              :       REAL(KIND=dp), DIMENSION(6, 6)                     :: amat
     917              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     918              : 
     919          702 :       idt = 1.0_dp/dt
     920          702 :       dtby2 = dt*.5_dp
     921         1404 :       DO iconst = 1, SIZE(g4x6_list)
     922          702 :          IF (g4x6_list(iconst)%restraint%active) CYCLE
     923          692 :          index_a = g4x6_list(iconst)%a + first_atom - 1
     924          692 :          index_b = g4x6_list(iconst)%b + first_atom - 1
     925          692 :          index_c = g4x6_list(iconst)%c + first_atom - 1
     926          692 :          index_d = g4x6_list(iconst)%d + first_atom - 1
     927         2768 :          v12(:) = vel(:, index_a) - vel(:, index_b)
     928         2768 :          v13(:) = vel(:, index_a) - vel(:, index_c)
     929         2768 :          v14(:) = vel(:, index_a) - vel(:, index_d)
     930         2768 :          v23(:) = vel(:, index_b) - vel(:, index_c)
     931         2768 :          v24(:) = vel(:, index_b) - vel(:, index_d)
     932         2768 :          v34(:) = vel(:, index_c) - vel(:, index_d)
     933              : 
     934         2768 :          r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
     935         2768 :          r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
     936         2768 :          r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
     937         2768 :          r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
     938         2768 :          r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
     939         2768 :          r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
     940          692 :          atomic_kind => particle_set(index_a)%atomic_kind
     941          692 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     942          692 :          imass1 = 1.0_dp/mass
     943          692 :          atomic_kind => particle_set(index_b)%atomic_kind
     944          692 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     945          692 :          imass2 = 1.0_dp/mass
     946          692 :          atomic_kind => particle_set(index_c)%atomic_kind
     947          692 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     948          692 :          imass3 = 1.0_dp/mass
     949          692 :          atomic_kind => particle_set(index_d)%atomic_kind
     950          692 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     951          692 :          imass4 = 1.0_dp/mass
     952         2768 :          lg4x6(iconst)%fa = -2.0_dp*r12
     953         2768 :          lg4x6(iconst)%fb = -2.0_dp*r13
     954         2768 :          lg4x6(iconst)%fc = -2.0_dp*r14
     955         2768 :          lg4x6(iconst)%fd = -2.0_dp*r23
     956         2768 :          lg4x6(iconst)%fe = -2.0_dp*r24
     957         2768 :          lg4x6(iconst)%ff = -2.0_dp*r34
     958              :          ! Check for fixed atom constraints
     959              :          CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
     960          692 :                                         index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
     961              :          ! construct matrix
     962         2768 :          amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, lg4x6(iconst)%fa)
     963         2768 :          amat(1, 2) = imass1*DOT_PRODUCT(r12, lg4x6(iconst)%fb)
     964         2768 :          amat(1, 3) = imass1*DOT_PRODUCT(r12, lg4x6(iconst)%fc)
     965         2768 :          amat(1, 4) = -imass2*DOT_PRODUCT(r12, lg4x6(iconst)%fd)
     966         2768 :          amat(1, 5) = -imass2*DOT_PRODUCT(r12, lg4x6(iconst)%fe)
     967          692 :          amat(1, 6) = 0.0_dp
     968              : 
     969         2768 :          amat(2, 1) = imass1*DOT_PRODUCT(r13, lg4x6(iconst)%fa)
     970         2768 :          amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, lg4x6(iconst)%fb)
     971         2768 :          amat(2, 3) = imass1*DOT_PRODUCT(r13, lg4x6(iconst)%fc)
     972         2768 :          amat(2, 4) = imass3*DOT_PRODUCT(r13, lg4x6(iconst)%fd)
     973          692 :          amat(2, 5) = 0.0_dp
     974         2768 :          amat(2, 6) = -imass3*DOT_PRODUCT(r13, lg4x6(iconst)%ff)
     975              : 
     976         2768 :          amat(3, 1) = imass1*DOT_PRODUCT(r14, lg4x6(iconst)%fa)
     977         2768 :          amat(3, 2) = imass1*DOT_PRODUCT(r14, lg4x6(iconst)%fb)
     978         2768 :          amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r14, lg4x6(iconst)%fc)
     979          692 :          amat(3, 4) = 0.0_dp
     980         2768 :          amat(3, 5) = imass4*DOT_PRODUCT(r14, lg4x6(iconst)%fe)
     981         2768 :          amat(3, 6) = imass4*DOT_PRODUCT(r14, lg4x6(iconst)%ff)
     982              : 
     983         2768 :          amat(4, 1) = -imass2*DOT_PRODUCT(r23, lg4x6(iconst)%fa)
     984         2768 :          amat(4, 2) = imass3*DOT_PRODUCT(r23, lg4x6(iconst)%fb)
     985          692 :          amat(4, 3) = 0.0_dp
     986         2768 :          amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r23, lg4x6(iconst)%fd)
     987         2768 :          amat(4, 5) = imass2*DOT_PRODUCT(r23, lg4x6(iconst)%fe)
     988         2768 :          amat(4, 6) = -imass3*DOT_PRODUCT(r23, lg4x6(iconst)%ff)
     989              : 
     990         2768 :          amat(5, 1) = -imass2*DOT_PRODUCT(r24, lg4x6(iconst)%fa)
     991          692 :          amat(5, 2) = 0.0_dp
     992         2768 :          amat(5, 3) = imass4*DOT_PRODUCT(r24, lg4x6(iconst)%fc)
     993         2768 :          amat(5, 4) = imass2*DOT_PRODUCT(r24, lg4x6(iconst)%fd)
     994         2768 :          amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r24, lg4x6(iconst)%fe)
     995         2768 :          amat(5, 6) = imass4*DOT_PRODUCT(r24, lg4x6(iconst)%ff)
     996              : 
     997          692 :          amat(6, 1) = 0.0_dp
     998         2768 :          amat(6, 2) = -imass3*DOT_PRODUCT(r34, lg4x6(iconst)%fb)
     999         2768 :          amat(6, 3) = imass4*DOT_PRODUCT(r34, lg4x6(iconst)%fc)
    1000         2768 :          amat(6, 4) = -imass3*DOT_PRODUCT(r34, lg4x6(iconst)%fd)
    1001         2768 :          amat(6, 5) = imass4*DOT_PRODUCT(r34, lg4x6(iconst)%fe)
    1002         2768 :          amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r34, lg4x6(iconst)%ff)
    1003              : 
    1004              :          ! construct solution vector
    1005         2768 :          bvec(1, 1) = DOT_PRODUCT(r12, v12)
    1006         2768 :          bvec(2, 1) = DOT_PRODUCT(r13, v13)
    1007         2768 :          bvec(3, 1) = DOT_PRODUCT(r14, v14)
    1008         2768 :          bvec(4, 1) = DOT_PRODUCT(r23, v23)
    1009         2768 :          bvec(5, 1) = DOT_PRODUCT(r24, v24)
    1010         2768 :          bvec(6, 1) = DOT_PRODUCT(r34, v34)
    1011         5536 :          bvec = -bvec*2.0_dp*idt
    1012              : 
    1013              :          ! get lambda
    1014          692 :          CALL solve_system(amat, 6, bvec)
    1015         4844 :          lg4x6(iconst)%lambda(:) = bvec(:, 1)
    1016              : 
    1017              :          fc1 = lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + &
    1018              :                lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb + &
    1019         2768 :                lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc
    1020              :          fc2 = -lg4x6(iconst)%lambda(1)*lg4x6(iconst)%fa + &
    1021              :                lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + &
    1022         2768 :                lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe
    1023              :          fc3 = -lg4x6(iconst)%lambda(2)*lg4x6(iconst)%fb - &
    1024              :                lg4x6(iconst)%lambda(4)*lg4x6(iconst)%fd + &
    1025         2768 :                lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff
    1026              :          fc4 = -lg4x6(iconst)%lambda(3)*lg4x6(iconst)%fc - &
    1027              :                lg4x6(iconst)%lambda(5)*lg4x6(iconst)%fe - &
    1028         2768 :                lg4x6(iconst)%lambda(6)*lg4x6(iconst)%ff
    1029         2768 :          vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
    1030         2768 :          vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
    1031         2768 :          vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
    1032         4172 :          vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
    1033              :       END DO
    1034          702 :    END SUBROUTINE rattle_4x6_low
    1035              : 
    1036              : ! **************************************************************************************************
    1037              : !> \brief ...
    1038              : !> \param fixd_list ...
    1039              : !> \param g4x6_list ...
    1040              : !> \param lg4x6 ...
    1041              : !> \param first_atom ...
    1042              : !> \param particle_set ...
    1043              : !> \param vel ...
    1044              : !> \param r_rattle ...
    1045              : !> \param dt ...
    1046              : !> \param veps ...
    1047              : !> \par History
    1048              : !>      none
    1049              : ! **************************************************************************************************
    1050         1024 :    SUBROUTINE rattle_roll_4x6_low(fixd_list, g4x6_list, lg4x6, first_atom, &
    1051         1024 :                                   particle_set, vel, r_rattle, dt, veps)
    1052              : 
    1053              :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
    1054              :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
    1055              :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
    1056              :       INTEGER, INTENT(IN)                                :: first_atom
    1057              :       TYPE(particle_type), POINTER                       :: particle_set(:)
    1058              :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
    1059              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_rattle
    1060              :       REAL(kind=dp), INTENT(in)                          :: dt
    1061              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: veps
    1062              : 
    1063              :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
    1064              :                                                             index_d
    1065              :       REAL(KIND=dp)                                      :: dtby2, idt, imass1, imass2, imass3, &
    1066              :                                                             imass4, mass
    1067              :       REAL(KIND=dp), DIMENSION(3)                        :: f_roll1, f_roll2, f_roll3, f_roll4, &
    1068              :                                                             f_roll5, f_roll6, fc1, fc2, fc3, fc4, &
    1069              :                                                             r12, r13, r14, r23, r24, r34, v12, &
    1070              :                                                             v13, v14, v23, v24, v34
    1071              :       REAL(KIND=dp), DIMENSION(6)                        :: lambda
    1072              :       REAL(KIND=dp), DIMENSION(6, 1)                     :: bvec
    1073              :       REAL(KIND=dp), DIMENSION(6, 6)                     :: amat
    1074              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
    1075              : 
    1076         1024 :       idt = 1.0_dp/dt
    1077         1024 :       dtby2 = dt*.5_dp
    1078         2048 :       DO iconst = 1, SIZE(g4x6_list)
    1079         1024 :          IF (g4x6_list(iconst)%restraint%active) CYCLE
    1080         1024 :          index_a = g4x6_list(iconst)%a + first_atom - 1
    1081         1024 :          index_b = g4x6_list(iconst)%b + first_atom - 1
    1082         1024 :          index_c = g4x6_list(iconst)%c + first_atom - 1
    1083         1024 :          index_d = g4x6_list(iconst)%d + first_atom - 1
    1084         4096 :          v12(:) = vel(:, index_a) - vel(:, index_b)
    1085         4096 :          v13(:) = vel(:, index_a) - vel(:, index_c)
    1086         4096 :          v14(:) = vel(:, index_a) - vel(:, index_d)
    1087         4096 :          v23(:) = vel(:, index_b) - vel(:, index_c)
    1088         4096 :          v24(:) = vel(:, index_b) - vel(:, index_d)
    1089         4096 :          v34(:) = vel(:, index_c) - vel(:, index_d)
    1090              : 
    1091         4096 :          r12(:) = particle_set(index_a)%r(:) - particle_set(index_b)%r(:)
    1092         4096 :          r13(:) = particle_set(index_a)%r(:) - particle_set(index_c)%r(:)
    1093         4096 :          r14(:) = particle_set(index_a)%r(:) - particle_set(index_d)%r(:)
    1094         4096 :          r23(:) = particle_set(index_b)%r(:) - particle_set(index_c)%r(:)
    1095         4096 :          r24(:) = particle_set(index_b)%r(:) - particle_set(index_d)%r(:)
    1096         4096 :          r34(:) = particle_set(index_c)%r(:) - particle_set(index_d)%r(:)
    1097         1024 :          atomic_kind => particle_set(index_a)%atomic_kind
    1098         1024 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1099         1024 :          imass1 = 1.0_dp/mass
    1100         1024 :          atomic_kind => particle_set(index_b)%atomic_kind
    1101         1024 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1102         1024 :          imass2 = 1.0_dp/mass
    1103         1024 :          atomic_kind => particle_set(index_c)%atomic_kind
    1104         1024 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1105         1024 :          imass3 = 1.0_dp/mass
    1106         1024 :          atomic_kind => particle_set(index_d)%atomic_kind
    1107         1024 :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
    1108         1024 :          imass4 = 1.0_dp/mass
    1109         4096 :          lg4x6(iconst)%fa = -2.0_dp*r12
    1110         4096 :          lg4x6(iconst)%fb = -2.0_dp*r13
    1111         4096 :          lg4x6(iconst)%fc = -2.0_dp*r14
    1112         4096 :          lg4x6(iconst)%fd = -2.0_dp*r23
    1113         4096 :          lg4x6(iconst)%fe = -2.0_dp*r24
    1114         4096 :          lg4x6(iconst)%ff = -2.0_dp*r34
    1115              :          ! Check for fixed atom constraints
    1116              :          CALL check_fixed_atom_cns_g4x6(imass1, imass2, imass3, imass4, &
    1117         1024 :                                         index_a, index_b, index_c, index_d, fixd_list, lg4x6(iconst))
    1118              :          ! roll the fc
    1119        13312 :          f_roll1 = MATMUL(r_rattle, lg4x6(iconst)%fa)
    1120        13312 :          f_roll2 = MATMUL(r_rattle, lg4x6(iconst)%fb)
    1121        13312 :          f_roll3 = MATMUL(r_rattle, lg4x6(iconst)%fc)
    1122        13312 :          f_roll4 = MATMUL(r_rattle, lg4x6(iconst)%fd)
    1123        13312 :          f_roll5 = MATMUL(r_rattle, lg4x6(iconst)%fe)
    1124        13312 :          f_roll6 = MATMUL(r_rattle, lg4x6(iconst)%ff)
    1125              :          ! construct matrix
    1126         4096 :          amat(1, 1) = (imass1 + imass2)*DOT_PRODUCT(r12, f_roll1)
    1127         4096 :          amat(1, 2) = imass1*DOT_PRODUCT(r12, f_roll2)
    1128         4096 :          amat(1, 3) = imass1*DOT_PRODUCT(r12, f_roll3)
    1129         4096 :          amat(1, 4) = -imass2*DOT_PRODUCT(r12, f_roll4)
    1130         4096 :          amat(1, 5) = -imass2*DOT_PRODUCT(r12, f_roll5)
    1131         1024 :          amat(1, 6) = 0.0_dp
    1132              : 
    1133         4096 :          amat(2, 1) = imass1*DOT_PRODUCT(r13, f_roll1)
    1134         4096 :          amat(2, 2) = (imass1 + imass3)*DOT_PRODUCT(r13, f_roll2)
    1135         4096 :          amat(2, 3) = imass1*DOT_PRODUCT(r13, f_roll3)
    1136         4096 :          amat(2, 4) = imass3*DOT_PRODUCT(r13, f_roll4)
    1137         1024 :          amat(2, 5) = 0.0_dp
    1138         4096 :          amat(2, 6) = -imass3*DOT_PRODUCT(r13, f_roll6)
    1139              : 
    1140         4096 :          amat(3, 1) = imass1*DOT_PRODUCT(r14, f_roll1)
    1141         4096 :          amat(3, 2) = imass1*DOT_PRODUCT(r14, f_roll2)
    1142         4096 :          amat(3, 3) = (imass1 + imass4)*DOT_PRODUCT(r14, f_roll3)
    1143         1024 :          amat(3, 4) = 0.0_dp
    1144         4096 :          amat(3, 5) = imass4*DOT_PRODUCT(r14, f_roll5)
    1145         4096 :          amat(3, 6) = imass4*DOT_PRODUCT(r14, f_roll6)
    1146              : 
    1147         4096 :          amat(4, 1) = -imass2*DOT_PRODUCT(r23, f_roll1)
    1148         4096 :          amat(4, 2) = imass3*DOT_PRODUCT(r23, f_roll2)
    1149         1024 :          amat(4, 3) = 0.0_dp
    1150         4096 :          amat(4, 4) = (imass3 + imass2)*DOT_PRODUCT(r23, f_roll4)
    1151         4096 :          amat(4, 5) = imass2*DOT_PRODUCT(r23, f_roll5)
    1152         4096 :          amat(4, 6) = -imass3*DOT_PRODUCT(r23, f_roll6)
    1153              : 
    1154         4096 :          amat(5, 1) = -imass2*DOT_PRODUCT(r24, f_roll1)
    1155         1024 :          amat(5, 2) = 0.0_dp
    1156         4096 :          amat(5, 3) = imass4*DOT_PRODUCT(r24, f_roll3)
    1157         4096 :          amat(5, 4) = imass2*DOT_PRODUCT(r24, f_roll4)
    1158         4096 :          amat(5, 5) = (imass4 + imass2)*DOT_PRODUCT(r24, f_roll5)
    1159         4096 :          amat(5, 6) = imass4*DOT_PRODUCT(r24, f_roll6)
    1160              : 
    1161         1024 :          amat(6, 1) = 0.0_dp
    1162         4096 :          amat(6, 2) = -imass3*DOT_PRODUCT(r34, f_roll2)
    1163         4096 :          amat(6, 3) = imass4*DOT_PRODUCT(r34, f_roll3)
    1164         4096 :          amat(6, 4) = -imass3*DOT_PRODUCT(r34, f_roll4)
    1165         4096 :          amat(6, 5) = imass4*DOT_PRODUCT(r34, f_roll5)
    1166         4096 :          amat(6, 6) = (imass3 + imass4)*DOT_PRODUCT(r34, f_roll6)
    1167              : 
    1168              :          ! construct solution vector
    1169        22528 :          bvec(1, 1) = DOT_PRODUCT(r12, v12 + MATMUL(veps, r12))
    1170        22528 :          bvec(2, 1) = DOT_PRODUCT(r13, v13 + MATMUL(veps, r13))
    1171        22528 :          bvec(3, 1) = DOT_PRODUCT(r14, v14 + MATMUL(veps, r14))
    1172        22528 :          bvec(4, 1) = DOT_PRODUCT(r23, v23 + MATMUL(veps, r23))
    1173        22528 :          bvec(5, 1) = DOT_PRODUCT(r24, v24 + MATMUL(veps, r24))
    1174        22528 :          bvec(6, 1) = DOT_PRODUCT(r34, v34 + MATMUL(veps, r34))
    1175         8192 :          bvec = -bvec*2.0_dp*idt
    1176              : 
    1177              :          ! get lambda
    1178         1024 :          CALL solve_system(amat, 6, bvec)
    1179         1024 :          lambda(:) = bvec(:, 1)
    1180         7168 :          lg4x6(iconst)%lambda(:) = lambda
    1181              : 
    1182              :          fc1 = lambda(1)*f_roll1 + &
    1183              :                lambda(2)*f_roll2 + &
    1184         4096 :                lambda(3)*f_roll3
    1185              :          fc2 = -lambda(1)*f_roll1 + &
    1186              :                lambda(4)*f_roll4 + &
    1187         4096 :                lambda(5)*f_roll5
    1188              :          fc3 = -lambda(2)*f_roll2 - &
    1189              :                lambda(4)*f_roll4 + &
    1190         4096 :                lambda(6)*f_roll6
    1191              :          fc4 = -lambda(3)*f_roll3 - &
    1192              :                lambda(5)*f_roll5 - &
    1193         4096 :                lambda(6)*f_roll6
    1194         4096 :          vel(:, index_a) = vel(:, index_a) + imass1*dtby2*fc1(:)
    1195         4096 :          vel(:, index_b) = vel(:, index_b) + imass2*dtby2*fc2(:)
    1196         4096 :          vel(:, index_c) = vel(:, index_c) + imass3*dtby2*fc3(:)
    1197         6144 :          vel(:, index_d) = vel(:, index_d) + imass4*dtby2*fc4(:)
    1198              :       END DO
    1199         1024 :    END SUBROUTINE rattle_roll_4x6_low
    1200              : 
    1201         6144 : END MODULE constraint_4x6
        

Generated by: LCOV version 2.0-1