CASSCF

TeraChem can carry out complete active space self-consistent field calculations, which allow the wavefunction to comprise multiple determinants. The CASSCF method is especially useful for treating multireference problems, such as diradicals, bond-breaking, and electronic excited states. The disadvantages of CASSCF are its lack of attention to dynamic electron correlation (this generally leads to too large S0-S1 excitation energies and to a disfavoring of charge transfer states) and the need to specify an active space. The active space is specified by giving a number of electrons and a number of orbitals, e.g. CAS(2/4) for an active space with 2 electrons in 4 orbitals. All configurations consistent with the active space are included in the calculation.

TeraChem implements CASSCF with a GPU-accelerated, atomic-orbital-based algorithm that exploits sparsity to reach systems of more than a thousand atoms, including analytic gradients and nonadiabatic coupling vectors. For a detailed description of the method and its implementation, see Theory and implementation. The empirical α-CASSCF correction can be layered on top of any state-averaged calculation to approximate MS-CASPT2 surfaces at SA-CASSCF cost. A closely related, cheaper method that skips the orbital optimization — CASCI and FOMO-CASCI — reuses most of the same keywords.

Configuration state functions and determinants

The CASSCF wavefunction is written as \(\Psi(\Theta) = \sum_I c_I \psi_I(\Theta)\) where \(\Theta\) is a shorthand for the set of molecular orbitals (which are optimized), \(\psi_I\) consists of either Slater determinants or spin-adapted combinations of Slater determinants (CSFs), and \(c_I\) are the configuration interaction coefficients corresponding to the many-electron basis functions.

TeraChem always works with a Slater determinant basis at the fundamental level. However, you may choose to transform the result to a spin-adapted CSF basis. The advantages of working with CSFs are that the basis set is smaller (so the calculation is faster and converges more readily) and there is no ambiguity about spin states. When using determinants, the code can find electronic states of arbitrary spin. When using CSFs, the user specifies the desired spin state (e.g. singlet or triplet) and only states with the desired spin are located.

State-averaging

When searching for excited states, it is often convenient to use state-averaging. The basic idea here is to find orbitals that minimize the weighted average of the energies of a set of electronic states. This ensures that the orbitals are not biased towards any of the given states. The primary advantage is not accuracy (it is widely believed that state-specific CASSCF is more accurate than state-averaged CASSCF, although this certainly depends on the property being computed), but rather stability when one is computing excited states. State-specific CASSCF for excited states is often difficult to converge and excited state solutions tend to "collapse" back to the ground state. When using state-averaging in TeraChem, you must specify how many states are included in the average (these will be selected as the lowest states) and what the weights are for each of the states. Often these weights are chosen to be equal, in order to ensure the correct treatment of conical intersections.

Dynamically-weighted state averaging

Choosing equal weights for all states in a state-averaged CASSCF calculation can be problematic. If one of the states is very high in energy and well-separated from the other electronic states, incuding it in the average might bias the orbitals to describe that stae (which is often not of any interest). One could instead limit the average to one fewer state, but one generally wants the active space and state-averaging to be consistent across the potential energy surface (e.g. when doing dynamics or computing a reaction path). It is often the case that N states are needed at some geometry (e.g. the Franck-Condon region for butadiene where there are two near-degenerate low-lying excited states) but after relaxation on the excited state, only N-1 of these states are relevant (e.g., after torsion of butadiene, one the near-degenerate states is strongly stabilized and intersects with the gorund state). The dynamic weighting procedure makes the weights depend on the energies of the individual states. This gradually removes high-lying states from the state-average and can provide a better global description than the usual insistence on a fixed number of states in an equally-weighted average. It also ensures that states which are near-degenerate (e.g. near a conical intersection) are assigned the same weight, as required.

