OpenQuantum

OpenQuantum is an ab initio molecular orbital (MO) calculation program written in Rust. It performs Hartree–Fock, Kohn–Sham DFT, and post–Hartree–Fock calculations within the Linear Combination of Atomic Orbitals (LCAO) framework. Supports RHF, UHF, ROHF, and RKS/UKS for closed- and open-shell molecules; MP2 and CCSD(T) correlation energies; density functional theory with GGA, meta-GGA, hybrid, and range-separated hybrid functionals; implicit solvation (PCM-family, ddCOSMO/ddPCM, SMD) with analytical gradients and Hessians; geometry optimization (BFGS, Berny RFO, primitive/DLC/TRIC internal coordinates) with analytic gradients; IRC and NEB pathways; harmonic frequency analysis with analytical Hessians and thermochemistry; and effective core potentials (ECPs).

Note: The project is under active development and not yet ready for production use.

Feature Highlights

Quick Navigation

• New to OpenQuantum? Start with Installation then Quick Start.
• Ready to run? See the Input File Format and Keyword Reference.
• Want to understand the algorithms? See Theoretical Background and Density Functional Theory.
• Running with solvent? See Implicit Solvent Models.
• Looking for the full API? See Developer Documentation.

Citation

If you use OpenQuantum in your research, please cite this preprint:

Pham, Le Nhan. OpenQuantum: A High-Performance Rust Implementation for Ab Initio Molecular Orbital Calculations. 2026.

Installation

Building from Source

Source code is not open yet. You can download binary files from https://github.com/lenhanpham/OpenQuantum-binary

# Clone the repository
git clone https://github.com/lenhanpham/OpenQuantum
cd OpenQuantum

# Build in release mode
cargo build --release

# Run all tests
cargo test

Note: The workspace contains multiple crates. The binary is named oquantum. When running the program, use -p oquantum to specify the package:

cargo run --release -p oquantum -- <input_file>

Requirements

• Rust 1.70+ (stable toolchain)
• Cargo (comes with Rust)
• Git (for cloning the repository)

Platform Support

OpenQuantum builds and runs on:

• Linux (x86_64, aarch64)
• macOS (Intel, Apple Silicon)
• Windows (x86_64)

Verifying the Installation

After building, test with a minimal input:

cat > test.inp << 'EOF'
task: RHF STO-3G

MOLECULE
  Charge        0
  Multiplicity  1
  Units         Angstrom
END

GEOMETRY
  H  0.0  0.0  0.0
  H  0.0  0.0  1.4
END
EOF

cargo run --release -p oquantum -- test.inp

You should see an RHF energy calculation output for H₂ with STO-3G basis.

Quick Start

OpenQuantum is a command-line program driven by plain text input files. The section-based format is the only supported format:

• Section-based: a task: header followed by MOLECULE, GEOMETRY, SCF, INT, DFT, SOLVATION, OPT, FREQ, MP2, CC, BASIS, SYMMETRY blocks terminated by END.

cargo run --release -p oquantum -- water_rhf.inp

Note: The binary is named oquantum. The -p oquantum is required because the workspace contains multiple crates.

Minimal RHF Energy on Water

task: RHF STO-3G

MOLECULE
  Charge        0
  Multiplicity  1
  Units         Angstrom
END

GEOMETRY
  O   0.000000   0.000000   0.117369
  H   0.756950   0.000000  -0.469476
  H  -0.756950   0.000000  -0.469476
END

Save this as water_rhf.inp and run:

cargo run --release -p oquantum -- water_rhf.inp

Anatomy of an Input File

Every input file consists of:

1. task: line — Required. Sets the run type, method, and basis set.

   task: RHF STO-3G
   

2. MOLECULE section — Optional. Charge, multiplicity, coordinate units.
3. GEOMETRY section — Required. Atom list with coordinates.
4. Optional sections — BASIS, SCF, INT, DFT, SOLVATION, OPT, FREQ, MP2, CC, SYMMETRY — only include when you need non-default settings.

