Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates the Exchange contribution in the BEEF-vdW functional
10 : ! (Norskov Group)
11 : !> \par History
12 : !> 02.2014 created based on xc_xbecke88.F [rkoitz]
13 : !> \author rkoitz
14 : ! **************************************************************************************************
15 : MODULE xc_xbeef
16 :
17 : USE bibliography, ONLY: Wellendorff2012,&
18 : cite_reference
19 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
20 : USE input_section_types, ONLY: section_vals_type,&
21 : section_vals_val_get
22 : USE kinds, ONLY: dp
23 : USE mathconstants, ONLY: pi
24 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
25 : deriv_norm_drhoa,&
26 : deriv_norm_drhob,&
27 : deriv_rho,&
28 : deriv_rhoa,&
29 : deriv_rhob
30 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
31 : xc_dset_get_derivative
32 : USE xc_derivative_types, ONLY: xc_derivative_get,&
33 : xc_derivative_type
34 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
35 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
36 : xc_rho_set_type
37 : #include "../base/base_uses.f90"
38 :
39 : IMPLICIT NONE
40 : PRIVATE
41 :
42 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbeef'
44 :
45 : PUBLIC :: xbeef_lda_info, xbeef_lsd_info, xbeef_lda_eval, xbeef_lsd_eval
46 :
47 : !for beef functional pulled out of GPAW source, Feb 2014
48 : REAL(kind=dp), PARAMETER :: a(0:29) = (/1.516501714304992365356_dp, 0.441353209874497942611_dp, -0.091821352411060291887_dp, &
49 : -0.023527543314744041314_dp, 0.034188284548603550816_dp, 0.002411870075717384172_dp, &
50 : -0.014163813515916020766_dp, 0.000697589558149178113_dp, 0.009859205136982565273_dp, &
51 : -0.006737855050935187551_dp, -0.001573330824338589097_dp, 0.005036146253345903309_dp, &
52 : -0.002569472452841069059_dp, -0.000987495397608761146_dp, 0.002033722894696920677_dp, &
53 : -0.000801871884834044583_dp, -0.000668807872347525591_dp, 0.001030936331268264214_dp, &
54 : -0.000367383865990214423_dp, -0.000421363539352619543_dp, 0.000576160799160517858_dp, &
55 : -0.000083465037349510408_dp, -0.000445844758523195788_dp, 0.000460129009232047457_dp, &
56 : -0.000005231775398304339_dp, -0.000423957047149510404_dp, 0.000375019067938866537_dp, &
57 : 0.000021149381251344578_dp, -0.000190491156503997170_dp, 0.000073843624209823442_dp/)
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief return various information on the functional
62 : !> \param reference string with the reference of the actual functional
63 : !> \param shortform string with the shortform of the functional name
64 : !> \param needs the components needed by this functional are set to
65 : !> true (does not set the unneeded components to false)
66 : !> \param max_deriv ...
67 : !> \par History
68 : !> 02.2014 created
69 : !> \author rkoitz
70 : ! **************************************************************************************************
71 :
72 23 : SUBROUTINE xbeef_lda_info(reference, shortform, needs, max_deriv)
73 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
74 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
75 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
76 :
77 23 : IF (PRESENT(reference)) THEN
78 1 : reference = "Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LDA}"
79 : END IF
80 23 : IF (PRESENT(shortform)) THEN
81 1 : shortform = "Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LDA}"
82 : END IF
83 23 : IF (PRESENT(needs)) THEN
84 22 : needs%rho = .TRUE.
85 22 : needs%rho_1_3 = .TRUE.
86 22 : needs%norm_drho = .TRUE.
87 : END IF
88 23 : IF (PRESENT(max_deriv)) max_deriv = 1
89 :
90 23 : END SUBROUTINE xbeef_lda_info
91 :
92 : ! **************************************************************************************************
93 : !> \brief return various information on the functional
94 : !> \param reference string with the reference of the actual functional
95 : !> \param shortform string with the shortform of the functional name
96 : !> \param needs the components needed by this functional are set to
97 : !> true (does not set the unneeded components to false)
98 : !> \param max_deriv ...
99 : !> \par History
100 : !> 02.2014 created
101 : !> \author rkoitz
102 : ! **************************************************************************************************
103 0 : SUBROUTINE xbeef_lsd_info(reference, shortform, needs, max_deriv)
104 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
105 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
106 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
107 :
108 0 : IF (PRESENT(reference)) THEN
109 0 : reference = "Wellendorff, J. et al., Phys. Rev. B 85, 235149 (2012) {LSD}"
110 : END IF
111 0 : IF (PRESENT(shortform)) THEN
112 0 : shortform = "Exchange Contribution to BEEF-vdW Functional (Wellendorff, 2012) {LSD}"
113 : END IF
114 0 : IF (PRESENT(needs)) THEN
115 0 : needs%rho_spin = .TRUE.
116 0 : needs%rho_spin_1_3 = .TRUE.
117 0 : needs%norm_drho_spin = .TRUE.
118 : END IF
119 0 : IF (PRESENT(max_deriv)) max_deriv = 1
120 :
121 0 : END SUBROUTINE xbeef_lsd_info
122 :
123 : ! **************************************************************************************************
124 : !> \brief evaluates the beef exchange functional for lda
125 : !> \param rho_set the density where you want to evaluate the functional
126 : !> \param deriv_set place where to store the functional derivatives (they are
127 : !> added to the derivatives)
128 : !> \param grad_deriv degree of the derivative that should be evaluated,
129 : !> if positive all the derivatives up to the given degree are evaluated,
130 : !> if negative only the given degree is calculated
131 : !> \param xbeef_params input parameters (scaling)
132 : !> \par History
133 : !> 02.2014 created
134 : !> \author rkoitz
135 : ! **************************************************************************************************
136 54 : SUBROUTINE xbeef_lda_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
137 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
138 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
139 : INTEGER, INTENT(in) :: grad_deriv
140 : TYPE(section_vals_type), POINTER :: xbeef_params
141 :
142 : CHARACTER(len=*), PARAMETER :: routineN = 'xbeef_lda_eval'
143 :
144 : INTEGER :: handle, npoints
145 : INTEGER, DIMENSION(2, 3) :: bo
146 : REAL(kind=dp) :: epsilon_rho, sx
147 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
148 18 : POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
149 18 : rho, rho_1_3
150 : TYPE(xc_derivative_type), POINTER :: deriv
151 :
152 18 : CALL timeset(routineN, handle)
153 :
154 18 : CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
155 :
156 18 : CALL cite_reference(Wellendorff2012)
157 :
158 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
159 18 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
160 18 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
161 :
162 18 : dummy => rho
163 :
164 18 : e_0 => dummy
165 18 : e_rho => dummy
166 18 : e_ndrho => dummy
167 :
168 18 : IF (grad_deriv >= 0) THEN
169 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
170 18 : allocate_deriv=.TRUE.)
171 18 : CALL xc_derivative_get(deriv, deriv_data=e_0)
172 : END IF
173 18 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
174 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
175 18 : allocate_deriv=.TRUE.)
176 18 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
177 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
178 18 : allocate_deriv=.TRUE.)
179 18 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
180 : END IF
181 18 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
182 0 : CPABORT("derivatives greater than 1 not implemented")
183 : END IF
184 :
185 : !$OMP PARALLEL DEFAULT(NONE) &
186 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
187 : !$OMP SHARED(e_ndrho) &
188 : !$OMP SHARED( grad_deriv, npoints) &
189 18 : !$OMP SHARED(epsilon_rho,sx)
190 : CALL xbeef_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
191 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
192 : grad_deriv=grad_deriv, &
193 : npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
194 : !$OMP END PARALLEL
195 :
196 18 : CALL timestop(handle)
197 18 : END SUBROUTINE xbeef_lda_eval
198 :
199 : ! **************************************************************************************************
200 : !> \brief evaluates the beef exchange functional for lda
201 : !> \param rho the density where you want to evaluate the functional
202 : !> \param rho_1_3 ...
203 : !> \param norm_drho ...
204 : !> \param e_0 ...
205 : !> \param e_rho ...
206 : !> \param e_ndrho ...
207 : !> \param grad_deriv degree of the derivative that should be evaluated,
208 : !> if positive all the derivatives up to the given degree are evaluated,
209 : !> if negative only the given degree is calculated
210 : !> \param npoints ...
211 : !> \param epsilon_rho ...
212 : !> \param sx scaling-parameter for exchange
213 : !> \par History
214 : !> 02.2014 created based on xc_xbecke88
215 : !> \author rkoitz
216 : ! **************************************************************************************************
217 18 : SUBROUTINE xbeef_lda_calc(rho, rho_1_3, norm_drho, &
218 18 : e_0, e_rho, e_ndrho, &
219 : grad_deriv, npoints, epsilon_rho, sx)
220 : INTEGER, INTENT(in) :: npoints, grad_deriv
221 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho, e_rho, e_0
222 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
223 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
224 :
225 : INTEGER, PARAMETER :: m = 30
226 :
227 : INTEGER :: i, ii
228 : REAL(kind=dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
229 : epsilon_rho43, k_lda, kf, my_rho, &
230 : my_rho_1_3, s, s2, t, t3
231 : REAL(kind=dp), DIMENSION(0:29) :: de_leg, e_leg
232 :
233 : !energies, first and second derivatives from legendre pol
234 :
235 18 : kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
236 18 : k_lda = -(2.0_dp/(2.0_dp**(4._dp/3._dp)))*(3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*pi))**(1.0_dp/3.0_dp)
237 : !exchange energy density of the uniform electron gas, non spin-polarized
238 :
239 18 : epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
240 :
241 : !$OMP DO
242 : DO ii = 1, npoints
243 :
244 1417176 : my_rho = rho(ii)
245 :
246 1417176 : IF (my_rho > epsilon_rho) THEN
247 1417176 : my_rho_1_3 = rho_1_3(ii)
248 :
249 1417176 : e_ueg = k_lda*my_rho*my_rho_1_3
250 1417176 : e_ueg_drho = (4.0_dp/3.0_dp)*k_lda*my_rho_1_3
251 :
252 1417176 : t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
253 :
254 1417176 : s = norm_drho(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
255 1417176 : s2 = s**2
256 1417176 : t = 2.0_dp*s2/(4.0_dp + s2) - 1.0_dp
257 :
258 1417176 : IF (grad_deriv >= 0) THEN !asking for pure e evaluation or also derivatives
259 1417176 : e_leg(0) = 1 !first legendre pol
260 1417176 : e_leg(1) = t !second legendre pol
261 : END IF
262 :
263 1417176 : IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
264 1417176 : de_leg(0) = 0
265 1417176 : de_leg(1) = 1
266 1417176 : dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
267 1417176 : ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
268 1417176 : ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
269 : END IF
270 :
271 41098104 : DO i = 2, m - 1 !LEGENDRE PART
272 39680928 : e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
273 : !taken from quantum espresso beef library.
274 :
275 41098104 : IF (ABS(grad_deriv) >= 1) THEN !first derivative
276 : !the zero-derivatives need to be available for the first deriv.
277 39680928 : de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
278 : END IF
279 : END DO
280 :
281 : !NO DERIVATIVE
282 1417176 : IF (grad_deriv >= 0) THEN
283 : !add the scaled legendre linear combination to e_0
284 43932456 : e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
285 : END IF
286 :
287 : !FIRST DERIVATIVE
288 1417176 : IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
289 86447736 : e_rho(ii) = e_rho(ii) + (SUM(e_leg*a)*e_ueg_drho + SUM(de_leg*a)*dt*ds_rho*e_ueg)*sx
290 43932456 : e_ndrho(ii) = e_ndrho(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
291 : END IF
292 :
293 : END IF
294 : END DO
295 :
296 : !$OMP END DO
297 :
298 18 : END SUBROUTINE xbeef_lda_calc
299 :
300 : ! **************************************************************************************************
301 : !> \brief evaluates the beef 88 exchange functional for lsd
302 : !> \param rho_set ...
303 : !> \param deriv_set ...
304 : !> \param grad_deriv ...
305 : !> \param xbeef_params ...
306 : !> \par History
307 : !> 2/2014 rkoitz [created based on Becke 88]
308 : !> \author rkoitz
309 : ! **************************************************************************************************
310 0 : SUBROUTINE xbeef_lsd_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
311 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
312 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
313 : INTEGER, INTENT(in) :: grad_deriv
314 : TYPE(section_vals_type), POINTER :: xbeef_params
315 :
316 : CHARACTER(len=*), PARAMETER :: routineN = 'xbeef_lsd_eval'
317 :
318 : INTEGER :: handle, i, ispin, npoints
319 : INTEGER, DIMENSION(2, 3) :: bo
320 : REAL(kind=dp) :: epsilon_rho, sx
321 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
322 0 : POINTER :: dummy, e_0
323 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, norm_drho, rho, rho_1_3
324 : TYPE(xc_derivative_type), POINTER :: deriv
325 :
326 0 : CALL timeset(routineN, handle)
327 :
328 0 : CALL cite_reference(Wellendorff2012)
329 :
330 0 : NULLIFY (deriv)
331 0 : DO i = 1, 2
332 0 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
333 : END DO
334 :
335 0 : CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
336 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
337 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
338 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
339 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
340 0 : local_bounds=bo)
341 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
342 :
343 0 : dummy => rho(1)%array
344 :
345 0 : e_0 => dummy
346 0 : DO i = 1, 2
347 0 : e_rho(i)%array => dummy
348 0 : e_ndrho(i)%array => dummy
349 : END DO
350 :
351 0 : IF (grad_deriv >= 0) THEN
352 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
353 0 : allocate_deriv=.TRUE.)
354 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
355 : END IF
356 0 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
357 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
358 0 : allocate_deriv=.TRUE.)
359 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
360 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
361 0 : allocate_deriv=.TRUE.)
362 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
363 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
364 0 : allocate_deriv=.TRUE.)
365 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
366 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
367 0 : allocate_deriv=.TRUE.)
368 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
369 : END IF
370 0 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
371 0 : CPABORT("derivatives greater than 1 not implemented")
372 : END IF
373 :
374 0 : DO ispin = 1, 2
375 :
376 : !$OMP PARALLEL DEFAULT(NONE) &
377 : !$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
378 : !$OMP SHARED(e_rho, e_ndrho) &
379 : !$OMP SHARED(grad_deriv, npoints) &
380 0 : !$OMP SHARED(epsilon_rho, sx)
381 :
382 : CALL xbeef_lsd_calc( &
383 : rho_spin=rho(ispin)%array, &
384 : rho_1_3_spin=rho_1_3(ispin)%array, &
385 : norm_drho_spin=norm_drho(ispin)%array, &
386 : e_0=e_0, e_rho_spin=e_rho(ispin)%array, &
387 : e_ndrho_spin=e_ndrho(ispin)%array, &
388 : grad_deriv=grad_deriv, npoints=npoints, &
389 : epsilon_rho=epsilon_rho, sx=sx)
390 :
391 : !$OMP END PARALLEL
392 :
393 : END DO
394 :
395 0 : CALL timestop(handle)
396 :
397 0 : END SUBROUTINE xbeef_lsd_eval
398 : ! **************************************************************************************************
399 : !> \brief low level calculation of the beef exchange functional for lsd
400 : !> \param rho_spin alpha or beta spin density
401 : !> \param rho_1_3_spin rho_spin**(1./3.)
402 : !> \param norm_drho_spin || grad rho_spin ||
403 : !> \param e_0 adds to it the local value of the functional
404 : !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
405 : !> named where the * is. Everything wrt. to the spin of the arguments.
406 : !> \param e_ndrho_spin ...
407 : !> \param grad_deriv ...
408 : !> \param npoints ...
409 : !> \param epsilon_rho ...
410 : !> \param sx scaling-parameter for exchange
411 : !> \par History
412 : !> 02.2014 created based on Becke88
413 : !> \author rkoitz
414 : ! **************************************************************************************************
415 0 : SUBROUTINE xbeef_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
416 : e_rho_spin, e_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
417 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
418 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin
419 : INTEGER, INTENT(in) :: grad_deriv, npoints
420 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
421 :
422 : INTEGER, PARAMETER :: m = 30
423 :
424 : INTEGER :: i, ii
425 : REAL(kind=dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
426 : epsilon_rho43, k_lsd, kf, &
427 : my_epsilon_rho, my_rho, my_rho_1_3, s, &
428 : s2, t, t3
429 : REAL(kind=dp), DIMENSION(0:29) :: de_leg, e_leg
430 :
431 : !energies and first derivatives from legendre pol
432 :
433 0 : kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
434 0 : k_lsd = (3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*pi))**(1.0_dp/3.0_dp)
435 : !exchange energy density of the uniform electron gas, spin-polarized
436 :
437 0 : my_epsilon_rho = 0.5_dp*epsilon_rho
438 0 : epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
439 :
440 0 : !$OMP DO
441 : DO ii = 1, npoints
442 0 : my_rho = rho_spin(ii)
443 :
444 0 : IF (my_rho > epsilon_rho) THEN
445 0 : my_rho_1_3 = rho_1_3_spin(ii)
446 :
447 0 : e_ueg = k_lsd*my_rho*my_rho_1_3
448 0 : e_ueg_drho = (4.0_dp/3.0_dp)*k_lsd*my_rho_1_3
449 :
450 0 : t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
451 :
452 0 : s = norm_drho_spin(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
453 0 : s2 = s**2
454 0 : t = 2.0_dp*s**2/(4.0_dp + s**2) - 1.0_dp
455 :
456 0 : IF (grad_deriv >= 0) THEN !asking for pure e evaluation or also derivatives
457 0 : e_leg(0) = 1 !first legendre pol
458 0 : e_leg(1) = t !second legendre pol
459 : END IF
460 :
461 0 : IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
462 0 : de_leg(0) = 0
463 0 : de_leg(1) = 1
464 0 : dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
465 0 : ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
466 0 : ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
467 : END IF
468 :
469 0 : DO i = 2, m - 1 !LEGENDRE PART
470 0 : e_leg(i) = 2.*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
471 : !taken from quantum espresso beef library.
472 :
473 0 : IF (ABS(grad_deriv) >= 1) THEN !first derivative
474 : !the zero-derivatives need to be available for the first deriv.
475 0 : de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
476 : END IF
477 :
478 : END DO
479 :
480 : !NO DERIVATIVE
481 0 : IF (grad_deriv >= 0) THEN
482 : !add the scaled legendre linear combination to e_0
483 0 : e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
484 : END IF
485 :
486 : !FIRST DERIVATIVE
487 0 : IF ((grad_deriv >= 1) .OR. (grad_deriv == -1)) THEN !asking for first derivative or higher
488 0 : e_rho_spin(ii) = e_rho_spin(ii) + (SUM(e_leg*a)*e_ueg_drho + SUM(de_leg*a)*dt*ds_rho*e_ueg)*sx
489 0 : e_ndrho_spin(ii) = e_ndrho_spin(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
490 : END IF
491 : END IF
492 : END DO
493 : !$OMP END DO
494 :
495 0 : END SUBROUTINE xbeef_lsd_calc
496 :
497 : END MODULE xc_xbeef
498 :
|
