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}\]
In COmputational Physics, Thijssen begins by looking at using non-orthogonal functions to find the eigen energies and eigen wave functions of the potential well. Solving this problem is solving the generalized eigenvecoter problem (GEP), where \(\mathbf{H}\) is a matrix of the Hamiltonian values of non-orthonogal functions (NOFs), \(\mathbf{S}\) is a matrix of the crossterms of the NOFs, and \(\mathbf{C}\) is the matrix of eigenvector of the energy \(E\), or \[ \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} = \text{ Eigenvector of } \mathbf{H} \text{ also, the compents of the linear addition of the NOFs of the eigen wave function} \\ E = \text{ Eigen energy or the eigenvalue associated with } \mathbf{C} \]
The variational method as stated it would be used in the chapter is in fact not used. Instead, the main disccusion of the chapter is the solution to the general eigenvalue problem (GEP): \(\mathbf{HC} = E\mathbf{SC}\).
In the GEP
There
The generalized eigenvalue eqaution \[ \beta_k \mathbf{A}\mathbf{x}_K = \alpha_k \mathbf{B}\mathbf{x}_K \\ k = 0,1,2,...\]
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 \]
The
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; tol = 1e-16)
sz = size(mat)
for i = [1:sz[1]]
for j = [1:sz[2]]
if ((mat[i,j] = NaN) || (mat[i,j] < tol))
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]