LCOV - code coverage report
Current view: top level - src - constraint_clv.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 235 249 94.4 %
Date: 2024-03-28 07:31:50 Functions: 16 18 88.9 %

          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 Module that handles the COLLECTIVE constraints
      10             : !> \par History
      11             : !>      none
      12             : !> \author Teodoro Laino [tlaino]
      13             : ! **************************************************************************************************
      14             : MODULE constraint_clv
      15             :    USE cell_types,                      ONLY: cell_type
      16             :    USE colvar_methods,                  ONLY: colvar_eval_mol_f
      17             :    USE colvar_types,                    ONLY: colvar_type,&
      18             :                                               diff_colvar
      19             :    USE input_section_types,             ONLY: section_vals_get,&
      20             :                                               section_vals_get_subs_vals,&
      21             :                                               section_vals_type,&
      22             :                                               section_vals_val_get,&
      23             :                                               section_vals_val_set
      24             :    USE kinds,                           ONLY: dp
      25             :    USE molecule_kind_types,             ONLY: colvar_constraint_type,&
      26             :                                               fixd_constraint_type,&
      27             :                                               get_molecule_kind,&
      28             :                                               molecule_kind_type
      29             :    USE molecule_types,                  ONLY: get_molecule,&
      30             :                                               global_constraint_type,&
      31             :                                               local_colvar_constraint_type,&
      32             :                                               molecule_type
      33             :    USE particle_types,                  ONLY: particle_type
      34             : #include "./base/base_uses.f90"
      35             : 
      36             :    IMPLICIT NONE
      37             : 
      38             :    PRIVATE
      39             :    PUBLIC :: shake_roll_colv_int, &
      40             :              rattle_roll_colv_int, &
      41             :              shake_colv_int, &
      42             :              rattle_colv_int, &
      43             :              shake_roll_colv_ext, &
      44             :              rattle_roll_colv_ext, &
      45             :              shake_colv_ext, &
      46             :              rattle_colv_ext, &
      47             :              shake_update_colv_int, &
      48             :              shake_update_colv_ext
      49             : 
      50             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'constraint_clv'
      51             : 
      52             : CONTAINS
      53             : 
      54             : ! **************************************************************************************************
      55             : !> \brief Intramolecular subroutine
      56             : !>      shake_colv algorithm for collective variables constraints
      57             : !>      updates the multiplier one molecule type at a time
      58             : !> \param molecule ...
      59             : !> \param particle_set ...
      60             : !> \param pos ...
      61             : !> \param vel ...
      62             : !> \param dt ...
      63             : !> \param ishake ...
      64             : !> \param cell ...
      65             : !> \param imass ...
      66             : !> \param max_sigma ...
      67             : !> \par History
      68             : !>      none
      69             : !> \author Teodoro Laino [tlaino]
      70             : ! **************************************************************************************************
      71       18698 :    SUBROUTINE shake_colv_int(molecule, particle_set, pos, vel, dt, ishake, &
      72       18698 :                              cell, imass, max_sigma)
      73             : 
      74             :       TYPE(molecule_type), POINTER                       :: molecule
      75             :       TYPE(particle_type), POINTER                       :: particle_set(:)
      76             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
      77             :       REAL(kind=dp), INTENT(in)                          :: dt
      78             :       INTEGER, INTENT(IN)                                :: ishake
      79             :       TYPE(cell_type), POINTER                           :: cell
      80             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
      81             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
      82             : 
      83       18698 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
      84       18698 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
      85       18698 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
      86             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
      87             : 
      88       18698 :       NULLIFY (fixd_list)
      89       18698 :       molecule_kind => molecule%molecule_kind
      90       18698 :       CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
      91       18698 :       CALL get_molecule(molecule, lcolv=lcolv)
      92             :       ! Real Shake
      93             :       CALL shake_colv_low(fixd_list, colv_list, lcolv, &
      94       18698 :                           particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
      95             : 
      96       18698 :    END SUBROUTINE shake_colv_int
      97             : 
      98             : ! **************************************************************************************************
      99             : !> \brief Intramolecular subroutine for updating the TARGET value of collective
     100             : !>        constraints
     101             : !> \param molecule ...
     102             : !> \param dt ...
     103             : !> \param motion_section ...
     104             : !> \author Teodoro Laino [tlaino] - University of Zurich
     105             : ! **************************************************************************************************
     106        2591 :    SUBROUTINE shake_update_colv_int(molecule, dt, motion_section)
     107             : 
     108             :       TYPE(molecule_type), POINTER                       :: molecule
     109             :       REAL(kind=dp), INTENT(in)                          :: dt
     110             :       TYPE(section_vals_type), POINTER                   :: motion_section
     111             : 
     112        2591 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     113             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     114             : 
     115        2591 :       molecule_kind => molecule%molecule_kind
     116        2591 :       CALL get_molecule_kind(molecule_kind, colv_list=colv_list)
     117             :       ! Real update of the Shake target
     118        2591 :       CALL shake_update_colv_low(colv_list, dt, motion_section)
     119             : 
     120        2591 :    END SUBROUTINE shake_update_colv_int
     121             : 
     122             : ! **************************************************************************************************
     123             : !> \brief Intramolecular subroutine
     124             : !>      rattle algorithm for collective variables constraints
     125             : !>      updates the multiplier one molecule type at a time
     126             : !> \param molecule ...
     127             : !> \param particle_set ...
     128             : !> \param vel ...
     129             : !> \param dt ...
     130             : !> \param irattle ...
     131             : !> \param cell ...
     132             : !> \param imass ...
     133             : !> \param max_sigma ...
     134             : !> \par History
     135             : !>      none
     136             : !> \author Teodoro Laino [tlaino]
     137             : ! **************************************************************************************************
     138        7833 :    SUBROUTINE rattle_colv_int(molecule, particle_set, vel, dt, irattle, &
     139        7833 :                               cell, imass, max_sigma)
     140             : 
     141             :       TYPE(molecule_type), POINTER                       :: molecule
     142             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     143             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     144             :       REAL(kind=dp), INTENT(in)                          :: dt
     145             :       INTEGER, INTENT(IN)                                :: irattle
     146             :       TYPE(cell_type), POINTER                           :: cell
     147             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     148             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     149             : 
     150        7833 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     151        7833 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     152        7833 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     153             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     154             : 
     155        7833 :       NULLIFY (fixd_list)
     156        7833 :       molecule_kind => molecule%molecule_kind
     157        7833 :       CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
     158        7833 :       CALL get_molecule(molecule, lcolv=lcolv)
     159             :       ! Real Rattle
     160             :       CALL rattle_colv_low(fixd_list, colv_list, lcolv, &
     161        7833 :                            particle_set, vel, dt, irattle, cell, imass, max_sigma)
     162             : 
     163        7833 :    END SUBROUTINE rattle_colv_int
     164             : 
     165             : ! **************************************************************************************************
     166             : !> \brief Intramolecular subroutine
     167             : !>      shake algorithm (box allowed to change) for collective variables constraints
     168             : !>      updates the multiplier one molecule type at a time
     169             : !> \param molecule ...
     170             : !> \param particle_set ...
     171             : !> \param pos ...
     172             : !> \param vel ...
     173             : !> \param r_shake ...
     174             : !> \param v_shake ...
     175             : !> \param dt ...
     176             : !> \param ishake ...
     177             : !> \param cell ...
     178             : !> \param imass ...
     179             : !> \param max_sigma ...
     180             : !> \par History
     181             : !>      none
     182             : !> \author Teodoro Laino [tlaino]
     183             : ! **************************************************************************************************
     184       15773 :    SUBROUTINE shake_roll_colv_int(molecule, particle_set, pos, vel, r_shake, v_shake, &
     185       15773 :                                   dt, ishake, cell, imass, max_sigma)
     186             : 
     187             :       TYPE(molecule_type), POINTER                       :: molecule
     188             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     189             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     190             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake, v_shake
     191             :       REAL(kind=dp), INTENT(in)                          :: dt
     192             :       INTEGER, INTENT(in)                                :: ishake
     193             :       TYPE(cell_type), POINTER                           :: cell
     194             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     195             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     196             : 
     197       15773 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     198       15773 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     199       15773 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     200             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     201             : 
     202       15773 :       NULLIFY (fixd_list)
     203       15773 :       molecule_kind => molecule%molecule_kind
     204       15773 :       CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
     205       15773 :       CALL get_molecule(molecule, lcolv=lcolv)
     206             :       ! Real Shake
     207             :       CALL shake_roll_colv_low(fixd_list, colv_list, lcolv, &
     208             :                                particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
     209       15773 :                                imass, max_sigma)
     210             : 
     211       15773 :    END SUBROUTINE shake_roll_colv_int
     212             : 
     213             : ! **************************************************************************************************
     214             : !> \brief Intramolecular subroutine
     215             : !>      rattle algorithm (box allowed to change) for collective variables constraints
     216             : !>      updates the multiplier one molecule type at a time
     217             : !> \param molecule ...
     218             : !> \param particle_set ...
     219             : !> \param vel ...
     220             : !> \param r_rattle ...
     221             : !> \param dt ...
     222             : !> \param irattle ...
     223             : !> \param veps ...
     224             : !> \param cell ...
     225             : !> \param imass ...
     226             : !> \param max_sigma ...
     227             : !> \par History
     228             : !>      none
     229             : !> \author Teodoro Laino [tlaino]
     230             : ! **************************************************************************************************
     231        5521 :    SUBROUTINE rattle_roll_colv_int(molecule, particle_set, vel, r_rattle, &
     232        5521 :                                    dt, irattle, veps, cell, imass, max_sigma)
     233             : 
     234             :       TYPE(molecule_type), POINTER                       :: molecule
     235             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     236             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     237             :       REAL(KIND=dp), INTENT(IN)                          :: r_rattle(:, :), dt
     238             :       INTEGER, INTENT(in)                                :: irattle
     239             :       REAL(KIND=dp), INTENT(IN)                          :: veps(:, :)
     240             :       TYPE(cell_type), POINTER                           :: cell
     241             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     242             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     243             : 
     244        5521 :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     245        5521 :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     246        5521 :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     247             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     248             : 
     249        5521 :       NULLIFY (fixd_list)
     250        5521 :       molecule_kind => molecule%molecule_kind
     251        5521 :       CALL get_molecule_kind(molecule_kind, colv_list=colv_list, fixd_list=fixd_list)
     252        5521 :       CALL get_molecule(molecule, lcolv=lcolv)
     253             :       ! Real Rattle
     254             :       CALL rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
     255             :                                 particle_set, vel, r_rattle, dt, irattle, veps, cell, &
     256        5521 :                                 imass, max_sigma)
     257             : 
     258        5521 :    END SUBROUTINE rattle_roll_colv_int
     259             : 
     260             : ! **************************************************************************************************
     261             : !> \brief Intermolecular subroutine
     262             : !>      shake_colv algorithm for collective variables constraints
     263             : !>      updates the multiplier one molecule type at a time
     264             : !> \param gci ...
     265             : !> \param particle_set ...
     266             : !> \param pos ...
     267             : !> \param vel ...
     268             : !> \param dt ...
     269             : !> \param ishake ...
     270             : !> \param cell ...
     271             : !> \param imass ...
     272             : !> \param max_sigma ...
     273             : !> \par History
     274             : !>      none
     275             : !> \author Teodoro Laino [tlaino]
     276             : ! **************************************************************************************************
     277        1719 :    SUBROUTINE shake_colv_ext(gci, particle_set, pos, vel, dt, ishake, &
     278        1719 :                              cell, imass, max_sigma)
     279             : 
     280             :       TYPE(global_constraint_type), POINTER              :: gci
     281             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     282             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     283             :       REAL(kind=dp), INTENT(in)                          :: dt
     284             :       INTEGER, INTENT(IN)                                :: ishake
     285             :       TYPE(cell_type), POINTER                           :: cell
     286             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     287             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     288             : 
     289             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     290             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     291             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     292             : 
     293        1719 :       colv_list => gci%colv_list
     294        1719 :       fixd_list => gci%fixd_list
     295        1719 :       lcolv => gci%lcolv
     296             :       ! Real Shake
     297             :       CALL shake_colv_low(fixd_list, colv_list, lcolv, &
     298        1719 :                           particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
     299             : 
     300        1719 :    END SUBROUTINE shake_colv_ext
     301             : 
     302             : ! **************************************************************************************************
     303             : !> \brief Intermolecular subroutine for updating the TARGET value for collective
     304             : !>        constraints
     305             : !> \param gci ...
     306             : !> \param dt ...
     307             : !> \param motion_section ...
     308             : !> \author Teodoro Laino [tlaino]
     309             : ! **************************************************************************************************
     310        1094 :    SUBROUTINE shake_update_colv_ext(gci, dt, motion_section)
     311             : 
     312             :       TYPE(global_constraint_type), POINTER              :: gci
     313             :       REAL(kind=dp), INTENT(in)                          :: dt
     314             :       TYPE(section_vals_type), POINTER                   :: motion_section
     315             : 
     316             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     317             : 
     318        1094 :       colv_list => gci%colv_list
     319             :       ! Real update of the Shake target
     320        1094 :       CALL shake_update_colv_low(colv_list, dt, motion_section)
     321             : 
     322        1094 :    END SUBROUTINE shake_update_colv_ext
     323             : 
     324             : ! **************************************************************************************************
     325             : !> \brief Intermolecular subroutine
     326             : !>      rattle algorithm for collective variables constraints
     327             : !>      updates the multiplier one molecule type at a time
     328             : !> \param gci ...
     329             : !> \param particle_set ...
     330             : !> \param vel ...
     331             : !> \param dt ...
     332             : !> \param irattle ...
     333             : !> \param cell ...
     334             : !> \param imass ...
     335             : !> \param max_sigma ...
     336             : !> \par History
     337             : !>      none
     338             : !> \author Teodoro Laino [tlaino]
     339             : ! **************************************************************************************************
     340        1442 :    SUBROUTINE rattle_colv_ext(gci, particle_set, vel, dt, irattle, &
     341        1442 :                               cell, imass, max_sigma)
     342             : 
     343             :       TYPE(global_constraint_type), POINTER              :: gci
     344             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     345             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     346             :       REAL(kind=dp), INTENT(in)                          :: dt
     347             :       INTEGER, INTENT(IN)                                :: irattle
     348             :       TYPE(cell_type), POINTER                           :: cell
     349             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     350             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     351             : 
     352             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     353             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     354             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     355             : 
     356        1442 :       colv_list => gci%colv_list
     357        1442 :       fixd_list => gci%fixd_list
     358        1442 :       lcolv => gci%lcolv
     359             :       ! Real Rattle
     360             :       CALL rattle_colv_low(fixd_list, colv_list, lcolv, &
     361        1442 :                            particle_set, vel, dt, irattle, cell, imass, max_sigma)
     362             : 
     363        1442 :    END SUBROUTINE rattle_colv_ext
     364             : 
     365             : ! **************************************************************************************************
     366             : !> \brief Intermolecular subroutine
     367             : !>      shake algorithm (box allowed to change) for collective variables constraints
     368             : !>      updates the multiplier one molecule type at a time
     369             : !> \param gci ...
     370             : !> \param particle_set ...
     371             : !> \param pos ...
     372             : !> \param vel ...
     373             : !> \param r_shake ...
     374             : !> \param v_shake ...
     375             : !> \param dt ...
     376             : !> \param ishake ...
     377             : !> \param cell ...
     378             : !> \param imass ...
     379             : !> \param max_sigma ...
     380             : !> \par History
     381             : !>      none
     382             : !> \author Teodoro Laino [tlaino]
     383             : ! **************************************************************************************************
     384           0 :    SUBROUTINE shake_roll_colv_ext(gci, particle_set, pos, vel, r_shake, v_shake, &
     385           0 :                                   dt, ishake, cell, imass, max_sigma)
     386             : 
     387             :       TYPE(global_constraint_type), POINTER              :: gci
     388             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     389             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     390             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake, v_shake
     391             :       REAL(kind=dp), INTENT(in)                          :: dt
     392             :       INTEGER, INTENT(in)                                :: ishake
     393             :       TYPE(cell_type), POINTER                           :: cell
     394             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     395             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     396             : 
     397             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     398             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     399             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     400             : 
     401           0 :       colv_list => gci%colv_list
     402           0 :       fixd_list => gci%fixd_list
     403           0 :       lcolv => gci%lcolv
     404             :       ! Real Shake
     405             :       CALL shake_roll_colv_low(fixd_list, colv_list, lcolv, &
     406             :                                particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
     407           0 :                                imass, max_sigma)
     408             : 
     409           0 :    END SUBROUTINE shake_roll_colv_ext
     410             : 
     411             : ! **************************************************************************************************
     412             : !> \brief Intermolecular subroutine
     413             : !>      rattle algorithm (box allowed to change) for collective variables constraints
     414             : !>      updates the multiplier one molecule type at a time
     415             : !> \param gci ...
     416             : !> \param particle_set ...
     417             : !> \param vel ...
     418             : !> \param r_rattle ...
     419             : !> \param dt ...
     420             : !> \param irattle ...
     421             : !> \param veps ...
     422             : !> \param cell ...
     423             : !> \param imass ...
     424             : !> \param max_sigma ...
     425             : !> \par History
     426             : !>      none
     427             : !> \author Teodoro Laino [tlaino]
     428             : ! **************************************************************************************************
     429           0 :    SUBROUTINE rattle_roll_colv_ext(gci, particle_set, vel, r_rattle, &
     430           0 :                                    dt, irattle, veps, cell, imass, max_sigma)
     431             : 
     432             :       TYPE(global_constraint_type), POINTER              :: gci
     433             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     434             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     435             :       REAL(KIND=dp), INTENT(IN)                          :: r_rattle(:, :), dt
     436             :       INTEGER, INTENT(in)                                :: irattle
     437             :       REAL(KIND=dp), INTENT(IN)                          :: veps(:, :)
     438             :       TYPE(cell_type), POINTER                           :: cell
     439             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     440             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     441             : 
     442             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     443             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     444             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     445             : 
     446           0 :       colv_list => gci%colv_list
     447           0 :       fixd_list => gci%fixd_list
     448           0 :       lcolv => gci%lcolv
     449             :       ! Real Rattle
     450             :       CALL rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
     451             :                                 particle_set, vel, r_rattle, dt, irattle, veps, cell, &
     452           0 :                                 imass, max_sigma)
     453             : 
     454           0 :    END SUBROUTINE rattle_roll_colv_ext
     455             : 
     456             : ! **************************************************************************************************
     457             : !> \brief Real Shake subroutine - Low Level
     458             : !>      shake_colv algorithm for collective variables constraints
     459             : !>      updates the multiplier one molecule type at a time
     460             : !> \param fixd_list ...
     461             : !> \param colv_list ...
     462             : !> \param lcolv ...
     463             : !> \param particle_set ...
     464             : !> \param pos ...
     465             : !> \param vel ...
     466             : !> \param dt ...
     467             : !> \param ishake ...
     468             : !> \param cell ...
     469             : !> \param imass ...
     470             : !> \param max_sigma ...
     471             : !> \par History
     472             : !>      none
     473             : !> \author Teodoro Laino [tlaino]
     474             : ! **************************************************************************************************
     475       20417 :    SUBROUTINE shake_colv_low(fixd_list, colv_list, lcolv, &
     476       20417 :                              particle_set, pos, vel, dt, ishake, cell, imass, max_sigma)
     477             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     478             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     479             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     480             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     481             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     482             :       REAL(kind=dp), INTENT(in)                          :: dt
     483             :       INTEGER, INTENT(IN)                                :: ishake
     484             :       TYPE(cell_type), POINTER                           :: cell
     485             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     486             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     487             : 
     488             :       INTEGER                                            :: iconst
     489             :       REAL(KIND=dp)                                      :: del_lam, dtby2, dtsqby2, fdotf_sum
     490             : 
     491       20417 :       dtsqby2 = dt*dt*.5_dp
     492       20417 :       dtby2 = dt*.5_dp
     493       20417 :       IF (ishake == 1) THEN
     494        9800 :          DO iconst = 1, SIZE(colv_list)
     495        6775 :             IF (colv_list(iconst)%restraint%active) CYCLE
     496             :             ! Update positions
     497             :             CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
     498             :                                  lambda=lcolv(iconst)%lambda, &
     499        5855 :                                  imass=imass)
     500             :             ! Update velocities
     501             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     502             :                                  lambda=lcolv(iconst)%lambda, &
     503        9800 :                                  imass=imass)
     504             :          END DO
     505             :       ELSE
     506      107106 :          DO iconst = 1, SIZE(colv_list)
     507       89714 :             IF (colv_list(iconst)%restraint%active) CYCLE
     508             :             ! Update colvar
     509             :             CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
     510       89680 :                                    pos=pos, fixd_list=fixd_list)
     511             :             lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
     512       89680 :                                               colv_list(iconst)%expected_value)
     513             :             fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar, &
     514       89680 :                                         lcolv(iconst)%colvar_old, imass=imass)
     515       89680 :             del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*dt*fdotf_sum)
     516       89680 :             lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
     517             : 
     518             :             ! Update positions
     519             :             CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
     520             :                                  lambda=del_lam, &
     521       89680 :                                  imass=imass)
     522             :             ! Update velocities
     523             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     524             :                                  lambda=del_lam, &
     525      107106 :                                  imass=imass)
     526             :          END DO
     527             :       END IF
     528             :       ! computing the constraint and value of tolerance
     529      116906 :       DO iconst = 1, SIZE(colv_list)
     530       96489 :          IF (colv_list(iconst)%restraint%active) CYCLE
     531             :          CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
     532       95535 :                                 pos=pos, fixd_list=fixd_list)
     533             :          lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
     534       95535 :                                            colv_list(iconst)%expected_value)
     535      116906 :          max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
     536             :       END DO
     537       20417 :    END SUBROUTINE shake_colv_low
     538             : 
     539             : ! **************************************************************************************************
     540             : !> \brief Real Shake subroutine - Low Level - for updating the TARGET value
     541             : !> \param colv_list ...
     542             : !> \param dt ...
     543             : !> \param motion_section ...
     544             : !> \date 02.2008
     545             : !> \author Teodoro Laino [tlaino] - University of Zurich
     546             : ! **************************************************************************************************
     547        7370 :    SUBROUTINE shake_update_colv_low(colv_list, dt, motion_section)
     548             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     549             :       REAL(kind=dp), INTENT(in)                          :: dt
     550             :       TYPE(section_vals_type), POINTER                   :: motion_section
     551             : 
     552             :       INTEGER                                            :: iconst, irep, n_rep
     553             :       LOGICAL                                            :: do_update_colvar, explicit
     554             :       REAL(KIND=dp)                                      :: clv_target, limit, new_clv_target, value
     555             :       TYPE(section_vals_type), POINTER                   :: collective_sections
     556             : 
     557             : ! Update globally for restart
     558             : 
     559        3685 :       collective_sections => section_vals_get_subs_vals(motion_section, "CONSTRAINT%COLLECTIVE")
     560        3685 :       CALL section_vals_get(collective_sections, n_repetition=n_rep)
     561        3685 :       IF (n_rep /= 0) THEN
     562       10798 :          DO irep = 1, n_rep
     563             :             CALL section_vals_val_get(collective_sections, "TARGET_GROWTH", r_val=value, &
     564        7813 :                                       i_rep_section=irep)
     565       10798 :             IF (value /= 0.0_dp) THEN
     566             :                CALL section_vals_val_get(collective_sections, "TARGET", r_val=clv_target, &
     567         200 :                                          i_rep_section=irep)
     568         200 :                new_clv_target = clv_target + value*dt
     569             :                ! Check limits..
     570             :                CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", explicit=explicit, &
     571         200 :                                          i_rep_section=irep)
     572         200 :                do_update_colvar = .TRUE.
     573         200 :                IF (explicit) THEN
     574             :                   CALL section_vals_val_get(collective_sections, "TARGET_LIMIT", r_val=limit, &
     575         100 :                                             i_rep_section=irep)
     576         100 :                   IF (value > 0.0_dp) THEN
     577          50 :                      IF (clv_target == limit) THEN
     578             :                         do_update_colvar = .FALSE.
     579          22 :                      ELSE IF (new_clv_target >= limit) THEN
     580           1 :                         new_clv_target = limit
     581             :                      END IF
     582             :                   ELSE
     583          50 :                      IF (clv_target == limit) THEN
     584             :                         do_update_colvar = .FALSE.
     585          39 :                      ELSE IF (new_clv_target <= limit) THEN
     586           1 :                         new_clv_target = limit
     587             :                      END IF
     588             :                   END IF
     589             :                END IF
     590             :                IF (do_update_colvar) THEN
     591             :                   CALL section_vals_val_set(collective_sections, "TARGET", r_val=new_clv_target, &
     592         161 :                                             i_rep_section=irep)
     593             :                END IF
     594             :             END IF
     595             :          END DO
     596             :       END IF
     597             : 
     598             :       ! Update locally the value to each processor
     599       13318 :       DO iconst = 1, SIZE(colv_list)
     600             :          ! Update local to each processor
     601        9633 :          IF (colv_list(iconst)%expected_value_growth_speed == 0.0_dp) CYCLE
     602             :          CALL section_vals_val_get(collective_sections, "TARGET", &
     603             :                                    r_val=colv_list(iconst)%expected_value, &
     604       13318 :                                    i_rep_section=colv_list(iconst)%inp_seq_num)
     605             :       END DO
     606             : 
     607        3685 :    END SUBROUTINE shake_update_colv_low
     608             : 
     609             : ! **************************************************************************************************
     610             : !> \brief Real Rattle - Low Level
     611             : !>      rattle algorithm for collective variables constraints
     612             : !>      updates the multiplier one molecule type at a time
     613             : !> \param fixd_list ...
     614             : !> \param colv_list ...
     615             : !> \param lcolv ...
     616             : !> \param particle_set ...
     617             : !> \param vel ...
     618             : !> \param dt ...
     619             : !> \param irattle ...
     620             : !> \param cell ...
     621             : !> \param imass ...
     622             : !> \param max_sigma ...
     623             : !> \par History
     624             : !>      none
     625             : !> \author Teodoro Laino [tlaino]
     626             : ! **************************************************************************************************
     627        9275 :    SUBROUTINE rattle_colv_low(fixd_list, colv_list, lcolv, &
     628        9275 :                               particle_set, vel, dt, irattle, cell, imass, max_sigma)
     629             : 
     630             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     631             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     632             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     633             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     634             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     635             :       REAL(kind=dp), INTENT(in)                          :: dt
     636             :       INTEGER, INTENT(IN)                                :: irattle
     637             :       TYPE(cell_type), POINTER                           :: cell
     638             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     639             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     640             : 
     641             :       INTEGER                                            :: iconst
     642             :       REAL(KIND=dp)                                      :: del_lam, dtby2, fdotf_sum
     643             : 
     644        9275 :       dtby2 = dt*.5_dp
     645        9275 :       IF (irattle == 1) THEN
     646        9812 :          DO iconst = 1, SIZE(colv_list)
     647        6781 :             IF (colv_list(iconst)%restraint%active) CYCLE
     648             :             ! Update colvar_old
     649             :             CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
     650        5861 :                                    particles=particle_set, fixd_list=fixd_list)
     651             :             ! Update velocities
     652             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     653             :                                  lambda=lcolv(iconst)%lambda, &
     654        9812 :                                  imass=imass)
     655             :          END DO
     656             :       ELSE
     657       35597 :          DO iconst = 1, SIZE(colv_list)
     658       29353 :             IF (colv_list(iconst)%restraint%active) CYCLE
     659       29333 :             lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel)
     660             :             fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar_old, &
     661       29333 :                                         lcolv(iconst)%colvar_old, imass=imass)
     662       29333 :             del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*fdotf_sum)
     663       29333 :             lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
     664             : 
     665             :             ! Update velocities
     666             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     667             :                                  lambda=del_lam, &
     668       35597 :                                  imass=imass)
     669             :          END DO
     670             :       END IF
     671             : 
     672       45409 :       DO iconst = 1, SIZE(colv_list)
     673       36134 :          IF (colv_list(iconst)%restraint%active) CYCLE
     674       35194 :          lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel)
     675       45409 :          max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
     676             :       END DO
     677             : 
     678        9275 :    END SUBROUTINE rattle_colv_low
     679             : 
     680             : ! **************************************************************************************************
     681             : !> \brief Real shake_roll - Low Level
     682             : !>      shake algorithm (box allowed to change) for collective variables constraints
     683             : !>      updates the multiplier one molecule type at a time
     684             : !> \param fixd_list ...
     685             : !> \param colv_list ...
     686             : !> \param lcolv ...
     687             : !> \param particle_set ...
     688             : !> \param pos ...
     689             : !> \param vel ...
     690             : !> \param r_shake ...
     691             : !> \param v_shake ...
     692             : !> \param dt ...
     693             : !> \param ishake ...
     694             : !> \param cell ...
     695             : !> \param imass ...
     696             : !> \param max_sigma ...
     697             : !> \par History
     698             : !>      none
     699             : !> \author Teodoro Laino [tlaino]
     700             : ! **************************************************************************************************
     701       15773 :    SUBROUTINE shake_roll_colv_low(fixd_list, colv_list, lcolv, &
     702       31546 :                                   particle_set, pos, vel, r_shake, v_shake, dt, ishake, cell, &
     703       15773 :                                   imass, max_sigma)
     704             : 
     705             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     706             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     707             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     708             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     709             :       REAL(KIND=dp), INTENT(INOUT)                       :: pos(:, :), vel(:, :)
     710             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: r_shake, v_shake
     711             :       REAL(kind=dp), INTENT(in)                          :: dt
     712             :       INTEGER, INTENT(in)                                :: ishake
     713             :       TYPE(cell_type), POINTER                           :: cell
     714             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     715             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     716             : 
     717             :       INTEGER                                            :: iconst
     718             :       REAL(KIND=dp)                                      :: del_lam, dtby2, dtsqby2, fdotf_sum
     719             : 
     720       15773 :       dtsqby2 = dt*dt*.5_dp
     721       15773 :       dtby2 = dt*.5_dp
     722       15773 :       IF (ishake == 1) THEN
     723        8664 :          DO iconst = 1, SIZE(colv_list)
     724        7312 :             IF (colv_list(iconst)%restraint%active) CYCLE
     725             :             ! Update positions
     726             :             CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
     727             :                                  lambda=lcolv(iconst)%lambda, &
     728        7312 :                                  roll=.TRUE., rmat=r_shake, imass=imass)
     729             :             ! Update velocities
     730             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     731             :                                  lambda=lcolv(iconst)%lambda, &
     732        8664 :                                  roll=.TRUE., rmat=v_shake, imass=imass)
     733             :          END DO
     734             :       ELSE
     735      100275 :          DO iconst = 1, SIZE(colv_list)
     736       85854 :             IF (colv_list(iconst)%restraint%active) CYCLE
     737             :             ! Update colvar
     738             :             CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
     739       85854 :                                    pos=pos, fixd_list=fixd_list)
     740             :             lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
     741       85854 :                                               colv_list(iconst)%expected_value)
     742             :             fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar, &
     743             :                                         lcolv(iconst)%colvar_old, roll=.TRUE., rmat=r_shake, &
     744       85854 :                                         imass=imass)
     745       85854 :             del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*dt*fdotf_sum)
     746       85854 :             lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
     747             : 
     748             :             ! Update positions
     749             :             CALL update_con_colv(pos, dtsqby2, lcolv(iconst), &
     750             :                                  lambda=del_lam, &
     751       85854 :                                  roll=.TRUE., rmat=r_shake, imass=imass)
     752             :             ! Update velocities
     753             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     754             :                                  lambda=del_lam, &
     755      100275 :                                  roll=.TRUE., rmat=v_shake, imass=imass)
     756             :          END DO
     757             :       END IF
     758             :       ! computing the constraint and value of tolerance
     759      108939 :       DO iconst = 1, SIZE(colv_list)
     760       93166 :          IF (colv_list(iconst)%restraint%active) CYCLE
     761             :          CALL colvar_eval_mol_f(lcolv(iconst)%colvar, cell, particles=particle_set, &
     762       93166 :                                 pos=pos, fixd_list=fixd_list)
     763             :          lcolv(iconst)%sigma = diff_colvar(lcolv(iconst)%colvar, &
     764       93166 :                                            colv_list(iconst)%expected_value)
     765      108939 :          max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
     766             :       END DO
     767             : 
     768       15773 :    END SUBROUTINE shake_roll_colv_low
     769             : 
     770             : ! **************************************************************************************************
     771             : !> \brief Real Rattle_roll - Low Level
     772             : !>      rattle algorithm (box allowed to change) for collective variables constraints
     773             : !>      updates the multiplier one molecule type at a time
     774             : !> \param fixd_list ...
     775             : !> \param colv_list ...
     776             : !> \param lcolv ...
     777             : !> \param particle_set ...
     778             : !> \param vel ...
     779             : !> \param r_rattle ...
     780             : !> \param dt ...
     781             : !> \param irattle ...
     782             : !> \param veps ...
     783             : !> \param cell ...
     784             : !> \param imass ...
     785             : !> \param max_sigma ...
     786             : !> \par History
     787             : !>      none
     788             : !> \author Teodoro Laino [tlaino]
     789             : ! **************************************************************************************************
     790        5521 :    SUBROUTINE rattle_roll_colv_low(fixd_list, colv_list, lcolv, &
     791       11042 :                                    particle_set, vel, r_rattle, dt, irattle, veps, cell, &
     792        5521 :                                    imass, max_sigma)
     793             : 
     794             :       TYPE(fixd_constraint_type), DIMENSION(:), POINTER  :: fixd_list
     795             :       TYPE(colvar_constraint_type), POINTER              :: colv_list(:)
     796             :       TYPE(local_colvar_constraint_type), POINTER        :: lcolv(:)
     797             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     798             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     799             :       REAL(KIND=dp), INTENT(IN)                          :: r_rattle(:, :), dt
     800             :       INTEGER, INTENT(in)                                :: irattle
     801             :       REAL(KIND=dp), INTENT(IN)                          :: veps(:, :)
     802             :       TYPE(cell_type), POINTER                           :: cell
     803             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     804             :       REAL(KIND=dp), INTENT(INOUT)                       :: max_sigma
     805             : 
     806             :       INTEGER                                            :: iconst
     807             :       REAL(KIND=dp)                                      :: del_lam, dtby2, fdotf_sum
     808             : 
     809        5521 :       dtby2 = dt*.5_dp
     810        5521 :       IF (irattle == 1) THEN
     811        7957 :          DO iconst = 1, SIZE(colv_list)
     812        6670 :             IF (colv_list(iconst)%restraint%active) CYCLE
     813             :             ! Update colvar_old
     814             :             CALL colvar_eval_mol_f(lcolv(iconst)%colvar_old, cell, &
     815        6670 :                                    particles=particle_set, fixd_list=fixd_list)
     816             :             ! Update velocities
     817             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     818             :                                  lambda=lcolv(iconst)%lambda, &
     819        7957 :                                  imass=imass)
     820             :          END DO
     821             :       ELSE
     822       29110 :          DO iconst = 1, SIZE(colv_list)
     823       24876 :             IF (colv_list(iconst)%restraint%active) CYCLE
     824             :             lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel, &
     825       24876 :                                                   roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set)
     826             :             fdotf_sum = eval_Jac_colvar(lcolv(iconst)%colvar_old, &
     827             :                                         lcolv(iconst)%colvar_old, roll=.TRUE., &
     828       24876 :                                         rmat=r_rattle, imass=imass)
     829       24876 :             del_lam = 2.0_dp*lcolv(iconst)%sigma/(dt*fdotf_sum)
     830       24876 :             lcolv(iconst)%lambda = lcolv(iconst)%lambda + del_lam
     831             :             ! Update velocities
     832             :             CALL update_con_colv(vel, dtby2, lcolv(iconst), &
     833             :                                  lambda=del_lam, &
     834       29110 :                                  roll=.TRUE., rmat=r_rattle, imass=imass)
     835             :          END DO
     836             :       END IF
     837             :       ! computing the constraint and value of the tolerance
     838       37067 :       DO iconst = 1, SIZE(colv_list)
     839       31546 :          IF (colv_list(iconst)%restraint%active) CYCLE
     840             :          lcolv(iconst)%sigma = rattle_con_eval(lcolv(iconst)%colvar_old, vel, &
     841       31546 :                                                roll=.TRUE., veps=veps, rmat=r_rattle, particles=particle_set)
     842       37067 :          max_sigma = MAX(ABS(lcolv(iconst)%sigma), max_sigma)
     843             :       END DO
     844             : 
     845        5521 :    END SUBROUTINE rattle_roll_colv_low
     846             : 
     847             : ! **************************************************************************************************
     848             : !> \brief Update position/velocities
     849             : !> \param wrk ...
     850             : !> \param fac ...
     851             : !> \param lcolv ...
     852             : !> \param lambda ...
     853             : !> \param roll ...
     854             : !> \param rmat ...
     855             : !> \param imass ...
     856             : !> \par History
     857             : !>      Teodoro Laino [teo] created 04.2006
     858             : !> \author Teodoro Laino [tlaino]
     859             : ! **************************************************************************************************
     860      444142 :    SUBROUTINE update_con_colv(wrk, fac, lcolv, lambda, roll, rmat, imass)
     861             :       REAL(KIND=dp), INTENT(INOUT)                       :: wrk(:, :)
     862             :       REAL(KIND=dp), INTENT(IN)                          :: fac
     863             :       TYPE(local_colvar_constraint_type), INTENT(IN)     :: lcolv
     864             :       REAL(KIND=dp), INTENT(IN)                          :: lambda
     865             :       LOGICAL, INTENT(in), OPTIONAL                      :: roll
     866             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     867             :          OPTIONAL                                        :: rmat
     868             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     869             : 
     870             :       INTEGER                                            :: iatm, ind
     871             :       LOGICAL                                            :: my_roll
     872             :       REAL(KIND=dp), DIMENSION(3)                        :: f_roll
     873             : 
     874      444142 :       my_roll = .FALSE.
     875      444142 :       IF (PRESENT(roll)) THEN
     876      211208 :          my_roll = roll
     877      211208 :          IF (my_roll) THEN
     878      211208 :             CPASSERT(PRESENT(rmat))
     879             :          END IF
     880             :       END IF
     881     1344082 :       DO iatm = 1, SIZE(lcolv%colvar_old%i_atom)
     882      899940 :          ind = lcolv%colvar_old%i_atom(iatm)
     883             :          !
     884      899940 :          IF (my_roll) THEN
     885             :             ! If ROLL rotate forces
     886     5520684 :             f_roll = MATMUL(rmat, lcolv%colvar_old%dsdr(:, iatm))
     887             :          ELSE
     888     1901088 :             f_roll = lcolv%colvar_old%dsdr(:, iatm)
     889             :          END IF
     890     4043902 :          wrk(:, ind) = wrk(:, ind) - imass(ind)*fac*lambda*f_roll
     891             :       END DO
     892      444142 :    END SUBROUTINE update_con_colv
     893             : 
     894             : ! **************************************************************************************************
     895             : !> \brief Evaluates the Jacobian of the collective variables constraints
     896             : !> \param colvar ...
     897             : !> \param colvar_old ...
     898             : !> \param roll ...
     899             : !> \param rmat ...
     900             : !> \param imass ...
     901             : !> \return ...
     902             : !> \par History
     903             : !>      Teodoro Laino [teo] created 04.2006
     904             : !> \author Teodoro Laino [tlaino]
     905             : ! **************************************************************************************************
     906      229743 :    FUNCTION eval_Jac_colvar(colvar, colvar_old, roll, rmat, imass) RESULT(res)
     907             :       TYPE(colvar_type), POINTER                         :: colvar, colvar_old
     908             :       LOGICAL, INTENT(IN), OPTIONAL                      :: roll
     909             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     910             :          OPTIONAL                                        :: rmat
     911             :       REAL(KIND=dp), DIMENSION(:)                        :: imass
     912             :       REAL(KIND=dp)                                      :: res
     913             : 
     914             :       INTEGER                                            :: i, iatom
     915             :       LOGICAL                                            :: my_roll
     916             :       REAL(KIND=dp), DIMENSION(3)                        :: tmp1, tmp2, tmp3
     917             : 
     918      229743 :       my_roll = .FALSE.
     919      229743 :       IF (PRESENT(roll)) THEN
     920      110730 :          my_roll = roll
     921      110730 :          IF (my_roll) THEN
     922      110730 :             CPASSERT(PRESENT(rmat))
     923             :          END IF
     924             :       END IF
     925             : 
     926      229743 :       res = 0.0_dp
     927      110730 :       IF (my_roll) THEN
     928      332989 :          DO i = 1, SIZE(colvar%i_atom)
     929      222259 :             iatom = colvar%i_atom(i)
     930      889036 :             tmp1 = colvar%dsdr(1:3, i)
     931      889036 :             tmp3 = colvar_old%dsdr(1:3, i)
     932     2889367 :             tmp2 = MATMUL(rmat, tmp3)
     933      999766 :             res = res + DOT_PRODUCT(tmp1, tmp2)*imass(iatom)
     934             :          END DO
     935             :       ELSE
     936      360221 :          DO i = 1, SIZE(colvar%i_atom)
     937      241208 :             iatom = colvar%i_atom(i)
     938      964832 :             tmp1 = colvar%dsdr(1:3, i)
     939      964832 :             tmp2 = colvar_old%dsdr(1:3, i)
     940     1083845 :             res = res + DOT_PRODUCT(tmp1, tmp2)*imass(iatom)
     941             :          END DO
     942             :       END IF
     943             : 
     944      229743 :    END FUNCTION eval_Jac_colvar
     945             : 
     946             : ! **************************************************************************************************
     947             : !> \brief Evaluates the constraint for the rattle scheme
     948             : !> \param colvar ...
     949             : !> \param vel ...
     950             : !> \param roll ...
     951             : !> \param veps ...
     952             : !> \param rmat ...
     953             : !> \param particles ...
     954             : !> \return ...
     955             : !> \par History
     956             : !>      Teodoro Laino [teo] created 04.2006
     957             : !> \author Teodoro Laino [tlaino]
     958             : ! **************************************************************************************************
     959      120949 :    FUNCTION rattle_con_eval(colvar, vel, roll, veps, rmat, particles) RESULT(res)
     960             :       TYPE(colvar_type), POINTER                         :: colvar
     961             :       REAL(KIND=dp), INTENT(INOUT)                       :: vel(:, :)
     962             :       LOGICAL, INTENT(IN), OPTIONAL                      :: roll
     963             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
     964             :          OPTIONAL                                        :: veps, rmat
     965             :       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
     966             :          POINTER                                         :: particles
     967             :       REAL(KIND=dp)                                      :: res
     968             : 
     969             :       INTEGER                                            :: iatm, ind
     970             :       LOGICAL                                            :: my_roll
     971             :       REAL(KIND=dp), DIMENSION(3)                        :: f_roll, pos, v_roll
     972             : 
     973      120949 :       my_roll = .FALSE.
     974      120949 :       IF (PRESENT(roll)) THEN
     975       56422 :          my_roll = roll
     976       56422 :          IF (my_roll) THEN
     977       56422 :             CPASSERT(PRESENT(rmat))
     978       56422 :             CPASSERT(PRESENT(veps))
     979       56422 :             CPASSERT(PRESENT(particles))
     980             :          END IF
     981             :       END IF
     982      120949 :       res = 0.0_dp
     983      368027 :       DO iatm = 1, SIZE(colvar%i_atom)
     984      247078 :          ind = colvar%i_atom(iatm)
     985      247078 :          IF (my_roll) THEN
     986      456848 :             pos = particles(ind)%r
     987     1484756 :             f_roll = MATMUL(rmat, colvar%dsdr(:, iatm))
     988     1827392 :             v_roll = vel(:, ind) + MATMUL(veps, pos)
     989             :          ELSE
     990      531464 :             f_roll = colvar%dsdr(:, iatm)
     991      531464 :             v_roll = vel(:, ind)
     992             :          END IF
     993     1109261 :          res = res + DOT_PRODUCT(f_roll, v_roll)
     994             :       END DO
     995             : 
     996      120949 :    END FUNCTION rattle_con_eval
     997             : 
     998             : END MODULE constraint_clv

Generated by: LCOV version 1.15