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 : 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 1212 : 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 1212 : CALL timeset(routineN, handle)
70 1212 : NULLIFY (Vloc, Vlocg, pw_pool_aux)
71 : !
72 : ! For Martyna-Tuckerman we set up an auxiliary pw_pool at an higher cutoff
73 : !
74 1212 : CALL cite_reference(Martyna1999)
75 1212 : IF (ASSOCIATED(super_ref_pw_grid)) THEN
76 1206 : CALL pw_pool_create(pw_pool_aux, pw_grid=super_ref_pw_grid)
77 : END IF
78 : NULLIFY (screen_function)
79 1212 : ALLOCATE (screen_function)
80 1212 : CALL pw_pool%create_pw(screen_function)
81 1212 : CALL pw_zero(screen_function)
82 2422 : SELECT CASE (method)
83 : CASE (MT0D)
84 1210 : NULLIFY (Vloc, Vlocg)
85 1210 : ALLOCATE (Vloc, Vlocg)
86 1210 : IF (ASSOCIATED(pw_pool_aux)) THEN
87 1204 : CALL pw_pool_aux%create_pw(Vloc)
88 1204 : 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 1210 : CALL mt0din(Vloc, alpha)
94 1210 : CALL pw_transfer(Vloc, Vlocg)
95 1210 : CALL pw_axpy(Vlocg, screen_function)
96 1210 : IF (ASSOCIATED(pw_pool_aux)) THEN
97 1204 : CALL pw_pool_aux%give_back_pw(Vloc)
98 1204 : 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 1210 : DEALLOCATE (Vloc, Vlocg)
104 : !
105 : ! Get rid of the analytical FT of the erf(a*r)/r
106 : !
107 1210 : alpha2 = alpha*alpha
108 81392745 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
109 81391535 : g2 = screen_function%pw_grid%gsq(ig)
110 81391535 : g3d = fourpi/g2
111 81392745 : screen_function%array(ig) = screen_function%array(ig) - g3d*EXP(-g2/(4.0E0_dp*alpha2))
112 : END DO
113 1210 : IF (screen_function%pw_grid%have_g0) &
114 611 : screen_function%array(1) = screen_function%array(1) + fourpi/(4.0E0_dp*alpha2)
115 : CASE (MT2D)
116 2 : iz = special_dimension ! iz is the direction with NO PBC
117 2 : zlength = slab_size ! zlength is the thickness of the cell
118 194401 : DO ig = screen_function%pw_grid%first_gne0, screen_function%pw_grid%ngpts_cut_local
119 194399 : gz = screen_function%pw_grid%g(iz, ig)
120 194399 : g2 = screen_function%pw_grid%gsq(ig)
121 194399 : gxy = SQRT(ABS(g2 - gz*gz))
122 194399 : g3d = fourpi/g2
123 194401 : screen_function%array(ig) = -g3d*COS(gz*zlength/2.0_dp)*EXP(-gxy*zlength/2.0_dp)
124 : END DO
125 2 : IF (screen_function%pw_grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
126 : CASE (MT1D)
127 0 : iz = special_dimension ! iz is the direction with PBC
128 0 : CALL mt1din(screen_function)
129 1212 : CPABORT("MT1D unimplemented")
130 : END SELECT
131 1212 : CALL pw_pool_release(pw_pool_aux)
132 1212 : CALL timestop(handle)
133 1212 : END SUBROUTINE MTin_create_screen_fn
134 :
135 : ! **************************************************************************************************
136 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
137 : !> according the scheme published on:
138 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 110, No. 6, 2810-2821
139 : !> \param Vloc ...
140 : !> \param alpha ...
141 : !> \author Teodoro Laino (09.03.2005)
142 : ! **************************************************************************************************
143 1210 : SUBROUTINE mt0din(Vloc, alpha)
144 : TYPE(pw_r3d_rs_type), POINTER :: Vloc
145 : REAL(KIND=dp), INTENT(in) :: alpha
146 :
147 : CHARACTER(len=*), PARAMETER :: routineN = 'mt0din'
148 :
149 : INTEGER :: handle, i, ii, j, jj, k, kk
150 1210 : INTEGER, DIMENSION(:), POINTER :: glb
151 1210 : INTEGER, DIMENSION(:, :), POINTER :: bo
152 : REAL(KIND=dp) :: dx, dy, dz, fact, omega, r, r2, x, y, &
153 : y2, z, z2
154 : REAL(KIND=dp), DIMENSION(3) :: box, box2
155 : TYPE(pw_grid_type), POINTER :: grid
156 :
157 1210 : CALL timeset(routineN, handle)
158 :
159 1210 : grid => Vloc%pw_grid
160 1210 : bo => grid%bounds_local
161 1210 : glb => grid%bounds(1, :)
162 220740493 : Vloc%array = 0.0_dp
163 4840 : box = REAL(grid%npts, kind=dp)*grid%dr
164 4840 : box2 = box/2.0_dp
165 4840 : omega = PRODUCT(box)
166 1210 : fact = omega
167 1210 : dx = grid%dr(1)
168 1210 : dy = grid%dr(2)
169 1210 : dz = grid%dr(3)
170 1210 : kk = bo(1, 3)
171 79922 : DO k = bo(1, 3), bo(2, 3)
172 78712 : z = REAL(k - glb(3), dp)*dz; IF (z > box2(3)) z = box(3) - z
173 78712 : z2 = z*z
174 78712 : jj = bo(1, 2)
175 5605928 : DO j = bo(1, 2), bo(2, 2)
176 5527216 : y = REAL(j - glb(2), dp)*dy; IF (y > box2(2)) y = box(2) - y
177 5527216 : y2 = y*y
178 5527216 : ii = bo(1, 1)
179 220660571 : DO i = bo(1, 1), bo(2, 1)
180 215133355 : x = REAL(i - glb(1), dp)*dx; IF (x > box2(1)) x = box(1) - x
181 215133355 : r2 = x*x + y2 + z2
182 215133355 : r = SQRT(r2)
183 215133355 : IF (r > 1.0E-10_dp) THEN
184 215132744 : Vloc%array(ii, jj, kk) = erf(alpha*r)/r*fact
185 : ELSE
186 611 : Vloc%array(ii, jj, kk) = 2.0_dp*alpha*oorootpi*fact
187 : END IF
188 220660571 : ii = ii + 1
189 : END DO
190 5605928 : jj = jj + 1
191 : END DO
192 79922 : kk = kk + 1
193 : END DO
194 1210 : CALL timestop(handle)
195 1210 : END SUBROUTINE Mt0din
196 :
197 : ! **************************************************************************************************
198 : !> \brief Calculates the Tuckerman Green's function in reciprocal space
199 : !> according the scheme published on:
200 : !> Martyna and Tuckerman, J. Chem. Phys. Vol. 121, No. 23, 11949
201 : !> \param screen_function ...
202 : !> \author Teodoro Laino (11.2005)
203 : ! **************************************************************************************************
204 0 : SUBROUTINE mt1din(screen_function)
205 : TYPE(pw_c1d_gs_type), POINTER :: screen_function
206 :
207 : CHARACTER(len=*), PARAMETER :: routineN = 'mt1din'
208 :
209 : INTEGER :: handle
210 : REAL(KIND=dp) :: dx, dy, dz, omega
211 : REAL(KIND=dp), DIMENSION(3) :: box, box2
212 : TYPE(pw_grid_type), POINTER :: grid
213 :
214 0 : CALL timeset(routineN, handle)
215 0 : grid => screen_function%pw_grid
216 0 : box = REAL(grid%npts, kind=dp)*grid%dr
217 0 : box2 = box/2.0_dp
218 0 : omega = PRODUCT(box)
219 0 : dx = grid%dr(1)
220 0 : dy = grid%dr(2)
221 0 : dz = grid%dr(3)
222 :
223 0 : CALL timestop(handle)
224 0 : END SUBROUTINE mt1din
225 :
226 : END MODULE mt_util
|
