Reviews in Computational Chemistry

164 Spin–Orbit Coupling in Molecules

energies Eðk1;aÞ are obtained by diagonalizing the interaction matrix in a basis spanned by the d components j ðk0;bÞi of j ðk0Þi, where b is an index from 1 to d.

0 h

ð0Þ

^

ð0Þ

h

ð0Þ

^

ð0Þ

h

ð0Þ

^

ð0Þ

1

ð0;Þj

^ j

ð0;

Þi

ð0;

Þj

^ j

ð0;Þi

. . .

ð0;Þ j

^ j

ð0;Þi

k 1

HSO

k 1

k 2

HSO

k 1

k d

HSO

k 1

C

B h k;1jH.SOj k;2i

h k;2jH.SOj k;2i

.. . . h k;djH.SOj k;2i

190

B

.

.

.

.

.

C

½ &

B

.

.

.

C

@

A

B

0

H^ SO

0

Þ

0

Þ

H^ SO

0

. . .

0

H^ SO

0

C

B

ð Þ

ð

h

ð

ð Þ

h

ð Þ

ð Þ

C

B h

k;1j

j

k;di

k;2j

j

k;di

k;dj

j

k;di

C

Only one of the matrix elements needs to be evaluated explicitly. All others can be obtained from the reduced matrix element by means of the Wigner– Eckart theorem, Eq. [166]. The eigenvectors

j kð0;aÞi ¼	d
	jckð0;aÞ;bj kð0;bÞi	½191&

X

b¼1

are zeroth-order wave functions adapted to this particular perturbation.

If the degeneracy is lifted completely in first-order or if at least secondorder effects do not introduce an additional splitting of degenerate levels, the

ð2Þ

second-order energy Ek;l can be expressed as

d0

ð0Þ ^

ð0Þ

ð0Þ ^

ð0Þ

ð2Þ

X X

1

h k;ajHSOj i;b0 ih i;b0 jHSOj k;ai

Ek;a

¼

i

k b0

¼

1

Eið0Þ

Eð0Þ

½192&

6¼

k

The corresponding first-order perturbed wave function reads

			d0			ð0Þ ^	ð0Þ
ð1Þ	X X							ð0Þ
	1					h i;b0 jHSOj k;ai				½193&
j k;ai ¼	i	k b0		¼	1	Eið0Þ	Eð0Þ	j i;b0	i
	¼6						k

If these conditions are not fulfilled, an auxiliary condition is necessary to deter-

ð0Þ

mine the proper linear combinations of the j k;ai: although each set of orthogonal eigenfunctions qualifies as a proper wave function as long as the degeneracy pertains, the wave function coefficients are required to change smoothly when the additional perturbation is switched on adiabatically. Symmetry properties of the involved states and operators may aid in predetermining these linear combinations. In general, this procedure is too complicated, however, and one resorts to quasi-degenerate perturbation theory.

Moreover, first and second orders of perturbation theory are not defined in a stringent way in a molecule. Consider, for instance, all potential energy curves of a diatomic molecule that correspond to a specific dissociation

Computational Aspects

165

channel. Near the equilibrium geometry, the low-lying electronic states are usually well separated. Their mutual interaction is then described properly by second-order perturbation theory. In the separated atom limit, first-order degenerate perturbation theory applies. In between, but still close to the dissociation limit, the states are nearly degenerate—or quasi-degenerate as we might say. Neither of the procedures is then strictly applicable and orders of perturbation theory are not well defined.

Quasi-Degenerate Perturbation Theory

An often chosen way out of this dilemma is to set up a so-called perturbation matrix in the basis of eigenvectors of a spin-free secular equation, multiplied by an appropriate spin function. In this basis, matrix elements of the spin-free Hamiltonian occur only in the diagonal, of course. Due to symmetry, diagonal spin–orbit matrix elements come along only with complex wave functions. In a Cartesian basis, the integrals of the electrostatic Hamiltonian

and the y component of ^ are real, whereas x and z integrals exhibit an

HSO