Dynamic weighting also cures a specific failure of equally-weighted SA-CASSCF. The SA-CASSCF energy is invariant to rotations among the averaged states only when those states carry equal weight; consequently, when a state that is not included in the average (zero weight) crosses one that is included, the state-averaged energy — and its gradient — become discontinuous. Simply adding more states to the average rarely helps, because molecular systems have an ever-increasing density of states at higher energy. By making the weights smooth functions of the state energies, dynamic weighting removes these discontinuities and yields energy-conserving excited-state dynamics in cases where SA-CASSCF fails.

TeraChem implements the dynamical-weighting-with-spline (DWS) scheme of Glover. The lowest "states of interest" are given equal weight, and each higher state's weight is tapered smoothly to zero with a cubic polynomial spline,

\[ w_i \propto \begin{cases} 1 & \Delta E_i \le 0 \\[2pt] 1 - 3\left(\Delta E_i/\alpha\right)^2 + 2\left(\Delta E_i/\alpha\right)^3 & 0 < \Delta E_i < \alpha \\[2pt] 0 & \Delta E_i \ge \alpha \end{cases} \]

where \(\Delta E_i\) is the state's energy measured from the highest state of interest and \(\alpha\) is the taper width (the dwbandwidth keyword, default ≈ 3 eV). The weights are recomputed self-consistently within the CASSCF optimization, and TeraChem provides the corresponding analytic gradients so that the scheme can be used for geometry optimization and ab initio molecular dynamics. Choosing \(\alpha\) too small makes the weights vary too rapidly for the CASSCF to converge; too large pulls in many high-lying states and the surface drifts away from the equally-weighted SA result.

Files produced by a CASSCF calculation

Molden files - There are a number of these, with different orbitals
Singlet.x.molden: The natural orbitals corresponding to the xth singlet state
state_averaged_natural_orbital.molden:

Effective number of unpaired electrons (ENUE)

The effective number of unpaired electrons (ENUE) is a measure of the polyradical character of a wavefunction — it counts how many electrons are effectively unpaired in a given electronic state. For a closed-shell single-determinant state ENUE is zero; for a perfect diradical it approaches two, and so on. Because CASSCF states are genuinely multiconfigurational, ENUE is a convenient diagnostic for how multireference a given state is (e.g. along a bond-dissociation coordinate or near a conical intersection).

TeraChem computes ENUE from the natural-orbital occupation numbers \(n_i\) of the state's one-particle reduced density matrix using the quadratic index of Staroverov and Davidson (equivalent to the "odd electron" distribution of Takatsuka, Fueno, and Yamaguchi),

\[ N_u = \sum_i n_i\,(2 - n_i), \]

where the \(n_i \in [0, 2]\) are the eigenvalues of the orthogonalized density. Each doubly-occupied (\(n_i = 2\)) or empty (\(n_i = 0\)) orbital contributes nothing, while a singly-occupied orbital (\(n_i = 1\)) contributes exactly one. The density is orthogonalized with Löwdin (\(S^{1/2}\)) population analysis before the trace is taken, so \(N_u = \mathrm{tr}\!\left(2P - P^2\right)\) with \(P = S^{1/2}\,D\,S^{1/2}\).

ENUE is enabled with cas_enue yes (default no) and is available for the CAS family of methods only: CASSCF, CASCI, and FOMO-CASCI. It is reported per electronic state and for every spin multiplicity in the calculation. The value is printed to the main output log as

<Multiplicity> state <n> ENUE is <value>

together with a per-atom Mulliken charge breakdown (this requires the default Mulliken charge scheme). The job also writes each state's AO-basis one-particle density matrix to a separate file opdm_AO_<Mult>_<n>.bin (e.g. opdm_AO_Singlet_1.bin) in the scratch directory (scrdir, default ./scr); see Binary Files for the layout.

Worked example: ethylene π/π* diradical character

A CAS(2,2) calculation on ethylene — two electrons in the π and π* orbitals, state-averaged with equal weights over the lowest two singlets — is the textbook illustration of ENUE. The input enables it with a single keyword on top of an otherwise ordinary state-averaged CASSCF job:

