C++ API

This section documents the C++ native library. The documentation is generated automatically from the source headers using Doxygen and Breathe.

Scheme Classes

class HF : public Logger
#include <hf.hpp>

Solver for the Hartree-Fock (HF) dielectric scheme.

Base class for all dielectric scheme solvers. Computes the ideal density response (IDR), the static structure factor (SSF), and the local field correction (LFC) on a wave-vector grid. Derived classes override the virtual compute methods to implement higher-level approximations.

Subclassed by Rpa

Public Functions

HF(const std::shared_ptr<const Input> &in_, const bool verbose_)

Construct with explicit verbosity flag.

Parameters
  • in_ – Shared pointer to the input parameters.

  • verbose_ – If false, solver output is suppressed.

inline explicit HF(const std::shared_ptr<const Input> &in_)

Construct with verbosity enabled.

Parameters

in_ – Shared pointer to the input parameters.

virtual ~HF() = default

Virtual destructor.

int compute()

Run the full solver pipeline.

Returns

0 on success.

inline const Vector2D &getIdr() const

Return the ideal density response grid (wave-vectors × Matsubara frequencies).

inline const Vector2D &getLfc() const

Return the local field correction grid (wave-vectors × Matsubara frequencies).

inline const std::vector<double> &getSsf() const

Return the static structure factor over the wave-vector grid.

inline const std::vector<double> &getWvg() const

Return the wave-vector grid.

inline double getChemicalPotential() const

Return the chemical potential.

Returns

Chemical potential (in units of the thermal energy).

std::vector<double> getSdr() const

Compute and return the static density response.

Returns

Vector of static density response values over the wave-vector grid.

double getUInt() const

Compute and return the interaction energy.

Returns

Dimensionless interaction energy per particle.

Protected Functions

inline const Input &in() const

Dereference the input shared pointer.

virtual void init()

Initialize the wave-vector grid and other basic derived quantities.

virtual void computeStructuralProperties()

Compute all structural properties (IDR, SSF, LFC).

virtual void computeSsf()

Dispatch the SSF calculation to the appropriate temperature branch.

virtual void computeSsfFinite()

Compute the static structure factor at finite temperature.

virtual void computeSsfGround()

Compute the static structure factor at zero temperature (ground state).

virtual void computeLfc()

Compute the local field correction (zero for bare HF).

Protected Attributes

const std::shared_ptr<const Input> inPtr

Shared pointer to the input parameters.

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

std::vector<double> wvg

Wave-vector grid.

Vector2D idr

Ideal density response (rows = wave-vectors, columns = Matsubara frequencies).

Vector2D lfc

Static local field correction (rows = wave-vectors, columns = Matsubara frequencies).

std::vector<double> ssf

Static structure factor over the wave-vector grid.

double mu

Chemical potential (in units of the thermal energy).

Private Functions

void buildWaveVectorGrid()

Construct the uniform wave-vector grid from the input parameters.

void computeChemicalPotential()

Compute the chemical potential from the normalization condition.

void computeIdr()

Dispatch the IDR calculation to the appropriate temperature branch.

void computeIdrFinite()

Compute the ideal density response at finite temperature.

void computeIdrGround()

Compute the ideal density response at zero temperature.

namespace HFUtil

Internal helpers for the Hartree-Fock ideal density response and SSF.

class Idr : public dimensionsUtil::DimensionsHandler
#include <hf.hpp>

Computes the ideal density response at finite temperature.

Evaluates the IDR for a given wave-vector x_ over all Matsubara frequencies by numerical integration over the auxiliary momentum.

Public Functions

inline Idr(const std::shared_ptr<const Input> in_, const double &x_, const double &mu_, const double &yMin_, const double &yMax_, std::shared_ptr<Integrator1D> itg_)

Construct for a finite-temperature IDR calculation.

Parameters
  • in_ – Shared pointer to the input parameters.

  • x_ – Wave-vector value.

  • mu_ – Chemical potential.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • itg_ – Shared pointer to a 1D integrator.

std::vector<double> get()

Compute and return the IDR for all Matsubara frequencies.

Returns

Vector of IDR values, one per Matsubara frequency.

Private Functions

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

double integrand(const double &y, const int &l) const

3D integrand for Matsubara frequency l.

Parameters
  • y – Auxiliary momentum variable.

  • l – Matsubara frequency index.

double integrand(const double &y) const

3D integrand for the zeroth Matsubara frequency.

Parameters

y – Auxiliary momentum variable.

double integrand2D(const double &y, const int &l) const

2D integrand for Matsubara frequency l.

Parameters
  • y – Auxiliary momentum variable.

  • l – Matsubara frequency index.

double integrand2D(const double &y) const

2D integrand for the zeroth Matsubara frequency.

Parameters

y – Auxiliary momentum variable.

Private Members

const std::shared_ptr<const Input> in

Input parameters.

const double x

Wave-vector.

const double mu

Chemical potential.

const double yMin

Lower integration limit.

const double yMax

Upper integration limit.

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

std::vector<double> res

Storage for the per-frequency integration results.

class IdrGround : public dimensionsUtil::DimensionsHandler
#include <hf.hpp>

Computes the ideal density response at zero temperature.

Public Functions

inline IdrGround(const dimensionsUtil::Dimension dim_, const double &x_, const double &Omega_)

Construct for a ground-state IDR calculation.

Parameters
  • dim_ – Dimensionality (2D or 3D).

  • x_ – Wave-vector value.

  • Omega_ – Frequency parameter.

double get()

Compute and return the ground-state IDR.

Returns

Analytic IDR value at (x_, Omega_).

Private Functions

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

Private Members

const dimensionsUtil::Dimension dim

Dimensionality.

const double x

Wave-vector.

const double Omega

Frequency parameter.

double res

Result of the ground-state IDR computation.

class Ssf : public dimensionsUtil::DimensionsHandler
#include <hf.hpp>

Computes the Hartree-Fock static structure factor at finite temperature.

Evaluates F_HF(x, 0) via numerical integration over the auxiliary momentum. For 3D systems, uses a single 1D integral. For 2D systems, uses a 2D integral (over y and angle p) plus the IDR contribution.

Public Functions

inline Ssf(const std::shared_ptr<const Input> in_, const double &x_, const double &mu_, const double &yMin_, const double &yMax_, std::shared_ptr<Integrator1D> itg_, const std::vector<double> &itgGrid_, std::shared_ptr<Integrator2D> itg2_, const double &idr0_)

Construct for a finite-temperature SSF calculation.

Parameters
  • in_ – Shared pointer to the input parameters.

  • x_ – Wave-vector value.

  • mu_ – Chemical potential.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • itg_ – Shared pointer to a 1D integrator.

  • itgGrid_ – Grid for 2D integration.

  • itg2_ – Shared pointer to a 2D integrator.

  • idr0_ – Ideal density response at l=0 for the current wave-vector.

double get()

Compute and return the HF static structure factor.

Returns

SSF value at the current wave-vector.

Private Functions

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

double integrand(const double &y) const

3D integrand over auxiliary momentum y.

Parameters

y – Auxiliary momentum variable.

double integrand2DOut(const double &y) const

Outer 2D integrand over auxiliary momentum y.

Parameters

y – Outer integration variable.

double integrand2DIn(const double &p) const

Inner 2D integrand (angle) over p.

Parameters

p – Inner integration variable (angle).

Private Members

const std::shared_ptr<const Input> in

Input parameters.

const double x

Wave-vector.

const double mu

Chemical potential.

const double yMin

Lower integration limit.

const double yMax

Upper integration limit.

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

const std::vector<double> &itgGrid

Grid for 2D integration.

const std::shared_ptr<Integrator2D> itg2

2D numerical integrator.

const double idr0

Ideal density response at l=0 for the current wave-vector.

double res

Result of the SSF computation.

class SsfGround : public dimensionsUtil::DimensionsHandler
#include <hf.hpp>

Computes the Hartree-Fock static structure factor at zero temperature.

Uses analytic ground-state expressions in both 2D and 3D.

Public Functions

inline SsfGround(const dimensionsUtil::Dimension dim_, const double &x_)

Construct for a zero-temperature SSF calculation.

Parameters
  • dim_ – Dimensionality (2D or 3D).

  • x_ – Wave-vector value.

double get()

Compute and return the ground-state HF static structure factor.

Returns

Ground-state SSF value at x_.

Private Functions

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

Private Members

const dimensionsUtil::Dimension dim

Dimensionality.

const double x

Wave-vector.

double res

Result of the ground-state SSF computation.

class Rpa : public HF
#include <rpa.hpp>

Solver for the Random Phase Approximation (RPA) dielectric scheme.

Extends the Hartree-Fock (HF) base class to compute the static structure factor (SSF) within the RPA. Both finite-temperature and zero-temperature (ground state) regimes are supported. Helper classes for the SSF computation are provided in the RpaUtil namespace.

Subclassed by ESA, Stls

Public Functions

Rpa(const std::shared_ptr<const Input> &in_, const bool verbose_)

Construct with explicit verbosity flag.

Parameters
  • in_ – Shared pointer to the input parameters.

  • verbose_ – If false, solver output is suppressed.

inline explicit Rpa(const std::shared_ptr<const Input> &in_)

Construct with verbosity enabled.

Parameters

in_ – Shared pointer to the input parameters.

Protected Functions

virtual void init() override

Initialize wave-vector grid and other basic properties.

virtual void computeSsfFinite() override

Compute the static structure factor at finite temperature.

virtual void computeSsfGround() override

Compute the static structure factor at zero temperature.

Protected Attributes

std::vector<double> ssfHF

Hartree-Fock static structure factor.

Private Functions

void computeSsfHF()

Compute the Hartree-Fock static structure factor.

virtual void computeLfc() override

Compute the local field correction (zero for RPA).

namespace RpaUtil

Internal helpers for the RPA static structure factor computation.

class SsfBase
#include <rpa.hpp>

Base class holding shared state for SSF helpers.

Subclassed by RpaUtil::Ssf, RpaUtil::SsfGround

Protected Functions

inline SsfBase(const double &x_, const double &ssfHF_, std::span<const double> lfc_, const std::shared_ptr<const Input> in_)

Construct the base with the quantities needed for SSF evaluation.

Parameters
  • x_ – Wave-vector value.

  • ssfHF_ – Hartree-Fock SSF at this wave-vector.

  • lfc_ – Span over the local field correction array.

  • in_ – Shared pointer to the input parameters.

double ip() const

Normalized interaction potential at the current wave-vector.

Protected Attributes

const double x

Wave-vector value.

const double ssfHF

Hartree-Fock SSF.

std::span<const double> lfc

Local field correction values.

const std::shared_ptr<const Input> in

Input parameters.

class Ssf : public RpaUtil::SsfBase, private dimensionsUtil::DimensionsHandler
#include <rpa.hpp>

Computes the finite-temperature RPA static structure factor.

Integrates over Matsubara frequencies using the ideal density response.

Public Functions

inline Ssf(const double &x_, const double &ssfHF_, std::span<const double> lfc_, const std::shared_ptr<const Input> in_, std::span<const double> idr_)

Construct for a finite-temperature calculation.

Parameters
  • x_ – Wave-vector value.

  • ssfHF_ – Hartree-Fock static structure factor at this wave-vector.

  • lfc_ – Span over the local field correction array.

  • in_ – Shared pointer to the input parameters.

  • idr_ – Span over the ideal density response array.

double get()

Compute and return the static structure factor.

Protected Attributes

const std::span<const double> idr

Ideal density response values over Matsubara frequencies.

Private Functions

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

double computeMatsubaraSummation() const

Compute the finite-temperature RPA SSF from Matsubara terms.

Returns

Static structure factor at the current wave-vector.

Private Members

double res

Stores the result of the frequency summation.

class SsfGround : public RpaUtil::SsfBase
#include <rpa.hpp>

Computes the zero-temperature RPA static structure factor.

Uses a 1D numerical integrator over real frequencies.

Subclassed by QstlsUtil::SsfGround

Public Functions

inline SsfGround(const double &x_, const double &ssfHF_, std::span<const double> lfc_, std::shared_ptr<Integrator1D> itg_, const std::shared_ptr<const Input> in_)

Construct for a zero-temperature calculation.

Parameters
  • x_ – Wave-vector value.

  • ssfHF_ – Hartree-Fock static structure factor at this wave-vector.

  • lfc_ – Span over the local field correction array.

  • itg_ – Shared pointer to a 1D integrator instance.

  • in_ – Shared pointer to the input parameters.

double get()

Compute and return the static structure factor.

Protected Functions

double integrand(const double &Omega) const

Integrand for the zero-temperature frequency integral.

Parameters

Omega – Real frequency value.

Returns

Value of the integrand at Omega.

Protected Attributes

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

class Stls : public Rpa
#include <stls.hpp>

Solver for the STLS (Singwi–Tosi–Land–Sjölander) dielectric scheme.

Extends the RPA solver by computing the static local field correction (SLFC) self-consistently via an iterative scheme. The SLFC accounts for short-range correlations beyond the random phase approximation.

Subclassed by Qstls, StlsIet, VSStlsManager, VSStlsWorker

Public Functions

Stls(const std::shared_ptr<const StlsInput> &in_, const bool verbose_)

Construct with explicit verbosity flag.

Parameters
  • in_ – Shared pointer to the STLS input parameters.

  • verbose_ – If false, solver output is suppressed.

inline explicit Stls(const std::shared_ptr<const StlsInput> &in_)

Construct with verbosity enabled.

Parameters

in_ – Shared pointer to the STLS input parameters.

~Stls() override = default

Virtual destructor.

inline double getError() const

Return the convergence error from the last iteration.

Returns

RMS difference between successive SSF iterates.

Protected Functions

virtual void computeStructuralProperties() override

Compute the SSF and SLFC self-consistently.

virtual void computeSsf() override

Compute the static structure factor using the current SLFC.

virtual void computeLfc() override

Compute the static local field correction from the current SSF.

virtual void initialGuess()

Set the initial guess for the iterative procedure.

virtual bool initialGuessFromInput()

Attempt to load the initial guess from the input parameters.

Returns

True if a valid guess was found in the input, false otherwise.

virtual double computeError() const

Compute the convergence error between successive iterates.

Returns

RMS difference between ssf and ssfOld.

virtual void updateSolution()

Mix the new and old SSF to advance the iteration.

Protected Attributes

std::vector<double> ssfOld

SSF from the previous iteration (used for mixing).

Private Functions

const StlsInput &in() const

Access the input as an StlsInput reference.

namespace StlsUtil

Internal helpers for the STLS static local field correction.

Functions

template<typename T>
T check_dynamic_cast_result(T ptr)

Perform a checked dynamic pointer cast between shared_ptr types.

Throws an MPI error if the cast result is null (i.e., the runtime type does not match TOut).

Template Parameters

T – Resulting shared_ptr type.

Parameters

ptr – Pointer to validate.

Returns

ptr unchanged.

template<typename TIn, typename TOut>
std::shared_ptr<const TOut> dynamic_pointer_cast(const std::shared_ptr<const TIn> &in)

Downcast a shared_ptr<const TIn> to shared_ptr<const TOut>.

Template Parameters
  • TIn – Source type.

  • TOut – Target type (must be derived from or base of TIn).

Parameters

in – Source pointer.

Returns

Downcast pointer; throws on failure.

template<typename TIn, typename TOut>
std::shared_ptr<TOut> dynamic_pointer_cast(const std::shared_ptr<TIn> &in)

Downcast a shared_ptr<TIn> to shared_ptr<TOut>.

Template Parameters
  • TIn – Source type.

  • TOut – Target type (must be derived from or base of TIn).

Parameters

in – Source pointer.

Returns

Downcast pointer; throws on failure.

class SlfcBase
#include <stls.hpp>

Base class holding shared state for SLFC helpers.

Subclassed by StlsIetUtil::Slfc, StlsUtil::Slfc

Protected Functions

inline SlfcBase(const double &x_, const double &yMin_, const double &yMax_, std::shared_ptr<Interpolator1D> ssfi_)

Construct with the quantities needed for SLFC evaluation.

Parameters
  • x_ – Wave-vector value.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • ssfi_ – Shared pointer to an SSF interpolator.

double ssf(const double &y) const

Evaluate the interpolated SSF at auxiliary momentum y.

Parameters

y – Auxiliary momentum value.

Returns

Interpolated SSF value.

Protected Attributes

const double x

Wave-vector value.

const double yMin

Lower integration limit.

const double yMax

Upper integration limit.

const std::shared_ptr<Interpolator1D> ssfi

Interpolator for the static structure factor.

class Slfc : public StlsUtil::SlfcBase, private dimensionsUtil::DimensionsHandler
#include <stls.hpp>

Computes the STLS static local field correction.

Evaluates the SLFC integral at a single wave-vector x_ by integrating over the auxiliary momentum with the SSF interpolator.

Public Functions

inline Slfc(const double &x_, const double &yMin_, const double &yMax_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Integrator1D> itg_, const std::shared_ptr<const Input> in_)

Construct for a SLFC calculation.

Parameters
  • x_ – Wave-vector value.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • itg_ – Shared pointer to a 1D integrator.

  • in_ – Shared pointer to the input parameters.

double get()

Compute and return the SLFC at the current wave-vector.

Returns

SLFC value.

Private Functions

double integrand(const double &y) const

3D integrand over auxiliary momentum y.

Parameters

y – Auxiliary momentum variable.

double integrand2D(const double &y) const

2D integrand over auxiliary momentum y.

Parameters

y – Auxiliary momentum variable.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

Private Members

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

const std::shared_ptr<const Input> in

Input parameters.

double res

Stores the result of the integration.

class ESA : public Rpa
#include <esa.hpp>

Solver for the ESA (Enhanced Screening Approximation) dielectric scheme.