Section names are case-insensitive. Each section ends with END (or end) on its own line. Lines inside a section are key value pairs (whitespace-separated, comment lines start with # or !, blank lines are ignored).

Next Steps

• Learn the Input Format in detail
• Browse Examples for common calculations
• Understand the Theoretical Background behind the methods

Environment Variables

Basis Search Path

When resolving an external or custom basis set (not built into the binary), OpenQuantum searches directories in this order:

1. Directory from the OPENQUANTUM_BASIS_PATH environment variable
2. basis/ subdirectory next to the running executable
3. basis/ subdirectory in the current working directory
4. Source-tree fallback: crates/basis/src/data/ (development only)

Files are named <basis_name_lowercase>.<ext> with the following extensions:

HPC Example

# SLURM job script
#SBATCH --job-name=oquantum
#SBATCH --ntasks=1

export OPENQUANTUM_SCRATCH=/scratch/$USER/$SLURM_JOB_ID
mkdir -p $OPENQUANTUM_SCRATCH

oquantum myjob.inp

The scratch directory is always printed at the start of each run:

  Scratch directory: /scratch/user/12345345

Debug Environment Variables

The Hartree–Fock Method

The Hartree–Fock (HF) method is a mean-field approximation to the many-electron Schrödinger equation. Each electron moves in an average field created by all other electrons, leading to a set of one-electron equations.

This chapter provides the theoretical foundation for the HF implementation in OpenQuantum, including the electronic Hamiltonian, the Roothaan–Hall equations, the Fock matrix construction, and total energy evaluation.

Overview of Chapters

• Electronic Hamiltonian — The Hamiltonian in atomic units
• Roothaan–Hall Equations — LCAO expansion and matrix equations
• Fock Matrix — Core Hamiltonian and two-electron contributions
• Total Energy — Energy expression and nuclear repulsion

UHF and ROHF

For open-shell systems, see Unrestricted & Restricted Open-Shell HF which covers separate α/β orbitals, spin contamination, and the ROHF formalism.

Key Equations

The central equations are summarized here for quick reference.

Electronic Hamiltonian (atomic units)

$$\hat{H}=\sum _i\hat{h}(i)+\sum _{i<j}\frac{1}{r_{ij}}$$

where:

• $\hat{h}(i)=-\frac{1}{2}{\nabla }_i^2-{\sum }_A\frac{Z_A}{r_{iA}}$ is the one-electron operator (kinetic + nuclear attraction)
• $1/r_{ij}$ is the electron–electron repulsion

Roothaan–Hall Equations

In the LCAO approximation, molecular orbitals are expanded in atomic orbitals:

$${\varphi }_i=\sum _{\mu }{C}_{\mu i}{\chi }_{\mu }$$

This leads to the matrix equations:

$$\mathbf{FC}=\mathbf{SC}\mathbf{\epsilon }$$

where:

• $\mathbf{F}$ is the Fock matrix
• $\mathbf{C}$ is the MO coefficient matrix
• $\mathbf{S}$ is the overlap matrix
• $\mathbf{\epsilon }$ is the diagonal matrix of orbital energies

Fock Matrix Elements

$$F_{\mu \nu }=H_{\mu \nu }+G_{\mu \nu }$$

where:

• $H_{\mu \nu }=T_{\mu \nu }+V_{\mu \nu }$ is the core Hamiltonian (kinetic + nuclear attraction)
• $G_{\mu \nu }={\sum }_{\lambda \sigma }{P}_{\lambda \sigma }\left[(\mu \nu |\lambda \sigma )-\frac{1}{2}(\mu \lambda |\nu \sigma )\right]$ is the two-electron contribution

The density matrix is:

$$P_{\mu \nu }=2\sum _{i(\text{occ})}{C}_{\mu i}{C}_{\nu i}$$

Total Energy

$$E=\frac{1}{2}\mathrm{Tr}[\mathbf{P}(\mathbf{H}+\mathbf{F})]+E_{\text{nuc}}$$

where $E_{\text{nuc}}={\sum }_{A<B}\frac{Z_A{Z}_B}{R_{AB}}$ is the nuclear repulsion energy.

The Electronic Hamiltonian

The electronic Hamiltonian in atomic units is:

$$\hat{H}=\sum _i\hat{h}(i)+\sum _{i<j}\frac{1}{r_{ij}}$$

where:

• $\hat{h}(i)=-\frac{1}{2}{\nabla }_i^2-{\sum }_A\frac{Z_A}{r_{iA}}$ is the one-electron operator (kinetic + nuclear attraction)
• $\frac{1}{r_{ij}}$ is the electron–electron repulsion

In the Born–Oppenheimer approximation, the nuclei are fixed and the electronic Hamiltonian acts only on electronic coordinates. The nuclear repulsion energy $E_{\text{nuc}}={\sum }_{A<B}\frac{Z_A{Z}_B}{R_{AB}}$ is a constant for a given geometry and is added to the electronic energy at the end.

Atomic Units

OpenQuantum uses atomic units throughout (Hartree for energy, Bohr for distance, electron mass for mass, $\hslash =1$). Conversion factors:

One-Electron Operator

The one-electron operator consists of kinetic energy and nuclear attraction:

$$\hat{h}(i)=\hat{T}(i)+\hat{V}(i)=-\frac{1}{2}{\nabla }_i^2-\sum _A\frac{Z_A}{r_{iA}}$$

In the LCAO basis, the matrix elements are:

$$h_{\mu \nu }=T_{\mu \nu }+V_{\mu \nu }$$

where $T_{\mu \nu }$ are kinetic energy integrals and $V_{\mu \nu }$ are nuclear attraction integrals, both computed via the Obara–Saika recurrence relations.

The Roothaan–Hall Equations

In the LCAO (Linear Combination of Atomic Orbitals) approximation, molecular orbitals are expanded in terms of atomic orbitals:

$${\varphi }_i=\sum _{\mu }{C}_{\mu i}{\chi }_{\mu }$$

This leads to the Roothaan–Hall matrix equations:

$$\mathbf{FC}=\mathbf{SC}\mathbf{\epsilon }$$

where:

• $\mathbf{F}$ is the Fock matrix
• $\mathbf{C}$ is the MO coefficient matrix
• $\mathbf{S}$ is the overlap matrix
• $\mathbf{\epsilon }$ is the diagonal matrix of orbital energies

Derivation

The HF energy for a single determinant is:

$$E=\sum _i⟨i|\hat{h}|i⟩+\frac{1}{2}\sum _{i,j}\left[⟨ij|ij⟩-⟨ij|ji⟩\right]$$

Minimizing $E$ with respect to orbital rotations subject to orthonormality constraints $⟨{\varphi }_i|{\varphi }_j⟩={\delta }_{ij}$ leads to the eigenvalue problem above.

Solution Procedure

The Roothaan–Hall equations are solved iteratively (SCF procedure):

1. Initial guess for $\mathbf{C}$ (core Hamiltonian, Hückel, SAD, or read from checkpoint)
2. Build density matrix: $P_{\mu \nu }=2{\sum }_{i(\text{occ})}{C}_{\mu i}{C}_{\nu i}$
3. Build Fock matrix using current density
4. Solve $\mathbf{FC}=\mathbf{SC}\mathbf{\epsilon }$ for new $\mathbf{C}$
5. Check convergence (density change, energy change, gradient norm)
6. Repeat from step 2 if not converged

Symmetric Orthogonalization

To solve the generalized eigenvalue problem, we transform to an orthonormal basis:

$$\mathbf{S}=\mathbf{Us}{\mathbf{U}}^{†}$$

$$\mathbf{X}=\mathbf{U}{\mathbf{s}}^{-1/2}{\mathbf{U}}^{†}$$

Then ${\mathbf{F}}^{\prime }={\mathbf{X}}^{†}\mathbf{FX}$ is diagonalized:

$${\mathbf{F}}^{\prime }{\mathbf{C}}^{\prime }={\mathbf{C}}^{\prime }\mathbf{\epsilon }$$

And the MO coefficients in the original basis are $\mathbf{C}=\mathbf{X}{\mathbf{C}}^{\prime }$.

Symmetry Adaptation

When point-group symmetry is enabled, the basis functions are symmetry-adapted, block-diagonalizing both $\mathbf{F}$ and $\mathbf{S}$ by irreducible representation. Each block is diagonalized separately, producing canonical-like orbital energies within each irrep. This is essential for post-HF calculations (MP2, CCSD(T)) which require symmetry-labeled orbitals.

The Fock Matrix

The Fock matrix elements are:

$$F_{\mu \nu }=H_{\mu \nu }+G_{\mu \nu }$$

where:

• $H_{\mu \nu }=T_{\mu \nu }+V_{\mu \nu }$ is the core Hamiltonian (kinetic + nuclear attraction)
• $G_{\mu \nu }={\sum }_{\lambda \sigma }{P}_{\lambda \sigma }\left[(\mu \nu |\lambda \sigma )-\frac{1}{2}(\mu \lambda |\nu \sigma )\right]$ is the two-electron contribution

The density matrix is:

$$P_{\mu \nu }=2\sum _{i(\text{occ})}{C}_{\mu i}{C}_{\nu i}$$

Core Hamiltonian

The core Hamiltonian matrix elements are computed directly from one-electron integrals:

$$H_{\mu \nu }=T_{\mu \nu }+V_{\mu \nu }$$

These integrals are evaluated using the Obara–Saika recurrence relations and are independent of the electron density.

Two-Electron Contribution

The two-electron part of the Fock matrix is built from the density matrix and electron repulsion integrals (ERIs):

$$G_{\mu \nu }=\sum _{\lambda \sigma }{P}_{\lambda \sigma }\left[(\mu \nu |\lambda \sigma )-\frac{1}{2}(\mu \lambda |\nu \sigma )\right]$$

The ERIs $(\mu \nu |\lambda \sigma )$ are computed using:

• McMurchie–Davidson algorithm (default, universal for all angular momenta)
• SP fast paths for S/P shells (up to 10× speedup)
• Rys quadrature for D/F+ shells
• Direct SCF mode for large systems (recomputes ERIs each iteration)

ERIs have 8-fold permutational symmetry:

$$(\mu \nu |\lambda \sigma )=(\nu \mu |\lambda \sigma )=(\mu \nu |\sigma \lambda )=(\nu \mu |\sigma \lambda )=(\lambda \sigma |\mu \nu )=(\sigma \lambda |\mu \nu )=(\lambda \sigma |\nu \mu )=(\sigma \lambda |\nu \mu )$$

Only unique ERIs are stored, and symmetry is exploited in all ERI-heavy computations (energy, gradients, Hessians).

Symmetry Acceleration

When point-group symmetry is enabled, only symmetry-unique shell quartets are computed. The detected point group is used to skip symmetry-equivalent ERI shell quartets during energy, gradient, and Hessian computations, and to assign irreducible representation labels (A₁, B₁, B₂, …) to molecular orbitals.

Total Energy

The total HF energy is:

$$E=\frac{1}{2}\mathrm{Tr}[\mathbf{P}(\mathbf{H}+\mathbf{F})]+E_{\text{nuc}}$$

where:

• $\mathbf{P}$ is the density matrix
• $\mathbf{H}$ is the core Hamiltonian matrix
• $\mathbf{F}$ is the Fock matrix
• $E_{\text{nuc}}={\sum }_{A<B}\frac{Z_A{Z}_B}{R_{AB}}$ is the nuclear repulsion energy

Derivation

Starting from the expectation value of the Hamiltonian for a single Slater determinant:

$$E=\sum _i⟨i|\hat{h}|i⟩+\frac{1}{2}\sum _{i,j}\left[⟨ij|ij⟩-⟨ij|ji⟩\right]$$

Substituting the LCAO expansion and density matrix definition yields the trace formula above.

Energy Components

The total energy is partitioned into:

The nuclear repulsion energy is a simple pairwise sum:

$$E_{\text{nuc}}=\sum _{A<B}\frac{Z_A{Z}_B}{R_{AB}}$$

where $Z_A$ are atomic numbers and $R_{AB}$ are inter-nuclear distances.

Energy Convergence

The SCF procedure monitors:

These can be tightened with the Conver keyword in the SCF section (Conver 10 for 10⁻¹⁰, etc.).

Unrestricted & Restricted Open-Shell HF

For open-shell systems, OpenQuantum supports both Unrestricted Hartree–Fock (UHF) and Restricted Open-Shell Hartree–Fock (ROHF).

UHF Equations

UHF uses separate spatial orbitals for α and β electrons:

$${\mathbf{F}}^{\alpha }=\mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\alpha }$$ $${\mathbf{F}}^{\beta }=\mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\beta }$$

where:

• $\mathbf{J}$ is the Coulomb matrix from total density ${\mathbf{P}}^T={\mathbf{P}}^{\alpha }+{\mathbf{P}}^{\beta }$
• ${\mathbf{K}}^{\alpha }$ and ${\mathbf{K}}^{\beta }$ are exchange matrices from respective spin densities

The density matrices are: $$P_{\mu \nu }^{\alpha }=\sum _{i(\text{occ}{)}^{\alpha }}{C}_{\mu i}^{\alpha }{C}_{\nu i}^{\alpha }$$ $$P_{\mu \nu }^{\beta }=\sum _{i(\text{occ}{)}^{\beta }}{C}_{\mu i}^{\beta }{C}_{\nu i}^{\beta }$$

The Fock matrices are coupled through the common Coulomb matrix $\mathbf{J}$.

UHF Energy

$$E=\sum _{\sigma =\alpha ,\beta }\mathrm{Tr}[{\mathbf{P}}^{\sigma }(\mathbf{H}+{\mathbf{F}}^{\sigma })]+E_{\text{nuc}}$$

Spin Contamination

The S² expectation value measures spin contamination:

$$⟨S^2⟩=S(S+1)+N_{\beta }-\sum _{i,j}|⟨{\varphi }_i^{\alpha }|{\varphi }_j^{\beta }⟩{|}^2$$

Ideal value: $S(S+1)$ where $S=(N_{\alpha }-N_{\beta })/2$

OpenQuantum reports $⟨S^2⟩$ in the output for every UHF calculation. Significant deviation from the ideal value indicates spin contamination (the determinant is not a pure spin eigenfunction).

ROHF

Restricted Open-Shell HF uses a single set of doubly-occupied orbitals plus singly-occupied orbitals. The Fock matrix is constructed to maintain orbital degeneracy within the open shell. ROHF avoids spin contamination but has a more complex Fock matrix construction.

At convergence, the occupied–occupied and virtual–virtual blocks of the Fock matrix are diagonalized separately (QCPsuC), producing canonical-like orbital energies required for subsequent post-HF calculations (MP2, CCSD(T)).

Unrestricted Hartree-Fock (UHF) Equations

