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