imaginary phase. (The spatial parts of the integrals are imaginary for all Cartesian components. It is the choice of an imaginary phase for S^y that makes the

matrix elements of the y component of ^ real.) As before, matrix elements

HSO

are computed only for one representative of a spin multiplet; all other matrix elements are generated by use of the WET.

Similar to first-order degenerate perturbation theory, perturbed energies and wave functions are obtained by matrix diagonalization. This approach has therefore been named quasi-degenerate perturbation theory. In spite of this designation, the procedure is also applicable to states with large energy separations. Other authors prefer to call it LS contracted spin–orbit configuration interaction (SOCI) in order to stress its relation to configuration interaction procedures in the limit of a complete set of eigenvectors.108 The perturbation matrix is in most cases small enough for the application of standard complex Hermitian eigenvalue and eigenvector solvers. Otherwise one can resort to the same methods as in SOCI.

One of the great advantages of quasi-degenerate and Rayleigh–Schro¨ dinger perturbation theories over SOCI procedures is that different levels of sophistication can easily be mixed. For reasons of efficiency or technical limitations, it is in general not possible to a perform a full CI calculation that yields the exact solution within a given basis set. Instead, configurations are selected according to some criterion such as excitation class or energy. Unfortunately, electron correlation contributions are slowly convergent. Furthermore, truncated CI expansions are not size extensive, that is, the correlation energy does not scale properly with the number of electrons. For spin-free states, long-standing experience exists on how to estimate correlation contributions from discarded configurations or excitation classes.109–115 These extrapolated energies or the eigenvalues of correspondingly dressed Hamiltonians can be taken as diagonal elements combined with spin–orbit matrix

166 Spin–Orbit Coupling in Molecules

elements of single excitation CI or smaller single and double excitation CI expansions.41,108,116

A weakness of these methods lies in the limited number of zeroth-order states that are used for an expansion of the first-order perturbed wave function. In particular, it has been demonstrated that probabilities of spin-forbid- den radiative transitions converge slowly with the length of the perturbation expansion.92

Variational Perturbation and Response Theory

As an alternative to sum-over-states methods, the perturbation equations can be solved directly. In the context of spin–orbit coupling, reviews on this

subject have recently been given by Yarkony	117	˚	118
		and by Agren et al.

Hess et al.119 utilized a Hamiltonian matrix approach to determine the spin–orbit coupling between a spin-free correlated wave function and the configuration state functions (CSFs) of the perturbing symmetries. Havriliak and Yarkony120 proposed to solve the matrix equation

Hð0Þ		Eð0Þ	1	Þ	ð1Þ	¼	H ð0Þ	½	194	&
ð		i			i		SO i

directly for ði1Þ in the basis of CSFs. The direct solution of such a perturbation

equation is generally known as Hylleraas (or variational) perturbation theory.121,122 Yarkony123 developed this approach further by introducing a sym-

bolic matrix element method, thereby extending the limits relating to the dimensionality of the Hamiltonian matrix.

Also in response theory the summation over excited states is effectively replaced by solving a system of linear equations. Spin–orbit matrix elements are obtained from linear response functions, whereas quadratic response functions can most elegantly be utilized to compute spin-forbidden radiative transition probabilities. We refrain from going into details here, because an

˚	118	While
excellent review on this subject has been published by Agren et al.

these authors focus on response theory and its application in the framework of CI and multiconfiguration self-consistent field (MCSCF) procedures, an analogous scheme using coupled-cluster electronic structure methods was presented lately by Christiansen et al.124

Variational Procedures

For compounds containing heavy atoms, spin–orbit and electron correlation energies are approximately of the same size, and one cannot expect these effects to be independent of each other. A variational approach that treats both interactions at the same level is then preferable to a perturbation expansion. Special care is advisable in the choice of the spin–orbit operator in this case. The variational determination of spin–orbit coupling requires a spin–orbit

Computational Aspects

167

operator that is bounded from below such as the no-pair Hamiltonian,103 the

zeroth-order regular approximation Hamiltonian (ZORA),125–127 or the ECP spin–orbit operators.39,40,52 Otherwise one has to constrain the basis set (e.g.,

by an energy selection) to prevent a variational collapse.