For open-shell systems, UHF uses separate spatial orbitals for alpha ($\alpha$) and beta ($\beta$) electrons. This breaks spin symmetry but allows different spatial distributions for each spin.

Pople-Nesbet Equations

The UHF Fock equations are:

$$\begin{array}{rl}{\mathbf{F}}^{\alpha }{\mathbf{C}}^{\alpha }& =\mathbf{S}{\mathbf{C}}^{\alpha }{\mathbf{\epsilon }}^{\alpha }\\ {\mathbf{F}}^{\beta }{\mathbf{C}}^{\beta }& =\mathbf{S}{\mathbf{C}}^{\beta }{\mathbf{\epsilon }}^{\beta }\end{array}$$

where ${\mathbf{F}}^{\alpha }$, ${\mathbf{F}}^{\beta }$ are the alpha and beta Fock matrices, ${\mathbf{C}}^{\alpha }$, ${\mathbf{C}}^{\beta }$ are the molecular orbital coefficient matrices, $\mathbf{S}$ is the overlap matrix, and ${\mathbf{\epsilon }}^{\alpha }$, ${\mathbf{\epsilon }}^{\beta }$ are the orbital energy vectors.

Fock Matrix Construction

The Fock matrices are:

$$\begin{array}{rl}{\mathbf{F}}^{\alpha }& ={\mathbf{H}}_{\text{core}}+\mathbf{J}-{\mathbf{K}}^{\alpha }\\ {\mathbf{F}}^{\beta }& ={\mathbf{H}}_{\text{core}}+\mathbf{J}-{\mathbf{K}}^{\beta }\end{array}$$

Components

The key difference from RHF:

• RHF: $\mathbf{F}={\mathbf{H}}_{\text{core}}+2\mathbf{J}-\mathbf{K}$ (factor 2 from spin summation)
• UHF: $\mathbf{J}$ uses total density (no factor 2); exchange uses respective spin density

Density Matrices

The alpha and beta density matrices are:

$$\begin{array}{rl}{P}_{\mu \nu }^{\alpha }& =\sum _i^{N_{\alpha }}{C}_{\mu i}^{\alpha }{C}_{\nu i}^{\alpha }\\ P_{\mu \nu }^{\beta }& =\sum _i^{N_{\beta }}{C}_{\mu i}^{\beta }{C}_{\nu i}^{\beta }\end{array}$$

where $N_{\alpha }$ and $N_{\beta }$ are the number of alpha and beta electrons (occupied orbitals).

Electronic Energy

The total electronic energy is the sum of alpha and beta contributions:

$$E_{\text{elec}}=\frac{1}{2}\text{Tr}\left[{\mathbf{P}}^{\alpha }({\mathbf{H}}_{\text{core}}+{\mathbf{F}}^{\alpha })\right]+\frac{1}{2}\text{Tr}\left[{\mathbf{P}}^{\beta }({\mathbf{H}}_{\text{core}}+{\mathbf{F}}^{\beta })\right]$$

Expanding:

$$E_{\text{elec}}=\frac{1}{2}\sum _{\mu \nu }\left[P_{\mu \nu }^{\alpha }(H_{\mu \nu }+F_{\mu \nu }^{\alpha })+P_{\mu \nu }^{\beta }(H_{\mu \nu }+F_{\mu \nu }^{\beta })\right]$$

The total energy includes nuclear repulsion:

$$E_{\text{total}}=E_{\text{elec}}+E_{\text{nuc}}$$

SCF Iteration Procedure

1. Initial guess: Build initial ${\mathbf{P}}^{\alpha }$, ${\mathbf{P}}^{\beta }$ from core Hamiltonian or Hückel/SAD
2. Build Fock matrices: ${\mathbf{F}}^{\alpha }=\mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\alpha }$, ${\mathbf{F}}^{\beta }=\mathbf{H}+\mathbf{J}-{\mathbf{K}}^{\beta }$
3. Solve eigenproblems: Diagonalize ${\mathbf{F}}^{\alpha }{\mathbf{C}}^{\alpha }=\mathbf{S}{\mathbf{C}}^{\alpha }{\mathbf{\epsilon }}^{\alpha }$ and similarly for $\beta$
4. Form new densities: ${\mathbf{P}}_{\text{new}}^{\alpha }$, ${\mathbf{P}}_{\text{new}}^{\beta }$ from occupied orbitals
5. Check convergence: RMS density change $\Delta P=\sqrt{\frac{1}{N}{\sum }_{\mu \nu }(\Delta P_{\mu \nu }^{\alpha }{)}^2+(\Delta P_{\mu \nu }^{\beta }{)}^2}$
6. Mix/damp/DIIS: Apply acceleration if enabled
7. Repeat until convergence

Convergence Criteria

OpenQuantum uses the combined RMS density change:

$$\Delta P_{\text{RMS}}=\sqrt{\frac{{\sum }_{\mu \nu }(P_{\mu \nu }^{\alpha ,\text{new}}-P_{\mu \nu }^{\alpha ,\text{old}}{)}^2+(P_{\mu \nu }^{\beta ,\text{new}}-P_{\mu \nu }^{\beta ,\text{old}}{)}^2}{N_{\text{basis}}^2}}$$

Default threshold: $10^{-8}$.

Acceleration Methods

Implementation in OpenQuantum

Fock Build (crates/scf/src/fock.rs::build_uhf_fock_matrices)

#![allow(unused)]
fn main() {
pub fn build_uhf_fock_matrices(
    h_core: &DMatrix<f64>,
    density_alpha: &DMatrix<f64>,
    density_beta: &DMatrix<f64>,
    eri: &(impl IntegralProvider + Sync),
) -> (DMatrix<f64>, DMatrix<f64>) {
    // Total density P^T = P^α + P^β
    let total_density = density_alpha + density_beta;

    // Coulomb matrix from total density (same for both spins)
    let j = build_coulomb_matrix(&total_density, eri);

    // Exchange matrices (different for each spin)
    let k_alpha = build_exchange_matrix(density_alpha, eri);
    let k_beta  = build_exchange_matrix(density_beta, eri);

    // F^α = H + J - K^α
    let fock_alpha = h_core + &j - &k_alpha;

    // F^β = H + J - K^β
    let fock_beta = h_core + &j - &k_beta;

    (fock_alpha, fock_beta)
}
}

SCF Driver (crates/scf/src/uhf.rs::uhf_scf)

The main loop handles:

• Separate alpha/beta diagonalization with generalized_eigensolve
• DIIS with combined error vectors
• Level shifting on both Fock matrices
• Spin contamination calculation via compute_s_squared
• Convergence checking on combined RMS density change

Initial Guess

OpenQuantum supports three initial guesses for UHF:

1. Core Hamiltonian (default): Diagonalize $H_{\text{core}}$ in orthogonal basis; same MOs for both spins, different occupations
2. Extended Hückel: Approximate MO coefficients from atomic orbital parameters
3. SAD (Superposition of Atomic Densities): Build initial density from atomic fragments

Comparison with RHF

ROHF Connection

Restricted open-shell HF (ROHF) uses a single set of orbitals for both spins but with different occupation numbers. The Fock matrix is a linear combination:

$${\mathbf{F}}^{\text{ROHF}}={\mathbf{F}}^{\text{closed}}+\frac{1}{2}({\mathbf{F}}^{\alpha }+{\mathbf{F}}^{\beta })-\frac{1}{2}({\mathbf{K}}^{\alpha }+{\mathbf{K}}^{\beta })$$

ROHF maintains spin purity ($⟨S^2⟩=S(S+1)$) but often gives higher energies than UHF for open-shell systems.

Further Reading

• Pople & Nesbet, "Self-consistent orbitals for radicals" (1954)
• Szabo & Ostlund, Modern Quantum Chemistry, Ch. 3
• Helgaker, Jørgensen & Olsen, Molecular Electronic-Structure Theory, Ch. 10

Spin Contamination in UHF

In unrestricted Hartree-Fock (UHF), the single determinant wavefunction is not an eigenfunction of the ${\hat{S}}^2$ operator. This leads to spin contamination — the wavefunction mixes in higher-spin states.

The $⟨S^2⟩$ Expectation Value

For a UHF determinant with $N_{\alpha }$ alpha and $N_{\beta }$ beta electrons, the expectation value of the total spin operator is:

$$⟨S^2⟩=S(S+1)+N_{\beta }-\sum _{i=1}^{N_{\alpha }}\sum _{j=1}^{N_{\beta }}|⟨{\varphi }_i^{\alpha }|{\varphi }_j^{\beta }⟩{|}^2$$

where:

• $S=\frac{N_{\alpha }-N_{\beta }}{2}$ is the expected spin quantum number
• ${\varphi }_i^{\alpha }$ and ${\varphi }_j^{\beta }$ are the occupied alpha and beta molecular orbitals
• The overlap integral is $⟨{\varphi }_i^{\alpha }|{\varphi }_j^{\beta }⟩={\sum }_{\mu \nu }{C}_{\mu i}^{\alpha }{S}_{\mu \nu }{C}_{\nu j}^{\beta }$

Ideal vs. Actual Values

The actual $⟨S^2⟩$ from a UHF calculation equals the ideal value only when the alpha and beta occupied orbitals are perfectly orthogonal ($⟨{\varphi }_i^{\alpha }|{\varphi }_j^{\beta }⟩=0$ for all $i,j$). Any non-zero overlap between occupied alpha and beta orbitals increases $⟨S^2⟩$ above the ideal value.

