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