Extends the RPA solver by replacing the zero local field correction with a neural-network parametrization of the static local field correction (SLFC). The parametrization satisfies the compressibility sum rule and uses machine-learned coefficients fitted to quantum Monte Carlo data.

Public Functions

inline explicit ESA(const std::shared_ptr<const Input> &in_)

Construct the ESA solver.

Parameters

in_ – Shared pointer to the input parameters.

Private Functions

virtual void computeLfc() override

Compute the neural-network-parametrized local field correction.

namespace ESAUtil

Internal helpers for the ESA static local field correction.

class Slfc
#include <esa.hpp>

Computes the ESA static local field correction for a given state point.

The SLFC is constructed from a neural network parametrization that blends the compressibility-sum-rule (CSR) result at long wavelengths with a machine-learned fit at shorter wavelengths via a smooth activation function.

Public Functions

inline explicit Slfc(const double &rs_, const double &theta_)

Construct for a given coupling and degeneracy parameter.

Parameters
  • rs_ – Coupling parameter (Wigner–Seitz radius).

  • theta_ – Degeneracy parameter (reduced temperature).

double get(const double &x)

Evaluate the SLFC at wave-vector x.

Parameters

x – Dimensionless wave-vector.

Returns

SLFC value at x.

double nn(const double &x) const

Neural-network parametrization of the SLFC.

Parameters

x – Dimensionless wave-vector.

Returns

Neural-network contribution to the SLFC at x.

double csr(const double &x) const

Compressibility-sum-rule contribution to the SLFC.

Parameters

x – Dimensionless wave-vector.

Returns

CSR contribution at x.

void computeCoefficients()

Compute and cache all SLFC coefficients.

void computeNNCoefficients()

Compute and cache the neural-network coefficients.

void computeCSRCoefficients()

Compute and cache the compressibility-sum-rule coefficient.

double onTop() const

On-top value of the radial distribution function.

Returns

g(0) at the current state point.

double activationFunction(const double &x) const

Activation function that blends the long- and short-wavelength limits.

Parameters

x – Dimensionless wave-vector.

Returns

Value of the activation function at x.

Dual22 freeEnergy() const

Evaluate the parametrized free energy (value only).

Returns

Free energy dual number at the current (rs, theta).

Dual22 freeEnergy(const Dual22 &rs, const Dual22 &theta) const

Evaluate the parametrized free energy at arbitrary (rs, theta).

Parameters
  • rs – Dual-number coupling parameter.

  • theta – Dual-number degeneracy parameter.

Returns

Free energy dual number carrying value and derivatives.

Public Members

const double rs

Coupling parameter.

const double theta

Degeneracy parameter.

Coefficients coeff

Cached coefficient set for the current state point.

struct Coefficients
#include <esa.hpp>

Cache of pre-computed SLFC expansion coefficients.

The coefficients are computed once and reused for all wave-vector evaluations at the same state point.

Public Members

bool valid = false

True if the coefficients are up to date and valid.

double lwl

Long-wavelength limit coefficient.

double afEta

Activation-function width parameter.

double afxm

Activation-function inflection point.

double nna

Neural-network coefficient a.

double nnb

Neural-network coefficient b.

double nnc

Neural-network coefficient c.

double nnd

Neural-network coefficient d.

double csr

Compressibility-sum-rule coefficient.

class Qstls : public Stls
#include <qstls.hpp>

Solver for the qSTLS (quantum STLS) dielectric scheme.

Extends the STLS solver by replacing the classical SLFC with a quantum local field correction derived from the auxiliary density response (ADR). The ADR is split into a “fixed” frequency-independent component (stored in an SQLite database) and a frequency-dependent part that is iterated to self-consistency together with the SSF.

Subclassed by QstlsIet, VSQstlsManager, VSQstlsWorker

Public Functions

Qstls(const std::shared_ptr<const QstlsInput> &in_, const bool verbose_)

Construct with explicit verbosity flag.

Parameters
  • in_ – Shared pointer to the qSTLS input parameters.

  • verbose_ – If false, solver output is suppressed.

inline explicit Qstls(const std::shared_ptr<const QstlsInput> &in_)

Construct with verbosity enabled.

Parameters

in_ – Shared pointer to the qSTLS input parameters.

inline const Vector3D &getAdrFixed() const

Return the fixed component of the auxiliary density response.

Returns

3D array indexed by (wave-vector, wave-vector, Matsubara frequency).

Protected Functions

virtual void init() override

Initialize wave-vector grid, chemical potential, and ADR fixed component.

virtual void computeLfc() override

Compute the quantum local field correction from the ADR.

void readAdrFixed(Vector3D &res, const std::string &name, int runId) const

Read the fixed ADR component from the database.

Parameters
  • res – Output array to populate.

  • name – Database file path.

  • runId – Run identifier in the database.

void writeAdrFixed(const Vector3D &res, const std::string &name) const

Write the fixed ADR component to the database.

Parameters
  • res – Array to persist.

  • name – Database file path.

Protected Attributes

Vector3D adrFixed

Fixed (frequency-independent) component of the ADR.

std::string adrFixedDatabaseName

Path to the SQLite database used to persist adrFixed.

Private Functions

inline const QstlsInput &in() const

Access the input as a QstlsInput reference.

void computeAdrFixed()

Compute and cache the fixed component of the ADR.

virtual void computeSsfGround() override

Compute the zero-temperature static structure factor.

namespace QstlsUtil

Internal helpers for the qSTLS auxiliary density response.

Functions

void deleteBlobDataOnDisk(const std::string &dbName, int runId)

Delete blob data associated with a specific run from disk.

Parameters
  • dbName – Path to the SQLite database file.

  • runId – Run identifier whose blobs should be deleted.

class AdrBase
#include <qstls.hpp>

Base class holding shared state for ADR helper classes.

Subclassed by QstlsIetUtil::AdrIet, QstlsUtil::Adr, QstlsUtil::AdrGround

Public Functions

inline AdrBase(const double &Theta_, const double &yMin_, const double &yMax_, const double &x_, std::shared_ptr<Interpolator1D> ssfi_)

Construct with thermodynamic and grid parameters.

Parameters
  • Theta_ – Degeneracy parameter.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • x_ – Wave-vector value.

  • ssfi_ – Shared pointer to an SSF interpolator.

Protected Functions

double ssf(const double &y) const

Evaluate the interpolated SSF at auxiliary momentum y.

Parameters

y – Auxiliary momentum.

Returns

Interpolated SSF value.

Protected Attributes

const double Theta

Degeneracy parameter.

const double yMin

Lower integration limit.

const double yMax

Upper integration limit.

const double x

Wave-vector.

const std::shared_ptr<Interpolator1D> ssfi

Interpolator for the static structure factor.

const double isc

Integrand scaling constant.

const double isc0

Integrand scaling constant for the zeroth Matsubara frequency.

class AdrFixedBase
#include <qstls.hpp>

Base class holding shared state for fixed-ADR helper classes.

Subclassed by QstlsIetUtil::AdrFixedIet, QstlsUtil::AdrFixed

Public Functions

inline AdrFixedBase(const double &Theta_, const double &qMin_, const double &qMax_, const double &x_, const double &mu_)

Construct with thermodynamic and grid parameters.

Parameters
  • Theta_ – Degeneracy parameter.

  • qMin_ – Lower integration limit over the auxiliary momentum.

  • qMax_ – Upper integration limit over the auxiliary momentum.

  • x_ – Wave-vector value.

  • mu_ – Chemical potential.

Protected Attributes

const double Theta

Degeneracy parameter.

const double qMin

Lower integration limit.

const double qMax

Upper integration limit.

const double x

Wave-vector.

const double mu

Chemical potential.

class Adr : public QstlsUtil::AdrBase
#include <qstls.hpp>

Computes the finite-temperature auxiliary density response.

Evaluates the ADR at wave-vector x_ and all Matsubara frequencies by integrating over the auxiliary momentum using the fixed ADR component and the current SSF.

Public Functions

inline Adr(const double &Theta_, const double &yMin_, const double &yMax_, const double &x_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Integrator1D> itg_)

Construct for a finite-temperature ADR calculation.

Parameters
  • Theta_ – Degeneracy parameter.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • x_ – Wave-vector value.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • itg_ – Shared pointer to a 1D integrator.

void get(const std::vector<double> &wvg, const Vector3D &fixed, Vector2D &res)

Compute and store the ADR for all Matsubara frequencies.

Parameters
  • wvg – Wave-vector grid.

  • fixed – Fixed ADR component (wave-vectors × wave-vectors × frequencies).

  • res – Output array (wave-vectors × frequencies) to accumulate into.

Private Functions

double fix(const double &y) const

Evaluate the fixed-ADR contribution at auxiliary momentum y.

Parameters

y – Auxiliary momentum.

Returns

Fixed-component value at y.

double integrand(const double &y) const

Integrand for the ADR frequency integral.

Parameters

y – Auxiliary momentum.

Returns

Integrand value at y.

Private Members

Interpolator1D fixi

Interpolator for the fixed ADR component.

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

class AdrFixed : public QstlsUtil::AdrFixedBase
#include <qstls.hpp>

Computes the fixed (frequency-independent) component of the ADR.

This component depends only on the Fermi-Dirac occupation numbers and does not change during the self-consistency iterations.

Public Functions

inline AdrFixed(const double &Theta_, const double &qMin_, const double &qMax_, const double &x_, const double &mu_, const std::vector<double> &itgGrid_, std::shared_ptr<Integrator2D> itg_)

Construct for a finite-temperature fixed-ADR calculation.

Parameters
  • Theta_ – Degeneracy parameter.

  • qMin_ – Lower integration limit.

  • qMax_ – Upper integration limit.

  • x_ – Wave-vector value.

  • mu_ – Chemical potential.

  • itgGrid_ – Grid for the 2D integration.

  • itg_ – Shared pointer to a 2D integrator.

void get(const std::vector<double> &wvg, Vector3D &res) const

Compute and store the fixed ADR component for all wave-vector pairs.

Parameters
  • wvg – Wave-vector grid.

  • res – Output 3D array (wave-vectors × wave-vectors × frequencies).

Private Functions

double integrand1(const double &q, const double &l) const

First integrand over auxiliary momentum q and frequency l.

Parameters
  • q – Auxiliary momentum.

  • l – Matsubara frequency value.

double integrand2(const double &t, const double &y, const double &l) const

Second integrand over variables t, y, and l.

Parameters
  • t – Intermediate momentum variable.

  • y – Auxiliary momentum variable.

  • l – Matsubara frequency value.

Private Members

const std::shared_ptr<Integrator2D> itg

2D numerical integrator.

const std::vector<double> &itgGrid

Grid for 2D integration.

class AdrGround : public QstlsUtil::AdrBase
#include <qstls.hpp>

Computes the auxiliary density response at zero temperature.

Evaluates the ground-state ADR by integrating over the real frequency and the auxiliary momentum.

Public Functions

inline AdrGround(const double &x_, const double &Omega_, std::shared_ptr<Interpolator1D> ssfi_, const double &yMax_, std::shared_ptr<Integrator2D> itg_)

Construct for a zero-temperature ADR calculation.

Parameters
  • x_ – Wave-vector value.

  • Omega_ – Real frequency value.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • yMax_ – Upper integration limit.

  • itg_ – Shared pointer to a 2D integrator.

double get()

Compute and return the ground-state ADR.

Returns

ADR value at (x_, Omega_).

Private Functions

double integrand1(const double &y) const

Outer integrand over auxiliary momentum y.

Parameters

y – Auxiliary momentum.

double integrand2(const double &t) const

Inner integrand over intermediate variable t.

Parameters

t – Intermediate integration variable.

Private Members

const double Omega

Real frequency.

const std::shared_ptr<Integrator2D> itg

2D numerical integrator.

class SsfGround : public RpaUtil::SsfGround
#include <qstls.hpp>

Computes the zero-temperature qSTLS static structure factor.

Extends the RPA ground-state SSF to include the quantum local field correction derived from the ADR.

Public Functions

inline SsfGround(const double &x_, const double &ssfHF_, const double &xMax_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Integrator1D> itg_, const std::shared_ptr<const Input> in_)

Construct for a zero-temperature qSTLS SSF calculation.

Parameters
  • x_ – Wave-vector value.

  • ssfHF_ – Hartree-Fock SSF at x_.

  • xMax_ – Upper integration limit for the wave-vector integral.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • itg_ – Shared pointer to a 1D integrator.

  • in_ – Shared pointer to the input parameters.

double get()

Compute and return the zero-temperature qSTLS SSF.

Returns

SSF value at x_.

Private Functions

double integrand(const double &Omega) const

Integrand for the zero-temperature frequency integral.

Parameters

Omega – Real frequency value.

Returns

Integrand value at Omega.

Private Members

const double xMax

Upper integration limit for the wave-vector integral.

const std::shared_ptr<Interpolator1D> ssfi

SSF interpolator.

class Iet : public Logger
#include <iet.hpp>

Computes the bridge function for IET-based dielectric schemes.

The IET (Integral Equation Theory) extension adds a bridge function to the STLS local field correction. Supported bridge-function theories are HNC (hypernetted-chain), IOI (Ichimaru–Ogata–Ichimaru), and LCT (Lucco Castello–Tolias). The quantum-to-classical mapping for the effective coupling parameter is selected via the input mapping string.

Public Functions

inline explicit Iet(const std::shared_ptr<const IetInput> &in_, const std::vector<double> &wvg_)

Construct the IET object.

Parameters
  • in_ – Shared pointer to the IET input parameters.

  • wvg_ – Wave-vector grid inherited from the parent scheme.

void init()

Compute the bridge function over the wave-vector grid.

inline const std::vector<double> &getBf() const

Return the bridge function values over the wave-vector grid.

bool initialGuessFromInput(Vector2D &lfc)

Attempt to load the initial LFC guess from the input parameters.

Parameters

lfc – Output LFC array to populate if a guess is available.

Returns

True if a valid guess was found, false otherwise.

Private Functions

inline const IetInput &in() const

Access the input as an IetInput reference.

inline const IterationInput &inRpa() const

Cast the input to an IterationInput reference for RPA parameters.

void computeBf()

Compute the bridge function for each wave-vector in wvg.

Private Members

const std::shared_ptr<const IetInput> inPtr

Shared pointer to the IET input parameters.

const std::vector<double> wvg

Wave-vector grid (shared with the parent scheme).

std::vector<double> bf

Bridge function values over the wave-vector grid.

namespace IetUtil

Internal helpers for the IET bridge function calculation.

class BridgeFunction
#include <iet.hpp>

Computes the bridge function at a single wave-vector.

Supports three bridge-function theories (HNC, IOI, LCT) and three quantum-to-classical mappings (coupling-parameter, energy, or pressure-based). The effective classical coupling parameter is computed first, then the corresponding bridge function is evaluated.

Public Functions

inline BridgeFunction(const std::string &theory_, const std::string &mapping_, const double &rs_, const double &Theta_, const double &x_, std::shared_ptr<Integrator1D> itg_)

Construct for a single wave-vector evaluation.

Parameters
  • theory_ – Name of the bridge-function theory (“HNC”, “IOI”, “LCT”).

  • mapping_ – Name of the quantum-classical mapping (“coupling”, …).

  • rs_ – Quantum coupling parameter (Wigner–Seitz radius).

  • Theta_ – Degeneracy parameter (reduced temperature).

  • x_ – Dimensionless wave-vector.

  • itg_ – Shared pointer to a 1D integrator (used by LCT).

double get() const

Compute and return the bridge function.

Returns

Bridge-function value at x_.

Private Functions

double hnc() const

Hypernetted-chain (HNC) bridge function.

Returns

HNC bridge-function value (zero by definition).

double ioi() const

Ichimaru–Ogata–Ichimaru (IOI) bridge function.

Returns

IOI bridge-function value at x_.

double lct() const

Lucco Castello–Tolias (LCT) bridge function.

Returns

LCT bridge-function value at x_.

double lctIntegrand(const double &r, const double &Gamma) const

LCT integrand as a function of real-space distance r and classical coupling parameter Gamma.

Parameters
  • r – Real-space distance.

  • Gamma – Classical coupling parameter.

Returns

Integrand value.

double couplingParameter() const

Compute the effective classical coupling parameter for the bridge function.

Returns

Classical coupling parameter \(\Gamma\).

Private Members

const std::string theory

Bridge-function theory identifier.

const std::string mapping

Quantum-classical mapping identifier.

const double rs

Quantum coupling parameter.

const double Theta

Degeneracy parameter.

const double x

Wave-vector.

const std::shared_ptr<Integrator1D> itg

1D numerical integrator (used by the LCT theory).

const double lambda = pow(4.0 / (9.0 * M_PI), 1.0 / 3.0)

Unit-conversion constant \(\lambda = (4/9\pi)^{1/3}\).

class StlsIet : public Stls
#include <stlsiet.hpp>

Solver for the STLS-IET dielectric scheme.

Combines the STLS iterative solver with an IET bridge function. The static local field correction includes both the STLS exchange-correlation term and a bridge-function contribution computed by the Iet helper. Convergence is driven by the same iterative mixing as for plain STLS.

Public Functions

explicit StlsIet(const std::shared_ptr<const StlsIetInput> &in_)

Construct the STLS-IET solver.

Parameters

in_ – Shared pointer to the STLS-IET input parameters.

inline const std::vector<double> &getBf() const

Return the bridge function values over the wave-vector grid.

Returns

Constant reference to the bridge function array.

Private Functions

inline const StlsIetInput &in() const

Access the input as a StlsIetInput reference.

virtual void init() override

Initialize wave-vector grid, chemical potential, and bridge function.

virtual void computeLfc() override

Compute the STLS-IET static local field correction.

virtual bool initialGuessFromInput() override

Attempt to load the initial LFC guess from the input parameters.

