1 2! Copyright (C) 2004, 2005 Rickard Armiento 3! This file is distributed under the terms of the GNU Lesser General Public 4! License. See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: xc_am05 8! !INTERFACE: 9subroutine xc_am05(n,rho,grho,g2rho,g3rho,ex,ec,vx,vc) 10! !INPUT/OUTPUT PARAMETERS: 11! n : number of density points (in,integer) 12! rho : charge density (in,real(n)) 13! grho : |grad rho| (in,real(n)) 14! g2rho : grad^2 rho (in,real(n)) 15! g3rho : (grad rho).(grad |grad rho|) (in,real(n)) 16! ex : exchange energy density (out,real(n)) 17! ec : correlation energy density (out,real(n)) 18! vx : spin-unpolarised exchange potential (out,real(n)) 19! vc : spin-unpolarised correlation potential (out,real(n)) 20! !DESCRIPTION: 21! Spin-unpolarised exchange-correlation potential and energy functional of 22! R. Armiento and A. E. Mattsson, {\it Phys. Rev. B} {\bf 72}, 085108 (2005). 23! 24! !REVISION HISTORY: 25! Created April 2005 (RAR); based on xc_pbe 26!EOP 27!BOC 28implicit none 29! arguments 30integer, intent(in) :: n 31real(8), intent(in) :: rho(n),grho(n),g2rho(n),g3rho(n) 32real(8), intent(out) :: ex(n),ec(n),vx(n),vc(n) 33! local variables 34integer i 35real(8), parameter :: pi=3.1415926535897932385d0 36real(8) r,kf,s,v,u 37real(8) grho_,g2rho_,g3rho_ 38do i=1,n 39 r=rho(i) 40 if (r.gt.1.d-12) then 41 grho_=grho(i) 42 g2rho_=g2rho(i) 43 g3rho_=g3rho(i) 44! exchange energy density and potential 45 kf=(r*3.d0*pi**2)**(1.d0/3.d0) 46 s=grho_/(2.d0*kf*r) 47 v=g2rho_/(r*(2.d0*kf)**2) 48 u=g3rho_/((r**2)*(2.d0*kf)**3) 49 call xc_am05_point(r,s,u,v,ex(i),ec(i),vx(i),vc(i),1) 50 else 51 ex(i)=0.d0 52 ec(i)=0.d0 53 vx(i)=0.d0 54 vc(i)=0.d0 55 end if 56end do 57end subroutine 58!EOC 59 60!BOP 61! !ROUTINE: xc_am05_point 62! !INTERFACE: 63subroutine xc_am05_point(rho,s,u,v,ex,ec,vx,vc,pot) 64! !INPUT/OUTPUT PARAMETERS: 65! rho : electron density (in,real) 66! s : gradient of n / (2 kF n) 67! u : grad n * grad | grad n | / (n**2 (2 kF)**3) 68! v : laplacian of density / (n**2 (2.d0*kf)**3) 69! ex : exchange energy density (out,real) 70! ec : correlation energy density (out,real) 71! vx : spin-unpolarised exchange potential (out,real) 72! vc : spin-unpolarised correlation potential (out,real) 73! !DESCRIPTION: 74! Calculate the spin-unpolarised exchange-correlation potential and energy for 75! the Armiento-Mattsson 05 functional for a single point. 76! 77! !REVISION HISTORY: 78! Created April 2005 (RAR) 79!EOP 80!BOC 81implicit none 82! arguments 83real(8), intent(in) :: rho, s, u, v 84integer, intent(in) :: pot 85real(8), intent(out) :: ex, ec, vx, vc 86! constants 87real(8), parameter :: pi=3.1415926535897932385d0 88real(8), parameter :: g = 0.8098d0 89real(8), parameter :: a = 2.804d0 90real(8), parameter :: c = 0.7168d0 91! local variables 92real(8) s2,exlda, vxlda, eclda, vclda, X, Xs, Xss 93real(8) F, Fs, Fss, Hx, Hxs, Hxss, Hc, Hcs, Hcss 94real(8) zb, zbs, zbss, w 95real(8) n0b, n0bs, n0bss 96real(8) ln0b, ln0bs, ln0bss 97real(8) zbb, zbbc, zbbs, zbbss 98real(8) fxb, fxbs, fxbss 99! cutoff 100if((rho .le. 1.0d-16)) then 101 ex = 0.0d0 102 ec = 0.0d0 103 vx = 0.0d0 104 vc = 0.0d0 105 return 106endif 107s2 = s**2 108! LDA correlation 109call xc_am05_ldapwc(rho,eclda,vclda) 110! LDA exchange 111call xc_am05_ldax(rho,exlda,vxlda) 112!------------------! 113! exchange ! 114!------------------! 115! interpolation index 116X = 1.0d0 - a*s2/(1.0d0 + a*s2) 117! Airy LAA refinement function 118call xc_am05_labertw(s**(3.0d0/2.0d0)/sqrt(24.0d0),w) 119zb = (3.0d0/2.0d0*w)**(2.0d0/3.0d0) 120n0b = w/(2.0d0*pi**2*s**3) 121ln0b = -3.0d0/(2.0d0*pi)*(3.0d0*pi**2*n0b)**(1.0d0/3.0d0) 122zbbc = ((4.0d0/3.0d0)**(1.0d0/3.0d0)*2.0d0*pi/3.0d0)**4 123zbb = (zbbc*zb**2 + zb**4)**(1.0d0/4.0d0) 124Fxb = -1.0d0/(ln0b*2.0d0*zbb) 125F = (c*s2 + 1.0d0)/(c*s2/Fxb + 1.0d0) 126! exchange refinement function 127Hx = X + (1.0d0 - X)*F 128! exchange energy per particle, Ex = Integrate[n*ex] 129ex = exlda*Hx 130!---------------------! 131! correlation ! 132!---------------------! 133! correlation refinement function 134Hc = X + g*(1.0d0 - X) 135! correlation energy per particle, Ec = Integrate[rho*ec] 136ec = eclda*Hc 137if (pot .eq. 0) return 138!----------------------------! 139! exchange potential ! 140!----------------------------! 141! interpolation index derivatives, dX/ds 142Xs = -2.0d0*a*s/(1.0d0 + a*s2)**2 143Xss = 2.0d0*a*(3.0d0*a*s2-1.0d0)/(1.0d0+a*s2)**3 144! airy LAA refinement function derivatives, dF/ds 145zbs = zb/(s + s*w) 146zbss = - zb*w*(5.0d0+2.0d0*w)/(2.0d0*s2*(1.0d0+w)**3) 147n0bs = sqrt(zb)*(-2.0d0*zb+s*zbs)/(2.0d0*pi**2*s2**2) 148n0bss = (16.0d0*zb**2+s**2*zbs**2+2.0d0*s*zb*(-6.0d0* & 149 zbs+s*zbss))/(4.0d0*pi**2*s**5*sqrt(zb)) 150ln0bs = -(3.0d0/pi)**(1.0d0/3.0d0)*n0bs/ & 151 (2.0d0*n0b**(2.0d0/3.0d0)) 152ln0bss = (2.0d0*n0bs**2-3.0d0*n0b*n0bss)/(2.0d0* & 153 3.0d0**(2.0d0/3.0d0)*pi**(1.0d0/3.0d0)*n0b**(5.0d0/3.0d0)) 154zbbs = zb*(zbbc+2*zb**2)*zbs/ & 155 (2.0d0*(zb**2*(zbbc+zb**2))**(3.0d0/4.0d0)) 156zbbss = zb**2*(-zbbc*(zbbc-2.0d0*zb**2)*zbs**2+ & 157 2.0d0*zb*(zbbc+zb**2)*(zbbc+2.0d0*zb**2)*zbss)/ & 158 (4.0d0*(zb**2*(zbbc+zb**2))**(7.0d0/4.0d0)) 159Fxbs = (zbb*ln0bs+ln0b*zbbs)/(2.0d0*ln0b**2*zbb**2) 160Fxbss = (-2.0d0*ln0b**2*zbbs**2+zbb**2*(-2.0d0*ln0bs**2 + & 161 ln0b*ln0bss)+ln0b*zbb*(-2.0d0*ln0bs*zbbs+ln0b*zbbss))/ & 162 (2.0d0*ln0b**3*zbb**3) 163Fs = (c*s*(2.0d0*(Fxb-1.0d0)*Fxb + s*(1.0d0+c*s2)*Fxbs))/ & 164 (c*s2 + Fxb)**2 165Fss = (c*(-2.0d0*(3.0d0*c*s2-Fxb)*(Fxb-1.0d0)*Fxb+ & 166 4.0d0*s*(-c*s2+Fxb+2.0d0*c*s2*Fxb)*Fxbs - & 167 2.0d0*s2*(1.0d0+c*s2)*Fxbs**2+s2*(1.0d0+c*s2)* & 168 (c*s2 + Fxb)*Fxbss))/(c*s2+Fxb)**3 169! GGA refinement function derivatives, dF/ds 170Hxs = - (X - 1.0d0)*Fs - (F - 1.0d0)*Xs 171Hxss = - 2.0d0*Fs*Xs - (X - 1.0d0)*Fss - (F - 1.0d0)*Xss 172! vx formula for gradient dependent functional, 173! generalized form of Eq. (24) in PRB 33, 8800 (1986) 174vx = vxlda*(Hx - s*Hxs) + & 175 exlda*((4.0d0/3.0d0*s-v/s)*Hxs - & 176 (u-4.0d0/3.0d0*s**3)*(Hxss/s-Hxs/s2)) 177!-------------------------------! 178! correlation potential ! 179!-------------------------------! 180! correlation refinement function derivatives, dF/ds 181Hcs = Xs - g*Xs 182Hcss = Xss - g*Xss 183! vc formula for gradient dependent functional, 184! generalized form of Eq. (24) in Phys. Rev. B 33, 8800 (1986) 185vc = vclda*(Hc - s*Hcs) + & 186 eclda*((4.0d0/3.0d0*s - v/s)*Hcs - & 187 (u - 4.0d0/3.0d0*s**3)*(Hcss/s - Hcs/s2)) 188end subroutine 189!EOC 190 191!BOP 192! !ROUTINE: xc_am05_ldax 193! !INTERFACE: 194subroutine xc_am05_ldax(n,ex,vx) 195! !INPUT/OUTPUT PARAMETERS: 196! n : electron density (in,real) 197! ex : exchange energy per electron (out,real) 198! vx : exchange potential (out,real) 199! !DESCRIPTION: 200! Local density approximation exchange. 201! 202! !REVISION HISTORY: 203! Created April 2005 (RAR) 204!EOP 205!BOC 206implicit none 207! arguments 208real(8), intent(in) :: n 209real(8), intent(out) :: ex 210real(8), intent(out) :: vx 211! constants 212real(8), parameter :: pi=3.1415926535897932385d0 213vx=-(3.d0*n/pi)**(1.d0/3.d0) 214ex=(3.d0/4.d0)*vx 215end subroutine 216!EOC 217 218!BOP 219! !ROUTINE: xc_am05_ldapwc 220! !INTERFACE: 221subroutine xc_am05_ldapwc(n,ec,vc) 222! !INPUT/OUTPUT PARAMETERS: 223! n : electron density (in,real) 224! ec : correlation energy per electron (out,real) 225! vc : correlation potential (out,real) 226! !DESCRIPTION: 227! Correlation energy and potential of the Perdew-Wang parameterisation of 228! the Ceperley-Alder electron gas {\it Phys. Rev. B} {\bf 45}, 13244 (1992) 229! and {\it Phys. Rev. Lett.} {\bf 45}, 566 (1980). This is a clean-room 230! implementation from paper. 231! 232! !REVISION HISTORY: 233! Created April 2005 (RAR) 234!EOP 235!BOC 236implicit none 237! arguments 238real(8), intent(in) :: n 239real(8), intent(out) :: ec 240real(8), intent(out) :: vc 241! constants 242real(8), parameter :: pi=3.1415926535897932385d0 243real(8), parameter :: a01 = 0.21370d0 244real(8), parameter :: b01 = 7.5957d0 245real(8), parameter :: b02 = 3.5876d0 246real(8), parameter :: b03 = 1.6382d0 247real(8), parameter :: b04 = 0.49294d0 248! paper actually use this: 249! real(8), parameter (A0 = 0.031091d0) 250! but routines now "defacto standard" was distributed using: 251real(8), parameter :: A0 = 0.0310907d0 252! local variables 253real(8) rsq 254real(8) Q0, Q1, Q1p, ecrs 255rsq = (3.0d0/(4.0d0*pi*n))**(1.0d0/6.0d0) 256ec = -2.0d0*A0*(1.0d0 + a01*rsq**2)* & 257 log(1.0d0 + 1.0d0/ & 258 (2.0d0*A0*rsq*(b01 + rsq*(b02 + rsq*(b03 + b04*rsq))))) 259Q0 = -2.0d0*A0*(1.0d0 + a01*rsq**2) 260Q1 = 2.0d0*A0*rsq*(b01 + rsq*(b02 + rsq*(b03 + b04*rsq))) 261Q1p = A0*(b01/rsq+2.0d0*b02+3.0d0*b03*rsq+4.0d0*b04*rsq**2) 262ecrs = -2.0d0*A0*a01*log(1.0d0 + 1.0d0/Q1)-Q0*Q1p/(Q1**2+Q1) 263vc = ec - rsq**2/3.0d0*ecrs 264end subroutine 265!EOC 266 267!BOP 268! !ROUTINE: xc_am05_labertw 269! !INTERFACE: 270subroutine xc_am05_labertw(z,val) 271! !INPUT/OUTPUT PARAMETERS: 272! z : function argument (in,real) 273! val : value of lambert W function of z (out,real) 274! !DESCRIPTION: 275! Lambert $W$-function using the method of Corless, Gonnet, Hare, Jeffrey and 276! Knuth, {\it Adv. Comp. Math.} {\bf 5}, 329 (1996). The approach is based 277! loosely on that in GNU Octave by N. N. Schraudolph, but this implementation 278! is for real values and the principal branch only. 279! 280! !REVISION HISTORY: 281! Created April 2005 (RAR) 282!EOP 283!BOC 284implicit none 285! arguments 286real(8), intent(in) :: z 287real(8), intent(out) :: val 288! local variables 289real(8) e,t,p 290integer i 291! if z too low, go with the first term of the power expansion, z 292if (z.lt.1.d-20) then 293 val=z 294 return 295end if 296e=exp(1.d0) 297! inital guess 298if (abs(z+1.d0/e).gt.1.45d0) then 299! asymptotic expansion at 0 and Inf 300 val=log(z) 301 val=val-log(val) 302else 303! series expansion about -1/e to first order 304 val=1.d0*sqrt(2.d0*e*z+2.d0)-1.d0 305end if 306! find val through iteration 307do i=1,10 308 p=exp(val) 309 t=val*p-z 310 if (val.ne.-1.d0) then 311 t=t/(p*(val+1.d0)-0.5d0*(val+2.d0)*t/(val+1.d0)) 312 else 313 t=0.d0 314 end if 315 val=val-t 316 if (abs(t).lt.(2.48d0*1.d-14)*(1.d0+abs(val))) return 317end do 318! this should never happen! 319write(*,*) 320write(*,'("Error(xc_am05_labertw): iteration limit reached")') 321write(*,'(" Likely cause: improper numbers (INFs, NaNs) in density")') 322write(*,*) 323stop 324end subroutine 325!EOC 326 327