1* 2* $Id$ 3* 4 5* *********************************************** 6* * * 7* * v_thomasfermi * 8* * * 9* *********************************************** 10 11* Computes the Thomas-Fermi potential and energy density 12 13 subroutine v_thomasfermi(n2ft3d,ispin,dn,xcp,xce) 14 implicit none 15 integer n2ft3d 16 integer ispin 17 real*8 dn(n2ft3d,2) 18 real*8 xcp(n2ft3d,2) 19 real*8 xce(n2ft3d) 20 21 22 integer ms,k 23 real*8 n,nup,ndn,eta,f 24 25 real*8 two3rd,five3rd,twotwothirds,CTF,ffac,dncut 26 parameter (two3rd=2.0d0/3.0d0,five3rd=5.0d0/3.0d0) 27 parameter (twotwothirds=1.587401052d0) 28 parameter (CTF=2.871234d0) 29 parameter (ffac=0.851207191959658d0) 30 parameter (dncut=1.0d-30) 31 32 33 call nwpw_timing_start(4) 34 35 call dcopy(n2ft3d,0.0d0,0,xce,1) 36 if (ispin.eq.1) then 37 do k=1,n2ft3d 38 xcp(k,1) = CTF*five3rd*(dn(k,1)+dn(k,ispin))**two3rd 39 xce(k) = xce(k) + CTF*(dn(k,1)+dn(k,ispin))**two3rd 40 end do 41 else 42 do k=1,n2ft3d 43 nup = dn(k,1)+0.5d0*dncut 44 ndn = dn(k,2)+0.5d0*dncut 45 n = nup+ndn 46 eta = (nup-ndn)/n 47 f = ffac*((1.0d0+eta)**five3rd+(1.0d0-eta)**five3rd-2.0d0) 48 xce(k) = xce(k)+(1.0d0+f*(twotwothirds-1.0d0))*CTF*n**two3rd 49 xcp(k,1) = twotwothirds*CTF*five3rd*nup**two3rd 50 xcp(k,2) = twotwothirds*CTF*five3rd*ndn**two3rd 51 end do 52 end if 53 54 call nwpw_timing_end(4) 55 56 return 57 end 58