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