Line data Source code
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 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 maximum degree of the derivative that should be evaluated
129 : !> \param xbeef_params input parameters (scaling)
130 : !> \par History
131 : !> 02.2014 created
132 : !> \author rkoitz
133 : ! **************************************************************************************************
134 54 : SUBROUTINE xbeef_lda_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
135 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
136 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
137 : INTEGER, INTENT(in) :: grad_deriv
138 : TYPE(section_vals_type), POINTER :: xbeef_params
139 :
140 : CHARACTER(len=*), PARAMETER :: routineN = 'xbeef_lda_eval'
141 :
142 : INTEGER :: handle, npoints
143 : INTEGER, DIMENSION(2, 3) :: bo
144 : REAL(kind=dp) :: epsilon_rho, sx
145 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
146 18 : POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
147 18 : rho, rho_1_3
148 : TYPE(xc_derivative_type), POINTER :: deriv
149 :
150 18 : CALL timeset(routineN, handle)
151 :
152 18 : CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
153 :
154 18 : CALL cite_reference(Wellendorff2012)
155 :
156 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
157 18 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
158 18 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
159 :
160 18 : dummy => rho
161 :
162 18 : e_0 => dummy
163 18 : e_rho => dummy
164 18 : e_ndrho => dummy
165 :
166 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
167 18 : allocate_deriv=.TRUE.)
168 18 : CALL xc_derivative_get(deriv, deriv_data=e_0)
169 18 : IF (grad_deriv >= 1) THEN
170 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
171 18 : allocate_deriv=.TRUE.)
172 18 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
173 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
174 18 : allocate_deriv=.TRUE.)
175 18 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
176 : END IF
177 18 : IF (grad_deriv > 1) THEN
178 0 : CPABORT("derivatives greater than 1 not implemented")
179 : END IF
180 :
181 : !$OMP PARALLEL DEFAULT(NONE) &
182 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
183 : !$OMP SHARED(e_ndrho) &
184 : !$OMP SHARED( grad_deriv, npoints) &
185 18 : !$OMP SHARED(epsilon_rho,sx)
186 : CALL xbeef_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
187 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
188 : grad_deriv=grad_deriv, &
189 : npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
190 : !$OMP END PARALLEL
191 :
192 18 : CALL timestop(handle)
193 18 : END SUBROUTINE xbeef_lda_eval
194 :
195 : ! **************************************************************************************************
196 : !> \brief evaluates the beef exchange functional for lda
197 : !> \param rho the density where you want to evaluate the functional
198 : !> \param rho_1_3 ...
199 : !> \param norm_drho ...
200 : !> \param e_0 ...
201 : !> \param e_rho ...
202 : !> \param e_ndrho ...
203 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
204 : !> \param npoints ...
205 : !> \param epsilon_rho ...
206 : !> \param sx scaling-parameter for exchange
207 : !> \par History
208 : !> 02.2014 created based on xc_xbecke88
209 : !> \author rkoitz
210 : ! **************************************************************************************************
211 18 : SUBROUTINE xbeef_lda_calc(rho, rho_1_3, norm_drho, &
212 18 : e_0, e_rho, e_ndrho, &
213 : grad_deriv, npoints, epsilon_rho, sx)
214 : INTEGER, INTENT(in) :: npoints, grad_deriv
215 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho, e_rho, e_0
216 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
217 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
218 :
219 : INTEGER, PARAMETER :: m = 30
220 :
221 : INTEGER :: i, ii
222 : REAL(kind=dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
223 : epsilon_rho43, k_lda, kf, my_rho, &
224 : my_rho_1_3, s, s2, t, t3
225 : REAL(kind=dp), DIMENSION(0:29) :: de_leg, e_leg
226 :
227 : !energies, first and second derivatives from legendre pol
228 :
229 18 : kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
230 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)
231 : !exchange energy density of the uniform electron gas, non spin-polarized
232 :
233 18 : epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
234 :
235 : !$OMP DO
236 : DO ii = 1, npoints
237 :
238 1417176 : my_rho = rho(ii)
239 :
240 1417176 : IF (my_rho > epsilon_rho) THEN
241 1417176 : my_rho_1_3 = rho_1_3(ii)
242 :
243 1417176 : e_ueg = k_lda*my_rho*my_rho_1_3
244 1417176 : e_ueg_drho = (4.0_dp/3.0_dp)*k_lda*my_rho_1_3
245 :
246 1417176 : t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
247 :
248 1417176 : s = norm_drho(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
249 1417176 : s2 = s**2
250 1417176 : t = 2.0_dp*s2/(4.0_dp + s2) - 1.0_dp
251 :
252 1417176 : e_leg(0) = 1 !first legendre pol
253 1417176 : e_leg(1) = t !second legendre pol
254 :
255 1417176 : IF ((grad_deriv >= 1)) THEN !asking for first derivative or higher
256 1417176 : de_leg(0) = 0
257 1417176 : de_leg(1) = 1
258 1417176 : dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
259 1417176 : ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
260 1417176 : ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
261 : END IF
262 :
263 41098104 : DO i = 2, m - 1 !LEGENDRE PART
264 39680928 : e_leg(i) = 2._dp*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
265 : !taken from quantum espresso beef library.
266 :
267 41098104 : IF (grad_deriv >= 1) THEN !first derivative
268 : !the zero-derivatives need to be available for the first deriv.
269 39680928 : de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
270 : END IF
271 : END DO
272 :
273 : !add the scaled legendre linear combination to e_0
274 43932456 : e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
275 :
276 : !FIRST DERIVATIVE
277 1417176 : IF ((grad_deriv >= 1)) THEN !asking for first derivative or higher
278 86447736 : e_rho(ii) = e_rho(ii) + (SUM(e_leg*a)*e_ueg_drho + SUM(de_leg*a)*dt*ds_rho*e_ueg)*sx
279 43932456 : e_ndrho(ii) = e_ndrho(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
280 : END IF
281 :
282 : END IF
283 : END DO
284 :
285 : !$OMP END DO
286 :
287 18 : END SUBROUTINE xbeef_lda_calc
288 :
289 : ! **************************************************************************************************
290 : !> \brief evaluates the beef 88 exchange functional for lsd
291 : !> \param rho_set ...
292 : !> \param deriv_set ...
293 : !> \param grad_deriv ...
294 : !> \param xbeef_params ...
295 : !> \par History
296 : !> 2/2014 rkoitz [created based on Becke 88]
297 : !> \author rkoitz
298 : ! **************************************************************************************************
299 0 : SUBROUTINE xbeef_lsd_eval(rho_set, deriv_set, grad_deriv, xbeef_params)
300 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
301 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
302 : INTEGER, INTENT(in) :: grad_deriv
303 : TYPE(section_vals_type), POINTER :: xbeef_params
304 :
305 : CHARACTER(len=*), PARAMETER :: routineN = 'xbeef_lsd_eval'
306 :
307 : INTEGER :: handle, i, ispin, npoints
308 : INTEGER, DIMENSION(2, 3) :: bo
309 : REAL(kind=dp) :: epsilon_rho, sx
310 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
311 0 : POINTER :: dummy, e_0
312 0 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, norm_drho, rho, rho_1_3
313 : TYPE(xc_derivative_type), POINTER :: deriv
314 :
315 0 : CALL timeset(routineN, handle)
316 :
317 0 : CALL cite_reference(Wellendorff2012)
318 :
319 0 : NULLIFY (deriv)
320 0 : DO i = 1, 2
321 0 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
322 : END DO
323 :
324 0 : CALL section_vals_val_get(xbeef_params, "scale_x", r_val=sx)
325 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
326 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
327 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
328 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
329 0 : local_bounds=bo)
330 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
331 :
332 0 : dummy => rho(1)%array
333 :
334 0 : e_0 => dummy
335 0 : DO i = 1, 2
336 0 : e_rho(i)%array => dummy
337 0 : e_ndrho(i)%array => dummy
338 : END DO
339 :
340 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
341 0 : allocate_deriv=.TRUE.)
342 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
343 0 : IF (grad_deriv >= 1) THEN
344 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
345 0 : allocate_deriv=.TRUE.)
346 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
347 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
348 0 : allocate_deriv=.TRUE.)
349 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
350 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
351 0 : allocate_deriv=.TRUE.)
352 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
353 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
354 0 : allocate_deriv=.TRUE.)
355 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
356 : END IF
357 0 : IF (grad_deriv > 1) THEN
358 0 : CPABORT("derivatives greater than 1 not implemented")
359 : END IF
360 :
361 0 : DO ispin = 1, 2
362 :
363 : !$OMP PARALLEL DEFAULT(NONE) &
364 : !$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
365 : !$OMP SHARED(e_rho, e_ndrho) &
366 : !$OMP SHARED(grad_deriv, npoints) &
367 0 : !$OMP SHARED(epsilon_rho, sx)
368 :
369 : CALL xbeef_lsd_calc( &
370 : rho_spin=rho(ispin)%array, &
371 : rho_1_3_spin=rho_1_3(ispin)%array, &
372 : norm_drho_spin=norm_drho(ispin)%array, &
373 : e_0=e_0, e_rho_spin=e_rho(ispin)%array, &
374 : e_ndrho_spin=e_ndrho(ispin)%array, &
375 : grad_deriv=grad_deriv, npoints=npoints, &
376 : epsilon_rho=epsilon_rho, sx=sx)
377 :
378 : !$OMP END PARALLEL
379 :
380 : END DO
381 :
382 0 : CALL timestop(handle)
383 :
384 0 : END SUBROUTINE xbeef_lsd_eval
385 : ! **************************************************************************************************
386 : !> \brief low level calculation of the beef exchange functional for lsd
387 : !> \param rho_spin alpha or beta spin density
388 : !> \param rho_1_3_spin rho_spin**(1./3.)
389 : !> \param norm_drho_spin || grad rho_spin ||
390 : !> \param e_0 adds to it the local value of the functional
391 : !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
392 : !> named where the * is. Everything wrt. to the spin of the arguments.
393 : !> \param e_ndrho_spin ...
394 : !> \param grad_deriv ...
395 : !> \param npoints ...
396 : !> \param epsilon_rho ...
397 : !> \param sx scaling-parameter for exchange
398 : !> \par History
399 : !> 02.2014 created based on Becke88
400 : !> \author rkoitz
401 : ! **************************************************************************************************
402 0 : SUBROUTINE xbeef_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
403 : e_rho_spin, e_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
404 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
405 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin
406 : INTEGER, INTENT(in) :: grad_deriv, npoints
407 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
408 :
409 : INTEGER, PARAMETER :: m = 30
410 :
411 : INTEGER :: i, ii
412 : REAL(kind=dp) :: ds_ndrho, ds_rho, dt, e_ueg, e_ueg_drho, &
413 : epsilon_rho43, k_lsd, kf, &
414 : my_epsilon_rho, my_rho, my_rho_1_3, s, &
415 : s2, t, t3
416 : REAL(kind=dp), DIMENSION(0:29) :: de_leg, e_leg
417 :
418 : !energies and first derivatives from legendre pol
419 :
420 0 : kf = (3.0_dp*pi**2)**(1.0_dp/3.0_dp) !only constants, without n^1/3
421 0 : k_lsd = (3.0_dp/2.0_dp)*(3.0_dp/(4.0_dp*pi))**(1.0_dp/3.0_dp)
422 : !exchange energy density of the uniform electron gas, spin-polarized
423 :
424 0 : my_epsilon_rho = 0.5_dp*epsilon_rho
425 0 : epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
426 :
427 0 : !$OMP DO
428 : DO ii = 1, npoints
429 0 : my_rho = rho_spin(ii)
430 :
431 0 : IF (my_rho > epsilon_rho) THEN
432 0 : my_rho_1_3 = rho_1_3_spin(ii)
433 :
434 0 : e_ueg = k_lsd*my_rho*my_rho_1_3
435 0 : e_ueg_drho = (4.0_dp/3.0_dp)*k_lsd*my_rho_1_3
436 :
437 0 : t3 = my_rho_1_3*my_rho*2*kf !reduced gradient, denominator
438 :
439 0 : s = norm_drho_spin(ii)/MAX(t3, epsilon_rho43) !reduced gradient finally
440 0 : s2 = s**2
441 0 : t = 2.0_dp*s**2/(4.0_dp + s**2) - 1.0_dp
442 :
443 0 : e_leg(0) = 1 !first legendre pol
444 0 : e_leg(1) = t !second legendre pol
445 :
446 0 : IF ((grad_deriv >= 1)) THEN !asking for first derivative or higher
447 0 : de_leg(0) = 0
448 0 : de_leg(1) = 1
449 0 : dt = 4.0_dp*s/(4.0_dp + s2) - 4.0_dp*s*s2/(4.0_dp + s2)**2
450 0 : ds_rho = -(4.0_dp*s)/(3.0_dp*MAX(my_rho, epsilon_rho))
451 0 : ds_ndrho = 1.0_dp/(MAX(t3, epsilon_rho43))
452 : END IF
453 :
454 0 : DO i = 2, m - 1 !LEGENDRE PART
455 0 : e_leg(i) = 2._dp*(t)*e_leg(i - 1) - e_leg(i - 2) - ((t)*e_leg(i - 1) - e_leg(i - 2))/(REAL(i, KIND=dp))
456 : !taken from quantum espresso beef library.
457 :
458 0 : IF (grad_deriv >= 1) THEN !first derivative
459 : !the zero-derivatives need to be available for the first deriv.
460 0 : de_leg(i) = e_leg(i - 1)*i + de_leg(i - 1)*(t)
461 : END IF
462 :
463 : END DO
464 :
465 : !NO DERIVATIVE
466 : !add the scaled legendre linear combination to e_0
467 0 : e_0(ii) = e_0(ii) + SUM(e_leg*a)*e_ueg*sx
468 :
469 : !FIRST DERIVATIVE
470 0 : IF ((grad_deriv >= 1)) THEN !asking for first derivative or higher
471 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
472 0 : e_ndrho_spin(ii) = e_ndrho_spin(ii) + (SUM(de_leg*a)*dt*ds_ndrho*e_ueg)*sx
473 : END IF
474 : END IF
475 : END DO
476 : !$OMP END DO
477 :
478 0 : END SUBROUTINE xbeef_lsd_calc
479 :
480 : END MODULE xc_xbeef
481 :
|
