From e650ed1e1e908e51c78c1b047bec0da7c4fea366 Mon Sep 17 00:00:00 2001 From: sotech117 Date: Sat, 27 Apr 2024 04:25:23 -0400 Subject: testing --- hw8/10-14-old.png | Bin 0 -> 17063 bytes hw8/10-14.jl | 82 ++++++++++++++++++++++++++++ hw8/10-14.png | Bin 0 -> 16696 bytes hw8/10-14t.png | Bin 0 -> 16696 bytes hw8/10-17.jl | 81 ++++++++++++++++++++++++++++ hw8/10-3.jl | 92 ++++++++++++++++++++++++++++++++ hw8/BoundStates.jl | 73 +++++++++++++++++++++++++ hw8/RadialBoundStates.jl | 59 ++++++++++++++++++++ hw8/TimeDependentSchrodingerEquation.jl | 81 ++++++++++++++++++++++++++++ hw8/tightbinding.jl | 24 +++++++++ 10 files changed, 492 insertions(+) create mode 100644 hw8/10-14-old.png create mode 100644 hw8/10-14.jl create mode 100644 hw8/10-14.png create mode 100644 hw8/10-14t.png create mode 100644 hw8/10-17.jl create mode 100644 hw8/10-3.jl create mode 100644 hw8/BoundStates.jl create mode 100644 hw8/RadialBoundStates.jl create mode 100644 hw8/TimeDependentSchrodingerEquation.jl create mode 100644 hw8/tightbinding.jl (limited to 'hw8') diff --git a/hw8/10-14-old.png b/hw8/10-14-old.png new file mode 100644 index 0000000..c1d5618 Binary files /dev/null and b/hw8/10-14-old.png differ diff --git a/hw8/10-14.jl b/hw8/10-14.jl new file mode 100644 index 0000000..27d1b5c --- /dev/null +++ b/hw8/10-14.jl @@ -0,0 +1,82 @@ +#!/usr/bin/env julia + +using Plots +using LinearAlgebra +using DifferentialEquations + +# useful functions: + +function H(psi, dx) # action of Hamiltonian on wavefunction + Hpsi = zeros(ComplexF64, size(psi)) + # -(1/2) * laplacian(psi) (m = hbar = 1) + Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx) + + # periodic boundary conditions + #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) + #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) + return Hpsi +end + +derivative(psi, dx) = -1.0im * H(psi, dx) + +function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0) + Delta2 = Delta^2 + return exp.(-(x .- x0) .^ 2 / Delta2) .* exp.(1.0im * k * x) +end + +function normalization(psi, dx) # normalization of wavefunction + n2 = dot(psi, psi) * dx + return sqrt(n2) +end + +prob(psi) = real(psi .* conj(psi)) + +# The actual simulation +N = 400 # number of lattice points +L = 20.0 # x runs from 0 to L +dx = L / N + +x = range(0.0, L, N) |> collect # lattice along x-axis +#println(x) + +# initial wavefunction has position (x0), width (Delta), and momentum (k) +psi0 = initialWavefunction(x, 10.0, 1.0, 1000.0) + +# normalize wavefunction +n = normalization(psi0, dx) +psi0 = psi0 / n +println("norm = ", normalization(psi0, dx)) + +# plot initial wavefunction and probability density +plot(prob(psi0)) + +# integrate forward in time +tf = 1.0 +dt = 0.1 +tspan = (0.0, tf) + +function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm + psi = psi0 + for t in range(0, stop = tf, step = dt) + psiMid = psi + 0.5 * dt * derivative(psi, dx) + psi = psi + dt * derivative(psiMid, dx) + end + return psi +end + +# tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver +# problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE +# sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy + +# compare initial and final wavefunction probabilities +psi = timeEvolve(psi0, tf, dt) +# psi = sol[:, end] +times = sol.t + +# check that normalization hasn't deviated too far from 1.0 +println("norm = ", normalization(psi, dx)) + +plot([prob(psi0), prob(psi)]) +savefig("10-14.png") + + diff --git a/hw8/10-14.png b/hw8/10-14.png new file mode 100644 index 0000000..3e18e5b Binary files /dev/null and b/hw8/10-14.png differ diff --git a/hw8/10-14t.png b/hw8/10-14t.png new file mode 100644 index 0000000..3e18e5b Binary files /dev/null and b/hw8/10-14t.png differ diff --git a/hw8/10-17.jl b/hw8/10-17.jl new file mode 100644 index 0000000..6dfff2b --- /dev/null +++ b/hw8/10-17.jl @@ -0,0 +1,81 @@ +#!/usr/bin/env julia + +using Plots +using LinearAlgebra +using DifferentialEquations + +# useful functions: + +function H(psi, dx) # action of Hamiltonian on wavefunction + Hpsi = zeros(ComplexF64, size(psi)) + # -(1/2) * laplacian(psi) (m = hbar = 1) + Hpsi[2:end-1] = 0.5 * (2.0 * psi[2:end-1] - psi[3:end] - psi[1:end-2]) / (dx * dx) + + # periodic boundary conditions + #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) + #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) + return Hpsi +end + +derivative(psi, dx) = -1.0im * H(psi, dx) + +function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0) + Delta2 = Delta^2 + return exp.(-(x .- x0) .^ 2 / Delta2) .* exp.(1.0im * k * x) +end + +function normalization(psi, dx) # normalization of wavefunction + n2 = dot(psi, psi) * dx + return sqrt(n2) +end + +prob(psi) = real(psi .* conj(psi)) + +# The actual simulation +N = 400 # number of lattice points +L = 20.0 # x runs from 0 to L +dx = L / N + +x = range(0.0, L, N) |> collect # lattice along x-axis +#println(x) + +# initial wavefunction has position (x0), width (Delta), and momentum (k) +psi0 = initialWavefunction(x, 10.0, 0.05, 700.0) + +# normalize wavefunction +n = normalization(psi0, dx) +psi0 = psi0 / n +println("norm = ", normalization(psi0, dx)) + +# plot initial wavefunction and probability density +plot(prob(psi0)) + +# integrate forward in time +tf = 1.0 +dt = 5e-7 +tspan = (0.0, tf) + +function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm + psi = psi0 + for t in range(0, stop = tf, step = dt) + psiMid = psi + 0.5 * dt * derivative(psi, dx) + psi = psi + dt * derivative(psiMid, dx) + end + return psi +end + +tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver +problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE +sol = solve(problem, Tsit5(), reltol = 1e-12, abstol = 1e-12) # solve using Tsit5 algorithm to specified accuracy + +# compare initial and final wavefunction probabilities +#psi = timeEvolve(psi0, tf, dt) +psi = sol[:, end] +times = sol.t + +# check that normalization hasn't deviated too far from 1.0 +println("norm = ", normalization(psi, dx)) + +plot([prob(psi0), prob(psi)]) + + diff --git a/hw8/10-3.jl b/hw8/10-3.jl new file mode 100644 index 0000000..8616a69 --- /dev/null +++ b/hw8/10-3.jl @@ -0,0 +1,92 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of 1-D quantum problems""" + +using LinearAlgebra +using Plots + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(x, n) + """ The potential energy""" + #return 0.0 # particle in a box + #return 0.5 * x^2 # SHO with the spring constant k = 1 + #return -6.0 * x^2 + 8.0 * x^6 # potential with zero ground state energy + #return 0.1 * x^4 - 2.0 * x^2 + 0.0 * x # double-well potential + # return 8 * x^6 - 8 * x^4 - 4 * x^2 + 1 # another double-well potential + return abs(x)^n +end + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + + +function map_n_to_energies(n) + N = 1000 # number of lattice points + L = 20.0 # x runs from -L/2 to L/2 + dx = L / N + + D = zeros(N, N) # discrete laplacian operator + V = zeros(N, N) # potential + + for i in 1:N + D[i, i] = -2.0 + end + + for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 + end + + + for i in 1:N + x = (i + 0.5) * dx - 0.5 * L # coordinates of lattice points + V[i, i] = potential(x, n) + end + + H = -0.5 * dx^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + + + e, v = eigen(H) # diagonalize Hamiltonian + + println("\n n = ", n) + println("Ground state energy = ", e[1]) + println("1st excited state energy = ", e[2]) + println("2nd excited state energy = ", e[3]) + println("3rd excited state energy = ", e[4]) + println("4th excited state energy = ", e[5], "\n") + + return e +end + +n_max = 18 +n_to_e = [map_n_to_energies(n) for n in 1:n_max] + +# plot e[0] for all N +eList = zeros(0) +for i in 1:n_max + push!(eList, n_to_e[i][1]) +end + +plot(eList) + +# gs(x) = exp(-0.5 * x^2) # Gaussian that is exact ground state of SHO + + +# plot(potential) + + +# plot(v[:, 1]) +#plot(v[:,2]) +# plot(gs) + +#= +eList = zeros(0) +for i in 1:20 + push!(eList, e[i]) +end + + +bar(eList, orientation = :horizontal) +=# diff --git a/hw8/BoundStates.jl b/hw8/BoundStates.jl new file mode 100644 index 0000000..22e608c --- /dev/null +++ b/hw8/BoundStates.jl @@ -0,0 +1,73 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of 1-D quantum problems""" + +using LinearAlgebra +using Plots + +N = 1000 # number of lattice points +L = 20.0 # x runs from -L/2 to L/2 +dx = L / N + +D = zeros(N, N) # discrete laplacian operator +V = zeros(N, N) # potential + +for i in 1:N + D[i, i] = -2.0 +end + +for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 +end + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(x) + """ The potential energy""" + #return 0.0 # particle in a box + #return 0.5 * x^2 # SHO with the spring constant k = 1 + #return -6.0 * x^2 + 8.0 * x^6 # potential with zero ground state energy + #return 0.1 * x^4 - 2.0 * x^2 + 0.0 * x # double-well potential + return 8 * x^6 - 8 * x^4 - 4 * x^2 + 1 # another double-well potential +end + +for i in 1:N + x = (i + 0.5) * dx - 0.5 * L # coordinates of lattice points + V[i, i] = potential(x) +end + +H = -0.5 * dx^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + +e, v = eigen(H) # diagonalize Hamiltonian + +println("\nGround state energy = ", e[1]) +println("\n1st excited state energy = ", e[2]) +println("\n2nd excited state energy = ", e[3]) +println("\n3rd excited state energy = ", e[4]) +println("\n4th excited state energy = ", e[5]) + + +gs(x) = exp(-0.5 * x^2) # Gaussian that is exact ground state of SHO + + +plot(potential) + + +plot(v[:, 1]) +#plot(v[:,2]) +#plot(gs) + +#= +eList = zeros(0) +for i in 1:20 + push!(eList, e[i]) +end + + +bar(eList, orientation = :horizontal) +=# diff --git a/hw8/RadialBoundStates.jl b/hw8/RadialBoundStates.jl new file mode 100644 index 0000000..ed7bcc9 --- /dev/null +++ b/hw8/RadialBoundStates.jl @@ -0,0 +1,59 @@ +#!/usr/bin/env julia + +"""Find eigenstates and eigenenergies of central potential problems""" + +using LinearAlgebra +using Plots + +N = 5000 # number of lattice points +L = 20.0 # r runs from 0 to L +dr = L / N + +D = zeros(N, N) # discrete radial 2nd derivative operator +V = zeros(N, N) # potential + +for i in 1:N + D[i, i] = -2.0 +end + +for i in 1:N-1 + D[i, i+1] = 1.0 + D[i+1, i] = 1.0 +end + +#println("\nLattice Laplacian operator") +#println(D) + +function potential(r, ℓ = 0) + """ The potential energy""" + #return 0.5 * ell * (ℓ+1.0) * pow(r, -2.0) # V=0: Free particle in spherical coordinates + + return -1.0 / r + 0.5 * ℓ * (ℓ + 1.0) * r^(-2.0) # Hydrogen atom + + #return -r^(-1.1) + 0.5 * ℓ * (ℓ+1.0) * r^(-2.0) # modified Coulomb potential +end + +for i in 1:N + r = (i + 0.5) * dr # radial coordinates of lattice points + V[i, i] = potential(r, 0) +end + +H = -0.5 * dr^(-2.0) * D + V # Hamiltonian. Here m = hbar = 1 + +#println("\nMatrix elements of Hamiltonian = ") +#println(H) + +e, v = eigen(H) # diagonalize Hamiltonian + +println("\nGround state energy = ", e[1]) +println("\n1st excited state energy = ", e[2]) +println("\n2nd excited state energy = ", e[3]) +println("\n3rd excited state energy = ", e[4]) +println("\n4th excited state energy = ", e[5]) + + +plot(potential) + + +plot(v[:, 1]) +#plot(v[:,2]) diff --git a/hw8/TimeDependentSchrodingerEquation.jl b/hw8/TimeDependentSchrodingerEquation.jl new file mode 100644 index 0000000..dacbced --- /dev/null +++ b/hw8/TimeDependentSchrodingerEquation.jl @@ -0,0 +1,81 @@ +#!/usr/bin/env julia + +using Plots +using LinearAlgebra +using DifferentialEquations + +# useful functions: + +function H(psi, dx) # action of Hamiltonian on wavefunction + Hpsi = zeros(ComplexF64, size(psi)) + # -(1/2) * laplacian(psi) (m = hbar = 1) + Hpsi[2:end-1] = 0.5 * (2.0*psi[2:end-1] - psi[3:end] - psi[1:end-2])/(dx*dx) + + # periodic boundary conditions + #Hpsi[1] = 0.5 * (2.0*psi[1] - psi[2] - psi[end])/(dx*dx) + #Hpsi[end] = 0.5 * (2.0*psi[end] - psi[end-1] - psi[1])/(dx*dx) + return Hpsi +end + +derivative(psi, dx) = -1.0im * H(psi, dx) + +function initialWavefunction(x::Vector{Float64}, x0 = 10.0, Delta = 1.0, k = 4.0) + Delta2 = Delta^2 + return exp.(- (x .- x0).^2 / Delta2 ) .* exp.(1.0im * k * x) +end + +function normalization(psi, dx) # normalization of wavefunction + n2 = dot(psi, psi) * dx + return sqrt(n2) +end + +prob(psi) = real(psi .* conj(psi)) + +# The actual simulation +N = 400 # number of lattice points +L = 20.0 # x runs from 0 to L +dx = L / N + +x = range(0.0, L, N) |> collect # lattice along x-axis +#println(x) + +# initial wavefunction has position (x0), width (Delta), and momentum (k) +psi0 = initialWavefunction(x, 10.0, 1.0, 0.0) + +# normalize wavefunction +n = normalization(psi0, dx) +psi0 = psi0 / n +println("norm = ", normalization(psi0, dx)) + +# plot initial wavefunction and probability density +plot(prob(psi0)) + +# integrate forward in time +tf = 1.0 +dt = 0.0001 +tspan = (0.0, tf) + +function timeEvolve(psi0, tf, dt) # second order Runge-Kutta algorithm + psi = psi0 + for t in range(0, stop=tf, step=dt) + psiMid = psi + 0.5 * dt * derivative(psi, dx) + psi = psi + dt * derivative(psiMid, dx) + end + return psi +end + +tendency(psi, dx, t) = derivative(psi, dx) # use ODE solver +problem = ODEProblem(tendency, psi0, tspan, dx) # specify ODE +sol = solve(problem, Tsit5(), reltol=1e-12, abstol=1e-12) # solve using Tsit5 algorithm to specified accuracy + +# compare initial and final wavefunction probabilities +#psi = timeEvolve(psi0, tf, dt) +psi = sol[:, end] +times = sol.t + +# check that normalization hasn't deviated too far from 1.0 +println("norm = ", normalization(psi, dx)) + +plot([prob(psi0), prob(psi)]) + + diff --git a/hw8/tightbinding.jl b/hw8/tightbinding.jl new file mode 100644 index 0000000..0d1cc8c --- /dev/null +++ b/hw8/tightbinding.jl @@ -0,0 +1,24 @@ +using LinearAlgebra + +# Define the tight binding Hamiltonian +function tight_binding_hamiltonian(n_sites::Int64, t::Float64) + H = zeros(n_sites, n_sites) + for i in 1:(n_sites-1) + H[i, i+1] = t + H[i+1, i] = t + end + H +end + +# Find the eigenenergies of the Hamiltonian +function eigenenergies(H::Array{Float64, 2}) + eigvals(H) +end + +n_sites = 10 # number of sites in the chain +t = 1.0 # hopping parameter +H = tight_binding_hamiltonian(n_sites, t) +energies = eigenenergies(H) +println(energies) + + -- cgit v1.2.3-70-g09d2