Returns

True if a valid guess was found, false otherwise.

Private Members

Iet iet

IET helper that computes the bridge function.

const std::shared_ptr<Integrator2D> itg2D

2D integrator for the SLFC double integral.

std::vector<double> itgGrid

Quadrature grid for the 2D integral.

namespace StlsIetUtil

Internal helpers for the STLS-IET static local field correction.

class Slfc : public StlsUtil::SlfcBase, private dimensionsUtil::DimensionsHandler
#include <stlsiet.hpp>

Computes the STLS-IET static local field correction.

Evaluates the SLFC integral that combines the standard STLS exchange- correlation term with a bridge-function contribution via a 2D quadrature.

Public Functions

inline Slfc(const double &x_, const double &yMin_, const double &yMax_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Interpolator1D> lfci_, std::shared_ptr<Interpolator1D> bfi_, const std::vector<double> &itgGrid_, std::shared_ptr<Integrator2D> itg_, const std::shared_ptr<const Input> in_)

Construct for a STLS-IET SLFC calculation.

Parameters
  • x_ – Wave-vector value.

  • yMin_ – Lower integration limit.

  • yMax_ – Upper integration limit.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • lfci_ – Shared pointer to an LFC interpolator.

  • bfi_ – Shared pointer to a bridge-function interpolator.

  • itgGrid_ – Grid for 2D integration.

  • itg_ – Shared pointer to a 2D integrator.

  • in_ – Shared pointer to the input parameters.

double get()

Compute and return the SLFC at the current wave-vector.

Returns

SLFC value including the bridge-function contribution.

Private Functions

double integrand1(const double &y) const

3D outer integrand over auxiliary momentum y.

Parameters

y – Outer integration variable.

double integrand2(const double &w) const

3D inner integrand over auxiliary momentum w.

Parameters

w – Inner integration variable.

double integrand1_2D(const double &y) const

2D outer integrand over auxiliary momentum y.

Parameters

y – Outer integration variable.

double integrand2_2D(const double &w) const

2D inner integrand over auxiliary momentum w.

Parameters

w – Inner integration variable.

double lfc(const double &x) const

Evaluate the interpolated LFC at wave-vector x.

Parameters

x – Wave-vector value.

Returns

Interpolated LFC value.

double bf(const double &x_) const

Evaluate the interpolated bridge function at wave-vector x_.

Parameters

x_ – Wave-vector value.

Returns

Interpolated bridge function value.

virtual void compute2D() override

Compute for a 2D system. Implemented by derived classes.

virtual void compute3D() override

Compute for a 3D system. Implemented by derived classes.

Private Members

const std::shared_ptr<Integrator2D> itg

2D numerical integrator.

const std::vector<double> itgGrid

Grid for 2D integration.

const std::shared_ptr<Interpolator1D> lfci

Interpolator for the static local field correction.

const std::shared_ptr<Interpolator1D> bfi

Interpolator for the bridge function.

double res

Result of the SLFC integration.

const std::shared_ptr<const Input> in

Input parameters.

class QstlsIet : public Qstls
#include <qstlsiet.hpp>

Solver for the qSTLS-IET (quantum STLS with Integral Equation Theory) scheme.

Extends the qSTLS solver by incorporating an IET bridge-function contribution into the local field correction. Both the quantum (ADR-derived) and the classical (bridge-function) parts of the LFC are iterated to self-consistency.

Public Functions

explicit QstlsIet(const std::shared_ptr<const QstlsIetInput> &in_)

Construct the qSTLS-IET solver.

Parameters

in_ – Shared pointer to the qSTLS-IET input parameters.

inline const std::vector<double> &getBf() const

Return the bridge function values over the wave-vector grid.

Returns

Constant reference to the bridge function array.

Private Functions

inline const QstlsIetInput &in() const

Access the input as a QstlsIetInput reference.

virtual void init() override

Initialize wave-vector grid, chemical potential, bridge function, and fixed ADR.

virtual void computeLfc() override

Compute the qSTLS-IET local field correction.

void computeAdrFixed()

Compute the fixed (frequency-independent) ADR component.

virtual bool initialGuessFromInput() override

Attempt to load the initial LFC guess from the input parameters.

Returns

True if a valid guess was found, false otherwise.

Private Members

Iet iet

IET helper that computes the bridge function.

Vector2D lfcIet

IET contribution to the local field correction.

std::vector<double> itgGrid

Quadrature grid for the 2D ADR integral.

namespace QstlsIetUtil

Internal helpers for the qSTLS-IET auxiliary density response.

class AdrIet : public QstlsUtil::AdrBase
#include <qstlsiet.hpp>

Computes the IET contribution to the finite-temperature ADR.

Extends the standard qSTLS ADR by adding bridge-function and dynamic LFC terms. The integration combines the fixed ADR component with the current SSF, LFC, and bridge function.

Public Functions

inline AdrIet(const double &Theta_, const double &qMin_, const double &qMax_, const double &x_, std::shared_ptr<Interpolator1D> ssfi_, std::vector<std::shared_ptr<Interpolator1D>> lfci_, std::shared_ptr<Interpolator1D> bfi_, const std::vector<double> &itgGrid_, std::shared_ptr<Integrator2D> itg_)

Construct for a finite-temperature IET-ADR calculation.

Parameters
  • Theta_ – Degeneracy parameter.

  • qMin_ – Lower integration limit.

  • qMax_ – Upper integration limit.

  • x_ – Wave-vector value.

  • ssfi_ – Shared pointer to an SSF interpolator.

  • lfci_ – Per-frequency LFC interpolators (one per Matsubara frequency).

  • bfi_ – Shared pointer to a bridge-function interpolator.

  • itgGrid_ – Grid for 2D integration.

  • itg_ – Shared pointer to a 2D integrator.

void get(const std::vector<double> &wvg, const Vector3D &fixed, Vector2D &res)

Compute and store the IET-ADR for all Matsubara frequencies.

Parameters
  • wvg – Wave-vector grid.

  • fixed – Fixed ADR component (wave-vectors × wave-vectors × frequencies).

  • res – Output array (wave-vectors × frequencies) to accumulate into.

Private Functions

double integrand1(const double &q, const int &l) const

Outer integrand over auxiliary momentum q and Matsubara index l.

Parameters
  • q – Auxiliary momentum.

  • l – Matsubara frequency index.

double integrand2(const double &y) const

Inner integrand over auxiliary momentum y.

Parameters

y – Auxiliary momentum.

double lfc(const double &y, const int &l) const

Evaluate the dynamic LFC at wave-vector y and frequency index l.

Parameters
  • y – Auxiliary momentum.

  • l – Matsubara frequency index.

double bf(const double &y) const

Evaluate the bridge-function contribution at wave-vector y.

Parameters

y – Wave-vector value.

double fix(const double &x, const double &y) const

Evaluate the fixed ADR component at (x, y).

Parameters
  • x – Primary wave-vector.

  • y – Auxiliary wave-vector.

Private Members

const double &qMin = yMin

Reference to lower integration limit (alias for yMin).

const double &qMax = yMax

Reference to upper integration limit (alias for yMax).

const std::shared_ptr<Integrator2D> itg

2D numerical integrator.

const std::vector<double> &itgGrid

Grid for 2D integration.

const std::vector<std::shared_ptr<Interpolator1D>> lfci

Per-frequency LFC interpolators.

const std::shared_ptr<Interpolator1D> bfi

Bridge function interpolator.

Interpolator2D fixi

2D interpolator for the fixed ADR component.

class AdrFixedIet : public QstlsUtil::AdrFixedBase
#include <qstlsiet.hpp>

Computes the IET fixed component of the ADR.

This component depends only on Fermi–Dirac occupation numbers and the Matsubara frequency index, not on the self-consistent SSF or LFC.

Public Functions

inline AdrFixedIet(const double &Theta_, const double &qMin_, const double &qMax_, const double &x_, const double &mu_, std::shared_ptr<Integrator1D> itg_)

Construct for a finite-temperature fixed-IET-ADR calculation.

Parameters
  • Theta_ – Degeneracy parameter.

  • qMin_ – Lower integration limit.

  • qMax_ – Upper integration limit.

  • x_ – Wave-vector value.

  • mu_ – Chemical potential.

  • itg_ – Shared pointer to a 1D integrator.

void get(int l, const std::vector<double> &wvg, Vector3D &res) const

Compute and store the fixed IET-ADR for a single Matsubara frequency.

Parameters
  • l – Matsubara frequency index.

  • wvg – Wave-vector grid.

  • res – Output 3D array to accumulate into.

Private Functions

double integrand(const double &t, const double &y, const double &q, const double &l) const

Integrand over variables t, y, q, and l.

Parameters
  • t – Intermediate momentum variable.

  • y – Auxiliary momentum variable.

  • q – Auxiliary momentum variable.

  • l – Matsubara frequency value.

Private Members

const double &tMin = qMin

Reference to lower integration limit (alias for qMin).

const double &tMax = qMax

Reference to upper integration limit (alias for qMax).

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

class VSStlsWorker : public VSWorker, public Stls
#include <vsstls.hpp>

STLS-based worker for a single grid point in the VS scheme.

Wraps a Stls solver to provide the VSWorker interface required by VSManager. Each instance is constructed for one (rs, Theta) point of the 3×3 variational-swarm grid and delegates all structural property computations to the underlying STLS solver.

Public Functions

inline VSStlsWorker(const std::shared_ptr<const VSStlsInput> &in)

Construct a worker for the given input parameters.

Parameters

in – Shared pointer to the VS-STLS input (verbosity disabled).

inline virtual const Vector2D &getLfc() const override

Return the local field correction (wave-vectors × Matsubara frequencies).

inline virtual const std::vector<double> &getWvg() const override

Return the wave-vector grid.

inline virtual const std::vector<double> &getSsf() const override

Return the static structure factor over the wave-vector grid.

inline virtual const Vector2D &getIdr() const override

Return the ideal density response.

inline virtual std::vector<double> getSdr() const override

Compute and return the static density response.

Returns

Vector of static density response values.

inline virtual double getUInt() const override

Return the interaction energy per particle.

Returns

Dimensionless interaction energy.

inline virtual double getChemicalPotential() const override

Return the chemical potential.

Returns

Chemical potential (in units of the thermal energy).

inline virtual double getQAdder() const override

Compute the Q-adder contribution (interaction energy at rs = 1).

Returns

Internal energy per particle evaluated at unit coupling strength.

inline virtual double getCoupling() const override

Return the coupling parameter for this grid point.

Returns

rs value.

inline virtual double getDegeneracy() const override

Return the degeneracy parameter for this grid point.

Returns

Theta value.

inline virtual void init() override

Initialize the worker (wave-vector grid, chemical potential, etc.).

inline virtual void initialGuess() override

Set the initial guess for the iterative procedure.

inline virtual void computeSsf() override

Compute the static structure factor using the current SLFC.

inline virtual void computeLfc() override

Compute the static local field correction from the current SSF.

inline virtual void applyLfcDiff(const Vector2D &v) override

Apply a LFC correction by subtracting v from the stored LFC.

Parameters

v – LFC correction to subtract.

inline virtual double computeError() const override

Compute the convergence error between successive iterates.

Returns

RMS difference between ssf and ssfOld.

inline virtual void updateSolution() override

Mix the new and old SSF to advance the iteration.

class VSStlsManager : public VSManager, public Stls
#include <vsstls.hpp>

VS manager for STLS-based schemes.

Combines VSManager (which drives the 3×3 grid iteration) with Stls (whose compute() runs the self-consistency loop for the central grid point). The nine VSStlsWorker objects are initialized in the constructor and coordinated through the VSManager interface.

Public Functions

explicit VSStlsManager(const std::shared_ptr<const VSStlsInput> &in)

Construct the manager and initialize all nine workers.

Parameters

in – Shared pointer to the VS-STLS input parameters.

inline virtual void init() override

Initialize wave-vector grid and other basic properties.

inline virtual void initialGuess() override

Set the initial guess for the iterative procedure.

inline virtual void computeSsf() override

Compute the static structure factor using the current SLFC.

inline virtual void computeLfc() override

Compute the static local field correction from the current SSF.

inline virtual double computeError() const override

Compute the convergence error between successive iterates.

Returns

RMS difference between ssf and ssfOld.

inline virtual void updateSolution() override

Mix the new and old SSF to advance the iteration.

inline virtual int compute() override

Run the STLS self-consistency loop as the manager’s compute method.

Returns

0 on success.

Private Functions

inline virtual const VSInput &inVS() const override

Return the VS-specific input parameters.

Returns

Constant reference to the VSInput.

inline virtual const Input &inScheme() const override

Return the scheme-specific base input parameters.

Returns

Constant reference to the Input.

Private Members

std::shared_ptr<const VSStlsInput> managerInPtr

Shared pointer to the input, used for inVS() and inScheme().

class VSStls : public VSBase
#include <vsstls.hpp>

Variational-swarm solver based on STLS structural properties.

Determines the optimal free parameter alpha by minimizing the exchange-correlation free energy over a 3×3 (rs, Theta) grid, using VSStlsManager to run the STLS solver at each grid point.

Public Functions

explicit VSStls(const std::shared_ptr<const VSStlsInput> &in)

Construct the VS-STLS solver.

Parameters

in – Shared pointer to the VS-STLS input parameters.

int compute()

Expose VSBase::compute() as the public entry point.

Private Functions

inline virtual VSManager &grid() override

Return the mutable VSManager (non-const).

Returns

Reference to the managed grid.

inline virtual const VSManager &grid() const override

Return the VSManager (const).

Returns

Constant reference to the managed grid.

virtual const VSInput &in() const override

Return the VS-specific input parameters.

Returns

Constant reference to VSInput.

virtual const Input &inScheme() const override

Return the scheme-specific base input parameters.

Returns

Constant reference to Input.

Private Members

std::shared_ptr<const VSStlsInput> inPtr

Shared pointer to the input parameters.

VSStlsManager grid_

3×3 grid manager.

class VSQstlsWorker : public VSWorker, public Qstls
#include <qvsstls.hpp>

qSTLS-based worker for a single grid point in the QVS scheme.

Wraps a Qstls solver to provide the VSWorker interface required by VSQstlsManager. Constructed for one (rs, Theta) grid point; all structural property queries delegate to the underlying qSTLS solver.

Public Types

using InputType = QVSStlsInput

Type alias for the input class.

Public Functions

VSQstlsWorker(const std::shared_ptr<const QVSStlsInput> &in, GridPoint p)

Construct a worker for the given input and grid point.

Parameters
  • in – Shared pointer to the QVS-STLS input parameters.

  • p – Grid point that determines the (rs, Theta) offset for this worker.

inline virtual const Vector2D &getLfc() const override

Return the local field correction (wave-vectors × Matsubara frequencies).

inline virtual const std::vector<double> &getWvg() const override

Return the wave-vector grid.

inline virtual const std::vector<double> &getSsf() const override

Return the static structure factor over the wave-vector grid.

inline virtual const Vector2D &getIdr() const override

Return the ideal density response.

inline virtual std::vector<double> getSdr() const override

Compute and return the static density response.

Returns

Vector of static density response values.

inline virtual double getUInt() const override

Return the interaction energy per particle.

Returns

Dimensionless interaction energy.

inline virtual double getChemicalPotential() const override

Return the chemical potential.

Returns

Chemical potential (in units of the thermal energy).

virtual double getQAdder() const override

Compute the Q-adder contribution for the quantum VS free-energy expression.

Returns

Q-adder value at this grid point.

inline virtual double getCoupling() const override

Return the coupling parameter for this grid point.

Returns

rs value.

inline virtual double getDegeneracy() const override

Return the degeneracy parameter for this grid point.

Returns

Theta value.

inline virtual void init() override

Initialize wave-vector grid, chemical potential, and ADR fixed component.

inline virtual void initialGuess() override

Set the initial guess for the iterative solver.

inline virtual void computeSsf() override

Compute the static structure factor using the current LFC.

inline virtual void computeLfc() override

Compute the quantum local field correction from the ADR.

inline virtual void applyLfcDiff(const Vector2D &v) override

Apply a LFC correction by subtracting v from the stored LFC.

Parameters

v – LFC correction to subtract.

inline virtual double computeError() const override

Compute the convergence error between successive SSF iterates.

Returns

RMS difference between the current and previous SSF.

inline virtual void updateSolution() override

Mix the new and old SSF to advance the iteration.

double computeQAdder(const std::shared_ptr<Integrator2D> &itg2D, const std::vector<double> &itgGrid) const

Compute the Q-adder using explicitly provided integrators.

Parameters
  • itg2D – Shared pointer to a 2D numerical integrator.

  • itgGrid – Grid for 2D integration.

Returns

Q-adder value.

class VSQstlsManager : public VSManager, public Qstls
#include <qvsstls.hpp>

VS manager for qSTLS-based schemes.

Combines VSManager (grid iteration) with Qstls (self-consistency loop) to implement the QVS-STLS variational-swarm scheme. Owns nine VSQstlsWorker objects and a shared 2D integrator for the Q-adder computations.

Public Functions

explicit VSQstlsManager(const std::shared_ptr<const QVSStlsInput> &in)

Construct the manager and initialize all nine workers.

Parameters

in – Shared pointer to the QVS-STLS input parameters.

inline virtual void init() override

Initialize wave-vector grid, chemical potential, and ADR fixed component.

inline virtual void initialGuess() override

Set the initial guess for the iterative procedure.

inline virtual void computeSsf() override

Dispatch the SSF calculation to the appropriate temperature branch.

inline virtual void computeLfc() override

Compute the quantum local field correction from the ADR.

inline virtual double computeError() const override

Compute the convergence error between successive iterates.

Returns

RMS difference between ssf and ssfOld.

inline virtual void updateSolution() override