Quantifying Spin Contamination

The degree of spin contamination is measured by the difference:

$$\Delta ⟨S^2⟩=⟨S^2{⟩}_{\text{UHF}}-⟨S^2{⟩}_{\text{ideal}}$$

• $\Delta ⟨S^2⟩\approx 0$: Pure spin state, minimal contamination
• $0<\Delta ⟨S^2⟩\lesssim 0.1$: Moderate contamination, usually acceptable
• $\Delta ⟨S^2⟩>0.1$: Significant contamination; results may be unreliable

In OpenQuantum, a warning is printed when $\Delta ⟨S^2⟩>0.01$:

Spin Information:
  N(alpha): 5
  N(beta):  4
  <S²>:     0.850000 (ideal: 0.750000)
  WARNING: Significant spin contamination detected (0.1000)

Physical Interpretation

Spin contamination arises because the UHF determinant allows alpha and beta electrons to occupy different spatial orbitals. This flexibility lowers the energy but mixes in higher-spin configurations:

• For a doublet ($S=1/2$), the contaminated wavefunction contains quartet ($S=3/2$) and higher-spin character
• For a triplet ($S=1$), it contains quintet ($S=2$) character
• The energy lowering is artificial — it comes from spin-symmetry breaking, not physical correlation

Consequences

1. Energies: UHF energies are variationally lower than RHF, but part of this lowering is due to spin contamination rather than physical correlation
2. Properties: Properties like dipole moments and gradients may be affected
3. Post-HF methods: MP2 and CCSD(T) built on a contaminated UHF reference inherit spin contamination
4. Geometry optimization: Contaminated forces can lead to incorrect geometries

Mitigation Strategies

Implementation in OpenQuantum

The computation is in crates/scf/src/fock.rs::compute_s_squared():

#![allow(unused)]
fn main() {
pub fn compute_s_squared(
    coefficients_alpha: &DMatrix<f64>,
    coefficients_beta: &DMatrix<f64>,
    overlap: &DMatrix<f64>,
    n_alpha: usize,
    n_beta: usize,
) -> (f64, f64) {
    // Ideal S²: S(S+1) where S = (N_α - N_β)/2
    let s = (n_alpha as f64 - n_beta as f64) / 2.0;
    let s_squared_ideal = s * (s + 1.0);

    // Overlap of occupied orbitals: S_αβ = C_occ^α^T · S · C_occ^β
    let c_alpha_occ = coefficients_alpha.columns(0, n_alpha);
    let c_beta_occ = coefficients_beta.columns(0, n_beta);
    let s_alpha_beta = c_alpha_occ.transpose() * overlap * c_beta_occ;

    // Sum of squared overlaps
    let overlap_sum: f64 = s_alpha_beta.iter().map(|x| x * x).sum();

    // ⟨S²⟩ = S_ideal + N_β - Σ|S_αβ|²
    let s_squared = s_squared_ideal + n_beta as f64 - overlap_sum;

    (s_squared, s_squared_ideal)
}
}

The result is returned in UhfResult and formatted in format_uhf_energy_output() with a warning for significant contamination.

Further Reading

• Szabo & Ostlund, Modern Quantum Chemistry, Ch. 3 (UHF theory)
• Krylov, "Spin-flip methods in quantum chemistry" (2006)
• Yamaguchi et al., "Spin-projected UHF energies" (1986)

Gaussian Basis Functions

OpenQuantum uses Cartesian Gaussian-type orbitals (GTOs):

$$\chi (\mathbf{r})=N{x}^{l_x}{y}^{l_y}{z}^{l_z}\exp (-\alpha r^2)$$

where:

• $N$ is the normalization constant
• $l_x+l_y+l_z=L$ is the angular momentum
• $\alpha$ is the Gaussian exponent

Contracted Gaussians

Contracted Gaussians are linear combinations of primitives:

$${\chi }_{\text{contracted}}=\sum _i{c}_i{\chi }_i({\alpha }_i)$$

The contraction coefficients $c_i$ and exponents ${\alpha }_i$ are read from basis set definition files (.obs format, compatible with EMSL BSE).

Normalization

Primitives are normalized so that:

$$\int |\chi (\mathbf{r}){|}^{2}d\mathbf{r}=1$$

The normalization constant for a Cartesian GTO is:

$$N={\left(\frac{2\alpha }{\pi }\right)}^{3/4}{\left[\frac{(2\alpha {)}^L(2l_x-1)!!(2l_y-1)!!(2l_z-1)!!}{(2L+1)!!}\right]}^{1/2}$$

Pure Spherical Harmonics

By default, OpenQuantum transforms Cartesian basis functions to pure spherical harmonics (5d, 7f, etc.). This reduces the number of basis functions for higher angular momentum:

Use 5d/7f (default) for pure harmonics, 6d/10f for Cartesian.

Embedded Basis Sets

OpenQuantum ships with over 100 basis sets via embedded .obs files:

See Basis Sets Reference for the complete list and usage.

Integral Evaluation

OpenQuantum computes all required molecular integrals using modern recurrence relations and efficient algorithms.

Overview of Chapters

• One-Electron Integrals — Overlap, kinetic, nuclear attraction
• Two-Electron Integrals — ERIs with MD, SP fast paths, Rys quadrature

One-Electron Integrals

One-electron integrals are evaluated using the Obara–Saika recurrence relations which provide an efficient and numerically stable way to compute integrals for arbitrary angular momentum.

Two-Electron Integrals

Electron repulsion integrals (ERIs):

$$(\mu \nu |\lambda \sigma )=\iint {\chi }_{\mu }({\mathbf{r}}_1){\chi }_{\nu }({\mathbf{r}}_1)\frac{1}{r_{12}}{\chi }_{\lambda }({\mathbf{r}}_2){\chi }_{\sigma }({\mathbf{r}}_2)d{\mathbf{r}}_{1}d{\mathbf{r}}_2$$

ERIs have 8-fold permutational symmetry. OpenQuantum provides multiple algorithms:

The selected ERI engine applies to all ERI computation paths: in-core integral build, direct SCF, analytical nuclear gradients, fully analytical Hessians, and semi-analytical (finite-difference of gradients) Hessians.

See Two-Electron Integrals for detailed algorithm descriptions.

One-Electron Integrals

OpenQuantum evaluates one-electron integrals using the Obara–Saika recurrence relations, which provide an efficient and numerically stable way to compute integrals for arbitrary angular momentum.

Overlap Integrals

$$S_{\mu \nu }=\int {\chi }_{\mu }(\mathbf{r}){\chi }_{\nu }(\mathbf{r})d\mathbf{r}$$

The overlap matrix $\mathbf{S}$ is used throughout the calculation for orthogonalization, density matrix formation, and energy evaluation.

Kinetic Energy Integrals

$$T_{\mu \nu }=\int {\chi }_{\mu }(\mathbf{r})\left(-\frac{1}{2}{\nabla }^2\right){\chi }_{\nu }(\mathbf{r})d\mathbf{r}$$

The kinetic energy operator in atomic units is $-\frac{1}{2}{\nabla }^2$. These integrals contribute to the core Hamiltonian $H_{\mu \nu }=T_{\mu \nu }+V_{\mu \nu }$.

Nuclear Attraction Integrals

$$V_{\mu \nu }=\int {\chi }_{\mu }(\mathbf{r})\sum _A\left(-\frac{Z_A}{|\mathbf{r}-{\mathbf{R}}_A|}\right){\chi }_{\nu }(\mathbf{r})d\mathbf{r}$$

The nuclear attraction potential is a sum over all nuclei $A$ with charge $Z_A$ at position ${\mathbf{R}}_A$. These integrals are computed efficiently using the Obara–Saika scheme with Rys quadrature for the $1/r$ operator.

Obara–Saika Recurrence

The Obara–Saika method generates higher angular momentum integrals from lower ones using recurrence relations. For a Gaussian product ${\chi }_a{\chi }_b$, the recurrence steps up angular momentum in each Cartesian direction:

$$\begin{array}{rl}(i+1)S_{i+1,j,k}& =R_{A,i}{S}_{i,j,k}+i\alpha S_{i-1,j,k}\\ & \text{etc. for }j,k\end{array}$$

where $R_A$ is the vector from the Gaussian center to nucleus $A$, and $\alpha$ is the Gaussian exponent.

This recurrence is numerically stable and avoids the catastrophic cancellation issues of earlier methods.

Two-Electron Integrals

Electron repulsion integrals (ERIs) are the computational bottleneck of Hartree–Fock and post-Hartree–Fock calculations:

$$(\mu \nu |\lambda \sigma )=\iint {\chi }_{\mu }({\mathbf{r}}_1){\chi }_{\nu }({\mathbf{r}}_1)\frac{1}{r_{12}}{\chi }_{\lambda }({\mathbf{r}}_2){\chi }_{\sigma }({\mathbf{r}}_2)d{\mathbf{r}}_{1}d{\mathbf{r}}_2$$

Permutational Symmetry

ERIs have 8-fold permutational symmetry:

$$(\mu \nu |\lambda \sigma )=(\nu \mu |\lambda \sigma )=(\mu \nu |\sigma \lambda )=(\nu \mu |\sigma \lambda )=(\lambda \sigma |\mu \nu )=(\sigma \lambda |\mu \nu )=(\lambda \sigma |\nu \mu )=(\sigma \lambda |\nu \mu )$$

