H atom

The following example shows calculations using the H atom. In this case the full wavefunction is analytically known and we thus can calculate the complete Hamiltonian. This example does not really use DFT, but shows how to calculate radial integrals from known wavefunctions.

HWaveFunctions.Quanty

-- The calculations are done using a set of basis spin-orbitals
-- These are very often given by some radial function times a spherical harmonic
-- (or given by sum of these)
-- Here I show for the very simple example of an H atom how one can calculate the 
-- appropriate integrals
 
-- a gnuplot script to make the plots
gnuplotInput = [[
set autoscale 
set xtic auto 
set ytic auto  
set style line  1 lt 1 lw 1 lc rgb "#FF0000"
set style line  2 lt 1 lw 1 lc rgb "#0000FF"
set style line  3 lt 1 lw 1 lc rgb "#00C000"
set style line  4 lt 1 lw 1 lc rgb "#FF00FF"
set style line  5 lt 1 lw 1 lc rgb "#00FFFF"
set style line  6 lt 1 lw 1 lc rgb "#000000"
 
set ylabel "Rnl" font "Times,12"
set xlabel "r(A)" font "Times,12"
 
set out 'HWaveFunctions.ps'
set size 1.0, 1.0
set terminal postscript portrait enhanced color  "Times" 8
 
plot "Radial" using 1:2  title "R_{1s}(r)" with lines ls 1
plot "Radial" using 1:3  title "R_{2s}(r)" with lines ls 1,\
     "Radial" using 1:4  title "R_{2p}(r)" with lines ls 2
plot "Radial" using 1:5  title "R_{3s}(r)" with lines ls 1,\
     "Radial" using 1:6  title "R_{3p}(r)" with lines ls 2,\
     "Radial" using 1:7  title "R_{3d}(r)" with lines ls 3
plot "Radial" using 1:8  title "R_{4s}(r)" with lines ls 1,\
     "Radial" using 1:9  title "R_{4p}(r)" with lines ls 2,\
     "Radial" using 1:10 title "R_{4d}(r)" with lines ls 3,\
     "Radial" using 1:11 title "R_{4f}(r)" with lines ls 4
plot "Radial" using 1:12 title "R_{5s}(r)" with lines ls 1,\
     "Radial" using 1:13 title "R_{5p}(r)" with lines ls 2,\
     "Radial" using 1:14 title "R_{5d}(r)" with lines ls 3,\
     "Radial" using 1:15 title "R_{5f}(r)" with lines ls 4,\
     "Radial" using 1:16 title "R_{5g}(r)" with lines ls 5
 
set ylabel "<R|j|R>" font "Times,12"
set xlabel "q(A^{-1})" font "Times,12"
 
i=1
do for [n1=1:3] {
  do for [l1=0:n1-1] {
    do for [n2=1:3] {
      do for [l2=0:n2-1] {
        do for [k=abs(l1-l2):l1+l2:2] {
          i=i+1
          plot "Bessel" using 1:(column(i)) title "<R_{".n1." ".l1."}(r)|j[qr]|R_{".n2." ".l2."}(r)" with lines ls 6
        }
      }
    }
  }
}
 
]]
 
-- bohr radius in units of Angstrom
a0 = 0.52917721092;
 
-- Hydrogen wave functions with effective Z
function Rnl (n, l, r, Zeff)
  return 2^(l+1) * e^(-r*Zeff/(a0 * n)) * n^(-2-l) * r^l * (Zeff/a0)^(l+3/2) * math.sqrt(math.fact(n-l-1)/math.fact(n+l)) * math.LaguerreL(n-l-1, 2*l+1, 2*r*Zeff/(a0*n))
end
 
-- write the radial functions to a file
file = io.open("Radial", "w")
dr=0.01
Nr=5000
for ir = 0, Nr, 1 do
  r = ir*dr
  file:write(string.format("%15.8E ",r))
  for n=1, 5, 1 do
    for l=0, n-1, 1 do
      file:write(string.format("%15.8E ",Rnl(n,l,r,1)))
    end
  end
  file:write("\n")
end
file:close()
 
