dictAntoineCoefficients

Antoine coefficients [A,B,C,D] for temperature ranges below and above the melting points from dictMeltingPoints. These coefficients are used in the Antoine equation to calculate the saturated vapor pressure p (in Pa) at temperature `T (in K),

\[\mathrm{log_{e}}p=A+B/T+C\mathrm{log_{10}}T+D\cdot T/1000.\]

Currently, only the coefficients for the metalic elements are implemented - see C. B. Alcock, V. P. Itkin and M. K. Horrigan, Canadian Metallurgical Quarterly, 23, 309 (1984).

Note that the melting point is taken to be pressure independent.

julia> dictAntoineCoefficients
Dict{Int64, Tuple{Any, Any, Tuple{Any, Any, Any}}} with 102 entries:
  5  => (nothing, nothing, (empty, 2349.0, empty))
  56 => ([40.0897, -22312.0, -5.27062, 0.0], [20.7526, -18796.0, 0.0, 0.0], (298.0, 1000.0, 1200.0))
  35 => (nothing, nothing, (empty, 265.8, empty))
  55 => ([22.3736, -9208.04, 0.0, 0.0], [21.1164, -8818.9, 0.0, 0.0], (298.0, 301.7, 550.0))
  60 => ([32.2402, -39751.8, -2.19183, 0.0], [22.8364, -36436.1, 0.0, 0.0], (298.0, 1297.0, 2000.0))
  30 => ([25.5765, -15602.3, 0.0, 0.0], [23.9094, -14474.0, 0.0, 0.0], (298.0, 692.68, 750.0))
  32 => (nothing, nothing, (empty, 1211.4, empty))
  ⋮  => ⋮

Examples:

To calculate the saturated vapor pressure of Li (in Pa) at T=623 K we use the Antoine equation implemented in thefunction svp.

julia> svp(55, 400)
0.39420306801845933

julia> svp("Cs", 400)
0.39420306801845933

dictAtomicNumbers

julia> dictAtomicNumbers
Dict{String, Int64} with 102 entries:
  "Pd" => 46
  "Si" => 14
  "C"  => 6
  "P"  => 15
  "Nb" => 41
    ⋮  =>  ⋮

Examples:

julia> Z = get(dictAtomicNumbers, "Rb", nothing)
37

julia> listElement(Z; fmt=Info)
Element: rubidium
symbol: Rb
atomic number: Z = 37
atomic weight (relative atomic mass): 85.468

dictBigConversion

Dictionary for conversion from Int-based types to BigInt-based types

Example:

julia> dictBigConversion
Dict{DataType, DataType} with 12 entries:
  Vector{Rational{Int64}}         => Vector{Rational{BigInt}}
  Int64                           => BigInt
  Vector{Float64}                 => Vector{BigFloat}
  ComplexF64                      => Complex{BigFloat}
  Vector{Vector{ComplexF64}}      => Vector{Vector{Complex{BigFloat}}}
  Vector{Int64}                   => Vector{BigInt}
  Vector{Vector{Int64}}           => Vector{Vector{BigInt}}
  Vector{Vector{Float64}}         => Vector{Vector{BigFloat}}
  Vector{Vector{Rational{Int64}}} => Vector{Vector{Rational{Int64}}}
  Float64                         => BigFloat
  Rational{Int64}                 => Rational{BigInt}
  Vector{ComplexF64}              => Vector{Complex{BigFloat}}

dictElements

Standard atomic weights of the elements 2021 - see IUPAC Technical Report

julia> dictElements
Dict{Int64, Tuple{String, String, Any}} with 102 entries:
  5  => ("boron", "B", 10.81)
  56 => ("barium", "Ba", 137.33)
  35 => ("bromine", "Br", 79.904)
  55 => ("caesium", "Cs", 132.91)
  60 => ("neodymium", "Nd", 144.24)
  30 => ("zinc", "Zn", 65.38)
   ⋮ =>  ⋮

Examples:

julia> get(dictElements, 37, nothing)
("rubidium", "Rb", 85.468)

julia> listElement("Rb", fmt=Info)
Element: rubidium
  symbol: Rb
  atomic number: Z = 37
  atomic weight (relative atomic mass): 85.468

dictIsotopes

Sources: AME2020, LINDC(NDS)-0794 and INDC(NDS)-0833

julia> dictIsotopes
Dict{Tuple{Int64, Int64}, Tuple{String, String, Int64, Int64, Int64, Float64, Float64, Real, Int64, Float64, Float64, Any, Any}} with 341 entries:
  (71, 175) => ("¹⁷⁵Lu", "lutetium", 71, 175, 104, 5.37, 174.941, 7//2, 1, 1.0e…
  (40, 92)  => ("⁹²Zr", "zirconium", 40, 92, 52, 4.3057, 91.905, 0, 1, 1.0e100,…
  (48, 111) => ("¹¹¹Cd", "cadmium", 48, 111, 63, 4.5845, 110.904, 1//2, 1, 1.0e…
  (72, 176) => ("¹⁷⁶Hf", "hafnium", 72, 176, 104, 5.3286, 175.941, 0, 1, 1.0e10…
  (30, 68)  => ("⁶⁸Zn", "zinc", 30, 68, 38, 3.9658, 67.9248, 0, 1, 1.0e100, 0.0…
  (76, 184) => ("¹⁸⁴Os", "osmium", 76, 184, 108, 5.3823, 183.952, 0, 1, 5.6e13,…
  (54, 129) => ("¹²⁹Xe", "xenon", 54, 129, 75, 4.7775, 128.905, 1//2, 1, 1.0e10…
      ⋮     =>                                ⋮

Example:

julia> get(dictIsotopes, (37,87), nothing)
("⁸⁷Rb", "rubidium", 37, 87, 50, 4.1989, 86.90918053, 3//2, -1, 4.97e10, 2.75129, 0.1335, 27.83)

julia> listIsotope(37, 87, fmt=Info)
Isotope: rubidium-87
  symbol: ⁸⁷Rb
  element: rubidium
  atomic number: Z = 37
  atomic mass number: A = 87
  neutron number: N = 50
  rms nuclear charge radius: R = 4.1989 fm
  atomic mass: M = 86.90918053 amu
  nuclear spin: I = 3/2 ħ
  parity of nuclear state: π = odd
  nuclear magnetic dipole moment: μI = 2.75129 μN
  nuclear electric quadrupole moment: Q = 0.1335 barn
  relative abundance: RA = 27.83%
  lifetime: 4.97e10 years

dictMeltingPoints

Melting points of the elements at standard pressure (1atm) - see Wikipedia

julia> dictMeltingPoints
Dict{Int64, Union{Nothing, Real}} with 102 entries:
  5  => 2349   
  56 => 1000   
  35 => 265.8  
  55 => 301.7  
  60 => 1297   
  30 => 692.68 
  32 => 1211.4 
  6  => nothing
  67 => 1734   
  ⋮  => ⋮

Examples:

julia> get(dictMeltingPoints, 3, nothing)
453.65

julia> melting_point("Li")
453.65

Adams{T}

• G: (:Vector{Matrix{T}})
• σ: (:Vector{Matrix{T}})
• Minv: (:Vector{Matrix{T}})
• Z: (:Vector{Complex{T}})

Atom(Z, A, Q, Zc, element, isotope)

Type with fields:

• .Z: atomic number (::Int)
• .A: atomic mass number in amu (::Int)
• .Q: ionic charge in a.u. (::Int)
• .Zc: Rydberg charge in a.u. (::Int)
• .element: (::Element)
• .isotope: (::Isotope)

The type Atom is best created with the function castAtom.

Codata

Object to hold the natural constants from CODATA. It is best created with the function castCodata

The fields are:

• .∆νCs: Cs hyperfine transition frequency (::Value)
• .c: speed of light in vacuum (::Value)
• .h: Planck constant (::Value)
• .ħ: Planck constant - reduced (::Value)
• .e: elementary charge (::Value)
• .kB: Boltzmann constant (::Value)
• .NA: Avogadro constant (::Value)
• .Kcd: Luminous efficacy (::Value)
• .me: electron rest mass (::Value)
• .R∞: Rydberg constant (::Value)
• .Ry: Rydberg frequency (::Value)
• .Eh: Hartree a.u. (::Value)
• .α: fine-structure constant (::Value)
• .μ0: magnetic permitivity of vacuum (::Value)
• .ε0: electric permitivity of vacuum (::Value)
• .KJ: Josephson constant (::Value)
• .RK: Von Klitzing constant (::Value)
• .R: Molar gas constant (::Value)
• .u: unified atomic mass unit (::Value)
• .matE: unit conversion matrix (Matrix{Float64})

Example:

codata = castCodata(2018)
codata.μ0
  Value(1.2566370621250601e-6, "N A⁻²")

codata.μ0.val
  1.2566370621250601e-6

Def(T, atom, orbit, pot, scr, o1, o2, o3, pos, epn, k, am, matLD)

Type with fields:

• .T: gridType (::Type)
• .atom: atom object (::Atom)
• .orbit: orbit object (::Orbit)
• .codata: codata object (::Codata)
• .pot: tabulated potential function (::Vector{T})
• .scr: tabulated screening function (::Vector{T})
• .o1: vector of zero-filled matrices (::Vector{Matrix{T}})
• .o2: vector of zero-filled matrices (::Vector{Matrix{T}})
• .o3: vector of unit-filled matrices (::Vector{Matrix{T}})
• .pos: object containing Na, Nlctp, Nmin, Nuctp, Nb, N and nodes (::Pos)
• .epn: number of endpoints trapezoidal correction - must be odd (::Int)
• .k: Adams-Moulton order (::Int)
• .am: Adams-Moulton weight coefficients (::Vector{T})
• .matLD: Lagrangian differentiation matrix (::Matrix{T})

The object Def is best created with the function castDef.

Element(name, symbol, weight)

Type with fields:

• .name: name of element (::String)
• .symbol: symbol of element (::String)
• .weight: relative atomic mass - atomic weight (::Float64)

The type Element is best created with the function castElement.

Grid(ID, name, T, N, r, r′, h, r0, epn, epw, k)

Type with fields:

• .ID: grid identifer name (::Int)
• .name: grid identifer name (::String)
• .T: gridType (::Type)
• .N: number of grid points (::Int)
• .r: tabulated grid function (::Vector{T})
• .r′: tabulated derivative of grid function (::Vector{T})
• .h : grid step multiplyer (::T)
• .r0: grid scale factor (::T)
• .epn: number of endpoints used for trapezoidal endpoint correction (must be odd) (::Int)
• .epw: trapezoidal endpoint weights for n=1:epn (::Vector{Vector{T}})
• .k: Adams-Moulton order (::Int)

The object Grid is best created with the function castGrid.

Isotope(symbol, name, Z, A, N, R, M, I, π, T½, mdm, eqm, ra)

Type with fields:

• .symbol: symbol (::String)
• .name: name (::String)
• .Z: atomic number (::Int)
• .A: atomic mass number in amu (::Int)
• .N: neutron number (::Int)
• .R: rms charge radius in Fermi (::Float64)
• .M: atomic mass in amu (::Float64)
• .I: nuclear spin in units of ħ (::Rational{Int})
• .π: parity of nuclear state (::Int)
• .T½: lifetime in years (::Float64)
• .mdm: nuclear magnetic dipole moment (::Float64)
• .eqm: nuclear electric quadrupole moment (::Float64)
• .ra: relative abundance in % (::Float64)

The type Isotope is best created with the function castIsotope.

NamedValue(val::Value, name::String, comment::String)

Object to hold a Value together with its symbolic name and a short description

The fields are:

• .val: Value (::Value)
• .name: symbolic name (::String)
• .comment: description (::String)

Named Value object The object NamedValue is best created using castNamedValue.

Example:

f = Value(1,"Hz")
  Value(1, "Hz", "frequency")

f.name
  "frequency"

Orbit(name, n, n′, ℓ, mℓ)

Type for specification of atomic orbitals with fields:

• .name: name
• .n: principal quantum number
• .n′: radial quantum number (number of nodes in radial wavefunction)
• .ℓ: orbital angular momentum valence electron
• .mℓ: orbital angular momentum projection valence electron

The type Orbit is best created with the function castOrbit.

Pos(Na::Int, Nlctp::Int, Nmin::Int, Nuctp::Int, Nb::Int, N::Int, nodes::Int, cWKB::Float64)

Type with fields:

• .Na: grid index of last leading point (::Int)
• .Nlctp: grid index of lower classical turning point (::Int)
• .Nmin: grid index of (screened) potential minimum (::Int)
• .Nuctp: grid index of upper classical turning point (::Int)
• .Nb: grid index first trailing point (::Int)
• .N: grid index last point (::Int)
• .nodes: number of nodes (::Int)
• .cWKB: WKB threshold level determining Na and Nb (::Float64)

Mutable struct to hold special grid indices as well as the number of nodes; Pos is one of the fields of the Def object

Examples:

pos = Pos(1, 2, 3, 4, 5, 6, 7, 8)
pos.Nuctp
    4

pos.Nuctp = 8
pos
    Pos(1, 2, 3, 9, 5, 6, 7, 8)

Spinorbital

Type for specification of atomic Spinorbitals with fields:

• .name: name
• .n: principal quantum number
• .n′: radial quantum number (number of nodes in radial wavefunction)
• .ℓ: orbital angular momentum valence electron
• .ms: spin magnetic quantum number

The type Spinorbital is best created with the function castSpinorbital.

Term(name::String, n::Int, ℓ::Int, S::Real, L::Int, J::Real)

Type for specification of atomic fine-structure Terms with fields:

• name: name
• .n: principal quantum number
• .n′: radial quantum number (number of nodes in wavefunction)
• .ℓ: orbital angular momentum valence electron
• .S: total electron spin in units of ħ
• .L: total orbital angular momentum in units of ħ
• .J: total electronic angular momentum in units of ħ

The type Term is best created with the function createTerm.

Value(val::Real, unit::String)

Object to hold a real numerical value together with a unit specifier.

The fields are:

• .val: numerical value (::Real)
• .unit: unit specifier (::String)

Example:

f = Value(1,"Hz")
  Value(1, "Hz")

f.val
  1

f.unit
  "Hz"

CGC(j1::Real, m1::Real, j2::Real, m2::Real, J::Real, M::Real; msg=false)

Clebsch-Gordan coefficient (CGC). This is a vector-coupling coefficient in Dirac notation. The CGCs are zero unless $Δ(j_{1},j_{2},j_{3})>0$ (triangle inequality holds) and $M=m_{1}+m_{2}$. The relation to the Wigner 3j symbols is given by:

\[\langle j_{1}m_{1};j_{2}m_{2}|JM\rangle\equiv (-1)^{j_{1}-j_{2}+M}\sqrt{2J+1}\left(\begin{array}{ccc} j_{1} & j_{2} & J\\ m_{1} & m_{2} & -M \end{array}\right)\]

Examples:

julia> j1=3; m1=0; j2=4; m2=-1; J=5; M=-1;

julia> a = CGC(j1, m1, j2, m2, J, M)
-0.36364052611670255269921486774521555203216489725107181148303161368088211274565

julia> b = (-1)^(j1-j2+M) * sqrt(2J+1) * threeJsymbol(j1, m1, j2, m2, J, -M)
-0.36364052611670255269921486774521555203216489725107181148303161368088211274565

julia> o = CGC(j1, m1, j2, m2, J, M; msg=true); println(" = $o")
-√(361/2730) = -0.36364052611670255269921486774521555203216489725107181148303161368088211274565

INSCH(E::T, grid::Grid{T}, def::Def{T}, adams::Adams{T}) where T<:Real

OUTSCH(E::T, grid::Grid{T}, def::Def{T}, σ::Vector{Matrix{T}}})) where T<:Real

Solution of the Schrödinger for the first $k$ points on the grid, where $k$ is the Adams-Moulton order. The WKB solution for energy E is used when the WKB approximation is valid (for nonzero angular momentum at distances below the inner classical turning point - ictp)

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0)
Z = OUTSCH(Ecal, grid, def, adams.σ)
println("\nZ: standard Ansatz for wavefunction (n < Na=$(def.pos.Na)))")
    Orbital: 1s
        principal quantum number: n = 1
        radial quantum number: n′ = 0 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 0
    Grid created: exponential, Float64, Rmax = 63.0 a.u., Ntot = 100, h = 0.1, r0 = 0.00286033
    Def created for hydrogen 1s on exponential grid

    Z: standard Ansatz for wavefunction (n < Na=8))

Ecal, grid, def, adams = demo_hydrogen(n=10, ℓ=5)
Z = OUTSCH(Ecal, grid, def, adams.σ);
println("\nZ: WKB Ansatz for wavefunction (n < Na=$(def.pos.Na)))")
    Orbital: 10h
        principal quantum number: n = 10
        radial quantum number: n′ = 4 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 5
    Grid created: exponential, Float64, Rmax = 360.0 a.u., Ntot = 550, h = 0.0181818, r0 = 0.0163447
    Def created for hydrogen 10h on exponential grid

    Z: WKB Ansatz for wavefunction (n < Na=70))

plot_wavefunction(Z, 1:def.pos.Na, grid, def; reduced=true)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

RH1s(Z::U, r::T) where {U <: Real, T <:Real}

Analytic expression for the hydrogenic 1s radial wavefunction and its derivative in the format $Z = \tilde{R} + i \tilde{R}^′$, where

\[ \tilde{R}_{1s}(ρ) = Z^{3/2} 2 e^{-Zρ}\]

is the radial wavefunction and

\[ \tilde{R}^′_{1s}(ρ) = -Z^{5/2} 2 e^{-Zρ}\]

its derivative, with $ρ$ the radial distance to the nucleus in a.u..

Example:

atom = castAtom(Z=1, A=1, Q=0; msg=false);
orbit = castOrbit(n=1, ℓ=0; msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=false);
def = castDef(grid, atom, orbit, codata);

RH1s_example = [RH1s(atom.Z, grid.r[n]) for n=1:grid.N];

plot_wavefunction(RH1s_example, 1:grid.N, grid, def; reduced=false)

The plot is made using CairomMakie. NB.: plot_function is not included in the CamiXon package.

RH2p(Z::U, r::T) where {U <: Real, T <:Real}

Analytic expression for the hydrogenic 1s reduced radial wavefunction and its derivative in the format $Z = \tilde{R} + i \tilde{R}^′$, where

\[ \tilde{R}_{2p}(ρ)=\left(Z/2\right)^{3/2}\sqrt{1/3}(Zρ/2)2e^{-Zρ/2}\]

is the radial wavefunction and

\[ \tilde{R}_{2p}(ρ)=\left(Z/2\right)^{3/2}\sqrt{1/3}(1-Zρ/2)2e^{-Zρ/2}\]

its derivative, with $ρ$ the radial distance to the nucleus in a.u..

Example:

atom = castAtom(Z=1, A=1, Q=0; msg=false);
orbit = castOrbit(n=2, ℓ=1; msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=false);
def = castDef(grid, atom, orbit, codata);

RH2p_example = [RH2p(atom.Z, grid.r[n]) for n=1:grid.N];

plot_wavefunction(RH2p_example, 1:grid.N, grid, def; reduced=false)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

RH2s(Z::U, r::T) where {U <: Real, T <:Real}

Analytic expression for the hydrogenic 1s reduced radial wavefunction and its derivative in the format $Z = \tilde{R} + i \tilde{R}^′$, where

\[ \tilde{R}_{2s}(ρ)=\left(Z/2\right)^{3/2}(1-Zρ/2)2e^{-Zρ/2}\]

is the radial wavefunction and

\[ \tilde{R}_{2s}(ρ)=-\left(Z/2\right)^{5/2}(2-Zρ/2)2e^{-Zρ/2}\]

its derivative, with $ρ$ the radial distance to the nucleus in a.u..

Example:

atom = castAtom(Z=1, A=1, Q=0; msg=false);
orbit = castOrbit(n=2, ℓ=0; msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=false);
def = castDef(grid, atom, orbit, codata; msg=false);

RH2s_example = [RH2s(atom.Z, grid.r[n]) for n=1:grid.N];

plot_wavefunction(RH2s_example, 1:grid.N, grid, def; reduced=false)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

UF(k::Int, P::Vector{T}, grid::Grid{V}) where {T<:Real, V<:Real}

Potential for directe screening,

\[U_{F}^{k}(\rho) =\frac{1}{\rho^{k+1}}\int_{0}^{\rho}\varrho^{k} \left[\tilde{R}_{nl}(\varrho)\right]^{2} \varrho^{2}d\varrho+\rho^{k}\int_{\rho}^{\infty} \frac{1}{\varrho^{k+1}} \left[\tilde{R}_{nl}(\varrho)\right]^{2}\varrho^{2}d\varrho.\]

Example:

atom = castAtom(Z=2, A=4, Q=0; msg=true)
orbit = castOrbit(n=1, ℓ=0; msg=true)
scr = nothing
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=true)
def = castDef(grid, atom, orbit, codata; scr)
E = initE(def1)
adams = castAdams(E, grid, def)
E, def, adams, Z = adams_moulton_master(E, grid, def, adams; Δν=Value(1,"kHz"), imax=50, msg=false);

P1 = real(Z)
UF0P1 = UF(0, P1, grid);
plot_function(scrUF0P1, 1:grid.N, grid; title="He4(1s,1s):  direct screening potential")

The plot is made using CairomMakie. NB.: plot_function is not included in the CamiXon package.

UG(k::Int, P1::Vector{T}, P2::Vector{T}, grid::Grid{V}) where {T<:Real, V<:Real}

Potential for exchange screening,

\[U_{G}^{k}(\rho) =\frac{1}{\rho^{k+1}}\int_{0}^{\rho}\varrho^{k}\tilde{R}_{nl}(\varrho) \tilde{R}_{n^{\prime}l^{\prime}}(\varrho)\, \varrho^{2}d\varrho+\rho^{k}\int_{\rho}^{\infty} \frac{1}{\varrho^{k+1}}\tilde{R}_{nl}(\varrho) \tilde{R}_{n^{\prime}l^{\prime}}(\varrho)\,\varrho^{2}d\varrho.\]

Example:

atom = castAtom(Z=2, A=4, Q=0; msg=true)
orbit1 = castOrbit(n=1, ℓ=0; msg=true)
orbit2 = castOrbit(n=2, ℓ=0; msg=true)
scr = nothing
grid = autoGrid(atom, [orbit1,orbit2], Float64; Nboost=1, msg=true)
def1 = castDef(grid, atom, orbit1, codata; scr)
E = initE(def1)
adams = castAdams(E, grid, def1)
E, def, adams, Z1 = adams_moulton_master(E, grid, def1, adams; Δν=Value(1,"kHz"), imax=50, msg=false);

def2 = castDef(grid, atom, orbit2, codata; scr)
E = initE(def2)
adams = castAdams(E, grid, def2)
E, def, adams, Z2 = adams_moulton_master(E, grid, def2, adams; Δν=Value(1,"kHz"), imax=50, msg=false);

P1 = real(Z1);
P2 = real(Z2);

UG0 = UG(0, P1, P2, grid);
plot_function(UG0, 1:grid.N, grid; title="He4(1s,2s):  exchange screening potential")

The plot is made using CairomMakie. NB.: plot_function is not included in the CamiXon package.

a_direct(k::Int, l::Int, ml::Int, l′::Int, ml′::Int)

Coulomb angular integral - direct part:

\[a^{k}(lm_{l};l^{\prime}m_{l^{\prime}})=(-)^{m_{l}+m_{l^{\prime}}} (2l+1)(2l^{\prime}+1)\left(\begin{array}{ccc} l & k & l\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} l & k & l\\ -m_{l} & 0 & m_{l} \end{array}\right)\left(\begin{array}{ccc} l^{\prime} & k & l^{\prime}\\ 0 & 0 & 0 \end{array}\right)\left(\begin{array}{ccc} l^{\prime} & k & l^{\prime}\\ -m_{l^{\prime}} & 0 & m_{l^{\prime}} \end{array}\right)\]

Example:

a_direct(2,1,1,2,2)
    2//35

a_direct(6,3,2,3,-1)
    -250//20449

adams_moulton_inward(E::T, grid::Grid{T}, def::Def{T}, adams::Adams{T}) where T<:Real

adams_moulton_iterate(init::NTuple{4,T}, grid::Grid{T}, def::Def{T}, adams::Adams{T}; imax=25, Δν=Value(1,"kHz")) where T<:Real

Solves the Schrödinger equation for an atom defined by def for energy E on grid the grid with the Adams-Moulton method defined by adams; E is adjusted in an iteration procedure until convergence is reached within the convergence goal Δν is reached (limited to a maximum of imax iterations).

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0);
    Def created for hydrogen 1s on exponential grid of 100 points

E = 1.5Ecal;
msg1, adams, init, Z = adams_moulton_prepare(E, grid, def, adams);
println("Ecal = $Ecal; E = $(init[2]); $(def.pos.nodes) nodes")
    Ecal = -0.5; E = -0.75; 0 nodes

msg2, adams, init, Z = adams_moulton_iterate(init, grid, def, adams; Δν=Value(1,"MHz"), imax=25)
println("Ecal = $Ecal; E = $(init[2]); $(def.pos.nodes) nodes")
    Ecal = -0.5; E = -0.49999997841850014; 0 nodes

plot_wavefunction(Z, 1:def.pos.N, grid, def; reduced=false)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

adams_moulton_master(E, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=false)

Solves the Schrödinger equation for an atom defined by def for energy E on grid the grid with the Adams-Moulton method defined by adams.

Δν: convergence goal

imax: maximum number of iterations

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0);
    Def created for hydrogen 1s on exponential grid of 100 points

E = 1.5Ecal;
E, def, adams, Z = adams_moulton_master(E, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=true);
plot_wavefunction(Z, 1:def.pos.N, grid, def; reduced=false)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

adams_moulton_normalized(Z::Vector{Complex{T}}, ΔQ::T, grid::Grid{T}, def::Def{T}) where T<:Real

adams_moulton_outward(def::Def{T}, adams::Adams{T}) where T<:Real

adams_moulton_patch(Z::Vector{Complex{T}}, def::Def{T}, adams::Adams{T}) where T<:Real

Correct first 2k points of Z.

adams_moulton_prepare(E::T, grid::Grid{T}, def::Def{T}, adams::Adams{T}) where T<:Real

Solves the Schrödinger equation for an atom defined by def for energy E on grid the grid with the Adams-Moulton method defined by adams. E is adjusted until the wavefunction has the correct number of n′ nodes.

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0);
    Def created for hydrogen 1s on exponential grid of 100 points

E = 1.5Ecal
msg, adams, init, Z = adams_moulton_prepare(E, grid, def, adams);
    Ecal = -0.5; E = -0.75; 0 nodes

plot_wavefunction(Z, 1:def.pos.N, grid, def; reduced=false)

The plot is made using CairomMakie. Note the discontinuity in the derivative. NB.: plot_wavefunction is not included in the CamiXon package.

adams_moulton_solve(E::T, grid::Grid{T}, def::Def{T}, adams::Adams) where T<:Real

Numerical solution of the 1D Schrödinger equation for the radial motion of a valence electron of energy E. Output: the improved Adams object, the energy convergence ΔE, and Z, where P = real(Z) is the reduced radial wavefunction and Q = imag(Z) its derivative.

Example:

atom = castAtom(Z=1, A=1, Q=0, msg=true)
orbit = castOrbit(n=1, ℓ=0)
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=true)
def = castDef(grid, atom, orbit, codata)
E = Ecal = convert(grid.T, bohrformula(atom.Z, orbit.n))
adams = castAdams(E, grid, def);

adams, ΔE, Z = adams_moulton_solve(E, grid, def, adams)
plot_wavefunction(Z, 1:grid.N, grid, def; reduced=true)

The plot is made using CairomMakie. NB.: plot_wavefunction is not part of the CamiXon package.

autoGrid(atom, orbit,  T; Nboost=1, epn=5, k=7, msg=true, p=0)
autoGrid(atom, orbits, T; Nboost=1, epn=5, k=7, msg=true, p=0)
autoGrid(atom, orbit,  T; Nboost=1, epn=5, k=7, msg=true, coords=[])
autoGrid(atom, orbits, T; Nboost=1, epn=5, k=7, msg=true, coords=[])

Automatic setting of grid parameters for a given orbit Orbit or an array of orbits - orbits = [orbit1, orbit2, ⋯]. Important cases:

• p == 0 (exponential radial grid)
• p == 1 (linear radial grid)
• p > 1 (quasi-exponential radial grid)
• coords=[] (free polynomial grid based on the coords)
• Nboost (multiplier to boost numerical precision)
• epn (endpoint number: odd number to be used for trapezoidal integration with endpoint correction)
• k (Adams-Moulton order to be used for k+1-point Adams-Moulton integration)

Example:

codata = castCodata(2018)
atom = castAtom(;Z=1, A=1, Q=0, msg=false)
orbit = castOrbit(n=75, ℓ=0, msg=false)
grid = autoGrid(atom, orbit, Float64);
    Grid created: exponential, Float64, Rmax = 16935.0 a.u., Ntot = 3800, h = 0.00263158, r0 = 0.768883

plot_gridfunction(grid, 1:grid.N; title="")

The plot is made using CairomMakie. NB.: plot_gridfunction is not part of the CamiXon package.

autoNtot(orbit::Orbit, Nboost=1)

Total number of gridpoints (rule of thumb value)

\[ N_{tot} = (70 + 50 * n) * N_{boost},\]

where $n$ is the principal quantum number and Nboost a multiplier to boost numerical precision

Example:

orbit = castOrbit(n=1, ℓ=0)
autoNtot(orbit)
 Orbit created: 1s - (n = 1, n′ = 0, ℓ = 0)

 100

autoPrecision(Rmax::T, orbit::Orbit) where T<:Real

Floating point precision (rule of thumb value)

Example:

atom = castAtom(Z=1)
orbit = castOrbit(n=1,ℓ=0)
Rmax = autoRmax(atom, orbit)
o = autoPrecision(Rmax, orbit); println("precision = $o")
    Element created: H, hydrogen, Z=1, weight=1.008
    Isotope created: ¹H, hydrogen, Z=1, A=1, N=0, R=0.8783, M=1.007825032, I=1/2⁺, μI=2.792847351, Q=0.0, RA=99.9855%, (stable)
    Atom created: hydrogen, neutral atom, ¹H, Z=1, A=1, Q=0, Zc=1
    Orbital: 1s
        principal quantum number: n = 1
        radial quantum number: n′ = 0 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 0
    precision = Float64

autoRmax(atom::Atom, orbit::Orbit)

Largest relevant radial distance in a.u. (rule of thumb value)

\[ R_{max} = (2n^2 + 20n + 62)/Zc,\]

where $n$ is the principal quantum number and $Z_c$ the Rydberg charge

Example:

codata = castCodata(2018)
atom = castAtom(Z=1, A=1, Q=0)
orbit = castOrbit(n=1, ℓ=0)
rmax = autoRmax(atom::Atom, orbit::Orbit); println("rmax = $(rmax) a.u.")
    Element created: H, hydrogen, Z=1, weight=1.008
    Isotope created: ¹H, hydrogen, Z=1, A=1, N=0, R=0.8783, M=1.007825032, I=1/2⁺, μI=2.792847351, Q=0.0, RA=99.9855%, (stable)
    Atom created: hydrogen, neutral atom, ¹H, Z=1, A=1, Q=0, Zc=1
    Orbital: 1s
        principal quantum number: n = 1
        radial quantum number: n′ = 0 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 0
    rmax = 63.0 a.u.

autoSteps(ID::Int, Ntot::Int, Rmax::T; p=5) where T<:Real
autoSteps(ID::Int, Ntot::Int, Rmax::T; coords=[0,1]) where T<:Real

Step size parameter (h) and range parameter (r0) (rule of thumb values).

Example:

(h, r0) = autoSteps(1, 100, 100)
    (0.1, 0.004540199100968777)

b_exchange(k::Int, l::Int, ml::Int, l′::Int, ml′::Int)

Coulomb angular integral - exchange part:

\[b^{k}(lm_{l};l^{\prime}m_{l^{\prime}})=(2l+1)(2l^{\prime}+1) \left(\begin{array}{ccc} l & k & l^{\prime}\\ 0 & 0 & 0 \end{array}\right)^{2}\left(\begin{array}{ccc} l & k & l^{\prime}\\ -m_{l} & (m_{l}-m_{l^{\prime}}) & m_{l^{\prime}} \end{array}\right)^{2}\]

Example:

b_exchange(1,1,1,2,2)
    2//5

b_exchange(6,3,2,3,-1)
    1050//20449

bigconvert(o)

Convert IntBased types to BigIntBased types for n > nc in accordance with dictBigConversion.

Example:

julia> o = [[1//1, 1//2, 1//3],[1//1, 1//2, 1//3]]
2-element Vector{Vector{Rational{Int64}}}:
 [1//1, 1//2, 1//3]
 [1//1, 1//2, 1//3]

julia> bigconvert(o)
2-element Vector{Vector{Rational{Int64}}}:
 [1//1, 1//2, 1//3]
 [1//1, 1//2, 1//3]

bohrformula(Z::Int, n::Int)

Hydrogenic energy (in Hartree a.u.) for atom with atomic number Z and principal quantum number n.

\[ E_n = - \frac{Z^2}{2n^2}\]

Example:

Z = 2
n = 4
bohrformula(Z,n)
 -0.125

calibrationReport(E, Ecal, codata::Codata; unitIn="Hartree", msg=true)

Comparison of energy E with calibration value Ecal

default input: Hartree

Example:

codata = castCodata(2018)
calibrationReport(1.1, 1.0, codata; unitIn="Hartree")
  calibration report (Float64):
  Ecal = 1.0 Hartree
  E = 1.1 Hartree
  absolute accuracy: ΔE = 0.1 Hartree (657.968 THz)
  relative accuracy: ΔE/E = 0.0909091

castAdams(E::T, grid::Grid{T}, def::Def{T}) where T<:Real

castAtom(;Z=1, A=1, Q=0, msg=false)
castAtom(elt::String; A=1, Q=0, msg=false)

Create Atom with fields:

• .Z: atomic number (::Int)
• .A: atomic mass number in amu (::Int)
• .Q: ionic charge in a.u. (::Int)
• .Zc: Rydberg charge in a.u. (::Int)
• .element: (::Element)
• .isotope: (::Isotope)

Examples:

castAtom("Rb"; A=87, Q=0) == castAtom(Z=37, A=87, Q=0)
  true

castAtom(Z=1, A=3, Q=0)
  Atom(1, 3, 0, 1, Element("hydrogen", "H", 1.008), Isotope("³T", "tritium",
  1, 3, 2, 1.7591, 3.016049281, 1//2, 1, 12.33, 2.97896246, 0.0, nothing))

atom = castAtom(Z=1, A=3, Q=0, msg=true);
  Element created: H, hydrogen, Z=1, weight=1.008
  Isotope created: ³T, tritium, Z=1, A=3, N=2, R=1.7591, M=3.016049281, I=1/2⁺, μI=2.97896246, Q=0.0, RA=trace, (radioactive)
  Atom created: tritium, neutral atom, ³T, Z=1, A=3, Q=0, Zc=1

atom
  Atom(1, 3, 0, 1, Element("hydrogen", "H", 1.008), Isotope("³T", "tritium",
  1, 3, 2, 1.7591, 3.016049281, 1//2, 1, 12.33, 2.97896246, 0.0, nothing))

atom.isotope.T½
  12.33

castCodata(year::Int)

Method to create the Codata object

Example:

codata = castCodata(2018)
strValue.([codata.∆νCs,codata.c,codata.h])
 3-element Vector{String}:
  "9192631770 Hz"
  "299792458 m s⁻¹"
  "6.62607e-34 J Hz⁻¹"

castDef(grid::Grid{T}, atom::Atom, orbit::Orbit, codata::Codata[; scr=nothing[, msg=true]]) where T <: Real

Create the Def object starting from the Grid object and the atomic properties of the objects Atom and Orbit. Optional: scr (supply screening array)

Example:

codata = castCodata(2018)
atom = castAtom(Z=1, A=1, Q=0)
orbit = castOrbit(n=7, ℓ=2)
grid = autoGrid(atom, orbit, Float64)
def = castDef(grid, atom, orbit, codata);
    Element created: H, hydrogen, Z=1, weight=1.008
    Isotope created: ¹H, hydrogen, Z=1, A=1, N=0, R=0.8783, M=1.007825032, I=1/2⁺, μI=2.792847351, Q=0.0, RA=99.9855%, (stable)
    Atom created: hydrogen, neutral atom, ¹H, Z=1, A=1, Q=0, Zc=1
    Orbital: 7d
        principal quantum number: n = 7
        radial quantum number: n′ = 4 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 2
    Grid created: exponential, Float64, Rmax = 207.0 a.u., Ntot = 400, h = 0.025, r0 = 0.00939821
    Def created for hydrogen 7d on exponential grid of 400 points

castElement(;Z=1, msg=true)
castElement(elt::String; msg=true)

Create Atom with fields

• .name: name of element
• .symbol: symbol of element
• .weight: relative atomic mass (atomic weight)

Example:

castElement("Rb"; msg=false) == castElement(Z=37, msg=false)
  true

element = castElement(;Z=1, msg=true)
element
  Element created: H, hydrogen, Z=1, weight=1.008

  Element("hydrogen", "H", 1.008)

castGrid(ID::Int, N::Int, T::Type; h=1, r0=1,  p=5, coords=[0,1], epn=5, k=7, msg=false)

Method to create the Grid object

ID = 1: exponential grid, ID = 2: quasi-exponential grid, ID = 3: linear grid ID = 4: polynomial grid

Examples:

h = 0.1
r0 = 1.0
grid = castGrid(1, 4, Float64; h, r0, msg=true)
grid.r
  create exponential Grid: Float64, Rmax = 0.491825 a.u., Ntot = 4, h = 0.1, r0 = 1.0
  [1.0e-100, 0.10517091807564771, 0.22140275816016985, 0.3498588075760032]

grid = castGrid(2, 4, Float64; p = 4, h, r0, msg=true))
grid.r
  create quasi-exponential Grid: Float64, Rmax = 0.491733 a.u., Ntot = 4, p = 4, h = 0.1, r0 = 1.0
  [1.0e-100, 0.10517083333333321, 0.22140000000000004, 0.3498375]

grid = castGrid(3, 4, Float64; coords=[0, 1, 1/2, 1/6, 1/24], h, r0, msg=true)
grid.r
  create polynomial Grid: Float64, Rmax = 0.491733 a.u., Ntot = 4, coords = [0.0, 1.0, 0.5, 0.166666, 0.0416666], h = 0.1, r0 = 1.0
  [1.0e-100, 0.10517083333333334, 0.2214, 0.3498375000000001]

grid = castGrid(4, 4, Float64; h, r0, msg=true)
grid.r
  create linear Grid: Float64, Rmax = 0.4 a.u., Ntot = 4, p = 1, h = 0.1, r0 = 1.0
  [1.0e-100, 0.1, 0.2, 0.3]

grid.r′
  [0.1, 0.1, 0.1, 0.1]

castIsotope(;Z=1, A=1, msg=false)
castIsotope(elt::String; A=1, msg=false)

Create Isotope with fields

• .symbol: symbol (::String)
• .name: symbol (::String)
• .Z: atomic number (::Int)
• .A: atomic mass number in amu (::Int)
• .N: neutron number (::Int)
• .R: rms charge radius in Fermi (::Float64)
• .M: atomic mass in amu (::Float64)
• .I: nuclear spin in units of ħ (::Rational{Int})
• .π: parity of nuclear state (::Int)
• .ra: relative abundance in % (::Float64)
• .mdm: nuclear magnetic dipole moment (::Float64)
• .eqm: nuclear electric quadrupole moment (::Float64)
• .T½: lifetime in years (::Float64)

Examples:

castIsotope("Rb"; A=87) == castIsotope(Z=37, A=87)
  true

isotope = castIsotope(Z=1, A=3)
  Isotope("³T", "tritium", 1, 3, 2, 1.7591, 3.016049281, 1//2, 1, 12.33, 2.97896246, 0, nothing)

isotope.T½
  12.33

castIsotope(Z=1, A=3, msg=true);
  Isotope created: tritium-3
      symbol: ³T
      element: tritium
      atomic number: Z = 1
      atomic mass number: A = 3
      neutron number: N = 2
      rms nuclear charge radius: R = 1.7591 fm
      atomic mass: M = 3.016049281 amu
      nuclear spin: I = 1/2 ħ
      parity of nuclear state: π = ⁺
      nuclear magnetic dipole moment: μI = 2.97896246μN
      nuclear electric quadrupole moment: Q = 0.0barn
      relative abundance: RA = trace
      lifetime: 12.33 years

castNamedValue(val::Value; name=" ", comment=" ")

Method to create a NamedValue object

Example

v = Value(1.602176634e-19, "C")
nv = castNamedValue(v; name="e")
nv.name * " = " * strValue2(nv.val)
  "e = 1.60218e-19 C"

castOrbit(;n=1, ℓ=0, mℓ=0, msg=true)

Create Orbit with fields:

• .name: name
• .n: principal quantum number
• .n′: radial quantum number (number of nodes in radial wavefunction)
• .ℓ: orbital angular momentum valence electron
• .mℓ: orbital angular momentum projection valence electron

Examples:

castOrbit(n=1, ℓ=0)
 Orbit created: 1s (n = 1, n′ = 0, ℓ = 0)
 Orbit("1s", 1, 0, 0)

castSpinorbital(o::Orbit; up=true, msg=true)

Specify Spinorbital with fields:

• .name: name
• .n: principal quantum number
• .n′: radial quantum number (number of nodes in radial wavefunction)
• .ℓ: orbital angular momentum valence electron
• .ms: spin magnetic quantum number

Examples:

s1s = castOrbit(1,0)
castSpinorbital(s1s; up=true)
  Spinorbital created: 1s↑ (n = 1, n′ = 0, ℓ = 0, ms = 1//2)
  Spinorbital("1s↑", 1, 0, 0, 1//2)

conditionalType(n::T, nc::T [; msg=true]) where T<:Integer

Convert type T to BigInt for n > nc.

Example:

julia> conditionalType(46, 46)
Int64

julia> conditionalType(47, 46)
BigInt

convertUnit(val, codata; unitIn="Hartree", unitOut="xHz")

Unit conversion between μHz,⋯ EHz, Hartree, Rydberg, Joule, and eV

default input: Hartree

default output: xHz ∈ {μHz, mHz, Hz, kHz, MHz, GHz, THz, PHz, EHz}

Example:

codata = castCodata(2018)
convertUnit(1, codata; unitIn="Hz", unitOut="Joule")
  6.62607015e-34

convertUnit(1, codata; unitIn="Hartree", unitOut="Hz")
  Value(6.57968392050182e15, "Hz")

f = convertUnit(1, codata) # default input (Hartree) and output (xHz)
strf = strValue(f)
  "6.57968 PHz"

count_nodes(Z::Vector{Complex{T}}, def::Def{T}) where T<:Real

Number of nodes (excluding the origin) of the reduced radial wavefunction χ(r) = real(Z).

Example:

atom = castAtom(Z=1, A=1, Q=0, msg=false);
orbit = castOrbit(n=3, ℓ=2, msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, epn=5, k=7, msg=false);
def = castDef(grid.T, atom, orbit, codata);
    Def created for hydrogen 3d on exponential grid of 200 points

E = convert(setT, bohrformula(atom.Z, orbit.n));
adams = castAdams(E, grid, def);
E, def, adams, Z = adams_moulton_master(E, codata, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=false);

o = count_nodes(Z, def); println("node count: $o nodes")
    node count: 0 nodes

createTerm(n::Int; ℓ=0, S=1//2, L=0, J=1//2, msg=true)

Specify Term in the Term notatation with fields:

• .n: principal quantum number
• .n′: radial quantum number (number of nodes - autogenerated)
• .ℓ: orbital angular momentum valence electron
• .S: total electron spin
• .L: total orbital angular momentum
• .J: total electronic angular momentum

Examples:

term_H1I = createTerm(1; ℓ=0, S=1//2, L=0, J=1//2)
 Term created: 1s ²S₁⸝₂, n = 1, n′ = 0, ℓ = 0, S = 1//2, L = 0, J = 1//2
 Term("1s ²S₁⸝₂", 1, 0, 0, 1//2, 0, 1//2)

create_adams_moulton_weights(k::Int [; rationalize=false [, devisor=false [, T=Int]]])

$k^{th}$-order Adams-Moulton weights vector,

\[y[n+1] = y[n] + \frac{1}{D}\sum_{j=0}^{k}a^k[j]f[n+1-k+j]\]

The weights are stored in the vector $a^k \equiv[a_k^k/D,⋯\ a_0^k/D]$ under the convention $a^k[j] \equiv a_{k-j}^k/D$, where $a_j^k$ are the Adams-Moulton weight coefficients and $D$ the corresponding Adams-Moulton divisor. By default the output is in Float64, optionally the output is rational, with or without specification of the gcd devisor.

Example:

[create_adams_moulton_weights(k; rationalize=true, devisor=true, T=Int) for k=1:8]
8-element Vector{Tuple{Int64, Int64, Vector{Int64}}}:
 (1, 2, [1, 1])
 (2, 12, [-1, 8, 5])
 (3, 24, [1, -5, 19, 9])
 (4, 720, [-19, 106, -264, 646, 251])
 (5, 1440, [27, -173, 482, -798, 1427, 475])
 (6, 60480, [-863, 6312, -20211, 37504, -46461, 65112, 19087])
 (7, 120960, [1375, -11351, 41499, -88547, 123133, -121797, 139849, 36799])
 (8, 3628800, [-33953, 312874, -1291214, 3146338, -5033120, 5595358, -4604594, 4467094, 1070017])

create_lagrange_differentiation_matrix(k::Int)

Lagrange differentiation matrix, $m[i,j]=s_{k-j}^k(i)$, for $k^{th}$-order lagrangian differentiation,

\[\frac{dy}{dx}[i]= \sum_{j=0}^{k}m[i,j]y[j],\]

Example:

k = 3
create_lagrange_differentiation_matrix(k)
 4×4 Matrix{Rational{Int64}}:
  -11//6   3//1  -3//2   1//3
   -1//3  -1//2   1//1  -1//6
    1//6  -1//1   1//2   1//3
   -1//3   3//2  -3//1  11//6

demo_hydrogen(; n=3, ℓ=2, codata=castCodata(2018))

Solves Schrödinger equation for hydrogen atom with principal quantum number n and rotational quantum number ℓ.

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0);
    Def created for hydrogen 1s on exponential grid of 100 points

E = 1.5Ecal
E, def, adams, Z = adams_moulton_master(E, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=true);

plot_wavefunction(Z, 1:def.pos.N, grid, def; reduced=false)

The plot is made using CairomMakie. Note the discontinuity in the derivative. NB.: plot_wavefunction is not included in the CamiXon package.

edges(px [, Δx[, x0]])

Heatmap range transformation from pixel coordinates to physical coordinates, with pixelsize Δx and offset x0, both in physical units.

Examples:

px = 1:5
Δx = 2.5
x0 = 2.5
edges(px)
 [0.5, 1.5, 2.5, 3.5, 4.5]

edges(px, Δx)
 [1.25, 3.75, 6.25, 8.75, 11.25]

edges(px, Δx, x0)
 [-1.25, 1.25, 3.75, 6.25, 8.75]

fdiff_adams_bashford_expansion_coeffs(k [; T=Int])

$(k+1)$-point Adams-Bashford expansion coefficients $B_p$.

\[-\frac{∇}{(1-∇)ln(1-∇)}=\sum_{p=0}^{\infty}B_p∇^p=1+\ \frac{1}{2}∇+\ \frac{5}{12}∇^2+\ ⋯.\]

The weights are stored in forward order: $[B_0^k,⋯\ B_k^k]$ - order of use in summation.

Examples:

k = 5
o = fdiff_adams_bashford_expansion_coeffs(k); println(o)
 Rational{Int64}[1//1, 1//2, 5//12, 3//8, 251//720, 95//288]

fdiff_adams_moulton_expansion_coeff(n::T; msg=true) where {T<:Integer}
fdiff_adams_moulton_expansion_coeffs(n::T; msg=true) where {T<:Integer}

Finite difference expansion coefficient vector $β ≡ [β_0(x),\ ⋯,\ β_p(x)]$ defining $k^{th}$-order Adams-Moulton expansion,

\[-\frac{∇}{ln(1-∇)} = \sum_{p=0}^{\infty}β_p∇^p = 1 - \frac{1}{2}∇ - \frac{1}{12}∇^2 - \frac{1}{24}∇^3 +⋯.\]

Examples:

julia> k = 5;
julia> β = fdiff_adams_moulton_expansion_coeffs(k); println(β)
Rational{Int64}[1//1, -1//2, -1//12, -1//24, -19//720, -3//160]

julia> D = denominator(gcd(β))
1440

julia> o = convert(Vector{Int},(β .* D)); println(o)
[1440, -720, -120, -60, -38, -27]

julia> k=20;
julia> fdiff_adams_moulton_expansion_coeff(k)
Integer-overflow protection: output converted to BigInt
-12365722323469980029//4817145976189747200000

fdiff_differentiation(f::Vector{T}, v::V; k=3) where {T<:Real, V<:Real}

$k^{th}$-order (default third order) lagrangian differentiation of the analytic function $f$, tabulated in forward order on a uniform grid.

Example:

f = [x^3 for x=-5:5]; println(round.(Int,f))
l = length(f)
f′ = [fdiff_differentiation(f, v) for v=1:l]; println(round.(Int,f′))
  [-125, -64, -27, -8, -1, 0, 1, 8, 27, 64, 125]
  [75, 48, 27, 12, 3, 0, 3, 12, 27, 48, 75]

f′= fdiff_differentiation(f, 1) ; println("f′(1) = $(f′))
  f′(1) = 75//1

f′= fdiff_differentiation(f, 6.5) ; println("f′(6.5) = $(f′)")
    f′(6.5) = 0.75

For a cubic function the third-order lagrangian differentiation is exact - see Figure below.

fdiff_differentiation_expansion_coeffs(ξ::T [, k=3 [, notation=bwd]]) where T<:Real

Finite-difference expansion coefficient vector defining $k^{th}$-order lagrangian differentiation of the tabulated analytic function $f[n]$ at offset $ξ$ (with respect to index position $n$), which is positive for increasing index and negative for decreasing index.

Forward difference notation (notation = fwd)

\[\frac{df}{dξ}[n+ξ]=\sum_{p=0}^kα_p(ξ)Δ^{p}f[n]\]

Offset convention: $ξ = -σ$ with respect to index $n$ in tabulated interval $f[n:n+k]$

Backward difference notation (notation = bwd)

\[\frac{df}{dξ}[n+ξ]=\sum_{p=0}^kβ_p(ξ)∇^{p}f[n]\]

where $β(ξ) ≡ [β_0(ξ),\ ⋯,\ β_p(ξ)]$

Offset convention: $ξ = σ$ with respect to index $n$ in tabulated interval $f[n-k:n]$

Example:

k = 2; ξ = 0
o = fdiff_differentiation_expansion_coeffs(ξ, k); println(o)
 [0.0, 1.0, -1.5]

fdiff_expansion(coeffs, f[, notation=bwd])

Finite difference expansion of the analytical function f(x) tabulated in forward order (growing index) at $k+1$ positions on a uniform grid. The expansion coefficients are specified by the vector coeffs. By default coeffs are assumed to be in backward-difference notation (bwd). For coeffs in forward-difference notation the third argument must be fwd.

Forward difference notation (notation = fwd)

\[\sum_{p=0}^{k}α_{p}Δ^{p}f[n] = F^{k} \cdot f[n:n+k],\]

where $f[n:n+k]$ are elements of the analytical function $f$ (tabulated in forward order) and $α ≡ [α_0,⋯\ α_k]$ is the vector coeffs, which has to be supplied to define the forward-difference expansion. The corresponding weights vector $F^{k}$ is internally generated.

Backward difference notation (notation = bwd)

\[\sum_{p=0}^{k}β_{p}∇^{p}f[n] = \bar{B}^k \cdot f[n-k:n].\]

where $f[n-k:n]$ are elements of the analytical function $f$ (tabulated in forward order) and $β ≡ [β_0,⋯\ β_k]$ is the vector coeffs, which has to be supplied to define the backward-difference expansion. The corresponding weights vector $\bar{B}^k$ is internally generated.

Examples:

Consider the function $f(x)=x^2$ and the expansions,

\[f(x-1)=(1+Δ)^{-1}=(1-Δ+Δ^2-Δ^3+⋯)f(x).\]

\[f(x+1)=(1-∇)^{-1}=(1+∇+∇^2+∇^3+⋯)f(x),\]

To fourth order (k=4) the forward- and backward-difference coefficient vectors are α=[1,-1,1,-1,1] and β=[1,1,1,1,1], respectively. We tabulate the function at $k+1$ points, f=[1,4,9,16,25].

α = [1,-1,1,-1,1]
Fk = fdiff_expansion_weights(α, fwd, reg); println("Fk = $(Fk)")
  Fk = [5, -10, 10, -5, 1]

β = [1,1,1,1,1]
revBk = fdiff_expansion_weights(β, bwd, rev); println("revBk = $(revBk)")
  revBk = [1, -5, 10, -10, 5]

f = [1,4,9,16,25]
o = fdiff_expansion(α, f, fwd); println("f[0] = $(o)")
  f[0] = 0

fdiff_expansion(α, f, fwd) == Fk ⋅ f == fdiff_interpolation(f, 0)
  true

o = fdiff_expansion(β, f, bwd); println("f[6] = $(o)")
  f[6] = 36

fdiff_expansion(β, f, bwd) == revBk ⋅ f == fdiff_interpolation(f, length(f)+1)
  true

In these cases the results are exact because the function is quadratic and the expansion is third order (based on the polynomial of $k^{th}$ degree running through the $k+1$ points of the tabulated function). Note the relation with fdiff_interpolation(f, v, k=3).

fdiff_expansion_weights(coeffs[, notation=bwd[, ordering=rev]])

Expansion weights corresponding to the expansion coefficients coeffs of a finite difference expansion.

Forward difference notation (notation = fwd)

Weight vector $F^k ≡ [F_k^k,⋯\ F_0^k]$ corresponding to the expansion coefficients $α ≡ [α_0^k,⋯\ α_k^k]$ of the $k^{th}$-order forward-difference expansion,

\[\sum_{p=0}^{k}α_{p}Δ^{p}f[n] =\sum_{j=0}^{k}F_{j}^{k}f[n+j] =F^{k} \cdot f[n:n+k],\]

where $f[n:n+k]$ are elements of the analytic function $f$ tabulated in forward order.

fdiff_expansion_weights(α, fwd, reg) $→ F^k ≡ [F_0^k,⋯\ F_k^k]$,

where $α ≡ [α_0,⋯\ α_k]$ has to be supplied in combination with fwd to indicate that the weights must be evaluated in forward-difference notation.

Backward difference notation (notation = bwd)

Weight vector $\bar{B}^{k} ≡ [B_k^k,⋯\ B_0^k]$ corresponding to the expansion coefficients $β ≡ [β_0,⋯\ β_k]$ of the $k^{th}$-order backward-difference expansion,

\[\sum_{p=0}^{k}β_{p}∇^{p}f[n] =\sum_{j=0}^{k}B_{j}^kf[n-j] =\bar{B}^k \cdot f[n-k:n].\]

where $f[n-k:n]$ are elements of the analytic function $f$ tabulated in forward order.

fdiff_expansion_weights(β, bwd, rev) $→ \bar{B}^{k} ≡ [B_k^k,⋯\ B_0^k]$,

where $β ≡ [β_0,⋯\ β_k]$ has to be supplied in combination with bwd to indicate that the weights must be evaluated in backward-difference notation.

Example:

Consider the expansions,

\[f[n-1]=(1+Δ)^{-1}f[n]=(1-Δ+Δ^2-Δ^3+⋯)f[n].\]

\[f[n+1]=(1-∇)^{-1}f[n]=(1+∇+∇^2+∇^3+⋯)f[n],\]

α = [1,-1,1,-1,1]
β = [1,1,1,1,1]
Fk = fdiff_expansion_weights(α, fwd, reg); println("Fk = $(Fk)")
  Fk = [5, -10, 10, -5, 1]

Bk = fdiff_expansion_weights(β, bwd, reg); println("Bk = $(Bk)")
  Bk = [5, -10, 10, -5, 1]

revFk = fdiff_expansion_weights(α, fwd, rev); println("revFk = $(revFk)")
  revFk = [1, -5, 10, -10, 5]

revBk = fdiff_expansion_weights(β, bwd, rev); println("revBk = $(revBk)")
  revBk = [1, -5, 10, -10, 5]

fdiff_interpolation(f::Vector{T}, v::V; k=3) where {T<:Real, V<:Real}

Finite difference lagrangian interpolation (by default third order) at real position v (in index units) with respect to the elements of the uniformly tabulated analytic function f[1:N]. The interpolation points lie on a Lagrange polynomial of degree $k$ (by default third degree) running through $k+1$ subsequenct points of the tabulated function. Outside the tabulated range, the method represents extrapolation on the lagrangian polynomial defined by the first/last $k+1$ tabulated points.

Beware that the interpolation becomes inaccurate if the tabulated function cannot be approximated by a polynomial of degree $k$.

Examples:

f = [1,2,3,4,5,6,7]
[fdiff_interpolation(f, v; k=3) for v=1:0.5:7]
  [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0]

f = [1,4,9,16,25,36,49]
[fdiff_interpolation(f, v; k=3) for v=1:0.5:7]
 [1.0, 2.25, 4.0, 6.25, 9.0, 12.25, 16.0, 20.25, 25.0, 30.25, 36.0, 42.25, 49.0]

 f = [x^3 for x=-4:2]
 f1(v) = fdiff_interpolation(f, v; k=1)
 f2(v) = fdiff_interpolation(f, v; k=2)
 f3(v) = fdiff_interpolation(f, v; k=3)
 [[f1(v),f2(v),f3(v)] for v=1:0.5:9]
   [[-64.0, -64.0, -64.0], [-45.5, -43.25, -42.875], [-27.0, -27.0, -27.0],
   [-17.5, -16.0, -15.625], [-8.0, -8.0, -8.0], [-4.5, -3.75, -3.375],
   [-1.0, -1.0, -1.0], [-0.5, -0.5, -0.125], [0.0, 0.0, 0.0],
   [0.5, -0.25, 0.125], [1.0, 1.0, 1.0], [4.5, 3.75, 3.375], [8.0, 8.0, 8.0],
   [11.5, 13.75, 15.625], [15.0, 21.0, 27.0], [18.5, 29.75, 42.875],
   [22.0, 40.0, 64.0]]

The result for f3(v) is exact because the function is cubic and the expansion is third order - see Figure below. The tabulated function is given by the black points. The interpolation and extrapolation points are red.

fdiff_interpolation_expansion_coeffs(ξ::T [, k=3 [, notation=bwd]]) where T<:Real

Finite-difference expansion coefficient vector defining the $k^{th}$-order (default third order) Lagrange-polynomial interpolation of a tabulated analytic function $f[n]$ at offset $ξ$ with respect to index position $n$, which is positive for increasing index and negative for decreasing index.

Forward difference notation (notation = fwd)

In this case we consider the tabulated interval $f[n:n+k]$. The interpolated value $f[n+ξ]$ is given by the forward-difference expansion

\[f[n+ξ] = \sum_{p=0}^k α_p(-ξ) Δ^p f[n] + ⋯,\]

where the expansion coefficients are given by

fdiff_interpolation_expansion_coeffs(ξ, k, fwd) $→ α(-ξ) ≡ [α_0(-ξ),⋯\ α_k(-ξ)]$. In this notation the range $0\leq ξ\leq k$ corresponds to interpolation and the ranges $ξ<0$ and $ξ>k$ to extrapolation.

Backward difference notation (notation = bwd)

In this case we consider the tabulated interval $f[n-k:n]$. The interpolated value $f[n+ξ]$ is given by the backward-difference expansion

\[f[n+ξ] = \sum_{p=0}^k β_p(ξ) ∇^p f[n] + ⋯,\]

where the expansion coefficients are given by

fdiff_interpolation_expansion_coeffs(ξ, k, bwd) $→ β(ξ) ≡ [β_0(ξ),⋯\ β_k(ξ)]$. In this notation the range $-k\leq ξ\leq0$ corresponds to interpolation and the ranges $ξ<-k$ and $ξ>0$ to extrapolation.

Examples:

k = 5
ξ = -1
α = fdiff_interpolation_expansion_coeffs(ξ, k, fwd); println("α = $α")
β = fdiff_interpolation_expansion_coeffs(ξ, k, bwd); println("β = $β")
  α = [1, 1, 0, 0, 0, 0]
  β = [1, 1, 1, 1, 1, 1]

ξ = 0
α = fdiff_interpolation_expansion_coeffs(ξ, k, fwd); println("α = $α")
β = fdiff_interpolation_expansion_coeffs(ξ, k, bwd); println("β = $β")
  α = [1, 0, 0, 0, 0, 0]
  β = [1, 0, 0, 0, 0, 0]

ξ = 1
α = fdiff_interpolation_expansion_coeffs(ξ, k, fwd); println("α = $α")
β = fdiff_interpolation_expansion_coeffs(ξ, k, bwd); println("β = $β")
  α = [1, -1, 1, -1, 1, -1]
  β = [1, -1, 0, 0, 0, 0]

fdiff_weight(k::Int, j::Int)

Finite difference weight coefficient,

\[c_{j}^{k}=(-1)^{k+j}\binom{k}{j}.\]

Example:

c(k,j) = fdiff_weight(k,j)

[[c(k,j) for j=0:k] for k=0:3] == [[1], [1, -1], [1, -2, 1], [1, -3, 3, -1]]
  true

[[c(k,k-j) for j=0:k] for k=0:3] == [[1], [-1, 1], [1, -2, 1], [-1, 3, -3, 1]]
  true

findIndex(rval::T, grid::Grid{T}) where T<:Number

The grid index corresponding to the position rval on the grid.

Example:

h = 0.1
r0 = 1.0
grid = castGrid(1, 4, Float64; h, r0)
r = grid.r; println("r[3] = $(r[3])")
  Grid created: exponential, Float64, Rmax = 0.491825 a.u., Ntot = 4, h = 0.1, r0 = 1.0
  r[3] = 0.22140275816016985

findIndex(0.222, grid)
  3

find_all(A [,a...]; count=false)

A: string/array of elements of the same type

default : Array containing the index (indices) of selected elements of A (default: all elements)

count=true: The number of indices found for selected elements of A (default: all elements)

Examples:

A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";
find_all(A) == find_all(B) == find_all(str)
true

find_all(A,:📌)
1-element Array{Array{Int64,1},1}:
 [2, 4]

find_all(str)
4-element Array{Array{Int64,1},1}:
 [1]
 [2, 4]
 [3]
 [5]

find_all(A; count=true)
4-element Array{Int64,1}:
 1
 2
 1
 1

str = "📑📌📢📌📞"
find_all(str,'📌')
1-element Array{Array{Int64,1},1}:
 [2, 4]

find_first(A [,a...]; dict=false)

The first index of selected Array element

A: string/array of elements of the same type

default : Array containing the first index (indices) of selected elements of A (default: all elements)

dict=true: Dict for the first index (indices) of selected elements of A (default: all elements)

Examples:

A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";

find_first(A) == find_first(B) == find_first(str)
true

find_first(A,:📌)
1-element Array{Array{Int64,1},1}:
 2

find_last(A,:📌; dict=true)
1-element Array{Pair{Symbol,Int64},1}:
 :📌 => 2

find_last(A; dict=true)
4-element Array{Pair{Symbol,Int64},1}:
 :📑 => 1
 :📌 => 2
 :📢 => 3
 :📞 => 5

find_first(str)
4-element Array{Int64,1}:
 1
 2
 3
 5

find_last(A [,a...]; dict=false)

The last index of selected Array element

A: string/array of elements of the same type

default : Array containing the lasst index (indices) of selected elements of A (default: all elements)

dict=true: Dict for the lasst index (indices) of selected elements of A (default: all elements)

Examples:

A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";
find_last(A) == find_first(B) == find_first(str)
true

find_last(A,:📌)
1-element Array{Array{Int64,1},1}:
 4

find_last(A,:📌; dict=true)
1-element Array{Pair{Symbol,Int64},1}:
 :📌 => 4

find_last(A; dict=true)
4-element Array{Pair{Symbol,Int64},1}:
 :📑 => 1
 :📌 => 4
 :📢 => 3
 :📞 => 5

find_last(str)
4-element Array{Int64,1}:
 1
 4
 3
 5

frac(i)

Fraction notation for rational numbers

Examples:

frac(-5//2)
 "-⁵/₂"

get_Na(Z::Vector{Complex{T}}, def::Def{T}) where T<:Real

Grid index of the starting point for outward numerical integration. This is k+1 or the point marking the end of the quasiclassical region below the lower classical turning point (lctp) as marked by the WKB threshold value (def.pos.cWKB).

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0)
E, def, adams, Z = adams_moulton_master(E, codata, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=false);
    Orbital: 1s
        principal quantum number: n = 1
        radial quantum number: n′ = 0 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 0
    Grid created: exponential, Float64, Rmax = 63.0 a.u., Ntot = 100, h = 0.1, r0 = 0.00286033
    Def created for hydrogen 1s on exponential grid

Na = get_Na(Z, def)
println("k + 1 = $(grid.k+1); Na = $Na")
    k + 1 = 8; Na = 8

Na == def.pos.Na
    true

get_Nb(Z::Vector{Complex{T}}, def::Def{T}) where T<:Real

Grid index of the stopping for outward numerical integration. This is N-k-1 or the point marking the start of the quasiclassical region above the upper classical turning point (Nuctp) as marked by the WKB threshold value (def.pos.cWKB).

Example:

Ecal, grid, def, adams = demo_hydrogen(n=1, ℓ=0)
E, def, adams, Z = adams_moulton_master(E, codata, grid, def, adams; Δν=Value(1,"kHz"), imax=25, msg=false);
    Orbital: 1s
        principal quantum number: n = 1
        radial quantum number: n′ = 0 (number of nodes in radial wavefunction)
        orbital angular momentum of valence electron: ℓ = 0
    Grid created: exponential, Float64, Rmax = 63.0 a.u., Ntot = 100, h = 0.1, r0 = 0.00286033
    Def created for hydrogen 1s on exponential grid

Nb = get_Nb(Z, def)
println("N - k - 1 = $(grid.N-grid.k-1); Nb = $Nb")
    N - k - 1 = 92; Nb = 92

Nb == def.pos.Nb
    true

get_Nlctp(E::T, def::Def{T}) where T<:Real

Grid index of the *lower classical turning point * of the screened potential curve. By definition get_Nlctp(E, def) = 2 for zero orbital angular momentum (ℓ=0).

get_Nmin(def::Def{T}) where T<:Real

Grid index of the minimum of the screened potential curve. By definition get_Nmin(def) = 1 for zero orbital angular momentum (ℓ=0).

get_Nuctp(E::T, def::Def{T}) where T<:Real

Grid index of the upper classical turning point of the screened potential curve. By definition get_Nuctp(E, def) = N-1 for zero orbital angular momentum ($ℓ=0$).

grid_differentiation(f::Vector{T}, grid::Grid{T}; k=3) where T<:Real

$k^{th}$-order lagrangian differentiation of the analytic function $f$, tabulated in forward order on a Grid of $n$ points, $f[1:n]$.

Example:

ID = 4 # linear grid
f = [0.0, 1.0, 4.0, 9.0, 16.0, 25.0]
grid = castGrid(ID, length(f), Float64; r0=1.0, h=1.0, k=3)  # linear grid
f′= grid_differentiation(f, grid; k=3); println("f′= $(f′)")
  Grid created: linear, Float64, Rmax = 6.0 a.u., Ntot = 6, p = 1, h = 1.0, r0 = 1.0
  f′= [0.0, 2.0, 4.0, 6.0, 7.999999999999998, 9.999999999999993]

grid_integration(f::Vector{T}, n1::Int, n2::Int, grid::Grid{V}) where {T<:Real, V<:Real}

Integral of the function $f=[f_0,⋯\ f_n]$ tabulated on a Grid using the trapezoidal rule optimized with endpoint correction by the weightsvector grid.epw,

\[ ∫_{0}^{r_n} f(r) dr = ∫_{0}^{n} f(x) r^{\prime}(x) dx,\]

where the latter integral corresponds to the optimized trapezoidal rule for a uniform grid (see trapezoidal_integration). The rule is exact for polynonials of degree $d=0,\ 1,⋯\ k-1$, where $k=$ grid.epn. For $k=1$ the rule reduces to the ordinary trapezoidal rule (weights = [1/2]).

Example:

f1s(r) = 2.0*r*exp(-r);  # hydrogen 1s wavefunction (reduced and unit normalized)
N = 1000;
grid = castGrid(1, N, Float64; h=0.01, r0=0.005)
    create exponential Grid: Float64, Rmax = 110.127 (a.u.), Ntot = 1000, h = 0.01, r0 = 0.005

r = grid.r;
f2 = [f1s(r[n])^2 for n=1:N];
grid_integration(f2, 1:N, grid) == grid_integration(f2, 1, N, grid)
    true

norm = grid_integration(f2, 1:N, grid)

    1.0

gridfunction(ID::Int, n::Int, h::T; p=5, coords=[0,1], deriv=0) where T <: Real

• ID = 1: exponential grid function,

\[ f[n] = \text{exp}(h(n-1)) - 1.0\]

• ID = 2: quasi-exponential grid function degree p (linear grid for p = 1),

\[ f[n] = h(n-1) + \frac{1}{2}(h(n-1))^2 + ⋯ + \frac{1}{p!}(h(n-1))^p\]

• ID = 3: linear grid function,

\[ f[n] = h(n-1)\]

• ID = 4: polynomial grid function of degree p = length(c) based on polynom $c = [c_1,c_2,⋯\ c_p]$,

\[ f[n] = c_1h(n-1) + c_2(h(n-1))^2 + ⋯ + c_p(h(n-1))^p\]

Examples:

h = 0.1
r = [gridfunction(1, n-1, h) for n=1:5]                            # exponential
 [0.0, 0.10517091807564771, 0.22140275816016985, 0.3498588075760032, 0.49182469764127035]

r = [gridfunction(2, n-1, h; p = 4) for n=1:5]  # quasi exponential (degree p=4)
 [0.0, 0.10517083333333321, 0.22140000000000004, 0.3498375, 0.49173333333333336]

r = [gridfunction(3, n-1, h) for n=1:5]              # linear
  [0.0, 0.1, 0.2, 0.3, 0.4]

r′= [gridfunction(3, n-1, h; deriv=1) for n=1:5]     # linear (first derivative)
   [0.1, 0.1, 0.1, 0.1, 0.1]

  r = [gridfunction(4, n-1, h; coords = [0,1,1/2,1/6,1/24]) for n=1:5]  # polynomial of degree 4)
   [0.0, 0.10517083333333334, 0.2214, 0.3498375000000001, 0.49173333333333336]

gridname(ID::Int)

Name corresponding to the grid ID.

Example:

n = gridname(2); println("The grid type with ID = 2 is called '$n'.")
  The grid type with ID = 2 is called 'quasi-exponential'.

hydrogenic_reduced_wavefunction(Zval, orbit::Orbit, grid::Grid)

Analytic expression for the hydrogenic wavefunction written in the format $Z = \tilde{χ} + i \tilde{χ}^′$, where $\tilde{χ}_{nℓ}(ρ)$ is the reduced radial wavefunction and $\tilde{χ}^′_{nℓ}(ρ)$ its derivative, with $ρ$ the radial distance to the nucleus in a.u.. The expression is evaluated for a given Atom in a given Orbit on a given Grid. The argument Def completes the definition of the problem.

\[ \tilde{\chi}_{nl}(\rho) =\mathcal{N}_{nl}^{-1/2}(2Z/n)^{l+3/2}\rho^{l+1}e^{-Zρ/n} L_{n-l-1}^{2l+1}(2Z\rho/n)\]

where $L_{n-l-1}^{2l+1}(2Z\rho/n)$ is the generalized Laguerre polynomial CamiMath.generalized_laguerreL and

\[ \mathcal{N}_{nl} = {\displaystyle \int\nolimits _{0}^{\infty}}x^{2l+2}e^{-x} \left[L_{n-l-1}^{2l+1}(x)\right]^{2}dx = \frac{2n\Gamma(n+l+1)}{\Gamma(n-l)}\]

is the norm of the wavefunction.

Example:

orbit = castOrbit(n=25, ℓ=10)
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=true)
def = castDef(grid, atom, orbit, codata)
Zval = 1
Z = hydrogenic_reduced_wavefunction(Zval, orbit, grid);
    Orbital: 25n
    principal quantum number: n = 25
    radial quantum number: n′ = 14 (number of nodes in radial wavefunction)
    orbital angular momentum of valence electron: ℓ = 10
    Grid created: exponential, Float64, Rmax = 1935.0 a.u., Ntot = 1300, h = 0.00769231, r0 = 0.0878529
    Def created for hydrogen 25n on exponential grid of 1300 points

plot_wavefunction(Z, 1:grid.N, grid, def)

The plot is made using CairomMakie. NB.: plot_wavefunction is not included in the CamiXon package.

initE(def::Def{T}) where T<:Real

Autogenerated seed value for the energy

Example:

codata = castCodata(2018)
atom = castAtom(Z=1, A=1, Q=0; msg=false)
orbit = castOrbit(n=1, ℓ=0; msg=false)
grid = autoGrid(atom, orbit, Float64; msg=false)
def = castDef(grid, atom, orbit, codata);
    Def created for hydrogen 1s on exponential grid of 100 points

E = initE(def); println("E = $E")
    E = -0.03508495857961283

function isforward(val)

Boolean status of val, with options: fwd (forward) and bwd (backward).

Example:

julia? isforward(fwd)
true

function isregular(val)

Boolean status of val, with options: reg (regular) and rev (reversed).

Example:

isregular(reg)
  true

latent_heat_vaporization(atomicnumber::Int, temp::Real)
latent_heat_vaporization(element::String, temp:Real)

Latent heat of vaporization, L(T) (in Joule/K), of a given element for temperature T (in K),

\[L(T) = -(B +C\cdot T \mathrm{log_{10}}T+D\cdot T^2/1000),\]

where A,B,C,D, are the Antoine coefficients collected in CamiXon.dictAntoineCoefficients.

Example:

julia> latent_heat_vaporization("Li", 623.0)
-18473.64020109123

latexIsotopeTable(Z1::Int, Z2::Int; continuation=false)
latexIsotopeTable(itrZ::UnitRange; continuation=false)

Isotope table for all isotopes with atomic number from Z1 to Z2.

Example:

o = latexIsotopeTable(1:3);
println(o)
  \setlength{\tabcolsep}{3pt}
  \renewcommand{\arraystretch}{1.2}
  \begin{table}[H]
    \centering
    \caption{\label{table:Isotopes-a-1}Properties of selected atomic isotopes. The Table is based on three databases: (a) AME2020 (atomic mass evaluation); (b) IAEA-INDC(NDS)-794 (magnetic dipole moments); (c) IAEA-INDC(NDS)-833 (electric quadrupole moments).}
    \begin{tabular}{r|lr|rrrr|r|r|r|r}
      \multicolumn{12}{r}\vspace{-18pt}\\
      \hline
      \hline
      $Z$ & element & symbol & $A$ & $N$ & radius & atomic mass & $I\,^\pi$ & $\mu_I $ & $Q$ & $RA$\\&  &  &  &  & (fm) & $(m_u)$ & $(\hbar)\ \ $ & $(\mu_N)$ & (barn) & (\%)\\\hline
      1 & hydrogen & $^{1}$H & 1\, & 0 & 0.8783 & 1.007825032 & 1/2$^+$ & 2.792847351 & 0.0 & 99.9855 \\
        &  & $^{2}$H & 2\, & 1 & 2.1421 & 2.014101778 & 1//1$^+$ & 0.857438231 & 0.0028578 & 0.0145 \\
        &  & $^{3}$H & 3$*\!\!$ & 2 & 1.7591 & 3.016049281 & 1/2$^+$ & 2.97896246 & 0.0 & trace \\
      \hline
      2 & helium & $^{3}$He & 3\, & 1 & 1.9661 & 3.016029322 & 1/2$^+$ & -2.12762531 & 0.0 & 0.0002 \\
        &  & $^{4}$He & 4\, & 2 & 1.6755 & 4.002603254 & 0//1$^+$ & 0.0 & 0.0 & 99.9998\% \\
      \hline
      3 & lithium & $^{6}$Li & 6\, & 3 & 2.589 & 6.015122887 & 1//1$^+$ & 0.822043 & -0.000806 & 4.85 \\
        &  & $^{7}$Li & 7\, & 4 & 2.444 & 7.016003434 & 3/2$^-$ & 3.256407 & -0.04 & 95.15 \\
      \hline
      \multicolumn{12}{l}{*radioactive }\\
    \end{tabular}
  \end{table}

The typeset result is shown in the figule below.

lc_eltype(o)

Lowest common eltype of a collection.

Examples:

julia> o = ([1//2, 1//3]; (1//4, 1//1, 1//6));
julia> lc_eltype(o)
Rational{Int64}

julia> o = ([1//2, 1//3]; (1//4, big(1)//big(5), 1//6));
julia> lc_eltype(o)
Rational

julia> o = ([1//2, 1//3]; (1//4, [big(1)//big(5)], 1//6));
julia> lc_eltype(o)
Any

julia> o = ([1/2, 1/3]; (1/4, 1/1, 1/6));
julia> lc_eltype(o)
Float64

lc_primitivetype(o::Any)

Lowest comon primitive type of Any Type

Examples:

julia> o = ([1//2, 1//3]; (1//4, 1//1, 1//6));
julia> lc_primitivetype(o)
Int64

listAtom(Z::Int, A::Int, Q::Int[; fmt=Object])

Properties of atom with atomic number Z, atomic mass number A, ionic charge Q.

Output options: fmt = Object (default), String, Info.

Example:

listAtom("H", 3, 0) == listAtom(1, 3, 0)
  true

listAtom(1, 3, 0; fmt=Info)
Element: hydrogen
    symbol: H
    element: tritium
    atomic number: Z = 1
    atomic weight (relative atomic mass): 1.008

listAtoms(Z1::Int, Z2::Int, Q::Int[; fmt=Object])

Properties of atoms with atomic number in the range Z1:Z3 and ionic charge Q.

Output options: fmt = Object (default), String, Info.

Example

listAtoms(1,3,0) == listAtoms(1:3,0)
  true

listAtoms(1:1, 0; fmt=Info);
  Atom: hydrogen, neutral atom
    symbol: ¹H
    atomic charge: Z = 1
    Rydberg charge: Zc = 1
  Atom: deuterium, neutral atom
    symbol: ²D
    atomic charge: Z = 1
    Rydberg charge: Zc = 1
  Atom: tritium, neutral atom
    symbol: ³T
    atomic charge: Z = 1
    Rydberg charge: Zc = 1

listCodata(codata::Codata)

Method to list the fields of Codata by their symbolic name

Example:

julia> codata = castCodata(2018);

julia> listCodata(codata)
∆νCs = 9192631770 Hz      - ¹³³Cs hyperfine transition frequency
   c = 299792458 m s⁻¹    - speed of light in vacuum
   h = 6.62607e-34 J Hz⁻¹ - Planck constant
   ħ = 1.05457e-34 J s    - Planck constant (reduced)
   e = 1.60218e-19 C      - elementary charge
  kB = 1.38065e-23 J K⁻¹  - Boltzmann constant
  NA = 6.02214e23 mol⁻¹   - Avogadro constant
 Kcd = 683 lm W⁻¹         - Luminous efficacy
  mₑ = 9.10938e-31 Kg     - electron rest mass
  R∞ = 1.09737e7 m⁻¹      - Rydberg constant
  Ry = 3.28984e15 Hz      - Rydberg frequency
  Eₕ = 4.35974e-18 Hartree a.u. - Hartree atomic unit
   α = 0.00729735         - fine-structure constant
  μ₀ = 1.25664e-6 N A⁻²   - magnetic permitivity of vacuum
  ε₀ = 8.85419e-12 F m⁻¹  - electric permitivity of vacuum
  KJ = 4.83598e14 Hz V⁻¹  - Josephson constant
  RK = 25812.8 Ω          - Von Klitzing constant
   R = 8.31446 J mol⁻¹K⁻¹ - Molar gas constant
   u = 1.66054e-27 kg     - unified atomic mass unit

julia> codata.u.val
1.6605390666e-27  

listElement(Z::Int[; fmt=Object])
listElement(elt::String[; fmt=Object])

Properties of element with atomic number Z.

Output options: fmt = Object (default), String, Info.

Example:

listElement("H") == listElement(1)
  true

listElement(1; fmt=Info)
  Element: hydrogen
    symbol: H
    atomic number: Z = 1
    atomic weight (relative atomic mass): 1.008

listElements(Z1::Int, Z2::Int[; fmt=Object])
listElements(itrZ::UnitRange{Int}; fmt=Object)

Properties of elements with atomic number in the range Z1:Z2.

Output options: fmt = Object (default), String, Info.

Example

listElements(1,3) == listElements(1:3)
  true

listElements(1:3; fmt=Info)
  Element: hydrogen
    symbol: H
    atomic number: Z = 1
    atomic weight (relative atomic mass): 1.008
  Element: helium
    symbol: He
    atomic number: Z = 2
    atomic weight (relative atomic mass): 4.0026
  Element: lithium
    symbol: Li
    atomic number: Z = 3
    atomic weight (relative atomic mass): 6.94

listIsotope(Z::Int, A::Int; fmt=Object)

Properties of isotopes with atomic number Z and atomic mass number A.

Output options: fmt = Object (default), String, Latex, Info.

Example:

listIsotope(1,3; fmt=Info)
  Isotope: tritium-3
    symbol: ³T
    element: tritium
    atomic number: Z = 1
    atomic mass number: A = 3
    neutron number: N = 2
    rms nuclear charge radius: R = 1.7591 fm
    atomic mass: M = 3.016049281 amu
    nuclear spin: I = 1/2 ħ
    parity of nuclear state: π = even
    nuclear magnetic dipole moment: μI = 2.97896246μN
    nuclear electric quadrupole moment: Q = 0.0barn
    relative abundance: RA = trace
    lifetime: 12.33 years

listIsotopes(Z1::Int, Z2::Int; fmt=Object)

All isotopes with atomic number from Z1 to Z2.

Output options: Object (default), String, Latex, Info.

Example:

listIsotopes(1,3) == listIsotopes(1:3)
 true

listIsotopes(1:1; fmt=Info)
 3-element Vector{Any}:
  Isotope("¹H", "hydrogen", 1, 1, 0, 0.8783, 1.007825032, 1//2, 1, 1.0e100, 2.792847351, 0.0, 99.9855)
  Isotope("²D", "deuterium", 1, 2, 1, 2.1421, 2.014101778, 1, 1, 1.0e100, 0.857438231, 0.0028578, 0.0145)
  Isotope("³T", "tritium", 1, 3, 2, 1.7591, 3.016049281, 1//2, 1, 12.33, 2.97896246, 0.0, nothing)

matG(E::T, grid::Grid{T}, def::Def{T}) where T<:Real

matMinv(E::T, grid::Grid{T}, def::Def{T}, amEnd::T) where T<:Real

matσ(E::T, grid::Grid{T}, def::Def{T}) where T<:Real

melting_point(atomicnumber::Int)
melting_point(element::String)

Melting point of a given element at standard pressure (1 atm).

Example:

julia> melting_point("Li")
453.65

primitivetype(T::Type)

The primitive type of a Type

Examples:

julia> T = Complex{Float16}};
julia> primitivetype(T)
Float16

julia> T = String
julia> primitivetype(T)
Char

reduce_wavefunction(Z::Vector{Complex{T}}, grid::Grid{V}) where {T<:Real, V<:Real}

Conversion from the ordinary radial wavefunction $\tilde{R}_{nl}(ρ)$ to the reduced radial wavefuntion

\[ \tilde{\chi}_{nl}(ρ) = ρ \tilde{R}_{nl}(ρ).\]

where $ρ$ is the radial distance to the nucleus in a.u..

Example:

atom = castAtom(Z=1, A=1, Q=0; msg=false);
orbit = castOrbit(n=1, ℓ=0; msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=false);
def = castDef(grid, atom, orbit, codata);

RH1s_example = [RH1s(atom.Z, grid.r[n]) for n=1:grid.N];
ZH1s_generic = hydrogenic_reduced_wavefunction(1, orbit, grid);

ZH1s_example = reduce_wavefunction(RH1s_example, grid);
RH1s_generic = restore_wavefunction(ZH1s_generic, grid);

ZH1s_example ≈ ZH1s_generic
    true

RH1s_example ≈ RH1s_generic
    true

restore_wavefunction(Z::Vector{Complex{T}}, grid::Grid{V}) where {T<:Real, V<:Real}

Conversion from the reduced radial wavefunction $\tilde{\chi}_{nl}(ρ)$ to the ordinary radial wavefuntion $\tilde{R}_{nl}(ρ)$,

\[ \tilde{R}_{nl}(ρ)=\tilde{\chi}_{nl}(ρ)/ρ,\]

where $ρ$ is the radial distance to the nucleus in a.u..

Example:

atom = castAtom(Z=1, A=1, Q=0; msg=false);
orbit = castOrbit(n=1, ℓ=0; msg=false);
grid = autoGrid(atom, orbit, Float64; Nboost=1, msg=false);
def = castDef(grid, atom, orbit, codata);

RH1s_example = [RH1s(atom.Z, grid.r[n]) for n=1:grid.N];
ZH1s_generic = hydrogenic_reduced_wavefunction(1, orbit, grid);

ZH1s_example = reduce_wavefunction(RH1s_example);
RH1s_generic = restore_wavefunction(ZH1s_generic);

RH1s_example ≈ RH1s_generic
    true

ZH1s_example ≈ ZH1s_generic
    true

f1 = real(ZH1s_example)
f2 = real(ZH1s_generic)

compare_functions(f1, f2, 1:grid.N, grid)

The plot is made using CairomMakie. NB.: compare_functions is not included in the CamiXon package.

select125(x)

Select elements of the collection x by index according to 1-2-5 scheme

Examples:

x = [1,2,4,6,8,10,13,16,18,20,40,60,80,100]
select125(x)
 [2, 6, 10, 16, 20, 60, 100]

x = string.(x)
select125(x)
 ["2", "6", "10", "16", "20", "60", "100"]

x = 1:100
select125(x)
 [20, 40, 60, 80, 100]

step125(x)

Step used for deviding the number x in steps according to 1-2-5 scheme

Examples:

step125.([5,10,21.3,50,100.1])
5-element Vector{Int64}:
  1
  2
  5
 10
 20

stepcenters(x)

Stepcenter positions for steplength specification vector x

Examples:

x = [4,2,6]
stepcenters(x)
 [2.0, 5.0, 9.0]

stepedges(x)

Stepedges for steplength specification vector x

Examples:

x = [4,2,6]
stepedges(x)
 [0, 4, 6, 12]

steps(x)

Heatmap range transformation for steplength specification vector x

Examples:

x = [4,2,6]
steps(x)
 [0, 4, 6, 12]

strRational(i)

Fraction notation for rational numbers and integers

Examples:

strRational(-5//2)
 "-5/2"

 strRational(-5//2)
  "-5"

strValue(f::Value)

String expression for a Value object in :compact => true representation

Example:

f = Value(1,"Hz")
strValue(f)
  "1 Hz"

sub(i::T) where T<:Real
sub(str::String)

Subscript notation for integers, rational numbers and a subset of lowercase characters ('a','e','h','k','l','m','n','o','p','r','s','t','x')

Examples:

'D' * sub(5//2)
 "D₅⸝₂"

"m" * sub("e")
 "mₑ"

sup(i::T) where T<:Real

Superscript notation for integers and rational numbers

Examples:

sup(3) * 'P'
 "³P"

svp(atomicnumber::Int, temp::Real)
svp(element::String, temp::Real)

Saturated vapor pressure, p (in Pa), of a given element for a given temperature T (in K),

\[\mathrm{log_{e}}p=A+B/T+C\mathrm{log_{10}}T+D\cdot T/1000,\]

where A,B,C,D, are the Antoine coefficients collected in CamiXon.dictAntoineCoefficients.

Examples:

julia> svp("Li", 623.0)
0.0015230367024569058

threeJsymbol(j1::Real, m1::Real, j2::Real, m2::Real, j3::Real, m3::Real; msg=false)

Wigner 3j symbol. This is a vector coupling coefficient with optimized symmetry properties. The 3j symbols are zero unless $Δ(j_{1},j_{2},j_{3})>0$ (triangle inequality holds) and $m_{1}+m_{2}+m_{3}=0$. The implementation is based on the Racah formula:

\[\left(\begin{array}{ccc} j_{1} & j_{2} & j_{3}\\ m_{1} & m_{2} & m_{3} \end{array}\right)= (-1)^{j_{1}-j_{2}-m_{3}}\sqrt{\Delta(j_{1}j_{2}J)}\\\times \sqrt{\left(j_{1}+m_{1}\right)! \left(j_{1}-m_{1}\right)! \left(j_{2}+m_{2}\right)! \left(j_{2}-m_{2}\right)! \left(j_{3}+m_{3}\right)! \left(j_{3}-m_{3}\right)!} \\\times\sum_{t}\frac{(-)^{t}}{t!(j_{3}-j_{2}+t+m_{1})! (j_{3}-j_{1}+t-m_{2})! (j_{1}+j_{2}-j_{3}-t)!(j_{1}-t-m_{1})!(j_{2}-t+m_{2})!}\]

Example:

julia> threeJsymbol(0, 0, 0, 0, 0, 0)
1.0

julia> threeJsymbol(3, 0, 4, -1, 5, 1)
-0.1096417439724123565166029917781360897459044055433631161836138910409772907333476

julia> o = threeJsymbol(3, 0, 4, -1, 5, 1; msg=true); println(" = $o")
-√(361/30030) = -0.1096417439724123565166029917781360897459044055433631161836138910409772907333476

trapezoidal_epw(k::Int [; rationalize=false [, devisor=false]])

Endpoint weights vector $a=[a_1,⋯\ a_k]$ of trapeziodal rule optimized for functions of polynomial form,

\[ ∫_0^n f(x) dx = a_1 (f_0+f_n) + ⋯ + a_k (f_{k-1}+f_{n-k+1}) + (f_k+⋯+f_{n-k}),\]

where $k$ is odd. The rule is exact for polynonials of degree $d=0,\ 1, ⋯,\ k-1$. For $k=1$ the rule reduces to the ordinary trapezoidal rule. By default the output is in Float64, optionally the output is rational, with or without specification of the gcd devisor.

Example:

[trapezoidal_epw(k; rationalize=true, devisor=true) for k=1:2:9]
5-element Vector{Tuple{Int64, Int64, Vector{Int64}}}:
  (1, 2, [1])
  (3, 24, [9, 28, 23])
  (5, 1440, [475, 1902, 1104, 1586, 1413])
  (7, 120960, [36799, 176648, 54851, 177984, 89437, 130936, 119585])
  (9, 7257600, [2082753, 11532470, 261166, 16263486, -1020160, 12489922,
                                                     5095890, 7783754, 7200319])

trapezoidal_integration(f, x1, x2, weights)

Integral of the tabulated function $f=[f_0,⋯\ f_n]$ over the domain $x1 ≤ x ≤ x2$ using the optimized trapezoidal rule with endpoint correction by the weights vector weights,

\[ ∫_0^n f(x) dx = a_1 (f_0+f_n) + ⋯ + a_k (f_{k-1}+f_{n-k+1}) + (f_k+⋯+f_{n-k}).\]

The rule is exact for polynonials of degree $d=0,\ 1,⋯\ k-1$. For $k=1$ the rule reduces to the ordinary trapezoidal rule (weights = [1/2]).

Examples::

p = 3
c = [1 for i=0:p]
pol = ImmutablePolynomial(c,:z)
Ipol = integrate(pol)
n = 10

x1=0.0
x2=5.0
x = collect(range(x1, x2, n))
f = pol.(x .-2.5)

w3 = trapezoidal_epw(3)
trapezoidal_integration(f, x1, x2, w3)
 15.416666666666673

Ipol(2.5)-Ipol(-2.5)
 15.41666666666666

updateAdams!(adams::Adams{T}, E, grid::Grid{T}, def::Def{T}) where T<:Real