(git:b195825)
xc_util.F
Go to the documentation of this file.
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: &
27  USE pw_types, ONLY: &
28  pw_c1d_gs_type, pw_r3d_rs_type
29  USE xc_input_constants, ONLY: &
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  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  (xc_deriv_id == xc_deriv_pw)
64  END FUNCTION
65 
66 ! **************************************************************************************************
67 !> \brief ...
68 !> \param pw_in ...
69 !> \param pw_out ...
70 !> \param xc_smooth_id ...
71 ! **************************************************************************************************
72  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  CALL timeset(routinen, handle)
82 
83  SELECT CASE (xc_smooth_id)
84  CASE (xc_rho_no_smooth)
85  CALL pw_copy(pw_in, pw_out)
87  CALL pw_zero(pw_out)
88  CALL pw_nn_smear_r(pw_in=pw_in, &
89  pw_out=pw_out, &
90  coeffs=spline2_coeffs)
92  CALL pw_zero(pw_out)
93  CALL pw_nn_smear_r(pw_in=pw_in, &
94  pw_out=pw_out, &
95  coeffs=spline3_coeffs)
96  CASE (xc_rho_nn10)
97  CALL pw_zero(pw_out)
98  CALL pw_nn_smear_r(pw_in=pw_in, &
99  pw_out=pw_out, &
100  coeffs=nn10_coeffs)
101  CASE (xc_rho_nn50)
102  CALL pw_zero(pw_out)
103  CALL pw_nn_smear_r(pw_in=pw_in, &
104  pw_out=pw_out, &
105  coeffs=nn50_coeffs)
106  CASE default
107  cpabort("Unsupported smoothing")
108  END SELECT
109 
110  CALL timestop(handle)
111 
112  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  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  DO idir = 1, 3
131  CALL pw_zero(gradient(idir))
132  CALL xc_pw_derive(pw_r, tmp_g, gradient(idir), idir, xc_deriv_method_id, pw_g=pw_g)
133  END DO
134 
135  END SUBROUTINE xc_pw_gradient
136 
137 ! **************************************************************************************************
138 !> \brief Calculates the Laplacian of pw
139 !> \param pw on input: pw of which the Laplacian shall be calculated, on output if pw_out is absent: Laplacian of input
140 !> \param pw_pool ...
141 !> \param deriv_method_id ...
142 !> \param pw_out if present, save the Laplacian of pw here
143 !> \param tmp_g scratch grid in reciprocal space, used instead of the internal grid if given explicitly to save memory
144 ! **************************************************************************************************
145  SUBROUTINE xc_pw_laplace_r3d_rs (pw, pw_pool, deriv_method_id, pw_out, tmp_g)
146  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw
147  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
148  INTEGER, INTENT(IN) :: deriv_method_id
149  TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: pw_out
150  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: tmp_g
151 
152  CHARACTER(len=*), PARAMETER :: routinen = 'xc_pw_laplace'
153 
154  INTEGER :: handle
155  LOGICAL :: owns_tmp_g
156  TYPE(pw_c1d_gs_type) :: my_tmp_g
157 
158  CALL timeset(routinen, handle)
159 
160  SELECT CASE (deriv_method_id)
161  CASE (xc_deriv_pw)
162 
163  IF (PRESENT(tmp_g)) my_tmp_g = tmp_g
164 
165  owns_tmp_g = .false.
166  IF (.NOT. ASSOCIATED(my_tmp_g%pw_grid)) THEN
167  CALL pw_pool%create_pw(my_tmp_g)
168  owns_tmp_g = .true.
169  END IF
170  CALL pw_zero(my_tmp_g)
171  CALL pw_transfer(pw, my_tmp_g)
172 
173  CALL pw_laplace(my_tmp_g)
174 
175  IF (PRESENT(pw_out)) THEN
176  CALL pw_transfer(my_tmp_g, pw_out)
177  ELSE
178  CALL pw_transfer(my_tmp_g, pw)
179  END IF
180  IF (owns_tmp_g) THEN
181  CALL pw_pool%give_back_pw(my_tmp_g)
182  END IF
183  CASE default
184  cpabort("Unsupported derivative method")
185  END SELECT
186 
187  CALL timestop(handle)
188 
189  END SUBROUTINE xc_pw_laplace_r3d_rs
190 ! **************************************************************************************************
191 !> \brief Calculates the Laplacian of pw
192 !> \param pw on input: pw of which the Laplacian shall be calculated, on output if pw_out is absent: Laplacian of input
193 !> \param pw_pool ...
194 !> \param deriv_method_id ...
195 !> \param pw_out if present, save the Laplacian of pw here
196 !> \param tmp_g scratch grid in reciprocal space, used instead of the internal grid if given explicitly to save memory
197 ! **************************************************************************************************
198  SUBROUTINE xc_pw_laplace_c1d_gs (pw, pw_pool, deriv_method_id, pw_out, tmp_g)
199  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: pw
200  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
201  INTEGER, INTENT(IN) :: deriv_method_id
202  TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: pw_out
203  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: tmp_g
204 
205  CHARACTER(len=*), PARAMETER :: routinen = 'xc_pw_laplace'
206 
207  INTEGER :: handle
208  LOGICAL :: owns_tmp_g
209  TYPE(pw_c1d_gs_type) :: my_tmp_g
210 
211  CALL timeset(routinen, handle)
212 
213  SELECT CASE (deriv_method_id)
214  CASE (xc_deriv_pw)
215 
216  IF (PRESENT(tmp_g)) my_tmp_g = tmp_g
217 
218  owns_tmp_g = .false.
219  IF (.NOT. ASSOCIATED(my_tmp_g%pw_grid)) THEN
220  CALL pw_pool%create_pw(my_tmp_g)
221  owns_tmp_g = .true.
222  END IF
223  CALL pw_zero(my_tmp_g)
224  CALL pw_transfer(pw, my_tmp_g)
225 
226  CALL pw_laplace(my_tmp_g)
227 
228  IF (PRESENT(pw_out)) THEN
229  CALL pw_transfer(my_tmp_g, pw_out)
230  ELSE
231  CALL pw_transfer(my_tmp_g, pw)
232  END IF
233  IF (owns_tmp_g) THEN
234  CALL pw_pool%give_back_pw(my_tmp_g)
235  END IF
236  CASE default
237  cpabort("Unsupported derivative method")
238  END SELECT
239 
240  CALL timestop(handle)
241 
242  END SUBROUTINE xc_pw_laplace_c1d_gs
243 
244 ! **************************************************************************************************
245 !> \brief Calculates the divergence of pw_to_deriv
246 !> \param xc_deriv_method_id ...
247 !> \param pw_to_deriv ...
248 !> \param tmp_g ...
249 !> \param vxc_g ...
250 !> \param vxc_r ...
251 ! **************************************************************************************************
252  SUBROUTINE xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_r)
253  INTEGER, INTENT(IN) :: xc_deriv_method_id
254  TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: pw_to_deriv
255  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: tmp_g, vxc_g
256  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
257 
258  CHARACTER(len=*), PARAMETER :: routinen = 'xc_pw_divergence'
259 
260  INTEGER :: handle, idir
261 
262  CALL timeset(routinen, handle)
263 
264  ! partial integration
265  IF (xc_deriv_method_id /= xc_deriv_pw) THEN
266  CALL pw_spline_scale_deriv(pw_to_deriv, transpose=.true.)
267  END IF
268 
269  IF (ASSOCIATED(vxc_g%pw_grid)) CALL pw_zero(vxc_g)
270 
271  DO idir = 1, 3
272  CALL xc_pw_derive(pw_to_deriv(idir), tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr=.false.)
273  IF (ASSOCIATED(tmp_g%pw_grid) .AND. ASSOCIATED(vxc_g%pw_grid)) CALL pw_axpy(tmp_g, vxc_g)
274  END DO
275 
276  IF (ASSOCIATED(vxc_g%pw_grid)) THEN
277  CALL pw_transfer(vxc_g, pw_to_deriv(1))
278  CALL pw_axpy(pw_to_deriv(1), vxc_r)
279  END IF
280 
281  CALL timestop(handle)
282 
283  END SUBROUTINE xc_pw_divergence
284 
285 ! **************************************************************************************************
286 !> \brief Calculates the derivative of a function on a planewave grid in a given direction
287 !> \param pw function to derive
288 !> \param tmp_g temporary grid in reciprocal space, only required if derivative method is pw or spline
289 !> \param vxc_r if tmp_g is not required, add derivative here
290 !> \param idir direction of derivative
291 !> \param xc_deriv_method_id ...
292 !> \param copy_to_vxcr ...
293 !> \param pw_g ...
294 ! **************************************************************************************************
295  SUBROUTINE xc_pw_derive_r3d_rs (pw, tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr, pw_g)
296  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
297  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: tmp_g
298  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
299  INTEGER, INTENT(IN) :: idir, xc_deriv_method_id
300  LOGICAL, INTENT(IN), OPTIONAL :: copy_to_vxcr
301  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: pw_g
302 
303  CHARACTER(len=*), PARAMETER :: routinen = 'xc_pw_derive'
304  INTEGER, DIMENSION(3, 3), PARAMETER :: nd = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
305 
306  INTEGER :: handle
307  LOGICAL :: my_copy_to_vxcr
308 
309  CALL timeset(routinen, handle)
310 
311  my_copy_to_vxcr = .true.
312  IF (PRESENT(copy_to_vxcr)) my_copy_to_vxcr = copy_to_vxcr
313 
314  IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
315 
316  IF (PRESENT(pw_g)) THEN
317  IF (ASSOCIATED(pw_g%pw_grid)) THEN
318  CALL pw_copy(pw_g, tmp_g)
319  ELSE
320  CALL pw_transfer(pw, tmp_g)
321  END IF
322  ELSE
323  CALL pw_transfer(pw, tmp_g)
324  END IF
325 
326  SELECT CASE (xc_deriv_method_id)
327  CASE (xc_deriv_pw)
328  CALL pw_derive(tmp_g, nd(:, idir))
329  CASE (xc_deriv_spline2)
331  CALL pw_spline2_deriv_g(tmp_g, idir=idir)
332  CASE (xc_deriv_spline3)
334  CALL pw_spline3_deriv_g(tmp_g, idir=idir)
335  CASE default
336  cpabort("Unsupported deriv method")
337  END SELECT
338 
339  IF (my_copy_to_vxcr) CALL pw_transfer(tmp_g, vxc_r)
340  ELSE
341  SELECT CASE (xc_deriv_method_id)
343  CALL pw_nn_deriv_r(pw_in=pw, &
344  pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
345  idir=idir)
347  CALL pw_nn_deriv_r(pw_in=pw, &
348  pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
349  idir=idir)
350  CASE (xc_deriv_nn10_smooth)
351  CALL pw_nn_deriv_r(pw_in=pw, &
352  pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
353  idir=idir)
354  CASE (xc_deriv_nn50_smooth)
355  CALL pw_nn_deriv_r(pw_in=pw, &
356  pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
357  idir=idir)
358  CASE default
359  cpabort("Unsupported derivative method")
360  END SELECT
361  END IF
362 
363  CALL timestop(handle)
364 
365  END SUBROUTINE xc_pw_derive_r3d_rs
366 ! **************************************************************************************************
367 !> \brief Calculates the derivative of a function on a planewave grid in a given direction
368 !> \param pw function to derive
369 !> \param tmp_g temporary grid in reciprocal space, only required if derivative method is pw or spline
370 !> \param vxc_r if tmp_g is not required, add derivative here
371 !> \param idir direction of derivative
372 !> \param xc_deriv_method_id ...
373 !> \param copy_to_vxcr ...
374 !> \param pw_g ...
375 ! **************************************************************************************************
376  SUBROUTINE xc_pw_derive_c1d_gs (pw, tmp_g, vxc_r, idir, xc_deriv_method_id, copy_to_vxcr, pw_g)
377  TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
378  TYPE(pw_c1d_gs_type), INTENT(INOUT) :: tmp_g
379  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vxc_r
380  INTEGER, INTENT(IN) :: idir, xc_deriv_method_id
381  LOGICAL, INTENT(IN), OPTIONAL :: copy_to_vxcr
382  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: pw_g
383 
384  CHARACTER(len=*), PARAMETER :: routinen = 'xc_pw_derive'
385  INTEGER, DIMENSION(3, 3), PARAMETER :: nd = reshape((/1, 0, 0, 0, 1, 0, 0, 0, 1/), (/3, 3/))
386 
387  INTEGER :: handle
388  LOGICAL :: my_copy_to_vxcr
389  TYPE(pw_r3d_rs_type) :: tmp_r
390 
391  CALL timeset(routinen, handle)
392 
393  my_copy_to_vxcr = .true.
394  IF (PRESENT(copy_to_vxcr)) my_copy_to_vxcr = copy_to_vxcr
395 
396  IF (xc_requires_tmp_g(xc_deriv_method_id)) THEN
397 
398  IF (PRESENT(pw_g)) THEN
399  IF (ASSOCIATED(pw_g%pw_grid)) THEN
400  CALL pw_copy(pw_g, tmp_g)
401  ELSE
402  CALL pw_transfer(pw, tmp_g)
403  END IF
404  ELSE
405  CALL pw_transfer(pw, tmp_g)
406  END IF
407 
408  SELECT CASE (xc_deriv_method_id)
409  CASE (xc_deriv_pw)
410  CALL pw_derive(tmp_g, nd(:, idir))
411  CASE (xc_deriv_spline2)
413  CALL pw_spline2_deriv_g(tmp_g, idir=idir)
414  CASE (xc_deriv_spline3)
416  CALL pw_spline3_deriv_g(tmp_g, idir=idir)
417  CASE default
418  cpabort("Unsupported deriv method")
419  END SELECT
420 
421  IF (my_copy_to_vxcr) CALL pw_transfer(tmp_g, vxc_r)
422  ELSE
423  CALL tmp_r%create(pw%pw_grid)
424  SELECT CASE (xc_deriv_method_id)
426  CALL pw_nn_deriv_r(pw_in=tmp_r, &
427  pw_out=vxc_r, coeffs=spline2_deriv_coeffs, &
428  idir=idir)
430  CALL pw_nn_deriv_r(pw_in=tmp_r, &
431  pw_out=vxc_r, coeffs=spline3_deriv_coeffs, &
432  idir=idir)
433  CASE (xc_deriv_nn10_smooth)
434  CALL pw_nn_deriv_r(pw_in=tmp_r, &
435  pw_out=vxc_r, coeffs=nn10_deriv_coeffs, &
436  idir=idir)
437  CASE (xc_deriv_nn50_smooth)
438  CALL pw_nn_deriv_r(pw_in=tmp_r, &
439  pw_out=vxc_r, coeffs=nn50_deriv_coeffs, &
440  idir=idir)
441  CASE default
442  cpabort("Unsupported derivative method")
443  END SELECT
444  CALL tmp_r%release()
445  END IF
446 
447  CALL timestop(handle)
448 
449  END SUBROUTINE xc_pw_derive_c1d_gs
450 
451 END MODULE xc_util
subroutine, public pw_laplace(pw)
Calculate the Laplacian of a plane wave vector.
Definition: pw_methods.F:10175
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
different utils that are useful to manipulate splines on the regular grid of a pw
subroutine, public pw_nn_deriv_r(pw_in, pw_out, coeffs, idir)
calculates a nearest neighbor central derivative. for the x dir: pw_outarray(i,j,k)=( pw_in(i+1,...
subroutine, public pw_spline3_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the cubic spline
real(kind=dp), dimension(4), parameter, public spline2_coeffs
real(kind=dp), dimension(3), parameter, public spline3_deriv_coeffs
subroutine, public pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
rescales the derivatives from gridspacing=1 to the real derivatives
real(kind=dp), dimension(3), parameter, public spline2_deriv_coeffs
subroutine, public pw_nn_smear_r(pw_in, pw_out, coeffs)
calculates the values of a nearest neighbor smearing
real(kind=dp), dimension(4), parameter, public spline3_coeffs
real(kind=dp), dimension(3), parameter, public nn50_deriv_coeffs
subroutine, public pw_spline3_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the2 cubic spline that interpolates the given values
subroutine, public pw_spline2_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the quadratic spline that interpolates the given values
real(kind=dp), dimension(3), parameter, public nn10_deriv_coeffs
subroutine, public pw_spline2_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the quadratic spline
real(kind=dp), dimension(4), parameter, public nn10_coeffs
real(kind=dp), dimension(4), parameter, public nn50_coeffs
input constants for xc
integer, parameter, public xc_deriv_spline2_smooth
integer, parameter, public xc_rho_spline2_smooth
integer, parameter, public xc_rho_spline3_smooth
integer, parameter, public xc_deriv_nn10_smooth
integer, parameter, public xc_deriv_pw
integer, parameter, public xc_deriv_spline2
integer, parameter, public xc_rho_nn10
integer, parameter, public xc_deriv_nn50_smooth
integer, parameter, public xc_deriv_spline3
integer, parameter, public xc_rho_nn50
integer, parameter, public xc_deriv_spline3_smooth
integer, parameter, public xc_rho_no_smooth
contains utility functions for the xc package
Definition: xc_util.F:14
subroutine, public xc_pw_divergence(xc_deriv_method_id, pw_to_deriv, tmp_g, vxc_g, vxc_r)
Calculates the divergence of pw_to_deriv.
Definition: xc_util.F:253
subroutine, public xc_pw_smooth(pw_in, pw_out, xc_smooth_id)
...
Definition: xc_util.F:73
elemental logical function, public xc_requires_tmp_g(xc_deriv_id)
...
Definition: xc_util.F:58
subroutine, public xc_pw_gradient(pw_r, pw_g, tmp_g, gradient, xc_deriv_method_id)
...
Definition: xc_util.F:123