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