GitHub

Solving for the ground state

The nonlinear eigenvalue problem is solved by a Self-Consistent Field (SCF) iteration: build an effective Hamiltonian from the current density, diagonalize it to get new orbitals, form the new density, and repeat until convergence.

SCF algorithm

The Optimal Damping Algorithm (ODA) is currently the supported SCF scheme. It mixes successive densities along an energy-minimizing line search rather than replacing them outright, which stabilizes convergence:

AtomicKohnSham.ODA — Type
ODA(; scftol, tinit, frozen_t=false, aufbau=OptimizedAufbau(),
     maxiter_ls=100, abstol_ls=…, reltol_ls=…)

Optimal Damping Algorithm (ODA) for self-consistent field (SCF) iterations.

ODA stabilizes SCF convergence by mixing successive densities using an energy-based line search along a linear interpolation path. At each iteration, the optimal mixing parameter is determined by minimizing a one-dimensional energy model including kinetic, Coulomb, and Hartree contributions, with an optional nonlinear correction (e.g. exchange–correlation).

The algorithm supports fractional occupations through a configurable Aufbau strategy and can optionally freeze the mixing parameter.

Keywords

  • scftol: Convergence tolerance for the SCF iteration.
  • tinit: Initial mixing (relaxation) parameter.
  • frozen_t: If true, the mixing parameter is kept fixed to tinit.
  • aufbau: Aufbau method used to fill orbital occupations.
  • maxiter_ls: Maximum number of iterations in the line search.
  • abstol_ls: Absolute tolerance for the line-search minimization.
  • reltol_ls: Relative tolerance for the line-search minimization.

ODA is designed to be robust in the presence of degeneracies and large density variations, and can be combined with acceleration techniques such as DIIS once the SCF iterations are stabilized.

source

Aufbau (occupation) scheme

ODA needs a strategy for turning orbital energies into occupation numbers at each iteration — the aufbau keyword:

AtomicKohnSham.OptimizedAufbau — Type
OptimizedAufbau(; max_degen=2, tol=1e-3, handle_degen_spin_1iter=true)

Aufbau occupation scheme with explicit treatment of low-dimensional degeneracies.

This scheme fills orbitals following the Aufbau principle and handles degenerate levels by solving a small variational problem to minimize the total energy when the last shell is only partially filled.

Parameters

  • max_degen::Int: Maximum degeneracy dimension handled explicitly (default: 2).
  • tol::Real: Energy tolerance to detect degeneracies.
  • handle_degen_spin_1iter::Bool: If true, spin-only degeneracies (same spatial orbital, opposite spins) are treated symmetrically at the first SCF iteration.
source

Two alternatives exist for special cases: FrozenAufbau(n::Dict) fixes the occupations to user-prescribed values (ignoring orbital energies entirely), and SmearedAufbau(Temp) is a thermal-smearing scheme.

Running the SCF loop

AtomicKohnSham.groundstate — Function
groundstate(model::KSEModel, discretization::KSEDiscretization, method::SCFAlgorithm;
           kwargs...) -> KSESolution

Compute the ground-state solution of a Kohn–Sham model using a given discretization and SCF algorithm.

Arguments

  • model::KSEModel: Kohn–Sham extended model.
  • discretization::KSEDiscretization: Discretization of the radial Kohn–Sham equations.
  • method::SCFAlgorithm: SCF algorithm.

Keyword arguments

  • maxiter::Int = 100: Maximum number of SCF iterations.
  • callback = nothing: Callback (or CallbackSet) invoked once per SCF iteration, e.g. a LogFileCallback for per-iteration diagnostics.

Returns

  • KSESolution: Ground-state solution including orbitals, energies, density, etc.
source

For finer control (inspecting the solver between iterations, custom callbacks — see Reports, logs & plots), the lower-level KSESolver/solve! pair that groundstate wraps is also part of the public API.

Example

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))

aufbau = OptimizedAufbau(max_degen = 2, tol = 1e-1)
alg = ODA(tinit = 0.6, aufbau = aufbau, scftol = 1e-9)

sol = groundstate(model, discretization, alg; maxiter = 100)
println(sol)
Name : Sodium
Success = SUCCESS
niter = 20
Stopping criteria = 5.569658855746759e-10
Occupation number = 
            1s : ε = -3.764758179762527e+01        n = 2.0
            2s : ε = -2.007737456192821e+00        n = 2.0
            2p : ε = -1.006028994917961e+00        n = 6.0
            3s : ε = -7.701608708314467e-02        n = 0.999999999991612
            3p : ε = -1.111637722765135e-02        n = 8.387956995647983e-12

Continue to Analyzing a solution to inspect sol in more detail.