How to compute SOCs and how to make sense out of them

Q-Chem can deal with spin-orbit effects by using a perturbative approach in which one first computes non-relativistic zero-order states (by CIS, EOM-CCSD, ADC, etc) and then computes matrix elements of the spin-orbit part of the Breit-Pauli Hamiltonian. This matrix (expressed in the basis of zero-order states) can then be treated as a perturbation that mixes non-relativistic states, giving rise to the spin-orbit perturbed states. Because there are different experimental manifestations of spin-orbit effects (i.e., spectroscopic splittings, intensity borrowing, inter-system crossing, magnetic anisotropies, to name a few), there are different ways in which spin-orbit couplings (SOCs) can be related to the experimental observables. Hence, the user should know what is it they want to compute and carry out the post-processing of the raw data (zero-order non-relativistic states, their energies and properties, and SOC matrix) accordingly.

Here are some pointers to help to navigate through the relativistic calculations. The key point to remember is that “coupling” and “splitting” are not the same. These two terms have distinct meanings. Couplings are the matrix elements of spin-orbit operator between non-relativistic states: <SM| Hso | S’M’>. This is precisely what is computed and printed in Q-Chem’s output. The definition of the splitting depends on a particular type of a transition; hence, the calculations of splittings are system-dependent and should be carried as a post-processing step.

Since spin-orbit operator depends on spin, the values of the matrix elements depend on the spin projection. The values of the individual matrix elements also depend on the orientation of coordinate system — this is perfectly normal.

In order to connect to the experimental observables, such as SOCC (spin-orbit coupling constant), one needs to have the SOC matrix elements for the entire multiplet. For example, when computing couplings between the singlet and triplet states, one needs to to obtain the couplings for all three components of the triplet (Ms=-1, 0, 1). These components are degenerate in non-relativistic calculations and usually only one component is computed (Ms=0).

In the EOM-CC implementation of SOC, the algorithm computes one spin-projection in a given calculation. Then it generates matrix elements for all other spin projections by using Wigner-Eckart’s theorem. The details of the theory are described here:
http://iopenshell.usc.edu/pubs/pdf/jcp-151-034106.pdf

SOCC stands for spin-orbit coupling constant, which is defined as the sum of the squared matrix elements between the states of all possible spin projections in a given transition A->B:
SOCC = \sqrt(\sum_{M, M’} |<B, S, M | H_so | A, S, M>|^2).
In contrast to SOC, SOCC is independent of molecular orientation. It is related to ISC rates via Fermi’s golden rule.

Breit-Pauli spin-orbit operator has one- and two-electron parts. The full two electron part is expensive, but it can be very accurately approximated by the spin-orbit mean-field (SOMF) treatment. This is what is implemented for EOM-EE. The output contains both one-electron part of SOCs and the full SOCs with two-electron part evaluated with SOMF approximation.

What does it mean when an EOM-CC spin-orbit coupling computation returns this for a singlet-singlet transition? Does this message simply indicate that a singlet-singlet spin-orbit coupling is zero by symmetry?

Clebsh-Gordan coefficient: <0.000,0.000;1.000,0.000|0.000,0.000> = 0.000
WARNING! Clebsh-Gordon coefficient is too small to find a rotationally invariant part of triplet density.
Please use another EOM approach with non-zero Clebsh-Gordon coefficients
I am going to skip this transition

The implemented algorithm can be presented as follows:

1. S^2, spin projections, and one-particle transition density matrix is computed between the pair of EOM states.
2. The Clebsh-Gordan coefficients are computed based on the spin projections and multiplicities. The multiplicities are recognized based on the S^2.
3. The triplet transition density matrix is extracted from the one-particle transition density matrix.
4. The spinless density matrix is obtained from the triplet density matrix and the Clebsh-Gordan coefficient, formed by the multiplicities and spin projections of the starting EOM states. This step involves division by the Clebsh-Gordan coefficient (see the paper for details). Therefore, this coefficient should not be zero. If it is zero, the algorithm checks for that and skips the transition.
5. Reduced matrix elements are formed as a contraction between integrals and the spinless density matrix.
6. All matrix elements are generated from the reduced matrix elements based on the corresponding Clebsh-Gordan coefficients.

For example, the triplet-triplet couplings between the states with zero Sz is always zero because the corresponding Clebsh-Gordan coefficient is zero. The algorithm cannot start from these states to generate all SO matrix elements between the triplets with all possible spin projections: -1, 0, +1 in this case. The warning is printed to alarm the user that this transition cannot be treated with the approach.
In order to compute the SOC matrix in the triplet-triplet case, the user needs to generate triplet with non-zero spin projections. For example, one can evaluate them between SF states from the closed-shell reference or as EE states, obtained from the triplet CCSD reference. The particular case of singlet-singlet couplings always gives

Singlet-singlet couplings are always zero, because the corresponding Clebsh-Gordan coefficient is zero. This transition is skipped by the algorithm based on that.

I agree that certain individual cases, where all SO matrix elements are zero, such as singlet-singlet couplings or singlet-quintet couplings, can be hard-coded to produce zero matrix in the output. At the moment this is not implemented and the transition is skipped based on the zero value of the initial Clebsh-Gordan coefficient.