(git:08825b9)
Loading...
Searching...
No Matches
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-2026 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
17
20 USE kinds, ONLY: dp
21 USE mathconstants, ONLY: pi
22 USE mathlib, ONLY: expint
23 USE xc_derivative_desc, ONLY: deriv_rho,&
33#include "../base/base_uses.f90"
34
35 IMPLICIT NONE
36
37 PRIVATE
38
39! *** Global parameters ***
40
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
53CONTAINS
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 maximum degree of the derivative that should be evaluated
119!> \param params input parameters (scaling, cutoff_radius)
120!> \par History
121!> 12.2008 created [Manuel Guidon]
122!> \author Manuel Guidon
123! **************************************************************************************************
124 SUBROUTINE xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
125
126 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
127 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
128 INTEGER, INTENT(IN) :: order
129 TYPE(section_vals_type), POINTER :: params
130
131 CHARACTER(len=*), PARAMETER :: routinen = 'xlda_hole_t_c_lr_lda_eval'
132
133 INTEGER :: handle, npoints
134 INTEGER, DIMENSION(2, 3) :: bo
135 REAL(kind=dp) :: epsilon_rho, r, sx
136 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
137 POINTER :: dummy, e_0, e_rho, rho
138 TYPE(xc_derivative_type), POINTER :: deriv
139
140 CALL timeset(routinen, handle)
141
142 CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
143 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
144
145 CALL xc_rho_set_get(rho_set, rho=rho, &
146 local_bounds=bo, rho_cutoff=epsilon_rho)
147 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
148
149 dummy => rho
150
151 e_0 => dummy
152 e_rho => dummy
153
154 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
155 allocate_deriv=.true.)
156 CALL xc_derivative_get(deriv, deriv_data=e_0)
157 IF (order >= 1) THEN
158 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
159 allocate_deriv=.true.)
160 CALL xc_derivative_get(deriv, deriv_data=e_rho)
161 END IF
162 IF (order > 1) THEN
163 cpabort("derivatives bigger than 1 not implemented")
164 END IF
165
166 IF (r == 0.0_dp) THEN
167 cpabort("Cutoff_Radius 0.0 not implemented")
168 END IF
169 CALL xlda_hole_t_c_lr_lda_calc(npoints, order, rho=rho, &
170 e_0=e_0, e_rho=e_rho, &
171 epsilon_rho=epsilon_rho, &
172 sx=sx, r=r)
173
174 CALL timestop(handle)
175
176 END SUBROUTINE xlda_hole_t_c_lr_lda_eval
177
178! **************************************************************************************************
179!> \brief Call low level routine
180!> \param npoints ...
181!> \param order ...
182!> \param rho ...
183!> \param e_0 ...
184!> \param e_rho ...
185!> \param epsilon_rho ...
186!> \param sx ...
187!> \param R ...
188!> \par History
189!> 12.2008 created [Manuel Guidon]
190!> \author Manuel Guidon
191! **************************************************************************************************
192 SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho, &
193 epsilon_rho, sx, R)
194
195 INTEGER, INTENT(in) :: npoints, order
196 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
197 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
198
199 INTEGER :: ip
200 REAL(dp) :: my_rho
201
202!$OMP PARALLEL DO DEFAULT(NONE) &
203!$OMP SHARED(npoints, rho, epsilon_rho, order, e_0, e_rho) &
204!$OMP SHARED(sx, r) &
205!$OMP PRIVATE(ip, my_rho)
206
207 DO ip = 1, npoints
208 my_rho = max(rho(ip), 0.0_dp)
209 IF (my_rho > epsilon_rho) THEN
210 CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip), &
211 sx, r)
212 END IF
213 END DO
214
215!$OMP END PARALLEL DO
216
217 END SUBROUTINE xlda_hole_t_c_lr_lda_calc
218
219! **************************************************************************************************
220!> \brief low level routine
221!> \param order ...
222!> \param rho ...
223!> \param e_0 ...
224!> \param e_rho ...
225!> \param sx ...
226!> \param R ...
227!> \par History
228!> 12.2008 created [Manuel Guidon]
229!> \author Manuel Guidon
230! **************************************************************************************************
231 SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, &
232 sx, R)
233 INTEGER, INTENT(IN) :: order
234 REAL(kind=dp), INTENT(IN) :: rho
235 REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho
236 REAL(kind=dp), INTENT(IN) :: sx, r
237
238 REAL(kind=dp) :: t1, t12, t14, t15, t19, t2, t22, t23, &
239 t24, t25, t3, t32, t33, t36, t4, t41, &
240 t46, t5, t6, t62, t64, t67, t68, t7, &
241 t82, t86, t9, t91, t95
242
243 t1 = rho**2
244 t2 = t1*pi
245 t3 = 3**(0.1e1_dp/0.3e1_dp)
246 t4 = pi**2
247 t5 = t4*rho
248 t6 = t5**(0.1e1_dp/0.3e1_dp)
249 t7 = t6**2
250 t9 = t3/t7
251 t12 = log(r*t3*t6)
252 t14 = r**2
253 t15 = t14**2
254 t19 = 0.1e1_dp/d
255 t22 = t3**2
256 t23 = t22*t7
257 t24 = d*t14*t23
258 t25 = exp(-t24)
259 t32 = 9 + 4*a*t14*t23
260 t33 = log(t32)
261 t36 = d**2
262 t41 = expint(1, t24)
263 t46 = 0.1e1_dp/t36
264 t62 = log(0.2e1_dp)
265 t64 = log(a)
266 t67 = a*t12 + 0.3e1_dp/0.2e1_dp*e*t15*t3*t6*t5*t19*t25 &
267 - a*t33/0.2e1_dp + e/t36/d*t25 + a*t41/0.2e1_dp + e*t14 &
268 *t22*t7*t46*t25 + b*t19*t25/0.2e1_dp + c*t46*t25/0.2e1_dp &
269 + c*t14*t22*t7*t19*t25/0.2e1_dp + a*t62 + a*t64 &
270 /0.2e1_dp
271 t68 = t9*t67
272 e_0 = e_0 + (0.2e1_dp/0.3e1_dp*t2*t68)*sx
273 IF (order >= 1) THEN
274 t82 = a/rho
275 t86 = t4**2
276 t91 = a**2
277 t95 = 0.1e1_dp/t6*t4
278 e_rho = e_rho + (0.4e1_dp/0.3e1_dp*rho*pi*t68 - 0.4e1_dp/0.9e1_dp*t1*t4*pi &
279 *t3/t7/t5*t67 + 0.2e1_dp/0.3e1_dp*t2*t9*(t82/0.3e1_dp - &
280 0.3e1_dp*e*t15*t14*t86*rho*t25 - 0.4e1_dp/0.3e1_dp*t91*t14 &
281 *t22*t95/t32 - t82*t25/0.3e1_dp - b*t14*t22*t95*t25 &
282 /0.3e1_dp - c*t15*t3*t6*t4*t25))*sx
283 END IF
284
285 END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0
286
287! **************************************************************************************************
288!> \brief evaluates the truncated lsd exchange hole. Calls the lda routine and
289!> applies spin scaling relation
290!> \param rho_set the density where you want to evaluate the functional
291!> \param deriv_set place where to store the functional derivatives (they are
292!> added to the derivatives)
293!> \param order maximum degree of the derivative that should be evaluated
294!> \param params input parameters (scaling, cutoff_radius)
295!> \par History
296!> 12.2008 created [Manuel Guidon]
297!> \author Manuel Guidon
298! **************************************************************************************************
299 SUBROUTINE xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
300
301 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
302 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
303 INTEGER, INTENT(IN) :: order
304 TYPE(section_vals_type), POINTER :: params
305
306 CHARACTER(len=*), PARAMETER :: routinen = 'xlda_hole_t_c_lr_lsd_eval'
307
308 INTEGER :: handle, npoints
309 INTEGER, DIMENSION(2, 3) :: bo
310 REAL(kind=dp) :: epsilon_rho, r, sx
311 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
312 POINTER :: dummy, e_0, e_rhoa, e_rhob, rhoa, rhob
313 TYPE(xc_derivative_type), POINTER :: deriv
314
315 CALL timeset(routinen, handle)
316
317 CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
318 CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=r)
319
320 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
321 local_bounds=bo, rho_cutoff=epsilon_rho)
322 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
323
324 dummy => rhoa
325
326 e_0 => dummy
327 e_rhoa => dummy
328 e_rhob => dummy
329
330 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
331 allocate_deriv=.true.)
332 CALL xc_derivative_get(deriv, deriv_data=e_0)
333 IF (order >= 1) THEN
334 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
335 allocate_deriv=.true.)
336 CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
337 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
338 allocate_deriv=.true.)
339 CALL xc_derivative_get(deriv, deriv_data=e_rhob)
340 END IF
341 IF (order > 1) THEN
342 cpabort("derivatives bigger than 2 not implemented")
343 END IF
344 IF (r == 0.0_dp) THEN
345 cpabort("Cutoff_Radius 0.0 not implemented")
346 END IF
347
348!$OMP PARALLEL DEFAULT(NONE) &
349!$OMP SHARED(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) &
350!$OMP SHARED(sx, r,rhob, e_rhob)
351
352 CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, &
353 e_0=e_0, e_rho=e_rhoa, &
354 epsilon_rho=epsilon_rho, &
355 sx=sx, r=r)
356
357 CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, &
358 e_0=e_0, e_rho=e_rhob, &
359 epsilon_rho=epsilon_rho, &
360 sx=sx, r=r)
361!$OMP END PARALLEL
362
363 CALL timestop(handle)
364
365 END SUBROUTINE xlda_hole_t_c_lr_lsd_eval
366
367! **************************************************************************************************
368!> \brief low level routine
369!> \param npoints ...
370!> \param order ...
371!> \param rho ...
372!> \param e_0 ...
373!> \param e_rho ...
374!> \param epsilon_rho ...
375!> \param sx ...
376!> \param R ...
377!> \par History
378!> 12.2008 created [Manuel Guidon]
379!> \author Manuel Guidon
380! **************************************************************************************************
381 SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho, &
382 epsilon_rho, sx, R)
383
384 INTEGER, INTENT(in) :: npoints, order
385 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
386 REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, r
387
388 INTEGER :: ip
389 REAL(dp) :: e_tmp, my_rho
390
391!$OMP DO
392
393 DO ip = 1, npoints
394 my_rho = 2.0_dp*max(rho(ip), 0.0_dp)
395 IF (my_rho > epsilon_rho) THEN
396 e_tmp = 0.0_dp
397 CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp, e_rho(ip), &
398 sx, r)
399 e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
400 END IF
401 END DO
402
403!$OMP END DO
404
405 END SUBROUTINE xlda_hole_t_c_lr_lsd_calc
406END MODULE xc_xlda_hole_t_c_lr
407
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
elemental impure real(dp) function, public expint(n, x, exp_mx)
computes the exponential integral En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,...
Definition mathlib.F:1531
Module with functions to handle derivative descriptors. derivative description are strings have the f...
integer, parameter, public deriv_rhob
integer, parameter, public deriv_rhoa
integer, parameter, public deriv_rho
represent a group ofunctional derivatives
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
contains the structure
contains the structure
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
Calculates the lda exchange hole in a truncated coulomb potential. Can be used as longrange correctio...
subroutine, public xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
evaluates the truncated lda exchange hole
subroutine, public xlda_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
returns various information on the functional
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
subroutine, public xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, sx, r)
low level routine
subroutine, public xlda_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
returns various information on the functional
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
represent a derivative of a functional
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation