LCOV - code coverage report
Current view: top level - src/motion - pint_gle.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:f856083) Lines: 61 61 100.0 %
Date: 2024-11-01 07:00:00 Functions: 3 3 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  Methods to apply GLE to PI runs.
      10             : !> \author michelec
      11             : !> \par History
      12             : !>      06.2010 created [michelec]
      13             : !> \note   trying to keep duplication at a minimum....
      14             : ! **************************************************************************************************
      15             : 
      16             : MODULE pint_gle
      17             :    USE gle_system_dynamics,             ONLY: gle_cholesky_stab
      18             :    USE gle_system_types,                ONLY: gle_type
      19             :    USE kinds,                           ONLY: dp
      20             :    USE pint_types,                      ONLY: pint_env_type
      21             : #include "../base/base_uses.f90"
      22             : 
      23             :    IMPLICIT NONE
      24             : 
      25             :    PRIVATE
      26             : 
      27             :    PUBLIC :: pint_gle_step, pint_gle_init, pint_calc_gle_energy
      28             : 
      29             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_gle'
      30             : 
      31             : CONTAINS
      32             : ! **************************************************************************************************
      33             : !> \brief ...
      34             : !> \param pint_env ...
      35             : ! **************************************************************************************************
      36           6 :    ELEMENTAL SUBROUTINE pint_calc_gle_energy(pint_env)
      37             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
      38             : 
      39             :       INTEGER                                            :: i
      40             : 
      41           6 :       pint_env%e_gle = 0._dp
      42           6 :       IF (ASSOCIATED(pint_env%gle)) THEN
      43       55302 :          DO i = 1, pint_env%gle%loc_num_gle
      44       55302 :             pint_env%e_gle = pint_env%e_gle + pint_env%gle%nvt(i)%thermostat_energy
      45             :          END DO
      46             :       END IF
      47           6 :    END SUBROUTINE
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief ...
      51             : !> \param pint_env ...
      52             : ! **************************************************************************************************
      53           2 :    SUBROUTINE pint_gle_init(pint_env)
      54             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
      55             : 
      56             :       INTEGER                                            :: i, ib, idim, imap, j
      57           4 :       REAL(dp) :: mf, rr(pint_env%gle%ndim), cc(pint_env%gle%ndim, pint_env%gle%ndim)
      58             : 
      59           2 :       CALL gle_cholesky_stab(pint_env%gle%c_mat, cc, pint_env%gle%ndim)
      60       18434 :       DO i = 1, pint_env%gle%loc_num_gle
      61       18432 :          imap = pint_env%gle%map_info%index(i)
      62       18432 :          ib = 1 + (imap - 1)/pint_env%ndim
      63       18432 :          idim = 1 + MOD(imap - 1, pint_env%ndim)
      64       18432 :          mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
      65      110592 :          DO j = 1, pint_env%gle%ndim
      66      110592 :             rr(j) = pint_env%gle%nvt(i)%gaussian_rng_stream%next()*mf
      67             :          END DO
      68      663554 :          pint_env%gle%nvt(i)%s = MATMUL(cc, rr)
      69             :       END DO
      70             : 
      71           2 :    END SUBROUTINE pint_gle_init
      72             : 
      73             : ! **************************************************************************************************
      74             : !> \brief ...
      75             : !> \param pint_env ...
      76             : ! **************************************************************************************************
      77           4 :    SUBROUTINE pint_gle_step(pint_env)
      78             :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
      79             : 
      80             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pint_gle_step'
      81             : 
      82             :       INTEGER                                            :: handle, iadd, ib, ideg, idim, imap, &
      83             :                                                             ndim, num
      84             :       REAL(dp)                                           :: alpha, beta, mf, rr
      85           4 :       REAL(dp), DIMENSION(:, :), POINTER                 :: a_mat, e_tmp, h_tmp, s_tmp
      86             :       TYPE(gle_type), POINTER                            :: gle
      87             : 
      88           4 :       CALL timeset(routineN, handle)
      89             : 
      90           4 :       gle => pint_env%gle
      91           4 :       ndim = gle%ndim
      92             : 
      93          16 :       ALLOCATE (s_tmp(ndim, gle%loc_num_gle))
      94      221188 :       s_tmp = 0.0_dp
      95          12 :       ALLOCATE (e_tmp(ndim, gle%loc_num_gle))
      96          12 :       ALLOCATE (h_tmp(ndim, gle%loc_num_gle))
      97             : 
      98       36868 :       DO ideg = 1, gle%loc_num_gle
      99       36864 :          imap = gle%map_info%index(ideg)
     100       36864 :          ib = 1 + (imap - 1)/pint_env%ndim
     101       36864 :          idim = 1 + MOD(imap - 1, pint_env%ndim)
     102             : 
     103       36864 :          gle%nvt(ideg)%s(1) = pint_env%uv_t(ib, idim)
     104             :          gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy &
     105       36864 :                                            + 0.5_dp*pint_env%mass_fict(ib, idim)*gle%nvt(ideg)%s(1)**2
     106       36864 :          s_tmp(1, imap) = gle%nvt(ideg)%s(1)
     107       36864 :          rr = gle%nvt(ideg)%gaussian_rng_stream%next()
     108       36864 :          mf = 1.0_dp/SQRT(pint_env%mass_fict(ib, idim))
     109       36864 :          e_tmp(1, imap) = rr*mf
     110      184324 :          DO iadd = 2, ndim
     111      147456 :             s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd)
     112      147456 :             rr = gle%nvt(ideg)%gaussian_rng_stream%next()
     113      184320 :             e_tmp(iadd, imap) = rr*mf
     114             :          END DO
     115             :       END DO
     116           4 :       num = gle%loc_num_gle
     117           4 :       a_mat => gle%gle_s
     118           4 :       alpha = 1.0_dp
     119           4 :       beta = 0.0_dp
     120             : 
     121           4 :       CALL DGEMM('N', 'N', ndim, num, ndim, alpha, a_mat(1, 1), ndim, e_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
     122             : 
     123           4 :       a_mat => gle%gle_t
     124           4 :       beta = 1.0_dp
     125           4 :       CALL DGEMM("N", "N", ndim, num, ndim, alpha, a_mat(1, 1), ndim, s_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim)
     126             : 
     127       36868 :       DO ideg = 1, gle%loc_num_gle
     128       36864 :          imap = gle%map_info%index(ideg)
     129             : 
     130      221184 :          DO iadd = 1, ndim
     131      221184 :             gle%nvt(ideg)%s(iadd) = h_tmp(iadd, imap)
     132             :          END DO
     133       36864 :          ib = 1 + (imap - 1)/pint_env%ndim
     134       36864 :          idim = 1 + MOD(imap - 1, pint_env%ndim)
     135       36864 :          pint_env%uv_t(ib, idim) = gle%nvt(ideg)%s(1)
     136             :          gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy &
     137       36868 :                                            - 0.5_dp*pint_env%mass_fict(ib, idim)*gle%nvt(ideg)%s(1)**2
     138             :       END DO
     139           4 :       pint_env%e_kin_t = 0.0_dp
     140           4 :       DEALLOCATE (e_tmp, s_tmp, h_tmp)
     141           4 :       CALL timestop(handle)
     142           4 :    END SUBROUTINE
     143             : END MODULE

Generated by: LCOV version 1.15