I have published my code here: GitHub - VikramD2k/Hydrogen-Atom-Simulator: This computes and plots the probability density of any orbital of the hydrogen atom
and I am also pasting my code down below
##
using Symbolics
using Plots
using Statistics
##
function Legrande(l::Int, cosθ::Num)
if l == 0
legrande = 1
else
D_cst_l = Differential(cosθ)^l
legrande = (1 / (2^l * factorial(l))) * D_cst_l((cosθ^2 - 1)^l)
end
return expand_derivatives(legrande)
end
##
function AssoLegrande(l::Int, m::Int, cst::Float64)
@variables cosθ
if m == 0
asslegrande = (-1)^m * (1 - cosθ^2)^(m / 2) * Legrande(l, cosθ)
else
D_cst_m = Differential(cosθ)^m
asslegrande = (-1)^m * (1 - cosθ^2)^(m / 2) * D_cst_m(Legrande(l, cosθ))
end
asslegrande = expand_derivatives(asslegrande)
asslegrande = substitute(asslegrande, cosθ => cst)
return asslegrande
end
##
function Laguerre(q::Int, x::Num)
if q == 0
leg = 1
else
Dx_q = Differential(x)^q
leg = exp(x) * Dx_q(exp(-x) * x^q) / factorial(q)
end
return expand_derivatives(leg)
end
##
function AssoLaguerre(p::Int, q::Int, k::Float64)
@variables ρ
if p == 0
asso_lag = (-1)^p * Laguerre(p + q, ρ)
else
Dρ_p = Differential(ρ)^p
asso_lag = (-1)^p * Dρ_p(Laguerre(p + q, ρ))
end
asso_lag = expand_derivatives(asso_lag)
asso_lag = substitute(asso_lag, ρ => k)
return asso_lag
end
##
function Radial(n::Int, l::Int, r::Float64)
a = 0.529 * 10^-10 # Radius of hydrogen atom in meters
r_val = 2 * r / (n * a)
radial = sqrt((2 / (n * a))^3 * factorial(n - l - 1) / (factorial(n + l) * 2 * n)) * exp(-r / (n * a)) * (2 * r / (n * a))^l
radial = radial * AssoLaguerre(2 * l + 1, n - l - 1, r_val)
radial = expand_derivatives(radial)
return radial
end
##
function Angular(l::Int, m::Int, cst::Float64, phi::Float64)
i = Complex{Float64}(0, 1)
angular = sqrt((2 * l + 1) * factorial(l - m) / (4 * π * factorial(l + m))) * AssoLegrande(l, m, cst) * exp(i * m * phi)
return angular
end
##
function get_prob_density(n, l, m, rₘᵢₙ, rₘₐₓ, r_steps, θ_steps, φ_steps)
x = zeros(Total_steps)
y = zeros(Total_steps)
z = zeros(Total_steps)
ψ² = zeros(Total_steps) # Initializing placeholder for |ψ|²
counter = 1
for r in range(rₘᵢₙ, rₘₐₓ, length=r_steps)
for θ in range(0, π, length=θ_steps)
for ϕ in range(0, 2 * π, length=φ_steps)
ψ = eval(Radial(n, l, r) * Angular(l, m, cos(θ), ϕ)) # Here we got probability amplitudes at the point (r,θ,φ)
print("ψ is ", ψ, " Type of ψ is ", typeof(ψ), "\n")
ψ²[counter] = abs2(ψ)
x[counter] = r * sin(θ) * cos(ϕ)
y[counter] = r * sin(θ) * sin(ϕ)
z[counter] = r * cos(θ)
counter = counter + 1
end
end
end
ψ²_sum = sum(ψ²) # This should be close to 1
ψ² = ψ² / ψ²_sum # Now it is normalized
return ψ², x, y, z
end
##
function plot_probs(Total_steps::Int, ψ², x, y, z)
X, Y, Z = [], [], []
c = 0
for i in 1:Total_steps
if ψ²[i] < 1.1 && ψ²[i] > 0.9
c += 1
push!(X, x[i])
push!(Y, y[i])
push!(Z, z[i])
end
end
Plots.scatter3d(X, Y, Z, markersize=2, color=:red)
title = "n is $n, l is $l, m is $m"
title!(title)
xlabel!("x")
ylabel!("y")
zlabel!("z")
end
##
# With the above all the required functions are defined
a = 0.529 * 10^-10 # Radius of hydrogen atom in meters
rₘᵢₙ = 0.5 * a
rₘₐₓ = 50 * a
r_steps = 10
θ_steps = 10
φ_steps = 10
Total_steps = r_steps * θ_steps * φ_steps
n = 1 # principle quantum number
l = 0 # azimuthal quantum number
m = 0 # magnetic quantum number
##
prob_density, x, y, z =get_prob_density(n, l, m, rₘᵢₙ, rₘₐₓ, r_steps, θ_steps, φ_steps)
plot_probs(Total_steps, prob_density, x, y, z)
##
the error I am encountering is the following
ψ is 8.893942444146948e14 Type of ψ is Complex{Num}
1-element ExceptionStack:
LoadError: MethodError: no method matching Float64(::Num)
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
Float64(::IrrationalConstants.Sqrt2)
@ IrrationalConstants C:\Users\vikra\.julia\packages\IrrationalConstants\vp5v4\src\macro.jl:112
...
Stacktrace:
[1] convert(::Type{Float64}, x::Num)
@ Base .\number.jl:7
[2] setindex!(A::Vector{Float64}, x::Num, i1::Int64)
@ Base .\array.jl:1021
[3] get_prob_density(n::Int64, l::Int64, m::Int64, rₘᵢₙ::Float64, rₘₐₓ::Float64, r_steps::Int64, θ_steps::Int64, φ_steps::Int64)
@ Main c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:78
[4] top-level scope
@ c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:123
[5] eval
@ .\boot.jl:385 [inlined]
[6] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
@ Base .\loading.jl:2076
[7] include_string(m::Module, txt::String, fname::String)
@ Base .\loading.jl:2086
[8] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
@ Base .\essentials.jl:892
[9] invokelatest(::Any, ::Any, ::Vararg{Any})
@ Base .\essentials.jl:889
[10] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:271
[11] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:181
[12] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:276
[13] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:179
[14] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:38
[15] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:150
[16] with_logstate(f::Function, logstate::Any)
@ Base.CoreLogging .\logging.jl:515
[17] with_logger
@ .\logging.jl:627 [inlined]
[18] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:263
[19] #invokelatest#2
@ .\essentials.jl:892 [inlined]
[20] invokelatest(::Any)
@ Base .\essentials.jl:889
[21] (::VSCodeServer.var"#64#65")()
@ VSCodeServer c:\Users\vikra\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:34
in expression starting at c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:123
I think the abs2(ψ)
is unable to convert the Complex{Num}
type ψ to a float.
I have no clue how to fix the above issue. Also tell me if I my analysis is correct.
additional info:
I am new to julia and I have migrated this code from MATLAB. I have gotten great suggestions from @mikmoore related to this code in a previous topic (Getting Stack Overflow error when running Derivative(x)^n) . I have incorporated all the changes they suggested.