Implicit Solvent Models

Implicit solvation treats the solvent as a structureless polarizable continuum rather than individual molecules. The solute occupies a cavity inside the dielectric, and the reaction field of the polarized solvent feeds back into the solute's Hamiltonian.

OpenQuantum supports seven electrostatic models across two discretization families.

Theoretical Framework

The total free energy in an implicit solvent is

G = E_gas + ΔG_solv

where ΔG_solv is the solvation free energy, split into electrostatic (ΔG_es) and non-electrostatic (ΔG_non‑es) contributions. The electrostatic part comes from polarizing the dielectric continuum; the non-electrostatic part covers cavitation, dispersion, and solvent-structure effects (CDS).

Apparent Surface Charge (ASC) Formalism

All PCM-family models represent the reaction field through apparent surface charges q on the cavity boundary. The charges solve a linear system

K·q = R·v

where v is the electrostatic potential (nuclear + electronic) evaluated at each cavity point, and matrices K and R define the electrostatic model. The reaction-field energy is

E_RF = ½ qᵀ · v

and the Fock matrix contribution is the AO potential of the Gaussian-smeared surface charges.

Cavity Surface: SWIG and ISWIG

Two smooth-surface discretizations are available:

• SWIG (Switching/Gaussian): Uses a polynomial switching function to smoothly attenuate surface elements that lie inside neighbouring atoms' spheres. The Gaussian charge exponent ξ controls the spatial extent of each surface element.

• ISWIG (Improved SWIG): Uses an error-function-based switching (erf) that provides smoother behaviour at the cost of slightly more computation.

Both methods produce a set of tesserae (surface points) per atom, each carrying a Gaussian charge density (η/π)^(3/2) exp(-η|r − C|²) with η = ξ².

Electrostatic Models

C-PCM (Conductor-like PCM)

The simplest ASC model. Assumes the solvent is a perfect conductor (ε → ∞), then scales back to finite ε:

K = S,    R = −f·I,    f = (ε − 1) / ε

S is the Coulomb interaction matrix between Gaussian surface charges. Suitable for high-dielectric solvents (water, DMSO).

COSMO (Conductor-like Screening Model)

Similar to C-PCM but uses a different scaling factor:

f = (ε − 1) / (ε + ½)

Often more accurate for low-dielectric solvents.

IEF-PCM (Integral Equation Formalism PCM)

The full dielectric model including the normal-derivative matrix D:

K = S − f/(2π)·D·A·S,    R = −f·[I − 1/(2π)·D·A]

where A is the diagonal area matrix. IEF-PCM is the most general ASC model and reduces to C-PCM or COSMO in the appropriate limits.

SS(V)PE (Surface and Simulation of Volume Polarization for Electrostatics)

A symmetrized variant of IEF-PCM:

K = S − f/(4π)·(D·A·S + S·A·Dᵀ),    R = −f·[I − 1/(2π)·D·A]

SS(V)PE is the model of choice for general quantum-chemical applications with implicit solvation; it matches the exact solution for spherical cavities.

ddCOSMO / ddPCM (Domain Decomposition)

Domain-decomposition variants expand the reaction potential in real spherical harmonics on atomic van der Waals spheres rather than on a discretized cavity surface. The L coupling matrix replaces S/D, and the system size scales as natm × (lmax+1)²:

L·X = φ,    E_RF = ½·f_ε·ψ·X

• ddCOSMO: Conductor-like (COSMO) variant.
• ddPCM: Full dielectric variant using the A matrix.

Both are well-suited for large systems because the number of degrees of freedom depends only on the number of atoms and the spherical-harmonic order (lmax = 6 by default).

SMD (Solvation Model based on Density)

SMD layers a non-electrostatic CDS term on IEF-PCM electrostatics:

ΔG_solv = ΔG_es(IEF-PCM) + ΔG_cds

The CDS term uses atomic surface tensions parametrized for 96 solvents plus water. SMD requires a named solvent (not just a dielectric constant).

Analytical Gradients and Hessians

All four PCM-family models (C-PCM, COSMO, IEF-PCM, SS(V)PE) provide:

• Gradient: Three contributions — nuclear–cavity (grad_nuc), cavity response (grad_solver), and electron–cavity integral derivative (grad_qv). The total solvent gradient is added to the gas-phase gradient.

• Hessian: Three contributions — nuclear–cavity curvature (hess_nuc), cavity-response curvature (hess_solver), and electron–cavity integral curvature (hess_qv). The total solvent Hessian is added to the gas-phase Hessian.

  The analytical Hessian with solvent is available for C-PCM, COSMO, IEF-PCM, and SS(V)PE via freq=(analytical). The Hessian uses the solvated density from rhf_scf_with_reaction_field or uhf_scf_with_reaction_field and includes the CPHF solvent response (AXPCM kernel).

• ddCOSMO/ddPCM gradient: The analytical gradient uses the adjoint- Lagrangian formulation (ζ = L⁻¹·ψ). The l=0,m=0 block of the L matrix uses the exact Coulomb formula −1/|R_ia−R_ka| instead of the singular quadrature, making the gradient physically meaningful.

DFT-solvent coupling

Kohn–Sham DFT with implicit solvation adds the XC potential and the PCM reaction-field potential to the same effective Fock matrix. The SCF coupling is handled by dedicated entry points in the scf crate:

• rks_scf_with_reaction_field — Restricted Kohn–Sham with any PCM-family, ddCOSMO/ddPCM, or SMD solvent.
• uks_scf_with_reaction_field — Unrestricted Kohn–Sham with solvent.

At each SCF iteration the Fock matrix is

$$F_{\mu \nu }=H_{\mu \nu }^{\text{core}}+J_{\mu \nu }+c_x{K}_{\mu \nu }+V_{\mu \nu }^{\text{xc}}+V_{\mu \nu }^{\text{RF}},$$

where $V_{\mu \nu }^{\text{xc}}$ is the XC potential from grid quadrature and $V_{\mu \nu }^{\text{RF}}$ is the reaction-field potential from the solvent surface charges. The exact-exchange coefficient $c_x$ is zero for pure functionals and non-zero for hybrids and range-separated hybrids.

The total free energy in a solvent DFT calculation is

$$G=E_{\text{KS}}[\rho ]+E_{\text{RF}}+G_{\text{non-es}},$$

with $E_{\text{RF}}=\frac{1}{2}{\mathbf{q}}^{⊺}\mathbf{v}$ for PCM models and $G_{\text{non-es}}$ the SMD CDS correction.

Geometry optimization and frequency analysis in solvent

Optimizations in solvent recompute the cavity at each geometry step because the cavity surface is a function of the nuclear positions. The analytical gradient carries the full PCM contribution (grad_nuc + grad_solver + grad_qv) alongside the DFT gradient.

Harmonic frequencies with solvent use the solvated electronic density and include the PCM curvature contributions. The analytical Hessian is available for RHF/UHF/RKS/UKS with C-PCM, COSMO, IEF-PCM, and SS(V)PE. For ddCOSMO/ddPCM, finite-difference Hessians (freq=(numerical)) are the fallback.

Supported task combinations

• Post-HF in solvent uses the frozen solvent approximation (Frozen true in the SOLVATION section) by default.