LCOV - code coverage report
Current view: top level - src/motion - pint_public.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 11 82 13.4 %
Date: 2024-04-20 06:29:22 Functions: 1 5 20.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  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 1.15