Table of Contents
Done
ToBeDone
Start with \[ E\left[ \psi \right] = \frac{\int dX\psi^{*}(X)H\psi(X)}{\int dX\psi^{*}(X)\psi(X)} = \frac{ \left\langle \psi \right| H \left| \psi \right\rangle}{\left\langle \psi | \psi \right\rangle}\]
Through the variational method, we can obtain approximately the correct energy of a Hamiltonian using any basis wave fucntions that satisfy the boundary conditions.
\[ \mathbf{H} = \left\langle \psi_m \right| H \left| \psi_n \right\rangle \\ \mathbf{S} = \left\langle \psi_m \mid \psi_n \right\rangle \\ \mathbf{C} = \ \]
The generalized eigenvalue eqaution \[ \beta_k \mathbf{A}\mathbf{x}_K = \alpha_k \mathbf{B}\mathbf{x}_K \\ k = 0,1,2,...\]
\(\mathbf{HC} = E\mathbf{SC}\) where \(\mathbf{C}\) is the matrix of basis vectors
To reach a Hamiltonian for which the eignstates can be found two approximation are made:
In \(\S4.3.1\), to find the wave function of a Helium atom, we start with the hamiltonian with Born-Oppenheimer Approximation, which is \[ \left[ -\frac{1}{2}\nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{ \left| \mathbf{r}_1 - \mathbf{r}_2 \right| } \right]\phi(\mathbf{r}_1)\phi(\mathbf{r}_2) = E\phi(\mathbf{r}_1)\phi(\mathbf{r}_2) \]
then multiply by \(\phi(\mathbf{r}_2)\) and integrate with respect to \(d\mathbf{r}_2\)
\[ \left[ -\frac{1}{2}\nabla_1^2 - \frac{2}{r_1} + \int d^3\mathbf{r}_2 \left|\phi(\mathbf{r}_2)\right|^2 \frac{1}{ \left| \mathbf{r}_1 - \mathbf{r}_2 \right|} \right] \phi(\mathbf{r}_1) = E'\phi(\mathbf{r}_1)\]
then \(\int d^3\mathbf{r}_2 \left|\phi(\mathbf{r}_2)\right|^2 \frac{1}{ \left| \mathbf{r}_1 - \mathbf{r}_2 \right|}\) is the Coulomb energy of particle 1 in the electric field generated by the charge density of particle 2
The wave function we have used is called uncorrelated because of the fact that the probablilty \(P(\mathbf{r}_1,\mathbf{r}_2)\) for finding an electron at \(\mathbf{r}_1\) and antoher one at \(\mathbf{r}_2\) is uncoorelated ... This does not mean that electrons do not feel each other: in the dertmination of the spatial function \(\phi\), the interaction term \(\frac{1}{\left|\mathbf{r}_1 - \mathbf{r}_2 \right|}\) have been taken into account. But this interaction has been an averaged way: it is not the actual position of \(\mathbf{r}_2\) that determines the wave funtion for electron 1, but the average charge distribution of electron 2.
Paremetrize \[\psi(\vec r) = \sum_1^4 C_p\chi_p (\vec r) \] leads directly to \[ \left[ -\frac{1}{2}\nabla_1^2 - \frac{2}{r_1} + \sum_{r,s=1}^4 C_r C_s \int d^3 r_2 \chi_r(\mathbf{r}_2) \chi_s(\mathbf{r}_2) \frac{1}{ \left| \mathbf{r}_1 - \mathbf{r}_2 \right|} \right] \sum_{q=1}^4 C_q \chi_q(\mathbf{r}_1) = E' \sum_{q=1}^4 C_q \chi_q(\mathbf{r}_1) \] multiply from the left by \(\chi_p(\mathbf{r}_1)\) and integrate by \(\mathbf{r}_1\) leads to \[ \sum_{pq} \left( h_pq + \sum_{rs} C_r C_s Q_{prqs} \right) C_q = E' \sum_{pq} S_{pq} C_q \]
eye(N)
eig( M )
inv( M )
M[1:,i]
quadgk( f, a, b)
#this is just here to make sure that I can plot every time I restart the ipython notebook
using PyPlot
#globals used in the Hydrogen 4Basis
H4B_ALPHA = [13.00773, 1.962079, 0.444529, 0.1219492]
pie3_2 = sqrt( pi^3 )
#compose two functions
function myCompose(f,g)
#given two functions of 1 variable each
return (x,y)->(f(x)*g(y))
end
function compose(f,g)
return (x)->(f(x)*g(x))
end
# used to clobber terms in a matrix with NaN
function squishMatrix(mat)
sz = size(mat)
for i = [1:sz[1]]
for j = [1:sz[2]]
if (mat[i,j] = NaN)
mat[i,j] = 0
end
end
end
return matm2_dif = squishMat(m2)
end
function printMatElem( f, n, m)
for i in [1:n]
for j in [1:m]
temp = f(i,j)
print("$(temp), ")
end
println()
end
end
function fillMat(n,m,f)
mat = Array(Float64,n,m)
for i in [1:n]
for j in [1:m]
mat[i,j] = f((i-1),(j-1))
end
end
return mat
end
# generates the nmth element of the S matrix from CP3.2
function S_inf(n,m)
if( (n+m)%2 == 0 )
return ( (2.0/(n+m+5)) - (4.0/(n+m+3)) + (2.0/(n+m+1)))
else
return 0
end
end
# generates the nmth element of the H matrix from CP3.2
function H_inf(n,m)
if((n+m)%2 == 0)
return ( -8.0 * ( (1-m-n-(2m*n)) / ((m+n+3)*(m+n+1)*(m+n-1)) ) )
else
return 0
end
end
function infWellEig(n)
return eig( inv(fillMat(n,n,S_inf))*fillMat(n,n,H_inf) )
end
function infWellStatesandEig(n)
mat =
(vals, vecs) = eig(mat)
return (mat, vals, vecs)
end
function infWellEnergies(n)
(energies,vecs) = infWellEig(n)
return sort(energies)
end
function infWellBasisFuncs(x,n)
return (x^(n+2))-(x^n)
end
#draws the wave functions
function infWellWavenVec(N; steps = 1000, x_start = -1.0, x_stop = 1.0, eig_vecs=nothing, n = 0)
#has the user passed any eigenvectors?
if (eig_vecs == nothing)
ev = eye(N)
else
ev = eig_vecs
end
x = linspace(x_start,x_stop,steps)
wave = 0*Array(Float64,steps)
for i = [1:steps]
for j = [1:N]
bfj_x = infWellBasisFuncs(x[i], (j-1))
if (bfj_x != NaN) #if the basis function j at x is NaN, skip it or it will create garbage
wave[i] += ev[j, (n+1)] * bfj_x #j-1 because the eigenfunctions starts an N=0,1,2,..
end
end
end
return (x,wave)
end
# infWellPlotn
function infWellPlotn(N; n = 0, steps = 1000, x_start = -1.0, x_stop = 1.0, eig_vecs=nothing)
(x,wave) = infWellWavenVec(N,n=n,steps=steps, x_start=x_start, x_stop=x_stop, eig_vecs=eig_vecs)
plot(x,wave)
end
# quantumWaveKVec
# description:
# generates a vector of the probability density based on the projection into an eigen base
# N - the number of basis functions, or the order of the approx Matrix
# domain
# basisFunction is a function of two values, position and wave number
# projectionEigenVector
function quantumWaveKVec(N, domain, basisFunctions, projectionEigenVector)
sz = length(domain)
wave = Array(Float64,sz)
bfk_x = Array(Float64,N)
for i = [1:sz]
for k = [1:N]
bfk_tmp = basisFunctions(domain[i], (k-1))
if (bfk_tmp != NaN)
wave[i] += ev[j, (n+1)] * bfj_x #j-1 because the eigenfunctions starts an N=0,1,2,..
end
end
end
return wave
end
# print the energies and Plot the prob densities for an inifite square well with N basis functions
function infSqrWellApproxPlotv1(N; steps = 1000, width = 2.0, center = 0.0, showWhich=-1)
#step on generate the eigen vectors and values
(vals, vecs) = infWellEig(N)
#find the order of sorting
sortOrder = sortperm(vals)
# create the domain and range for the wave functions
domain = linspace( ((-width/2.0)-center), ((width/2.0)+center), steps)
projectionEnergy = 0*Array(Float64,N)
projectionRange = 0*Array(Float64,steps,N) # an Nx1000 matrix for simpelist use
for i in [1:N]
projectionEnergy[i] = vals[sortOrder[i]]
projectionRange[1:,i] = quantumWaveKVec(N,domain, infWellBasisFuncs, vecs[1:,sortOrder[i]])
end
if (showWhich == -2)
for i in [1:N]
plot(domain,projectionRange[1:,i])
end
elseif (showWhich == -1) #show the first five (default)
for i in [1:5]
plot(domain,projectionRange[1:,i], label="Wave $(i-1) with Energy $(round(projectionEnergy[i],4))")
end
else
plot(domain,projectionRange[1:,showWhich])
end
if ( sign(projectionRange[ int(steps/2), 1]) == -1)
legend(loc="lower right")
else
legend(loc="upper right")
end
end
# infWellBasisPlot plots the first 5 orders of the basis wave functions
function infWellBasisPlot(;steps = 1000)
domain = linspace(-1.0,1.0,steps)
pRange = 0*Array(Float64,steps,5)
for i = [1:steps]
for j = [1:5]
pRange[i,j] = infWellBasisFuncs( domain[i], (j-1) )
end
end
for j = [1:5]
plot(domain, pRange[1:,j], label="Eqn order $(j)")
end
legend()
end
# calc the prob densities for an inifite square well with N basis functions
function infSqrWellApprox(N; steps = 1000, width = 2.0, center = 0.0)
#step on generate the eigen vectors and values
(vals, vecs) = infWellEig(N)
#find the order of sorting
sortOrder = sortperm(vals)
# create the domain and range for the wave functions
domain = linspace( ((-width/2.0)-center), ((width/2.0)+center), steps)
projectionEnergy = 0*Array(Float64,N)
projectionRange = 0*Array(Float64,steps,N) # an Nx1000 matrix for simpelist use
for i in [1:N]
projectionEnergy[i] = vals[sortOrder[i]]
projectionRange[1:,i] = quantumWaveKVec(N,domain, infWellBasisFuncs, vecs[1:,sortOrder[i]])
end
return (domain, projectionRange, projectionEnergy)
end
# plot the infSqrWellWaveApprox
function infSqrWellApproxPlot(N)
(d, pR, pE) = infSqrWellApprox(N)
for i in [1:5]
plot(domain,projectionRange[1:,i], label="Wave $(i-1) with Energy $(round(projectionEnergy[i],4))")
end
mid = int(size(pR)[1] /2.0)
if ( sign(pR[ int(steps/2), 1]) == -1)
legend(loc="lower right")
else
legend(loc="upper right")
end
end
function hydrogen4Basis(r,a)
return exp( -1 * H4B_ALPHA[a] * r^2)
end
# the matrix element generation functions have their indicies increased by one to work with fillMat funciton
function hydrogen4BOverlap(p,q)
return sqrt( (pi / ( H4B_ALPHA[p+1] + H4B_ALPHA[q+1] ))^3 )
end
function hydrogen4BKinetic(p,q)
return (3.0 * (H4B_ALPHA[q+1]*H4B_ALPHA[p+1]*pie3_2) / sqrt( (H4B_ALPHA[p+1]+H4B_ALPHA[q+1])^5 ))
end
function hydrogen4BCoulomb(p,q)
return (-2pi / (H4B_ALPHA[p+1]+H4B_ALPHA[q+1]))
end
function hydrogen4BHamiltonian(p,q)
return hydrogen4BKinetic(p,q) + hydrogen4BCoulomb(p,q)
end
#return the hamilton and the overlap matricies
function hydrogen_4BH(N)
return fillMat(N,N,hydrogen4BHamiltonian)
end
function hydrogen_4BS(N)
return fillMat(N,N,hydrogen4BOverlap)
end
function hydrogen_4BEigen(;N=4)
return eig( inv(hydrogen_4BS(N))*hydrogen_4BH(N) )
end
#return a domain in r and a range as a LinComb of gaussian waves
function hydrogen_ExactWaves(;r=6, giveMeEig = hydrogen_4BEigen )
itr = [1:4]
(h_vals,h_vecs) = giveMeEig()
h_r = linspace(0.0,r,1000)
#sort the eigen values and vectors
h_perm = sortperm(h_vals)
#create a matrix to store the values of the LinComb
h_wpd = zeros(1000,4)
h_rpd = zeros(4)
for i_r = [1:1000]
# for each eigen vector
for i_e in itr # iterate over the eigenvector values
for a in itr # iterate over basis functions
tmp = h_vecs[a, h_perm[i_e] ] * hydrogen4Basis( h_r[i_r] , a )
if ( (tmp == NaN) || (tmp < 1e-12) ) # squish values that create garbage
h_rpd[a] = 0.0
else
h_rpd[a] = tmp
end
end
h_wpd[i_r, i_e] = dot( h_rpd, h_vecs[1:,i_e] )
end
end
return (h_r, h_wpd, h_vals, h_vecs, h_perm)
end
# plot the "exact" solutions to the hydrogen atom wave functions using a Gaussian Type Orbital (GTO)
function hydrogen_PlotExactWaves(;r=6, giveMeEig= hydrogen_4BEigen )
(h_r, h_wpd, h_vals, h_vecs, h_perm) = hydrogen_ExactWaves(r=r, giveMeEig = giveMeEig)
#plot each exact wave funciton
for i = [1:4]
plot(h_r, h_wpd[1:,h_perm[i]], label="Wave $(i)")
end
legend()
end
function hydrogen_GEPEig()
#overlap matrix S
hS = hydrogen_4BS(4)
# diagonalize the S matrix
hEV = eigvecs(hS)
hs = squishMat(conj(transpose(hEV))*hS*hEV)
hV = hEV* hs^(-1/2)
hH = hydrogen_4BH(4)
# hEO is the value of Energies and Waves ("orbitals") of hydrogen atom using the method in sec 3.3
(h_vals,h_vecs) = eig( transpose(hV)*hH*hV )
return (h_vals, hV*h_vecs)
end
he_a = [0.298073, 1.242567, 5.782948, 38.474970]
function he_basis(r, a)
return exp( -1* he_a[a] * r^2 )
end
# the matrix element generation functions have their indicies increased by one to work with fillMat funciton
function he_s(p,q)
return sqrt( (pi / ( he_a[p] + he_a[q] ))^3 )
end
function he_k(p,q)
return (3.0 * (he_a[q+1]*he_a[p+1]*pie3_2) / sqrt( (he_a[p+1]+he_a[q+1])^5 ))
end
function he_c(p,q)
return (-4pi / (he_a[p+1]+he_a[q+1]))
end
function he_h(p,q)
return he_k(p,q) + he_c(p,q)
end
function he_q(p,r,q,s)
return ( 2*pi^(5/2) ) / ( (he_a[p] + he_a[q])*(he_a[r]+he_a[s])*sqrt( (he_a[p]+he_a[r]+he_a[q]+he_a[s]) ) )
end
#return the hamilton and the overlap matricies
function he_H()
return fillMat(4,4,he_h)
end
function he_S(N)
return fillMat(4,4,he_s)
end
function he_Q()
Q = zeros(4,4,4,4)
itr = [1:4]
for i = itr
for j = itr
for k = itr
for l = itr
Q[i,j,k,l] = he_q(i,j,k,l)
end
end
end
end
return Q
end
# a function that returns a 4x4 matrix from a 4x4x4x4 matrix
function he_subQ(Q,r,s)
tmp = zeros(4,4)
for i in [1:4]
tmp[:,i] = Q[:,r,:,s][:,:,i]
end
return tmp
end
function he_F(;C=eye(4))
Q = he_Q()
H = he_H()
Q_pq = zeros(4,4)
itr = [1:4]
for r = itr
for s = itr
Q_pq += (he_subQ(Q,r,s)*C[r,s])
end
end
return squishMat(H + Q_pq)
end
ffo = myCompose(x->x, x->x)
println(ffo(2,2))
println(ffo(2,4))
println(quadgk(sin,0,pi))
println(quadgk(compose(sin,cos),0,pi))
# test that the H_inf returns values for H matrix
printMatElem( H_inf, 5, 5)
# test that the S_inf returns values for H matrix
printMatElem( S_inf, 5, 5)
# test fillMat for the S matrix for the infitite sqaure well
S_temp = fillMat(5,5,S_inf)
# test fillMat for the H matrix for the infitite sqaure well
H_temp = fillMat(5,5,H_inf)
# do the eigen values match those in the book in table 3.1
eig(inv(S_temp)*H_temp )
# YES THEY DO! but they are out of order
(m1,v1,e1) = infWellStatesandEig(5)
infWellEnergies(8)
infWellEnergies(12)[1:5]
infWellEnergies(16)[1:5]
(p,q) = infWellEig(8)
# test ploting and
x_tmp = linspace(-1.0,1.0,100)
w_tmp = 0*Array(Float64,100)
ev_tmp = eye(5)
for i in [1:100]
for j in [1:5]
bf = infWellBasisFuncs(x_tmp[i],(j-1))
if( bf != NaN )
w_tmp[i] += (ev_tmp[j,1] * bf)
end
end
end
plot(x_tmp,w_tmp)
infWellWavenVec(5, n=0)
infWellPlotn(5)
infSqrWellApproxPlotv1(5)
infSqrWellApproxPlotv1(8)
infSqrWellApproxPlotv1(12)
infSqrWellApproxPlotv1(16)
infSqrWellApproxPlotv1(20)
20 is the last value before the the values start becomming imaginary.
infWellBasisPlot()
hfoo = hydrogen_4BEigen()
dot(hfoo[2][1:,1],hfoo[2][1:,2])
heig = hfoo[2]
hS = hydrogen_4BS(4)
hEV = eigvecs(hS)
hs = squishMat(conj(transpose(hEV))*hS*hEV)
hV = hEV* hs^(-1/2)
hId = squishMat( transpose(hV)*hS*hV )
hH = hydrogen_4BH(4)
# hEO_1 is the value of Energies and Waves ("orbitals") of hydrogen atom using the method in sec 3.3
hEO_1 = eig( transpose(hV)*hH*hV )
# this is my method for finding the Enegies and Orbitals by finding the eigens of the inverse of S multiplied by H
hEO_2 = eig( inv(hydrogen_4BS(4)) * hydrogen_4BH(4) )
hV*hEO_1[2]
inv(hV)*hEO_1[2]
#now find the difference between the two methods
hDiff = hEO_2[2]-hEO_1[2]
#reduce r range ( r =100 )
hydrogen_PlotExactWaves(r=100)
# r = 10
hydrogen_PlotExactWaves(r=10)
# changed hyrodgen_PlotExactWaves and hydrogen_ExactWaves to dynamically change r
hydrogen_PlotExactWaves(r=6)
hydrogen_PlotExactWaves(r=4)
hydrogen_PlotExactWaves(r=4,giveMeEig=hydrogen_GEPEig)
qaz = zeros(4,4,4,4)
qaz[1,1,1,1]
qaz[:,1,:,1]
qaz = he_Q()
qaz[:,:,1,1]
qaz[:,1,:,1]
qaz[:,:,1,2]
qaz[:,1,:,2]
qaz[:,1,:,2][:,:,1]
qaz[:,1,:,2][:,:,2]
qua = zeros(4,4)
for i in [1:4]
qua[:,i] = qaz[:,1,:,2][:,:,i]
end
qua
qaz[:,1,:,1]