Skip to content

Theory and implementation

This page gives a more detailed description of the complete active space self-consistent field (CASSCF) method and of the GPU-accelerated, atomic-orbital (AO) based formulation used in TeraChem. For the user-facing keywords, see the CASSCF overview. The empirical α-CASSCF correction is described on its own page.

The CASSCF energy

CASSCF is defined by the variational minimization of a configuration-interaction (CI) wavefunction

\[ |\Psi_{\mathrm{CI}}\rangle = \sum_I c_I |\Phi_I\rangle , \]

with respect to both the CI coefficients \(c_I\) and the molecular orbital (MO) coefficients. The \(|\Phi_I\rangle\) are Slater determinants (or spin-adapted combinations of determinants — see Configuration state functions and determinants). The MOs are parameterized through an antihermitian orbital-rotation operator \(\hat{\kappa}\),

\[ |\Psi_{\mathrm{CASSCF}}\rangle = e^{-\hat{\kappa}} |\Psi_{\mathrm{CI}}\rangle , \qquad \hat{\kappa} = \sum_{p<q} \kappa_{pq}\left(\hat{E}_{pq} - \hat{E}_{qp}\right) . \]

The orbitals are partitioned into three spaces:

  • Inactive (core) orbitals — doubly occupied in every determinant,
  • Active orbitals — variably occupied; all configurations consistent with the active space are generated, and
  • Virtual orbitals — unoccupied in every determinant.

In terms of the one-particle density matrix (OPDM) \(\boldsymbol{\gamma}\) and the two-particle density matrix (TPDM) \(\boldsymbol{\Gamma}\),

\[ \gamma_{pq} = \sum_{IJ} c_I c_J \langle\Phi_I|\hat{E}_{pq}|\Phi_J\rangle , \qquad \Gamma_{pqrs} = \tfrac{1}{2}\sum_{IJ} c_I c_J \langle\Phi_I|\hat{E}_{pq}\hat{E}_{rs} - \delta_{qr}\hat{E}_{ps}|\Phi_J\rangle , \]

the CASSCF energy takes the compact form

\[ E = \sum_{pq} \gamma_{pq}\,(p|\hat{h}|q) + \sum_{pqrs} \Gamma_{pqrs}\,(pq|rs) . \]

The orbital optimization is driven by the orbital gradient \(g_{pq} = \partial E / \partial \kappa_{pq} = 2\,(X_{pq} - X_{qp})\), where \(X_{pq} = \sum_r \gamma_{pr}(q|\hat{h}|r) + 2\sum_{rst}\Gamma_{prst}(qr|st)\) is the Lagrangian (generalized Fock) matrix.

AO-based, GPU-accelerated formulation

State-of-the-art algorithms for single-reference methods (HF, KS-DFT) have long exploited sparsity in the AO basis. TeraChem brings the same machinery to CASSCF: the entire orbital optimization is recast in terms of Coulomb (J) and exchange (K) matrices,

\[ J^{(\mathbf{D})}_{\mu\nu} = \sum_{\rho\sigma}(\mu\nu|\rho\sigma)D_{\rho\sigma} , \qquad K^{(\mathbf{D})}_{\mu\nu} = \sum_{\rho\sigma}(\mu\rho|\nu\sigma)D_{\rho\sigma} , \]

exactly the Fock-like contractions that the GPU HF/DFT code already builds efficiently. The inactive-core operator is built from the core density, and the active Fock matrix is built by back-transforming the active OPDM to the AO basis. The remaining half-transformed integrals \((\mu\nu|tu)\) — which would otherwise require an expensive integral transformation — are evaluated as a sum of Coulomb matrices over pseudodensities \(\mathbf{P}^{(tu)}_{\mu\nu} = C_{\mu t}C_{\nu u}\).

Two distinct sources of sparsity are exploited simultaneously:

  1. Sparsity in the AO electron-repulsion integrals, from the exponential decay of basis-function overlap, and
  2. Sparsity in the density and pseudodensity matrices, which becomes pronounced when the active orbitals are localized.

Together these reduce the formal \(O(N^4)\) scaling of the rate-limiting steps to roughly \(O(N^2)\) — and as low as near-linear for the half-transformed integrals of a localized chromophore in a large solvent bath — without introducing any new approximation to the theory or to the Coulomb operator. The Fock-like intermediates can additionally be built incrementally from the density change between macroiterations; in practice these incremental quantities are rebuilt in full every 8–10 iterations to flush accumulated screening errors. This is what allows TeraChem to run CASSCF on molecular systems of more than one thousand atoms.

Orbital optimization

The default optimizer uses an approximate, diagonal-Hessian update and neglects the coupling between the CI and orbital degrees of freedom. The approximate Hessian elements are assembled from quantities already formed while building the orbital gradient. Convergence behavior is controlled by the casscfmaxiter, casscfconvthre, and related keywords documented in the overview.

Analytic gradients and nonadiabatic couplings

Because the CASSCF energy is variational in both the CI and MO coefficients, analytic nuclear gradients are obtained by contracting the (already available) density matrices with derivative one- and two-electron integrals; the two-electron contribution is evaluated as derivative J and K matrices, so the gradient inherits the same GPU acceleration and AO-basis sparsity as the energy and scales approximately \(O(N^2)\).

For state-averaged CASSCF the orbitals minimize the state-averaged energy, not the energy of any individual state. Computing a state-specific gradient (or a nonadiabatic coupling vector) therefore requires the response of the wavefunction, obtained by solving the coupled-perturbed SA-CASSCF (CP-SA-CASSCF / CPMCSCF) equations. TeraChem solves these with a preconditioned conjugate-gradient (PCG) solver in which the preconditioner is an approximate diagonal Hessian; the relevant controls are the cpsacasscf* keywords in the overview. The nonadiabatic coupling (NAC) vector between states \(\Theta\) and \(\Lambda\) follows from essentially the same machinery, with the state-specific density replaced by a transition density and the orbital interstate coupling gradient replacing the orbital gradient. Each CPMCSCF iteration also scales approximately \(O(N^2)\).

These low-scaling analytic gradients and NAC vectors make minimum-energy conical intersection (MECI) searches and nonadiabatic molecular dynamics tractable for systems with hundreds of atoms.

Notation

Throughout the TeraChem documentation and the references below, a state-averaged CASSCF wavefunction is denoted SA-\(N\)-CAS(\(m\),\(n\)), where \(N\) is the number of states included in the average, \(m\) is the number of active electrons, and \(n\) is the number of active orbitals. For example, SA-2-CAS(2,2) averages two states over an active space of two electrons in two orbitals.

References

The GPU-accelerated, AO-based CASSCF method in TeraChem is described in:

  1. E. G. Hohenstein, N. Luehr, I. S. Ufimtsev, and T. J. Martínez, "An atomic orbital-based formulation of the complete active space self-consistent field method on graphical processing units," J. Chem. Phys. 142, 224103 (2015). doi:10.1063/1.4921956
  2. J. W. Snyder Jr., E. G. Hohenstein, N. Luehr, and T. J. Martínez, "An atomic orbital-based formulation of analytical gradients and nonadiabatic coupling vector elements for the state-averaged complete active space self-consistent field method on graphical processing units," J. Chem. Phys. 143, 154107 (2015). doi:10.1063/1.4932613