Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Out-source n!/k! to Combinatorics.jl #7

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Out-source n!/k! to Combinatorics.jl #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
50387d5
Out-source n!/k! to Combinatorics.jl
jagot Jul 8, 2020
abf9519
Implemented radial moments for hydrogenic orbitals
jagot Dec 4, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ version = "0.1.0"
[deps]
AngularMomentumAlgebra = "0bb1b1ac-0216-11e9-361b-bb110c659ce1"
AtomicLevels = "10933b4c-d60f-11e8-1fc6-bd9035a249a1"
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
EnergyExpressions = "f4b57a2e-27c7-11e9-0cb0-edd185b706f6"
HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -15,6 +16,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
SpecialPolynomials = "a25cea48-d430-424a-8ee7-0d3ad3742e9e"

[compat]
Combinatorics = "1.0"
HypergeometricFunctions = "0.2, 0.3"
Polynomials = "1"
SpecialFunctions = "0.10"
Expand Down
2 changes: 1 addition & 1 deletion src/Hydrogen.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ module Hydrogen

using LinearAlgebra
using SparseArrays
using Combinatorics

using AtomicLevels
using AngularMomentumAlgebra
Expand All @@ -13,7 +14,6 @@ using SpecialPolynomials: Laguerre
using SpecialFunctions: gamma

include("sparse_builder.jl")
include("factorial_ratio.jl")

include("constants.jl")
include("energies.jl")
Expand Down
2 changes: 1 addition & 1 deletion src/dipoles.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ function radial_dipole_moment(a::Orbital, b::Orbital)
n == n′ && return 3/2*n*√(n^2-ℓ^2)

A = powneg1(n′-ℓ)/(4*factorial(2ℓ-1))
B = √(factorial_ratio(n+ℓ,n-ℓ-1)*factorial_ratio(n′+ℓ-1,n′-ℓ))
B = √(factorial(n+ℓ,n-ℓ-1)*factorial(n′+ℓ-1,n′-ℓ))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't the gamma function from SpecialFunctions be better here, for the high n/ell cases?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am actually unsure, when talking about ratios of factorials. Should be investigated. This is of course related to #2, however, a numerically stable evaluation of the formula will not help for arbitrarily large values of n/ell, anyway, since the dipole moments (in the length gauge) are analytically divergent. This is physically reasonable, since the matrix element is ~ <n' ell' | r | n ell>, which of course will increase, the larger the values of the quantum numbers are.

C = (4*n*n′)^(ℓ+1)*(n-n′)^(n+n′-2ℓ-2)/((n+n′)^(n+n′))

nᵣ = n-ℓ-1
Expand Down
25 changes: 0 additions & 25 deletions src/factorial_ratio.jl
View file
Open in desktop

This file was deleted.

37 changes: 36 additions & 1 deletion src/orbitals.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,4 +43,39 @@ function non_relativistic_orbital(n::Integer, ℓ::Integer, Z = 1)
r -> C(r) * exp(q*r)
end

export non_relativistic_orbital
"""
radial_moment(k, n, ℓ[, Z = 1])

Compute the radial moment ``⟨rᵏ⟩`` for a hydrogenic orbital ``(n,ℓ)``
in a Coulomb potential of charge ``Z``. The formulas are taken from Table 2⁵ of

- Condon, E. U., & Shortley, G. H. (1951). The theory of atomic
spectra. Cambridge: Cambridge University Press.

and are available for ``k ∈ -4:4``.
"""
function radial_moment(k::Integer, n::Integer, ℓ::Integer, Z = 1)
if k == 1
(3n^2 - ℓ*(ℓ+1))/2Z
elseif k == 2
n^2*(5n^2 + 1 - 3ℓ*(ℓ+1))/2Z^2
elseif k == 3
n^2*(35n^2*(n^2-1) - 30n^2*(ℓ+2)*(ℓ-1) + 3*(ℓ+2)*(ℓ+1)*ℓ*(ℓ-1))/8Z^3
elseif k == 4
n^4*(63n^4 - 35n^2*(2ℓ^2+2ℓ-3) + 5ℓ*(ℓ+1)*(3ℓ^2+3ℓ-10) + 12)/8Z^4
elseif k == -1
Z/n^2
elseif k == -2
Z^2/(n^3*(ℓ+1//2))
elseif k == -3
Z^3/(n^3*(ℓ+1)*(ℓ+1//2)*ℓ)
elseif k == -4
Z^4/2*(3n^2 - ℓ*(ℓ+1))/(n^5*(ℓ+3//2)*(ℓ+1)*(ℓ+1//2)*ℓ*(ℓ-1//2))
elseif k == 0
one(Z)
else
throw(ArgumentError("Radial moment not available for k = $(k)"))
end
end

export non_relativistic_orbital, radial_moment
8 changes: 0 additions & 8 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,14 +31,6 @@ function lyman_balmer²(o, n, ℓ)
end

@testset "Hydrogen.jl" begin
@testset "Factorial ratios" begin
for a = 1:6
for b = 1:6
@test Hydrogen.factorial_ratio(a,b) == factorial(a)/factorial(b)
end
end
end

@testset "Radial dipole moments" begin
# We cannot compare with Table 13 of Bethe and Salpeter 1977,
# since it provides so few figures and some of the values
Expand Down