(git:374b731)
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-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
16MODULE 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
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 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
416END 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
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:49
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78