1*
2* $Id$
3*
4
5      subroutine nwpw_vosko(n2ft3d,ng,ispin,dn,xce,xcp,x)
6      implicit none
7
8      integer n2ft3d
9      integer ng
10      integer ispin
11      real*8 dn(n2ft3d,2)
12      real*8 xce(n2ft3d,2)
13      real*8 xcp(n2ft3d,2)
14      real*8 x(n2ft3d)
15
16      double precision one3rd,one6th,dncut,for3rd
17      parameter (one3rd=1.0d0/3.0d0,for3rd=4.0d0/3.0d0)
18      parameter (one6th=1.0d0/6.0d0)
19      parameter (dncut=1.0d-30)
20
21*---- parameters given by vosko et al -----------------*
22      double precision ap,x0p,bp,cp,af,x0f,bf,cf
23      parameter (ap = 3.109070d-02, af = 1.554530d-02)
24      parameter (x0p=-1.049800d-01, x0f=-3.250000d-01)
25      parameter (bp = 3.727440d+00, bf = 7.060420d+00)
26      parameter (cp = 1.293520d+01, cf = 1.805780d+01)
27*------------------------------------------------------*
28
29*     constants calculated from vosko's parameters
30      double precision xp,qp,xx0p,cp1,cp2,cp3,cp4,cp5,cp6
31      double precision dp1,dp2,dp3,dp4,dp5,dp6,dp7,fc,crs
32      double precision xf,qf,xx0f,cf1,cf2,cf3,cf4,cf5,cf6
33      double precision df1,df2,df3,df4,df5,df6,df7,fd
34      parameter (xp  = -4.581653d-01,  xf  = -5.772521d-01)
35      parameter (qp  =  6.151991d+00,  qf  =  4.730927d+00)
36      parameter (xx0p=  1.255491d+01,  xx0f=  1.586879d+01)
37      parameter (cp1 =  3.109070d-02,  cf1 =  1.554530d-02)
38      parameter (cp2 =  9.690228d-04,  cf2 =  2.247860d-03)
39      parameter (cp3 =  1.049800d-01,  cf3 =  3.250000d-01)
40      parameter (cp4 =  3.878329d-02,  cf4 =  5.249122d-02)
41      parameter (cp5 =  3.075995d+00,  cf5 =  2.365463d+00)
42      parameter (cp6 =  1.863720d+00,  cf6 =  3.530210d+00)
43      parameter (dp1 =  6.218140d-02,  df1 =  3.109060d-02)
44      parameter (dp2 =  1.938045d-03,  df2 =  4.495720d-03)
45      parameter (dp3 =  1.049800d-01,  df3 =  3.250000d-01)
46      parameter (dp4 = -3.205972d-02,  df4 = -1.779316d-02)
47      parameter (dp5 = -1.192972d-01,  df5 = -1.241661d-01)
48      parameter (dp6 =  1.863720d+00,  df6 =  3.530210d+00)
49      parameter (dp7 =  9.461748d+00,  df7 =  5.595417d+00)
50      parameter (fc  =  1.923661d+00,  fd  =  2.564881d+00)
51      parameter (crs =  7.876233d-01)
52
53      double precision df,pi,rho,zdw,zup,f,xx
54
55      integer k
56
57      pi=4.0d0*datan(1.0d0)
58
59*     square root of wigner radius
60!$OMP DO
61      do 100 k=1,ng
62c        rho=dn(k,1)+dn(k,ispin)+dncut
63        rho=dn(k,1)+dn(k,ispin)+dncut
64        x(k)=crs/rho**one6th
65  100 end do
66!$OMP END DO
67
68*     paramagnetic correlation energy & potential
69!$OMP DO
70      do 110 k=1,ng
71        xx=1.0d0/(x(k)*(x(k)+bp)+cp)
72        xce(k,1)=cp1*dlog(xx*x(k)**2)+cp2*dlog(xx*(x(k)+cp3)**2)
73     &          +cp4*datan(cp5/(x(k)+cp6))
74        xcp(k,1)=xce(k,1)-one6th*x(k)*(
75     &           dp1/x(k)+dp2/(x(k)+dp3)+dp4*xx*(2.0d0*x(k)+bp)
76     &          +dp5/((x(k)+dp6)**2+dp7) )
77  110 end do
78!$OMP END DO
79
80*     paramagnetic exchange energy & potential
81!$OMP DO
82      do 120 k=1,ng
83        xce(k,1)=xce(k,1)+(xp/x(k)**2)
84        xcp(k,1)=xcp(k,1)+for3rd*(xp/x(k)**2)
85  120 end do
86!$OMP END DO
87
88*     return if spin-restricted lda
89*      if(ispin.eq.1) go to 200
90      if(ispin.eq.1) then
91!$OMP DO
92        do k=1,ng
93          xcp(k,2) = xcp(k,1)
94        end do
95!$OMP END DO
96        return
97      end if
98
99*     ferromagnetic correlation energy & potential
100!$OMP DO
101      do 130 k=1,ng
102        xx=1.0d0/(x(k)*(x(k)+bf)+cf)
103        xce(k,2)=cf1*dlog(xx*x(k)**2)+cf2*dlog(xx*(x(k)+cf3)**2)
104     &          +cf4*datan(cf5/(x(k)+cf6))
105        xcp(k,2)=xce(k,2)-one6th*x(k)*(
106     &           df1/x(k)+df2/(x(k)+df3)+df4*xx*(2.0d0*x(k)+bf)
107     &          +df5/((x(k)+df6)**2+df7) )
108  130 end do
109!$OMP END DO
110
111*     ferromagnetic exchange-energy & potential
112!$OMP DO
113      do 140 k=1,ng
114        xce(k,2)=xce(k,2)+(xf/x(k)**2)
115        xcp(k,2)=xcp(k,2)+for3rd*(xf/x(k)**2)
116  140 end do
117!$OMP END DO
118
119*     spin polarized exchange-correlation potential
120!$OMP DO
121      do 150 k=1,ng
122        rho=dn(k,1)+dn(k,ispin)+dncut
123        zup=2.0d0*dn(k,1)/rho
124        zdw=2.0d0*dn(k,2)/rho
125        f=(zup*(zup**one3rd)+zdw*(zdw**one3rd)-2.0d0)*fc
126        xcp(k,1)=(1.0d0-f)*xcp(k,1)+f*xcp(k,2)
127        df=(zup**one3rd-zdw**one3rd)*(xce(k,2)-xce(k,1))*fd
128        xcp(k,2)=xcp(k,1)-zup*df
129        xcp(k,1)=xcp(k,1)+zdw*df
130        xce(k,1)=xce(k,1)+f*(xce(k,2)-xce(k,1))
131  150 end do
132!$OMP END DO
133
134  200 continue
135
136      return
137      end
138