Mix the new and old SSF to advance the iteration.

inline virtual int compute() override

Run the qSTLS self-consistency loop as the manager’s compute method.

Returns

0 on success.

Private Functions

inline virtual const VSInput &inVS() const override

Return the VS-specific input parameters.

Returns

Constant reference to the VSInput.

inline virtual const Input &inScheme() const override

Return the scheme-specific base input parameters.

Returns

Constant reference to the Input.

Private Members

std::shared_ptr<const QVSStlsInput> managerInPtr_

Shared pointer to the input, used for inVS() and inScheme().

std::shared_ptr<Integrator2D> itg2D

Shared 2D integrator for Q-adder computations.

std::vector<double> itgGrid

Grid for 2D integration of the Q-adder.

class QVSStls : public VSBase
#include <qvsstls.hpp>

Variational-swarm solver based on qSTLS structural properties.

Determines the optimal free parameter alpha by minimizing the exchange-correlation free energy over a 3×3 (rs, Theta) grid, using VSQstlsManager to run the qSTLS solver at each grid point.

Public Functions

explicit QVSStls(const std::shared_ptr<const QVSStlsInput> &in)

Construct the QVS-STLS solver.

Parameters

in – Shared pointer to the QVS-STLS input parameters.

int compute()

Expose VSBase::compute() as the public entry point.

Private Functions

inline virtual VSManager &grid() override

Return the mutable VSManager (non-const).

Returns

Reference to the managed grid.

inline virtual const VSManager &grid() const override

Return the VSManager (const).

Returns

Constant reference to the managed grid.

virtual const VSInput &in() const override

Return the VS-specific input parameters.

Returns

Constant reference to VSInput.

virtual const Input &inScheme() const override

Return the scheme-specific base input parameters.

Returns

Constant reference to Input.

Private Members

std::shared_ptr<const QVSStlsInput> inPtr

Shared pointer to the input parameters.

VSQstlsManager grid_

3×3 grid manager.

class QAdder
#include <qvsstls.hpp>

Computes the Q-adder contribution for the quantum VS free-energy expression.

The Q-adder is a functional of the static structure factor that enters the self-consistency equation for the free parameter alpha in the quantum VS (QVS) scheme. It combines 1D and 2D integrations of the SSF with the Fermi–Dirac kernel.

Public Functions

inline QAdder(const double &Theta_, const double &mu_, const double &limitMin, const double &limitMax, const std::vector<double> &itgGrid_, std::shared_ptr<Integrator1D> itg1_, std::shared_ptr<Integrator2D> itg2_, std::shared_ptr<Interpolator1D> interp_)

Construct the Q-adder calculator.

Parameters
  • Theta_ – Degeneracy parameter.

  • mu_ – Chemical potential.

  • limitMin – Lower integration limit.

  • limitMax – Upper integration limit.

  • itgGrid_ – Grid for 2D integration.

  • itg1_ – Shared pointer to a 1D integrator.

  • itg2_ – Shared pointer to a 2D integrator.

  • interp_ – Shared pointer to an SSF interpolator.

double get() const

Compute and return the Q-adder value.

Returns

Q-adder at the current state point.

Private Functions

double ssf(const double &y) const

Evaluate the interpolated SSF at wave-vector y.

Parameters

y – Wave-vector value.

double integrandDenominator(const double q) const

Denominator integrand over q.

Parameters

q – Auxiliary momentum.

double integrandNumerator1(const double q) const

First numerator integrand over q.

Parameters

q – Auxiliary momentum.

double integrandNumerator2(const double w) const

Second numerator integrand over w.

Parameters

w – Auxiliary frequency variable.

void getIntDenominator(double &res) const

Compute the denominator integral.

Parameters

res – Output variable for the result.

Private Members

const double Theta

Degeneracy parameter.

const double mu

Chemical potential.

const std::pair<double, double> limits

Integration limits [min, max].

const std::vector<double> &itgGrid

Grid for 2D integration.

const std::shared_ptr<Integrator1D> itg1

1D numerical integrator.

const std::shared_ptr<Integrator2D> itg2

2D numerical integrator.

const std::shared_ptr<Interpolator1D> interp

Interpolator for the static structure factor.

Input Classes

Variables

constexpr double DEFAULT_DOUBLE = numUtil::NaN

Default value for uninitialized double parameters (signaling NaN).

constexpr int DEFAULT_INT = numUtil::iNaN

Default value for uninitialized integer parameters (signaling NaN).

constexpr bool DEFAULT_BOOL = false

Default value for uninitialized boolean parameters.

class Input
#include <input.hpp>

Base input class for all dielectric scheme solvers.

Stores the thermodynamic state point (coupling parameter rs and degeneracy parameter Theta), numerical grid parameters, and metadata such as the theory name and database information. All derived input classes extend this base with scheme-specific settings.

Subclassed by IterationInput

Public Functions

inline explicit Input()

Construct with all parameters set to their sentinel defaults.

virtual ~Input() = default

Virtual destructor.

void setCoupling(const double &rs)

Set the coupling parameter (Wigner–Seitz radius).

Parameters

rs – Non-negative coupling parameter.

void setDatabaseInfo(const databaseUtil::DatabaseInfo &dbInfo)

Set the database metadata for result storage.

Parameters

dbInfo – Database information struct.

void setDimension(const dimensionsUtil::Dimension &dimension)

Set the spatial dimension of the system.

Parameters

dimension – Dimension enum value (D2 or D3).

void setDegeneracy(const double &Theta)

Set the degeneracy parameter (reduced temperature \(\Theta = k_B T / E_F\)).

Parameters

Theta – Non-negative degeneracy parameter.

void setInt2DScheme(const std::string &int2DScheme)

Set the 2D integration quadrature scheme.

Parameters

int2DScheme – Name of the scheme (e.g., “segregated”).

void setIntError(const double &intError)

Set the relative error tolerance for numerical integrals.

Parameters

intError – Positive relative error tolerance.

void setNThreads(const int &nThreads)

Set the number of OpenMP threads for parallel sections.

Parameters

nThreads – Positive thread count.

void setTheory(const std::string &theory)

Set the name of the dielectric theory to solve.

Parameters

theory – Theory identifier string (e.g., “STLS”, “qSTLS”).

void setChemicalPotentialGuess(const std::vector<double> &muGuess)

Set the initial bracket for the chemical potential root-finder.

Parameters

muGuess – Two-element vector [lower, upper] bracketing the solution.

void setNMatsubara(const int &nMatsubara)

Set the number of Matsubara frequencies used in finite-temperature sums.

Parameters

nMatsubara – Positive integer number of frequencies.

void setWaveVectorGridRes(const double &dx)

Set the wave-vector grid resolution.

Parameters

dx – Positive grid spacing.

void setWaveVectorGridCutoff(const double &xmax)

Set the wave-vector grid cutoff.

Parameters

xmax – Positive maximum wave-vector value.

void setFrequencyCutoff(const double &OmegaMax)

Set the frequency cutoff (ground-state calculations only).

Parameters

OmegaMax – Positive maximum real frequency.

inline double getCoupling() const

Return the coupling parameter.

inline databaseUtil::DatabaseInfo getDatabaseInfo() const

Return the database metadata.

inline dimensionsUtil::Dimension getDimension() const

Return the spatial dimension.

inline double getDegeneracy() const

Return the degeneracy parameter.

inline std::string getInt2DScheme() const

Return the 2D integration quadrature scheme name.

inline double getIntError() const

Return the relative integral error tolerance.

inline int getNThreads() const

Return the number of OpenMP threads.

inline std::string getTheory() const

Return the theory name.

inline bool isClassic() const

Return true if the theory is a classical (non-quantum) approximation.

inline std::vector<double> getChemicalPotentialGuess() const

Return the chemical potential initial guess bracket.

inline int getNMatsubara() const

Return the number of Matsubara frequencies.

inline double getWaveVectorGridRes() const

Return the wave-vector grid resolution.

inline double getWaveVectorGridCutoff() const

Return the wave-vector grid cutoff.

inline double getFrequencyCutoff() const

Return the frequency cutoff.

Protected Attributes

double intError

Relative error tolerance for numerical integrals.

double rs

Quantum coupling parameter (Wigner–Seitz radius).

double Theta

Degeneracy parameter \(\Theta = k_B T / E_F\).

int nThreads

Number of OpenMP threads for parallel sections.

bool isClassicTheory

True if the theory is a classical (non-quantum) approximation.

bool isQuantumTheory

True if the theory includes quantum (exchange-correlation) effects.

std::string int2DScheme

2D quadrature scheme name.

std::string theory

Theory identifier string.

databaseUtil::DatabaseInfo dbInfo

Database metadata for result storage.

dimensionsUtil::Dimension dimension

Spatial dimension (default: 3D).

double dx

Wave-vector grid resolution.

double xmax

Wave-vector grid cutoff.

double OmegaMax

Frequency cutoff (relevant only for ground-state calculations).

int nl

Number of Matsubara frequencies.

std::vector<double> muGuess

Initial bracket for the chemical potential root-finder.

struct Guess
#include <input.hpp>

Carries the initial guess for iterative schemes.

Holds a pre-computed wave-vector grid, static structure factor, and local field correction that can be used to warm-start the iterative solver.

Public Members

std::vector<double> wvg

Wave-vector grid.

std::vector<double> ssf

Static structure factor values.

Vector2D lfc

Local field correction values (wave-vectors × Matsubara frequencies).

class IterationInput : public Input
#include <input.hpp>

Input class for iteratively solved dielectric schemes.

Extends Input with convergence criteria (minimum error, maximum iterations) and a linear-mixing parameter for the iterative updates. An optional initial guess can be provided to warm-start the iteration.

Subclassed by StlsInput

Public Functions

inline explicit IterationInput()

Construct with all iteration parameters set to their sentinel defaults.

void setErrMin(const double &errMin)

Set the convergence threshold.

Parameters

errMin – Positive minimum RMS error for stopping the iteration.

void setGuess(const Guess &guess)

Set the initial guess for the iterative solver.

Parameters

guess – Struct containing wvg, ssf, and lfc arrays.

void setMixingParameter(const double &aMix)

Set the linear-mixing parameter.

Parameters

aMix – Mixing fraction in (0, 1].

void setNIter(const int &nIter)

Set the maximum number of iterations.

Parameters

nIter – Positive integer maximum.

inline double getErrMin() const

Return the convergence threshold.

inline Guess getGuess() const

Return the initial guess struct.

inline double getMixingParameter() const

Return the linear-mixing parameter.

inline int getNIter() const

Return the maximum number of iterations.

Protected Attributes

double aMix

Linear-mixing parameter for SSF updates.

double errMin

Minimum RMS error for convergence.

int nIter

Maximum number of self-consistency iterations.

Guess guess

Optional warm-start initial guess.

class QuantumInput
#include <input.hpp>

Mixin input class for quantum (qSTLS and qSTLS-IET) schemes.

Carries the run identifier used to load a pre-computed fixed component of the auxiliary density response from the database.

Subclassed by QstlsInput

Public Functions

void setFixedRunId(const int &fixedRunId)

Set the run ID for the fixed ADR component.

Parameters

fixedRunId – Integer run identifier in the database.

inline int getFixedRunId() const

Return the fixed ADR run identifier.

Protected Attributes

int fixedRunId

Database run ID for the fixed ADR component.

class IetInput
#include <input.hpp>

Mixin input class for IET-based schemes.

Carries the name of the quantum-to-classical mapping used when evaluating the bridge function (e.g., “HNC”, “IOI”, “LCT”).

Subclassed by QstlsIetInput, StlsIetInput

Public Functions

virtual ~IetInput() = default

Virtual destructor.

void setMapping(const std::string &mapping)

Set the quantum-classical mapping for the bridge function.

Parameters

mapping – Mapping identifier string.

inline std::string getMapping() const

Return the mapping identifier.

Protected Attributes

std::string mapping

Quantum-classical mapping identifier for the IET bridge function.

class StlsInput : public IterationInput
#include <input.hpp>

Input class for the STLS scheme.

Inherits all iteration parameters from IterationInput with no additional STLS-specific settings.

Subclassed by QstlsInput, StlsIetInput, VSStlsInput

Public Functions

explicit StlsInput() = default

Construct with default parameters.

class StlsIetInput : public StlsInput, public IetInput
#include <input.hpp>

Input class for the STLS-IET scheme.

Combines STLS iteration parameters with the IET mapping setting.

Public Functions

explicit StlsIetInput() = default

Construct with default parameters.

class QstlsInput : public StlsInput, public QuantumInput
#include <input.hpp>

Input class for the qSTLS scheme.

Combines STLS iteration parameters with the fixed ADR run identifier.

Subclassed by QVSStlsInput, QstlsIetInput

Public Functions

explicit QstlsInput() = default

Construct with default parameters.

class QstlsIetInput : public QstlsInput, public IetInput
#include <input.hpp>

Input class for the qSTLS-IET scheme.

Combines qSTLS parameters with the IET mapping setting.

Public Functions

explicit QstlsIetInput() = default

Construct with default parameters.

class VSInput
#include <input.hpp>

Mixin input class for variational-swarm (VS) schemes.

Carries the free-parameter guess, finite-difference grid resolutions, convergence settings for the alpha iterations, and optionally a pre-computed free-energy integrand to warm-start the VS optimization.

Subclassed by QVSStlsInput, VSStlsInput

Public Functions

inline explicit VSInput()

Construct with all VS parameters set to their sentinel defaults.

virtual ~VSInput() = default

Virtual destructor.

void setAlphaGuess(const std::vector<double> &alphaGuess)

Set the initial bracket for the free parameter alpha.

Parameters

alphaGuess – Two-element vector bracketing alpha.

void setCouplingResolution(const double &drs)

Set the coupling-parameter grid resolution for finite differences.

Parameters

drs – Positive grid spacing in rs.

void setDegeneracyResolution(const double &dTheta)

Set the degeneracy-parameter grid resolution for finite differences.

Parameters

dTheta – Positive grid spacing in Theta.

void setErrMinAlpha(const double &errMinAlpha)

Set the convergence threshold for the alpha iterations.

Parameters

errMinAlpha – Positive minimum error for stopping.

void setNIterAlpha(const int &nIterAlpha)

Set the maximum number of alpha iterations.

Parameters

nIterAlpha – Positive integer maximum.

void setFreeEnergyIntegrand(const FreeEnergyIntegrand &freeEnergyIntegrand)

Provide a pre-computed free-energy integrand for warm-starting.

Parameters

freeEnergyIntegrand – Struct containing grid and integrand values.

inline std::vector<double> getAlphaGuess() const

Return the alpha guess bracket.

inline double getCouplingResolution() const

Return the coupling-parameter grid resolution.

inline double getDegeneracyResolution() const

Return the degeneracy-parameter grid resolution.

inline double getErrMinAlpha() const

Return the alpha convergence threshold.

inline double getNIterAlpha() const

Return the maximum number of alpha iterations.

inline FreeEnergyIntegrand getFreeEnergyIntegrand() const

Return the pre-computed free-energy integrand.

Private Members

std::vector<double> alphaGuess

Initial bracket for the free parameter alpha.

double drs

Coupling-parameter grid resolution (finite-difference step in rs).

double dTheta

Degeneracy-parameter grid resolution (finite-difference step in Theta).

double errMinAlpha

Convergence threshold for the alpha iterations.

int nIterAlpha

Maximum number of alpha iterations.

FreeEnergyIntegrand fxcIntegrand

Pre-computed free-energy integrand for warm-starting.

struct FreeEnergyIntegrand
#include <input.hpp>

Pre-computed free-energy integrand for warm-starting.

The grid field is the coupling-parameter grid and integrand contains the per-theta integrand values.

Public Members

std::vector<double> grid

Coupling-parameter grid.

std::vector<std::vector<double>> integrand

Free-energy integrand values (indexed by theta, then rs).

class VSStlsInput : public VSInput, public StlsInput
#include <input.hpp>

Input class for the VS-STLS scheme.

Combines VS optimization parameters with STLS iteration parameters.

Public Functions

explicit VSStlsInput() = default

Construct with default parameters.

class QVSStlsInput : public VSInput, public QstlsInput
#include <input.hpp>

Input class for the QVS-STLS (quantum VS-STLS) scheme.

Combines VS optimization parameters with qSTLS iteration parameters (including the fixed ADR run identifier).

Public Functions

explicit QVSStlsInput() = default

Construct with default parameters.

VS

class VSBase : public Logger
#include <vsbase.hpp>

Base class orchestrating the variational-swarm (VS) free-energy optimization.

VSBase drives the outer loop that determines the free variational parameter alpha by minimizing the exchange-correlation free energy. For a given alpha it delegates the self-consistency iterations to a VSManager, which runs nine (rs, Theta) grid points in parallel. The free-energy integrand is accumulated from the interaction energy at each grid point and integrated over the coupling-parameter grid.

Derived classes implement the abstract interface methods (in(), inScheme(), grid()) to supply the scheme-specific input and manager objects.

Subclassed by QVSStls, VSStls

Public Functions

inline explicit VSBase()

Default constructor.

virtual ~VSBase() = default

Virtual destructor.

int compute()

Run the full VS optimization pipeline.

Returns

0 on success.

inline double getAlpha() const

Return the converged value of the free parameter alpha.

Returns

Optimized alpha value.

const std::vector<std::vector<double>> &getFreeEnergyIntegrand() const

Return the free-energy integrand accumulated during the optimization.

Returns

2D array indexed by [theta-row][rs-column].

const std::vector<double> &getFreeEnergyGrid() const

Return the coupling-parameter grid used for the free-energy integration.

Returns

Coupling-parameter grid vector.

const std::vector<double> &getSsf() const

Return the SSF at the central (target) grid point.

Returns

Constant reference to the SSF vector.

const Vector2D &getLfc() const

