Chapter 13 - The H₂ molecule

References :

1. Principles of quantum mechanics; nonrelativistic wave mechanics with illustrative applications.  

3.  http://hyperphysics.phy-astr.gsu.edu/hbase/molecule/hmol.html

4.   The hydrogen molecule ion      http://farside.ph.utexas.edu/teaching/qmech/lectures/node129.html

5.  Computational Physics by Paul Devries  , exercise 4.24 page 195.

6.  http://en.wikipedia.org/wiki/Monte_Carlo_integration

7.  S. Weinzierl, Introduction to Monte Carlo methods

8. Schaum's Outline of Numerical Analysis by Francis J. Scheid , chapter 30.

9. http://en.wikipedia.org/wiki/Dihydrogen_cation

10. http://xbeams.chem.yale.edu/~batista/vvv/node32.html

Fig 1. Position vectors for H₂⁺ .

Before considering the H₂ molecule, we will study the H₂ ⁺ molecular ion.  

The hamiltonian is

                                      H = -(1/2) ∆   - 1/ r_A - 1/r_B  + 1 / R                                     .      (1)

Let the origin be at A , and the position of  atom B is along the X-axis ,

                                  r_A = x i +y j + z k     , r_AB = D i    ,

                            r_B = (x-D) +y j  + z k  = r_{A -}  r_AB        .                                          (2)

Assume that the wave function is the symmetric linear combination ot two 1s orbitals

                         Ψ ( r_A ) = φ_1s (r_A )  +   φ_1s (r_B )

                                   =   φ_1s (r_A )  +   φ_1s ( r_{A -}  r_AB )                               ,    (3)

where  φ_1s (r) = (Z³ /π)^1/2 exp(-Zr)  . φ_1s (r_A ) is a 1s orbital centered on atom A and φ_1s (r_B ) is a 1s orbital centered on atom B.

The energy expectation value is  a function of the atoms separation D

 E (D) = < φ_1s (r_A ) + φ_1s ( r_{A -}  r_AB ) / -(1/2) ∆  - 1/ r_A - 1/r_B  + 1 / R  / φ_1s (r_A ) + φ_1s ( r_{A -}  r_AB ) > /< Ψ /Ψ > .  (4)

All the integrals in(4) will be evaluated using the Monte Carlo method of random sampling.

1. A brief review of the method is given. It is taken from  "The hydrogen atom by Monte Carlo method" ,

http://www1.uprh.edu/rbaretti/HydrogenMonteCarlo22ene2010.htm .

In one dimension the average of a function f(x) in the range   a ≤  x ≤ b is given by

                                     f _ave = (1/(b-a)) ∫ f(x) dx                       .    (4)

The average f _ave can also be found by randomly sampling N values  of x in the range a ≤  x ≤ b. One has

                                        f _ave = (1/N) { ∑^N _i=1 f(x_i ) }             .  (5)

So the definite integral is approximated by

                                       ∫ f(x) dx ≈ (b-a) (1/N) { ∑^N _i=1 f(x_i ) }     .   (6) 

The generalization of (6) to three dimensions (x,y,z)  is 

                           I ≈ (x₂-x₁)   (y₂-y₁)   (z₂-z₁) (1/N) { ∑^N _i=1 f(x_i , y_i , z_i ) }      . (7)

 If importance sampling is used an auxiliary function g(x,y,z) is introduced so that g(x,y,z) resembles f(x,y,z). A change of variables is also introduced and the summation is modified to

                               I ≈ (x₂-x₁)   (y₂-y₁)   (z₂-z₁) (1/N) { ∑^N _i=1 f(x_i , y_i , z_i )/g( x_i , y_i , z_i ) }    . (8)

The variables x,y,z have a different range in eq (8) than in  (7).

For example one may have ,  in one dimension , a function f(x) = exp(-x²)  and  the definite integral

                                                              I = ∫^∞ _-∞  exp(-x²) dx = π^1/2   .                                    (9)

A simple Monte Carlo integration would be

                                               I_MC =  (x_final - x_initial ) (1/N) { ∑^N _i=1 f(x_i ) }                       .     (10)

Here x_final = - x_initial and they have to be "large" to represent an infinite limit.

To have importance sampling introduce  g(x)= 1/(1+x²)    and the new variable y

                                          y = ∫ g(x) dx = tan ^-1( x)                                                 (11)

     thus                              dy = g(x) dx = dx/(1+x² )                                   ,             (12)

                                         x= tan(y)  ,       - ∞ ≤ x ≤ ∞     ,   -π/2   ≤ y ≤ π/2      .       (13)

Integral (9) is now in terms of the variable (y) equal to

                                            I = ∫^π/2 _-π/2  exp[ -tan ²(y) ] /(1+tan² (y) ) } dy       ,       (14)

                                              = ∫^π/2 _-π/2 { f(y) /g(y) } dy                                .     (15)

Its Monte Carlo integration has incorporated importance sampling.

                             I_MC =  [ π/2 - ( - π/2   ) ] (1/N) { ∑^N _i=1 f(y_i )/g(y_i)  }        .     (16)

The random sampling in (16) is over the range - π/2  to + π/2 . 

2.  A first step is the integration of  < Ψ /Ψ > = ∫ ( φ_1s (r_A )  +   φ_1s ( r_{A -}  r_AB ) )² dτ = A ( the norm value)   (17)

   One should verify that the Monte Carlo integration assumes the limiting values

                                      lim _{RAB →∞}  < Ψ /Ψ > = 2                             ,                                  (18)

                                      lim _{RAB →0} < Ψ /Ψ >  = 4                                          .                    (19)    

The following two tests bear out the expected result. When D= 10.0 au  ,the normalization constant is 1.98 . When D=0.0 the normalization constant is 3.96

N=16000 sample points and a simple Monte Carlo integration is carried without importance sampling.

nstep,D,norm= 16000  0.10000E+02      0.19803E+01

nstep,D,norm= 16000 0.00000E+00       0.39618E+01


 

Fortran code for normalization test