jj Coupling Methods

Two-component approaches employing one-particle spinors as basis functions are referred to as jj coupling methods. They are typical one-step procedures that take electron correlation and spin–orbit coupling into account simultaneously. Kramers-restricted two-component calculations at the Har-

tree–Fock, second-order Møller-Plesset, singles and doubles CI, and coupled cluster levels have been reported for atoms and diatomic molecules.128,129

Technically Lee and co-workers utilized modified versions of four-component programs.130 Two-component wave function approaches that include SOC right away at the orbital optimization stage are rarely applied to molecules with more than two atoms. Possible reasons are the complications that arise from a formulation of subsequent correlation treatments due to complex-valued orbital coefficients. Among the most widely used twocomponent methods in the field of chemistry are spin density functional meth- ods.131–133 In the Kohn–Sham approximation, complex orbital coefficients emerge, but this is not problematic due to the restriction to a single determinant.

Intermediate Coupling Procedures

Most of the variational treatments of spin–orbit interaction utilize onecomponent MOs as the one-particle basis. The SOC is then introduced at the CI level. A so-called SOCI can be realized either as a oneor two-step procedure. Evidently, one-step methods determine spin–orbit coupling and electron correlation simultaneously. In two-step procedures, typically different matrix representations of the electrostatic and magnetic Hamiltonians are chosen.

In a one-step SOCI, a Hamiltonian matrix is set up in the basis of CSFs or determinants. It contains simultaneously matrix elements of the usual electrostatic Hamiltonian and of the spin–orbit operator and is complex in general. Since spin–orbit coupling mixes spin and spatial degrees of freedom, there are fewer possibilities of blocking the Hamiltonian matrix than for a spin-inde-

pendent operator. To this end, one makes use of double group symmetries and transforms the Cartesian basis functions to so-called Kramers pairs.52,134

By a special choice of phases, the matrix can even be made real in some cases. In general, however, very large complex eigenvalue problems have to be solved in a SOCI.

Due to its dimensionality, a complete diagonalization of the CI matrix is not feasible. One usually looks only for a few roots by means of an iterative procedure. For quantum chemical purposes, the Davidson algorithm or modifications thereof have proven to be well suited.135–137 In its original

168 Spin–Orbit Coupling in Molecules

formulation, the Davidson method is applicable to real symmetric matrices only. With some care, it can nevertheless be applied also to a complex eigenvalue problem. Every complex Hermitian matrix C can be written as

C ¼ A þ iB

½195&

where A and B are real matrices. Since C is Hermitian, A is a symmetric

matrix, and B is skew symmetric.	g				~
		Similarly, the eigenvectors Zi of C can be
~				~	~	itself is real.
divided into a real part Ui and an imaginary part iVi, where Vi
The complex eigenvalue problem
	~		~			½196&
CZi ¼ Ei Zi

can then be replaced by a real symmetric eigenvalue problem of two times its original dimension.138

A	B		~	! ¼ Ei	~	!
			Ui		Ui
B	A		V~i		V~i		½197&

All eigenvalues Ei of the real eigenvalue problem [197] are doubly degenerate.

This degeneracy is of purely technical origin and should not be confused with Kramers degeneracy,139,140 which may occur in addition. (For instance, four

degenerate roots are obtained for a doublet state, i.e., two for each Kramers

~ ~ t

level.) In addition to the transposed eigenvector ðUiViÞ , a second one with

~ ~ t

structure ð ViUiÞ is obtained. One of these solutions can be discarded. Alter-

natively, root-homing procedures can be employed to avoid its evaluation right from the beginning.136,137

Older versions of SOCI programs are very I/O intensive because they

used to store the Hamiltonian matrix on disk and read it in every iteration step.52,141 Integral-driven direct methods for spin–orbit coupling came up in the mid 1980s123,142 following the original fomulation of direct CI methods

for spin-independent Hamiltonians.143–146 Modern direct SOCI programs can easily handle several million determinants.108,147–151