Return the LFC at the central (target) grid point.

Returns

Constant reference to the LFC array.

const std::vector<double> &getWvg() const

Return the wave-vector grid at the central (target) grid point.

Returns

Constant reference to the wave-vector grid.

const Vector2D &getIdr() const

Return the IDR at the central (target) grid point.

Returns

Constant reference to the IDR array.

std::vector<double> getSdr() const

Compute and return the static density response at the central grid point.

Returns

Static density response vector.

double getUInt() const

Return the interaction energy at the central grid point.

Returns

Dimensionless interaction energy per particle.

double getChemicalPotential() const

Return the chemical potential at the central grid point.

Returns

Chemical potential (in units of the thermal energy).

double getError() const

Return the convergence error from the last iteration.

Returns

RMS difference between successive SSF iterates.

Protected Functions

virtual const VSInput &in() const = 0

Return the VS-specific input parameters.

Returns

Constant reference to VSInput.

virtual const Input &inScheme() const = 0

Return the scheme-specific base input parameters.

Returns

Constant reference to Input.

virtual VSManager &grid() = 0

Return the mutable VSManager (non-const).

Returns

Reference to the managed grid.

virtual const VSManager &grid() const = 0

Return the VSManager (const).

Returns

Constant reference to the managed grid.

int runGrid()

Run one full self-consistency sweep over all grid points.

Returns

0 on success.

double getCoupling(GridPoint p) const

Return the coupling parameter at a given grid point.

Parameters

p – Grid point identifier.

Returns

rs value.

double getDegeneracy(GridPoint p) const

Return the degeneracy parameter at a given grid point.

Parameters

p – Grid point identifier.

Returns

Theta value.

double getFxcIntegrandValue(GridPoint p) const

Return the free-energy integrand value at a given grid point.

Parameters

p – Grid point identifier.

Returns

Free-energy integrand value.

double getQAdder(GridPoint p) const

Return the Q-adder (free-energy derivative) at a given grid point.

Parameters

p – Grid point identifier.

Returns

Q-adder value.

void setRsGrid()

Build the coupling-parameter grid for the free-energy integration.

void setFxcIntegrand()

Initialize the free-energy integrand array to the correct shape.

void updateFxcIntegrand()

Update the free-energy integrand values from the current grid state.

double computeFreeEnergy(GridPoint p, bool normalize) const

Compute the exchange-correlation free energy at a grid point.

Parameters
  • p – Grid point identifier.

  • normalize – If true, divide by \(r_s^2\).

Returns

Free-energy value.

std::vector<double> getFreeEnergyData() const

Collect the free-energy integrand column for the target Theta row.

Returns

Vector of integrand values over the rs grid.

GridPoint getOutputGridPoint() const

Return the grid point corresponding to the target state point.

Returns

GridPoints::CENTER.

Protected Attributes

double alpha

Converged value of the free variational parameter alpha.

std::vector<double> rsGrid

Coupling-parameter grid for the free-energy integral.

std::vector<std::vector<double>> fxcIntegrand

Free-energy integrand values.

Outer index: theta row (DOWN = 0, CENTER = 1, UP = 2). Inner index: rs column (same ordering as rsGrid).

Private Functions

double computeAlpha()

Find the optimal alpha by minimizing the free energy.

Returns

Converged alpha value.

std::vector<double> computeQData()

Compute the Q-data vector used in the alpha minimization.

Returns

Q-data values for all rs grid points.

void doIterations()

Drive the alpha iteration loop to convergence.

double alphaDifference(const double &alphaTmp)

Evaluate the residual of the alpha self-consistency equation.

Parameters

alphaTmp – Trial alpha value.

Returns

Residual value.

class VSManager
#include <vsmanager.hpp>

Manages the 3×3 grid of workers for the variational-swarm scheme.

VSManager owns nine VSWorker instances arranged on a 3×3 grid in the (rs, Theta) plane. It coordinates the self-consistency iterations by driving each worker through the init / guess / SSF / LFC cycle and then applying finite-difference LFC derivatives across the grid to implement the variational-swarm local-field-correction correction.

Concrete subclasses supply input accessors and, via compute(), the per-scheme initialization of the workers array.

Subclassed by VSQstlsManager, VSStlsManager

Public Functions

inline void setAlpha(double alpha_)

Set the free variational parameter alpha.

Parameters

alpha_ – New value of alpha.

inline double getAlpha() const

Return the current value of the free parameter alpha.

Returns

Alpha value.

inline double getError() const

Return the convergence error from the last iteration.

Returns

RMS difference between successive SSF iterates (maximum over all workers).

virtual int compute() = 0

Run one full self-consistency sweep over all grid points.

Implemented by derived classes to initialize workers and drive the iterative loop until convergence.

Returns

0 on success.

const std::vector<double> &getSsf(GridPoint p) const

Return the SSF at a given grid point.

Parameters

p – Grid point identifier.

Returns

Constant reference to the SSF vector.

const Vector2D &getLfc(GridPoint p) const

Return the LFC at a given grid point.

Parameters

p – Grid point identifier.

Returns

Constant reference to the LFC array.

const std::vector<double> &getWvg(GridPoint p) const

Return the wave-vector grid at a given grid point.

Parameters

p – Grid point identifier.

Returns

Constant reference to the wave-vector grid.

const Vector2D &getIdr(GridPoint p) const

Return the ideal density response at a given grid point.

Parameters

p – Grid point identifier.

Returns

Constant reference to the IDR array.

std::vector<double> getSdr(GridPoint p) const

Compute and return the static density response at a given grid point.

Parameters

p – Grid point identifier.

Returns

Static density response vector.

double getCoupling(GridPoint p) const

Return the coupling parameter at a given grid point.

Parameters

p – Grid point identifier.

Returns

rs value.

double getDegeneracy(GridPoint p) const

Return the degeneracy parameter at a given grid point.

Parameters

p – Grid point identifier.

Returns

Theta value.

double getUInt(GridPoint p) const

Return the interaction energy at a given grid point.

Parameters

p – Grid point identifier.

Returns

Dimensionless interaction energy per particle.

double getChemicalPotential(GridPoint p) const

Return the chemical potential at a given grid point.

Parameters

p – Grid point identifier.

Returns

Chemical potential (in units of the thermal energy).

double getQAdder(GridPoint p) const

Return the free-energy integrand contribution (Q-adder) at a grid point.

Parameters

p – Grid point identifier.

Returns

Q-adder value.

double getFxcIntegrandValue(GridPoint p) const

Return the free-energy integrand value at a given grid point.

Parameters

p – Grid point identifier.

Returns

Free-energy integrand value.

Protected Functions

inline explicit VSManager()

Construct with alpha initialized to infinity and initDone = false.

void setupDerivativeData()

Populate rsDerivData and thetaDerivData from the grid layout.

void init()

Initialize all workers and set up derivative metadata.

void computeLfc()

Compute the LFC for all workers, then update LFC derivatives.

void computeSsf()

Compute the SSF for all workers.

double computeError() const

Compute the maximum convergence error across all workers.

Returns

Maximum per-worker RMS difference between successive SSF iterates.

void updateSolution()

Mix the new and old SSF for all workers.

void initialGuess()

Set the initial SSF guess for all workers.

virtual const VSInput &inVS() const = 0

Return the VS-specific input parameters.

Returns

Constant reference to the VSInput.

virtual const Input &inScheme() const = 0

Return the scheme-specific base input parameters.

Returns

Constant reference to the Input.

Protected Attributes

std::array<std::unique_ptr<VSWorker>, N> workers

The nine worker objects, one per grid point.

std::array<Vector2D, N> lfcDerivatives

LFC finite-difference derivatives, one per grid point.

std::array<DerivativeData, N> rsDerivData

Finite-difference metadata along the rs axis.

std::array<DerivativeData, N> thetaDerivData

Finite-difference metadata along the Theta axis.

std::array<double, N> rsValues

rs values for each of the nine grid points.

std::array<double, N> thetaValues

Theta values for each of the nine grid points.

double alpha

Current value of the free variational parameter alpha.

bool initDone

True after the first call to init(); prevents re-initialization.

Protected Static Attributes

static constexpr int N = 9

Number of grid points (3 × 3 = 9).

Private Functions

void computeLfcDerivatives()

Compute and store LFC finite-difference derivatives for all grid points.

double derivative(const Vector2D &f, int l, size_t i, DerivativeData::Type t) const

Compute the finite difference of LFC column l at grid index i.

Parameters
  • f – LFC array.

  • l – Matsubara frequency column index.

  • i – Flat grid index.

  • t – Derivative type (centered / forward / backward).

Returns

Finite-difference value.

double derivative(double f0, double f1, double f2, DerivativeData::Type t) const

Compute a scalar finite difference from three adjacent values.

Parameters
  • f0 – “Down” value.

  • f1 – “Center” value.

  • f2 – “Up” value.

  • t – Derivative type.

Returns

Finite-difference value.

struct DerivativeData
#include <vsmanager.hpp>

Metadata for finite-difference derivative computation at one grid point.

Public Types

enum class Type

Derivative type for the finite-difference stencil.

Values:

enumerator CENTERED
enumerator FORWARD
enumerator BACKWARD

Public Members

Type type

Stencil type to use at this grid point.

size_t upIdx

Flat index of the “up” neighbour in the stencil.

size_t downIdx

Flat index of the “down” neighbour in the stencil.

class VSWorker
#include <vsworker.hpp>

Abstract interface for a single-point variational-swarm worker.

Each VSWorker computes structural and thermodynamic properties for one (rs, Theta) point in the variational-swarm 3×3 grid. VSManager holds an array of nine workers and coordinates the synchronized local-field- correction update that mixes derivatives across grid points.

Concrete workers are created by the scheme-specific VSManager subclass (e.g., one wrapping a Stls solver, another wrapping a Qstls solver).

Subclassed by VSQstlsWorker, VSStlsWorker

Public Functions

virtual ~VSWorker() = default

Virtual destructor.

virtual const Vector2D &getLfc() const = 0

Return the local field correction (wave-vectors × Matsubara frequencies).

virtual const std::vector<double> &getWvg() const = 0

Return the wave-vector grid.

virtual const std::vector<double> &getSsf() const = 0

Return the static structure factor over the wave-vector grid.

virtual const Vector2D &getIdr() const = 0

Return the ideal density response.

virtual std::vector<double> getSdr() const = 0

Compute and return the static density response.

Returns

Vector of static density response values.

virtual double getUInt() const = 0

Return the interaction energy per particle.

Returns

Dimensionless interaction energy.

virtual double getChemicalPotential() const = 0

Return the chemical potential.

Returns

Chemical potential (in units of the thermal energy).

virtual double getQAdder() const = 0

Return the free-energy integrand contribution (Q-adder).

Returns

Q-adder value at this grid point.

virtual double getCoupling() const = 0

Return the coupling parameter for this grid point.

Returns

rs value.

virtual double getDegeneracy() const = 0

Return the degeneracy parameter for this grid point.

Returns

Theta value.

virtual void init() = 0

Initialize the worker (wave-vector grid, chemical potential, etc.).

virtual void initialGuess() = 0

Set the initial guess for the iterative solver.

virtual void computeSsf() = 0

Compute the static structure factor using the current LFC.

virtual void computeLfc() = 0

Compute the local field correction using the current SSF.

virtual void applyLfcDiff(const Vector2D &v) = 0

Apply a corrected LFC difference to update the worker’s LFC.

Used by VSManager to propagate the finite-difference-corrected LFC back to each worker after the synchronized derivative computation.

Parameters

v – LFC correction array to add.

virtual double computeError() const = 0

Compute the convergence error between successive SSF iterates.

Returns

RMS difference between the current and previous SSF.

virtual void updateSolution() = 0

Mix the new and old SSF to advance the iteration.

struct GridPoint
#include <grid_point.hpp>

Strong type for addressing a point in the 3×3 (rs, Theta) state-point grid.

The variational-swarm optimization evaluates the free energy on a 3×3 grid of (rs, Theta) values centred on the target state point. GridPoint encodes a grid position using typed enumerators for the coupling-parameter (Rs) and degeneracy-parameter (Theta) axes, and provides a mapping to the flat 0–8 index used internally by VSManager.

Public Types

enum class Rs

Position along the coupling-parameter (rs) axis.

Values:

enumerator DOWN
enumerator CENTER
enumerator UP
enum class Theta

Position along the degeneracy-parameter (Theta) axis.

Values:

enumerator DOWN
enumerator CENTER
enumerator UP

Public Functions

inline constexpr size_t toIndex() const

Map this grid point to a flat index in [0, 8].

The mapping uses theta as the outer (slow) index and rs as the inner (fast) index, matching the legacy StructIdx convention.

Returns

Flat index in [0, 8].

Public Members

Rs rs

Rs-axis position of this grid point.

Theta theta

Theta-axis position of this grid point.

namespace GridPoints

Named constants for all nine grid points in the 3×3 (rs, Theta) grid.

The naming convention is RS_{rs-position}_THETA_{theta-position}, where positions are DOWN, (CENTER omitted), or UP. The centre point is simply CENTER.

Variables

constexpr GridPoint RS_DOWN_THETA_DOWN = {GridPoint::Rs::DOWN, GridPoint::Theta::DOWN}

Grid point at (rs–, Theta–).

constexpr GridPoint RS_THETA_DOWN = {GridPoint::Rs::CENTER, GridPoint::Theta::DOWN}

Grid point at (rs, Theta–).

constexpr GridPoint RS_UP_THETA_DOWN = {GridPoint::Rs::UP, GridPoint::Theta::DOWN}

Grid point at (rs+, Theta–).

constexpr GridPoint RS_DOWN_THETA = {GridPoint::Rs::DOWN, GridPoint::Theta::CENTER}

Grid point at (rs–, Theta).

constexpr GridPoint CENTER = {GridPoint::Rs::CENTER, GridPoint::Theta::CENTER}

Grid point at the target state point (rs, Theta).

constexpr GridPoint RS_UP_THETA = {GridPoint::Rs::UP, GridPoint::Theta::CENTER}

Grid point at (rs+, Theta).

constexpr GridPoint RS_DOWN_THETA_UP = {GridPoint::Rs::DOWN, GridPoint::Theta::UP}

Grid point at (rs–, Theta+).

constexpr GridPoint RS_THETA_UP = {GridPoint::Rs::CENTER, GridPoint::Theta::UP}

Grid point at (rs, Theta+).

constexpr GridPoint RS_UP_THETA_UP = {GridPoint::Rs::UP, GridPoint::Theta::UP}

Grid point at (rs+, Theta+).

Thermodynamic Properties

namespace thermoUtil

High-level utility functions for computing thermodynamic properties.

These free functions provide a convenient interface for evaluating the most commonly requested thermodynamic quantities directly from the arrays returned by the dielectric scheme solvers.

Functions

double computeInternalEnergy(const std::vector<double> &wvg, const std::vector<double> &ssf, const double &coupling, const dimensionsUtil::Dimension &dim)

Compute the interaction energy per particle.

Integrates \(v(k)\,[S(k) - 1]\) over the wave-vector grid using the dimension-appropriate Coulomb potential.

Parameters
  • wvg – Wave-vector grid.

  • ssf – Static structure factor values (same length as wvg).

  • coupling – Coupling parameter \(r_s\).

  • dim – Spatial dimension.

Returns

Dimensionless interaction energy per particle.

double computeFreeEnergy(const std::vector<double> &grid, const std::vector<double> &rsu, const double &coupling)

Compute the exchange-correlation free energy (without normalization).

Integrates \(r_s \cdot u(r_s) / r_s\) over the coupling-parameter grid up to coupling.

Parameters
  • grid – Coupling-parameter grid.

  • rsu – Pre-computed \(r_s \cdot u(r_s)\) values (same length as grid).

  • coupling – Target coupling parameter (upper integration limit).

Returns

Exchange-correlation free energy.

double computeFreeEnergy(const std::vector<double> &grid, const std::vector<double> &rsu, const double &coupling, const bool normalize)

Compute the exchange-correlation free energy with optional normalization.

Parameters
  • grid – Coupling-parameter grid.

  • rsu – Pre-computed \(r_s \cdot u(r_s)\) values.

  • coupling – Target coupling parameter.

  • normalize – If true, divide the result by \(r_s^2\).

Returns

Exchange-correlation free energy (optionally normalized).

std::vector<double> computeRdf(const std::vector<double> &r, const std::vector<double> &wvg, const std::vector<double> &ssf, const dimensionsUtil::Dimension &dim)

Compute the radial distribution function over a real-space grid.

Fourier-transforms the static structure factor to obtain g(r) at each point in r.

Parameters
  • r – Real-space distance grid.

  • wvg – Wave-vector grid.

  • ssf – Static structure factor values.

  • dim – Spatial dimension.

Returns

Vector of g(r) values, one per entry in r.

Vector2D computeItcf(const std::shared_ptr<const Input> &in, const std::vector<double> &wvg, const std::vector<double> &tauValues, const double mu, const Vector2D &idr, const Vector2D &lfc)

Compute the imaginary-time correlation function.

First computes the non-interacting (HF) contribution. For HF theory (in->getTheory() == “HF”) the result is returned immediately. For all other theories the interacting correction is applied using the local field correction.

Parameters
  • in – Shared pointer to the input parameters.

  • wvg – Wave-vector grid.

  • tauValues – Imaginary time values in [0, 1] (normalised by beta).

  • mu – Chemical potential.

  • idr – Ideal density response (rows = wave-vectors, columns = Matsubara frequencies).

  • lfc – Local field correction (rows = wave-vectors, columns = Matsubara frequencies or 1 for static LFC). Ignored for HF theory.

Returns

2D vector containing ITCF values (rows = wave-vectors, columns = tau values).

class ChemicalPotential : public dimensionsUtil::DimensionsHandler
#include <chemical_potential.hpp>

