Line data Source code
1 : !!****m* m_paw_exactDC/m_paw_exactDC
2 : !! NAME
3 : !! m_paw_exactDC
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines related to the exact formula for the double counting.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2025-2026 ABINIT group
10 : !! These routines were inspired by K. Haule routines in embedded DMFT.
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : MODULE m_paw_exactDC
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 :
29 : implicit none
30 :
31 : private
32 :
33 : public :: compute_exactDC
34 :
35 : CONTAINS !========================================================================================
36 : !!***
37 :
38 : !----------------------------------------------------------------------
39 :
40 : !!****f* m_paw_exactDC/compute_exactDC
41 : !! NAME
42 : !! compute_exactDC
43 : !!
44 : !! FUNCTION
45 : !!
46 : !! Compute the exact formula for the double counting.
47 : !! See Physical review letters, Haule, K. (2015), 115(19), 196403 for formula.
48 : !!
49 : !! INPUTS
50 : !! lpawu = angular momentum for correlated species
51 : !! pawtab <type(pawtab_type)>=paw tabulated starting data
52 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data
53 : !! occ(2*lpawu+1,2*lpawu+1) = occupation matrix (summed over spins) in
54 : !! the real spherical harmonics basis, with the convention
55 : !! occ(i,j) = <c_j^dagger c_i>
56 : !! ixc = index of the XC functional
57 : !!
58 : !! OUTPUT
59 : !! vdc(2*lpawu+1,2*lpawu+1) = double counting potential
60 : !! edc = double counting energy
61 : !! edcdc = integral of Vdc(r)*rho_loc(r)
62 : !!
63 : !! SOURCE
64 :
65 0 : subroutine compute_exactDC(lpawu,pawtab,pawrad,occ,vdc,edc,edcdc,ixc)
66 :
67 : use m_pawtab, only : pawtab_type
68 : use m_pawrad, only : pawrad_type,simp_gen
69 : use m_paw_sphharm, only : ylmc,ylmcd
70 : use m_splines, only : spline2
71 :
72 : !Arguments ------------------------------------
73 : integer, intent(in) :: ixc,lpawu
74 : type(pawtab_type), intent(in) :: pawtab
75 : type(pawrad_type), intent(in) :: pawrad
76 : complex(dp), intent(in) :: occ(2*lpawu+1,2*lpawu+1)
77 : complex(dp), intent(inout) :: vdc(2*lpawu+1,2*lpawu+1)
78 : real(dp), intent(out) :: edc,edcdc
79 : !Local variables ------------------------------
80 : integer :: i,ir,j,k,l,m,m1,mm,minusm,meshsz,ndim
81 : logical :: need_gradient
82 : real(dp) :: cphi,ctheta,ecc,excint,exx,grad,gradphi,gradr,gradth
83 : real(dp) :: onemsqrt2,phi,rad,rho,rhor,sexc,sphi,stheta,svxc
84 : real(dp) :: theta,vcc,vxcf,vxx
85 : real(dp) :: kcart(3)
86 : complex(dp) :: dphi,dth,ylm_c
87 : integer, parameter :: ln = 13
88 0 : real(dp), allocatable :: exc(:,:,:),exci(:),phis(:),proj_dr(:),rho_angle(:,:)
89 0 : real(dp), allocatable :: rho_angle_dphi(:,:),rho_angle_dtheta(:,:)
90 0 : real(dp), allocatable :: slm_dphi(:,:,:),slm_dtheta(:,:,:),thetas(:),tweights(:)
91 0 : real(dp), allocatable :: vxc(:,:,:),vxci(:,:,:),ylm(:,:,:)
92 : !************************************************************************
93 :
94 : ! VERY IMPORTANT: This routine assumes that we work in the real spherical harmonic basis.
95 : ! If you want to generalize that to the complex case for some reason, the formulas
96 : ! below are not valid, so you would need to add some complex conjugates.
97 : ! Also, we work with the convention occ(i,j) = <c_j^dagger c_i>.
98 :
99 0 : if (ixc /= 7 .and. ixc /= 11 .and. ixc /= -1012 .and. ixc /= -101130) &
100 0 : & ABI_ERROR("Only PW92 and PBE are handled!")
101 :
102 0 : need_gradient = (ixc == 11 .or. ixc == -101130)
103 :
104 0 : meshsz = size(pawtab%proj2(:)) ! Size of radial mesh
105 0 : ndim = 2*lpawu + 1
106 :
107 0 : ABI_MALLOC(exc,(ln+1,2*ln+1,meshsz))
108 0 : ABI_MALLOC(exci,(meshsz))
109 0 : ABI_MALLOC(phis,(2*ln+1))
110 0 : ABI_MALLOC(rho_angle,(ln+1,2*ln+1))
111 0 : ABI_MALLOC(thetas,(ln+1))
112 0 : ABI_MALLOC(tweights,(ln+1))
113 0 : ABI_MALLOC(vxc,(ln+1,2*ln+1,meshsz))
114 0 : ABI_MALLOC(vxci,(meshsz,ndim,ndim))
115 0 : ABI_MALLOC(ylm,(ln+1,2*ln+1,ndim))
116 :
117 0 : if (need_gradient) then
118 0 : ABI_MALLOC(proj_dr,(meshsz))
119 0 : ABI_MALLOC(rho_angle_dphi,(ln+1,2*ln+1))
120 0 : ABI_MALLOC(rho_angle_dtheta,(ln+1,2*ln+1))
121 0 : ABI_MALLOC(slm_dphi,(ln+1,2*ln+1,ndim))
122 0 : ABI_MALLOC(slm_dtheta,(ln+1,2*ln+1,ndim))
123 : call spline2(pawrad%rad(1:meshsz),pawtab%proj2(:)/max(epsilon(one),pawrad%rad(1:meshsz)**2), &
124 0 : & meshsz,proj_dr(:),zero,zero,3,3)
125 0 : rho_angle_dphi(:,:) = zero
126 0 : rho_angle_dtheta(:,:) = zero
127 : end if ! gradient
128 :
129 0 : vdc(:,:) = czero
130 0 : edc = zero
131 :
132 : ! Hartree contribution
133 :
134 0 : do l=1,ndim
135 0 : do j=1,ndim
136 0 : do k=1,ndim
137 0 : do i=1,ndim
138 : ! In the complex case, you should use conjg(pawtab%vee(i,k,j,l))
139 0 : vdc(i,j) = vdc(i,j) + pawtab%vee(i,k,j,l)*occ(k,l)
140 0 : edc = edc + pawtab%vee(i,k,j,l)*dble(occ(i,j)*occ(k,l))
141 : end do ! i
142 : end do ! k
143 : end do ! j
144 : end do ! l
145 :
146 0 : edc = half * edc
147 :
148 : ! XC contribution
149 :
150 : ! Prepare angular mesh
151 0 : call angular_mesh(thetas(:),phis(:),tweights(:),ln)
152 :
153 : ! Compute the real spherical harmonics on the angular mesh
154 0 : do m=0,lpawu
155 0 : mm = m + lpawu + 1 ; minusm = - m + lpawu + 1
156 0 : onemsqrt2 = (-one)**m * sqrt2
157 0 : do j=1,2*ln+1
158 0 : phi = phis(j) ; cphi = cos(phi) ; sphi = sin(phi)
159 0 : do i=1,ln+1
160 0 : theta = thetas(i) ; ctheta = cos(theta) ; stheta = sin(theta)
161 0 : kcart(1) = cphi * stheta ; kcart(2) = sphi * stheta ; kcart(3) = ctheta
162 0 : ylm_c = ylmc(lpawu,m,kcart(:))
163 : ! Convert complex harmonics to real harmonics
164 0 : if (m == 0) then
165 0 : ylm(i,j,mm) = dble(ylm_c)
166 : else
167 0 : ylm(i,j,mm) = onemsqrt2 * dble(ylm_c)
168 0 : ylm(i,j,minusm) = onemsqrt2 * aimag(ylm_c)
169 : end if ! m=0
170 0 : if (need_gradient) then
171 0 : call ylmcd(lpawu,m,kcart(:),dth,dphi)
172 0 : if (m == 0) then
173 0 : slm_dphi(i,j,mm) = dble(dphi)
174 0 : slm_dtheta(i,j,mm) = dble(dth)
175 : else
176 0 : slm_dphi(i,j,mm) = onemsqrt2 * dble(dphi)
177 0 : slm_dtheta(i,j,mm) = onemsqrt2 * dble(dth)
178 0 : slm_dphi(i,j,minusm) = onemsqrt2 * aimag(dphi)
179 0 : slm_dtheta(i,j,minusm) = onemsqrt2 * aimag(dth)
180 : end if ! m=0
181 : end if ! gradient
182 : end do ! i
183 : end do ! j
184 : end do ! m
185 :
186 0 : rho_angle(:,:) = zero
187 :
188 0 : do m1=1,ndim
189 0 : do m=1,ndim
190 : ! In the complex case, you should use ylm(:,:,m)*conjg(ylm(:,:,m1))
191 0 : rho_angle(:,:) = rho_angle(:,:) + dble(occ(m,m1))*ylm(:,:,m)*ylm(:,:,m1)
192 0 : if (need_gradient) then
193 : ! In the complex case, there should be a conjg each time there is a m1
194 0 : rho_angle_dphi(:,:) = rho_angle_dphi(:,:) + dble(occ(m,m1))*(slm_dphi(:,:,m)*ylm(:,:,m1)+ylm(:,:,m)*slm_dphi(:,:,m1))
195 0 : rho_angle_dtheta(:,:) = rho_angle_dtheta(:,:) + dble(occ(m,m1))*(slm_dtheta(:,:,m)*ylm(:,:,m1)+ylm(:,:,m)*slm_dtheta(:,:,m1))
196 : end if ! gradient
197 : end do ! m
198 : end do ! m1
199 :
200 0 : if (maxval(rho_angle(:,:)-abs(rho_angle(:,:))) > tol10) ABI_WARNING("WARNING: the density is negative !")
201 0 : rho_angle(:,:) = abs(rho_angle(:,:))
202 :
203 0 : do ir=1,meshsz
204 0 : rad = pawrad%rad(ir)
205 0 : rhor = pawtab%proj2(ir) / max(rad**2,epsilon(one))
206 0 : do j=1,2*ln+1
207 0 : do i=1,ln+1
208 :
209 0 : rho = rho_angle(i,j) * rhor
210 0 : grad = zero
211 :
212 0 : if (need_gradient) then
213 0 : stheta = sin(thetas(i))
214 0 : gradr = rho_angle(i,j) * proj_dr(ir)
215 0 : gradth = (rho_angle_dtheta(i,j)*rhor) / max(rad,epsilon(one))
216 0 : gradphi = (rho_angle_dphi(i,j)*rhor) / (max(rad,epsilon(one))*stheta)
217 0 : grad = gradr*gradr + gradth*gradth + gradphi*gradphi
218 : end if ! gradient
219 :
220 0 : call exchange_yukawa(exx,vxx,rho,pawtab%lambda,pawtab%eps,grad,ixc)
221 0 : exc(i,j,ir) = exx ; vxc(i,j,ir) = vxx
222 :
223 0 : call correlation_yukawa(ecc,vcc,rho,pawtab%lambda,pawtab%eps,grad,ixc)
224 0 : exc(i,j,ir) = exc(i,j,ir) + ecc ; vxc(i,j,ir) = vxc(i,j,ir) + vcc
225 :
226 : end do ! i
227 : end do ! j
228 : end do ! ir
229 :
230 0 : exci(:) = zero
231 :
232 : ! Integrals over theta and phi
233 :
234 0 : do m1=1,ndim
235 0 : do m=1,ndim
236 0 : do ir=1,meshsz
237 : svxc = zero
238 : sexc = zero
239 0 : do j=1,2*ln+1
240 : ! In the complex case, you should use ylm(:,j,m)*conjg(ylm(:,j,m1))
241 0 : svxc = svxc + sum(ylm(:,j,m)*ylm(:,j,m1)*tweights(:)*vxc(:,j,ir))
242 0 : sexc = sexc + sum(ylm(:,j,m)*ylm(:,j,m1)*tweights(:)*exc(:,j,ir))
243 : end do ! j
244 0 : vxci(ir,m,m1) = two_pi * svxc / dble(2*ln+1)
245 0 : exci(ir) = exci(ir) + two_pi * sexc * dble(occ(m,m1)) / dble(2*ln+1)
246 : end do ! ir
247 : end do ! m
248 : end do ! m1
249 :
250 : ! Integrals over r
251 :
252 0 : call simp_gen(excint,exci(:)*pawtab%proj2(:),pawrad,r_for_intg=pawrad%rad(meshsz))
253 0 : edc = edc + excint
254 0 : edcdc = zero
255 :
256 : ! Careful, vdc(i,j) as defined here is the i,j-th matrix element of the TRANSPOSE of vdc !
257 :
258 0 : do m1=1,ndim
259 0 : do m=1,ndim
260 0 : call simp_gen(vxcf,vxci(:,m,m1)*pawtab%proj2(:),pawrad,r_for_intg=pawrad%rad(meshsz))
261 0 : vdc(m,m1) = vdc(m,m1) + cmplx(vxcf,zero,kind=dp)
262 0 : edcdc = edcdc + dble(vdc(m,m1)*occ(m,m1))
263 : end do ! m
264 : end do ! m1
265 :
266 0 : ABI_FREE(exc)
267 0 : ABI_FREE(exci)
268 0 : ABI_FREE(phis)
269 0 : ABI_FREE(rho_angle)
270 0 : ABI_FREE(thetas)
271 0 : ABI_FREE(tweights)
272 0 : ABI_FREE(vxc)
273 0 : ABI_FREE(vxci)
274 0 : ABI_FREE(ylm)
275 0 : ABI_SFREE(proj_dr)
276 0 : ABI_SFREE(rho_angle_dphi)
277 0 : ABI_SFREE(rho_angle_dtheta)
278 0 : ABI_SFREE(slm_dphi)
279 0 : ABI_SFREE(slm_dtheta)
280 :
281 0 : end subroutine compute_exactDC
282 : !!***
283 :
284 : !----------------------------------------------------------------------
285 :
286 : !!****f* m_paw_exactDC/angular_mesh
287 : !! NAME
288 : !! angular_mesh
289 : !!
290 : !! FUNCTION
291 : !!
292 : !! Prepare the angular mesh for the integration over thetas and phis.
293 : !!
294 : !! INPUTS
295 : !! ln = controls the number of integration points
296 : !!
297 : !! OUTPUT
298 : !! thetas(ln+1) = Gauss-Legendre integration points for theta grid
299 : !! phis(2*ln+1) = uniform points for phi grid
300 : !! tweights(ln+1) = Gauss-Legendre weights for theta grid
301 : !!
302 : !! SOURCE
303 :
304 0 : subroutine angular_mesh(thetas,phis,tweights,ln)
305 :
306 : use m_numeric_tools, only : coeffs_gausslegint
307 :
308 : !Arguments ------------------------------------
309 : integer, intent(in) :: ln
310 : real(dp), intent(inout) :: phis(2*ln+1),thetas(ln+1),tweights(ln+1)
311 : !Local variables ------------------------------
312 : integer :: i,lg
313 : real(dp) :: dsum,phi
314 : real(dp), parameter :: fake_shift = exp(-4.0_dp) ! small shift such that we do not start at phi=0
315 : !************************************************************************
316 :
317 0 : lg = ln + 1
318 0 : call coeffs_gausslegint(-1.0_dp,1.0_dp,thetas(:),tweights(:),lg)
319 0 : do i=1,lg
320 0 : thetas(i) = acos(thetas(i))
321 : end do
322 :
323 0 : do i=0,2*lg-2 ! 2*lg-1 points in phi direction
324 0 : phi = pi * (two*dble(i)/dble(2*lg-1)+fake_shift)
325 0 : phis(i+1) = phi
326 : end do ! i
327 :
328 0 : dsum = sum(tweights(:))
329 0 : tweights(:) = tweights(:) * 2.0_dp / dsum
330 :
331 0 : end subroutine angular_mesh
332 : !!***
333 :
334 : !----------------------------------------------------------------------
335 :
336 : !!****f* m_paw_exactDC/exchange_yukawa
337 : !! NAME
338 : !! exchange_yukawa
339 : !!
340 : !! FUNCTION
341 : !!
342 : !! Compute the exchange contribution with a Yukawa potential
343 : !!
344 : !! INPUTS
345 : !! rho = density at current point
346 : !! lambda = parameter for Yukawa potential (inverse screening length)
347 : !! eps = parameter for Yukawa potential (dielectric constant)
348 : !! grad = square of the gradient of the density
349 : !! ixc = index of the XC functional
350 : !!
351 : !! OUTPUT
352 : !! ex = exchange energy per particle
353 : !! vx = exchange potential
354 : !!
355 : !! SOURCE
356 :
357 0 : subroutine exchange_yukawa(ex,vx,rho,lambda,eps,grad,ixc)
358 :
359 : !Arguments ------------------------------------
360 : integer, intent(in) :: ixc
361 : real(dp), intent(in) :: grad,lambda,eps,rho
362 : real(dp), intent(out) :: ex,vx
363 : !Local variables ------------------------------
364 : logical :: islambda
365 : real(dp) :: dfdkappa,dfdmu,dfdrho,dfdss,dfx,div,div2,dkappadx,dlogf,dmudx,dssdrho
366 : real(dp) :: dxdrho,fx,fx_pbe,kappa,mu,rhothird,rhotwothird,rsinv,ss,x,x2
367 : real(dp), parameter :: c0 = (9.0_dp/(4.0_dp*(pi**2)))**third * three_quarters
368 : real(dp), parameter :: c1 = 1.804_dp
369 : real(dp), parameter :: kf_fac = (3.0_dp*(pi**2))**third
370 : real(dp), parameter :: mu0 = 0.2195149727645171_dp
371 : real(dp), parameter :: rsinv_fac = (4.0_dp*pi/3.0_dp)**third
372 : real(dp), parameter :: twotwothird = (2.0_dp)**(2.0_dp*third)
373 : real(dp), parameter :: x_fac = (9.0_dp*pi/4.0_dp)**third
374 : !************************************************************************
375 :
376 0 : rhothird = rho**third ; rsinv = rsinv_fac * rhothird
377 0 : islambda = (abs(lambda) > tol10)
378 :
379 : ! LDA exchange
380 0 : if (islambda) then
381 0 : x = x_fac * rsinv / lambda
382 0 : call fexchange(x,fx,dfx)
383 : else
384 0 : fx = 1.0_dp
385 0 : dfx = 0.0_dp
386 : end if
387 :
388 0 : ex = - c0 * fx * rsinv / eps
389 0 : vx = 4.0_dp*ex/3.0_dp - c0*x*dfx*rsinv/(3.0_dp*eps)
390 :
391 0 : if (ixc == 7 .or. ixc == -1012) return
392 :
393 : ! PBE exchange
394 0 : if (abs(rho) < tol30) then
395 0 : ex = zero ; vx = zero
396 0 : return
397 : end if
398 :
399 0 : if (islambda) then
400 0 : x2 = x * x
401 0 : kappa = x2 / (x2 + twotwothird)
402 : else
403 : kappa = 1.0_dp
404 : end if
405 :
406 0 : kappa = c1*kappa/fx - 1.0_dp
407 0 : rhotwothird = rhothird * rhothird
408 0 : ss = grad / (4.0_dp*(rho**2)*rhotwothird*(kf_fac**2))
409 0 : mu = mu0 / fx
410 0 : div = 1.0_dp+mu*ss/kappa
411 0 : div2 = div * div
412 0 : fx_pbe = 1.0_dp + kappa*(1.0_dp-1.0_dp/div)
413 0 : dssdrho = -8.0_dp * third * ss / rho
414 0 : dfdss = mu / div2
415 0 : dfdrho = dfdss * dssdrho
416 0 : if (islambda) then
417 0 : dlogf = dfx / fx
418 0 : dxdrho = third * x / rho
419 0 : dkappadx = (kappa+1.0_dp) * (2.0_dp*twotwothird/(x*(x2+twotwothird))-dlogf)
420 0 : dfdkappa = 1.0_dp - (1.0_dp+2.0_dp*mu*ss/kappa)/div2
421 0 : dmudx = -mu * dlogf
422 0 : dfdmu = ss / div2
423 0 : dfdrho = dfdrho + (dfdkappa*dkappadx+dfdmu*dmudx)*dxdrho
424 : end if
425 0 : vx = vx*fx_pbe + rho*ex*dfdrho
426 0 : ex = ex * fx_pbe
427 :
428 : end subroutine exchange_yukawa
429 : !!***
430 :
431 : !----------------------------------------------------------------------
432 :
433 : !!****f* m_paw_exactDC/fexchange
434 : !! NAME
435 : !! fexchange
436 : !!
437 : !! FUNCTION
438 : !!
439 : !! Compute the function in the HEG exchange energy with Yukawa potential
440 : !!
441 : !! INPUTS
442 : !! x = input of the function (=(9*pi/4)**(1/3) / (lambda*rs)
443 : !! with lambda the Yukawa parameter and rs the Wigner-Seitz radius)
444 : !!
445 : !! OUTPUT
446 : !! fx = value of the function at x
447 : !! dfx = derivative of the function at x
448 : !!
449 : !! SOURCE
450 :
451 0 : subroutine fexchange(x,fx,dfx)
452 :
453 : !Arguments ------------------------------------
454 : real(dp), intent(in) :: x
455 : real(dp), intent(out) :: fx,dfx
456 : !Local variables ------------------------------
457 : real(dp) :: at2,lg2,x2,x3,x4,x5
458 : !************************************************************************
459 :
460 0 : x2 = x * x ; x3 = x2 * x ; x4 = x3 * x ; x5 = x4 * x
461 0 : if (x < tol2) then ! Taylor expansion
462 0 : fx = 4.0_dp * x2 * (1.0_dp/9.0_dp-2.0_dp*x2/15.0_dp+8.0_dp*x4/35.0_dp)
463 0 : dfx = 8.0_dp * x * (1.0_dp/9.0_dp-4.0_dp*x2/15.0_dp+24.0_dp*x4/35.0_dp)
464 : else
465 0 : at2 = atan(2.0_dp*x) ; lg2 = log(1.0_dp+4.0_dp*x2)
466 : fx = 1.0_dp - 1.0_dp/(6.0_dp*x2) - 4.0_dp*at2/(3.0_dp*x) + &
467 0 : & (1.0_dp+12.0_dp*x2)*lg2/(24.0_dp*x4)
468 0 : dfx = (8.0_dp*x3*at2+4.0_dp*x2-(1.0_dp+6.0_dp*x2)*lg2)/(6.0_dp*x5)
469 : end if
470 :
471 0 : end subroutine fexchange
472 : !!***
473 :
474 : !!****f* m_paw_exactDC/correlation_yukawa
475 : !! NAME
476 : !! correlation_yukawa
477 : !!
478 : !! FUNCTION
479 : !!
480 : !! Compute the correlation contribution with a Yukawa potential
481 : !!
482 : !! INPUTS
483 : !! rho = density at current point
484 : !! lambda = parameter for screened potential (inverse screening length)
485 : !! eps = parameter for screened potential (dielectric constant)
486 : !! grad = square of the gradient of the density
487 : !! ixc = index of the XC functional
488 : !!
489 : !! OUTPUT
490 : !! ec = correlation energy per particle
491 : !! vc = correlation potential
492 : !!
493 : !! SOURCE
494 :
495 0 : subroutine correlation_yukawa(ec,vc,rho,lambda,eps,grad,ixc)
496 :
497 : !Arguments ------------------------------------
498 : integer, intent(in) :: ixc
499 : real(dp), intent(in) :: eps,grad,lambda,rho
500 : real(dp), intent(out) :: ec,vc
501 : !Local variables ------------------------------
502 : logical :: islambda
503 : real(dp) :: aa_pbe,alb,alb_pow,arg_log,daadec,decdrho,den,df,dhdaa
504 : real(dp) :: dhdrho,dhdtt,div,div2,dttdrho,eps2,eps3,exp_pbe,fx,h_pbe,lg,pade,pow,q0
505 : real(dp) :: q1,q1p,rhothird,rs,sqr_rs,tt,xx
506 : real(dp), parameter :: aa = 0.031091_dp,a1 = 0.21370_dp
507 : real(dp), parameter :: b1 = 7.5957_dp,b2 = 3.5876_dp
508 : real(dp), parameter :: b3 = 1.6382_dp,b4 = 0.49294_dp
509 : real(dp), parameter :: a = 0.47808102_dp,b = 0.84449703_dp
510 : real(dp), parameter :: c = 1.30089155_dp,d = 0.02949437_dp,beta = 1.34835105_dp
511 : real(dp), parameter :: kf_fac = (3.0_dp*(pi**2))**third
512 : real(dp), parameter :: rs_fac = (3.0_dp/(4.0_dp*pi))**third
513 : real(dp), parameter :: beta_pbe = 0.066725_dp,gamma_pbe = (1.0_dp-log(2.0_dp))/(pi**2)
514 : !************************************************************************
515 :
516 0 : if (abs(rho) < tol30) then
517 0 : ec = 0.0_dp
518 0 : vc = 0.0_dp
519 0 : return
520 : end if
521 :
522 0 : rhothird = rho**third
523 0 : rs = rs_fac / (rhothird*eps) ! Scaling law for r_s
524 0 : sqr_rs = sqrt(rs)
525 :
526 0 : eps2 = eps * eps
527 :
528 : ! LDA correlation for rescaled Coulomb potential (PW91 parametrization)
529 0 : q0 = -2.0_dp * aa * (1.0_dp+a1*rs)
530 0 : q1 = 2.0_dp * aa * (b1*sqr_rs+b2*rs+b3*rs*sqr_rs+b4*rs*rs)
531 0 : q1p = aa * (b1/sqr_rs+2._dp*b2+3._dp*b3*sqr_rs+4._dp*b4*rs)
532 0 : den = 1.0_dp / (q1*q1+q1)
533 0 : lg = -log(q1*q1*den)
534 0 : ec = q0 * lg
535 0 : vc = -2.0_dp*aa*a1*lg - q0*q1p*den
536 0 : vc = ec - rs*vc/3.0_dp
537 :
538 0 : islambda = (abs(lambda) > tol10)
539 :
540 : ! Correction factor for Yukawa potential
541 : ! Parametrization based on data from Savin, "Beyond the Kohn-Sham Determinant"
542 0 : pow = c + d*log(1.0_dp+rs)
543 0 : alb = a*lambda*eps*(rs**b) ; alb_pow = alb**pow ! Scaling law for lambda
544 0 : fx = (1.0_dp+alb_pow)**(-beta-1.0_dp)
545 0 : if (islambda) then
546 0 : df = -beta * fx * alb_pow * (pow*b/rs+d*log(alb)/(1.0_dp+rs))
547 : else
548 : df = zero
549 : end if
550 0 : fx = fx * (1.0_dp+alb_pow)
551 :
552 0 : vc = (vc*fx-ec*rs*df/3.0_dp) / eps2 ! Scaling law
553 0 : ec = ec * fx / eps2
554 :
555 0 : if (ixc == 7 .or. ixc == -1012) return
556 :
557 0 : eps3 = eps2 * eps
558 0 : tt = grad * pi / ((rho**2)*16.0_dp*kf_fac*rhothird)
559 0 : exp_pbe = exp(-ec*eps2/gamma_pbe)
560 0 : if (abs(exp_pbe-1.0_dp) < tol30) return
561 :
562 0 : aa_pbe = eps * beta_pbe / (gamma_pbe*(exp_pbe-1.0_dp))
563 0 : daadec = eps3 * beta_pbe * exp_pbe / (gamma_pbe*(exp_pbe-1.0_dp))**2
564 0 : decdrho = (vc-ec) / rho
565 0 : xx = aa_pbe * tt
566 0 : div = 1.0_dp + xx + xx**2
567 0 : div2 = div * div
568 0 : pade = (1.0_dp+xx) / div
569 0 : arg_log = 1.0_dp + eps*beta_pbe*tt*pade/gamma_pbe
570 0 : h_pbe = gamma_pbe * log(arg_log) / eps2
571 0 : dhdaa = -beta_pbe * (tt**2) * xx * (2.0_dp+xx) / (eps*div2*arg_log)
572 0 : dttdrho = -7.0_dp * tt / (3.0_dp*rho)
573 0 : dhdtt = beta_pbe * (pade-(xx**2)*(2.0_dp+xx)/div2) / (arg_log*eps)
574 0 : dhdrho = dhdtt*dttdrho + dhdaa*daadec*decdrho
575 0 : vc = vc + h_pbe + rho*dhdrho
576 0 : ec = ec + h_pbe
577 :
578 : end subroutine correlation_yukawa
579 : !!***
580 :
581 : END MODULE m_paw_exactDC
582 : !!***
583 :
|
