Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE mt_util
10 : USE bibliography, ONLY: Martyna1999,&
11 : cite_reference
12 : USE kinds, ONLY: dp
13 : USE mathconstants, ONLY: fourpi,&
14 : oorootpi,&
15 : pi
16 : USE pw_grid_types, ONLY: pw_grid_type
17 : USE pw_methods, ONLY: pw_axpy,&
18 : pw_transfer,&
19 : pw_zero
20 : USE pw_pool_types, ONLY: pw_pool_create,&
21 : pw_pool_release,&
22 : pw_pool_type
23 : USE pw_types, ONLY: pw_c1d_gs_type,&
24 : pw_r3d_rs_type
25 : #include "../base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 :
29 : PRIVATE
30 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mt_util'
32 :
33 : INTEGER, PARAMETER, PUBLIC :: MT2D = 1101, &
34 : MT1D = 1102, &
35 : MT0D = 1103
36 :
37 : PUBLIC :: MTin_create_screen_fn, mt2d_kernel_strain
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief MT2D Coulomb kernel and its affine strain derivative in the periodic plane.
42 : !> \param g Cartesian reciprocal vector, transformed by the inverse transpose of the strain.
43 : !> \param special_dimension Cartesian direction normal to the slab.
44 : !> \param slab_size Cell thickness in the nonperiodic direction, fixed for in-plane strain.
45 : !> \param kernel Reciprocal Coulomb kernel, including the existing MT2D G=0 convention.
46 : !> \param strain Kernel derivative with respect to in-plane affine strain.
47 : ! **************************************************************************************************
48 378918 : SUBROUTINE mt2d_kernel_strain(g, special_dimension, slab_size, kernel, strain)
49 : REAL(KIND=dp), INTENT(IN) :: g(3)
50 : INTEGER, INTENT(IN) :: special_dimension
51 : REAL(KIND=dp), INTENT(IN) :: slab_size
52 : REAL(KIND=dp), INTENT(OUT) :: kernel, strain(3, 3)
53 :
54 : INTEGER :: alpha, beta
55 : REAL(KIND=dp) :: f, g2, k, parallel(3), q, radius
56 :
57 378918 : CPASSERT(special_dimension >= 1 .AND. special_dimension <= 3)
58 378918 : CPASSERT(slab_size > 0.0_dp)
59 378918 : strain = 0.0_dp
60 1515672 : g2 = SUM(g**2)
61 378918 : IF (g2 == 0.0_dp) THEN
62 : ! This kernel is independent of in-plane strain. Its volume term is
63 : ! already included in the Hartree energy contribution to the virial.
64 78 : kernel = pi*slab_size**2/2.0_dp
65 906 : RETURN
66 : END IF
67 378840 : parallel = g
68 378840 : parallel(special_dimension) = 0.0_dp
69 1515360 : k = NORM2(parallel)
70 378840 : radius = slab_size/2.0_dp
71 378840 : q = COS(g(special_dimension)*radius)*EXP(-k*radius)
72 378840 : f = 1.0_dp - q
73 378840 : kernel = fourpi/g2*f
74 : ! For k=0 the reciprocal vector is unchanged by every in-plane strain.
75 : ! This also avoids dividing by k at the even-normal modes where K=0.
76 378840 : IF (k == 0.0_dp) RETURN
77 1512048 : DO alpha = 1, 3
78 1134036 : IF (alpha == special_dimension) CYCLE
79 3024084 : DO beta = alpha, 3
80 1890048 : IF (beta == special_dimension) CYCLE
81 1134036 : strain(alpha, beta) = fourpi/g2*g(alpha)*g(beta)*(2.0_dp*f/g2 - q*radius/k)
82 3024084 : strain(beta, alpha) = strain(alpha, beta)
83 : END DO
84 : END DO
85 : END SUBROUTINE mt2d_kernel_strain
86 :
87 : ! **************************************************************************************************
88 : !> \brief Initialize the Martyna && Tuckerman Poisson Solver
89 : !> \param screen_function ...
90 : !> \param pw_pool ...
91 : !> \param method ...
92 : !> \param alpha ...
93 : !> \param special_dimension ...
94 : !> \param slab_size ...
95 : !> \param super_ref_pw_grid ...
96 : !> \author Teodoro Laino (16.06.2004)
97 : ! **************************************************************************************************
98 1630 : SUBROUTINE MTin_create_screen_fn(screen_function, pw_pool, method, alpha, &
99 : special_dimension, slab_size, super_ref_pw_grid)
100 : TYPE(pw_c1d_gs_type), POINTER :: screen_function
101 : TYPE(pw_pool_type), POINTER :: pw_pool
102 : INTEGER, INTENT(IN) :: method
103 : REAL(KIND=dp), INTENT(in) :: alpha
104 : INTEGER, INTENT(IN) :: special_dimension
105 : REAL(KIND=dp), INTENT(in) :: slab_size
106 : TYPE(pw_grid_type), POINTER :: super_ref_pw_grid
107 :
108 : CHARACTER(len=*), PARAMETER :: routineN = 'MTin_create_screen_fn'
109 :
110 : INTEGER :: handle, ig, iz
111 : REAL(KIND=dp) :: alpha2, g2, g3d, gxy, gz, zlength
112 : TYPE(pw_c1d_gs_type), POINTER :: Vlocg
113 : TYPE(pw_pool_type), POINTER :: pw_pool_aux
114 : TYPE(pw_r3d_rs_type), POINTER :: Vloc
115 :
116 1630 : CALL timeset(routineN, handle)
117 1630 : NULLIFY (Vloc, Vlocg, pw_pool_aux)
118 : !
119 : ! For Martyna-Tuckerman we set up an auxiliary pw_pool at an higher cutoff
120 : !
121 1630 : CALL cite_reference(Martyna1999)
122 1630 : IF (ASSOCIATED(super_ref_pw_grid)) THEN
123 1620 : CALL pw_pool_create(pw_pool_aux, pw_grid=super_ref_pw_grid)
124 : END IF
125 : NULLIFY (screen_function)
126 1630 : ALLOCATE (screen_function)
127 1630 : CALL pw_pool%create_pw(screen_function)
128 1630 : CALL pw_zero(screen_function)
129 3190 : SELECT CASE (method)
130 : CASE (MT0D)
131 1560 : NULLIFY (Vloc, Vlocg)
132 1560 : ALLOCATE (Vloc, Vlocg)
133 1560 : IF (ASSOCIATED(pw_pool_aux)) THEN
134 1554 : CALL pw_pool_aux%create_pw(Vloc)
135 1554 : CALL pw_pool_aux%create_pw(Vlocg)
136 : ELSE
137 6 : CALL pw_pool%create_pw(Vloc)
138 6 : CALL pw_pool%create_pw(Vlocg)
139 : END IF
140 1560 : CALL mt0din(Vloc, alpha)
141 1560 : CALL pw_transfer(Vloc, Vlocg)
142 1560 : CALL pw_axpy(Vlocg, screen_function)
143 1560 : IF (ASSOCIATED(pw_pool_aux)) THEN
144 1554 : CALL pw_pool_aux%give_back_pw(Vloc)
145 1554 : CALL pw_pool_aux%give_back_pw(Vlocg)
146 : ELSE
147 6 : CALL pw_pool%give_back_pw(Vloc)
148 6 : CALL pw_pool%give_back_pw(Vlocg)
149 : END IF
150 1560 : DEALLOCATE (Vloc, Vlocg)
151 : !
152 : ! Get rid of the analytical FT of the erf(a*r)/r
153 : !
154 1560 : alpha2 = alpha*alpha
155 116994904 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
156 116993344 : g2 = screen_function%pw_grid%gsq(ig)
157 116993344 : g3d = fourpi/g2
158 116994904 : screen_function%array(ig) = screen_function%array(ig) - g3d*EXP(-g2/(4.0E0_dp*alpha2))
159 : END DO
160 1560 : IF (screen_function%pw_grid%have_g0) THEN
161 786 : screen_function%array(1) = screen_function%array(1) + fourpi/(4.0E0_dp*alpha2)
162 : END IF
163 : CASE (MT2D)
164 70 : iz = special_dimension ! iz is the direction with NO PBC
165 70 : zlength = slab_size ! zlength is the thickness of the cell
166 1189951 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
167 1189881 : gz = screen_function%pw_grid%g(iz, ig)
168 1189881 : g2 = screen_function%pw_grid%gsq(ig)
169 1189881 : gxy = SQRT(ABS(g2 - gz*gz))
170 1189881 : g3d = fourpi/g2
171 1189951 : screen_function%array(ig) = -g3d*COS(gz*zlength/2.0_dp)*EXP(-gxy*zlength/2.0_dp)
172 : END DO
173 70 : IF (screen_function%pw_grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
174 : CASE (MT1D)
175 0 : iz = special_dimension ! iz is the direction with PBC
176 0 : CALL mt1din(screen_function)
177 1630 : CPABORT("MT1D unimplemented")
178 : END SELECT
179 1630 : CALL pw_pool_release(pw_pool_aux)
180 1630 : CALL timestop(handle)
181 1630 : END SUBROUTINE MTin_create_screen_fn
182 :
183 : ! **************************************************************************************************
184 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
185 : !> according the scheme published on:
186 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 110, No. 6, 2810-2821
187 : !> \param Vloc ...
188 : !> \param alpha ...
189 : !> \author Teodoro Laino (09.03.2005)
190 : ! **************************************************************************************************
191 1560 : SUBROUTINE mt0din(Vloc, alpha)
192 : TYPE(pw_r3d_rs_type), POINTER :: Vloc
193 : REAL(KIND=dp), INTENT(in) :: alpha
194 :
195 : CHARACTER(len=*), PARAMETER :: routineN = 'mt0din'
196 :
197 : INTEGER :: handle, i, ii, j, jj, k, kk
198 1560 : INTEGER, DIMENSION(:), POINTER :: glb
199 1560 : INTEGER, DIMENSION(:, :), POINTER :: bo
200 : REAL(KIND=dp) :: dx, dy, dz, fact, omega, r, r2, x, y, &
201 : y2, z, z2
202 : REAL(KIND=dp), DIMENSION(3) :: box, box2
203 : TYPE(pw_grid_type), POINTER :: grid
204 :
205 1560 : CALL timeset(routineN, handle)
206 :
207 1560 : grid => Vloc%pw_grid
208 1560 : bo => grid%bounds_local
209 1560 : glb => grid%bounds(1, :)
210 324165498 : Vloc%array = 0.0_dp
211 6240 : box = REAL(grid%npts, kind=dp)*grid%dr
212 6240 : box2 = box/2.0_dp
213 6240 : omega = PRODUCT(box)
214 1560 : fact = omega
215 1560 : dx = grid%dr(1)
216 1560 : dy = grid%dr(2)
217 1560 : dz = grid%dr(3)
218 1560 : kk = bo(1, 3)
219 107262 : DO k = bo(1, 3), bo(2, 3)
220 105702 : z = REAL(k - glb(3), dp)*dz; IF (z > box2(3)) z = box(3) - z
221 105702 : z2 = z*z
222 105702 : jj = bo(1, 2)
223 7885588 : DO j = bo(1, 2), bo(2, 2)
224 7779886 : y = REAL(j - glb(2), dp)*dy; IF (y > box2(2)) y = box(2) - y
225 7779886 : y2 = y*y
226 7779886 : ii = bo(1, 1)
227 324058236 : DO i = bo(1, 1), bo(2, 1)
228 316278350 : x = REAL(i - glb(1), dp)*dx; IF (x > box2(1)) x = box(1) - x
229 316278350 : r2 = x*x + y2 + z2
230 316278350 : r = SQRT(r2)
231 316278350 : IF (r > 1.0E-10_dp) THEN
232 316277564 : Vloc%array(ii, jj, kk) = erf(alpha*r)/r*fact
233 : ELSE
234 786 : Vloc%array(ii, jj, kk) = 2.0_dp*alpha*oorootpi*fact
235 : END IF
236 324058236 : ii = ii + 1
237 : END DO
238 7885588 : jj = jj + 1
239 : END DO
240 107262 : kk = kk + 1
241 : END DO
242 1560 : CALL timestop(handle)
243 1560 : END SUBROUTINE Mt0din
244 :
245 : ! **************************************************************************************************
246 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
247 : !> according the scheme published on:
248 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 121, No. 23, 11949
249 : !> \param screen_function ...
250 : !> \author Teodoro Laino (11.2005)
251 : ! **************************************************************************************************
252 0 : SUBROUTINE mt1din(screen_function)
253 : TYPE(pw_c1d_gs_type), POINTER :: screen_function
254 :
255 : CHARACTER(len=*), PARAMETER :: routineN = 'mt1din'
256 :
257 : INTEGER :: handle
258 : REAL(KIND=dp) :: dx, dy, dz, omega
259 : REAL(KIND=dp), DIMENSION(3) :: box, box2
260 : TYPE(pw_grid_type), POINTER :: grid
261 :
262 0 : CALL timeset(routineN, handle)
263 0 : grid => screen_function%pw_grid
264 0 : box = REAL(grid%npts, kind=dp)*grid%dr
265 0 : box2 = box/2.0_dp
266 0 : omega = PRODUCT(box)
267 0 : dx = grid%dr(1)
268 0 : dy = grid%dr(2)
269 0 : dz = grid%dr(3)
270 :
271 0 : CALL timestop(handle)
272 0 : END SUBROUTINE mt1din
273 :
274 : END MODULE mt_util
|
