API reference

numerov(inputFileName::String)

Run a full Numerov calculation: reads the Numerov input file, solves the 1D, 2D or 3D time-independent Schrödinger equation on the grid potential defined therein and writes the resulting eigenvalues, eigenvectors and frequencies - and, if requested, the band structure - to .dat files together with a log file (default: Numerov.out).

Arguments

• inputFileName::String: path to the Numerov input file.

solve_schrodinger(V, coords; kwargs...) -> SchrodingerSolution

Solve the time-independent Schrödinger equation for the potential V given on a regular grid and return the lowest eigenpairs, without reading or writing any files.

Arguments

• V::AbstractArray{<:Real}: the potential on the grid, as a 1-, 2- or 3-dimensional array. V[i, j, k] is the value at (coords[1][i], coords[2][j], coords[3][k]).
• coords: the grid axes — one AbstractVector/AbstractRange per dimension of V (a single vector may be passed directly for 1D). Each axis must be ascending and equally spaced.

Keyword arguments

• mass = 1.0: reduced mass, a scalar or one value per dimension.
• periodic = false: periodic boundary conditions, a scalar or one flag per dimension.
• n_eigenvalues = 5: number of eigenpairs to compute.
• k = nothing: solve at this k-point (a vector in inverse coordinate units); requires at least one periodic dimension. nothing solves the real, non-periodic problem.
• stencil = 9: finite-difference stencil size (3, 5, 7, 9, 11 or 13); stencil_laplace and stencil_nabla override it individually.
• solver = :arpack: eigensolver backend (:arpack, :krylov, :lobpcg or :lu); :lobpcg is recommended for large non-periodic 3D problems.
• potential_unit = u"hartree", coord_unit = u"bohr", mass_unit = u"m_e": units of the inputs; energies are returned in potential_unit.

The grid spacings must be equal across dimensions after mass weighting, i.e. step(coords[i]) * sqrt(mass[i]) must be the same for every dimension — this is a restriction of the underlying discretization.

Example

x = range(-5.0, 5.0; length = 201)
result = solve_schrodinger(0.5 .* x .^ 2, x; n_eigenvalues = 3)
result.energies  # ≈ [0.5, 1.5, 2.5]

solve_schrodinger(V::AbstractArray{<:Unitful.Energy}, coords; kwargs...)

Convenience method accepting a potential of Unitful energy quantities and coordinate axes of length quantities; the units are inferred from the inputs.

band_structure(V, coords; n_kpoints, kwargs...) -> BandStructure

Compute the band structure of a periodic system along the path through the high-symmetry points of the Brillouin zone, without reading or writing any files.

Takes the same potential/coordinate arguments and keyword arguments as solve_schrodinger (except k), plus:

• n_kpoints::Integer: number of k-points per path segment (must be > 1).

At least one dimension must be periodic.

SchrodingerSolution

Result of solve_schrodinger.

Fields

• energies::Vector{Float64}: the lowest eigenvalues in the unit the potential was given in (absolute, i.e. referenced to the same zero as the input potential).
• states::Matrix: the corresponding eigenvectors as columns, in the same grid ordering as the flattened input potential array. Normalized so that sum(abs2, state) * prod(spacings) == 1 with the spacings expressed in the unit the coordinates were given in. Real for non-periodic calculations, complex for periodic ones.
• kpoint::Union{Nothing, Vector{Float64}}: the k-point the equation was solved at (in inverse coordinate units), or nothing for non-periodic calculations.

BandStructure

Result of band_structure.

Fields

• kpoints::Vector{Vector{Float64}}: the k-points along the path through the high-symmetry points of the Brillouin zone, in inverse coordinate units.
• kpath::Vector{Float64}: cumulative distance along the k-path (inverse coordinate units), useful as the x axis of a band-structure plot.
• energies::Matrix{Float64}: band energies, one row per k-point and one column per band, in the unit the potential was given in (absolute).

KineticPreconditioner

Tensor-product preconditioner for the shifted kinetic operator, used to precondition the LOBPCG eigensolver.

