(git:34ef472)
xc_xlda_hole_t_c_lr.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 Calculates the lda exchange hole in a truncated coulomb potential.
10 !> Can be used as longrange correction for truncated hfx calculations
11 !> \par History
12 !> Manuel Guidon (12.2008) : created
13 !> \author Manuel Guidon (06.2008)
14 ! **************************************************************************************************
15 
16 MODULE xc_xlda_hole_t_c_lr
17 
18  USE input_section_types, ONLY: section_vals_type,&
19  section_vals_val_get
20  USE kinds, ONLY: dp
21  USE mathconstants, ONLY: pi
22  USE mathlib, ONLY: expint
23  USE xc_derivative_desc, ONLY: deriv_rho,&
24  deriv_rhoa,&
25  deriv_rhob
26  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
27  xc_dset_get_derivative
28  USE xc_derivative_types, ONLY: xc_derivative_get,&
29  xc_derivative_type
30  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
31  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
32  xc_rho_set_type
33 #include "../base/base_uses.f90"
34 
35  IMPLICIT NONE
36 
37  PRIVATE
38 
39 ! *** Global parameters ***
40 
41  PUBLIC :: xlda_hole_t_c_lr_lda_eval, xlda_hole_t_c_lr_lda_info, &
42  xlda_hole_t_c_lr_lsd_eval, xlda_hole_t_c_lr_lsd_info, &
43  xlda_hole_t_c_lr_lda_calc_0
44 
45  REAL(KIND=dp), PARAMETER :: a = 1.0161144_dp, &
46  b = -0.37170836_dp, &
47  c = -0.077215461_dp, &
48  d = 0.57786348_dp, &
49  e = -0.051955731_dp
50 
51  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xlda_hole_t_c_lr'
52 
53 CONTAINS
54 
55 ! **************************************************************************************************
56 !> \brief returns various information on the functional
57 !> \param reference string with the reference of the actual functional
58 !> \param shortform string with the shortform of the functional name
59 !> \param needs the components needed by this functional are set to
60 !> true (does not set the unneeded components to false)
61 !> \param max_deriv controls the number of derivatives
62 !> \par History
63 !> 12.2008 created [mguidon]
64 !> \author mguidon
65 ! **************************************************************************************************
66  SUBROUTINE xlda_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
67  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
68  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
69  INTEGER, INTENT(out), OPTIONAL :: max_deriv
70 
71  IF (PRESENT(reference)) THEN
72  reference = "{LDA version}"
73  END IF
74  IF (PRESENT(shortform)) THEN
75  shortform = "{LDA}"
76  END IF
77  IF (PRESENT(needs)) THEN
78  needs%rho = .true.
79  END IF
80  IF (PRESENT(max_deriv)) max_deriv = 1
81 
82  END SUBROUTINE xlda_hole_t_c_lr_lda_info
83 
84 ! **************************************************************************************************
85 !> \brief returns various information on the functional
86 !> \param reference string with the reference of the actual functional
87 !> \param shortform string with the shortform of the functional name
88 !> \param needs the components needed by this functional are set to
89 !> true (does not set the unneeded components to false)
90 !> \param max_deriv controls the number of derivatives
91 !> \par History
92 !> 12.2008 created [mguidon]
93 !> \author mguidon
94 ! **************************************************************************************************
95  SUBROUTINE xlda_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
96  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
97  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
98  INTEGER, INTENT(out), OPTIONAL :: max_deriv
99 
100  IF (PRESENT(reference)) THEN
101  reference = "{LSD version}"
102  END IF
103  IF (PRESENT(shortform)) THEN
104  shortform = "{LSD}"
105  END IF
106  IF (PRESENT(needs)) THEN
107  needs%rho_spin = .true.
108  END IF
109  IF (PRESENT(max_deriv)) max_deriv = 1
110 
111  END SUBROUTINE xlda_hole_t_c_lr_lsd_info
112 
113 ! **************************************************************************************************
114 !> \brief evaluates the truncated lda exchange hole
115 !> \param rho_set the density where you want to evaluate the functional
116 !> \param deriv_set place where to store the functional derivatives (they are
117 !> added to the derivatives)
118 !> \param order degree of the derivative that should be evaluated,
119 !> if positive all the derivatives up to the given degree are evaluated,
120 !> if negative only the given degree is calculated
121 !> \param params input parameters (scaling, cutoff_radius)
122 !> \par History
123 !> 12.2008 created [Manuel Guidon]
124 !> \author Manuel Guidon
125 ! **************************************************************************************************
126  SUBROUTINE xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
127 
128  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
129  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
130  INTEGER, INTENT(IN) :: order
131  TYPE(section_vals_type), POINTER :: params
132 
133  CHARACTER(len=*), PARAMETER :: routinen = 'xlda_hole_t_c_lr_lda_eval'
134 
135  INTEGER :: handle, npoints
136  INTEGER, DIMENSION(2, 3) :: bo
137  REAL(kind=dp) :: epsilon_rho, r, sx
138  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
139  POINTER :: dummy, e_0, e_rho, rho
140  TYPE(xc_derivative_type), POINTER :: deriv
141 
142  CALL timeset(routinen, handle)
143 
144  CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
145  CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
146 
147  CALL xc_rho_set_get(rho_set, rho=rho, &
148  local_bounds=bo, rho_cutoff=epsilon_rho)
149  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
150 
151  dummy => rho
152 
153  e_0 => dummy
154  e_rho => dummy
155 
156  IF (order >= 0) THEN
157  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
158  allocate_deriv=.true.)
159  CALL xc_derivative_get(deriv, deriv_data=e_0)
160  END IF
161  IF (order >= 1 .OR. order == -1) THEN
162  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
163  allocate_deriv=.true.)
164  CALL xc_derivative_get(deriv, deriv_data=e_rho)
165  END IF
166  IF (order > 1 .OR. order < -1) THEN
167  cpabort("derivatives bigger than 1 not implemented")
168  END IF
169 
170  IF (r == 0.0_dp) THEN
171  cpabort("Cutoff_Radius 0.0 not implemented")
172  END IF
173  CALL xlda_hole_t_c_lr_lda_calc(npoints, order, rho=rho, &
174  e_0=e_0, e_rho=e_rho, &
175  epsilon_rho=epsilon_rho, &
176  sx=sx, r=r)
177 
178  CALL timestop(handle)
179 
180  END SUBROUTINE xlda_hole_t_c_lr_lda_eval
181 
182 ! **************************************************************************************************
183 !> \brief Call low level routine
184 !> \param npoints ...
185 !> \param order ...
186 !> \param rho ...
187 !> \param e_0 ...
188 !> \param e_rho ...
189 !> \param epsilon_rho ...
190 !> \param sx ...
191 !> \param R ...
192 !> \par History
193 !> 12.2008 created [Manuel Guidon]
194 !> \author Manuel Guidon
195 ! **************************************************************************************************
196  SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho, &
197  epsilon_rho, sx, R)
198 
199  INTEGER, INTENT(in) :: npoints, order
200  REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
201  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
202 
203  INTEGER :: ip
204  REAL(dp) :: my_rho
205 
206 !$OMP PARALLEL DO DEFAULT(NONE) &
207 !$OMP SHARED(npoints, rho, epsilon_rho, order, e_0, e_rho) &
208 !$OMP SHARED(sx, r) &
209 !$OMP PRIVATE(ip, my_rho)
210 
211  DO ip = 1, npoints
212  my_rho = max(rho(ip), 0.0_dp)
213  IF (my_rho > epsilon_rho) THEN
214  CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip), &
215  sx, r)
216  END IF
217  END DO
218 
219 !$OMP END PARALLEL DO
220 
221  END SUBROUTINE xlda_hole_t_c_lr_lda_calc
222 
223 ! **************************************************************************************************
224 !> \brief low level routine
225 !> \param order ...
226 !> \param rho ...
227 !> \param e_0 ...
228 !> \param e_rho ...
229 !> \param sx ...
230 !> \param R ...
231 !> \par History
232 !> 12.2008 created [Manuel Guidon]
233 !> \author Manuel Guidon
234 ! **************************************************************************************************
235  SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, &
236  sx, R)
237  INTEGER, INTENT(IN) :: order
238  REAL(kind=dp), INTENT(IN) :: rho
239  REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho
240  REAL(kind=dp), INTENT(IN) :: sx, r
241 
242  REAL(kind=dp) :: t1, t12, t14, t15, t19, t2, t22, t23, &
243  t24, t25, t3, t32, t33, t36, t4, t41, &
244  t46, t5, t6, t62, t64, t67, t68, t7, &
245  t82, t86, t9, t91, t95
246 
247  IF (order >= 0) THEN
248  t1 = rho**2
249  t2 = t1*pi
250  t3 = 3**(0.1e1_dp/0.3e1_dp)
251  t4 = pi**2
252  t5 = t4*rho
253  t6 = t5**(0.1e1_dp/0.3e1_dp)
254  t7 = t6**2
255  t9 = t3/t7
256  t12 = log(r*t3*t6)
257  t14 = r**2
258  t15 = t14**2
259  t19 = 0.1e1_dp/d
260  t22 = t3**2
261  t23 = t22*t7
262  t24 = d*t14*t23
263  t25 = exp(-t24)
264  t32 = 9 + 4*a*t14*t23
265  t33 = log(t32)
266  t36 = d**2
267  t41 = expint(1, t24)
268  t46 = 0.1e1_dp/t36
269  t62 = log(0.2e1_dp)
270  t64 = log(a)
271  t67 = a*t12 + 0.3e1_dp/0.2e1_dp*e*t15*t3*t6*t5*t19*t25 &
272  - a*t33/0.2e1_dp + e/t36/d*t25 + a*t41/0.2e1_dp + e*t14 &
273  *t22*t7*t46*t25 + b*t19*t25/0.2e1_dp + c*t46*t25/0.2e1_dp &
274  + c*t14*t22*t7*t19*t25/0.2e1_dp + a*t62 + a*t64 &
275  /0.2e1_dp
276  t68 = t9*t67
277  e_0 = e_0 + (0.2e1_dp/0.3e1_dp*t2*t68)*sx
278  END IF
279  IF (order >= 1 .OR. order == -1) THEN
280  t82 = a/rho
281  t86 = t4**2
282  t91 = a**2
283  t95 = 0.1e1_dp/t6*t4
284  e_rho = e_rho + (0.4e1_dp/0.3e1_dp*rho*pi*t68 - 0.4e1_dp/0.9e1_dp*t1*t4*pi &
285  *t3/t7/t5*t67 + 0.2e1_dp/0.3e1_dp*t2*t9*(t82/0.3e1_dp - &
286  0.3e1_dp*e*t15*t14*t86*rho*t25 - 0.4e1_dp/0.3e1_dp*t91*t14 &
287  *t22*t95/t32 - t82*t25/0.3e1_dp - b*t14*t22*t95*t25 &
288  /0.3e1_dp - c*t15*t3*t6*t4*t25))*sx
289  END IF
290 
291  END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0
292 
293 ! **************************************************************************************************
294 !> \brief evaluates the truncated lsd exchange hole. Calls the lda routine and
295 !> applies spin scaling relation
296 !> \param rho_set the density where you want to evaluate the functional
297 !> \param deriv_set place where to store the functional derivatives (they are
298 !> added to the derivatives)
299 !> \param order degree of the derivative that should be evaluated,
300 !> if positive all the derivatives up to the given degree are evaluated,
301 !> if negative only the given degree is calculated
302 !> \param params input parameters (scaling, cutoff_radius)
303 !> \par History
304 !> 12.2008 created [Manuel Guidon]
305 !> \author Manuel Guidon
306 ! **************************************************************************************************
307  SUBROUTINE xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
308 
309  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
310  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
311  INTEGER, INTENT(IN) :: order
312  TYPE(section_vals_type), POINTER :: params
313 
314  CHARACTER(len=*), PARAMETER :: routinen = 'xlda_hole_t_c_lr_lsd_eval'
315 
316  INTEGER :: handle, npoints
317  INTEGER, DIMENSION(2, 3) :: bo
318  REAL(kind=dp) :: epsilon_rho, r, sx
319  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
320  POINTER :: dummy, e_0, e_rhoa, e_rhob, rhoa, rhob
321  TYPE(xc_derivative_type), POINTER :: deriv
322 
323  CALL timeset(routinen, handle)
324 
325  CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
326  CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
327 
328  CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
329  local_bounds=bo, rho_cutoff=epsilon_rho)
330  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
331 
332  dummy => rhoa
333 
334  e_0 => dummy
335  e_rhoa => dummy
336  e_rhob => dummy
337 
338  IF (order >= 0) THEN
339  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
340  allocate_deriv=.true.)
341  CALL xc_derivative_get(deriv, deriv_data=e_0)
342  END IF
343  IF (order >= 1 .OR. order == -1) THEN
344  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
345  allocate_deriv=.true.)
346  CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
347  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
348  allocate_deriv=.true.)
349  CALL xc_derivative_get(deriv, deriv_data=e_rhob)
350  END IF
351  IF (order > 1 .OR. order < -1) THEN
352  cpabort("derivatives bigger than 2 not implemented")
353  END IF
354  IF (r == 0.0_dp) THEN
355  cpabort("Cutoff_Radius 0.0 not implemented")
356  END IF
357 
358 !$OMP PARALLEL DEFAULT(NONE) &
359 !$OMP SHARED(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) &
360 !$OMP SHARED(sx, r,rhob, e_rhob)
361 
362  CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, &
363  e_0=e_0, e_rho=e_rhoa, &
364  epsilon_rho=epsilon_rho, &
365  sx=sx, r=r)
366 
367  CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, &
368  e_0=e_0, e_rho=e_rhob, &
369  epsilon_rho=epsilon_rho, &
370  sx=sx, r=r)
371 !$OMP END PARALLEL
372 
373  CALL timestop(handle)
374 
375  END SUBROUTINE xlda_hole_t_c_lr_lsd_eval
376 
377 ! **************************************************************************************************
378 !> \brief low level routine
379 !> \param npoints ...
380 !> \param order ...
381 !> \param rho ...
382 !> \param e_0 ...
383 !> \param e_rho ...
384 !> \param epsilon_rho ...
385 !> \param sx ...
386 !> \param R ...
387 !> \par History
388 !> 12.2008 created [Manuel Guidon]
389 !> \author Manuel Guidon
390 ! **************************************************************************************************
391  SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho, &
392  epsilon_rho, sx, R)
393 
394  INTEGER, INTENT(in) :: npoints, order
395  REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
396  REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
397 
398  INTEGER :: ip
399  REAL(dp) :: e_tmp, my_rho
400 
401 !$OMP DO
402 
403  DO ip = 1, npoints
404  my_rho = 2.0_dp*max(rho(ip), 0.0_dp)
405  IF (my_rho > epsilon_rho) THEN
406  e_tmp = 0.0_dp
407  CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp, e_rho(ip), &
408  sx, r)
409  e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
410  END IF
411  END DO
412 
413 !$OMP END DO
414 
415  END SUBROUTINE xlda_hole_t_c_lr_lsd_calc
416 END MODULE xc_xlda_hole_t_c_lr
417 
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
mathlib::expint
elemental impure real(dp) function, public expint(n, x)
computes the exponential integral En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,...
Definition: mathlib.F:1404
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition: xc_derivative_desc.F:26
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition: xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition: xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition: xc_derivative_set_types.F:66
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition: xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition: xc_derivative_types.F:103
xc_rho_cflags_types
contains the structure
Definition: xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition: xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition: xc_rho_set_types.F:281
xc_xlda_hole_t_c_lr
Calculates the lda exchange hole in a truncated coulomb potential. Can be used as longrange correctio...
Definition: xc_xlda_hole_t_c_lr.F:16
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lda_eval
subroutine, public xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
evaluates the truncated lda exchange hole
Definition: xc_xlda_hole_t_c_lr.F:127
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lda_info
subroutine, public xlda_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
returns various information on the functional
Definition: xc_xlda_hole_t_c_lr.F:67
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lsd_eval
subroutine, public xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
evaluates the truncated lsd exchange hole. Calls the lda routine and applies spin scaling relation
Definition: xc_xlda_hole_t_c_lr.F:308
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lda_calc_0
subroutine, public xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, sx, R)
low level routine
Definition: xc_xlda_hole_t_c_lr.F:237
xc_xlda_hole_t_c_lr::xlda_hole_t_c_lr_lsd_info
subroutine, public xlda_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
returns various information on the functional
Definition: xc_xlda_hole_t_c_lr.F:96