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