The convergence of the iterative determination of eigenvalues and eigenvectors is accelerated appreciably if spin–orbit CI and quasi-degenerate perturbation theory procedures are combined. To this end, the perturbation matrix is set up in the basis of the most important LS contracted CI vectors j ðm0;1Þi. The solutions of this small eigenvalue problem

j kð1;1Þi ¼	sf CI states
	X	ckmð1;1Þj mð0;1Þi	½198&

m

are then used to start up the Davidson iteration.

gA matrix B is said to be skew-symmetric if Bij ¼ Bji.

Computational Aspects

169

As the number of electrons increases, the dimension of the CI space becomes increasingly large, even if excitations are restricted to single and double replacements. Many approximate schemes to restrict the number of configurations in a SOCI have been devised. The selected intermediate coupling CI (SICCI) code is an extension of the original COLUMBUS DGCI program and

utilizes electrostatic and magnetic interactions as a selection criterion.152 The SPDIAG program141 within the BNSOC package153 is based on the MRD-CI

approach. The latter makes use of a correlation energy criterion for configuration selection and estimates the contribution of the discarded configurations to the spin-free correlated energy Eðm0;1Þ by means of Epstein–Nesbet perturbation theory.111 If one assumes that the expansion [198] represents a decent approx-

imation to the SOCI solution, the MRD-CI extrapolation scheme can easily be extended to the spin–orbit coupled case.154,155 Also other approaches toward

a balanced treatment of spin–orbit interaction and electron correlation are based on a manipulation of the spin-free energies and wave functions. A pure shift of excitation energies in the Davidson135 start-up matrix is not suf-

^

ficient, because the original eigenvalues of Hð0Þ are restored during the iterative process.141 In the so-called spin-free state shifted (SFSS) SOCI method, this problem is circumvented by introducing a projector on the set of ðm0;1Þ.156

sf CIXstates

^ sfss

^

Emj

ð0;1Þ

ih

ð0;1Þ

j

½

199

&

HSO

¼ HSO þ

m

m

m

with

Em ¼ ð

E

m

E

XÞ ð

Eð0;1Þ

Eð0;1Þ

Þ

½

200

&

m

X

The Em quantities shift the spin-free excitation energies ðEðm0;1Þ EðX0;1ÞÞ, calculated at a lower level of correlation treatment, to the exact values ðEm EXÞ or at least to higher accuracy estimates. Herein, X denotes a common reference state, in general the electronic ground state.

Due to the structure of the spin–orbit Hamiltonian, off-diagonal blocks of the CI matrix are dominated by single excitations. On the other hand, for electron correlation effects (diagonal blocks) double and higher excitations are decisive. In many approximate schemes, this fact is exploited. One of the first of these approaches was the relativistic CI (RCI) algorithm by Balasubramanian.157 This is a two-step procedure. In a first step, large-scale multiconfiguration SCF and multireference singles and doubles CI (MRSDCI) calculations are carried out for a spin-free Hamiltonian to determine a set of natural orbitals (NOs). In the second step, spin–orbit interaction integrals are transformed to the NO basis and added to the one-electron Hamiltonian matrix elements of a smaller, truncated MRSDCI in relativistic symmetries. DiLabio and Christiansen158 studied the separability of spin–orbit and correlation energies for

170 Spin–Orbit Coupling in Molecules

sixth-row main group compounds. They propose to determine the energy shift DSOs due to spin–orbit coupling as the difference between spin-free and intermediate coupling single excitation CI calculations. The energy of the correlated spin–orbit coupled states is then estimated by summing DSOs into the spin-free correlated energies Elssd . Rakowitz et al.159 employed a SFSS SOCI approach on the Irþ ion. They demonstrated that the heavily spin–orbit perturbed spectrum of this ion can be obtained in good accord with experiment at the single excitation level, if higher level correlated electrostatic energies are used to determine the shifts Em. Recently Vallet et al.160 extended the effective Hamiltonian approach, originally implemented in the CIPSO code at the level of quasi-degenerate perturbation theory,41 to be used in intermediate coupling