Computes the chemical potential from the density normalization condition.

The chemical potential is found by solving the equation \(\int_0^\infty f(\epsilon; \mu)\, g(\epsilon)\, d\epsilon = n\) where \(f\) is the Fermi–Dirac distribution and \(g\) is the density of states. Both 2D and 3D geometries are supported via the DimensionsHandler dispatch mechanism.

Public Functions

inline explicit ChemicalPotential(const std::shared_ptr<const Input> in_)

Construct with the simulation input parameters.

Parameters

in_ – Shared pointer to the input parameters.

inline double get() const

Return the computed chemical potential.

Returns

Chemical potential in units of the Fermi energy.

Private Functions

virtual void compute2D() override

Compute the chemical potential for a 2D system.

virtual void compute3D() override

Compute the chemical potential for a 3D system.

double normalizationCondition(const double &mu) const

Evaluate the normalization condition at a trial chemical potential.

Used as the root-finding objective function.

Parameters

mu – Trial chemical potential value.

Returns

Residual: computed density minus target density.

Private Members

const std::shared_ptr<const Input> in

Shared pointer to the input parameters.

double mu = DEFAULT_DOUBLE

Computed chemical potential (initialized to the sentinel NaN).

class InternalEnergy
#include <internal_energy.hpp>

Computes the interaction (internal) energy per particle.

Evaluates the coupling-constant integral \(u(r_s) = \frac{1}{2} \int_0^\infty v(k)\, [S(k) - 1]\, dk / (2\pi)^d\) where \(v(k)\) is the Coulomb potential in \(d\) dimensions and \(S(k)\) is the static structure factor. Both 2D and 3D geometries are supported.

Public Functions

inline InternalEnergy(const double &rs_, const double &yMin_, const double &yMax_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Integrator1D> itg_, const dimensionsUtil::Dimension &dim_)

Construct with the integration data and parameters.

Parameters
  • rs_ – Coupling parameter (Wigner–Seitz radius).

  • yMin_ – Lower integration limit (minimum wave-vector).

  • yMax_ – Upper integration limit (maximum wave-vector).

  • ssfi_ – Shared pointer to an SSF interpolator.

  • itg_ – Shared pointer to a 1D numerical integrator.

  • dim_ – Spatial dimension of the system.

double get() const

Compute and return the interaction energy per particle.

Returns

Dimensionless interaction energy \(u / (e^2 / a_0)\).

Private Functions

double integrand(const double &y) const

Evaluate the integrand at wave-vector y.

Parameters

y – Wave-vector value.

Returns

Integrand value \(v(y)\,[S(y) - 1]\) at y.

double ssf(const double &y) const

Evaluate the interpolated SSF at wave-vector y.

Parameters

y – Wave-vector value.

Returns

Interpolated SSF value.

Private Members

const double rs

Coupling parameter.

const double yMin

Lower integration limit (minimum wave-vector).

const double yMax

Upper integration limit (maximum wave-vector).

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

const std::shared_ptr<Interpolator1D> ssfi

Interpolator for the static structure factor.

const dimensionsUtil::Dimension dim

Spatial dimension.

class FreeEnergy
#include <free_energy.hpp>

Computes the exchange-correlation free energy by integrating over the coupling parameter.

