LCOV - code coverage report
Current view: top level - src - constraint_util.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:c7c47bf) Lines: 264 266 99.2 %
Date: 2024-12-29 23:06:28 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Contains routines useful for the application of constraints during MD
      10             : !> \par History
      11             : !>      none
      12             : ! **************************************************************************************************
      13             : MODULE constraint_util
      14             :    USE cell_types,                      ONLY: cell_type
      15             :    USE colvar_methods,                  ONLY: colvar_eval_mol_f
      16             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      17             :    USE kinds,                           ONLY: dp
      18             :    USE message_passing,                 ONLY: mp_comm_type
      19             :    USE molecule_kind_types,             ONLY: colvar_constraint_type,&
      20             :                                               fixd_constraint_type,&
      21             :                                               g3x3_constraint_type,&
      22             :                                               g4x6_constraint_type,&
      23             :                                               get_molecule_kind,&
      24             :                                               molecule_kind_type
      25             :    USE molecule_types,                  ONLY: get_molecule,&
      26             :                                               global_constraint_type,&
      27             :                                               local_colvar_constraint_type,&
      28             :                                               local_constraint_type,&
      29             :                                               local_g3x3_constraint_type,&
      30             :                                               local_g4x6_constraint_type,&
      31             :                                               molecule_type
      32             :    USE particle_types,                  ONLY: particle_type
      33             :    USE virial_types,                    ONLY: virial_type
      34             : #include "./base/base_uses.f90"
      35             : 
      36             :    IMPLICIT NONE
      37             : 
      38             :    PRIVATE
      39             :    PUBLIC :: getold, &
      40             :              pv_constraint, &
      41             :              check_tol, &
      42             :              get_roll_matrix, &
      43             :              restore_temporary_set, &
      44             :              update_temporary_set
      45             : 
      46             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_util'
      47             : 
      48             : CONTAINS
      49             : 
      50             : ! **************************************************************************************************
      51             : !> \brief saves all of the old variables
      52             : !> \param gci ...
      53             : !> \param local_molecules ...
      54             : !> \param molecule_set ...
      55             : !> \param molecule_kind_set ...
      56             : !> \param particle_set ...
      57             : !> \param cell ...
      58             : !> \par History
      59             : !>      none
      60             : ! **************************************************************************************************
      61       16934 :    SUBROUTINE getold(gci, local_molecules, molecule_set, molecule_kind_set, &
      62             :                      particle_set, cell)
      63             : 
      64             :       TYPE(global_constraint_type), POINTER              :: gci
      65             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
      66             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
      67             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
      68             :       TYPE(particle_type), POINTER                       :: particle_set(:)
      69             :       TYPE(cell_type), POINTER                           :: cell
      70             : 
      71             :       INTEGER                                            :: first_atom, i, ikind, imol, n3x3con, &
      72             :                                                             n4x6con, nkind, nmol_per_kind
      73       16934 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
      74       16934 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
      75       16934 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
      76       16934 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
      77       16934 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
      78             :       TYPE(local_constraint_type), POINTER               :: lci
      79       16934 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
      80       16934 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
      81             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      82             :       TYPE(molecule_type), POINTER                       :: molecule
      83             : 
      84       16934 :       NULLIFY (fixd_list)
      85       16934 :       nkind = SIZE(molecule_kind_set)
      86             :       ! Intramolecular constraints
      87       84474 :       MOL: DO ikind = 1, nkind
      88       67540 :          nmol_per_kind = local_molecules%n_el(ikind)
      89      396005 :          DO imol = 1, nmol_per_kind
      90      311531 :             i = local_molecules%list(ikind)%array(imol)
      91      311531 :             molecule => molecule_set(i)
      92      311531 :             molecule_kind => molecule%molecule_kind
      93             :             CALL get_molecule_kind(molecule_kind, ng3x3=n3x3con, ng4x6=n4x6con, &
      94             :                                    colv_list=colv_list, g3x3_list=g3x3_list, g4x6_list=g4x6_list, &
      95      311531 :                                    fixd_list=fixd_list)
      96      311531 :             CALL get_molecule(molecule, lci=lci)
      97      311531 :             IF (.NOT. ASSOCIATED(lci)) CYCLE
      98             :             CALL get_molecule(molecule, first_atom=first_atom, &
      99      311531 :                               lcolv=lcolv, lg3x3=lg3x3, lg4x6=lg4x6)
     100             :             CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
     101      690602 :                             lcolv, lg3x3, lg4x6, first_atom, particle_set, cell)
     102             :          END DO
     103             :       END DO MOL
     104             :       ! Intermolecular constraints
     105       16934 :       IF (gci%ntot > 0) THEN
     106        1150 :          n3x3con = gci%ng3x3
     107        1150 :          n4x6con = gci%ng4x6
     108        1150 :          colv_list => gci%colv_list
     109        1150 :          g3x3_list => gci%g3x3_list
     110        1150 :          g4x6_list => gci%g4x6_list
     111        1150 :          fixd_list => gci%fixd_list
     112        1150 :          lcolv => gci%lcolv
     113        1150 :          lg3x3 => gci%lg3x3
     114        1150 :          lg4x6 => gci%lg4x6
     115             :          CALL getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
     116        1150 :                          lcolv, lg3x3, lg4x6, 1, particle_set, cell)
     117             :       END IF
     118       16934 :    END SUBROUTINE getold
     119             : 
     120             : ! **************************************************************************************************
     121             : !> \brief saves all of the old variables - Low Level
     122             : !> \param n3x3con ...
     123             : !> \param n4x6con ...
     124             : !> \param colv_list ...
     125             : !> \param g3x3_list ...
     126             : !> \param g4x6_list ...
     127             : !> \param fixd_list ...
     128             : !> \param lcolv ...
     129             : !> \param lg3x3 ...
     130             : !> \param lg4x6 ...
     131             : !> \param first_atom ...
     132             : !> \param particle_set ...
     133             : !> \param cell ...
     134             : !> \par History
     135             : !>      none
     136             : ! **************************************************************************************************
     137      312681 :    SUBROUTINE getold_low(n3x3con, n4x6con, colv_list, g3x3_list, g4x6_list, fixd_list, &
     138             :                          lcolv, lg3x3, lg4x6, first_atom, particle_set, cell)
     139             : 
     140             :       INTEGER, INTENT(IN)                                :: n3x3con, n4x6con
     141             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     142             :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     143             :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     144             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     145             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     146             :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     147             :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     148             :       INTEGER, INTENT(IN)                                :: first_atom
     149             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     150             :       TYPE(cell_type), POINTER                           :: cell
     151             : 
     152             :       INTEGER                                            :: iconst, index
     153             : 
     154      312681 :       IF (ASSOCIATED(colv_list)) THEN
     155             :          ! Collective constraints
     156       43227 :          DO iconst = 1, SIZE(colv_list)
     157             :             CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
     158       43227 :                                    particles=particle_set, fixd_list=fixd_list)
     159             :          END DO
     160             :       END IF
     161             :       ! 3x3 constraints
     162      494908 :       DO iconst = 1, n3x3con
     163      182227 :          index = g3x3_list(iconst)%a + first_atom - 1
     164     1457816 :          lg3x3(iconst)%ra_old = particle_set(index)%r
     165      182227 :          index = g3x3_list(iconst)%b + first_atom - 1
     166     1457816 :          lg3x3(iconst)%rb_old = particle_set(index)%r
     167      182227 :          index = g3x3_list(iconst)%c + first_atom - 1
     168     1770497 :          lg3x3(iconst)%rc_old = particle_set(index)%r
     169             :       END DO
     170             :       ! 4x6 constraints
     171      313708 :       DO iconst = 1, n4x6con
     172        1027 :          index = g4x6_list(iconst)%a + first_atom - 1
     173        8216 :          lg4x6(iconst)%ra_old = particle_set(index)%r
     174        1027 :          index = g4x6_list(iconst)%b + first_atom - 1
     175        8216 :          lg4x6(iconst)%rb_old = particle_set(index)%r
     176        1027 :          index = g4x6_list(iconst)%c + first_atom - 1
     177        8216 :          lg4x6(iconst)%rc_old = particle_set(index)%r
     178        1027 :          index = g4x6_list(iconst)%d + first_atom - 1
     179      320897 :          lg4x6(iconst)%rd_old = particle_set(index)%r
     180             :       END DO
     181             : 
     182      312681 :    END SUBROUTINE getold_low
     183             : 
     184             : ! **************************************************************************************************
     185             : !> \brief ...
     186             : !> \param gci ...
     187             : !> \param local_molecules ...
     188             : !> \param molecule_set ...
     189             : !> \param molecule_kind_set ...
     190             : !> \param particle_set ...
     191             : !> \param virial ...
     192             : !> \param group ...
     193             : !> \par History
     194             : !>      none
     195             : ! **************************************************************************************************
     196       21178 :    SUBROUTINE pv_constraint(gci, local_molecules, molecule_set, molecule_kind_set, &
     197             :                             particle_set, virial, group)
     198             : 
     199             :       TYPE(global_constraint_type), POINTER              :: gci
     200             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     201             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     202             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     203             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     204             :       TYPE(virial_type), INTENT(INOUT)                   :: virial
     205             : 
     206             :       CLASS(mp_comm_type), INTENT(IN)                     :: group
     207             : 
     208             :       INTEGER                                            :: first_atom, i, ikind, imol, ng3x3, &
     209             :                                                             ng4x6, nkind, nmol_per_kind
     210             :       REAL(KIND=dp)                                      :: pv(3, 3)
     211       21178 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     212       21178 :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     213       21178 :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     214       21178 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     215             :       TYPE(local_constraint_type), POINTER               :: lci
     216       21178 :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     217       21178 :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     218             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     219             :       TYPE(molecule_type), POINTER                       :: molecule
     220             : 
     221       21178 :       pv = 0.0_dp
     222       21178 :       nkind = SIZE(molecule_kind_set)
     223             :       ! Intramolecular Constraints
     224       92988 :       MOL: DO ikind = 1, nkind
     225       71810 :          nmol_per_kind = local_molecules%n_el(ikind)
     226      660848 :          DO imol = 1, nmol_per_kind
     227      567860 :             i = local_molecules%list(ikind)%array(imol)
     228      567860 :             molecule => molecule_set(i)
     229      567860 :             molecule_kind => molecule%molecule_kind
     230             :             CALL get_molecule_kind(molecule_kind, ng3x3=ng3x3, &
     231             :                                    ng4x6=ng4x6, g3x3_list=g3x3_list, g4x6_list=g4x6_list, &
     232      567860 :                                    colv_list=colv_list)
     233      567860 :             CALL get_molecule(molecule, lci=lci)
     234      567860 :             IF (.NOT. ASSOCIATED(lci)) CYCLE
     235             :             CALL get_molecule(molecule, first_atom=first_atom, lg3x3=lg3x3, &
     236      567860 :                               lg4x6=lg4x6, lcolv=lcolv)
     237             :             CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
     238     1207530 :                                    first_atom, lg3x3, lg4x6, lcolv, particle_set, pv)
     239             :          END DO
     240             :       END DO MOL
     241             :       ! Intermolecular constraints
     242       21178 :       IF (gci%ntot > 0) THEN
     243        1120 :          ng3x3 = gci%ng3x3
     244        1120 :          ng4x6 = gci%ng4x6
     245        1120 :          colv_list => gci%colv_list
     246        1120 :          g3x3_list => gci%g3x3_list
     247        1120 :          g4x6_list => gci%g4x6_list
     248        1120 :          lcolv => gci%lcolv
     249        1120 :          lg3x3 => gci%lg3x3
     250        1120 :          lg4x6 => gci%lg4x6
     251             :          CALL pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
     252        1120 :                                 1, lg3x3, lg4x6, lcolv, particle_set, pv)
     253             :       END IF
     254       21178 :       CALL group%sum(pv)
     255             :       ! Symmetrize PV
     256       21178 :       pv(1, 2) = 0.5_dp*(pv(1, 2) + pv(2, 1))
     257       21178 :       pv(2, 1) = pv(1, 2)
     258       21178 :       pv(1, 3) = 0.5_dp*(pv(1, 3) + pv(3, 1))
     259       21178 :       pv(3, 1) = pv(1, 3)
     260       21178 :       pv(3, 2) = 0.5_dp*(pv(3, 2) + pv(2, 3))
     261       21178 :       pv(2, 3) = pv(3, 2)
     262             :       ! Store in virial type
     263      275314 :       virial%pv_constraint = pv
     264             : 
     265       21178 :    END SUBROUTINE pv_constraint
     266             : 
     267             : ! **************************************************************************************************
     268             : !> \brief ...
     269             : !> \param ng3x3 ...
     270             : !> \param ng4x6 ...
     271             : !> \param g3x3_list ...
     272             : !> \param g4x6_list ...
     273             : !> \param colv_list ...
     274             : !> \param first_atom ...
     275             : !> \param lg3x3 ...
     276             : !> \param lg4x6 ...
     277             : !> \param lcolv ...
     278             : !> \param particle_set ...
     279             : !> \param pv ...
     280             : !> \par History
     281             : !>      none
     282             : ! **************************************************************************************************
     283      568980 :    SUBROUTINE pv_constraint_low(ng3x3, ng4x6, g3x3_list, g4x6_list, colv_list, &
     284             :                                 first_atom, lg3x3, lg4x6, lcolv, particle_set, pv)
     285             : 
     286             :       INTEGER, INTENT(IN)                                :: ng3x3, ng4x6
     287             :       TYPE(g3x3_constraint_type), POINTER                :: g3x3_list(:)
     288             :       TYPE(g4x6_constraint_type), POINTER                :: g4x6_list(:)
     289             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     290             :       INTEGER, INTENT(IN)                                :: first_atom
     291             :       TYPE(local_g3x3_constraint_type), POINTER          :: lg3x3(:)
     292             :       TYPE(local_g4x6_constraint_type), POINTER          :: lg4x6(:)
     293             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     294             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     295             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: pv
     296             : 
     297             :       INTEGER                                            :: iconst, index_a, index_b, index_c, &
     298             :                                                             index_d
     299             :       REAL(KIND=dp)                                      :: fc1(3), fc2(3), fc3(3), fc4(3), &
     300             :                                                             lambda_3x3(3), lambda_4x6(6)
     301             : 
     302      568980 :       IF (ASSOCIATED(colv_list)) THEN
     303             :          ! Colvar Constraints
     304       62284 :          DO iconst = 1, SIZE(colv_list)
     305       62284 :             CALL pv_colv_eval(pv, lcolv(iconst), particle_set)
     306             :          END DO
     307             :       END IF
     308             :       ! 3x3
     309     1000052 :       DO iconst = 1, ng3x3
     310             :          !  pv gets updated with FULL multiplier
     311     1724288 :          lambda_3x3 = lg3x3(iconst)%lambda
     312             : 
     313             :          fc1 = lambda_3x3(1)*lg3x3(iconst)%fa + &
     314     1724288 :                lambda_3x3(2)*lg3x3(iconst)%fb
     315             :          fc2 = -lambda_3x3(1)*lg3x3(iconst)%fa + &
     316     1724288 :                lambda_3x3(3)*lg3x3(iconst)%fc
     317             :          fc3 = -lambda_3x3(2)*lg3x3(iconst)%fb - &
     318     1724288 :                lambda_3x3(3)*lg3x3(iconst)%fc
     319      431072 :          index_a = g3x3_list(iconst)%a + first_atom - 1
     320      431072 :          index_b = g3x3_list(iconst)%b + first_atom - 1
     321      431072 :          index_c = g3x3_list(iconst)%c + first_atom - 1
     322             : 
     323             :          !pv(1,1)
     324      431072 :          pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
     325      431072 :          pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
     326      431072 :          pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
     327             :          !pv(1,2)
     328      431072 :          pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
     329      431072 :          pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
     330      431072 :          pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
     331             :          !pv(1,3)
     332      431072 :          pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
     333      431072 :          pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
     334      431072 :          pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
     335             :          !pv(2,1)
     336      431072 :          pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
     337      431072 :          pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
     338      431072 :          pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
     339             :          !pv(2,2)
     340      431072 :          pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
     341      431072 :          pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
     342      431072 :          pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
     343             :          !pv(2,3)
     344      431072 :          pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
     345      431072 :          pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
     346      431072 :          pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
     347             :          !pv(3,1)
     348      431072 :          pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
     349      431072 :          pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
     350      431072 :          pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
     351             :          !pv(3,2)
     352      431072 :          pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
     353      431072 :          pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
     354      431072 :          pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
     355             :          !pv(3,3)
     356      431072 :          pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
     357      431072 :          pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
     358     1000052 :          pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
     359             :       END DO
     360             : 
     361             :       ! 4x6
     362      572437 :       DO iconst = 1, ng4x6
     363             :          !  pv gets updated with FULL multiplier
     364       24199 :          lambda_4x6 = lg4x6(iconst)%lambda
     365             : 
     366             :          fc1 = lambda_4x6(1)*lg4x6(iconst)%fa + &
     367             :                lambda_4x6(2)*lg4x6(iconst)%fb + &
     368       13828 :                lambda_4x6(3)*lg4x6(iconst)%fc
     369             :          fc2 = -lambda_4x6(1)*lg4x6(iconst)%fa + &
     370             :                lambda_4x6(4)*lg4x6(iconst)%fd + &
     371       13828 :                lambda_4x6(5)*lg4x6(iconst)%fe
     372             :          fc3 = -lambda_4x6(2)*lg4x6(iconst)%fb - &
     373             :                lambda_4x6(4)*lg4x6(iconst)%fd + &
     374       13828 :                lambda_4x6(6)*lg4x6(iconst)%ff
     375             :          fc4 = -lambda_4x6(3)*lg4x6(iconst)%fc - &
     376             :                lambda_4x6(5)*lg4x6(iconst)%fe - &
     377       13828 :                lambda_4x6(6)*lg4x6(iconst)%ff
     378        3457 :          index_a = g4x6_list(iconst)%a + first_atom - 1
     379        3457 :          index_b = g4x6_list(iconst)%b + first_atom - 1
     380        3457 :          index_c = g4x6_list(iconst)%c + first_atom - 1
     381        3457 :          index_d = g4x6_list(iconst)%d + first_atom - 1
     382             : 
     383             :          !pv(1,1)
     384        3457 :          pv(1, 1) = pv(1, 1) + fc1(1)*particle_set(index_a)%r(1)
     385        3457 :          pv(1, 1) = pv(1, 1) + fc2(1)*particle_set(index_b)%r(1)
     386        3457 :          pv(1, 1) = pv(1, 1) + fc3(1)*particle_set(index_c)%r(1)
     387        3457 :          pv(1, 1) = pv(1, 1) + fc4(1)*particle_set(index_d)%r(1)
     388             :          !pv(1,2)
     389        3457 :          pv(1, 2) = pv(1, 2) + fc1(1)*particle_set(index_a)%r(2)
     390        3457 :          pv(1, 2) = pv(1, 2) + fc2(1)*particle_set(index_b)%r(2)
     391        3457 :          pv(1, 2) = pv(1, 2) + fc3(1)*particle_set(index_c)%r(2)
     392        3457 :          pv(1, 2) = pv(1, 2) + fc4(1)*particle_set(index_d)%r(2)
     393             :          !pv(1,3)
     394        3457 :          pv(1, 3) = pv(1, 3) + fc1(1)*particle_set(index_a)%r(3)
     395        3457 :          pv(1, 3) = pv(1, 3) + fc2(1)*particle_set(index_b)%r(3)
     396        3457 :          pv(1, 3) = pv(1, 3) + fc3(1)*particle_set(index_c)%r(3)
     397        3457 :          pv(1, 3) = pv(1, 3) + fc4(1)*particle_set(index_d)%r(3)
     398             :          !pv(2,1)
     399        3457 :          pv(2, 1) = pv(2, 1) + fc1(2)*particle_set(index_a)%r(1)
     400        3457 :          pv(2, 1) = pv(2, 1) + fc2(2)*particle_set(index_b)%r(1)
     401        3457 :          pv(2, 1) = pv(2, 1) + fc3(2)*particle_set(index_c)%r(1)
     402        3457 :          pv(2, 1) = pv(2, 1) + fc4(2)*particle_set(index_d)%r(1)
     403             :          !pv(2,2)
     404        3457 :          pv(2, 2) = pv(2, 2) + fc1(2)*particle_set(index_a)%r(2)
     405        3457 :          pv(2, 2) = pv(2, 2) + fc2(2)*particle_set(index_b)%r(2)
     406        3457 :          pv(2, 2) = pv(2, 2) + fc3(2)*particle_set(index_c)%r(2)
     407        3457 :          pv(2, 2) = pv(2, 2) + fc4(2)*particle_set(index_d)%r(2)
     408             :          !pv(2,3)
     409        3457 :          pv(2, 3) = pv(2, 3) + fc1(2)*particle_set(index_a)%r(3)
     410        3457 :          pv(2, 3) = pv(2, 3) + fc2(2)*particle_set(index_b)%r(3)
     411        3457 :          pv(2, 3) = pv(2, 3) + fc3(2)*particle_set(index_c)%r(3)
     412        3457 :          pv(2, 3) = pv(2, 3) + fc4(2)*particle_set(index_d)%r(3)
     413             :          !pv(3,1)
     414        3457 :          pv(3, 1) = pv(3, 1) + fc1(3)*particle_set(index_a)%r(1)
     415        3457 :          pv(3, 1) = pv(3, 1) + fc2(3)*particle_set(index_b)%r(1)
     416        3457 :          pv(3, 1) = pv(3, 1) + fc3(3)*particle_set(index_c)%r(1)
     417        3457 :          pv(3, 1) = pv(3, 1) + fc4(3)*particle_set(index_d)%r(1)
     418             :          !pv(3,2)
     419        3457 :          pv(3, 2) = pv(3, 2) + fc1(3)*particle_set(index_a)%r(2)
     420        3457 :          pv(3, 2) = pv(3, 2) + fc2(3)*particle_set(index_b)%r(2)
     421        3457 :          pv(3, 2) = pv(3, 2) + fc3(3)*particle_set(index_c)%r(2)
     422        3457 :          pv(3, 2) = pv(3, 2) + fc4(3)*particle_set(index_d)%r(2)
     423             :          !pv(3,3)
     424        3457 :          pv(3, 3) = pv(3, 3) + fc1(3)*particle_set(index_a)%r(3)
     425        3457 :          pv(3, 3) = pv(3, 3) + fc2(3)*particle_set(index_b)%r(3)
     426        3457 :          pv(3, 3) = pv(3, 3) + fc3(3)*particle_set(index_c)%r(3)
     427      572437 :          pv(3, 3) = pv(3, 3) + fc4(3)*particle_set(index_d)%r(3)
     428             :       END DO
     429             : 
     430      568980 :    END SUBROUTINE pv_constraint_low
     431             : 
     432             : ! **************************************************************************************************
     433             : !> \brief ...
     434             : !> \param pv ...
     435             : !> \param lcolv ...
     436             : !> \param particle_set ...
     437             : !> \par History
     438             : !>      none
     439             : ! **************************************************************************************************
     440       25701 :    SUBROUTINE pv_colv_eval(pv, lcolv, particle_set)
     441             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: pv
     442             :       TYPE(local_colvar_constraint_type), INTENT(IN)     :: lcolv
     443             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     444             : 
     445             :       INTEGER                                            :: i, iatm, ind, j
     446             :       REAL(KIND=dp)                                      :: lambda, tmp
     447             :       REAL(KIND=dp), DIMENSION(3)                        :: f
     448             : 
     449       81829 :       DO iatm = 1, SIZE(lcolv%colvar_old%i_atom)
     450       56128 :          ind = lcolv%colvar_old%i_atom(iatm)
     451      224512 :          f = -lcolv%colvar_old%dsdr(:, iatm)
     452             :          !  pv gets updated with FULL multiplier
     453       56128 :          lambda = lcolv%lambda
     454      250213 :          DO i = 1, 3
     455      168384 :             tmp = lambda*particle_set(ind)%r(i)
     456      729664 :             DO j = 1, 3
     457      673536 :                pv(j, i) = pv(j, i) + f(j)*tmp
     458             :             END DO
     459             :          END DO
     460             :       END DO
     461             : 
     462       25701 :    END SUBROUTINE pv_colv_eval
     463             : 
     464             : ! **************************************************************************************************
     465             : !> \brief ...
     466             : !> \param roll_tol ...
     467             : !> \param iroll ...
     468             : !> \param char ...
     469             : !> \param matrix ...
     470             : !> \param veps ...
     471             : !> \par History
     472             : !>      none
     473             : ! **************************************************************************************************
     474        3620 :    SUBROUTINE check_tol(roll_tol, iroll, char, matrix, veps)
     475             : 
     476             :       REAL(KIND=dp), INTENT(OUT)                         :: roll_tol
     477             :       INTEGER, INTENT(INOUT)                             :: iroll
     478             :       CHARACTER(LEN=*), INTENT(IN)                       :: char
     479             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     480             :          OPTIONAL                                        :: matrix, veps
     481             : 
     482             :       REAL(KIND=dp)                                      :: local_tol
     483             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: diff_rattle, diff_shake
     484             :       REAL(KIND=dp), DIMENSION(3, 3), SAVE               :: matrix_old, veps_old
     485             : 
     486        1826 :       SELECT CASE (char)
     487             :       CASE ('SHAKE')
     488        1826 :          IF (iroll == 1) THEN
     489        8814 :             matrix_old = matrix
     490         678 :             roll_tol = -1.E10_dp
     491             :          ELSE
     492             :             roll_tol = 0.0_dp
     493       14924 :             diff_shake = ABS(matrix_old - matrix)
     494       14924 :             local_tol = MAXVAL(diff_shake)
     495        1148 :             roll_tol = MAX(roll_tol, local_tol)
     496       14924 :             matrix_old = matrix
     497             :          END IF
     498        1826 :          iroll = iroll + 1
     499             :       CASE ('RATTLE')
     500        1794 :          IF (iroll == 1) THEN
     501        8814 :             veps_old = veps
     502         678 :             roll_tol = -1.E+10_dp
     503             :          ELSE
     504             :             roll_tol = 0.0_dp
     505             :             ! compute tolerance on veps
     506       14508 :             diff_rattle = ABS(veps - veps_old)
     507       14508 :             local_tol = MAXVAL(diff_rattle)
     508        1116 :             roll_tol = MAX(roll_tol, local_tol)
     509       14508 :             veps_old = veps
     510             :          END IF
     511        5414 :          iroll = iroll + 1
     512             :       END SELECT
     513             : 
     514        3620 :    END SUBROUTINE check_tol
     515             : 
     516             : ! **************************************************************************************************
     517             : !> \brief ...
     518             : !> \param char ...
     519             : !> \param r_shake ...
     520             : !> \param v_shake ...
     521             : !> \param vector_r ...
     522             : !> \param vector_v ...
     523             : !> \param u ...
     524             : !> \par History
     525             : !>      none
     526             : ! **************************************************************************************************
     527        3620 :    SUBROUTINE get_roll_matrix(char, r_shake, v_shake, vector_r, vector_v, u)
     528             : 
     529             :       CHARACTER(len=*), INTENT(IN)                       :: char
     530             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
     531             :          OPTIONAL                                        :: r_shake, v_shake
     532             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: vector_r, vector_v
     533             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     534             :          OPTIONAL                                        :: u
     535             : 
     536             :       INTEGER                                            :: i
     537             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: diag
     538             : 
     539       25532 :       IF (PRESENT(r_shake)) r_shake = 0.0_dp
     540       47060 :       IF (PRESENT(v_shake)) v_shake = 0.0_dp
     541        3620 :       diag = 0.0_dp
     542             : 
     543        1826 :       SELECT CASE (char)
     544             :       CASE ('SHAKE')
     545        1826 :          IF (PRESENT(u) .AND. PRESENT(vector_v) .AND. &
     546        1794 :              PRESENT(vector_r)) THEN
     547          20 :             diag(1, 1) = vector_r(1)
     548          20 :             diag(2, 2) = vector_r(2)
     549          20 :             diag(3, 3) = vector_r(3)
     550        2660 :             r_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
     551          20 :             diag(1, 1) = vector_v(1)
     552          20 :             diag(2, 2) = vector_v(2)
     553          20 :             diag(3, 3) = vector_v(3)
     554        2660 :             v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
     555         800 :             diag = MATMUL(r_shake, v_shake)
     556         260 :             r_shake = diag
     557        1806 :          ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v) .AND. &
     558             :                  PRESENT(vector_r)) THEN
     559        7224 :             DO i = 1, 3
     560        5418 :                r_shake(i, i) = vector_r(i)*vector_v(i)
     561        7224 :                v_shake(i, i) = vector_v(i)
     562             :             END DO
     563             :          ELSE
     564           0 :             CPABORT("Not sufficient parameters")
     565             :          END IF
     566             :       CASE ('RATTLE')
     567        3620 :          IF (PRESENT(u) .AND. PRESENT(vector_v)) THEN
     568          20 :             diag(1, 1) = vector_v(1)
     569          20 :             diag(2, 2) = vector_v(2)
     570          20 :             diag(3, 3) = vector_v(3)
     571        2660 :             v_shake = MATMUL(MATMUL(u, diag), TRANSPOSE(u))
     572        1774 :          ELSEIF (.NOT. PRESENT(u) .AND. PRESENT(vector_v)) THEN
     573        7096 :             DO i = 1, 3
     574        7096 :                v_shake(i, i) = vector_v(i)
     575             :             END DO
     576             :          ELSE
     577           0 :             CPABORT("Not sufficient parameters")
     578             :          END IF
     579             :       END SELECT
     580             : 
     581        3620 :    END SUBROUTINE get_roll_matrix
     582             : 
     583             : ! **************************************************************************************************
     584             : !> \brief ...
     585             : !> \param particle_set ...
     586             : !> \param local_particles ...
     587             : !> \param pos ...
     588             : !> \param vel ...
     589             : !> \par History
     590             : !>      Teodoro Laino [tlaino] 2007
     591             : ! **************************************************************************************************
     592        3345 :    SUBROUTINE restore_temporary_set(particle_set, local_particles, pos, vel)
     593             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     594             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     595             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     596             :          OPTIONAL                                        :: pos, vel
     597             : 
     598             :       INTEGER                                            :: iparticle, iparticle_kind, &
     599             :                                                             iparticle_local, nparticle_local
     600        3345 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: wrk
     601             : 
     602       10035 :       ALLOCATE (wrk(SIZE(particle_set)))
     603      187160 :       wrk = .TRUE.
     604        8555 :       DO iparticle_kind = 1, SIZE(local_particles%n_el)
     605        5210 :          nparticle_local = local_particles%n_el(iparticle_kind)
     606      100608 :          DO iparticle_local = 1, nparticle_local
     607       92053 :             iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     608       97263 :             wrk(iparticle) = .FALSE.
     609             :          END DO
     610             :       END DO
     611        3345 :       IF (PRESENT(vel)) THEN
     612      187160 :          DO iparticle = 1, SIZE(particle_set)
     613      187160 :             IF (wrk(iparticle)) THEN
     614      367048 :                vel(:, iparticle) = 0.0_dp
     615             :             END IF
     616             :          END DO
     617             :       END IF
     618        3345 :       IF (PRESENT(pos)) THEN
     619       98242 :          DO iparticle = 1, SIZE(particle_set)
     620       98242 :             IF (wrk(iparticle)) THEN
     621      192504 :                pos(:, iparticle) = 0.0_dp
     622             :             END IF
     623             :          END DO
     624             :       END IF
     625        3345 :       DEALLOCATE (wrk)
     626        3345 :    END SUBROUTINE restore_temporary_set
     627             : 
     628             : ! **************************************************************************************************
     629             : !> \brief ...
     630             : !> \param group ...
     631             : !> \param pos ...
     632             : !> \param vel ...
     633             : !> \par History
     634             : !>      Teodoro Laino [tlaino] 2007
     635             : ! **************************************************************************************************
     636        3345 :    SUBROUTINE update_temporary_set(group, pos, vel)
     637             :       CLASS(mp_comm_type), INTENT(IN)                     :: group
     638             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
     639             :          OPTIONAL                                        :: pos, vel
     640             : 
     641        3345 :       IF (PRESENT(pos)) THEN
     642      773035 :          CALL group%sum(pos)
     643             :       END IF
     644        3345 :       IF (PRESENT(vel)) THEN
     645     1473865 :          CALL group%sum(vel)
     646             :       END IF
     647        3345 :    END SUBROUTINE update_temporary_set
     648             : 
     649          60 : END MODULE constraint_util

Generated by: LCOV version 1.15