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
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Initialize the Martyna && Tuckerman Poisson Solver
42 : !> \param screen_function ...
43 : !> \param pw_pool ...
44 : !> \param method ...
45 : !> \param alpha ...
46 : !> \param special_dimension ...
47 : !> \param slab_size ...
48 : !> \param super_ref_pw_grid ...
49 : !> \author Teodoro Laino (16.06.2004)
50 : ! **************************************************************************************************
51 1542 : SUBROUTINE MTin_create_screen_fn(screen_function, pw_pool, method, alpha, &
52 : special_dimension, slab_size, super_ref_pw_grid)
53 : TYPE(pw_c1d_gs_type), POINTER :: screen_function
54 : TYPE(pw_pool_type), POINTER :: pw_pool
55 : INTEGER, INTENT(IN) :: method
56 : REAL(KIND=dp), INTENT(in) :: alpha
57 : INTEGER, INTENT(IN) :: special_dimension
58 : REAL(KIND=dp), INTENT(in) :: slab_size
59 : TYPE(pw_grid_type), POINTER :: super_ref_pw_grid
60 :
61 : CHARACTER(len=*), PARAMETER :: routineN = 'MTin_create_screen_fn'
62 :
63 : INTEGER :: handle, ig, iz
64 : REAL(KIND=dp) :: alpha2, g2, g3d, gxy, gz, zlength
65 : TYPE(pw_c1d_gs_type), POINTER :: Vlocg
66 : TYPE(pw_pool_type), POINTER :: pw_pool_aux
67 : TYPE(pw_r3d_rs_type), POINTER :: Vloc
68 :
69 1542 : CALL timeset(routineN, handle)
70 1542 : NULLIFY (Vloc, Vlocg, pw_pool_aux)
71 : !
72 : ! For Martyna-Tuckerman we set up an auxiliary pw_pool at an higher cutoff
73 : !
74 1542 : CALL cite_reference(Martyna1999)
75 1542 : IF (ASSOCIATED(super_ref_pw_grid)) THEN
76 1536 : CALL pw_pool_create(pw_pool_aux, pw_grid=super_ref_pw_grid)
77 : END IF
78 : NULLIFY (screen_function)
79 1542 : ALLOCATE (screen_function)
80 1542 : CALL pw_pool%create_pw(screen_function)
81 1542 : CALL pw_zero(screen_function)
82 3018 : SELECT CASE (method)
83 : CASE (MT0D)
84 1476 : NULLIFY (Vloc, Vlocg)
85 1476 : ALLOCATE (Vloc, Vlocg)
86 1476 : IF (ASSOCIATED(pw_pool_aux)) THEN
87 1470 : CALL pw_pool_aux%create_pw(Vloc)
88 1470 : CALL pw_pool_aux%create_pw(Vlocg)
89 : ELSE
90 6 : CALL pw_pool%create_pw(Vloc)
91 6 : CALL pw_pool%create_pw(Vlocg)
92 : END IF
93 1476 : CALL mt0din(Vloc, alpha)
94 1476 : CALL pw_transfer(Vloc, Vlocg)
95 1476 : CALL pw_axpy(Vlocg, screen_function)
96 1476 : IF (ASSOCIATED(pw_pool_aux)) THEN
97 1470 : CALL pw_pool_aux%give_back_pw(Vloc)
98 1470 : CALL pw_pool_aux%give_back_pw(Vlocg)
99 : ELSE
100 6 : CALL pw_pool%give_back_pw(Vloc)
101 6 : CALL pw_pool%give_back_pw(Vlocg)
102 : END IF
103 1476 : DEALLOCATE (Vloc, Vlocg)
104 : !
105 : ! Get rid of the analytical FT of the erf(a*r)/r
106 : !
107 1476 : alpha2 = alpha*alpha
108 109551449 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
109 109549973 : g2 = screen_function%pw_grid%gsq(ig)
110 109549973 : g3d = fourpi/g2
111 109551449 : screen_function%array(ig) = screen_function%array(ig) - g3d*EXP(-g2/(4.0E0_dp*alpha2))
112 : END DO
113 1476 : IF (screen_function%pw_grid%have_g0) THEN
114 744 : screen_function%array(1) = screen_function%array(1) + fourpi/(4.0E0_dp*alpha2)
115 : END IF
116 : CASE (MT2D)
117 66 : iz = special_dimension ! iz is the direction with NO PBC
118 66 : zlength = slab_size ! zlength is the thickness of the cell
119 1189761 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
120 1189695 : gz = screen_function%pw_grid%g(iz, ig)
121 1189695 : g2 = screen_function%pw_grid%gsq(ig)
122 1189695 : gxy = SQRT(ABS(g2 - gz*gz))
123 1189695 : g3d = fourpi/g2
124 1189761 : screen_function%array(ig) = -g3d*COS(gz*zlength/2.0_dp)*EXP(-gxy*zlength/2.0_dp)
125 : END DO
126 66 : IF (screen_function%pw_grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
127 : CASE (MT1D)
128 0 : iz = special_dimension ! iz is the direction with PBC
129 0 : CALL mt1din(screen_function)
130 1542 : CPABORT("MT1D unimplemented")
131 : END SELECT
132 1542 : CALL pw_pool_release(pw_pool_aux)
133 1542 : CALL timestop(handle)
134 1542 : END SUBROUTINE MTin_create_screen_fn
135 :
136 : ! **************************************************************************************************
137 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
138 : !> according the scheme published on:
139 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 110, No. 6, 2810-2821
140 : !> \param Vloc ...
141 : !> \param alpha ...
142 : !> \author Teodoro Laino (09.03.2005)
143 : ! **************************************************************************************************
144 1476 : SUBROUTINE mt0din(Vloc, alpha)
145 : TYPE(pw_r3d_rs_type), POINTER :: Vloc
146 : REAL(KIND=dp), INTENT(in) :: alpha
147 :
148 : CHARACTER(len=*), PARAMETER :: routineN = 'mt0din'
149 :
150 : INTEGER :: handle, i, ii, j, jj, k, kk
151 1476 : INTEGER, DIMENSION(:), POINTER :: glb
152 1476 : INTEGER, DIMENSION(:, :), POINTER :: bo
153 : REAL(KIND=dp) :: dx, dy, dz, fact, omega, r, r2, x, y, &
154 : y2, z, z2
155 : REAL(KIND=dp), DIMENSION(3) :: box, box2
156 : TYPE(pw_grid_type), POINTER :: grid
157 :
158 1476 : CALL timeset(routineN, handle)
159 :
160 1476 : grid => Vloc%pw_grid
161 1476 : bo => grid%bounds_local
162 1476 : glb => grid%bounds(1, :)
163 303849202 : Vloc%array = 0.0_dp
164 5904 : box = REAL(grid%npts, kind=dp)*grid%dr
165 5904 : box2 = box/2.0_dp
166 5904 : omega = PRODUCT(box)
167 1476 : fact = omega
168 1476 : dx = grid%dr(1)
169 1476 : dy = grid%dr(2)
170 1476 : dz = grid%dr(3)
171 1476 : kk = bo(1, 3)
172 100938 : DO k = bo(1, 3), bo(2, 3)
173 99462 : z = REAL(k - glb(3), dp)*dz; IF (z > box2(3)) z = box(3) - z
174 99462 : z2 = z*z
175 99462 : jj = bo(1, 2)
176 7396784 : DO j = bo(1, 2), bo(2, 2)
177 7297322 : y = REAL(j - glb(2), dp)*dy; IF (y > box2(2)) y = box(2) - y
178 7297322 : y2 = y*y
179 7297322 : ii = bo(1, 1)
180 303748264 : DO i = bo(1, 1), bo(2, 1)
181 296450942 : x = REAL(i - glb(1), dp)*dx; IF (x > box2(1)) x = box(1) - x
182 296450942 : r2 = x*x + y2 + z2
183 296450942 : r = SQRT(r2)
184 296450942 : IF (r > 1.0E-10_dp) THEN
185 296450198 : Vloc%array(ii, jj, kk) = erf(alpha*r)/r*fact
186 : ELSE
187 744 : Vloc%array(ii, jj, kk) = 2.0_dp*alpha*oorootpi*fact
188 : END IF
189 303748264 : ii = ii + 1
190 : END DO
191 7396784 : jj = jj + 1
192 : END DO
193 100938 : kk = kk + 1
194 : END DO
195 1476 : CALL timestop(handle)
196 1476 : END SUBROUTINE Mt0din
197 :
198 : ! **************************************************************************************************
199 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
200 : !> according the scheme published on:
201 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 121, No. 23, 11949
202 : !> \param screen_function ...
203 : !> \author Teodoro Laino (11.2005)
204 : ! **************************************************************************************************
205 0 : SUBROUTINE mt1din(screen_function)
206 : TYPE(pw_c1d_gs_type), POINTER :: screen_function
207 :
208 : CHARACTER(len=*), PARAMETER :: routineN = 'mt1din'
209 :
210 : INTEGER :: handle
211 : REAL(KIND=dp) :: dx, dy, dz, omega
212 : REAL(KIND=dp), DIMENSION(3) :: box, box2
213 : TYPE(pw_grid_type), POINTER :: grid
214 :
215 0 : CALL timeset(routineN, handle)
216 0 : grid => screen_function%pw_grid
217 0 : box = REAL(grid%npts, kind=dp)*grid%dr
218 0 : box2 = box/2.0_dp
219 0 : omega = PRODUCT(box)
220 0 : dx = grid%dr(1)
221 0 : dy = grid%dr(2)
222 0 : dz = grid%dr(3)
223 :
224 0 : CALL timestop(handle)
225 0 : END SUBROUTINE mt1din
226 :
227 : END MODULE mt_util
|
