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 contains utility functions for the xc package
10 : !> \par History
11 : !> 03.2022 created [F. Stein]
12 : !> \author Frederick Stein
13 : ! **************************************************************************************************
14 : MODULE xc_util
15 : USE pw_methods, ONLY: pw_axpy, &
16 : pw_copy, &
17 : pw_derive, &
18 : pw_laplace, &
19 : pw_transfer, &
20 : pw_zero
21 : USE pw_pool_types, ONLY: pw_pool_type
22 : USE pw_spline_utils, ONLY: &
23 : nn10_coeffs, nn10_deriv_coeffs, nn50_coeffs, nn50_deriv_coeffs, pw_nn_deriv_r, &
24 : pw_nn_smear_r, pw_spline2_deriv_g, pw_spline2_interpolate_values_g, pw_spline3_deriv_g, &
25 : pw_spline3_interpolate_values_g, pw_spline_scale_deriv, spline2_coeffs, &
26 : spline2_deriv_coeffs, spline3_coeffs, spline3_deriv_coeffs
27 : USE pw_types, ONLY: &
28 : pw_c1d_gs_type, pw_r3d_rs_type
29 : USE xc_input_constants, ONLY: &
30 : xc_deriv_nn10_smooth, xc_deriv_nn50_smooth, xc_deriv_pw, xc_deriv_spline2, &
31 : xc_deriv_spline2_smooth, xc_deriv_spline3, xc_deriv_spline3_smooth, xc_rho_nn10, &
32 : xc_rho_nn50, xc_rho_no_smooth, xc_rho_spline2_smooth, xc_rho_spline3_smooth
33 : #include "../base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 :
37 : PRIVATE
38 :
39 : PUBLIC :: xc_pw_smooth, xc_pw_laplace, xc_pw_divergence, xc_pw_derive, xc_requires_tmp_g, xc_pw_gradient
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_util'
41 :
42 : INTERFACE xc_pw_derive
43 : MODULE PROCEDURE xc_pw_derive_r3d_rs, xc_pw_derive_c1d_gs
44 : END INTERFACE
45 :
46 : INTERFACE xc_pw_laplace
47 : MODULE PROCEDURE xc_pw_laplace_r3d_rs, xc_pw_laplace_c1d_gs
48 : END INTERFACE
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief ...
54 : !> \param xc_deriv_id ...
55 : !> \return ...
56 : ! **************************************************************************************************
57 646057 : ELEMENTAL FUNCTION xc_requires_tmp_g(xc_deriv_id) RESULT(requires)
58 : INTEGER, INTENT(IN) :: xc_deriv_id
59 : LOGICAL :: requires
60 :
61 : requires = (xc_deriv_id == xc_deriv_spline2) .OR. &
62 : (xc_deriv_id == xc_deriv_spline3) .OR. &
63 646057 : (xc_deriv_id == xc_deriv_pw)
64 646057 : END FUNCTION
65 :
66 : ! **************************************************************************************************
67 : !> \brief ...
68 : !> \param pw_in ...
69 : !> \param pw_out ...
70 : !> \param xc_smooth_id ...
71 : ! **************************************************************************************************
72 132147 : SUBROUTINE xc_pw_smooth(pw_in, pw_out, xc_smooth_id)
73 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
74 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
75 : INTEGER, INTENT(IN) :: xc_smooth_id
76 :
77 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_pw_smooth'
78 :
79 : INTEGER :: handle
80 :
81 132147 : CALL timeset(routineN, handle)
82 :
83 263324 : SELECT CASE (xc_smooth_id)
84 : CASE (xc_rho_no_smooth)
85 131177 : CALL pw_copy(pw_in, pw_out)
86 : CASE (xc_rho_spline2_smooth)
87 0 : CALL pw_zero(pw_out)
88 : CALL pw_nn_smear_r(pw_in=pw_in, &
89 : pw_out=pw_out, &
90 0 : coeffs=spline2_coeffs)
91 : CASE (xc_rho_spline3_smooth)
92 0 : CALL pw_zero(pw_out)
93 : CALL pw_nn_smear_r(pw_in=pw_in, &
94 : pw_out=pw_out, &
95 0 : coeffs=spline3_coeffs)
96 : CASE (xc_rho_nn10)
97 752 : CALL pw_zero(pw_out)
98 : CALL pw_nn_smear_r(pw_in=pw_in, &
99 : pw_out=pw_out, &
100 752 : coeffs=nn10_coeffs)
101 : CASE (xc_rho_nn50)
102 218 : CALL pw_zero(pw_out)
103 : CALL pw_nn_smear_r(pw_in=pw_in, &
104 : pw_out=pw_out, &
105 218 : coeffs=nn50_coeffs)
106 : CASE default
107 132147 : CPABORT("Unsupported smoothing")
108 : END SELECT
109 :
110 132147 : CALL timestop(handle)
111 :
112 132147 : END SUBROUTINE xc_pw_smooth
113 :
114 : ! **************************************************************************************************
115 : !> \brief ...
116 : !> \param pw_r ...
117 : !> \param pw_g ...
118 : !> \param tmp_g ...
119 : !> \param gradient ...
120 : !> \param xc_deriv_method_id ...
121 : ! **************************************************************************************************
122 102575 : SUBROUTINE xc_pw_gradient(pw_r, pw_g, tmp_g, gradient, xc_deriv_method_id)
123 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_r
124 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw_g, tmp_g
125 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: gradient
126 : INTEGER, INTENT(IN) :: xc_deriv_method_id
127 :
128 : INTEGER :: idir
129 :
130 410300 : DO idir = 1, 3
131 307725 : CALL pw_zero(gradient(idir))
132 410300 : CALL xc_pw_derive(pw_r, tmp_g, gradient(idir), idir, xc_deriv_method_id, pw_g=pw_g)
133 : END DO
134 :
135 102575 : END SUBROUTINE xc_pw_gradient
136 :
137 : #:for kind in ["r3d_rs", "c1d_gs"]
138 : ! **************************************************************************************************
139 : !> \brief Calculates the Laplacian of pw
140 : !> \param pw on input: pw of which the Laplacian shall be calculated, on output if pw_out is absent: Laplacian of input
141 : !> \param pw_pool ...
142 : !> \param deriv_method_id ...
143 : !> \param pw_out if present, save the Laplacian of pw here
144 : !> \param tmp_g scratch grid in reciprocal space, used instead of the internal grid if given explicitly to save memory
145 : ! **************************************************************************************************
146 2406 : SUBROUTINE xc_pw_laplace_${kind}$ (pw, pw_pool, deriv_method_id, pw_out, tmp_g)
147 : TYPE(pw_${kind}$_type), INTENT(INOUT) :: pw
148 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
149 : INTEGER, INTENT(IN) :: deriv_method_id
150 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: pw_out
151 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: tmp_g
152 :
153 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_pw_laplace'
154 :
155 : INTEGER :: handle
156 : LOGICAL :: owns_tmp_g
157 : TYPE(pw_c1d_gs_type) :: my_tmp_g
158 :
159 2406 : CALL timeset(routineN, handle)
160 :
161 4812 : SELECT CASE (deriv_method_id)
162 : CASE (xc_deriv_pw)
163 :
164 2406 : IF (PRESENT(tmp_g)) my_tmp_g = tmp_g
165 :
166 2406 : owns_tmp_g = .FALSE.
167 2406 : IF (.NOT. ASSOCIATED(my_tmp_g%pw_grid)) THEN
168 1074 : CALL pw_pool%create_pw(my_tmp_g)
169 1074 : owns_tmp_g = .TRUE.
170 : END IF
171 2406 : CALL pw_zero(my_tmp_g)
172 2406 : CALL pw_transfer(pw, my_tmp_g)
173 :
174 2406 : CALL pw_laplace(my_tmp_g)
175 :
176 2406 : IF (PRESENT(pw_out)) THEN
177 1332 : CALL pw_transfer(my_tmp_g, pw_out)
178 : ELSE
179 1074 : CALL pw_transfer(my_tmp_g, pw)
180 : END IF
181 2406 : IF (owns_tmp_g) THEN
182 1074 : CALL pw_pool%give_back_pw(my_tmp_g)
183 : END IF
184 : CASE default
185 2406 : CPABORT("Unsupported derivative method")
186 : END SELECT
187 :
188 2406 : CALL timestop(handle)
189 :
190 2406 : END SUBROUTINE xc_pw_laplace_${kind}$
191 : #:endfor
192 :
193 : ! **************************************************************************************************
194 : !> \brief Calculates the divergence of pw_to_deriv
195 : !> \param xc_deriv_method_id ...
196 : !> \param pw_to_deriv ...
197 : !> \param tmp_g ...
198 : !> \param vxc_g ...
199 : !> \param vxc_r ...
200 : ! **************************************************************************************************
201 86541 : SUBROUTINE xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_r)
202 : INTEGER, INTENT(IN) :: xc_deriv_method_id
203 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: pw_to_deriv
204 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: tmp_g, vxc_g
205 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
206 :
207 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_pw_divergence'
208 :
209 : INTEGER :: handle, idir
210 :
211 86541 : CALL timeset(routineN, handle)
212 :
213 : ! partial integration
214 86541 : IF (xc_deriv_method_id /= xc_deriv_pw) THEN
215 8212 : CALL pw_spline_scale_deriv(pw_to_deriv, transpose=.TRUE.)
216 : END IF
217 :
218 86541 : IF (ASSOCIATED(vxc_g%pw_grid)) CALL pw_zero(vxc_g)
219 :
220 346164 : DO idir = 1, 3
221 259623 : CALL xc_pw_derive(pw_to_deriv(idir), tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr=.FALSE.)
222 346164 : IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(vxc_g%pw_grid)) CALL pw_axpy(tmp_g, vxc_g)
223 : END DO
224 :
225 86541 : IF (ASSOCIATED(vxc_g%pw_grid)) THEN
226 84111 : CALL pw_transfer(vxc_g, pw_to_deriv(1))
227 84111 : CALL pw_axpy(pw_to_deriv(1), vxc_r)
228 : END IF
229 :
230 86541 : CALL timestop(handle)
231 :
232 86541 : END SUBROUTINE xc_pw_divergence
233 :
234 : #:for kind in ["r3d_rs", "c1d_gs"]
235 : ! **************************************************************************************************
236 : !> \brief Calculates the derivative of a function on a planewave grid in a given direction
237 : !> \param pw function to derive
238 : !> \param tmp_g temporary grid in reciprocal space, only required if derivative method is pw or spline
239 : !> \param vxc_r if tmp_g is not required, add derivative here
240 : !> \param idir direction of derivative
241 : !> \param xc_deriv_method_id ...
242 : !> \param copy_to_vxcr ...
243 : !> \param pw_g ...
244 : ! **************************************************************************************************
245 567348 : SUBROUTINE xc_pw_derive_${kind}$ (pw, tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr, pw_g)
246 : TYPE(pw_${kind}$_type), INTENT(IN) :: pw
247 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: tmp_g
248 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
249 : INTEGER, INTENT(IN) :: idir, xc_deriv_method_id
250 : LOGICAL, INTENT(IN), OPTIONAL :: copy_to_vxcr
251 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: pw_g
252 :
253 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_pw_derive'
254 : INTEGER, DIMENSION(3, 3), PARAMETER :: nd = RESHAPE((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
255 :
256 : INTEGER :: handle
257 : LOGICAL :: my_copy_to_vxcr
258 : #:if kind=="c1d_gs"
259 : TYPE(pw_r3d_rs_type) :: tmp_r
260 : #:endif
261 :
262 567348 : CALL timeset(routineN, handle)
263 :
264 567348 : my_copy_to_vxcr = .TRUE.
265 567348 : IF (PRESENT(copy_to_vxcr)) my_copy_to_vxcr = copy_to_vxcr
266 :
267 567348 : IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
268 :
269 551946 : IF (PRESENT(pw_g)) THEN
270 299613 : IF (ASSOCIATED(pw_g%pw_grid)) THEN
271 299613 : CALL pw_copy(pw_g, tmp_g)
272 : ELSE
273 0 : CALL pw_transfer(pw, tmp_g)
274 : END IF
275 : ELSE
276 252333 : CALL pw_transfer(pw, tmp_g)
277 : END IF
278 :
279 1064964 : SELECT CASE (xc_deriv_method_id)
280 : CASE (xc_deriv_pw)
281 513018 : CALL pw_derive(tmp_g, nd(:, idir))
282 : CASE (xc_deriv_spline2)
283 38532 : CALL pw_spline2_interpolate_values_g(tmp_g)
284 38532 : CALL pw_spline2_deriv_g(tmp_g, idir=idir)
285 : CASE (xc_deriv_spline3)
286 396 : CALL pw_spline3_interpolate_values_g(tmp_g)
287 396 : CALL pw_spline3_deriv_g(tmp_g, idir=idir)
288 : CASE default
289 551946 : CPABORT("Unsupported deriv method")
290 : END SELECT
291 :
292 551946 : IF (my_copy_to_vxcr) CALL pw_transfer(tmp_g, vxc_r)
293 : ELSE
294 : #:if kind=="r3d_rs"
295 26802 : SELECT CASE (xc_deriv_method_id)
296 : CASE (xc_deriv_spline2_smooth)
297 : CALL pw_nn_deriv_r(pw_in=pw, &
298 : pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
299 11400 : idir=idir)
300 : CASE (xc_deriv_spline3_smooth)
301 : CALL pw_nn_deriv_r(pw_in=pw, &
302 : pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
303 1014 : idir=idir)
304 : CASE (xc_deriv_nn10_smooth)
305 : CALL pw_nn_deriv_r(pw_in=pw, &
306 : pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
307 1230 : idir=idir)
308 : CASE (xc_deriv_nn50_smooth)
309 : CALL pw_nn_deriv_r(pw_in=pw, &
310 : pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
311 1758 : idir=idir)
312 : CASE default
313 15402 : CPABORT("Unsupported derivative method")
314 : END SELECT
315 : #:else
316 0 : CALL tmp_r%create(pw%pw_grid)
317 0 : SELECT CASE (xc_deriv_method_id)
318 : CASE (xc_deriv_spline2_smooth)
319 : CALL pw_nn_deriv_r(pw_in=tmp_r, &
320 : pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
321 0 : idir=idir)
322 : CASE (xc_deriv_spline3_smooth)
323 : CALL pw_nn_deriv_r(pw_in=tmp_r, &
324 : pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
325 0 : idir=idir)
326 : CASE (xc_deriv_nn10_smooth)
327 : CALL pw_nn_deriv_r(pw_in=tmp_r, &
328 : pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
329 0 : idir=idir)
330 : CASE (xc_deriv_nn50_smooth)
331 : CALL pw_nn_deriv_r(pw_in=tmp_r, &
332 : pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
333 0 : idir=idir)
334 : CASE default
335 0 : CPABORT("Unsupported derivative method")
336 : END SELECT
337 0 : CALL tmp_r%release()
338 : #:endif
339 : END IF
340 :
341 567348 : CALL timestop(handle)
342 :
343 567348 : END SUBROUTINE xc_pw_derive_${kind}$
344 : #:endfor
345 :
346 : END MODULE xc_util
|