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