Discretization

Assuming spherical symmetry, each Kohn–Sham orbital factors into a radial part u(r) and a spherical harmonic. AtomicKohnSham.jl discretizes the radial part u(r) with a finite element method (FEM): a mesh cuts the radial domain [0, Rmax] into cells, and a polynomial basis is built on top of it.

Mesh

Several radial mesh types are available; expmesh (exponentially graded, denser near the origin where orbitals vary fastest) is recommended for atomic problems:

AtomicKohnSham.Mesh — Type

Mesh(points::AbstractVector{T}, name::AbstractString = "", params::NamedTuple = NamedTuple{}()) where T <: Real

A container for a one-dimensional mesh used in numerical discretizations.

Arguments

points::AbstractVector{T}: The coordinates of the mesh points, typically sorted in ascending order.
name::AbstractString = "": An optional name or identifier for the mesh.
params::NamedTuple = NamedTuple{}(): A named tuple storing the parameters used to generate the mesh (e.g., domain bounds, number of points, transformation options, etc.).

Fields

name::AbstractString: Name of the mesh (for logging or reconstruction).
points::AbstractVector{T}: Coordinates of the mesh points.
params::NamedTuple: Parameters used to build the mesh, useful for reproducibility.

source

AtomicKohnSham.expmesh — Function

expmesh(a::Real, b::Real, n::Int; T = Float64, s::Real)

Build an exponentially graded Mesh of n points on [a, b], denser near a for s > 1. Recommended for atomic radial problems, with s between 1 and 2, to resolve the orbitals' fast variation near the nucleus while still covering a large Rmax.

source

linmesh, geometricmesh, polynomialmesh and explinmesh follow the same (a, b, n; kwargs...) -> Mesh convention for, respectively, a uniform, geometric, polynomially-graded, or exponential-then-linear mesh.

FEM basis

Currently, the only implemented basis is a hierarchical integrated Legendre polynomial basis (P1 vertex functions plus higher-order "bubble" functions per cell, up to ordermax):

AtomicKohnSham.P1IntLegendreBasis — Function

P1IntLegendreBasis(mesh::Mesh, T::Type = Float64; ordermax::Int)

Build the hierarchical FEM basis on mesh: a P1 "hat" function at each mesh node plus, within each cell, integrated Legendre polynomial "bubble" functions up to degree ordermax, all continuous but with (in general) discontinuous derivatives across cell boundaries.

source

Integration method

FEM matrix assembly and energy integrals need numerical quadrature:

AtomicKohnSham.GaussLegendre — Type

GaussLegendre(basis::FEMBasis, npoints::Int = 1000)

Gauss–Legendre quadrature method for FEM matrix assembly and energy integrals, with npoints points per mesh cell.

Precomputes and caches the quadrature nodes/weights (rescaled to each reference cell and to the full radial domain [0, Rmax]), together with the FEM generator polynomials evaluated at the nodes, so that repeated assemblies during the SCF loop avoid recomputing them.

source

Putting it together

AtomicKohnSham.KSEDiscretization — Type

KSEDiscretization

Finite element discretization of the radial Kohn–Sham equations with spherical symmetry.

This structure defines the numerical discretization used to represent radial Kohn–Sham orbitals in a basis of finite element functions combined with spherical harmonics. It stores the FEM operators, Kohn–Sham Hamiltonian components, and associated workspaces required for SCF computations.

The discretization is independent of the SCF algorithm and can be reused with different solvers.

Mathematical form

The radial part of a Kohn–Sham orbital is expanded as

u(r,θ,φ) = ∑_{l,m,n} u_{nlm} Qₙ(r) / r · Yₗᵐ(θ,φ),

where (Qₙ) are radial FEM basis functions and Yₗᵐ are spherical harmonics.

Fields

lₕ: Maximum angular momentum quantum number.
nₕ: Number of orbitals per angular momentum channel.
Nₕ: Number of radial basis functions.
nspin: Number of spin components.
N: Total number of electrons.
Rmax: Radial cutoff.
basis: Radial FEM basis.
femops: Finite element operators.
ksham: Kohn–Sham Hamiltonian components.
cache: Preallocated workspaces.
fem_integration_method: Quadrature method for FEM integrals.

Use init_cache! to assemble all discretization-dependent operators.

source

lh (the angular momentum cutoff) should cover every occupied shell's l — e.g. lh = 1 for an atom occupying up to a p shell. nh bounds how many radial solutions per channel are kept.

Example

This continues the running example from The physical model (each manual page is self-contained, so we redefine model here too):

using AtomicKohnSham
using Libxc

model = KSEModel(Z = 11, N = 11, ex = Functional(:lda_x; n_spin = 1), ec = NoFunctional(1))

mesh = expmesh(0, 500, 60; s = 1.2)
basis = P1IntLegendreBasis(mesh; ordermax = 10)
discretization = KSEDiscretization(basis, model; lh = 2, nh = 10,
    fem_integration_method = GaussLegendre(basis, 2000))
discretization.Nₕ

Continue to Solving for the ground state to run the SCF loop.