T̃, the Kronecker SUM of the per-dimension 1D operators t_d = -Δ_d / (2 Δq_d²), replaces the true (possibly non-separable) kinetic operator: T̃'s eigendecomposition factorizes into the per-dimension eigenpairs, so (T̃ + σI)⁻¹ x is applied exactly - with one small dense eigenbasis transform per dimension - without ever factorizing a full-size operator, so there is no fill-in.

T̃ equals the true production kinetic operator exactly for 1D problems (trivially) and for the 2D 5-point stencil (buildLaplace_2d's 5-point Laplacian is a pure row/column "cross" pattern, which is separable: dividing by 2^(dimension-1), as solve() does, recovers the Kronecker sum exactly). Every other case - 2D with stencil 3, 7, 9 (the package default) or 11, and 3D with any stencil (buildLaplace_3d always uses a more elaborate, non-separable stencil) - is only an APPROXIMATION; this does not compromise correctness, only convergence speed, since every LOBPCG result is independently verified against the true Hamiltonian's residual in solveWrapper (and re-solved with Arpack if that check fails) regardless of how good an approximation T̃ is.

KineticPreconditioner(potential, system; σ = 1.0)

Build the preconditioner for the given problem. σ > 0 keeps the operator safely positive definite; its exact value only affects convergence speed.

accumulate_stencil_blocks!(Is, Js, Vs, block_size, n_outer, total_points,
                            periodic_outer, sub_matrix_indices, block_triplets)

Shared core of build_2d_stencil and build_3d_stencil: given the (I, J, V) nonzero triplets of each stencil offset's sub-block, plus the row/col block bookkeeping the caller works out, append the final triplets for one "outer axis" level of block placement to Is/Js/Vs.

Building the result via COO triplet accumulation followed by a single sparse(...) call (instead of the equivalent-but-much-slower repeated matrix[range_i, range_j] = block assignment into an already-existing large sparse matrix - each such assignment can touch and reallocate a large fraction of the matrix's internal CSC storage) is what makes this fast: it turns an operation that scaled worse than quadratically in the grid size into one that is close to linear.

For each row-block j, a small per-j Dict (at most stencil_size entries, discarded every iteration) resolves which stencil offset wins the placement for a given target column-block before any triplets are appended, exactly reproducing the original's matrix[...] = block last-write-wins overwrite semantics - this matters only when a periodic wraparound collides with a non-wrapped placement, which can only happen for grids smaller than the stencil width, but is reproduced correctly regardless.

Apply (T̃ + σI)⁻¹ to a flattened grid vector via per-dimension eigenbasis transforms.

build_1d_stencil(system, stencilCoefficients, stencil_size)

Builds a 1D stencil matrix for the given system.

Arguments

- `system::System`: The system to build the stencil matrix for.
- `stencilCoefficients::Array{Float64,1}`: The stencil coefficients.
- `stencil_size::Int`: The size of the stencil.

Returns

- `matrix::SparseMatrixCSC{Float64,Int}`: The 1D stencil matrix.

build_2d_stencil(system, n_datapoints, stencil, stencil_size)

Builds a 2D stencil matrix for the given system.

Arguments

- `system::System`: The system to build the stencil matrix for.
- `n_datapoints::Tuple{Int,Int}`: The number of datapoints in each dimension.
- `stencil::Array{Array{Float64,1},2}`: The stencil coefficients.
- `stencil_size::Int`: The size of the stencil.

Returns

- `matrix::SparseMatrixCSC{Float64,Int}`: The 2D stencil matrix.

build_3d_stencil(system, n_datapoints, stencil, stencil_size)

Builds a 3D stencil matrix for the given system.

Arguments

- `system::System`: The system to build the stencil matrix for.
- `n_datapoints::Tuple{Int,Int,Int}`: The number of datapoints in each dimension.
- `stencil::Array{Array{Array{Float64,1},2},3}`: The stencil coefficients.
- `stencil_size::Int`: The size of the stencil.

Returns

- `matrix::SparseMatrixCSC{Float64,Int}`: The 3D stencil matrix.

buildΔ(system::System, potential::Potential)

Builds the Laplace operator matrix Δ for the given system and potential.

Arguments

- `system::System`: The system to build the Laplace operator for.
- `potential::Potential`: The potential to build the Laplace operator for.

Returns

  - `nothing`

buildΔ_1D(system::System)

Build the Laplace operator for a 1D system. The Laplace operator is built using a finite difference stencil of the order system.stencilΔ. The Laplace operator is stored in system.Δ. Stencil size options are 3, 5, 7, 9, 11 and 13 points.

Arguments

- `system::System`: The system to build the Laplace operator for.

buildΔ_2D(system::System)

Build the Laplace operator for a 2D system. The Laplace operator is built using a finite difference stencil of the order system.stencilΔ. The Laplace operator is stored in system.Δ. Stencil size options are 3, 5, 7, 9 and 11 points.

Arguments

- `system::System`: The system to build the Laplace operator for.

Returns

- `nothing`

buildΔ_3D(system::System)

This file contains the functions to build the 3D Laplacian matrix for the system. The Laplacian matrix is built using a finite difference stencil of the order system.stencilΔ. The Laplacian matrix is stored in system.Δ. Stencil size options are 3, 5, 7, 9 and 11 points.

Arguments

- `system::System`: The system to build the Laplacian matrix for.

Returns

- `nothing`

build∇(system::System, potential::Potential)

Builds the gradient operator matrix ∇ for the given system and potential.

Arguments

- `system::System`: The system to build the gradient operator for.
- `potential::Potential`: The potential to build the gradient operator for.

Returns

- `nothing`

build∇_1D(system::System)

Builds the gradient operator matrix ∇ for the given system in 1D.

Arguments

- `system::System`: The system to build the gradient operator for.

Returns

- `nothing`

build∇_2D(system::System)

Builds the gradient operator matrix ∇ for the given system in 2D.

Arguments

• system::System: The system to build the gradient operator for.

Returns

- `nothing`

build∇_3D(system::System)

Builds the gradient operator matrix ∇ for the given system in 3D.

Arguments

- `system::System`: The system to build the gradient operator for.

Returns

- `nothing`

build_1d_stencil(system::System, stencil_1d::Array{Float64, 1}, stencil_size::Int64)

Builds the 1D gradient operator matrix ∇ for the given system.

Arguments

- `system::System`: The system to build the gradient operator for.
- `stencil_1d::Array{Float64, 1}`: The 1D stencil kernel array.
- `stencil_size::Int64`: The size of the stencil kernel.

Returns

- `∇::SparseMatrixCSC{Float64, Int64}`: The 1D gradient operator matrix.

Per-point column of axis d for a grid flattened with the last dimension varying fastest.

Generate the internal (mass-weighted) k-points for a band-structure run, mirroring readPotential.jl:121-168.

get_1d_stencil(system::System)

Returns the 1D stencil kernel array for the given system. The stencil kernel array is determined by the stencil∇ attribute of the system. Available stencil kernels are 3, 5, 7, 9, and 11 points.

Arguments

- `system::System`: The system to get the 1D stencil kernel array for.

Returns

- `stencil_1d::Array{Float64, 1}`: The 1D stencil kernel array.

1D kinetic matrix -Δ/(2 Δq²) for one dimension, built with the same stencil machinery as the full operator.

Largest relative eigenpair residual max ‖Hx - λx‖ / (‖x‖ max(1, |λ|)). Degenerate (near-zero) eigenvectors count as infinitely loose, so silently collapsed solver output can never pass verification.

Validate the array inputs and build fully initialized Potential, System, Output and Files structs, mirroring exactly what readInputFile + checkInput + readPotential + setupSystem produce for the file pipeline (see those files for the unit conventions). Also builds the stencil operator matrices.

Arpack shift-invert about a small negative σ. The potential is shifted so that min(V) = 0, making the Hamiltonian positive (semi)definite: its smallest eigenvalues are the ones closest to σ ≈ 0, so shift-invert converges in a few iterations where the plain :SM mode needs thousands of restarts. σ sits slightly BELOW zero so that H - σI stays safely invertible even when H itself is exactly singular.

Solve at one internal (mass-weighted) k-point, replicating the shift dance of the main.jl k-loop, and return (energies in hartree with the potential shift restored, states as normalized matrix columns).

Multiply every mode of X by op(Q_d) for its dimension's eigenbasis.