Skip to content

Commit

Permalink
Change InfinitePotentialWell3D
Browse files Browse the repository at this point in the history
  • Loading branch information
@ohno
ohno committed Jul 2, 2024
1 parent bf57669 commit ab6dacc
Show file tree
Hide file tree
Showing 5 changed files with 454 additions and 99 deletions.
48 changes: 34 additions & 14 deletions docs/src/InfinitePotentialWell3D.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ Antique.InfinitePotentialWell3D

#### Potential
```@docs; canonical=false
Antique.V(::InfinitePotentialWell3D, ::Any, ::Any, ::Any)
Antique.V(::InfinitePotentialWell3D, ::Any...)
```

#### Eigenvalues
Expand All @@ -24,7 +24,7 @@ Antique.E(::InfinitePotentialWell3D)

#### Eigenfunctions
```@docs; canonical=false
Antique.ψ(::InfinitePotentialWell3D, ::Any, ::Any, ::Any)
Antique.ψ(::InfinitePotentialWell3D, ::Any...)
```

## Usage & Examples
Expand All @@ -33,31 +33,51 @@ Antique.ψ(::InfinitePotentialWell3D, ::Any, ::Any, ::Any)

```@example IPW3D
using Antique
IPW3D = InfinitePotentialWell3D(Lx=1.0, Ly=1.0, Lz=1.0, m=1.0, ℏ=1.0)
IPW3D = InfinitePotentialWell3D(L=[1.0,1.0,1.0], m=1.0, ℏ=1.0)
; #hide
```

Parameters:

```@repl IPW3D
IPW3D.Lx
IPW3D.Ly
IPW3D.Lz
IPW3D.L
IPW3D.m
IPW3D.ℏ
```

Eigenvalues:

```@repl IPW3D
E(IPW3D, nx=1, ny=1, nz=1)
E(IPW3D, nx=2, ny=1, nz=1)
E(IPW3D, nx=1, ny=2, nz=1)
E(IPW3D, nx=1, ny=1, nz=2)
E(IPW3D, nx=2, ny=2, nz=1)
E(IPW3D, nx=2, ny=1, nz=2)
E(IPW3D, nx=1, ny=2, nz=2)
E(IPW3D, nx=2, ny=2, nz=2)
E(IPW3D, n=[1,1,1])
E(IPW3D, n=[2,1,1])
E(IPW3D, n=[1,2,1])
E(IPW3D, n=[1,1,2])
E(IPW3D, n=[2,2,1])
E(IPW3D, n=[2,1,2])
E(IPW3D, n=[1,2,2])
E(IPW3D, n=[2,2,2])
```

Wave functions:

```@example IPW3D
using CairoMakie
# settings
f = Figure()
ax = Axis(f[1,1], xlabel=L"$x$", ylabel=L"$\psi(x,0.5,0.5)$")
# plot
w1 = lines!(ax, 0..1, x -> ψ(IPW3D, x,0.5,0.5, n=[1,1,1]))
w2 = lines!(ax, 0..1, x -> ψ(IPW3D, x,0.5,0.5, n=[2,1,1]))
w3 = lines!(ax, 0..1, x -> ψ(IPW3D, x,0.5,0.5, n=[1,2,1]))
w4 = lines!(ax, 0..1, x -> ψ(IPW3D, x,0.5,0.5, n=[3,1,1]))
w5 = lines!(ax, 0..1, x -> ψ(IPW3D, x,0.5,0.5, n=[4,1,1]))
# legend
axislegend(ax, [w1, w2, w3, w4, w5], [L"n=[1,1,1]", L"n=[2,1,1]", L"n=[1,2,1]", L"n=[3,1,1]", L"n=[4,1,1]"], position=:lb)
f
```

## Testing
Expand Down
67 changes: 35 additions & 32 deletions src/InfinitePotentialWell3D.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,43 +2,35 @@ export InfinitePotentialWell3D, V, E, ψ

# parameters
@kwdef struct InfinitePotentialWell3D
Lx = 1.0
Ly = 1.0
Lz = 1.0
L = [1.0, 1.0, 1.0]
m = 1.0
ℏ = 1.0
end

# potential
function V(model::InfinitePotentialWell3D, x,y,z)
Lx = model.Lx
Ly = model.Ly
Lz = model.Lz
return (0<x<Lx&&0<y<Ly&&0<z<Lz) ? 0 : Inf
function V(model::InfinitePotentialWell3D, x...)
L = model.L
return prod(@. 0<x<L) ? 0 : Inf
end

