Self-consistent reaction-field engine — owns the cavity, the COSMO/CPCM pre-factored response, and (optionally) the SMD CDS cavity. More...

#include <reaction_field.h>

Public Member Functions
	ReactionFieldEngine (Options opts={})

void	initialize (const Mat3N &atom_positions_bohr, const IVec &atomic_numbers)
	Build the cavity, pre-factor the COSMO A matrix, and (if `include_cds`) pre-evaluate the CDS branch.

void	solve_asc (const Vec &phi_at_cavity, std::optional< double > total_solute_charge=std::nullopt)
	Eulerian: solve `A σ = −f(ε)·φ_cavity`.

void	update_from_atom_charges (const Vec &atom_charges)
	Atom-resolved: from input atomic charges, compute φ = B·q, σ = G·q, V_atom = J_solv·q, E_es = ½ q·V_atom.

void	update_from_atom_multipoles (const Vec &atom_charges, const Mat3N &dipoles, const Mat &quadrupoles)
	Atom-resolved including the anisotropic part of the solute density.

void	set_multipole_damping (const Vec &rco_bohr, double kdmp3, double kdmp5)
	Short-range damping for the atom→cavity dipole and quadrupole kernels (see `detail::MultipoleDamping`).

double	energy () const
	Total cached solvation energy: E_es + E_cds.

double	energy_es () const

double	energy_cds () const

const Vec &	atom_potential () const
	Cached per-atom screening potential V_atom (Hartree, length = natom).

const Mat3N &	dipole_potential () const
	Conjugates of the atomic dipoles (3 × natom) and quadrupoles (6 × natom, in the input layout) at the last update.

const Mat &	quadrupole_potential () const

const Vec &	damping_radius_gradient () const
	∂E_es/∂R_co per atom at the last update (empty when the kernels are undamped).

const Vec &	surface_charges () const
	Cached apparent surface charges σ on the ES cavity.

SolvationSurfaces	surfaces () const
	Per-element decomposition of the latest update.

Mat3N	gradient () const
	Frozen-cavity analytical gradient (Hartree/Bohr, 3 × natom).

const occ::solvent::surface::Surface &	es_cavity () const

const occ::solvent::surface::Surface &	cds_cavity () const

size_t	num_es_surface_points () const

size_t	num_cds_surface_points () const

const Mat &	B () const

const Mat &	G () const

const Mat &	J_solv () const

double	dielectric () const

double	f_epsilon () const

const Options &	options () const

const occ::solvent::SMDSolventParameters &	smd_parameters () const

const Vec &	es_atom_radii () const
	ES cavity radii (Bohr) used at initialize(). Empty before initialize().

const Vec &	cds_energy_elements () const
	Per-element CDS energy (Hartree).

Detailed Description

Self-consistent reaction-field engine — owns the cavity, the COSMO/CPCM pre-factored response, and (optionally) the SMD CDS cavity.

Provides both an atom-resolved API (q → V_atom, used by xTB) and an Eulerian API (φ at cavity → σ, used by HF/DFT).

Lifecycle:

construct with Options,
initialize(atoms, Z) to build the cavity and factor A,
per-iteration: update_from_atom_charges(q) OR solve_asc(phi),
energy(), surfaces(), gradient() reflect the latest update.