-- normalization
io.write("\n<Rnl|Rnl>\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
    sum=0
    for ir = 0, Nr, 1 do
      r = ir*dr
      sum = sum + dr * r^2 * Rnl(n,l,r,1)^2
    end
    io.write(string.format("%15.8E ",sum))
  end
  io.write("\n")
end
 
-- average distance
io.write("\n<Rnl|r|Rnl> [units of A]\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
    sum=0
    for ir = 0, Nr, 1 do
      r = ir*dr
      sum = sum + dr * r^2 * r * Rnl(n,l,r,1)^2
    end
    io.write(string.format("%15.8E ",sum))
  end
io.write("\n")
end
 
-- average distance
io.write("\n<Rnl|r|Rnl> [units of a0]\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
  sum=0
    for ir = 0, Nr, 1 do
    r = ir*dr
    sum = sum + dr * r^2 * r * Rnl(n,l,r,1)^2 / a0
  end
  io.write(string.format("%15.8E ",sum))
end
io.write("\n")
end
 
-- dipole moments
io.write("\n<Rn1l1|r|Rn2l2> [units of a0]\n");
dr=0.01
Nr=5000
io.write("        R(1s)           R(2s)           R(2p)           R(3s)           R(3p)            R(3d)\n")
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    io.write(string.format("R(%i %i) ",n1,l1))
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        sum=0
        for ir = 0, Nr, 1 do
          r = ir*dr
          sum = sum + dr * r^2 * Rnl(n1,l1,r,1) * r * Rnl(n2,l2,r,1) / a0
        end
        io.write(string.format("%15.8E ",sum))
      end
    end
    io.write("\n")
  end
end
 
-- Slater Integrals
io.write("\n<Rn1l1(r1) Rn2(r2)l2(r)| e^2/(|r1-r2|) |Rn1l1(r1) Rn2(r2)l2(r)> and <Rn1l1(r1) Rn2(r)l2(r2)| e^2/(|r1-r2|) |Rn1l1(r2) Rn2(r)l2(r1)> [units of Rydberg]\n");
dr=0.05
Nr=1000
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        io.write(string.format("R1(%i %i) R2(%i %i) ",n1,l1,n2,l2))
        for k=0, 2*math.min(l1,l2), 2 do
          io.write(string.format("F[%i]= ",k))
          sum=0
          for ir1 = 1, Nr, 1 do
            r1 = ir1*dr
            for ir2 = 1, Nr, 1 do
              r2 = ir2*dr
              sum = sum + a0 * dr^2* r1^2 * r2^2 * Rnl(n1,l1,r1,1)^2 * Rnl(n2,l2,r2,1)^2 * math.min(r1,r2)^k / math.max(r1,r2)^(k+1)
            end
          end
          io.write(string.format("%15.8E ",sum))
        end
        io.write("\n")
        io.write(string.format("                "))
        for k=math.abs(l1-l2), l1+l2, 2 do
          io.write(string.format("G[%i]= ",k))
          sum=0
          for ir1 = 1, Nr, 1 do
            r1 = ir1*dr
            for ir2 = 1, Nr, 1 do
              r2 = ir2*dr
                 sum = sum + a0 * dr^2* r1^2 * r2^2 * Rnl(n1,l1,r1,1) * Rnl(n1,l1,r2,1) * Rnl(n2,l2,r1,1) * Rnl(n2,l2,r2,1) * math.min(r1,r2)^k / math.max(r1,r2)^(k+1)
              end
            end
            io.write(string.format("%15.8E ",sum))
          end
          io.write("\n")
      end
    end
  end
end
 
-- Expectation value of bessel functions
-- q in units of inverse bohr radii
file = io.open("Bessel", "w")
file:write("<Rn1l1|j_k(qr)|Rn2l2> ");
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        for k=math.abs(l1-l2), l1+l2, 2 do
          file:write(string.format("%2i %2i %2i %2i %2i         ",n1,l1,n2,l2,k));
        end
      end
    end
  end
end
file:write("\n")
dr=0.05
Nr=1000
dq=0.1
Nq=100
for iq=0, Nq, 1 do
  q = iq*dq
  file:write(string.format("%15.8E ",q))
  for n1=1, 3, 1 do
    for l1=0, n1-1, 1 do
      for n2=1, 3, 1 do
        for l2=0, n2-1, 1 do
          for k=math.abs(l1-l2), l1+l2, 2 do
            sum=0
            for ir = 0, Nr, 1 do
              r = ir*dr
              sum = sum + dr * r^2 * Rnl(n1,l1,r,1) * math.SphericalBesselJ(k,q*r) * Rnl(n2,l2,r,1)
            end
            file:write(string.format("%15.8E ",sum))
          end
        end
      end
    end
  end
file:write("\n")
end
 
-- write the gnuplot script to a file
file = io.open("HWaveFunctions.gnuplot", "w")
file:write(gnuplotInput)
file:close()
 
-- call gnuplot to execute the script
os.execute("gnuplot HWaveFunctions.gnuplot ; ps2pdf HWaveFunctions.ps")

The text output of the example is:

HWaveFunctions.out

<Rnl|Rnl>
 9.99999991E-01 
 9.99999999E-01  1.00000000E+00 
 1.00000000E+00  1.00000000E+00  1.00000000E+00 
 1.00000000E+00  1.00000000E+00  1.00000000E+00  1.00000000E+00 
 9.99996618E-01  9.99997469E-01  9.99998634E-01  9.99999518E-01  9.99999915E-01 
 
<Rnl|r|Rnl> [units of A]
 7.93765819E-01 
 3.17506327E+00  2.64588605E+00 
 7.14389235E+00  6.61471514E+00  5.55636071E+00 
 1.27002531E+01  1.21710758E+01  1.11127214E+01  9.52518980E+00 
 1.98439701E+01  1.93148370E+01  1.82565430E+01  1.66690572E+01  1.45523689E+01 
 
<Rnl|r|Rnl> [units of a0]
 1.50000000E+00 
 6.00000000E+00  5.00000000E+00 
 1.35000000E+01  1.25000000E+01  1.05000000E+01 
 2.40000000E+01  2.30000000E+01  2.10000000E+01  1.80000000E+01 
 3.74996687E+01  3.64997521E+01  3.44998663E+01  3.14999528E+01  2.74999917E+01 
 
<Rn1l1|r|Rn2l2> [units of a0]
        R(1s)           R(2s)           R(2p)           R(3s)           R(3p)            R(3d)
R(1 0)  1.50000000E+00 -5.58701653E-01  1.29026620E+00 -2.43569644E-01  5.16689243E-01  3.85117423E-01 
R(2 0) -5.58701653E-01  6.00000000E+00 -5.19615242E+00 -1.85110879E+00  3.06481541E+00 -5.93938418E+00 
R(2 1)  1.29026620E+00 -5.19615242E+00  5.00000000E+00  9.38404238E-01 -1.76947200E+00  4.74799161E+00 
R(3 0) -2.43569644E-01 -1.85110879E+00  9.38404238E-01  1.35000000E+01 -1.27279221E+01  9.48683298E+00 
R(3 1)  5.16689243E-01  3.06481541E+00 -1.76947200E+00 -1.27279221E+01  1.25000000E+01 -1.00623059E+01 
R(3 2)  3.85117423E-01 -5.93938418E+00  4.74799161E+00  9.48683298E+00 -1.00623059E+01  1.05000000E+01 
 
<Rn1l1(r1) Rn2(r2)l2(r)| e^2/(|r1-r2|) |Rn1l1(r1) Rn2(r2)l2(r)> and <Rn1l1(r1) Rn2(r)l2(r2)| e^2/(|r1-r2|) |Rn1l1(r2) Rn2(r)l2(r1)> [units of Rydberg]
R1(1 0) R2(1 0) F[0]=  6.25361997E-01 
                G[0]=  6.25361997E-01 
R1(1 0) R2(2 0) F[0]=  2.09911370E-01 
                G[0]=  2.19839069E-02 
R1(1 0) R2(2 1) F[0]=  2.42809270E-01 
                G[1]=  5.12484496E-02 
R1(1 0) R2(3 0) F[0]=  9.94970763E-02 
                G[0]=  5.77815972E-03 
R1(1 0) R2(3 1) F[0]=  1.08828908E-01 
                G[1]=  1.36070257E-02 
R1(1 0) R2(3 2) F[0]=  1.11022542E-01 
                G[2]=  1.23687003E-03 
R1(2 0) R2(1 0) F[0]=  2.09911370E-01 
                G[0]=  2.19839069E-02 
R1(2 0) R2(2 0) F[0]=  1.50397569E-01 
                G[0]=  1.50397569E-01 
R1(2 0) R2(2 1) F[0]=  1.62113568E-01 
                G[1]=  8.79037050E-02 
R1(2 0) R2(3 0) F[0]=  8.41155296E-02 
                G[0]=  7.47769953E-03 
R1(2 0) R2(3 1) F[0]=  9.00128174E-02 
                G[1]=  8.49566045E-03 
R1(2 0) R2(3 2) F[0]=  1.03579160E-01 
                G[2]=  3.65275093E-02 
R1(2 1) R2(1 0) F[0]=  2.42809270E-01 
                G[1]=  5.12484496E-02 
R1(2 1) R2(2 0) F[0]=  1.62113568E-01 
                G[1]=  8.79037050E-02 
R1(2 1) R2(2 1) F[0]=  1.81647891E-01 F[2]=  8.79269507E-02 
                G[0]=  1.81647891E-01 G[2]=  8.79269507E-02 
R1(2 1) R2(3 0) F[0]=  8.67567735E-02 
                G[1]=  9.42538747E-03 
R1(2 1) R2(3 1) F[0]=  9.39457829E-02 F[2]=  2.58889132E-02 
                G[0]=  9.91050614E-03 G[2]=  1.13319343E-02 
R1(2 1) R2(3 2) F[0]=  1.06331648E-01 F[2]=  3.71302391E-02 
                G[1]=  3.73743204E-02 G[3]=  2.18070620E-02 
R1(3 0) R2(1 0) F[0]=  9.94970763E-02 
                G[0]=  5.77815972E-03 
R1(3 0) R2(2 0) F[0]=  8.41155296E-02 
                G[0]=  7.47769953E-03 
R1(3 0) R2(2 1) F[0]=  8.67567735E-02 
                G[1]=  9.42538747E-03 
R1(3 0) R2(3 0) F[0]=  6.64069248E-02 
                G[0]=  6.64069248E-02 
R1(3 0) R2(3 1) F[0]=  6.87938354E-02 
                G[1]=  4.23190720E-02 
R1(3 0) R2(3 2) F[0]=  7.31340175E-02 
                G[2]=  2.27882524E-02 
R1(3 1) R2(1 0) F[0]=  1.08828908E-01 
                G[1]=  1.36070257E-02 
R1(3 1) R2(2 0) F[0]=  9.00128174E-02 
                G[1]=  8.49566045E-03 
R1(3 1) R2(2 1) F[0]=  9.39457829E-02 F[2]=  2.58889132E-02 
                G[0]=  9.91050614E-03 G[2]=  1.13319343E-02 
R1(3 1) R2(3 0) F[0]=  6.87938354E-02 
                G[1]=  4.23190720E-02 
R1(3 1) R2(3 1) F[0]=  7.18684241E-02 F[2]=  3.59914255E-02 
                G[0]=  7.18684241E-02 G[2]=  3.59914255E-02 
R1(3 1) R2(3 2) F[0]=  7.69318422E-02 F[2]=  3.61349053E-02 
                G[1]=  3.41809433E-02 G[3]=  2.40553028E-02 
R1(3 2) R2(1 0) F[0]=  1.11022542E-01 
                G[2]=  1.23687003E-03 
R1(3 2) R2(2 0) F[0]=  1.03579160E-01 
                G[2]=  3.65275093E-02 
R1(3 2) R2(2 1) F[0]=  1.06331648E-01 F[2]=  3.71302391E-02 
                G[1]=  3.73743204E-02 G[3]=  2.18070620E-02 
R1(3 2) R2(3 0) F[0]=  7.31340175E-02 
                G[2]=  2.27882524E-02 
R1(3 2) R2(3 1) F[0]=  7.69318422E-02 F[2]=  3.61349053E-02 
                G[1]=  3.41809433E-02 G[3]=  2.40553028E-02 
R1(3 2) R2(3 2) F[0]=  8.60467604E-02 F[2]=  4.54247742E-02 F[4]=  2.96291766E-02 
                G[0]=  8.60467604E-02 G[2]=  4.54247742E-02 G[4]=  2.96291766E-02 

The created pdf is: HWaveFunctions.pdf