The free energy is obtained from the coupling-constant integration \(f_{xc}(r_s) = \int_0^{r_s} \frac{u(r_s')}{r_s'}\, dr_s'\) where \(u(r_s')\) is the interaction energy per particle at coupling parameter \(r_s'\). The integrand \(r_s' u(r_s')\) is provided via a 1D spline interpolator.

Public Functions

inline FreeEnergy(const double &rs_, std::shared_ptr<Interpolator1D> rsui_, std::shared_ptr<Integrator1D> itg_, const bool normalize_)

Construct with the integration data and options.

Parameters
  • rs_ – Target coupling parameter (upper integration limit).

  • rsui_ – Interpolator for the integrand \(r_s \cdot u(r_s)\).

  • itg_ – Shared pointer to a 1D numerical integrator.

  • normalize_ – If true, divide the result by \(r_s^2\).

double get() const

Compute and return the free energy.

Returns

Exchange-correlation free energy (optionally normalized by \(r_s^2\)).

Private Functions

double integrand(const double y) const

Evaluate the coupling-constant integrand at y.

Parameters

y – Coupling parameter value.

Returns

Integrand value \((r_s' \cdot u(r_s')) / r_s'\) at y.

Private Members

const double rs

Target coupling parameter (upper integration limit).

const std::shared_ptr<Integrator1D> itg

1D numerical integrator.

const std::shared_ptr<Interpolator1D> rsui

Interpolator for the integrand \(r_s \cdot u(r_s)\).

const bool normalize

If true, the result is divided by \(r_s^2\).

class Rdf : public dimensionsUtil::DimensionsHandler
#include <rdf.hpp>

Computes the radial distribution function (RDF) at a single real-space distance.

The RDF is obtained by a Fourier transform of the static structure factor: \(g(r) = 1 + \frac{1}{n} \int \frac{d^d k}{(2\pi)^d}\, e^{i\mathbf{k}\cdot\mathbf{r}} [S(k) - 1]\). Both 2D and 3D geometries are supported via the DimensionsHandler dispatch.

Public Functions

inline Rdf(const double &r_, const double &cutoff_, std::shared_ptr<Interpolator1D> ssfi_, std::shared_ptr<Integrator1D> itg_, std::shared_ptr<Integrator1D> itgf_, const dimensionsUtil::Dimension &dim_)

Construct for evaluation at a single real-space point.

Parameters
  • r_ – Real-space distance at which to evaluate g(r).

  • cutoff_ – Wave-vector cutoff (upper limit of the SSF grid).

  • ssfi_ – Shared pointer to an SSF interpolator.

  • itg_ – Shared pointer to a 1D numerical integrator (standard integrals).

  • itgf_ – Shared pointer to a 1D Fourier integrator.

  • dim_ – Spatial dimension of the system.

double get()

Compute and return g(r) at the configured distance.

Returns

Radial distribution function value at r.

Private Functions

double integrand(const double &y) const

3D integrand \(k^2 [S(k) - 1] \sin(kr) / (kr)\).

Parameters

y – Wave-vector value.

double integrand2D(const double &y) const

2D integrand \(k [S(k) - 1] J_0(kr)\).

Parameters

y – Wave-vector value.

double ssf(const double &y) const

Evaluate the interpolated SSF at wave-vector y.

Parameters

y – Wave-vector value.

Returns

Interpolated SSF value.

virtual void compute2D() override

Compute the RDF for a 2D system.

virtual void compute3D() override

Compute the RDF for a 3D system.

Private Members

const double r

Real-space distance.

const double cutoff

Wave-vector cutoff (upper integration limit).

const std::shared_ptr<Integrator1D> itgf

Fourier-type 1D integrator.

const std::shared_ptr<Integrator1D> itg

Standard 1D integrator.

const std::shared_ptr<Interpolator1D> ssfi

Interpolator for the static structure factor.

const dimensionsUtil::Dimension dim

Spatial dimension of the system.

double res

Result of the RDF computation.

Utility Classes

class Logger
#include <logger.hpp>

Provides controlled stdout output with optional verbosity.

Base class for components that need to print progress or diagnostic messages. Output can be suppressed by constructing with verbose_ = false.

Subclassed by HF, Iet, VSBase

Protected Functions

inline Logger(const bool &verbose_)

Construct with explicit verbosity flag.

Parameters

verbose_ – If false, all print calls are silenced.

inline Logger()

Construct with verbosity enabled.

void print(const std::string &msg) const

Print a message without a trailing newline.

Parameters

msg – Message to print.

void println(const std::string &msg) const

Print a message followed by a newline.

Parameters

msg – Message to print.

Private Functions

bool log_message() const

Returns true if the message should be logged.

Private Members

const bool verbose

Whether messages should be printed.

class Interpolator1D
#include <numerics.hpp>

1D cubic-spline interpolator backed by GSL.

Wraps gsl_spline with caching of the last evaluation point via a gsl_interp_accel. Values outside the data range are extrapolated by holding the boundary value (constant extrapolation beyond cutoff).

Public Functions

Interpolator1D(const std::vector<double> &x, const std::vector<double> &y)

Construct from separate x and y data vectors.

Parameters
  • x – Strictly increasing abscissa values.

  • y – Ordinate values (same length as x).

Interpolator1D(const double &x, const double &y, const size_t n_)

Construct from raw pointers.

Parameters
  • x – Pointer to the first abscissa value.

  • y – Pointer to the first ordinate value.

  • n_ – Number of data points.

explicit Interpolator1D()

Construct an empty (invalid) interpolator.

~Interpolator1D()

Destructor — releases GSL resources.

bool isValid() const

Return true if the interpolator holds valid data.

Returns

False for default-constructed or reset-but-not-populated objects.

void reset(const double &x, const double &y, const size_t n_)

Reinitialize with new raw-pointer data.

Parameters
  • x – Pointer to the first abscissa value.

  • y – Pointer to the first ordinate value.

  • n_ – Number of data points.

double eval(const double &x) const

Evaluate the spline at x.

Parameters

x – Evaluation point.

Returns

Interpolated value (constant extrapolation beyond the data range).

Private Functions

void setup(const double &x, const double &y, const size_t n_)

Shared initialization helper.

Parameters
  • x – Pointer to the first abscissa value.

  • y – Pointer to the first ordinate value.

  • n_ – Number of data points.

Private Members

const gsl_interp_type *TYPE = gsl_interp_cspline

GSL spline type (cubic).

gsl_spline *spline

GSL spline object.

gsl_interp_accel *acc

GSL interpolation accelerator.

double cutoff

Maximum abscissa value of the data range.

size_t n

Number of data points.

class Interpolator2D
#include <numerics.hpp>

2D bicubic-spline interpolator backed by GSL.

Wraps gsl_spline2d. The data must be on a regular rectangular grid with nx abscissa values and ny ordinate values.

Public Functions

Interpolator2D(const double &x, const double &y, const double &z, const int nx_, const int ny_)

Construct from raw-pointer grid data.

Parameters
  • x – Pointer to the first x abscissa value.

  • y – Pointer to the first y abscissa value.

  • z – Pointer to the first grid value (row-major: z[i*ny + j]).

  • nx_ – Number of x grid points.

  • ny_ – Number of y grid points.

explicit Interpolator2D(const Interpolator2D &it)

Copy constructor (deep copy of the underlying GSL objects).

Parameters

it – Source interpolator.

explicit Interpolator2D()

Construct an empty (invalid) interpolator.

~Interpolator2D()

Destructor — releases GSL resources.

bool isValid() const

Return true if the interpolator holds valid data.

Returns

False for default-constructed objects.

void reset(const double &x, const double &y, const double &z, const int szx_, const int szy_)

Reinitialize with new grid data.

Parameters
  • x – Pointer to the first x abscissa.

  • y – Pointer to the first y abscissa.

  • z – Pointer to the first grid value.

  • szx_ – Number of x grid points.

  • szy_ – Number of y grid points.

double eval(const double &x, const double &y) const

Evaluate the bicubic spline at (x, y).

Parameters
  • x – x-coordinate.

  • y – y-coordinate.

Returns

Interpolated value.

Private Functions

void setup(const double &x, const double &y, const double &z, const int nx_, const int ny_)

Shared initialization helper.

Parameters
  • x – Pointer to the first x abscissa.

  • y – Pointer to the first y abscissa.

  • z – Pointer to the first grid value.

  • nx_ – Number of x grid points.

  • ny_ – Number of y grid points.

Private Members

const gsl_interp2d_type *TYPE = gsl_interp2d_bicubic

GSL 2D spline type (bicubic).

gsl_spline2d *spline

GSL 2D spline object.

gsl_interp_accel *xacc

GSL accelerator for the x axis.

gsl_interp_accel *yacc

GSL accelerator for the y axis.

size_t nx

Number of x grid points.

size_t ny

Number of y grid points.

class RootSolverBase
#include <numerics.hpp>

Base class for scalar root solvers.

Tracks iteration count, solver status, and the solution in a uniform way for all concrete root-solver implementations.

Subclassed by BrentRootSolver, SecantSolver

Public Functions

inline double getSolution() const

Return the last computed solution.

Returns

Root value found by the solver.

Protected Functions

inline RootSolverBase(const double &relErr_, const int maxIter_)

Construct with explicit tolerances.

Parameters
  • relErr_ – Relative error tolerance.

  • maxIter_ – Maximum iteration count.

inline explicit RootSolverBase()

Construct with default tolerances (1e-10, 1000).

Protected Attributes

const double relErr

Relative error tolerance for convergence.

const int maxIter

Maximum number of iterations.

int iter

Current iteration count.

int status

GSL solver status code (GSL_CONTINUE until converged).

double sol

Root found by the solver.

class BrentRootSolver : public RootSolverBase
#include <numerics.hpp>

Root solver using the GSL Brent bracketing algorithm.

Requires an initial bracket \([a, b]\) with \(f(a) \cdot f(b) < 0\).

Public Functions

explicit BrentRootSolver()

Construct the Brent solver with default tolerances.

~BrentRootSolver()

Destructor — releases GSL resources.

void solve(const std::function<double(double)> &func, const std::vector<double> &guess)

Find the root of func within the initial bracket guess.

Parameters
  • func – Function whose root is sought.

  • guess – Two-element vector \([a, b]\) bracketing the root.

Private Members

gsl_function *F

GSL function wrapper.

const gsl_root_fsolver_type *rst

GSL root-solver algorithm type.

gsl_root_fsolver *rs

GSL root-solver workspace.

class SecantSolver : public RootSolverBase
#include <numerics.hpp>

Root solver using the secant (chord) method.

Does not require a bracket but needs two initial guesses that are sufficiently close to the root.

Public Functions

inline SecantSolver(const double relErr_, const int maxIter_)

Construct with explicit tolerances.

Parameters
  • relErr_ – Relative error tolerance.

  • maxIter_ – Maximum iteration count.

inline explicit SecantSolver()

Construct with default tolerances.

void solve(const std::function<double(double)> &func, const std::vector<double> &guess)

Find the root of func starting from the two guesses in guess.

Parameters
  • func – Function whose root is sought.

  • guess – Two-element vector of initial approximations.

class Integrator1D
#include <numerics.hpp>

1D numerical integrator wrapping GSL quadrature routines.

Supports three quadrature modes (CQUAD for general integrands, QAWO for oscillatory Fourier-type integrals, QAGS for integrands with endpoint singularities) selected via the Type enum at construction time.

Public Types

enum Type

Quadrature mode selection.

  • DEFAULT: GSL CQUAD (general-purpose doubly adaptive).

  • FOURIER: GSL QAWO (weighted oscillatory integrals).

  • SINGULAR: GSL QAGS (endpoint-singularity handling).

Values:

enumerator DEFAULT
enumerator FOURIER
enumerator SINGULAR

Public Functions

Integrator1D(const Type &type, const double &relErr)

Construct with an explicit quadrature type and relative error.

Parameters
  • type – Quadrature mode.

  • relErr – Relative error tolerance.

inline explicit Integrator1D(const double &relErr)

Construct with the default (CQUAD) quadrature type.

Parameters

relErr – Relative error tolerance.

inline Integrator1D(const Integrator1D &other)

Copy constructor.

Parameters

other – Source integrator; type and accuracy are copied.

void compute(const std::function<double(double)> &func, const Param &param) const

Evaluate the integral of func over the domain specified by param.

Parameters
  • func – Integrand function \(f : \mathbb{R} \to \mathbb{R}\).

  • param – Integration parameters (limits or Fourier radius).

double getSolution() const

Return the result of the last compute() call.

Returns

Integral value.

inline double getAccuracy() const

Return the relative error tolerance.

Returns

Tolerance value.

inline Type getType() const

Return the quadrature type.

Returns

Type enum value.

Private Members

std::unique_ptr<Base> gslIntegrator

Pointer to the active GSL integrator backend.

class Base

Abstract base class for GSL-backed integrators.

Public Functions

inline Base(const Type &type_, const size_t &limit_, const double &relErr_)

Construct with workspace limit and accuracy.

Parameters
  • type_ – Quadrature mode.

  • limit_ – GSL workspace limit.

  • relErr_ – Relative error tolerance.

virtual ~Base() = default

Virtual destructor.

inline double getSolution() const

Return the integral value from the last compute.

Returns

Solution value.

inline double getAccuracy() const

Return the relative error tolerance.

inline Type getType() const

Return the quadrature type.

virtual void compute(const std::function<double(double)> &func, const Param &param) = 0

Compute the integral.

Parameters
  • func – Integrand.

  • param – Integration parameters.

Protected Attributes

gsl_function *F

GSL function wrapper (populated by subclasses).

const Type type

Quadrature type.

const size_t limit

GSL workspace limit.

const double relErr

Relative error tolerance.

double err

Residual absolute error from the last compute.

double sol

Integral value from the last compute.

class CQUAD : public Integrator1D::Base

GSL CQUAD doubly-adaptive quadrature.

Public Functions

explicit CQUAD(const double &relErr_)

Construct with a relative error tolerance.

Parameters

relErr_ – Relative error tolerance.

inline CQUAD(const CQUAD &other)

Copy constructor.

Parameters

other – Source CQUAD integrator.

~CQUAD()

Destructor — releases the GSL workspace.

void compute(const std::function<double(double)> &func, const Param &param) override

Compute the integral.

Parameters
  • func – Integrand.

  • param – Integration parameters.

Private Members

gsl_integration_cquad_workspace *wsp

GSL CQUAD workspace.

size_t nEvals

Number of integrand evaluations from the last compute.

class Param
#include <numerics.hpp>

Integration parameters.

Encapsulates the integration limits for finite intervals or the Fourier radius for QAWO-type integrals.

Subclassed by Integrator2D::Param

Public Functions

inline Param(const double &xMin_, const double &xMax_)

Construct for a finite-interval integral.

Parameters
  • xMin_ – Lower limit.

  • xMax_ – Upper limit.

inline explicit Param(const double &fourierR_)

Construct for a Fourier-type integral.

Parameters

fourierR_ – Fourier radius.

Public Members

const double xMin = numUtil::NaN

Lower integration limit (NaN for FOURIER type).

const double xMax = numUtil::NaN

Upper integration limit (NaN for FOURIER type).

const double fourierR = numUtil::NaN

Fourier integration radius (NaN for finite-interval types).

class QAGS : public Integrator1D::Base

GSL QAGS adaptive quadrature with singularity handling.

Public Functions

explicit QAGS(const double &relErr_)

Construct with a relative error tolerance.

Parameters

relErr_ – Relative error tolerance.

inline QAGS(const QAGS &other)

Copy constructor.

Parameters

other – Source QAGS integrator.

~QAGS()

Destructor — releases the GSL workspace.

void compute(const std::function<double(double)> &func, const Param &param) override

Compute the integral.

Parameters
  • func – Integrand.

  • param – Integration parameters.

Private Members

gsl_integration_workspace *wsp

GSL QAGS workspace.

class QAWO : public Integrator1D::Base

GSL QAWO oscillatory quadrature for Fourier-type integrals.

Public Functions

explicit QAWO(const double &relErr_)

Construct with a relative error tolerance.

Parameters

relErr_ – Relative error tolerance.

inline QAWO(const QAWO &other)

Copy constructor.

Parameters

other – Source QAWO integrator.

~QAWO()

Destructor — releases the GSL workspaces.

void compute(const std::function<double(double)> &func, const Param &param) override

Compute the integral.

Parameters
  • func – Integrand.

  • param – Integration parameters.

Private Members

gsl_integration_workspace *wsp

Primary GSL integration workspace.

gsl_integration_workspace *wspc

Cycle workspace for the QAWO algorithm.

gsl_integration_qawo_table *qtab

QAWO table storing the oscillation parameters.

class Integrator2D
#include <numerics.hpp>

2D numerical integrator using nested 1D GSL quadrature.

Computes \(\int_{x_{\min}}^{x_{\max}} f_1(x) \left[ \int_{y_{\min}(x)}^{y_{\max}(x)} f_2(y)\, dy \right] dx\) by applying an outer 1D integrator over \(x\) and an inner 1D integrator over \(y\), where the inner limits may depend on \(x\).

Public Types

using Type = Integrator1D::Type

Alias for the 1D integrator type enum.

using Param1D = Integrator1D::Param

Alias for the 1D parameter struct.

Public Functions

inline Integrator2D(const Type &type1, const Type &type2, const double &relErr)

Construct with independent types for the outer and inner integrators.

Parameters
  • type1 – Quadrature type for the outer (x) integral.

  • type2 – Quadrature type for the inner (y) integral.

  • relErr – Relative error tolerance for both integrators.

inline Integrator2D(const Type &type, const double &relErr)

Construct with the same quadrature type for both levels.

Parameters
  • type – Quadrature type.

  • relErr – Relative error tolerance.

inline Integrator2D(const double &relErr)

Construct using the default (CQUAD) quadrature type.

Parameters

relErr – Relative error tolerance.

void compute(const std::function<double(double)> &func1, const std::function<double(double)> &func2, const Param &param, const std::vector<double> &xGrid)

Compute the double integral.

Parameters
  • func1 – Outer integrand \(f_1(x)\).

  • func2 – Inner integrand \(f_2(y)\).

  • param – 2D integration parameters.

  • xGrid – Grid for grid-based outer integration.

inline double getX() const

Return the current outer-integration variable.

Returns

Most recent value of \(x\) used in the inner integral.

inline double getSolution() const

Return the result of the last compute() call.

Returns

Double-integral value.

Private Members

Integrator1D itg1

Outer (level 1) integrator over \(x\).

Integrator1D itg2

Inner (level 2) integrator over \(y\).

double x

Temporary storage for the current outer variable \(x\).

double sol

Integral result from the last compute() call.

class Param : public Integrator1D::Param
#include <numerics.hpp>

2D integration parameters.

Extends the 1D parameter with per-x lower and upper y-limits (which may be constant or x-dependent functions).

Public Types

using Func = std::function<double(double)>

Function type for x-dependent integration limits.

Public Functions

inline Param(const double &xMin_, const double &xMax_, const Func &yMin_, const Func &yMax_)

Construct with x-dependent y-limits.

Parameters
  • xMin_ – Lower x-limit.

  • xMax_ – Upper x-limit.

  • yMin_ – Lower y-limit function of x.

  • yMax_ – Upper y-limit function of x.

inline Param(const double &xMin_, const double &xMax_, const double &yMin_, const double &yMax_)

Construct with constant y-limits.

Parameters
  • xMin_ – Lower x-limit.

  • xMax_ – Upper x-limit.

  • yMin_ – Constant lower y-limit.

  • yMax_ – Constant upper y-limit.

inline Param(const double &fourierR_)

Construct for a Fourier-type 2D integral.

Parameters

fourierR_ – Fourier radius.

Public Members

const std::vector<double> xGrid

Optional x-grid for grid-based 2D integration.

Private Members

const double yMinNum = numUtil::NaN

Numeric constant for the lower y-limit.

const double yMaxNum = numUtil::NaN

Numeric constant for the upper y-limit.

namespace GslWrappers

Thin C++ wrappers around GSL objects.

Provides type-safe adapters for gsl_function (function objects) and error-checked wrappers for GSL allocation and call patterns.

Functions

template<typename Func, typename ...Args>
void callGSLFunction(Func &&gslFunction, Args&&... args)

Call a GSL function and propagate any GSL error as a C++ exception.

Template Parameters
  • Func – GSL function type.

  • Args – Argument types.

Parameters
  • gslFunction – GSL function to call.

  • args – Arguments forwarded to gslFunction.

template<typename Ptr, typename Func, typename ...Args>
void callGSLAlloc(Ptr &ptr, Func &&gslFunction, Args&&... args)

Allocate a GSL resource and propagate any allocation error.

Template Parameters
  • Ptr – Pointer type to receive the allocation result.

  • Func – GSL allocation function type.

  • Args – Argument types.

Parameters
  • ptr – Pointer to populate with the allocated resource.

  • gslFunction – GSL allocation function.

  • args – Arguments forwarded to gslFunction.

template<typename T>
class GslFunctionWrap : public gsl_function
#include <numerics.hpp>

Adapts a callable object to a gsl_function.

Template Parameters

T – Type of the callable (lambda, functor, or function pointer).

Public Functions

inline explicit GslFunctionWrap(const T &func_)

Construct the wrapper.

Parameters

func_ – Callable to wrap.

Private Members

const T &func

The wrapped callable.

Private Static Functions

static inline double invoke(double x, void *params)

Static trampoline invoked by GSL.

Parameters
  • x – Evaluation point.

  • params – Pointer to the GslFunctionWrap instance.

Returns

\(f(x)\).

namespace SpecialFunctions

Wrappers for commonly needed GSL special functions.

Functions

double fermiDirac12(const double &x)

Fermi–Dirac integral of order 1/2.

Parameters

x – Argument.

Returns

\(F_{1/2}(x) = \frac{1}{\Gamma(3/2)} \int_0^\infty \frac{t^{1/2}}{e^{t-x} + 1} dt\).

double fermiDiracm12(const double &x)

Fermi–Dirac integral of order -1/2.

Parameters

x – Argument.

Returns

\(F_{-1/2}(x) = \frac{1}{\Gamma(1/2)} \int_0^\infty \frac{t^{-1/2}}{e^{t-x} + 1} dt\).

double coth(const double &x)

Hyperbolic cotangent.

Parameters

x – Argument.

Returns

\(\coth(x)\).

double ellipticK(const double &x)

Complete elliptic integral of the first kind.

Parameters

x – Modulus \(k\).

Returns

\(K(k)\).

double ellipticE(const double &x)

Complete elliptic integral of the second kind.

Parameters

x – Modulus \(k\).

Returns

\(E(k)\).

double besselJ0(const double &x)

Bessel function of the first kind of order zero.

Parameters

x – Argument.

Returns

\(J_0(x)\).

class Vector2D
#include <vector2D.hpp>

Flat-storage 2D array with row-major layout.

Stores an \(s_1 \times s_2\) matrix as a contiguous std::vector<double> in row-major (C) order. Provides element access via operator()(i, j), row-span access via operator[](i), and arithmetic helpers used throughout the self-consistency iterations.

Public Functions

inline Vector2D(const size_t s1_, const size_t s2_)

Construct a zero-initialized matrix of size s1_ × s2_.

Parameters
  • s1_ – Number of rows.

  • s2_ – Number of columns.

inline explicit Vector2D()

Construct a 0 × 0 empty array.

explicit Vector2D(const std::vector<std::vector<double>> &v_)

Construct from a vector-of-vectors.

Parameters

v_ – Source 2D data; inner vectors must all have the same length.

explicit Vector2D(const std::vector<double> &v_)

Construct a 1 × n array from a flat vector.

Parameters

v_ – Source data.

size_t size() const

Return the total number of stored elements ( \(s_1 \times s_2\)).

Returns

Total element count.

size_t size(const size_t i) const

Return the size along dimension i.

Parameters

i – 0 for rows, 1 for columns.

Returns

Size of the requested dimension.

bool empty() const

Return true if the array has no elements.

Returns

True if \(s_1 = 0\) or \(s_2 = 0\).

void resize(const size_t s1_, const size_t s2_)

Resize the array, discarding existing data.

Parameters
  • s1_ – New row count.

  • s2_ – New column count.

double &operator()(const size_t i, const size_t j)

Mutable element access.

Parameters
  • i – Row index.

  • j – Column index.

Returns

Reference to element (i, j).

const double &operator()(const size_t i, const size_t j) const

Immutable element access.

Parameters
  • i – Row index.

  • j – Column index.

Returns

Const reference to element (i, j).

std::span<double> operator[](const size_t i)

Return a mutable span over row i.

Parameters

i – Row index.

Returns

Span covering the s2 elements of row i.

std::span<const double> operator[](const size_t i) const

Return an immutable span over row i.

Parameters

i – Row index.

Returns

Const span covering the s2 elements of row i.

bool operator==(const Vector2D &other) const

Element-wise equality comparison.

Parameters

other – Array to compare against.

Returns

True if all elements are equal.

std::vector<double>::iterator begin()

Return an iterator to the first element.

std::vector<double>::iterator end()

Return an iterator past the last element.

std::vector<double>::const_iterator begin() const

Return a const iterator to the first element.

std::vector<double>::const_iterator end() const

Return a const iterator past the last element.

double *data()

Return a mutable pointer to the underlying contiguous storage.

Returns

Pointer to the first element.

const double *data() const

Return an immutable pointer to the underlying contiguous storage.

Returns

Const pointer to the first element.

void fill(const double &num)

Set every element to num.

Parameters

num – Fill value.

void fill(const size_t i, const double &num)

Set every element of row i to num.

Parameters
  • i – Row index.

  • num – Fill value.

void fill(const size_t i, const std::vector<double> &num)

Set row i to a copy of num.

Parameters
  • i – Row index.

  • num – Source vector (must have length s2).

void sum(const Vector2D &v_)

Add v_ element-wise to this array (in-place).

Parameters

v_ – Array to add (must have the same shape).

void diff(const Vector2D &v_)

Subtract v_ element-wise from this array (in-place).

Parameters

v_ – Array to subtract (must have the same shape).

void mult(const Vector2D &v_)

Multiply element-wise by v_ (in-place).

Parameters

v_ – Array to multiply by (must have the same shape).

void mult(const double &num)

Multiply every element by the scalar num (in-place).

Parameters

num – Scaling factor.

void div(const Vector2D &v_)

Divide element-wise by v_ (in-place).

Parameters

v_ – Divisor array (must have the same shape).

void linearCombination(const Vector2D &v_, const double &num_)

Accumulate \(\text{num} \cdot v\_\) into this array (in-place).

Parameters
  • v_ – Array to scale and add.

  • num_ – Scaling factor.

Private Members

std::vector<double> v

Flat storage for all \(s_1 \times s_2\) elements (row-major).

size_t s1

Number of rows.

size_t s2

Number of columns.

class Vector3D
#include <vector3D.hpp>

Flat-storage 3D array with row-major layout.

Stores an \(s_1 \times s_2 \times s_3\) array as a contiguous std::vector<double> in row-major (C) order. Elements are accessed via operator()(i, j, k). Used primarily to store the fixed component of the qSTLS auxiliary density response.

Public Functions

inline Vector3D(const size_t s1_, const size_t s2_, const size_t s3_)

Construct a zero-initialized array of size s1_ × s2_ × s3_.

Parameters
  • s1_ – Size of the first dimension.

  • s2_ – Size of the second dimension.

  • s3_ – Size of the third dimension.

inline explicit Vector3D()

Construct a 0 × 0 × 0 empty array.

size_t size() const

Return the total number of stored elements.

Returns

\(s_1 \times s_2 \times s_3\).

size_t size(const size_t i) const

Return the size along dimension i.

Parameters

i – 0, 1, or 2 for the first, second, or third dimension.

Returns

Size of the requested dimension.

bool empty() const

Return true if the array has no elements.

Returns

True if any dimension size is zero.

void resize(const size_t s1_, const size_t s2_, const size_t s3_)

Resize the array, discarding existing data.

Parameters
  • s1_ – New size of the first dimension.

  • s2_ – New size of the second dimension.

  • s3_ – New size of the third dimension.

double &operator()(const size_t i, const size_t j, const size_t k)

Mutable element access.

Parameters
  • i – First-dimension index.

  • j – Second-dimension index.

  • k – Third-dimension index.

Returns

Reference to element (i, j, k).

const double &operator()(const size_t i, const size_t j, const size_t k) const

Immutable element access.

Parameters
  • i – First-dimension index.

  • j – Second-dimension index.

  • k – Third-dimension index.

Returns

Const reference to element (i, j, k).

bool operator==(const Vector3D &other) const

Element-wise equality comparison.

Parameters

other – Array to compare against.

Returns

True if all elements are equal.

std::vector<double>::iterator begin()

Return an iterator to the first element.

std::vector<double>::iterator end()

Return an iterator past the last element.

std::vector<double>::const_iterator begin() const

Return a const iterator to the first element.

std::vector<double>::const_iterator end() const

Return a const iterator past the last element.

double *data()

Return a mutable pointer to the underlying contiguous storage.

Returns

Pointer to the first element.

const double *data() const

Return an immutable pointer to the underlying contiguous storage.

Returns

Const pointer to the first element.

void fill(const double &num)

Set every element to num.

Parameters

num – Fill value.

void fill(const size_t i, const double &num)

Set every element with first index i to num.

Parameters
  • i – First-dimension index.

  • num – Fill value.

void fill(const size_t i, const size_t j, const double &num)

Set every element with indices (i, j, *) to num.

Parameters
  • i – First-dimension index.

  • j – Second-dimension index.

  • num – Fill value.

void fill(const size_t i, const size_t j, const std::vector<double> &num)

Set the row (i, j, *) to a copy of num.

Parameters
  • i – First-dimension index.

  • j – Second-dimension index.

  • num – Source vector (must have length s3).

void sum(const Vector3D &v_)

Add v_ element-wise to this array (in-place).

Parameters

v_ – Array to add (must have the same shape).

void diff(const Vector3D &v_)

Subtract v_ element-wise from this array (in-place).

Parameters

v_ – Array to subtract (must have the same shape).

void mult(const Vector3D &v_)

Multiply element-wise by v_ (in-place).

Parameters

v_ – Array to multiply by (must have the same shape).

void mult(const double &num)

Multiply every element by the scalar num (in-place).

Parameters

num – Scaling factor.

void div(const Vector3D &v_)

Divide element-wise by v_ (in-place).

Parameters

v_ – Divisor array (must have the same shape).

void linearCombination(const Vector3D &v_, const double &num_)

Accumulate \(\text{num} \cdot v\_\) into this array (in-place).

Parameters
  • v_ – Array to scale and add.

  • num_ – Scaling factor.

Private Members

std::vector<double> v

Flat storage for all \(s_1 \times s_2 \times s_3\) elements (row-major).

size_t s1

Size of the first dimension.

size_t s2

Size of the second dimension.

size_t s3

Size of the third dimension.

Functions

template<int Order>
Dual<Order> operator+(const Dual<Order> &dual1, const Dual<Order> &dual2)

Add two dual numbers.

Parameters
  • dual1 – Left operand.

  • dual2 – Right operand.

Returns

Element-wise sum.

template<int Order>
Dual<Order> operator+(const Dual<Order> &dual, const double &scalar)

Add a dual number and a scalar.

Parameters
  • dualDual operand.

  • scalar – Scalar to add to the primal.

Returns

Dual number with scalar added to the primal.

template<int Order>
Dual<Order> operator+(const double &scalar, const Dual<Order> &dual)

Add a scalar and a dual number.

Parameters
  • scalar – Scalar to add to the primal.

  • dualDual operand.

Returns

Dual number with scalar added to the primal.

template<int Order>
Dual<Order> operator-(const Dual<Order> &dual1, const Dual<Order> &dual2)

Subtract dual2 from dual1.

Parameters
  • dual1 – Left operand.

  • dual2 – Right operand.

Returns

Element-wise difference.

template<int Order>
Dual<Order> operator-(const Dual<Order> &dual, double scalar)

Subtract a scalar from a dual number.

Parameters
  • dual – Left operand.

  • scalar – Scalar to subtract from the primal.

Returns

Dual with scalar subtracted from the primal.

template<int Order>
Dual<Order> operator-(const double &scalar, const Dual<Order> &dual)

Subtract a dual number from a scalar.

Parameters
  • scalar – Left operand.

  • dual – Right operand.

Returns

Negated dual with scalar added to the primal.

template<int Order>
Dual<Order> operator*(const Dual<Order> &dual1, const Dual<Order> &dual2)

Multiply two dual numbers using the product rule.

Parameters
  • dual1 – Left operand.

  • dual2 – Right operand.

Returns

Product dual number.

template<int Order>
Dual<Order> operator*(const Dual<Order> &dual, const double &scalar)

Multiply a dual number by a scalar.

Parameters
  • dualDual operand.

  • scalar – Scalar multiplier.

Returns

Scaled dual number.

template<int Order>
Dual<Order> operator*(const double &scalar, const Dual<Order> &dual)

Multiply a scalar by a dual number.

Parameters
  • scalar – Scalar multiplier.

  • dualDual operand.

Returns

Scaled dual number.

template<int Order>
Dual<Order> operator/(const Dual<Order> &dual1, const Dual<Order> &dual2)

Divide dual1 by dual2 using the quotient rule.

Parameters
  • dual1 – Numerator.

  • dual2 – Denominator.

Returns

Quotient dual number.

template<int Order>
Dual<Order> operator/(const Dual<Order> &dual, const double &scalar)

Divide a dual number by a scalar.

Parameters
  • dual – Numerator.

  • scalar – Denominator scalar.

Returns

Scaled dual number.

template<int Order>
Dual<Order> operator/(const double &scalar, const Dual<Order> &dual)

Divide a scalar by a dual number.

Parameters
  • scalar – Numerator scalar.

  • dual – Denominator dual number.

Returns

Quotient dual number.

template<int Order>
Dual<Order> exp(const Dual<Order> &x)

Dual exponential \(e^x\).

Parameters

x – Argument.

Returns

Dual number with primal \(e^{x.func}\) and propagated gradient.

template<int Order>
Dual<Order> log(const Dual<Order> &x)

Dual natural logarithm \(\ln x\).

Parameters

x – Argument.

Returns

Dual number with primal \(\ln(x.func)\) and propagated gradient.

template<int Order>
Dual<Order> sqrt(const Dual<Order> &x)

Dual square root \(\sqrt{x}\).

Parameters

x – Argument.

Returns

Dual number with primal \(\sqrt{x.func}\) and propagated gradient.

template<int Order>
Dual<Order> tanh(const Dual<Order> &x)

Dual hyperbolic tangent \(\tanh x\).

Parameters

x – Argument.

Returns

Dual number with primal \(\tanh(x.func)\) and propagated gradient.

template<int Order>
class Dual
#include <dual.hpp>

Recursive template class for forward-mode automatic differentiation.

A Dual<Order> number carries a primal value and a gradient vector of Dual<Order-1> entries, enabling computation of derivatives up to the specified order. The base case Dual<1> stores a double primal and a std::vector<double> gradient.

Arithmetic operators (+, -, *, /) and common transcendental functions (exp, log, sqrt, tanh) are overloaded to propagate the chain rule automatically.

Convenience wrappers (Dual11, Dual21, Dual12, Dual22) are provided for the most common combinations of order and number of variables.

Template Parameters

Order – Differentiation order (≥ 1).

Public Functions

inline Dual(const Dual<Order - 1> &func_, const int nvar, const int index)

Construct from a lower-order dual primal with a seed index.

Parameters
  • func_ – Primal value.

  • nvar – Number of independent variables.

  • index – Index of the variable seeded with 1; -1 for no seeding.

inline Dual(const double &func_, const int nvar, const int index)

Construct from a scalar primal with a seed index.

Parameters
  • func_ – Scalar primal value.

  • nvar – Number of independent variables.

  • index – Index of the variable seeded with 1; -1 for no seeding.

Public Members

Dual<Order - 1> func

Primal value (itself a Dual<Order-1>).

std::vector<Dual<Order - 1>> grad

Gradient vector (one entry per independent variable).

template<>
class Dual<1>
#include <dual.hpp>

Base case: first-order dual number.

Subclassed by Dual11, Dual12

Public Functions

inline Dual(const double &func_, const int nvar, const int index)

Construct a first-order dual number.

Parameters
  • func_ – Scalar primal value.

  • nvar – Number of independent variables.

  • index – Index of the seeded variable; -1 for no seed.

Public Members

double func

Primal (function) value.

std::vector<double> grad

Gradient vector over all independent variables.

class Dual11 : public Dual<1>
#include <dual.hpp>

First-order dual number for a function of one variable.

Provides named accessors val() and dx().

Public Functions

inline explicit Dual11(const double func_, const int index = -1)

Construct a Dual11 number.

Parameters
  • func_ – Scalar primal value.

  • index – 0 to seed \(dx\); -1 for a constant.

inline Dual11(const Dual<1> &other)

Construct from a generic Dual<1>.

Parameters

other – Source dual number.

inline const double &val() const

Return the primal value \(f\).

inline const double &dx() const

Return the first derivative \(\partial f / \partial x\).

class Dual21 : public Dual<2>
#include <dual.hpp>

Second-order dual number for a function of one variable.

Provides named accessors val(), dx(), and dxx().

Public Functions

inline explicit Dual21(const double func_, const int index = -1)

Construct a Dual21 number.

Parameters
  • func_ – Scalar primal value.

  • index – 0 to seed \(dx\); -1 for a constant.

inline Dual21(const Dual<2> &other)

Construct from a generic Dual<2>.

Parameters

other – Source dual number.

inline const double &val() const

Return the primal value \(f\).

inline const double &dx() const

Return the first derivative \(\partial f / \partial x\).

inline const double &dxx() const

Return the second derivative \(\partial^2 f / \partial x^2\).

class Dual12 : public Dual<1>
#include <dual.hpp>

First-order dual number for a function of two variables.

Provides named accessors val(), dx(), and dy().

Public Functions

inline explicit Dual12(const double func_, const int index = -1)

Construct a Dual12 number.

Parameters
  • func_ – Scalar primal value.

  • index – 0 to seed \(dx\), 1 to seed \(dy\); -1 for a constant.

inline Dual12(const Dual<1> &other)

Construct from a generic Dual<1>.

Parameters

other – Source dual number.

inline const double &val() const

Return the primal value \(f\).

inline const double &dx() const

Return \(\partial f / \partial x\).

inline const double &dy() const

Return \(\partial f / \partial y\).

class Dual22 : public Dual<2>
#include <dual.hpp>

Second-order dual number for a function of two variables.

Provides named accessors val(), dx(), dy(), dxx(), dxy(), dyx(), and dyy().

Public Functions

inline explicit Dual22(const double &func_, const int index = -1)

Construct a Dual22 number.

Parameters
  • func_ – Scalar primal value.

  • index – 0 to seed \(dx\), 1 to seed \(dy\); -1 for a constant.

inline Dual22(const Dual<2> &other)

Construct from a generic Dual<2>.

Parameters

other – Source dual number.

inline const double &val() const

Return the primal value \(f\).

inline const double &dx() const

Return \(\partial f / \partial x\).

inline const double &dy() const

Return \(\partial f / \partial y\).

inline const double &dxx() const

Return \(\partial^2 f / \partial x^2\).

inline const double &dxy() const

Return \(\partial^2 f / \partial x \partial y\).

inline const double &dyx() const

Return \(\partial^2 f / \partial y \partial x\).

inline const double &dyy() const

Return \(\partial^2 f / \partial y^2\).

namespace MPIUtil

Utility functions for optional MPI-based parallel calculations.

Provides a thin abstraction layer over MPI that compiles to no-op stubs when MPI is disabled (USE_MPI not defined). The isUsed flag can be queried at runtime to determine whether MPI support is active.

Typedefs

using MPIParallelForData = std::vector<std::pair<int, int>>

Per-rank loop index range (inclusive start, exclusive end).

Used to distribute a loop of loopSize iterations across MPI ranks.

Functions

void init()

Initialize the MPI runtime.

void finalize()

Finalize the MPI runtime and release resources.

bool isInitialized()

Return true if MPI_Init has been called.

Returns

True if MPI is initialized.

int rank()

Return the rank of the calling process.

Returns

Rank in [0, numberOfRanks() - 1].

int numberOfRanks()

Return the total number of MPI processes in MPI_COMM_WORLD.

Returns

Process count.

void barrier()

Synchronize all MPI processes at a barrier.

bool isRoot()

Return true if the calling process is the root (rank 0).

Returns

True for rank 0.

bool isSingleProcess()

Return true if only one MPI process is running.

Returns

True if numberOfRanks() == 1.

void throwError(const std::string &errMsg)

Throw a runtime error with the given message.

Parameters

errMsg – Human-readable error description.

void abort()

Abort the MPI job immediately.

double timer()

Return the wall-clock time in seconds.

Returns

Elapsed wall time.

bool isEqualOnAllRanks(const int &myNumber)

Check that all ranks hold the same integer value.

Parameters

myNumber – Local integer value to compare across ranks.

Returns

True if all ranks agree on myNumber.

std::pair<int, int> getLoopIndexes(const int loopSize, const int thisRank)

Compute the loop index range for the calling rank.

Parameters
  • loopSize – Total number of loop iterations.

  • thisRank – Rank of the calling process.

Returns

Pair {start, end} for this rank’s slice.

MPIParallelForData getAllLoopIndexes(const int loopSize)

Compute loop index ranges for all ranks.

Parameters

loopSize – Total number of loop iterations.

Returns

Vector of {start, end} pairs, one per rank.

MPIParallelForData parallelFor(const std::function<void(int)> &loopFunc, const int loopSize, const int ompThreads)

Execute a parallel for loop across all MPI ranks.

Distributes loopSize iterations across ranks, invokes loopFunc on each rank’s slice (optionally parallelized with OpenMP), and returns the per-rank index ranges.

Parameters
  • loopFunc – Function to call with the loop index.

  • loopSize – Total number of iterations.

  • ompThreads – Number of OpenMP threads per rank.

Returns

Per-rank index ranges (needed for gatherLoopData).

void gatherLoopData(double *dataToGather, const MPIParallelForData &loopData, const int countsPerLoop)

Gather and synchronize data produced by a parallel for loop.

Each rank has computed countsPerLoop values for its slice of the loop. After this call dataToGather is consistent on all ranks.

Parameters
  • dataToGather – Pointer to the data buffer (size = loopSize × countsPerLoop).

  • loopData – Per-rank index ranges returned by parallelFor.

  • countsPerLoop – Number of values produced per loop iteration.

Variables

constexpr bool isUsed = false

True if the library was compiled with MPI support.

namespace vecUtil

Free-function utilities for element-wise operations on std::vector<double>.

All binary functions require the two input vectors to have the same length and return a new vector of that length.

Functions

std::vector<double> sum(const std::vector<double> &v1, const std::vector<double> &v2)

Compute the element-wise sum of two vectors.

Parameters
  • v1 – First operand.

  • v2 – Second operand (same length as v1).

Returns

Vector \(v_1 + v_2\).

std::vector<double> diff(const std::vector<double> &v1, const std::vector<double> &v2)

Compute the element-wise difference \(v_1 - v_2\).

Parameters
  • v1 – Minuend.

  • v2 – Subtrahend (same length as v1).

Returns

Vector \(v_1 - v_2\).

std::vector<double> mult(const std::vector<double> &v1, const std::vector<double> &v2)

Compute the element-wise product \(v_1 \cdot v_2\).

Parameters
  • v1 – First operand.

  • v2 – Second operand (same length as v1).

Returns

Vector \(v_1 \odot v_2\).

std::vector<double> div(const std::vector<double> &v1, const std::vector<double> &v2)

Compute the element-wise quotient \(v_1 / v_2\).

Parameters
  • v1 – Dividend.

  • v2 – Divisor (same length as v1; no zero-division check).

Returns

Vector \(v_1 \oslash v_2\).

std::vector<double> mult(const std::vector<double> &v, const double a)

Scale a vector by a scalar.

Parameters
  • v – Vector to scale.

  • a – Scalar factor.

Returns

Vector \(a \cdot v\).

std::vector<double> linearCombination(const std::vector<double> &v1, const double a, const std::vector<double> &v2, const double b)

Compute the linear combination \(a \cdot v_1 + b \cdot v_2\).

Parameters
  • v1 – First vector.

  • a – Coefficient for v1.

  • v2 – Second vector (same length as v1).

  • b – Coefficient for v2.

Returns

Linear combination vector.

double rms(const std::vector<double> &v1, const std::vector<double> &v2, const bool normalize)

Compute the root-mean-square difference between two vectors.

Parameters
  • v1 – First vector.

  • v2 – Second vector (same length as v1).

  • normalize – If true, divide the RMS by the RMS of v1 (relative error).

Returns

RMS (or relative RMS) difference.

void fill(std::vector<double> &v, const double &num)

Fill every element of v with the constant num.

Parameters
  • v – Vector to fill (modified in place).

  • num – Fill value.

namespace dimensionsUtil

Dimension-aware dispatch infrastructure for 2D and 3D calculations.

Provides an enumeration of supported spatial dimensions and a base class whose compute() method dispatches to dimension-specific virtual overrides.

Enums

enum class Dimension

Supported spatial dimensions. Default maps to 3D.

Values:

enumerator D3
enumerator D2
enumerator Default
class DimensionsHandler
#include <dimensions_util.hpp>

Base class that dispatches computations based on the spatial dimension.

Derived classes implement compute2D() and compute3D() to provide dimension-specific logic. Call compute(dim) to invoke the correct override at runtime.

Subclassed by ChemicalPotential, HFUtil::Idr, HFUtil::IdrGround, HFUtil::Ssf, HFUtil::SsfGround, Rdf, RpaUtil::Ssf, StlsIetUtil::Slfc, StlsUtil::Slfc, thermoUtil::Itcf, thermoUtil::ItcfNonInteracting, thermoUtil::ItcfNonInteractingGround

Public Functions

virtual ~DimensionsHandler() noexcept = default

Virtual destructor.

void compute(const Dimension &dim)

Dispatch to the appropriate dimension-specific compute method.

Parameters

dim – Spatial dimension to use for the computation.

Protected Functions

virtual void compute2D() = 0

Compute for a 2D system. Implemented by derived classes.

virtual void compute3D() = 0

Compute for a 3D system. Implemented by derived classes.

namespace databaseUtil

Utilities for SQLite database interaction used by the solvers.

Provides a metadata struct for identifying a run’s database entry and helper functions for managing blob data stored alongside the database.

Functions

void deleteBlobDataOnDisk(const std::string &dbInfo, int runId)

Delete all blob data on disk associated with a given run.

Parameters
  • dbInfo – Path to the SQLite database file.

  • runId – Run identifier whose blobs should be deleted.

bool blobDataTableExists(const std::string &dbName)

Check whether the blob-data table exists in a database.

Parameters

dbName – Path to the SQLite database file.

Returns

True if the table exists, false otherwise.

struct DatabaseInfo
#include <database.hpp>

Metadata identifying a simulation run in the SQLite database.

Each solver run that stores results in the database is assigned a unique integer runId within the table named runTableName. Large binary data (e.g., the fixed ADR component) may be stored outside the database in a directory given by blobStorage.

Public Functions

inline DatabaseInfo()

Construct with runId set to the sentinel integer NaN.

Public Members

std::string name

Path to the SQLite database file.

std::string blobStorage

Directory used to store large binary (blob) data on disk.

int runId

Integer run identifier within runTableName.

std::string runTableName

Name of the table that stores per-run metadata.

namespace numUtil

Numeric constants and comparison utilities for double-precision values.

Provides sentinel values for uninitialized parameters, a dimensionless constant used throughout the dielectric scheme calculations, and tolerance- based comparison helpers.

Functions

bool isZero(const double &x)

Test whether x is zero within the tolerance dtol.

Parameters

x – Value to test.

Returns

True if \(|x| \le \text{dtol}\).

bool equalTol(const double &x, const double &y)

Test whether two doubles are equal within the tolerance dtol.

Parameters
  • x – First value.

  • y – Second value.

Returns

True if \(|x - y| \le \text{dtol}\).

bool largerThan(const double &x, const double &y)

Test whether x is strictly greater than y (within dtol).

Parameters
  • x – Left-hand operand.

  • y – Right-hand operand.

Returns

True if \(x > y + \text{dtol}\).

Variables

constexpr double Inf = std::numeric_limits<double>::infinity()

Positive infinity sentinel.

constexpr double NaN = std::numeric_limits<double>::signaling_NaN()

Signaling NaN sentinel for uninitialized double parameters.

constexpr double iNaN = std::numeric_limits<int>::signaling_NaN()

Signaling NaN cast to int for uninitialized integer parameters.

constexpr double dtol = 1e-10

Tolerance used in floating-point comparisons.

constexpr double lambda3 = 4.0 / (9.0 * M_PI)

Cube of the Wigner–Seitz constant \(\lambda^3 = 4/(9\pi)\).

const double lambda = pow(lambda3, 1.0 / 3.0)

Wigner–Seitz constant \(\lambda = (4/(9\pi))^{1/3}\).

namespace formatUtil

Platform-portable string formatting wrapper.

On macOS, delegates to the {fmt} library because Apple’s C++ standard library does not ship std::format. On other platforms, std::format is used directly. The calling code can use formatUtil::format with the same syntax as std::format regardless of platform.

Functions

template<typename ...Args>
std::string format(std::format_string<Args...> fmt_str, Args&&... args)

Format a string using std::format.

Template Parameters

Args – Argument types to format.

Parameters
  • fmt_str – Format string (must be a compile-time constant).

  • args – Arguments to substitute into fmt_str.

Returns

Formatted std::string.