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
|
