Branch data 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 : :
|
