GROUND STATE ENERGY OF LITHIUM

By Reinaldo Baretti Machin
Reibaretti2004@yahoo.com    URL www.oocities.org/reibaretti2004
www.oocities.org/serienumerica






c   lithium variation calculation-with hydrogenic wavefunctions
c and two parameters  alfa , alfa2
      dimension phi1s(0:2000),phi2s(0:2000),vx1s2s(0:2000)
      equivalence (dx,h)
      real j1s1s,j1s2s ,jcoul
      psi1s(x)=sqrt(alfa**3/pi)*exp(-alfa*x)
      psi2s(x)=sqrt(alfa2**3/(32.*pi))*(2.-alfa2*x)*exp(-alfa2*x/2.)
      rho(x)=f1s*psi1s(x)**2+f2s*psi2s(x)**2
      rho1s(x)=psi1s(x)**2
      rho2s(x)=psi2s(x)**2
      rho1s2s(x)=psi1s(x)*psi2s(x)
      pi=2.*asin(1.)
      z=3.
      f1s=2.
      f2s=1.
      alfa=2.7
      nstep=2000
      xlim=7.0
      dx=xlim/float(nstep)
      do 100 iter=1,1
      alfa2=.7*alfa
c       potencial phi mediante ecuacion de Poisson
c  calculo phicero, phi uno calculo de phicero=potencial en el origen
      phicero=0.
      phi2s0=0.
      vx0=0.
      do 10 i=1,nstep
      x=dx*float(i)
      phicero=phicero+rho1s(x)*x+rho1s(x-dx)*(x-dx)
      if(f2s.ge.1.)then
      phi2s0=phi2s0+rho2s(x)*x+rho2s(x-dx)*(x-dx)
      vx0=vx0+rho1s2s(x)*x+rho1s2s(x-dx)*(x-dx)
      endif
10    continue
      phicero=phicero*4.*pi*dx/2.
      phi2s0 =phi2s0*4.*pi*dx/2.
      vx0=vx0*4.*pi*dx/2.
      phi1s(0)=phicero
      phi1s(1)=phi1s(0)
      phi2s(0)=phi2s0
      phi2s(1)=phi2s(0)
      vx1s2s(0)=vx0
      vx1s2s(1)=vx0
c calculo del potencial  resolviendo ec de poisson
c LLeva  factor 4.*Pi ec de Poisson
      do 20 i=2,nstep
      x=dx*float(i)
      phi1s(i)=-dx**2*4.*pi*rho1s(x-dx)-(2.*dx/(x-dx))*
     $ (phi1s(i-1)-phi1s(i-2))+2.*phi1s(i-1)-phi1s(i-2)
      if(f2s.ge.1.)then
      phi2s(i)=-dx**2*4.*pi*rho2s(x-dx)-(2.*dx/(x-dx))*
     $ (phi2s(i-1)-phi2s(i-2))+2.*phi2s(i-1)-phi2s(i-2)
      vx1s2s(i)=-dx**2*4.*pi*rho1s2s(x-dx)-(2.*dx/(x-dx))*
     $ (vx1s2s(i-1)-vx1s2s(i-2))+2.*vx1s2s(i-1)-vx1s2s(i-2)
      endif
20    continue
      J1s1s=0.
      j1s2s=0.
      vnuc=0.
      tk1s=0.
      tk2s=0.
      ex1s2s=0.
      vn2s=0.
      do 30 i=1,nstep
      x=dx*float(i)
      j1s1s=j1s1s+rho1s(x)*x**2*phi1s(i)+ rho1s(x-dx)*(x-dx)**2
     $*phi1s(i-1)
      j1s2s=j1s2s+rho1s(x)*x**2*phi2s(i)+ rho1s(x-dx)*(x-dx)**2
     $*phi2s(i-1)
      vnuc=vnuc+rho(x)*x + rho(x-dx)*(x-dx)
      vn2s=vn2s+ psi2s(x)**2*x + psi2s(x-dx)**2*(x-dx)
c  if psi is not a hydrogenic wavefunction the following step are
c  necessary to calculate kinetic energy
c      tk1s=tk1s-.5*psi1s(x)*x**2*(psi1s(x+dx)-2.*psi1s(x)
c     $+psi1s(x-dx))/h**2- x*psi1s(x)*(psi1s(x+dx)-psi1s(x-dx))/(2.*h)
c      tk2s=tk2s-.5*psi2s(x)*x**2*(psi2s(x+dx)-2.*psi2s(x)
c     $+psi2s(x-dx))/h**2- x*psi2s(x)*(psi2s(x+dx)-psi2s(x-dx))/(2.*h)
      ex1s2s=ex1s2s+rho1s2s(x)*vx1s2s(i)*x**2 +
     $rho1s2s(x-dx)*vx1s2s(i-1)*(x-dx)**2
c      tkele=tkele+.5*psit(x)**2*al*(al+1.)
30    continue
      tk1s=alfa**2/2.
      tk2s=alfa2**2/8.
      j1s1s=j1s1s*4.*pi*dx/2.
      j1s2s=j1s2s*4.*pi*dx/2.
      jcoul=j1s1s+f1s*f2s*j1s2s
      vnuc=-z*vnuc*4.*pi*dx/2.
      vn2s=-z*vn2s*4.*pi*dx/2.
c      tk1s=tk1s*4.*pi*dx
c      tk2s=tk2s*4.*pi*dx
      tkin=f1s*tk1s+f2s*tk2s
      ex1s2s=ex1s2s*4.*pi*dx/2.
      ehtree=tkin+vnuc+j1s1s+f1s*f2s*j1s2s
      ehf=ehtree-ex1s2s
      e2s=tk2s + vn2s + 2.*j1s2s-ex1s2s
c      print*,'ekin1s,tk1s,ekin2s,tk2s=',alfa**2/2.,tk1s, alfa2**2/8.
c     $,tk2s
      print*,'alfa,alfa2,ex12s=', alfa,alfa2,ex1s2s
      print*,'ecin,jcoul,vnuc=', tkin ,jcoul,vnuc
      print*,'e2s,ehtree, etot=',e2s,ehtree, ehf
      print*,'  '
      alfa=alfa+.1
100   continue
      stop
      end 

RUN
alfa,alfa2,ex12s=  2.70000005  1.88999999  0.0758056119
 ecin,jcoul,vnuc=  7.73651314  2.50965858 -17.6142197
 e2s,ehtree, etot= -0.215423375 -7.36804771 -7.44385338