Differences

This shows you the differences between two versions of the page.

Link to this comparison view

documentation:tutorials:connection_to_dft:h_atom [2016/10/07 21:04] – created Maurits W. Haverkortdocumentation:tutorials:connection_to_dft:h_atom [2016/10/10 09:41] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +{{indexmenu_n>1}}
 +====== 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.
 +###
 +
 +===== Input =====
 +<code Quanty 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")
 +</code>
 +
 +===== Result =====
 +
 +The text output of the example is:
 +<file Quanty_Output 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 
 +</file>
 +
 +The created pdf is: {{ :documentation:tutorials:connection_to_dft:hwavefunctions.pdf |HWaveFunctions.pdf}}
 +
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » Connection to DFT and radial integrals » H atom

Print/export