Only unique ERIs are stored and computed. This reduces storage and computation by a factor of 8.

ERI Engine Selection

OpenQuantum provides multiple ERI algorithms, selectable via the INT section:

INT
  Eri  auto  # auto, md, sp, rys
END

The selected engine applies to all ERI computation paths:

• In-core integral build
• Direct SCF
• Analytical nuclear gradients
• Fully analytical Hessians
• Semi-analytical (FD) Hessians

McMurchie–Davidson Algorithm

The MD algorithm uses Hermite Gaussian intermediates and is universal for all angular momenta. It builds ERIs by:

1. Computing Hermite integrals over Gaussian primitives
2. Contracting to shell-pair level
3. Applying permutational symmetry

Complexity scales as $O(L^6)$ for maximum angular momentum $L$.

SP Fast Paths

For S and P shells (angular momentum $L\le 1$), analytical formulas give exact ERIs with far fewer operations. These are ~10× faster than MD for organic molecules. For higher angular momentum shells, MD is used as fallback.

Rys Quadrature

Rys quadrature is a numerical integration method that becomes competitive for D/F+ shells. OpenQuantum implements the Phenix-style Rys algorithm with precomputed roots and weights. MD is used as fallback when the Rys path is unavailable.

Schwarz Screening

For direct SCF and large systems, ERIs are pre-screened using the Cauchy–Schwarz inequality:

$$(\mu \nu |\lambda \sigma )\le \sqrt{(\mu \nu |\mu \nu )}\sqrt{(\lambda \sigma |\lambda \sigma )}$$

Only ERIs above the threshold Acc2E (default 10⁻¹²) are computed. This dramatically reduces the number of ERIs for large systems.

Analytic Nuclear Gradients

The RHF analytic gradient is:

$$\frac{\partial E}{\partial X_A}=\sum _{\mu \nu }{P}_{\mu \nu }\frac{\partial h_{\mu \nu }}{\partial X_A}+\frac{1}{2}\sum _{\mu \nu \lambda \sigma }{P}_{\mu \nu }{P}_{\lambda \sigma }\frac{\partial (\mu \nu |\lambda \sigma )}{\partial X_A}-\sum _{\mu \nu }{W}_{\mu \nu }\frac{\partial S_{\mu \nu }}{\partial X_A}+\frac{\partial E_{\text{nuc}}}{\partial X_A}$$

where $W=-{\mathbf{C}}_{\text{occ}}{\mathbf{\epsilon }}_{\text{occ}}{\mathbf{C}}_{\text{occ}}^{†}$ is the energy-weighted density matrix.

The UHF gradient has separate α and β components for the density and energy-weighted density matrices.

Gradient Components

Symmetry Acceleration

Analytic gradients use symmetry-accelerated ERI derivatives (skips symmetry-equivalent shell quartets). The detected point group is used to skip symmetry-equivalent ERI shell quartets during gradient computations, matching the GRAD2E SymShl approach.

This provides near-proportional speed-ups based on symmetry order (e.g., 2× for C₂, 4× for D₂ₕ) for quadratically-scaling steps.

Implementation

Gradient integrals are evaluated using the Obara–Saika recurrence for derivatives. The derivative of a Gaussian integral with respect to nuclear position is computed by differentiating the recurrence relations.

The gradient is used by:

• BFGS optimizer (Cartesian)
• Berny optimizer (internal coordinates)
• IRC path following
• NEB pathway optimization
• Geometry-space DIIS (GDIIS/GEDIIS)

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.

Geometry Optimization Algorithms

OpenQuantum provides a comprehensive suite of geometry optimization algorithms for minima and transition-state searches.

Overview of Chapters

• Primitive Internal Coordinates — Topology, B-matrix, Lindh Hessian, TRM step
• Delocalized & TRIC Coordinates — SVD rank detection, rigid-body removal
• GeomeTRIC Optimizer — Native TRIC quasi-Newton minimizer / RS-P-RFO TS engine
• IRC — Intrinsic Reaction Coordinate — Gonzalez–Schlegel mass-weighted predictor-corrector
• Transition-State Search — P-RFO, Bofill update, TS-specific trust radius
• Quasi-Newton Hessian Updates — MSP, BFGS, TS-BFGS, SR1, DFP, Bofill, PSB, MS, BFGS/Powell
• NEB — Nudged Elastic Band — Henkelman–Jónsson with climbing image
• Davidson Partial Eigensolver — Rayleigh–Ritz for lowest eigenvalues
• Minimum Mode Following — Alternative TS step with |λ| denominator
• Sella-Style Optimizer — TS-BFGS, sigma-trust, adaptive eigenvalue correction

Optimizer Selection

The optimizer is selected via the OPT section:

OPT
  Algorithm  bfgs | berny | rberny | sella | geometric
  Coord      cartesian | primitive | dlc | hdlc | tric | tric-p
  ...
END

Coordinate Models

Transition-State Search

Activate with TransitionState true or the TS task token. Two TS step engines are available, depending on the optimizer backend:

• P-RFO (Berny / Sella) partitions the Hessian into a TS mode (maximize) and minimization modes.
• RS-P-RFO (GeomeTRIC) is the restricted-step partitioned RFO variant; it applies a positive metric scaling and a damped Hebden iteration to keep the step inside the trust radius. See GeomeTRIC Optimizer for details.

GEDIIS Acceleration

Geometry-space DIIS (GDIIS/GEDIIS) extrapolates from previous geometry/gradient points. Enable with Diis true in the OPT section. Most effective for difficult optimizations near flat PES regions.

Orbital Warm-Start

Between optimization cycles the molecular geometry changes only slightly, so the converged orbitals from one cycle are an excellent initial guess for the SCF at the next. By default the optimizer carries the previous cycle's density matrix forward and uses it as the SCF guess, which typically reduces the SCF iteration count markedly after the first cycle. The first cycle has no prior orbitals and uses the configured SCF guess (core Hamiltonian by default).

This warm-start applies to RHF, UHF, and ROHF and to every optimizer backend and run mode (minimization, TS, IRC, NEB). Disable it with MOGuess false in the OPT section to force a fresh guess at every geometry.

Convergence Criteria

5-criterion Gaussian-style check:

• Energy change ΔE
• Gradient RMS
• Gradient max
• Displacement RMS
• Displacement max

Presets: Gauss default, Gauss loose, Gauss tight, or customize with Tighten N (tightens all by 10⁻ᴺ).

Primitive Internal Coordinates (IC)

When Coord primitive, Coord dlc, or Coord tric is requested, geometry steps are taken in internal-coordinate (IC) space using:

1. Topology Detection

• Covalent-radius bond detection (38 elements, Alvarez 2008 radii)
• Bond graph traversal for angles, dihedrals, out-of-plane bending (sp² centers), and supplementary linear-angle components for near-linear triplets

2. Wilson B-Matrix

Analytic Cartesian derivatives for all primitive types:

• Bonds, angles: standard Wilson formulae
• Dihedrals / out-of-plane: analytic 4-atom gradient (cross-product formula)
• Linear-angle: stateless perpendicular-axis construction + finite-difference gradient

The B-matrix relates Cartesian displacements to IC displacements: $$\Delta \mathbf{q}=\mathbf{B}\Delta \mathbf{x}$$

3. Lindh Model Hessian

Initial diagonal estimate in IC space:

• Bonds: $0.45\times \exp (-(r-r_0{)}^2)$ per Lindh
• Angles / linear-angle: $0.15$
• Dihedrals / out-of-plane: $0.005$

4. IC Back-Transform

Iterative Newton method $({\mathbf{B}}^T{\mathbf{G}}^{-1}\Delta \mathbf{q}$ loop) with periodic-angle wrapping; returns (new_xyz, bork) convergence flag.

5. Hessian Update

MSP (mixed-mode) update blending BFGS and Murtagh–Sargent with parameter $\varphi =1-(\Delta \mathbf{y}\cdot \mathbf{\Xi }{)}^2/((\Delta \mathbf{y}\cdot \Delta \mathbf{y})(\mathbf{\Xi }\cdot \mathbf{\Xi }))$; falls back to pure BFGS.

6. Step Engine

Trust-radius method (TRM) with level-shifted Newton step; Brent bisection to satisfy $\Vert \text{step}(\lambda )\Vert =\text{trust}$.

7. Convergence

5-criterion Gaussian-style check:

• Energy change ΔE
• Gradient RMS/max
• Displacement RMS/max

Thresholds from ConvergenceCriteria (gaussian_default, gaussian_loose, gaussian_tight).

8. Trust-Radius Update

ρ-based acceptance ratio; shrink on $\rho <0.25$, grow on $\rho >0.75$.

Delocalized & TRIC Coordinates

Delocalized Internal Coordinates (DLC)

DLC forms a non-redundant basis from the redundant primitive ICs:

1. SVD rank detection on the primitive B-matrix builds a non-redundant delocalized basis
2. Projection and back-transform operate in the reduced space

The DLC basis matrix $\mathbf{U}$ (rank × m) is obtained via SVD of the B-matrix: $$\mathbf{B}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{†}$$

Gradients and Hessians are transformed: $${\mathbf{g}}_{\text{DLC}}=\mathbf{U}{\mathbf{g}}_{\text{prim}}$$ $${\mathbf{H}}_{\text{DLC}}=\mathbf{U}{\mathbf{H}}_{\text{prim}}{\mathbf{U}}^{†}$$