# eigenvalues
function E(model::InfinitePotentialWell3D; nx::Int=1, ny::Int=1, nz::Int=1)
if !(1 ≤ nx && 1 ≤ ny && 1 ≤ nz)
throw(DomainError("(nx,ny,nz) = ($nx,$ny,$nz)", "This function is defined for 1 ≤ nx, 1 ≤ ny and 1 ≤ nz."))
function E(model::InfinitePotentialWell3D; n::Vector{Int64}=[1,1,1])
if !(prod(1 .≤ n))
throw(DomainError("n = $n", "This function is defined for 1 .≤ n."))
end
Lx = model.Lx
Ly = model.Ly
Lz = model.Lz
L = model.L
m = model.m
ℏ = model.ℏ
return (ℏ^2*π^2) / (2*m) * (nx^2/Lx^2 + ny^2/Ly^2 + nz^2/Lz^2)
return sum(@. (ℏ^2*n^2*π^2) / (2*m*L^2))
end

# eigenfunctions
function ψ(model::InfinitePotentialWell3D, x,y,z; nx::Int=1, ny::Int=1, nz::Int=1)
if !(1 ≤ nx && 1 ≤ ny && 1 ≤ nz)
throw(DomainError("(nx,ny,nz) = ($nx,$ny,$nz)", "This function is defined for 1 ≤ nx, 1 ≤ ny and 1 ≤ nz."))
function ψ(model::InfinitePotentialWell3D, x...; n::Vector{Int64}=[1,1,1])
if !(prod(1 .≤ n))
throw(DomainError("n = $n", "This function is defined for 1 .≤ n."))
end
Lx = model.Lx
Ly = model.Ly
Lz = model.Lz
return (0<x<Lx&&0<y<Ly&&0<z<Lz) ? sqrt(8/(Lx*Ly*Lz))*sin(nx*π*x/Lx)*sin(ny*π*y/Ly)*sin(nz*π*z/Lz) : 0
L = model.L
return prod(@. 0<x<L) ? prod(@. sqrt(2/L) * sin(n*π*x/L)) : 0
end