gradient() is the frozen-cavity analytical gradient (plus FD CDS when include_cds). It requires an atom-resolved update (it goes through the SCC's Mulliken charges); calling it after a pure Eulerian solve returns an empty matrix.

Constructor & Destructor Documentation

◆ ReactionFieldEngine()

occ::scrf::ReactionFieldEngine::ReactionFieldEngine ( Options opts = {} )

explicit

Member Function Documentation

◆ atom_potential()

const Vec & occ::scrf::ReactionFieldEngine::atom_potential ( ) const

inline

Cached per-atom screening potential V_atom (Hartree, length = natom).

Empty until the first atom-resolved update.

◆ B()

const Mat & occ::scrf::ReactionFieldEngine::B ( ) const

inline

◆ cds_cavity()

const occ::solvent::surface::Surface & occ::scrf::ReactionFieldEngine::cds_cavity ( ) const

inline

◆ cds_energy_elements()

const Vec & occ::scrf::ReactionFieldEngine::cds_energy_elements ( ) const

inline

Per-element CDS energy (Hartree).

Length = ncav_cds; empty if include_cds == false. Stable across update* calls (geometry-only).

◆ damping_radius_gradient()

const Vec & occ::scrf::ReactionFieldEngine::damping_radius_gradient ( ) const

inline

∂E_es/∂R_co per atom at the last update (empty when the kernels are undamped).

Cut-off radii that depend on the geometry contribute to the nuclear gradient through this; the caller owns that chain because it owns the radii.

◆ dielectric()

double occ::scrf::ReactionFieldEngine::dielectric ( ) const

inline

◆ dipole_potential()

const Mat3N & occ::scrf::ReactionFieldEngine::dipole_potential ( ) const

inline

Conjugates of the atomic dipoles (3 × natom) and quadrupoles (6 × natom, in the input layout) at the last update.

Both are empty unless that update went through update_from_atom_multipoles.

◆ energy()

double occ::scrf::ReactionFieldEngine::energy ( ) const

inline

Total cached solvation energy: E_es + E_cds.

(E_cds = 0 when include_cds == false.)

◆ energy_cds()

double occ::scrf::ReactionFieldEngine::energy_cds ( ) const

inline

◆ energy_es()

double occ::scrf::ReactionFieldEngine::energy_es ( ) const

inline

◆ es_atom_radii()

const Vec & occ::scrf::ReactionFieldEngine::es_atom_radii ( ) const

inline

ES cavity radii (Bohr) used at initialize(). Empty before initialize().

◆ es_cavity()

const occ::solvent::surface::Surface & occ::scrf::ReactionFieldEngine::es_cavity ( ) const

inline

◆ f_epsilon()

double occ::scrf::ReactionFieldEngine::f_epsilon ( ) const

inline

◆ G()

const Mat & occ::scrf::ReactionFieldEngine::G ( ) const

inline

◆ gradient()

Mat3N occ::scrf::ReactionFieldEngine::gradient ( ) const

Frozen-cavity analytical gradient (Hartree/Bohr, 3 × natom).

Requires a prior atom-resolved update — returns a zero matrix otherwise. When include_cds, adds an FD CDS contribution on top.

This is the explicit derivative at frozen moments. After update_from_atom_multipoles the dipoles and quadrupoles are included in the source field, but their own response to the geometry (∂μ/∂R, ∂Θ/∂R) belongs to the caller's density-Pulay chain, exactly as it does for the charges.

◆ initialize()

void occ::scrf::ReactionFieldEngine::initialize	(	const Mat3N &	atom_positions_bohr,
		const IVec &	atomic_numbers
	)

Build the cavity, pre-factor the COSMO A matrix, and (if include_cds) pre-evaluate the CDS branch.

Idempotent on geometry change — call again whenever atoms move.

◆ J_solv()

const Mat & occ::scrf::ReactionFieldEngine::J_solv ( ) const

inline

◆ num_cds_surface_points()

size_t occ::scrf::ReactionFieldEngine::num_cds_surface_points ( ) const

inline

◆ num_es_surface_points()

size_t occ::scrf::ReactionFieldEngine::num_es_surface_points ( ) const

inline

◆ options()

const Options & occ::scrf::ReactionFieldEngine::options ( ) const

inline

◆ quadrupole_potential()

const Mat & occ::scrf::ReactionFieldEngine::quadrupole_potential ( ) const

inline

◆ set_multipole_damping()

void occ::scrf::ReactionFieldEngine::set_multipole_damping	(	const Vec &	rco_bohr,
		double	kdmp3,
		double	kdmp5
	)

Short-range damping for the atom→cavity dipole and quadrupole kernels (see detail::MultipoleDamping).

rco_bohr is one cut-off radius per atom; pass it empty to switch damping off. The radii belong to whatever model produced the moments — for GFN2 they are its own CN-dependent multipole radii — so the engine takes them rather than choosing them.

◆ smd_parameters()

const occ::solvent::SMDSolventParameters & occ::scrf::ReactionFieldEngine::smd_parameters ( ) const

inline

◆ solve_asc()

void occ::scrf::ReactionFieldEngine::solve_asc	(	const Vec &	phi_at_cavity,
		std::optional< double >	total_solute_charge = `std::nullopt`
	)

Eulerian: solve A σ = −f(ε)·φ_cavity.

Caches σ and φ so energy_es() / surfaces() reflect the result. Does NOT compute V_atom — use the atom-resolved path for that.

With total_solute_charge set, the solve is constrained to Σσ = −f(ε)·q, which for a conductor is exact by Gauss's law. The unconstrained solve falls short by the charge lying outside the cavity. Enforced with a Lagrange multiplier on the existing factorisation:

σ = σ_free − λ A⁻¹1,   λ = (1ᵀσ_free − σ_target) / (1ᵀA⁻¹1)

This restores the sum rule but does not redistribute σ; the shape error needs the Klamt–Jonas double-cavity correction, which is not implemented.

◆ surface_charges()

const Vec & occ::scrf::ReactionFieldEngine::surface_charges ( ) const

inline

Cached apparent surface charges σ on the ES cavity.

Length = ncav_es. Empty before the first update.

◆ surfaces()

SolvationSurfaces occ::scrf::ReactionFieldEngine::surfaces ( ) const

Per-element decomposition of the latest update.

Coulomb branch is populated whenever there has been an update; CDS branch is populated when include_cds (and reflects the cavity at initialize(), i.e. geometry-only).

◆ update_from_atom_charges()

void occ::scrf::ReactionFieldEngine::update_from_atom_charges ( const Vec & atom_charges )

Atom-resolved: from input atomic charges, compute φ = B·q, σ = G·q, V_atom = J_solv·q, E_es = ½ q·V_atom.

Updates all cached state; subsequent energy(), surfaces(), and gradient() reflect this update.

◆ update_from_atom_multipoles()

void occ::scrf::ReactionFieldEngine::update_from_atom_multipoles	(	const Vec &	atom_charges,
		const Mat3N &	dipoles,
		const Mat &	quadrupoles
	)

Atom-resolved including the anisotropic part of the solute density.

The cavity potential comes from the atom-centred multipole expansion

φ_i = Σ_A [ q_A/d + μ_A·d/d³ + Σ_αβ Θ^A_αβ d_α d_β/d⁵ ], d = r_i − R_A

so that a solute with no net atomic charges — benzene, say — still polarises the continuum through its quadrupoles. Caches the conjugate of every moment, each being that moment's own kernel contracted with σ:

V_atom = ∂E_es/∂q, V_dipole = ∂E_es/∂μ, V_quad = ∂E_es/∂Θ

dipoles is 3 × natom. quadrupoles is 6 × natom, traceless, ordered (xx, xy, yy, xz, yz, zz) with the off-diagonal components entering φ twice — the layout and normalisation of occ::xtb::CammMoments::qp, so the potentials come back conjugate to the stored components directly.

Costs one LU solve per call: the pre-solved response G only covers the charge column, and widening it to the dipoles and quadrupoles would make it nine times larger for no saving on the systems this is used for.

The documentation for this class was generated from the following file:

/home/runner/work/occ/occ/include/occ/scrf/reaction_field.h

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ ReactionFieldEngine()

Member Function Documentation

◆ atom_potential()

◆ B()

◆ cds_cavity()

◆ cds_energy_elements()

◆ damping_radius_gradient()

◆ dielectric()

◆ dipole_potential()

◆ energy()

◆ energy_cds()

◆ energy_es()

◆ es_atom_radii()

◆ es_cavity()

◆ f_epsilon()

◆ G()

◆ gradient()

◆ initialize()

◆ J_solv()

◆ num_cds_surface_points()

◆ num_es_surface_points()

◆ options()

◆ quadrupole_potential()

◆ set_multipole_damping()

◆ smd_parameters()

◆ solve_asc()

◆ surface_charges()

◆ surfaces()

◆ update_from_atom_charges()

◆ update_from_atom_multipoles()