By Reinaldo Baretti Machin
e-mail reibaretti2004@yahoo.com
URL www.oocities.org/reibareti2004

c     atomo de berilio-the potentials are phi1s,phi2s,vx1s2s
      dimension phi1s(0:2000),phi2s(0:2000),vx1s2s(0:2000)
      equivalence (dx,h)
      real j1s1s,j1s2s,j2s2s
      psi1s(x)=sqrt(alfa**3/pi)*exp(-alfa*x)
      psi2s(x)=sqrt(alfa2**3/(32.*pi))*(2.-alfa2*x)*exp(-alfa2*x/2.)
c      psi1(x)=anorm*((x+1.)/(x+1.-c1))*(1.+5.*log(1.+x))*exp(-alfa*x)
c      psi1(x)=anorm*(1.-x/6.+x**2/36.)
c     $ *exp(-alfa*x)
      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=4.
      f1s=2.
      f2s=2.
      alfa=3.7
      nstep=2000
      xlim=8.2
      dx=xlim/float(nstep)
      do 100 iter=1,1
      alfa2=.6*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.
      j2s2s=0.
      vnuc=0.
      tk1s=0.
      tk2s=0.
      ex1s2s=0.
      vn2s=0.
      do 30 i=1,nstep
      x=dx*float(i)
      j1s2s=j1s2s+rho1s(x)*x**2*phi2s(i)+ rho1s(x-dx)*(x-dx)**2
     $*phi2s(i-1)
      j1s1s=j1s1s+rho1s(x)*x**2*phi1s(i)+ rho1s(x-dx)*(x-dx)**2
     $*phi1s(i-1)
      j2s2s=j2s2s+rho2s(x)*x**2*phi2s(i)+ rho2s(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      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
      j1s1s=j1s1s*4.*pi*dx/2.
      j1s2s=j1s2s*4.*pi*dx/2.
      j2s2s=j2s2s*4.*pi*dx/2.
      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
      tk1s=alfa**2/2.
      tk2s=alfa2**2/8.
      tkin=f1s*tk1s+f2s*tk2s
      ex1s2s=ex1s2s*4.*pi*dx/2.
      ehtree=tkin+vnuc+j1s1s+f1s*f2s*j1s2s +j2s2s
      ehf=ehtree-2.*ex1s2s
      e2s=tk2s + vn2s + 2.*j1s2s+j2s2s-ex1s2s
      print*,'alfa,alfa2=', alfa,alfa2
      Print*,' j1s1s,j1s2s,j2s2s,ex1s2s=',j1s1s,j1s2s,j2s2s,ex1s2s
      print*,'ecin,vnuc=', tkin ,vnuc
      print*,'e2s,ehtree,ehf=',e2s,ehtree, ehf
      print*,'  '
      alfa=alfa+.1
100   continue
      stop
      end