eff	^
calculations. Instead of representing Hel	þ HSO in the basis of LS contracted

states, the new effective and polarized spin–orbit CI (EPCISO) operates in a determinantal model space. This model space is chosen as the union of the separate reference spaces in nonrelativistic symmetries augmented by the most significant configurations contributing to either the correlation energy or spin–orbit coupling.

COMPARISON OF FINE-STRUCTURE SPLITTINGS WITH EXPERIMENT

Spectroscopic parameters of a molecule are derived from experimentally determined spectra by fitting term values to a properly chosen model Hamiltonian.161 Usually, the model Hamiltonian is an effective one-state Hamiltonian that incorporates the interactions with other electronic states parametrically. In rare cases, experimentalists have used a multistate ansatz like the supermultiplet approach162 to fit the rovibronic spectra of strongly interacting near-degenerate electronic states. The safest way of comparing theoretical data to experiment is to compute the spectrum and to fit the calculated term energies to the same model Hamiltonian as the experimentalists use.

The vibrational dependence of an effective molecular parameter A is usually expressed as

	X	ak vi þ	d		k
Avi ¼ Ae þ			i		½201&
	k		2

where Ae is the value of the property at the equilibrium distance, and di is the degeneracy of the vibrational mode vi. The vibrational dependence is caused by the anharmonicity of the potential energy surface and by the variation of A with the vibrational coordinate Q. Such information can in practice be obtained from a quantum chemical treatment only for diatomic and triatomic molecules. To this end, A is fitted to a function of Q—mostly a polynomial or

Comparison of Fine-Structure Splittings with Experiment

171

a cubic spline—and vibrationally averaged, that is, AðQÞ is weighted with the probability density of a particular vibrational state wvi ðQÞ and integrated over the vibrational coordinate Q

Avi ¼ hwvi ðQÞjAðQÞjwvi ðQÞi

½202&

Alternatively, expectation values computed at the equilibrium geometry of a given electronic state can be directly compared with experimental parameters extrapolated from rovibrational branches.

In many instances, the ‘‘fine-structure’’ splitting caused by (first-order) spin–orbit coupling is considerably larger than the energy separation between adjacent rotational levels—at least for low values of the total angular momentum J—and rotational excitations can thus be neglected in a first approximation. Hydrides are exceptional in this respect because of their low reduced masses and the resulting large rotational constants B. Methods for computing rovibrational spectra of diatomic molecules from ab initio data including spin– orbit and rotational coupling were proposed among others by Yarkony117 and by the author163 employing Hund’s case ðbÞ and ðaÞ basis functions, respectively.161 For a review on the theoretical determination of rovibronic spectra of triatomic molecules, please refer to the work of Peric´ et al.164 A general expression for the rotational dependence of a spectroscopic parameter cannot be given. Its functional form varies with the type of basis functions chosen for describing the actual rotation of the nuclear frame; the choice of basis functions in turn depends on the order in which the angular momenta are coupled (Hund’s cases in linear molecules) and the type of rotor [spherical top, symmetric top (prolate, oblate), asymmetric top]. Looking at these particular cases in detail goes far beyond the scope of the present chapter. For additional reading, see, for example, Refs. 165–167.

First-Order Spin–Orbit Splitting

Only spatially degenerate states exhibit a first-order zero-field splitting. This condition restricts the phenomenon to atoms, diatomics, and highly symmetric polyatomic molecules. For a comparison with experiment, computed matrix elements of one or the other microscopic spin–orbit Hamiltonian have to be equated with those of a phenomenological operator. One has to be aware of the fact, however, that experimentally determined parameters are effective ones and may contain second-order contributions. Second-order SOC may be large, particularly in heavy element compounds. As discussed in the next section, it is not always distinguishable from first-order effects.

A phenomenological spin–orbit Hamiltonian, formulated in terms of tensor operators, was presented already in the subsection on tensor operators. Few experimentalists utilize an effective Hamiltonian of this form (see Eq. [159]). Instead, shift operators are used to represent space and spin angular

172	Spin–Orbit Coupling in Molecules
					Ω = 3/2
			ASO