coordinates   ethylene.xyz
basis         6-31g*
method        hf
run           energy
charge        0
spinmult      1

casscf        yes
closed        7
active        2
cassinglets   2
casweights    [1.0,1.0]
castargetmult 1
castarget     1

cas_enue      yes
end

For planar ethylene this prints

Performing effective number of unpaired electrons analysis (ENUE).
Computing ENUE using Lowdin population analysis

Singlet state 1 ENUE is   0.29558367015999

Singlet state 2 ENUE is   2.00000000000001

Singlet state Mulliken charges for ENUE:

  Atom   Root  1     Root  2
------------------------
 C      0.1478      1.0000
 C      0.1478      1.0000
 H      0.0000      0.0000
 ...

The ground state (N) is predominantly closed-shell, with only a small residual \(N_u \approx 0.30\), while the excited ππ* state (V) has two effectively unpaired electrons, \(N_u \approx 2.0\). If one methylene is rotated 90° about the C=C axis, the π and π* orbitals become degenerate and the ordering flips: the ground state becomes a perfect diradical with \(N_u \approx 2.0\) — the signature of twisted ethylene — while the second state collapses to \(N_u \approx 0\).

Geometry	State	Character	\(N_u\)
Planar	S₀ (N)	closed-shell ground state	0.30
Planar	S₁ (V)	open-shell ππ*	2.00
90° twisted	S₀	perfect diradical	2.00
90° twisted	S₁	ionic / closed-shell	0.00

These two cases ship as the integration tests tests/tests/CASSCF/cas_ethylene_enue_planar and cas_ethylene_enue_twist, which make good starting points for your own ENUE calculations.

Parameters for CASSCF jobs

A CASSCF calculation is enabled by casscf yes. The active space is defined by closed (number of doubly-occupied orbitals) and active (number of orbitals in the active space). The number of CSF/determinant configurations to compute is set by the per-multiplicity keywords cassinglets, castriplets, etc. For example, closed 7 / active 2 / cassinglets 2 defines a CAS(2,2) with two singlet states.

Active space and states

Keyword	Type	Default	Description
`casscf`	bool	no	Master switch: enable CASSCF orbital optimization
`casci`	bool	no	Enable CASCI (no orbital optimization); often paired with `fon yes` to give FOMO-CASCI
`closed`	int	—	Number of doubly-occupied orbitals (orbitals below the active space). Required for any CAS method.
`active`	int	—	Number of active orbitals. Required for any CAS method.
`activeorbs`	string	not set	Comma-separated list specifying which orbitals are placed in the active space (overrides the default contiguous choice)
`cassinglets`	int	0	Number of singlet states to include in the CI / state average
`casdoublets`	int	0	Number of doublet states
`castriplets`	int	0	Number of triplet states
`casquartets`, `casquintets`, `cassextets`, `casseptets`	int	0	Higher multiplicities
`castarget`	int or `sa`	0	Target state for gradient calculations (0-based, within the chosen multiplicity). Set to `sa` to compute a state-averaged-energy gradient.
`castargetmult`	int	1	Spin multiplicity of the target state

State averaging

Keyword	Type	Default	Description
`casweights`	string	equal weights	Comma-separated list of weights for each state in the state-averaged energy
`dynamicweights`	bool	no	Use energy-dependent dynamic weighting via the cubic-spline (DWS) scheme (see "Dynamically-weighted state averaging" above)
`dwbandwidth`	float	≈ 0.110 a.u. (3 eV)	Energy range \(\alpha\) (a.u.) over which the spline tapers a state's weight from 1 to 0, measured from the highest state of interest

Convergence

