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