# docstrings
Expand All @@ -57,18 +49,18 @@ and the Hamiltonian
Parameters are specified with the following struct:
```
IPW3D = InfinitePotentialWell3D(Lx=1.0, Ly=1.0, Lz=1.0, m=1.0, ℏ=1.0)
IPW3D = InfinitePotentialWell3D(L=[1.0,1.0,1.0], m=1.0, ℏ=1.0)
```
``L_x,L_y,L_z`` are the lengths of the box in ``x``,``y``,``z``-direction, ``m`` is the mass of particle and ``\hbar`` is the reduced Planck constant (Dirac's constant).
``L`` is a vector of the lengths of the box in ``x``,``y``,``z``-direction, ``m`` is the mass of particle and ``\hbar`` is the reduced Planck constant (Dirac's constant).
## References
- [D. A. McQuarrie, J. D. Simon, _Physical chemistry : a molecular approach_ (University Science Books, 1997)](https://uscibooks.aip.org/books/physical-chemistry-a-molecular-approach/) p.90, 3-9. The Problem of a Particle in a Three-Dimensional Box Is a Simple Extension of the One-Dimensional Case
""" InfinitePotentialWell3D

@doc raw"""
`V(model::InfinitePotentialWell3D, x,y,z)`
`V(model::InfinitePotentialWell3D, x...)`
```math
V(x,y,z) =
Expand All @@ -79,22 +71,33 @@ V(x,y,z) =
\end{array}
\right.
```
""" V(model::InfinitePotentialWell3D, x,y,z)
""" V(model::InfinitePotentialWell3D, x...)

@doc raw"""
`E(model::InfinitePotentialWell3D; nx::Int=1, ny::Int=1, nz::Int=1)`
`E(model::InfinitePotentialWell3D; n::Vector{Int64}=[1,1,1])`
```math
E_{n_x,n_y,n_z} = \frac{\hbar^2 \pi^2}{2 m} \left(\frac{n_x^2}{L_x^2} + \frac{n_y^2}{L_y^2} + \frac{n_z^2}{L_z^2}\right)
E_{n_x,n_y,n_z}
= \frac{\hbar^2 n_x^2 \pi^2}{2 m L_x^2}
+ \frac{\hbar^2 n_y^2 \pi^2}{2 m L_y^2}
+ \frac{\hbar^2 n_z^2 \pi^2}{2 m L_z^2}
```
""" E(model::InfinitePotentialWell3D; nx::Int=1, ny::Int=1, nz::Int=1)
""" E(model::InfinitePotentialWell3D; n::Vector{Int64}=[1,1,1])

@doc raw"""
`ψ(model::InfinitePotentialWell3D, x,y,z; nx::Int=1, ny::Int=1, nz::Int=1)`
`ψ(model::InfinitePotentialWell3D, x...; n::Vector{Int64}=[1,1,1])`
The wave functions can be expressed as products of wave functions in a one-dimensional box.
```math
\psi_{n_x,n_y,n_z}(x,y,z) = \psi_{n_x}(x)\psi_{n_y}(y)\psi_{n_z}(z) = \sqrt{\frac{8}{L_xL_yL_z}} \sin\left(\frac{n_x\pi x}{L_x}\right) \sin\left(\frac{n_y\pi y}{L_y}\right) \sin\left(\frac{n_z\pi z}{L_z}\right)
\begin{aligned}
\psi_{n_x,n_y,n_z}(x,y,z)
&= \psi_{n_x}(x)
\times \psi_{n_y}(y)
\times \psi_{n_z}(z) \\
&= \sqrt{\frac{2}{L_x}} \sin \frac{n_x \pi x}{L_x}
\times \sqrt{\frac{2}{L_y}} \sin \frac{n_y \pi y}{L_y}
\times \sqrt{\frac{2}{L_z}} \sin \frac{n_z \pi z}{L_z}
\end{aligned}
```
""" ψ(model::InfinitePotentialWell3D, x,y,z; nx::Int=1, ny::Int=1, nz::Int=1)
""" ψ(model::InfinitePotentialWell3D, x...; n::Vector{Int64}=[1,1,1])
116 changes: 66 additions & 50 deletions test/InfinitePotentialWell3D.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
IPW3D = InfinitePotentialWell3D(Lx=1.0, Ly=1.1, Lz=1.2, m=1.0, ℏ=1.0)


# <ψᵢ|ψⱼ> = ∫ψₙ*ψₙdx = δᵢⱼ
Expand All @@ -14,31 +13,40 @@ println(raw"""
```""")

@testset "<ψᵢ|ψⱼ> = ∫ψₙ*ψₙdx = δᵢⱼ" begin
println("ix | iy | iz | jx | jy | jz | analytical | numerical ")
println("-- | -- | -- | -- | -- | -- | ----------------- | ----------------- ")
# for ix in 1:2
# for iy in 1:2
# for iz in 1:2
# for jx in 1:2
# for jy in 1:2
# for jz in 1:2
# analytical = ((ix==jx && iy==jy && iz==jz) ? 1 : 0)
# numerical = quadgk(x ->
# quadgk(y ->
# quadgk(z ->
# conj(ψ(IPW3D, x,y,z, nx=ix, ny=iy, nz=iz)) * ψ(IPW3D, x,y,z, nx=jx, ny=jy, nz=jz)
# , 0.0, IPW3D.Lz, maxevals=10)[1]
# , 0.0, IPW3D.Ly, maxevals=10)[1]
# , 0.0, IPW3D.Lz, maxevals=10)[1]
# acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-2) : isapprox(analytical, numerical, rtol=1e-2)
# @printf("%2d | %2d | %2d | %2d | %2d | %2d | %17.12f | %17.12f %s\n", ix, iy, iz, jx, jy, jz, analytical, numerical, acceptance ? "✔" : "✗")
# @test acceptance
# end
# end
# end
# end
# end
# end
for IPW3D in [
InfinitePotentialWell3D(L=[1.0,1.0,1.0], m=1.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=1.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=2.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=1.0, ℏ=2.0)
]
@show IPW3D
println("ix | iy | iz | jx | jy | jz | analytical | numerical ")
println("-- | -- | -- | -- | -- | -- | ----------------- | ----------------- ")
for ix in 1:2
for iy in 1:2
for iz in 1:2
for jx in 1:2
for jy in 1:2
for jz in 1:2
analytical = ((ix==jx && iy==jy && iz==jz) ? 1 : 0)
numerical = quadgk(x ->
quadgk(y ->
quadgk(z ->
conj(ψ(IPW3D, x,y,z, n=[ix,iy,iz])) * ψ(IPW3D, x,y,z, n=[jx,jy,jz])
, 0.0, IPW3D.L[3], maxevals=10)[1]
, 0.0, IPW3D.L[2], maxevals=10)[1]
, 0.0, IPW3D.L[1], maxevals=10)[1]
acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5)
@printf("%2d | %2d | %2d | %2d | %2d | %2d | %17.12f | %17.12f %s\n", ix, iy, iz, jx, jy, jz, analytical, numerical, acceptance ? "✔" : "✗")
@test acceptance
end
end
end
end
end
end
println()
end
end

println("""```
Expand Down Expand Up @@ -113,33 +121,41 @@ are given by the sum of 2 Taylor series:
```
```""")

ψTψ(IPW3D, x,y,z; nx=0,ny=0,nz=0, Δx=0.01,Δy=0.01,Δz=0.01) = -IPW3D.ℏ^2/(2*IPW3D.m) * conj(ψ(IPW3D,x,y,z,nx=nx,ny=ny,nz=nz)) * (
( ψ(IPW3D,x+Δx,y,z,nx=nx,ny=ny,nz=nz) -2*ψ(IPW3D,x,y,z,nx=nx,ny=ny,nz=nz) + ψ(IPW3D,x-Δx,y,z,nx=nx,ny=ny,nz=nz) ) / Δx^2 +
( ψ(IPW3D,x,y+Δy,z,nx=nx,ny=ny,nz=nz) -2*ψ(IPW3D,x,y,z,nx=nx,ny=ny,nz=nz) + ψ(IPW3D,x,y-Δy,z,nx=nx,ny=ny,nz=nz) ) / Δy^2 +
( ψ(IPW3D,x,y,z+Δz,nx=nx,ny=ny,nz=nz) -2*ψ(IPW3D,x,y,z,nx=nx,ny=ny,nz=nz) + ψ(IPW3D,x,y,z-Δz,nx=nx,ny=ny,nz=nz) ) / Δz^2
ψTψ(IPW3D, x, y, z; n=[1,1,1], Δx=0.01, Δy=0.01, Δz=0.01) = -IPW3D.ℏ^2/(2*IPW3D.m) * conj(ψ(IPW3D,x,y,z,n=n)) * (
( ψ(IPW3D,x+Δx,y,z,n=n) -2*ψ(IPW3D,x,y,z,n=n) + ψ(IPW3D,x-Δx,y,z,n=n) ) / Δx^2 +
( ψ(IPW3D,x,y+Δy,z,n=n) -2*ψ(IPW3D,x,y,z,n=n) + ψ(IPW3D,x,y-Δy,z,n=n) ) / Δy^2 +
( ψ(IPW3D,x,y,z+Δz,n=n) -2*ψ(IPW3D,x,y,z,n=n) + ψ(IPW3D,x,y,z-Δz,n=n) ) / Δz^2
)

@testset "<ψₙ|H|ψₙ> = ∫ψₙ*Tψₙdx = Eₙ" begin
println(" nx | ny | nz | analytical | numerical ")
println(" -- | --- | --- | ----------------- | ----------------- ")
# for nx in [1,2]
# for ny in [1,2]
# for nz in [1,2]
# IPW3D = InfinitePotentialWell3D(Lx=1.0,Ly=2.0,Lz=3.0)
# analytical = E(IPW3D,nx=nx,ny=ny,nz=nz)
# numerical = quadgk(x ->
# quadgk(y ->
# quadgk(z ->
# ψTψ(IPW3D, x, y, z, nx=nx, ny=ny, nz=nz, Δx=IPW3D.Lx*0.0001, Δy=IPW3D.Ly*0.0001, Δz=IPW3D.Lz*0.0001)
# , 0, IPW3D.Lz, maxevals=5)[1]
# , 0, IPW3D.Ly, maxevals=5)[1]
# , 0, IPW3D.Lz, maxevals=5)[1]
# acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-1) : isapprox(analytical, numerical, rtol=1e-1)
# @test acceptance
# @printf(" %2d | %3d | %3d | %17.12f | %17.12f %s\n", nx, ny, nz, numerical, analytical, acceptance ? "✔" : "✗")
# end
# end
# end
for IPW3D in [
InfinitePotentialWell3D(L=[1.0,1.0,1.0], m=1.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=1.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=2.0, ℏ=1.0)
InfinitePotentialWell3D(L=[1.2,3.4,4.5], m=1.0, ℏ=2.0)
]
@show IPW3D
println(" nx | ny | nz | analytical | numerical ")
println(" -- | -- | -- | ----------------- | ----------------- ")
for nx in [1,2]
for ny in [1,2]
for nz in [1,2]
analytical = E(IPW3D,n=[nx,ny,nz])
numerical = quadgk(x ->
quadgk(y ->
quadgk(z ->
ψTψ(IPW3D, x, y, z, n=[nx,ny,nz], Δx=IPW3D.L[1]*0.0001, Δy=IPW3D.L[2]*0.0001, Δz=IPW3D.L[3]*0.0001)
, 0, IPW3D.L[3], maxevals=5)[1]
, 0, IPW3D.L[2], maxevals=5)[1]
, 0, IPW3D.L[1], maxevals=5)[1]
acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5)
@test acceptance
@printf(" %2d | %2d | %2d | %17.12f | %17.12f %s\n", nx, ny, nz, numerical, analytical, acceptance ? "✔" : "✗")
end
end
end
println()
end
end

println("""```""")
Loading

0 comments on commit ab6dacc

Please sign in to comment.