Figure 14 Spin–orbit splitting pattern of a					Ω = 1/2
(regular) 2 r state. (The dashed line marks
the position of the unperturbed state.)			ASO > 0

momentum operators in their tensorial forms. Employing shift operator conventions and making use of the Russell–Saunders (LS) coupling scheme, the phenomenological spin–orbit Hamiltonian reads

~	~	^ ^		1		^	^ ^
^	^				^
ASOL S ¼ akL0S0			þ	2	a?ðLþS þ L SþÞ			½203&

For the sign of the spin–orbit parameter ASO, the following conventions apply:

*A positive ASO denotes a regular state (Figure 14). In regular states the sublevel with smallest J or quantum number is lowest in energy.

*A negative ASO denotes an inverted state (Figure 15). In inverted states the sublevel with largest J or quantum number is lowest in energy.

Typically, states with less than half-filled shells (particle states) are regular, whereas states with more than half-filled shells (hole states) are inverted.

Atoms: The Lande´ Interval Rule

To determine the first-order splitting pattern of an atomic state in terms of the phenomenological spin–orbit parameter ASO [Eq. 203], we utilize

				Ω = 1
		2 ASO
				Ω = 2
		2 ASO
	Figure 15 Spin–orbit splitting pattern of an			Ω = 3
		ASO < 0
	(inverted) 3 i state.

Comparison of Fine-Structure Splittings with Experiment

173

Russell–Saunders coupling, that

is,

~

~

~

~2

~

~

2

^

^

^

^

^

^

¼

~ 2

^

~

J

¼ ðL þ S Þ. From J

¼ ðL þ S Þ

^

^ ^

L

þ 2L S þ S, we obtain

1

~

~

2

2

2

^

^

^

^

^

½204&

L S ¼

2

ðJ

L

S Þ

Exploiting the fact that in first-order perturbation theory all sublevels exhibit

identical ^ 2 and S^2 eigenvalues, respectively, yields a spin–orbit splitting of

L

EðJÞ EðJ 1Þ ¼ ASOJ

½205&

for two neighboring atomic fine-structure levels with quantum numbers J and J 1. This expression is the famous Lande´ interval rule.

As an example where this rule can be applied favorably to determine an atomic spin–orbit parameter, consider the first excited state of atomic copper, 2Dgðd9s2Þ. The 2Dg state of copper is well separated from other electronic states that are allowed by symmetry to couple in second order. A 2D state gives rise to two fine-structure levels with total angular momentum quantum numbers J ¼ 32 and J ¼ 52. The d9s2 configuration corresponds to a hole state, and we therefore expect an inverted splitting pattern (Figure 16). From Lande´’s interval rule, we determine a splitting of 52 ASO. ASO may thus be calculated directly from the experimentally determined splitting168 of 2042 cm 1 by multiplying with 25, yielding ASO ¼ 816:8 cm 1; this value is in good agreement

with the theoretically determined spin–orbit parameter of ASO ¼ 802:4 cm 1.169

Things are not as easy for the 3Dgðd9s1Þ ground state of nickel. From a 3D atomic state, three fine-structure levels originate with total angular momentum quantum numbers J ¼ 1, J ¼ 2, and J ¼ 3 (Figure 17). The ratio of the measured fine-structure splittings ½Eð2Þ Eð1Þ&=½Eð3Þ Eð2Þ& amounts to 1.234 instead of 23 as expected from the Lande´ interval rule. The main cause for this deviation is an energy shift of the J ¼ 2 level due to the interaction with neighboring states of the same angular momentum quantum number. The J ¼ 3 and J ¼ 1 levels of 3Dg are only weakly perturbed and their energy

separation, EðJ ¼ 3Þ EðJ ¼ 1Þ ¼ 5ASO, may be taken to extract an experimental value for ASO ¼ 301:6 cm 1. The experimental spin–orbit

constant is in excellent agreement with the theoretical value of 301 cm 1.

J = 3/2

5/2 ASO

Figure 16 Spin–orbit splitting pattern of an

J = 5/2 inverted 2D atomic state.