The primitive Hessian ${\mathbf{H}}_{\text{IC}}$ is retained in full redundant IC space and updated via MSP as usual, providing stable curvature propagation across steps.

Translation-Rotation Internal Coordinates (TRIC)

TRIC adds rigid-body translational and rotational modes via weighted internal coordinates:

1. Modified Gram-Schmidt removes translational/rotational modes before and after DLC projection
2. Rigid-body components suppressed to machine precision

TRIC is required for:

• Rigid true option
• RemoveTr true option (also works with DLC/HDLC)

Native TRIC Coordinate System

The GeomeTRIC optimizer backend builds an explicit translation-rotation internal-coordinate system rather than removing TR modes from a DLC basis. For each disconnected fragment it appends six rigid-body coordinates in the order $[T_x,T_y,T_z,R_x,R_y,R_z]$:

• Translation coordinates are the fragment centroid components.
• Rotation coordinates are an exponential-map representation of the optimal rotation relative to the geometry captured when the system is built. The exponential map avoids gimbal-lock singularities, and its analytic derivatives supply the rotational rows of the Wilson $\mathbf{B}$-matrix.

Dihedral and rotation coordinates are periodic: internal-coordinate differences for these wrap into $(-\pi ,\pi ]$ during gradient projection and the internal-to-Cartesian back-transform. See the GeomeTRIC Optimizer chapter for the full step engine, initial-Hessian strategy, and trust-radius control.

HDLC (Hybrid Delocalized)

HDLC adds Cartesian components to the DLC basis, providing a complete spanning set while retaining the benefits of delocalization for the internal degrees of freedom.

Usage

OPT
  Coord  dlc      # or tric, hdlc, tric-p
  ...
END

TRIC is recommended for systems with floppy modes or when exact translational/ rotational separation is needed.

GeomeTRIC Optimizer

The GeomeTRIC backend is a quasi-Newton geometry optimizer that works natively in translation-rotation internal coordinates (TRIC). It targets both energy minima and first-order saddle points, drives an adaptive trust-radius loop, and back-transforms each internal-coordinate step to Cartesian space iteratively.

Select it from the OPT section:

OPT
  Algorithm  geometric      # aliases: geomet, tric
  Coord      tric | tric-p  # tric-p uses primitives only
  ...
END

Coordinate System

The optimizer builds a [TricSystem] from molecular topology:

1. Topology primitives — bonds, angles, dihedrals, and out-of-plane terms inferred from the connectivity graph. A linear/chain fallback is used when the topology yields no primitives.
2. Per-fragment external coordinates — for the tric model, each disconnected fragment contributes six rigid-body coordinates ordered as $[T_x,T_y,T_z,R_x,R_y,R_z]$:
   • Translation is the fragment centroid.
   • Rotation is encoded as an exponential map relative to the geometry captured when the system was built.

The tric-p model uses topology primitives only (no rigid-body coordinates), matching the Primitive coordinate kind.

Rotation Coordinates

Rotational internal coordinates use a quaternion/exponential-map algebra:

• The optimal rotation quaternion between the current fragment geometry and its reference is obtained from the eigenvector of a $4\times 4$ correlation matrix.
• The quaternion is mapped to a three-component exponential map that is free of gimbal singularities over the working range.
• Analytic derivatives of the exponential map with respect to Cartesian coordinates supply the rotational rows of the Wilson $\mathbf{B}$-matrix.

B-Matrix and Gradient Projection

The Wilson $\mathbf{B}$-matrix collects $\partial q_i/\partial x$ for every internal coordinate. Primitive rows come from the topology set; the centroid translation rows are constant ($1/N_{\text{frag}}$), and the rotation rows come from the analytic exponential-map derivatives.

Cartesian gradients are projected into internals through the Wilson $\mathbf{G}$-matrix:

$${\mathbf{g}}_q={\mathbf{G}}^{-1}\mathbf{B}{\mathbf{g}}_x,\mathbf{G}=\mathbf{B}{\mathbf{B}}^{†}$$

with ${\mathbf{G}}^{-1}$ formed as an SVD pseudo-inverse to handle the redundancy of the internal set.

Initial Hessian

initial_ic_hessian seeds the internal-coordinate Hessian:

• From a Cartesian Hessian (recommended for TS): transformed into internals via ${\mathbf{H}}_q={\mathbf{G}}^{-1}\mathbf{B}{\mathbf{H}}_x{\mathbf{B}}^{†}{\mathbf{G}}^{-1}$. For minimization it is forced positive-definite; for TS the single negative eigenvalue (the reaction coordinate) is preserved.
• Model guess (fallback): the Lindh diagonal for primitives plus weak positive curvature ($0.05$ a.u.) on each translation/rotation coordinate. This fallback also covers degenerate B-matrices from near-linear fragments.

Step Engine

Each cycle proposes a step with the trust-radius engine:

• Minimization uses the trust-radius Newton-Raphson (TRM) step $\Delta y=-(\mathbf{H}+v\mathbf{I}{)}^{-1}\mathbf{g}$ with BFGS Hessian updates by default.
• Transition states use restricted-step partitioned RFO (RS-P-RFO), which maximizes along the lowest Hessian eigenmode (transition vector) and minimizes along the rest, with Bofill Hessian updates by default.

The trust_step driver scales the proposed step to the current trust radius via a damped Hebden iteration. The iteration is bounded (oscillation detection plus a hard MAX_ITER cap) and returns the best step found. The RS-P-RFO metric uses a strictly positive scaling parameter; non-finite or non-positive values are clamped to a small positive floor so the symmetric eigensolver always converges.

Internal-to-Cartesian Back-Transform

ic_to_cartesian recovers Cartesian coordinates that reproduce a target internal displacement through the iteration

$$\Delta x={\mathbf{B}}^{†}{\mathbf{G}}^{-1}\Delta q$$

with periodic-aware residuals (dihedral and rotation coordinates wrap into $(-\pi ,\pi ]$), per-atom step capping, and adaptive damping. A bork flag is raised when the residual stays large, which the optimizer treats as a rejected step.

Trust-Radius Control

The step quality factor $\rho =\Delta E_{\text{actual}}/\Delta E_{\text{predicted}}$ drives an adaptive trust radius:

• $\rho >0.75$ at the trust boundary → grow trust (up to TMax).
• $\rho <0.25$ → shrink trust (down to the floor).
• Bad steps are rejected and retried with a halved trust radius; once the floor is reached the step is accepted to keep making progress.

Minimization rejects energy-raising steps against a confident prediction; TS searches legitimately raise the energy along the reaction coordinate and only reject pathological quality factors or failed back-transforms.

Convergence

Uses the same 5-criterion Gaussian-style test as the other backends (energy change, gradient RMS/max, displacement RMS/max). See Geometry Optimization Algorithms for the shared convergence presets.

IRC — Intrinsic Reaction Coordinate

The Gonzalez–Schlegel (1990) algorithm for mass-weighted IRC following:

Algorithm

1. Mass-weight the input Hessian: $${\mathbf{H}}_{\text{mw}}={\mathbf{M}}^{-1/2}\mathbf{H}{\mathbf{M}}^{-1/2}$$

2. Diagonalize to find the mode of largest imaginary frequency

3. Displace along that mode (forward and reverse) to get two starting points

4. For each direction: predictor (mass-weighted steepest descent) + corrector (constrained optimization orthogonal to the path tangent)

   • The corrector uses a proper Gonzalez–Schlegel corrector sub-loop (up to 10 iterations per IRC step)
   • The gradient is projected orthogonal to the radial (path tangent) direction in mass-weighted coordinates
   • A constrained steepest-descent step is taken on the hypersphere
   • Iterates until the perpendicular gradient norm is below $10^{-6}$

5. Convergence by max-gradient norm; output is a list of IrcPoint (arc, geometry, energy)

Bidirectional IRC

The IRC true directive triggers bidirectional IRC from the provided geometry. Restrict to one direction with IRCDir +1 (forward) or IRCDir -1 (reverse).

Requirements

• Starting geometry should be a pre-converged transition state (imaginary frequency expected)
• The initial Hessian must have at least one negative eigenvalue for the reaction mode

Usage

task: RHF STO-3G

GEOMETRY
  ...
END

OPT
  IRC     true
  IRCDir  +1
END

The output includes per-step energies, geometries, and a formatted path summary.

Transition-State Search (Berny Optimizer)

The Berny optimizer supports transition-state (TS) searches via the TS task token or by setting TransitionState true in the OPT section. The implementation uses Partitioned Rational Function Optimization (P-RFO) to locate first-order saddle points.

See also: the GeomeTRIC Optimizer provides an alternative restricted-step partitioned RFO (RS-P-RFO) TS engine in translation-rotation internal coordinates. Select it with Algorithm geometric plus TransitionState true.

Theory

A transition state is a first-order saddle point on the potential energy surface (PES) — a stationary point with exactly one negative eigenvalue of the Hessian matrix. The negative eigenvalue corresponds to the reaction coordinate (the direction connecting reactant and product).

The P-RFO method partitions the Hessian eigenvectors into two subspaces:

• TS mode (mode 1, lowest eigenvalue): maximize energy along this direction
• Minimisation modes (all other modes): minimize energy along these directions

The step is computed as:

$$\Delta \mathbf{x}=\sum _i{c}_i{\mathbf{v}}_i$$

where the coefficients $c_i$ are:

$$c_i=-\frac{g_i}{\lambda -{\mu }_i}\text{for }i<\text{ModMin}[\text{TS mode: upper RFO root}]$$ $$c_i=-\frac{g_i}{\lambda -{\mu }_i}\text{for }i\ge \text{ModMin}[\text{minimisation: lower RFO root}]$$

The RFO roots are:

$${\mu }_{\text{TS}}=\frac{1}{2}\left({\lambda }_{\text{max}}+\sqrt{{\lambda }_{\text{max}}^2+4g_{\text{max}}^2}\right)[\text{upper root for TS mode}]$$ $${\mu }_{\text{min}}\text{ found by bisection to satisfy: }\sum \frac{g_i^2}{\lambda -{\mu }_i}=\mu [\text{lower root for min modes}]$$

Implementation Details

The TS search in the Berny optimizer includes:

1. Initial TS Hessian (init_ts_hessian):

   • When no analytical Hessian is provided, the Lindh model Hessian is used
   • The softest IC mode (lowest eigenvalue) is set to −0.2 a.u.
   • This ensures P-RFO has a well-defined ascent direction from the first step

2. P-RFO Step (prfo_step in step_engine.rs):

   • Eigendecomposes the Hessian: $\mathbf{H}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{†}$
   • Sorts eigenvalues/eigenvectors by ascending eigenvalue
   • Projects gradient onto eigenvectors: ${\tilde{g}}_i={\mathbf{v}}_i\cdot \mathbf{g}$
   • Computes RFO roots for TS mode (upper) and minimisation modes (lower)
   • Forms step in eigenbasis, transforms back to IC space
   • Scales to trust radius if step is too large

3. Hessian Eigenvalue Check (ensure_correct_negative_eigenvalues):

   • Runs every 3 cycles during TS optimization
   • Checks that the Hessian has exactly 1 negative eigenvalue
   • If too few: flips smallest positive eigenvalues to −0.2 a.u.
   • If too many: flips least-negative eigenvalues to +0.2 a.u.

4. Bofill Hessian Update (update_hessian_bofill):

   • Blends PSB (Powell-Symmetric-Broyden) and MS (Murtagh-Sargent) updates
   • Weight: $\varphi =1-({\mathbf{s}}^{†}\cdot \mathbf{y}{)}^2/({\mathbf{s}}^{†}\cdot \mathbf{s}\cdot {\mathbf{y}}^{†}\cdot \mathbf{y})$
   • When $\varphi \to 1$: pure PSB (preserves curvature info)
   • When $\varphi \to 0$: pure MS (rank-1 correction)
   • Automatically adapts based on local surface curvature

5. TS-Specific Trust Radius:

   • Initial trust: 0.01 Å (vs 0.3 Å for minima)
   • Maximum trust: 0.03 Å (vs 0.3 Å for minima)
   • Conservative values prevent large steps near saddle points

6. TS-Specific DIIS Weighting:

   • Reaction coordinate gradient component scaled by factor 2.0
   • Minimum subspace size: 3 points (vs 2 for minima)
   • Cosine threshold: 0.5 (vs 0.0 for minima)

7. TS-Specific Output:

   • Prints Hessian negative eigenvalue count at each cycle
   • Header shows "Transition-State Optimization (Berny IC-space, P-RFO)"

TS Search Algorithm Flow

berny_optimize_ic_with_controls()
│
├─ Set TS-specific trust defaults (0.01, 0.03)
│
├─ Initialize Hessian
│   ├─ If cart_hessian provided: use it
│   └─ Else: guess_hessian_lindh()
│       └─ If transition_state: init_ts_hessian(h, 1, -0.2)
│
└─ Main loop
   │
   ├─ Compute IC gradient
   │
   ├─ TS diagnostics (print negative eigenvalue count)
   │
   ├─ Select step method
   │   ├─ If transition_state: prfo_step()
   │   └─ Else: trm_step() or gediis()
   │
   ├─ IC → Cartesian back-transform
   │
   ├─ Evaluate energy and gradient
   │
   ├─ Convergence check
   │
   ├─ Trust-radius update
   │
   ├─ Hessian update
   │   ├─ If transition_state: update_hessian_bofill()
   │   └─ Else: update_hessian_msp()
   │
   ├─ Periodic eigenvalue check (every 3 cycles)
   │   └─ ensure_correct_negative_eigenvalues(h, 1, 0.2)
   │
   └─ Accept step

Verification

A genuine TS shows exactly one imaginary frequency (printed as a negative value):

Frequencies (cm⁻¹):  -1247.3   1652.1   3825.4

Verify with a subsequent FREQ job on the converged TS geometry.

Quasi-Newton Hessian Update Methods

Four update formulas are available via HessianUpdateMethod, selected with the Update keyword in the OPT section:

OPT
  Update  bfgs | psb | ms | bofill | sr1 | dfp | bfgs_powell | ts-bfgs
END

MSP (Mixed Symmetric Powell — default for minimisations)

Adaptively blends BFGS and Murtagh–Sargent updates via a mixing parameter $\varphi =1-(\Delta \mathbf{y}\cdot \mathbf{\Xi }{)}^2/((\Delta \mathbf{y}\cdot \Delta \mathbf{y})(\mathbf{\Xi }\cdot \mathbf{\Xi }))$ where $\mathbf{\Xi }=\mathbf{H}\cdot \mathbf{s}-\mathbf{y}$:

$${\mathbf{H}}_{\text{new}}=\mathbf{H}+\frac{\mathbf{y}{\mathbf{y}}^{†}}{\mathbf{y}\cdot \mathbf{s}}-\frac{(\mathbf{Hs})(\mathbf{Hs}{)}^{†}}{{\mathbf{s}}^{†}\mathbf{Hs}}\text{(BFGS part)}$$ $$+\varphi \cdot \frac{\mathbf{\Xi }{\mathbf{\Xi }}^{†}}{\mathbf{\Xi }\cdot \mathbf{s}}\text{(SR1 correction)}$$

BFGS (Broyden–Fletcher–Goldfarb–Shanno)

The classic rank-2 update:

$$\Delta \mathbf{H}=\frac{\mathbf{y}{\mathbf{y}}^{†}}{\mathbf{y}\cdot \mathbf{s}}-\frac{(\mathbf{Hs})(\mathbf{Hs}{)}^{†}}{{\mathbf{s}}^{†}\mathbf{Hs}}$$

Preserves positive-definiteness when $\mathbf{y}\cdot \mathbf{s}>0$; unsuitable for TS searches.

TS-BFGS (Transition-State BFGS — default for saddle-point searches)

Standard BFGS drives the Hessian positive-definite, which destroys the negative curvature direction needed for TS optimization. TS-BFGS replaces the curvature metric ${\mathbf{s}}^{†}\mathbf{Hs}$ with ${\mathbf{s}}^{†}|\mathbf{H}|\mathbf{s}$ (using the absolute-value Hessian $|\mathbf{H}|$, obtained by diagonalizing $\mathbf{H}$ and flipping negative eigenvalues to their absolute values) so that the negative eigenvalue directions are maintained across steps.

The rank-2 single-step formula:

$$\mathbf{u}=\frac{(\mathbf{s}\cdot \mathbf{y})\mathbf{y}+(\mathbf{s}\cdot |\mathbf{H}|\mathbf{s})|\mathbf{H}|\mathbf{s}}{(\mathbf{s}\cdot \mathbf{y}{)}^2+(\mathbf{s}\cdot |\mathbf{H}|\mathbf{s}{)}^2}$$ $$\mathbf{J}=\mathbf{y}-\mathbf{Hs}$$ $$\Delta \mathbf{H}=\mathbf{u}{\mathbf{J}}^{†}+\mathbf{J}{\mathbf{u}}^{†}-(\mathbf{J}\cdot \mathbf{s})\mathbf{u}{\mathbf{u}}^{†}$$ $${\mathbf{H}}_{\text{new}}=\mathbf{H}+\Delta \mathbf{H}$$

where $\mathbf{s}={\mathbf{q}}_{n+1}-{\mathbf{q}}_n$ is the IC-space displacement and $\mathbf{y}={\mathbf{g}}_{n+1}-{\mathbf{g}}_n$ is the gradient change.

Properties:

• The quasi-Newton secant condition ${\mathbf{H}}_{\text{new}}\mathbf{s}\approx \mathbf{y}$ is approximately satisfied
• Negative eigenvalues are preserved, not destroyed
• Reduces exactly to standard BFGS when $\mathbf{H}$ is positive-definite

SR1 (Symmetric Rank-1)

$$\mathbf{\xi }=\mathbf{y}-\mathbf{Hs}$$ $$\Delta \mathbf{H}=\frac{\mathbf{\xi }{\mathbf{\xi }}^{†}}{\mathbf{\xi }\cdot \mathbf{s}}[\text{skipped when }|\mathbf{\xi }\cdot \mathbf{s}|<10^{-8}]$$

SR1 can build indefinite Hessian approximations naturally (no positive-definite constraint). Can be unstable if $\mathbf{\xi }\cdot \mathbf{s}$ is small, so a skip condition is applied.

DFP (Davidon–Fletcher–Powell)

$$\Delta \mathbf{H}=\frac{\mathbf{y}{\mathbf{y}}^{†}}{\mathbf{y}\cdot \mathbf{s}}-\frac{(\mathbf{Hs})(\mathbf{Hs}{)}^{†}}{{\mathbf{s}}^{†}\mathbf{Hs}}$$