c FORTRAN code (run for normalization test)
c hydrogen molecular ion H2 + using hydrogen atom by Monte Carlo
dimension x1a(30000),y1a(30000),z1a(30000)
data ndim,nstep/3,16000/
c D = internuclear distance along +x axis
data znuc,x1, D/1.,-3.5, 10./
equivalence(dx,dy,dz)
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(x,y,z)=psi1s(x,y,z)+psi1s(xb,yb,zb)
psicart(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
g(u1,u2,u3)=(1.+tan(u1)**2)*(1.+tan(u2)**2)*
$ (1.+tan(u3)**2)
c f(r1,r2)=psi1s(r1)**2*psi1s(r2)**2*(1./r12)*g(u1,u2,u3,u4,u5,u6)
c f(r)=(-Znuc/r)*psi1s(r)**2*g(u1,u2,u3)
c d2psi(x,y,z)=(psicart(x+dx,y,z)-2.*psicart(x,y,z)+psicart(x-dx,
cc $y,z) + psicart(x,y+dx,z)-2.*psicart(x,y,z)+psicart(x,y-dx,z) +
c $ psicart(x,y,z+dx)-2.*psicart(x,y,z) +psicart(x,y,z-dx) )/dx**2
c f(u) =sin(u)
c d2f(u)=(f(u+delta)-2.*f(u)+f(u-delta))/delta**2
pi=2.*asin(1.)
dx=1.e-3
nm=100000
x2= D + abs(x1)
alength=x2-x1
slopex=x2-x1
bx=x1
y1=x1
y2=-y1
z1=x1
z2=-x1
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
slope2=4.5
b2=-slope2/2.
alength2=slope2
r=22.6781
s= 14.53631
call arand(pi,x1a,r,s,nstep,nm)
r=23.6781
s= 18.53631
call arand(pi,y1a,r,s,nstep,nm)
r=21.8951
s= 15.53631
call arand(pi,z1a,r,s,nstep,nm)
sumpot=0.
sumn=0.
sumkin=0.
do 10 i=1,nstep
c u1= slope*x1a(i) + b
c u2= slope*y1a(i) + b
c u3= slope*z1a(i) + b
c x= slope2*x1a(i) + b2
c y= slope2*y1a(i) + b2
c z= slope2*z1a(i) + b2
x= slopex*x1a(i) + bx
y= slopey*y1a(i) + by
z= slopez*z1a(i) + bz
xb=x-D
yb=y
zb=z
c ra=sqrt(tan(u1)**2 + tan(u2)**2 +tan(u3)**2 )
c rb=sqrt((tan(u1)-D)**2 + tan(u2)**2 +tan(u3)**2 )
ra=sqrt(x**2 + y**2 + z**2 )
rb=sqrt((x-D)**2 + y**2 + z**2 )
c sumpot=sumpot+f(r)
c sumn=sumn+psi1s(r)**2*g(u1,u2,u3)
sumn=sumn+ psih2(x,y,z)**2
c sumkin=sumkin+(-1./2.)*psicart(x,y,z)*d2psi(x,y,z)
10 continue
c defpot=alength**ndim*sumpot/float(nstep)
defnorm=vol*sumn/float(nstep)
c defkin= alength2**ndim*sumkin/float(nstep)
c print*,'nstep=',nstep
c print*,'nstep, pot energy =',nstep, defpot
c print*,'Kin energy =',defkin
c print*,'nstep,anorm =',defnorm
c e=(defkin+defpot)/defnorm
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
print 110 ,nstep, d, defnorm
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('nstep,D,norm= ',i5,2(3x,e12.5))
stop
end

subroutine arand (pi,x,r,s,nstep,nm)
dimension x(30000)
x0=1.
anm=float(nm)
do 10 i=1,nstep
x1=r*x0 + s
c if(i.le.int(float(nstep)/2.))x1=r*x0 + s
c if(i.gt.int(float(nstep)/2.))x1=22.6781*x0 + 14.53631
nx1=int(x1)
j= mod(nx1,nm)
c print*,'j=', j
if(j.lt.0)then
print*,' error j<0 '
goto 20
endif
x(i)= abs(float(j)/anm)
c print*,'i,j=mod(nx1/nm),aj,=',i,j,abs(float(j)/anm)
x0=float(j)
10 continue
20 return
end



3. Integration of the attractive potential         <V> =  < Ψ /-Z/r_A  - Z/r_B /Ψ >  / < Ψ /Ψ > .

Its limits are            lim _{D →∞}  < V > = 2 < φ_1s  / -Z/r_A /  φ_1s > / (2) 

                                                       =  -Z  =  -1                                                                  .     (20)

  Also                     lim _{D →0}  < V > = 2 < 2φ_1s  / -Z/r_A /  2φ_1s > / (4) 

                                                      = -2Z=  -2                                                                .  (21)

The results show  <V>= -1.06 when the inter nuclear separation is  D=12.0 au   and <V>= -1.96 when D=0.

nstep,D,norm,V= 16000 0.12000E+02 0.20314E+01 -0.10630E+01

nstep,D,norm,V= 16000 0.00000E+00 0.38220E+01 -0.19607E+01

Fortran code for part 3.

c FORTRAN code (run for normalization test)
c hydrogen molecular ion H2 + using hydrogen atom by Monte Carlo
dimension x1a(30000),y1a(30000),z1a(30000)
data ndim,nstep/3,16000/
c D = internuclear distance along +x axis
data znuc,x1, D/1.,-2.5, 0./
equivalence(dx,dy,dz)
v(ra,rb)= -znuc/ra -znuc/rb
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(x,y,z)=psi1s(x,y,z)+psi1s(xb,yb,zb)
c psicart(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
c g(u1,u2,u3)=(1.+tan(u1)**2)*(1.+tan(u2)**2)*
c $ (1.+tan(u3)**2)
c f(r1,r2)=psi1s(r1)**2*psi1s(r2)**2*(1./r12)*g(u1,u2,u3,u4,u5,u6)
c f(r)=(-Znuc/r)*psi1s(r)**2*g(u1,u2,u3)
c d2psi(x,y,z)=(psicart(x+dx,y,z)-2.*psicart(x,y,z)+psicart(x-dx,
cc $y,z) + psicart(x,y+dx,z)-2.*psicart(x,y,z)+psicart(x,y-dx,z) +
c $ psicart(x,y,z+dx)-2.*psicart(x,y,z) +psicart(x,y,z-dx) )/dx**2
c f(u) =sin(u)
c d2f(u)=(f(u+delta)-2.*f(u)+f(u-delta))/delta**2
pi=2.*asin(1.)
dx=1.e-3
nm=100000
x2= D + abs(x1)
alength=x2-x1
slopex=x2-x1
bx=x1
y1=x1
y2=-y1
z1=x1
z2=-x1
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
slope2=4.5
b2=-slope2/2.
alength2=slope2
r=22.6781
s= 14.53631
call arand(pi,x1a,r,s,nstep,nm)
r=23.6781
s= 18.53631
call arand(pi,y1a,r,s,nstep,nm)
r=21.8951
s= 15.53631
call arand(pi,z1a,r,s,nstep,nm)
sumpot=0.
sumn=0.
sumkin=0.
do 10 i=1,nstep
c u1= slope*x1a(i) + b
c u2= slope*y1a(i) + b
c u3= slope*z1a(i) + b
c x= slope2*x1a(i) + b2
c y= slope2*y1a(i) + b2
c z= slope2*z1a(i) + b2
x= slopex*x1a(i) + bx
y= slopey*y1a(i) + by
z= slopez*z1a(i) + bz
xb=x-D
yb=y
zb=z
c ra=sqrt(tan(u1)**2 + tan(u2)**2 +tan(u3)**2 )
c rb=sqrt((tan(u1)-D)**2 + tan(u2)**2 +tan(u3)**2 )
ra=sqrt(x**2 + y**2 + z**2 )
rb=sqrt((x-D)**2 + y**2 + z**2 )
sumpot=sumpot+v(ra,rb)*psih2(x,y,z)**2
c sumn=sumn+psi1s(r)**2*g(u1,u2,u3)
sumn=sumn+ psih2(x,y,z)**2
c sumkin=sumkin+(-1./2.)*psicart(x,y,z)*d2psi(x,y,z)
10 continue
anorm=vol*sumn/float(nstep)
sumpot=(1./anorm)*vol*sumpot/float(nstep)
c defkin= alength2**ndim*sumkin/float(nstep)
c print*,'nstep=',nstep
c print*,'nstep, pot energy =',nstep, anorm
c print*,'Kin energy =',defkin
c print*,'nstep,anorm =',defnorm
c e=(defkin+defpot)/defnorm
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
print 110 ,nstep, d, anorm,sumpot
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('nstep,D,norm,V= ',i5,3(3x,e12.5))
stop
end

subroutine arand (pi,x,r,s,nstep,nm)
dimension x(30000)
x0=1.
anm=float(nm)
do 10 i=1,nstep
x1=r*x0 + s
c if(i.le.int(float(nstep)/2.))x1=r*x0 + s
c if(i.gt.int(float(nstep)/2.))x1=22.6781*x0 + 14.53631
nx1=int(x1)
j= mod(nx1,nm)
c print*,'j=', j
if(j.lt.0)then
print*,' error j<0 '
goto 20
endif
x(i)= abs(float(j)/anm)
c print*,'i,j=mod(nx1/nm),aj,=',i,j,abs(float(j)/anm)
x0=float(j)
10 continue
20 return
end



4. Tests of the kinetic energy average   < Ψ / -(1/2) ∆  / Ψ > / < Ψ / Ψ >

If D becomes large  , < Ψ / Ψ > ≈ 2 ,  <Ψ / -(1/2) ∆  / Ψ > ≈  2 <  φ_1s/ -(1/2) ∆  /   φ_1s> ≈ 2(1/2) =1. Therefore

                                   < Ψ / -(1/2) ∆  / Ψ > / < Ψ / Ψ > ≈ 1/2        .                          (22)

If D goes to zero   , < Ψ / Ψ > ≈ 4 ,  <Ψ / -(1/2) ∆  / Ψ > ≈   <2  φ_1s/ -(1/2) ∆  / 2 φ_1s> ≈ 4(1/2) =2.  Therefore

                                  < Ψ / -(1/2) ∆  / Ψ > / < Ψ / Ψ > ≈ 1/2   .                        (23)

The results obtained with the modified FORTRAN code given below are ,

 nstep,D= 16000 0.
norm,V, ekin= 0.38220E+01 -0.19607E+01  0.48145E+00

 nstep,D= 16000 12.
norm,V, ekin= 0.18397E+01 -0.11435E+01   0.55467E+00   .

 

FORTRAN code for part 4.

c FORTRAN code (run for normalization test)
c hydrogen molecular ion H2 + using hydrogen atom by Monte Carlo
dimension x1a(30000),y1a(30000),z1a(30000)
data ndim,nstep/3,16000/
c D = internuclear distance along +x axis
data znuc,x1, D/1.,-2.5, 0./
equivalence(dx,dy,dz)
v(ra,rb)= -znuc/ra -znuc/rb
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(x,y,z)=psi1s(x,y,z)+psi1s(x-d,y,z)
c psicart(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
c g(u1,u2,u3)=(1.+tan(u1)**2)*(1.+tan(u2)**2)*
c $ (1.+tan(u3)**2)
c f(r1,r2)=psi1s(r1)**2*psi1s(r2)**2*(1./r12)*g(u1,u2,u3,u4,u5,u6)
c f(r)=(-Znuc/r)*psi1s(r)**2*g(u1,u2,u3)
d2psih2(x,y,z)=(psih2(x+dx,y,z)-2.*psih2(x,y,z)+psih2(x-dx,
$y,z) + psih2(x,y+dx,z)-2.*psih2(x,y,z)+psih2(x,y-dx,z) +
$ psih2(x,y,z+dx)-2.*psih2(x,y,z) +psih2(x,y,z-dx) )/dx**2
pi=2.*asin(1.)
dx=1.e-3
nm=100000
x2= D + abs(x1)
alength=x2-x1
slopex=x2-x1
bx=x1
y1=x1
y2=-y1
z1=x1
z2=-x1
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
slope2=4.5
b2=-slope2/2.
alength2=slope2
r=22.6781
s= 14.53631
call arand(pi,x1a,r,s,nstep,nm)
r=23.6781
s= 18.53631
call arand(pi,y1a,r,s,nstep,nm)
r=21.8951
s= 15.53631
call arand(pi,z1a,r,s,nstep,nm)
sumpot=0.
sumn=0.
sumkin=0.
do 10 i=1,nstep
c u1= slope*x1a(i) + b
c u2= slope*y1a(i) + b
c u3= slope*z1a(i) + b
c x= slope2*x1a(i) + b2
c y= slope2*y1a(i) + b2
c z= slope2*z1a(i) + b2
x= slopex*x1a(i) + bx
y= slopey*y1a(i) + by
z= slopez*z1a(i) + bz
xb=x-D
yb=y
zb=z
c ra=sqrt(tan(u1)**2 + tan(u2)**2 +tan(u3)**2 )
c rb=sqrt((tan(u1)-D)**2 + tan(u2)**2 +tan(u3)**2 )
ra=sqrt(x**2 + y**2 + z**2 )
rb=sqrt((x-D)**2 + y**2 + z**2 )
sumpot=sumpot+v(ra,rb)*psih2(x,y,z)**2
c sumn=sumn+psi1s(r)**2*g(u1,u2,u3)
sumn=sumn+ psih2(x,y,z)**2
sumkin=sumkin+(-1./2.)*psih2(x,y,z)*d2psih2(x,y,z)
10 continue
anorm=vol*sumn/float(nstep)
sumpot=(1./anorm)*vol*sumpot/float(nstep)
ekin= (1./anorm)*vol*sumkin/float(nstep)
c print*,'nstep=',nstep
c print*,'nstep, pot energy =',nstep, anorm
c print*,'Kin energy =',defkin
c print*,'nstep,anorm =',defnorm
c e=(defkin+defpot)/defnorm
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
print*,'nstep,D=',nstep, d
print 110 ,anorm , sumpot, ekin
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('norm,V, ekin= ',3(3x,e12.5))
stop
end

subroutine arand (pi,x,r,s,nstep,nm)
dimension x(30000)
x0=1.
anm=float(nm)
do 10 i=1,nstep
x1=r*x0 + s
c if(i.le.int(float(nstep)/2.))x1=r*x0 + s
c if(i.gt.int(float(nstep)/2.))x1=22.6781*x0 + 14.53631
nx1=int(x1)
j= mod(nx1,nm)
c print*,'j=', j
if(j.lt.0)then
print*,' error j<0 '
goto 20
endif
x(i)= abs(float(j)/anm)
c print*,'i,j=mod(nx1/nm),aj,=',i,j,abs(float(j)/anm)
x0=float(j)
10 continue
20 return
end


 

5. Energy of  H₂ ⁺

To find the energy E(D) -see eq.(4), we use Fortran code #5.

The reader should compare our treatment with those given in Chapter XII of  reference 2 and also with Ref. 9 .

Fig 2. Ground state energy curve of H₂⁺   . The minima of energy is about  -.582 -(-.5) = -.082 au = -2.23 eV at D=2.40 au .

Fig 3. Plot of Ψ² along X axis  of  the ground state wave function, D=2.4 au.

Fig 4. Plot of Ψ² along X axis  of  the anti bonding wave function, D=2.4 au. Notice that the probability density is

zero at midpoint.

Fortran code for plots of  Ψ²

c hydrogen molecular ion H2 + using hydrogen atom by Monte Carlo
c D = internuclear distance along +x axis
data znuc,nx,x1/1.,70,-2.5 /
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(x,y,z)=psi1s(x,y,z)-psi1s(x-d,y,z)
pi=2.*asin(1.)
y=0.
z=0.
d=2.4
x2=abs(x1)+d
dx=(x2-x1)/float(nx)
do 10 i=0,nx
x=x1+dx*float(i)
print 100,x , psih2(x,y,z)**2
10 continue
100 format('x, psih2^2=',2(3x,e11.4))
stop
end
 

Fig 5.   Energies for H₂ ^{+ ,}  taken from reference  9. The lowest curve corresponds to the ground state energy.

FORTRAN  code (#5)   for Part 5.

c hydrogen molecular ion H2 + using hydrogen atom by Monte Carlo
dimension x1a(30000),y1a(30000),z1a(30000)
data ndim,nstep/3,16000/
c D = internuclear distance along +x axis
data znuc,x1/1.,-2.5 /
data nd,d1,d2/20,0.5,5.5/
equivalence(dx,dy,dz)
v(ra,rb)= -znuc/ra -znuc/rb
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(x,y,z)=psi1s(x,y,z)+psi1s(x-d,y,z)
c psicart(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
c g(u1,u2,u3)=(1.+tan(u1)**2)*(1.+tan(u2)**2)*
c $ (1.+tan(u3)**2)
c f(r1,r2)=psi1s(r1)**2*psi1s(r2)**2*(1./r12)*g(u1,u2,u3,u4,u5,u6)
c f(r)=(-Znuc/r)*psi1s(r)**2*g(u1,u2,u3)
c second derivatives of psih2
d2psih2(x,y,z)=(psih2(x+dx,y,z)-2.*psih2(x,y,z)+psih2(x-dx,
$y,z) + psih2(x,y+dx,z)-2.*psih2(x,y,z)+psih2(x,y-dx,z) +
$ psih2(x,y,z+dx)-2.*psih2(x,y,z) +psih2(x,y,z-dx) )/dx**2
pi=2.*asin(1.)
deltad=(d2-d1)/float(nd)
dx=1.e-3
nm=100000
do 50 id=0,nd
d=d1+deltad*float(id)
x2= D + abs(x1)
slopex=x2-x1
bx=x1
y1=x1
y2=-y1
z1=x1
z2=-x1
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
r=22.6781
s= 14.53631
call arand(pi,x1a,r,s,nstep,nm)
r=23.6781
s= 18.53631
call arand(pi,y1a,r,s,nstep,nm)
r=21.8951
s= 15.53631
call arand(pi,z1a,r,s,nstep,nm)
sumpot=0.
sumn=0.
sumkin=0.
do 10 i=1,nstep
x= slopex*x1a(i) + bx
y= slopey*y1a(i) + by
z= slopez*z1a(i) + bz
xb=x-D
yb=y
zb=z
c ra=sqrt(tan(u1)**2 + tan(u2)**2 +tan(u3)**2 )
c rb=sqrt((tan(u1)-D)**2 + tan(u2)**2 +tan(u3)**2 )
ra=sqrt(x**2 + y**2 + z**2 )
rb=sqrt((x-D)**2 + y**2 + z**2 )
sumpot=sumpot+v(ra,rb)*psih2(x,y,z)**2
c sumn=sumn+psi1s(r)**2*g(u1,u2,u3)
sumn=sumn+ psih2(x,y,z)**2
sumkin=sumkin+(-1./2.)*psih2(x,y,z)*d2psih2(x,y,z)
10 continue
anorm=vol*sumn/float(nstep)
sumpot=(1./anorm)*vol*sumpot/float(nstep)
ekin= (1./anorm)*vol*sumkin/float(nstep)
Vnuclear=1./d
c print*,'nstep=',nstep
c print*,'nstep, pot energy =',nstep, anorm
c print*,'Kin energy =',defkin
c print*,'nstep,anorm =',defnorm
c e=(defkin+defpot)/defnorm
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
c print*,'nstep,anorm,D=',nstep,anorm, d
c print 110 ,vnuclear, sumpot, ekin
et=ekin+sumpot+vnuclear
print 115,d,et
50 continue
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('vnuclear,V, ekin= ',3(3x,e12.5))
115 format('d,etotal= ',2(3x,e11.4))
stop
end

subroutine arand (pi,x,r,s,nstep,nm)
dimension x(30000)
x0=1.
anm=float(nm)
do 10 i=1,nstep
x1=r*x0 + s
c if(i.le.int(float(nstep)/2.))x1=r*x0 + s
c if(i.gt.int(float(nstep)/2.))x1=22.6781*x0 + 14.53631
nx1=int(x1)
j= mod(nx1,nm)
c print*,'j=', j
if(j.lt.0)then
print*,' error j<0 '
goto 20
endif
x(i)= abs(float(j)/anm)
c print*,'i,j=mod(nx1/nm),aj,=',i,j,abs(float(j)/anm)
x0=float(j)
10 continue
20 return
end



6. Anti bonding state of H₂⁺

An anti bonding state results from the wave function

                                         Ψ ( r_A ) = φ_1s (r_A )  -   φ_1s (r_B )                         .   (24)

Fig 6. Anti bonding energy curve . Compare with the lowest green curve in Fig 5.



7. Ground state wave function of   H₂  in the Heitler -London ( valence bond method) approach (ref. 10 ).

Fig 7.   Position vectors for the hydrogen molecule.

The hamiltonian is

                         H = -(1/2) ∆₁   - 1/ r_A1 - 1/r_B1   -(1/2) ∆₂   - 1/ r_A2 - 1/r_B2 + 1/r₁₂ + 1 / R     ,    (25-a )

                            ≡  H₀₁ + H₀₂  + 1/r₁₂  +  1/R                                                               ,        (25-b)

where    H₀₁  =  -(1/2) ∆₁   - 1/ r_A1 - 1/r_B1  and  H₀₂  = -(1/2) ∆₂   - 1/ r_A2 - 1/r_B2 .

      One  plausible approach to the problem is to let the ground state be represented by the spatially symmetric  wave function ,                                  

       Ψ (1,2) ~  φ_1s (r_A1 ) φ_1s (r_B2 ) + φ_1s (r_B1 ) φ_1s (r_A2 )                                 .  ( 26  )

In the first term electron number #1 "belongs " to atom A and in the second term electron #1 is attached to atom B. The opposite can be said about electron number 2.

The complete wave function has to be anti symmetric. This implies that (26 ) must  be multiplied by an anti symmetric spinor wave function. The un normalized wave function is ,

       Ψ (1,2) = [ φ_1s (r_A1 ) φ_1s (r_B2 ) + φ_1s (r_B1 ) φ_1s (r_A2 )] {α (1)β(2) -   β(1)α(2) }        (Valence bond)      (  27  )

In the molecular orbit approximation there are more terms in the spatial part of the wave function

Ψ (1,2) = {φ_1s (r_A1 )  + φ_1s (r_B1 ) } {φ_1s (r_A2 ) + φ_1s (r_B2 ) } [α (1)β(2) -   β(1)α(2) ]   (MO-molecular orbit)    (28) 

           E(R) = 2 < Ψ /H₀₁ /  Ψ   >/ <Ψ /Ψ >    +   < Ψ /(1/r₁₂) /Ψ  > / <Ψ /Ψ >    +      1/R              .   (29)

FORTRAN code  for H₂ molecule (Fig 8) - see below-  was used to integrate (29) by Monte Carlo methods. The result is shown in Fig. 8.  E_min ≈  -1.17 au  at R ≈ 1.5 au. The asymptotic state energy is that of two separated atoms, each with energy -0.5au .

Fig 8. Energy curve of H₂ using simple Monte Carlo integration.

Fig.  9 from http://hyperphysics.phy-astr.gsu.edu/hbase/molecule/hmol.html

In Fig 9 the energy minima is -4.52 eV=  -0.166 au and 0.074 nm= 1.40 au.

A second Monte Carlo code with importance sampling is employed  ( see FORTRAN code with importance sampling below).

The energy plot is shown in Fig 10.

Fig 10. For what it is worth. The energy minimum is E = -1.17 au   at R= 1.60 au .

We reproduce a recent abstract of a paper published in the Interantional Journal of Quantum Chemistry dealing with the

energy of  H₂ using the method of basis sets.

Lastly we investigate the energy curve for the wave functions with anti symmetric spatial part but symmetric spinor part.

It could be anyone of the following triplet

Ψ (1,2) = [ φ_1s (r_A1 ) φ_1s (r_B2 ) - φ_1s (r_B1 ) φ_1s (r_A2 )] {α (1)β(2) +   β(1)α(2) }   ,   (30)

or          [ φ_1s (r_A1 ) φ_1s (r_B2 ) - φ_1s (r_B1 ) φ_1s (r_A2 )] {α (1)α(2) }       ,                 (31)

or           [ φ_1s (r_A1 ) φ_1s (r_B2 ) - φ_1s (r_B1 ) φ_1s (r_A2 )] {(β1) β(2) }  .                    (32)

As will be shown next these are anti binding states.

For any two hydrogen atoms interacting in their ground states, the probability of forming a bound state is  1/4  while the probability of not forming a stable molecule is 3/4.

Fig 11. Anti binding state.

FORTRAN code  for H₂ molecule (Fig 8)

c FORTRAN code for H2  by Monte Carlo
dimension x1aa(60000),y1aa(60000),z1aa(60000)
dimension x2aa(60000),y2aa(60000),z2aa(60000)
data nstep/60000/
c D = internuclear distance along + X axis
c data znuc,x1/1.,-2.5 /
data znuc,x1/1.,-2.5 /
data nd,d1,d2/20,.5,5./
equivalence(dx,dy,dz)
v(r1a,r1b,r2a,r2b,r12)= -znuc/r1a -znuc/r1b -znuc/r2a -znuc/r2b
$ + 1./r12
c vee(r12)=1./r12
psi(r)=sqrt(znuc**3/pi)*exp(-znuc*r)
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(xa,ya,za,xb,yb,zb)=psi1s(xa,ya,za)*psi1s(xb-d,yb,zb)+
$ psi1s(xa-d,ya,za)*psi1s(xb,yb,zb)
c second derivatives of psih2
d2psih2(xa,ya,za,xb,yb,zb)=(psih2(xa+dx,ya,za,xb,yb,zb)
$ -6.*psih2(xa,ya,za,xb,yb,zb)+
$ psih2(xa-dx,ya,za,xb,yb,zb) +
$ psih2(xa,ya+dy,za,xb,yb,zb)
$ +psih2(xa,ya-dy,za,xb,yb,zb) +
$ psih2(xa,ya,za+dx,xb,yb,zb) +
$ psih2(xa,ya,za-dx,xb,yb,zb) )/dx**2
c d2psih2(x,y,z)=(psih2(x+dx,y,z)-2.*psih2(x,y,z)+psih2(x-dx,
c $y,z) + psih2(x,y+dx,z)-2.*psih2(x,y,z)+psih2(x,y-dx,z) +
c $ psih2(x,y,z+dx)-2.*psih2(x,y,z) +psih2(x,y,z-dx) )/dx**2
pi=2.*asin(1.)
deltad=(d2-d1)/float(nd)
dx=1.e-3
nm=100000
do 50 id=1,nd
d=d1+deltad*float(id-1)
x2= D + abs(x1)
slopex=x2-x1
bx=x1
y1=x1
y2=-y1
z1=x1
z2=-x1
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)**2*(y2-y1)**2*(z2-z1)**2
r=5.6781
s= 14.53631
if(id.eq.1)call arand(pi,x1aa,r,s,nstep,nm)
r=21.3781
s= 14.98763
IF(ID.EQ.1)call arand(pi,y1aa,r,s,nstep,nm)
r=7.951
s= 5.53631
if(id.eq.1)call arand(pi,z1aa,r,s,nstep,nm)
r=24.8951
s= 17.53631
if(id.eq.1)call arand(pi,x2aa,r,s,nstep,nm)
r=3.6517
s= 9.6316
if(id.eq.1)call arand(pi,y2aa,r,s,nstep,nm)
r=2.3148
s=6.7316
if(id.eq.1)call arand(pi,z2aa,r,s,nstep,nm)
sumpot=0.
sumn=0.
sumkin=0.
do 10 i=1,nstep
x1a= slopex*x1aa(i) + bx
y1a= slopey*y1aa(i) + by
z1a= slopez*z1aa(i) + bz
x1b=x1a-d
y1b=y1a
z1b=z1a
x2a= slopex*x2aa(i) + bx
y2a= slopey*y2aa(i) + by
z2a= slopez*z2aa(i) + bz
x2b=x2a-d
y2b=y2a
z2b=z2a
r12=((x2a-x1a)**2+(y2a-y1a)**2 +(z2a-z1a)**2 )**.5
r1a=sqrt(x1a**2 + y1a**2 + z1a**2 )
r1b=sqrt(x1b**2 + y1b**2 + z1b**2 )
r2a=sqrt(x2a**2 + y2a**2 + z2a**2 )
r2b=sqrt(x2b**2 + y2b**2 + z2b**2 )
sumpot=sumpot+ v(r1a,r1b,r2a,r2b,r12)*
$ psih2(x1a,y1a,z1a,x2a,y2a,z2a)**2
c sumn=sumn+psi1s(r)**2*g(u1,u2,u3)
c sumn=sumn+ psih2(x1a,y1a,z1a,x2a,y2a,z2a)**2
sumn=sumn+(psi(r1a)*psi(r2b)+psi(r1b)*psi(r2a))**2
sumkin=sumkin+ 2.*(-1./2.)*psih2(x1a,y1a,z1a,x2a,y2a,z2a)*
$ d2psih2(x1a,y1a,z1a,x2a,y2a,z2a)
10 continue
anorm=vol*sumn/float(nstep)
sumpot=(1./anorm)*vol*sumpot/float(nstep)
ekin= (1./anorm)*vol*sumkin/float(nstep)
Vnuclear=1./d
c print*,'nstep,anorm,potenergy=',nstep ,anorm,sumpot
c print*,'Vnuclear,Kinetic energy =',vnuclear, ekin
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
c print*,'nstep,anorm,D=',nstep,anorm, d
c print 110 ,vnuclear, sumpot, ekin
et=ekin+sumpot+vnuclear
print 115,d,et
50 continue
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('vnuclear,V, ekin= ',3(3x,e12.5))
115 format('d,etotal= ',2(3x,e11.4))
stop
end

subroutine arand (pi,x,r,s,nstep,nm)
dimension x(60000)
x0=1.
anm=float(nm)
do 10 i=1,nstep
x1=r*x0 + s
c if(i.le.int(float(nstep)/2.))x1=r*x0 + s
c if(i.gt.int(float(nstep)/2.))x1=22.6781*x0 + 14.53631
nx1=int(x1)
j= mod(nx1,nm)
c print*,'j=', j
if(j.lt.0)then
print*,' error j<0 '
goto 20
endif
x(i)= abs(float(j)/anm)
c print*,'i,j=mod(nx1/nm),aj,=',i,j,abs(float(j)/anm)
x0=float(j)
10 continue
20 return
end



FORTRAN code for H₂ with importance sampling (fig 10.)

c hydrogen molecule H2 tests with Monte Carlo 5 febrero 2010
dimension x1aa(60000),y1aa(60000),z1aa(60000)
dimension x2aa(60000),y2aa(60000),z2aa(60000)
data nstep/40000/
c D = internuclear distance along + X axis
c data znuc,x1/1.,-2.5 /
data znuc/1. /
data nd,d1,d2/50,.5,5.5/
equivalence(dx,dy,dz)
v(r1a,r1b,r2a,r2b,r12)= -znuc/r1a -znuc/r1b -znuc/r2a -znuc/r2b
$ + 1./r12
vee(r12)=1./r12
psi(r)=sqrt(znuc**3/pi)*exp(-znuc*r)
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(u1a,u2a,u3a,v1a,v2a,v3a)=psi1s(u1a,u2a,u3a)*
$ psi1s(v1a-d,v2a,v3a )+
$ psi1s(u1a-d,u2a,u3a)*psi1s(v1a,v2a,v3a )
c second derivatives of psih2

c d2psih2(xa,ya,za,xb,yb,zb)=(psih2(xa+dx,ya,za,xb,yb,zb)
c $ + psih2(xa-dx,ya,za,xb,yb,zb)
c $ + psih2(xa,ya+dy,za,xb,yb,zb)
c $ +psih2(xa,ya-dy,za,xb,yb,zb)
c $ + psih2(xa,ya,za+dx,xb,yb,zb)
c $ + psih2(xa,ya,za-dx,xb,yb,zb)
c $ -6.*psih2(xa,ya,za,xb,yb,zb) )/dx**2
d2psih2(u1a,u2a,u3a,v1a,v2a,v3a)=
$ ( psih2(u1a+dx,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a-dx,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a+dy,u3a,v1a,v2a,v3a)
$ +psih2(u1a,u2a-dy,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a+dx,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a-dx,v1a,v2a,v3a )
$ -6.*psih2(u1a,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a+dx,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a-dx,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a,v2a+dy,v3a)
$ +psih2(u1a,u2a,u3a,v1a,v2a-dy,v3a)
$ + psih2(u1a,u2a,u3a,v1a,v2a,v3a+dz)
$ + psih2(u1a,u2a,u3a,v1a,v2a,v3a-dz)
$ -6.*psih2(u1a,u2a,u3a,v1a,v2a,v3a) )/dx**2

pi=2.*asin(1.)
deltad=(d2-d1)/float(nd)
dx=1.e-3
do 50 id=1,nd
d=d1+deltad*float(id-1)
x1=-pi/2.
x2= pi/2.
slopex=x2-x1
bx=x1
y1=x1
y2=x2
z1=x1
z2=x2
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
vol2=vol**2
pwr=1./2.
seed=sqrt(pi)
if(id.eq.1) call aleat(nstep,pwr,seed,pi,x1aa)
seed=2.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,y1aa)
seed=3.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,z1aa)
seed=5.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,x2aa)
seed=7.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,y2aa)
seed=11.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,z2aa)
sumpot=0.
sumn=0.
sumkin=0.
ajee=0.
do 10 i=1,nstep
c r1a position of electron #1 with respect to atom A at origin
x1a= slopex*x1aa(i) + bx
y1a= slopey*y1aa(i) + by
z1a= slopez*z1aa(i) + bz
u1a=tan(x1a)
u2a=tan(y1a)
u3a=tan(z1a)
c print*,'i,x1a,u1a=',i,x1a,u1a
c vector r1b
x1b=x1a-d
y1b=y1a
z1b=z1a
c r2a position of electron #2 with respect to atom A at origin
x2a= slopex*x2aa(i) + bx
y2a= slopey*y2aa(i) + by
z2a= slopez*z2aa(i) + bz
v1a=tan(x2a)
v2a=tan(y2a)
v3a=tan(z2a)
c vector r2b
x2b=x2a-d
y2b=y2a
z2b=z2a
r12=sqrt((u1a-v1a)**2 +(u2a-v2a)**2 +(u3a-v3a)**2 )
r1a=sqrt(u1a**2 + u2a**2 + u3a**2 )
r1b=sqrt((u1a-d)**2 + u2a**2 + u3a**2 )
r2a=sqrt(v1a**2 + v2a**2 + v3a**2 )
r2b=sqrt((v1a-d)**2 + v2a**2 + v3a**2 )
sumpot=sumpot + v(r1a,r1b,r2a,r2b,r12)*(psi(r1a)*psi(r2b)+
$psi(r1b)*psi(r2a))**2*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
$(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)
c sumpot=sumpot + v(r1a,r1b,r2a,r2b,r12)*
c $ (psi(r1a)*psi(r2b)+psi(r1b)*psi(r2a))**2
c sumpot=sumpot +(-znuc/r1a)*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)
c $ *psi(r1a)**2
c sumn=sumn+ psih2(x1a,y1a,z1a,x2a,y2a,z2a)**2
c sumn=sumn+(psi(r1a)*psi(r2b)+psi(r1b)*psi(r2a))**2
sumn=sumn+(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
$(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)*(psi(r1a)*psi(r2b)+
$psi(r1b)*psi(r2a))**2
sumkin=sumkin + (-1./2.)*psih2(u1a,u2a,u3a,v1a,v2a,v3a)*
$ d2psih2(u1a,u2a,u3a,v1a,v2a,v3a)*(1.+u1a**2)*(1.+u2a**2)*
$(1.+u3a**2)*(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)

c sumkin=sumkin+(-1./2.)*psi(r1a)*d2psih2(x1a,y1a,z1a)
c ajee=ajee+vee(r12)*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
c $(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)*psi(r1a)**2*psi(r2a)**2
10 continue
anorm=vol2*sumn/float(nstep)
c sumpot=(1./anorm)*vol*sumpot/float(nstep)
sumpot= (1./anorm)*vol2*sumpot/float(nstep)
ajee=(vol2/anorm)*ajee/float(nstep)
c ekin= (1./anorm)*vol*sumkin/float(nstep)
ekin=(1./anorm)*vol2*sumkin/float(nstep)
Vnuclear=1./d
c print*,'nstep,potenergy=',nstep ,sumpot
c print*,'Kinetic energy =', ekin
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
c print*,'nstep,d,anorm=',nstep,d, anorm
c print*,'ajee,5./8.=', ajee ,5./8.
et=ekin + sumpot + vnuclear
print 110 ,d, sumpot, ekin,et
c print 115,d,et
50 continue
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('d,V,ek,et=',4(3x,e10.3))
115 format('d,etotal=',2(3x,e11.4))
stop
end

subroutine aleat(nstep,pwr,seed,pi,x)
dimension x(60000)
x0=seed
do 10 i=1,nstep
y=x0**pwr +sqrt(5.)
rnd=y-int(y)
x(i)=rnd
c print*,'rnd=',rnd
x0=rnd*pi*1.e4
10 continue
return
end


Fortran code for anti binding state (Fig. 11)

c hydrogen molecule H2 tests with Monte Carlo 5 febrero 2010
dimension x1aa(60000),y1aa(60000),z1aa(60000)
dimension x2aa(60000),y2aa(60000),z2aa(60000)
data nstep/40000/
c D = internuclear distance along + X axis
c data znuc,x1/1.,-2.5 /
data znuc/1. /
data nd,d1,d2/50,.5,5.5/
equivalence(dx,dy,dz)
v(r1a,r1b,r2a,r2b,r12)= -znuc/r1a -znuc/r1b -znuc/r2a -znuc/r2b
$ + 1./r12
vee(r12)=1./r12
psi(r)=sqrt(znuc**3/pi)*exp(-znuc*r)
psi1s(x,y,z)=sqrt(znuc**3/pi)*exp(-znuc*sqrt(x**2+y**2+z**2))
psih2(u1a,u2a,u3a,v1a,v2a,v3a)=psi1s(u1a,u2a,u3a)*
$ psi1s(v1a-d,v2a,v3a )-psi1s(u1a-d,u2a,u3a)*psi1s(v1a,v2a,v3a )
c second derivatives of psih2

c d2psih2(xa,ya,za,xb,yb,zb)=(psih2(xa+dx,ya,za,xb,yb,zb)
c $ + psih2(xa-dx,ya,za,xb,yb,zb)
c $ + psih2(xa,ya+dy,za,xb,yb,zb)
c $ +psih2(xa,ya-dy,za,xb,yb,zb)
c $ + psih2(xa,ya,za+dx,xb,yb,zb)
c $ + psih2(xa,ya,za-dx,xb,yb,zb)
c $ -6.*psih2(xa,ya,za,xb,yb,zb) )/dx**2
d2psih2(u1a,u2a,u3a,v1a,v2a,v3a)=
$ ( psih2(u1a+dx,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a-dx,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a+dy,u3a,v1a,v2a,v3a)
$ +psih2(u1a,u2a-dy,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a+dx,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a-dx,v1a,v2a,v3a )
$ -6.*psih2(u1a,u2a,u3a,v1a,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a+dx,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a-dx,v2a,v3a)
$ + psih2(u1a,u2a,u3a,v1a,v2a+dy,v3a)
$ +psih2(u1a,u2a,u3a,v1a,v2a-dy,v3a)
$ + psih2(u1a,u2a,u3a,v1a,v2a,v3a+dz)
$ + psih2(u1a,u2a,u3a,v1a,v2a,v3a-dz)
$ -6.*psih2(u1a,u2a,u3a,v1a,v2a,v3a) )/dx**2

pi=2.*asin(1.)
deltad=(d2-d1)/float(nd)
dx=1.e-3
do 50 id=1,nd
d=d1+deltad*float(id-1)
x1=-pi/2.
x2= pi/2.
slopex=x2-x1
bx=x1
y1=x1
y2=x2
z1=x1
z2=x2
slopey=y2-y1
by=y1
slopez=z2-z1
bz=y1
vol=(x2-x1)*(y2-y1)*(z2-z1)
vol2=vol**2
pwr=1./2.
seed=sqrt(pi)
if(id.eq.1) call aleat(nstep,pwr,seed,pi,x1aa)
seed=2.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,y1aa)
seed=3.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,z1aa)
seed=5.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,x2aa)
seed=7.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,y2aa)
seed=11.
if(id.eq.1)call aleat(nstep,pwr,seed,pi,z2aa)
sumpot=0.
sumn=0.
sumkin=0.
ajee=0.
do 10 i=1,nstep
c r1a position of electron #1 with respect to atom A at origin
x1a= slopex*x1aa(i) + bx
y1a= slopey*y1aa(i) + by
z1a= slopez*z1aa(i) + bz
u1a=tan(x1a)
u2a=tan(y1a)
u3a=tan(z1a)
c print*,'i,x1a,u1a=',i,x1a,u1a
c vector r1b
x1b=x1a-d
y1b=y1a
z1b=z1a
c r2a position of electron #2 with respect to atom A at origin
x2a= slopex*x2aa(i) + bx
y2a= slopey*y2aa(i) + by
z2a= slopez*z2aa(i) + bz
v1a=tan(x2a)
v2a=tan(y2a)
v3a=tan(z2a)
c vector r2b
x2b=x2a-d
y2b=y2a
z2b=z2a
r12=sqrt((u1a-v1a)**2 +(u2a-v2a)**2 +(u3a-v3a)**2 )
r1a=sqrt(u1a**2 + u2a**2 + u3a**2 )
r1b=sqrt((u1a-d)**2 + u2a**2 + u3a**2 )
r2a=sqrt(v1a**2 + v2a**2 + v3a**2 )
r2b=sqrt((v1a-d)**2 + v2a**2 + v3a**2 )
sumpot=sumpot + v(r1a,r1b,r2a,r2b,r12)*(psi(r1a)*psi(r2b)-
$psi(r1b)*psi(r2a))**2*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
$(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)
c sumpot=sumpot + v(r1a,r1b,r2a,r2b,r12)*
c $ (psi(r1a)*psi(r2b)+psi(r1b)*psi(r2a))**2
c sumpot=sumpot +(-znuc/r1a)*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)
c $ *psi(r1a)**2
c sumn=sumn+ psih2(x1a,y1a,z1a,x2a,y2a,z2a)**2
c sumn=sumn+(psi(r1a)*psi(r2b)+psi(r1b)*psi(r2a))**2
sumn=sumn+(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
$(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)*(psi(r1a)*psi(r2b)-
$psi(r1b)*psi(r2a))**2
sumkin=sumkin + (-1./2.)*psih2(u1a,u2a,u3a,v1a,v2a,v3a)*
$ d2psih2(u1a,u2a,u3a,v1a,v2a,v3a)*(1.+u1a**2)*(1.+u2a**2)*
$(1.+u3a**2)*(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)

c sumkin=sumkin+(-1./2.)*psi(r1a)*d2psih2(x1a,y1a,z1a)
c ajee=ajee+vee(r12)*(1.+u1a**2)*(1.+u2a**2)*(1.+u3a**2)*
c $(1.+v1a**2)*(1.+v2a**2)*(1.+v3a**2)*psi(r1a)**2*psi(r2a)**2
10 continue
anorm=vol2*sumn/float(nstep)
c sumpot=(1./anorm)*vol*sumpot/float(nstep)
sumpot= (1./anorm)*vol2*sumpot/float(nstep)
ajee=(vol2/anorm)*ajee/float(nstep)
c ekin= (1./anorm)*vol*sumkin/float(nstep)
ekin=(1./anorm)*vol2*sumkin/float(nstep)
Vnuclear=1./d
c print*,'nstep,potenergy=',nstep ,sumpot
c print*,'Kinetic energy =', ekin
c print*,'e=(defkin+defpot)/defnorm=',(defkin+defpot)/defnorm
c print*,'nstep,d,anorm=',nstep,d, anorm
c print*,'ajee,5./8.=', ajee ,5./8.
et=ekin + sumpot + vnuclear
print 110 ,d, sumpot, ekin,et
c print 115,d,et
50 continue
c print 100,nstep,defpot,defkin,defnorm,e
100 format(I6,4(3x,e10.3))
110 format('d,V,ek,et=',4(3x,e10.3))
115 format('d,etotal=',2(3x,e11.4))
stop
end

subroutine aleat(nstep,pwr,seed,pi,x)
dimension x(60000)
x0=seed
do 10 i=1,nstep
y=x0**pwr + sqrt(5.)
rnd=y-int(y)
x(i)=rnd
c print*,'rnd=',rnd
x0=rnd*pi*1.e4
10 continue
return
end