LCOV - code coverage report
Current view: top level - src/motion - pint_public.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 13.4 % 82 11
Test Date: 2025-07-25 12:55:17 Functions: 20.0 % 5 1

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief  Public path integral routines that can be called from other modules
      10              : !> \author Lukasz Walewski
      11              : !> \date   2009-07-24
      12              : !> \note   Avoiding circular dependencies: please design new members of this
      13              : !>         module in such a way that they use pint_types module only.
      14              : ! **************************************************************************************************
      15              : MODULE pint_public
      16              : 
      17              :    USE kinds,                           ONLY: dp
      18              :    USE parallel_rng_types,              ONLY: rng_stream_type
      19              :    USE pint_types,                      ONLY: pint_env_type
      20              : #include "../base/base_uses.f90"
      21              : 
      22              :    IMPLICIT NONE
      23              : 
      24              :    PRIVATE
      25              : 
      26              :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      27              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pint_public'
      28              : 
      29              :    PUBLIC :: pint_com_pos
      30              :    PUBLIC :: pint_levy_walk
      31              :    PUBLIC :: pint_calc_centroid
      32              : 
      33              : CONTAINS
      34              : 
      35              : ! ***************************************************************************
      36              : !> \brief  Return the center of mass of the PI system
      37              : !> \param pint_env ...
      38              : !> \return ...
      39              : !> \date   2009-07-24
      40              : !> \par    History
      41              : !>           2009-11-30 fixed serious bug in pint_env%x indexing [lwalewski]
      42              : !> \author Lukasz Walewski
      43              : ! **************************************************************************************************
      44          168 :    PURE FUNCTION pint_com_pos(pint_env) RESULT(com_r)
      45              : 
      46              :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      47              :       REAL(kind=dp), DIMENSION(3)                        :: com_r
      48              : 
      49              :       INTEGER                                            :: ia, ib, ic
      50              :       REAL(kind=dp)                                      :: tmass
      51              : 
      52          168 :       tmass = 0.0_dp
      53          672 :       com_r(:) = 0.0_dp
      54          672 :       DO ia = 1, pint_env%ndim/3
      55         4464 :          DO ib = 1, pint_env%p
      56        15672 :             DO ic = 1, 3
      57              :                com_r(ic) = com_r(ic) + &
      58        11376 :                            pint_env%x(ib, (ia - 1)*3 + ic)*pint_env%mass((ia - 1)*3 + ic)
      59        15168 :                tmass = tmass + pint_env%mass((ia - 1)*3 + ic)
      60              :             END DO
      61              :          END DO
      62              :       END DO
      63              :       ! pint_env%mass is REAL, DIMENSION(NDIM) which means that each atom
      64              :       ! has its mass defined three times - here we hope that all three
      65              :       ! values are equal
      66          168 :       tmass = tmass/3.0_dp
      67          672 :       com_r(:) = com_r(:)/tmass
      68          168 :    END FUNCTION pint_com_pos
      69              : 
      70              : ! ***************************************************************************
      71              : !> \brief  Return the center of geometry of the PI system
      72              : !> \param pint_env ...
      73              : !> \return ...
      74              : !> \date   2009-11-30
      75              : !> \author Lukasz Walewski
      76              : ! **************************************************************************************************
      77            0 :    PURE FUNCTION pint_cog_pos(pint_env) RESULT(cntrd_r)
      78              : 
      79              :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
      80              :       REAL(kind=dp), DIMENSION(3)                        :: cntrd_r
      81              : 
      82              :       INTEGER                                            :: ia, ib, ic, natoms
      83              : 
      84            0 :       cntrd_r(:) = 0.0_dp
      85            0 :       natoms = pint_env%ndim/3
      86            0 :       DO ia = 1, natoms
      87            0 :          DO ib = 1, pint_env%p
      88            0 :             DO ic = 1, 3
      89            0 :                cntrd_r(ic) = cntrd_r(ic) + pint_env%x(ib, (ia - 1)*3 + ic)
      90              :             END DO
      91              :          END DO
      92              :       END DO
      93            0 :       cntrd_r(:) = cntrd_r(:)/REAL(pint_env%p, dp)/REAL(natoms, dp)
      94            0 :    END FUNCTION pint_cog_pos
      95              : 
      96              : ! ***************************************************************************
      97              : !> \brief  Given the number of beads n and the variance t returns the
      98              : !>         positions of the beads for a non-interacting particle.
      99              : !> \param n ...
     100              : !> \param t ...
     101              : !> \param rng_gaussian ...
     102              : !> \param x ...
     103              : !> \param nout ...
     104              : !> \date   2010-12-13
     105              : !> \author Lukasz Walewski
     106              : !> \note  This routine implements Levy argorithm (see e.g. Rev. Mod. Phys.
     107              : !>         67 (1995) 279, eq. 5.35) and requires that n is a power of 2. The
     108              : !>         resulting bead positions are centered around (0,0,0).
     109              : ! **************************************************************************************************
     110            0 :    SUBROUTINE pint_free_part_bead_x(n, t, rng_gaussian, x, nout)
     111              : !
     112              : !TODO this routine gives wrong spread of the particles, please fix before usage.
     113              : !
     114              :       INTEGER, INTENT(IN)                                :: n
     115              :       REAL(kind=dp), INTENT(IN)                          :: t
     116              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_gaussian
     117              :       REAL(kind=dp), DIMENSION(:), POINTER               :: x
     118              :       INTEGER, INTENT(OUT)                               :: nout
     119              : 
     120              :       INTEGER                                            :: dl, i1, i2, ib, ic, il, ip, j, nlevels, &
     121              :                                                             np
     122              :       REAL(kind=dp)                                      :: rtmp, tcheck, vrnc, xc
     123              :       REAL(kind=dp), DIMENSION(3)                        :: cntrd_r
     124              : 
     125            0 :       nout = 0
     126              : 
     127            0 :       IF (n < 1) THEN
     128            0 :          RETURN
     129              :       END IF
     130              : 
     131              :       ! if number of beads is not a power of 2 return
     132            0 :       nlevels = NINT(LOG(REAL(n, KIND=dp))/LOG(2.0_dp))
     133            0 :       rtmp = 2**nlevels
     134            0 :       tcheck = ABS(REAL(n, KIND=dp) - rtmp)
     135            0 :       IF (tcheck > 100.0_dp*EPSILON(0.0_dp)) THEN
     136              :          RETURN
     137              :       END IF
     138              : 
     139              :       ! generate at least first point at (0,0,0)
     140            0 :       DO ic = 1, 3
     141            0 :          x(ic) = 0.0_dp
     142              :       END DO
     143            0 :       nout = 1
     144              : 
     145              :       ! loop over Levy levels
     146            0 :       vrnc = 2.0_dp*t
     147            0 :       DO il = 0, nlevels - 1
     148              : 
     149            0 :          np = 2**il ! number of points to be generated at this level
     150            0 :          dl = n/(2*np) ! interval betw points (in index numbers)
     151            0 :          vrnc = vrnc/2.0_dp; ! variance at this level (=t at level 0)
     152              : 
     153              :          ! loop over points added in this level
     154            0 :          DO ip = 0, np - 1
     155              : 
     156            0 :             j = (2*ip + 1)*dl ! index of currently generated point
     157              : 
     158              :             ! indices of two points betw which to generate a new point
     159            0 :             i1 = 2*dl*ip
     160            0 :             i2 = 2*dl*(ip + 1)
     161            0 :             IF (i2 .EQ. n) THEN
     162            0 :                i2 = 0
     163              :             END IF
     164              : 
     165              :             ! generate new point and save it under j
     166            0 :             DO ic = 1, 3
     167            0 :                xc = (x(3*i1 + ic) + x(3*i2 + ic))/2.0
     168            0 :                xc = xc + rng_gaussian%next(variance=vrnc)
     169            0 :                x(3*j + ic) = xc
     170              :             END DO
     171            0 :             nout = nout + 1
     172              : 
     173              :          END DO
     174              :       END DO
     175              : 
     176              :       ! translate the centroid to the origin
     177            0 :       cntrd_r(:) = 0.0_dp
     178            0 :       DO ib = 1, n
     179            0 :          DO ic = 1, 3
     180            0 :             cntrd_r(ic) = cntrd_r(ic) + x((ib - 1)*3 + ic)
     181              :          END DO
     182              :       END DO
     183            0 :       cntrd_r(:) = cntrd_r(:)/REAL(n, dp)
     184            0 :       DO ib = 1, n
     185            0 :          DO ic = 1, 3
     186            0 :             x((ib - 1)*3 + ic) = x((ib - 1)*3 + ic) - cntrd_r(ic)
     187              :          END DO
     188              :       END DO
     189              : 
     190              :    END SUBROUTINE pint_free_part_bead_x
     191              : 
     192              : ! ***************************************************************************
     193              : !> \brief  Perform a Brownian walk of length n around x0 with the variance v.
     194              : !> \param x0 ...
     195              : !> \param n ...
     196              : !> \param v ...
     197              : !> \param x ...
     198              : !> \param rng_gaussian ...
     199              : !> \date   2011-01-06
     200              : !> \author Lukasz Walewski
     201              : !> \note  This routine implements Levy argorithm (Phys. Rev. 143 (1966) 58)
     202              : ! **************************************************************************************************
     203            0 :    SUBROUTINE pint_levy_walk(x0, n, v, x, rng_gaussian)
     204              : 
     205              :       REAL(kind=dp), DIMENSION(3), INTENT(IN)            :: x0
     206              :       INTEGER, INTENT(IN)                                :: n
     207              :       REAL(kind=dp), INTENT(IN)                          :: v
     208              :       REAL(kind=dp), DIMENSION(:), POINTER               :: x
     209              :       TYPE(rng_stream_type), INTENT(INOUT)               :: rng_gaussian
     210              : 
     211              :       INTEGER                                            :: ib, ic
     212              :       REAL(kind=dp)                                      :: r, tau_i, tau_i1
     213              :       REAL(kind=dp), DIMENSION(3)                        :: cntrd_r
     214              : 
     215            0 :       x(1) = x0(1)
     216            0 :       x(2) = x0(2)
     217            0 :       x(3) = x0(3)
     218            0 :       DO ib = 1, n - 1
     219            0 :          DO ic = 1, 3
     220            0 :             r = rng_gaussian%next(variance=1.0_dp)
     221            0 :             tau_i = (REAL(ib, dp) - 1.0_dp)/REAL(n, dp)
     222            0 :             tau_i1 = (REAL(ib + 1, dp) - 1.0_dp)/REAL(n, dp)
     223              :             x(ib*3 + ic) = (x((ib - 1)*3 + ic)*(1.0_dp - tau_i1) + &
     224              :                             x(ic)*(tau_i1 - tau_i))/ &
     225              :                            (1.0_dp - tau_i) + &
     226              :                            r*v*SQRT( &
     227              :                            (tau_i1 - tau_i)* &
     228              :                            (1.0_dp - tau_i1)/ &
     229              :                            (1.0_dp - tau_i) &
     230            0 :                            )
     231              :          END DO
     232              :       END DO
     233              : 
     234              :       ! translate the centroid to the origin
     235            0 :       cntrd_r(:) = 0.0_dp
     236            0 :       DO ib = 1, n
     237            0 :          DO ic = 1, 3
     238            0 :             cntrd_r(ic) = cntrd_r(ic) + x((ib - 1)*3 + ic)
     239              :          END DO
     240              :       END DO
     241            0 :       cntrd_r(:) = cntrd_r(:)/REAL(n, dp)
     242            0 :       DO ib = 1, n
     243            0 :          DO ic = 1, 3
     244            0 :             x((ib - 1)*3 + ic) = x((ib - 1)*3 + ic) - cntrd_r(ic)
     245              :          END DO
     246              :       END DO
     247              : 
     248            0 :    END SUBROUTINE pint_levy_walk
     249              : 
     250              : ! ***************************************************************************
     251              : !> \brief  Calculate the centroid
     252              : !> \param  pint_env path integral environment
     253              : !> \date   2013-02-10
     254              : !> \author lwalewski
     255              : ! **************************************************************************************************
     256            0 :    PURE SUBROUTINE pint_calc_centroid(pint_env)
     257              : 
     258              :       TYPE(pint_env_type), INTENT(INOUT)                 :: pint_env
     259              : 
     260              :       INTEGER                                            :: ia, ib
     261              :       REAL(KIND=dp)                                      :: invp
     262              : 
     263            0 :       invp = 1.0_dp/pint_env%p
     264            0 :       pint_env%centroid(:) = 0.0_dp
     265            0 :       DO ia = 1, pint_env%ndim
     266            0 :          DO ib = 1, pint_env%p
     267            0 :             pint_env%centroid(ia) = pint_env%centroid(ia) + pint_env%x(ib, ia)
     268              :          END DO
     269            0 :          pint_env%centroid(ia) = pint_env%centroid(ia)*invp
     270              :       END DO
     271              : 
     272            0 :    END SUBROUTINE pint_calc_centroid
     273              : 
     274              : END MODULE pint_public
        

Generated by: LCOV version 2.0-1