DFP is the inverse-Hessian analogue of BFGS. It preserves positive-definiteness but can be slow to converge on ill-conditioned surfaces.

Bofill (Bofill Weighted Update — default for TS searches)

The Bofill update blends the Powell-Symmetric-Broyden (PSB) and Murtagh-Sargent (MS) updates with a weight $\varphi$ that measures how well the secant condition is satisfied:

$$\varphi =1-\frac{({\mathbf{s}}^{†}\mathbf{y}{)}^2}{{\mathbf{s}}^{†}\mathbf{s}\cdot {\mathbf{y}}^{†}\mathbf{y}}$$ $${\mathbf{H}}_{\text{new}}=\varphi \cdot \text{PSB}(\mathbf{H},\mathbf{s},\mathbf{y})+(1-\varphi )\cdot \text{MS}(\mathbf{H},\mathbf{s},\mathbf{y})$$

Properties:

• When $\varphi \to 1$: pure PSB (good for TS, preserves curvature information)
• When $\varphi \to 0$: pure MS (good for minima, rank-1 correction)
• Automatically adapts between PSB and MS based on the local surface curvature

PSB (Powell-Symmetric-Broyden)

$$\mathbf{\xi }=\mathbf{y}-\mathbf{Hs}$$ $${\mathbf{H}}_{\text{new}}=\mathbf{H}+\frac{\mathbf{\xi }{\mathbf{s}}^{†}+\mathbf{s}{\mathbf{\xi }}^{†}}{{\mathbf{s}}^{†}\mathbf{s}}-\frac{({\mathbf{\xi }}^{†}\mathbf{s})\mathbf{s}{\mathbf{s}}^{†}}{({\mathbf{s}}^{†}\mathbf{s}{)}^2}$$

PSB is a symmetric rank-2 update that does not preserve positive-definiteness. It is used as part of the Bofill update and is suitable for TS searches where the Hessian must remain indefinite.

MS (Murtagh-Sargent)

$$\mathbf{\xi }=\mathbf{y}-\mathbf{Hs}$$ $${\mathbf{H}}_{\text{new}}=\mathbf{H}+\frac{\mathbf{\xi }{\mathbf{\xi }}^{†}}{{\mathbf{\xi }}^{†}\mathbf{s}}$$

MS is a rank-1 update that can produce indefinite Hessians naturally. It is used as part of the Bofill update and is the simplest update that allows negative eigenvalues.

BFGS/Powell Mixed

if is_ts:
    H_new = Bofill(H,s,y)    # φ·PSB + (1-φ)·MS
else:
    H_new = BFGS(H,s,y)      # Standard BFGS

This method automatically selects BFGS for minima and Bofill for TS searches, providing a seamless transition between the two regimes.

Summary

NEB — Nudged Elastic Band

Henkelman–Jónsson (2000) NEB with climbing-image (CI-NEB):

Algorithm

1. Linear interpolation of $N$ images between start and end endpoints

2. Energy-weighted tangent selection: uphill images use max-energy neighbor, downhill use min-energy neighbor; blended at energy crossings

3. Spring force along tangent: $$F_{\text{spring}}=k_{\text{spring}}\times (|R_{i+1}-R_i|-|R_i-R_{i-1}|)\times \hat{\tau }$$

4. CI-NEB activates after $N/3$ iterations: highest-energy image climbs by inverting the gradient component along the tangent

5. Gradient-descent relaxation of the full image band; convergence by max force across all moveable images

6. Output includes per-image energies, geometries, and a formatted path summary

Usage

task: RHF STO-3G

GEOMETRY
  # Start geometry
  O   0.000000   0.000000   0.100000
  H   0.750000   0.000000  -0.450000
  H  -0.750000   0.000000  -0.450000

  # End geometry
  O   0.000000   0.000000   0.100000
  H   0.750000   0.100000  -0.450000
  H  -0.750000  -0.100000  -0.450000
END

OPT
  NEB     true
  Images  8
END

The start and end geometries are provided as consecutive blocks in the GEOMETRY section separated by a blank line. Images N controls the number of interpolating images (default 8).

Davidson Partial Eigensolver

A Rayleigh–Ritz iterative procedure for finding the k lowest eigenvalues and eigenvectors of a large matrix $\mathbf{A}$ without forming or storing it explicitly. Only matrix-vector products $\mathbf{Av}$ are required (supplied as a closure or via finite-difference gradient evaluations).

Algorithm

1. Initialise subspace $\mathbf{V}$ with $k$ unit vectors (or a user-supplied guess)
2. Build projected matrix $\tilde{\mathbf{A}}={\mathbf{V}}^{†}\mathbf{AV}$ (size $m\times m$, grows each step)
3. Diagonalise $\tilde{\mathbf{A}}\to$ Ritz values ${\theta }_i$ and Ritz vectors ${\mathbf{v}}_i=\mathbf{V}{\mathbf{e}}_i$
4. Check residuals ${\mathbf{r}}_i=\mathbf{A}{\mathbf{v}}_i-{\theta }_i{\mathbf{v}}_i$ for all $i<k$:
   • Converged when $\Vert {\mathbf{r}}_i\Vert \le \text{tol}\cdot \max (|{\theta }_i|,1)$ for all $i$
5. Expand: Orthogonalise the worst residual against $\mathbf{V}$ via Modified Gram-Schmidt (MGS); append the new column to $\mathbf{V}$ and $\mathbf{AV}$
6. Restart when $m\ge \text{max\_subspace}$: retain the $k$ best Ritz vectors

Modified Gram-Schmidt Orthogonalisation

Uses a sequential projection:

$$\mathbf{v}\leftarrow \mathbf{v}-(\mathbf{v}\cdot {\mathbf{w}}_j){\mathbf{w}}_j\text{for each basis column }{\mathbf{w}}_j$$ $$\mathbf{v}\leftarrow \mathbf{v}/\Vert \mathbf{v}\Vert$$

Returns None when the projected norm falls below tol (linear dependence).

Usage in OpenQuantum

The Davidson solver is used internally by refine_ts_hessian_via_davidson to determine an initial TS-mode direction from finite-difference Hessian-vector products, before the TS-BFGS update takes over.

It is also used by the Sella optimizer for adaptive eigenvalue correction.

Minimum Mode Following (MMF) Step

An alternative to P-RFO for saddle-point and TS searches that uses $|\lambda |$ in the denominator instead of RFO roots. This makes the step well-conditioned even when the Hessian eigenvalues have not yet converged to the correct sign pattern during the early steps of a TS search.

Formula

Given the eigendecomposition $\mathbf{H}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{†}$ (eigenvalues sorted ascending), the MMF step is:

$$c_i=+\frac{{\mathbf{v}}_i^{†}\mathbf{g}}{|{\lambda }_i|+\alpha }\text{for }i<\text{order}[\text{ascent along TS modes}]$$ $$c_i=-\frac{{\mathbf{v}}_i^{†}\mathbf{g}}{|{\lambda }_i|+\alpha }\text{for }i\ge \text{order}[\text{descent along min modes}]$$

$$\mathbf{s}=\sum _i{c}_i{\mathbf{v}}_i$$

The Lagrange multiplier $\alpha \ge 0$ is determined by bisection to satisfy $\Vert \mathbf{s}(\alpha )\Vert =\text{trust}$. When the unconstrained step already satisfies the trust radius, $\alpha =0$ is used directly.

Comparison with P-RFO

Usage

Enable with Step mmf in the OPT section:

OPT
  TransitionState true
  Step     mmf
  ...
END

When the Hessian eigenvalue signs are not yet correct (e.g., very early steps on a flat surface), step=mmf is more robust than P-RFO because MMF uses $|\lambda |$ in the denominator.

Sella-Style Optimizer

The Sella optimizer (Algorithm sella or Sella task token) provides a quasi-Newton geometry optimization backend with explicit Cartesian and internal-coordinate execution paths, combining:

1. TS-BFGS Hessian update (maintains order negative eigenvalues)

2. P-RFO or MMF step computation

3. Sigma-based trust-radius schedule: multiplicative update driven by the model-quality ratio $\rho =\Delta E_{\text{actual}}/\Delta E_{\text{predicted}}$:

   Condition	Action
   $\rho \ge \text{rho\_inc}$	$\delta \times ={\sigma }_{\text{inc}}$ (grow)
   $\rho <1/\text{rho\_dec}$	$\delta \times ={\sigma }_{\text{dec}}$ (shrink)
   otherwise	unchanged

4. Adaptive Hessian eigenvalue correction: every nsteps_per_diag steps, the Hessian eigenspectrum is inspected and corrected if the number of negative eigenvalues deviates from order:

   • Too few negative eigenvalues → flip the smallest positive ones to −0.2 a.u.
   • Too many negative eigenvalues → flip the least-negative ones to +0.2 a.u.

5. Initial TS Hessian: TS startup uses Davidson partial eigensolver refinement when available; if refinement is unavailable, fallback initialization sets the order softest modes to −0.2 a.u. so that P-RFO/MMF has a well-defined ascent direction from the first step.

6. Coordinate-space support:

   • Coord cartesian: native Cartesian Sella engine path
   • Coord primitive: native internal-coordinate Sella path
   • Coord dlc / Coord tric: reduced-space internal step path

Default Parameters

Relationship to Existing Optimizers

Usage

Sella minimum optimisation:

task: SELLA RHF STO-3G

GEOMETRY
  ...
END

Sella TS search with MMF step (more robust on flat surfaces):