Keyword	Type	Default	Description
`casscfmaxiter`	int	as set internally	Maximum number of CASSCF macroiterations
`casscfconvthre`	float	as set internally	Gradient-norm convergence threshold
`casscfenergyconvthre`	float	as set internally	Energy-change convergence threshold
`casscforbnriter`	int	as set internally	Maximum number of orbital Newton-Raphson microiterations
`casscfnriter`	int	as set internally	Maximum number of Newton-Raphson iterations in the orbital step
`casscftrustmaxiter`	int	as set internally	Maximum number of trust-region subproblem iterations
`casscftrustconvthre`	float	as set internally	Convergence threshold for the trust-region subproblem
`casscflambdainit`	float	as set internally	Initial trust-region Lagrange multiplier
`casscflambdamin` / `casscflambdamax`	float	—	Bounds on the trust-region multiplier

CP-SA-CASSCF (response equations for gradients)

Keyword	Type	Default	Description
`cpsacasscfmaxiter`	int	as set internally	Maximum iterations in the coupled-perturbed state-averaged CASSCF equations
`cpsacasscfconvthre`	float	as set internally	Convergence threshold for the CP-SA-CASSCF response equations
`cpsacasscfsolver`	string	as set internally	Solver: `direct`, `pcg`, `inc_diis`, or `inc_block_diis`

Initial guesses and orbital control

Keyword	Type	Default	Description
`casguess`	filename	not set	Use the molecular orbitals in this file as an initial guess for the CASSCF orbitals
`casaopguess`	filename	not set	Rotate the starting orbitals to align with the active space defined in this file
`zvecguess`	filename	not set	Initial guess for the Z-vector used in gradient evaluation

CI solver

Keyword	Type	Default	Description
`ci_solver`	string	`hard_coded`	CI solver: `hard_coded`, `direct`, `rank_reduced`, or `dmrg`

Fractional orbital occupation (FOMO)

When combined with casci yes, FOMO produces FOMO-CASCI — see CASCI and FOMO-CASCI for what FOMO-CASCI does and why. With casscf yes, it gives the related FOMO-CASSCF method. The fractional-occupation procedure itself is part of the SCF; see FOMO-SCF for the distributions, temperature, and the full set of fon* keywords.

Keyword	Type	Default	Description
`fon`	bool	no	Enable fractional orbital occupations
`fon_temperature`	float	—	Electronic temperature (a.u.) controlling the spread of fractional occupations
`fon_method`	string	`gaussian`	Distribution: `gaussian`, `fermi`, `marzari-vanderbilt` (alias `mv`), or `constant`

Output and analysis

Keyword	Type	Default	Description
`casgradprint`	bool	no	Print the CAS gradient in the output
`cascharges`	bool	no	Compute and print state-specific atomic charges
`cas_enue`	bool	no	Compute the effective number of unpaired electrons (ENUE) per state (see "Effective number of unpaired electrons" above)
`cas_ntos`	bool	no	Compute natural transition orbitals between states
`casrelaxeddipoles`	bool	no	Compute orbital-relaxed dipole moments
`cassoc`	bool	no	Compute spin-orbit couplings between CAS states

References

The GPU-accelerated, AO-based CASSCF method, its analytic gradients, and nonadiabatic coupling vectors are described on the Theory and implementation page.
Dynamically-weighted state averaging (the cubic-spline DWS scheme and its analytic gradients): W. J. Glover, "Communication: Smoothing out excited-state dynamics: Analytical gradients for dynamically weighted complete active space self-consistent field," J. Chem. Phys. 141, 171102 (2014). doi:10.1063/1.4901328
Effective number of unpaired electrons (the ENUE index used by cas_enue): V. N. Staroverov and E. R. Davidson, "Distribution of effectively unpaired electrons," Chem. Phys. Lett. 330, 161–168 (2000). doi:10.1016/S0009-2614(00)01088-5 The same index was introduced earlier as the "odd electron" distribution by K. Takatsuka, T. Fueno, and K. Yamaguchi, Theor. Chim. Acta 48, 175 (1978).