Spin-orbit-coupling
Introduction
Spin-orbit matrix elements and eigenstates can be computed using either the Breit-Pauli (BP) operator or spin-orbit pseudopotentials (ECPs). The state-interacting method is employed, which means that the spin-orbit eigenstates are obtained by diagonalizing $\hat H_{el} + \hat H_{SO}$ in a basis of eigenfunctions of $\hat H_{el}$. The full Breit-Pauli SO-operator can be used only for MCSCF wavefunctions. For MRCI wavefunctions, the full BP operator is used for computing the matrix elements between internal configurations (no electrons in external orbitals), while for contributions of external configurations a mean-field one-electron fock operator is employed. The error caused by this approximation is usually smaller than 1 cm$^{-1}$.
The program allows either the computation of individual spin-orbit matrix elements for a given pair of states, or the automatic setting-up and diagonalization of the whole matrix for a given set of electronic states. In the latter case, matrix elements over one-electron operators are also computed and transformed to the spin-orbit eigenstates (by default, the dipole matrix elements are computed; other operators can be specified on the GEXPEC
or EXPEC
cards, see section One-electron operators and expectation values (GEXPEC)). Since it may be often sufficient to compute the spin-orbit matrix elements in a smaller basis than the energies, it is possible to replace the energy eigenvalues by precomputed values, which are passed to the spin-orbit program by the MOLPRO
variable HLSDIAG
.
Calculation of SO integrals
The one-and two-electron spin-orbit integrals over the BP Hamiltonian can be precomputed and stored on disk using the command
LSINT
[,X
] [,Y
] [,Z
] [,ONECENTER
] [;TWOINT
,twoint;] [;PREFAC
,prefac;]
X
, Y
, and Z
specify the components to be computed. If none of these is given, all three are evaluated. The advantage of precomputing the integrals is that they can then be used in any number of subsequent SO calculations, but this may require a large amount of disk space (note that there are 6 times as many integrals as in an energy calculation). If the LSINT
card is not given, the integrals are computed whenever needed. The keyword ONECENTER
activates the one-center approximation for one- and two- electron spin-orbit integrals. This can reduce drastically the computing time for large molecules. TWOINT
and PREFAC
can be used to control the accuracy of spin-orbit integrals. These thresholds are similar to TWOINT
and PREFAC
for standard integrals. The default value for PREFAC
is TWOINT/100
, and the default value for TWOINT
is $10^{-7}$. In the case when no integrals are precomputed, these thresholds can be specified as options for HLSMAT
or TRANLS
cards, see below.
The input for spin-orbit ECPs is described in section effective core potentials. Of course, in ECP-LS calculations the LSINT
card is not needed.
Calculation of individual SO matrix elements
Individual spin-orbit matrix elements can be computed within the MRCI program using
TRANLS
,record1.file, record2.file, bra2ms, ket2ms, lsop;
where
- record1.file Record holding the bra-wavefunction.
- record2.file Record holding the ket-wavefunction. Both records must have been generated using the
SAVE
directive of the MRCI program. - bra2ms $2 \times M_S$ value of the bra-wavefunction.
- ket2ms $2 \times M_S$ value of the ket-wavefunction.
- lsop Cartesian component of the Spin-orbit Hamiltonian.
This can be one of ${\tt LSX}$, ${\tt LSY}$, or ${\tt LSZ}$ in all electron calculations, and ${\tt ECPLSX}$, ${\tt ECPLSY}$, or ${\tt ECPLSZ}$ in ECP calculations. Starting from the MOLPRO
version 2008.1, more types are available which control the approximation level. These are described in section approximations used in calculating spin-orbit integrals and matrix elements.
Since the spin-orbit program is part of the MRCI program, the TRANLS
card must be preceded by a [MR]CI
card. If in the MRCI step several states of the same symmetry are computed simultaneously using the STATE
directive, the matrix elements are computed for all these states. Note that the OCC
and CLOSED
cards must be the same for all states used in a TRANLS
calculation. TRANLS
can also be used with records saved by the MULTI
program (e.g. multi,cisave=record1.file
).
The selection rules for the $M_S$ values are $\Delta M_S = \pm 1$ for the LSX
and LSY
operators, and $\Delta M_S=0$ for the LSZ
operator. Note that $2M_S$ has to be specified, and so the selection rules applying to the difference of the input values are 0 or 2.
In all-electron SO calculations the value of the calculated spin-orbit matrix element is saved (in atomic units) in the MOLPRO
variables TRLSX
, TRLSY
and TRLSZ
for the $x$, $y$, and $z$ components respectively. For ECP-LS calculations the variables TRECPLSX
, TRECPLSY
, and TRECPLSZ
are used. Note that for imaginary matrix elements (i.e., for the $x$ and $z$ components of the SO Hamiltonian) the matrix elements are imaginary and the stored real values have to be multiplied by $i$. If matrix elements for several states are computed, all values are stored in the respective variable-arrays with the bra-states running fastest.
Approximations used in calculating spin-orbit integrals and matrix elements
Recently, more sophisticated approximations were introduced to simplify spin-orbit calculations for larger molecules. These are controlled by specifying the spin-orbit operator type lsop as follows (we omit suffixes X, Y, Z
which specify the component):
LS
Standard spin-orbit calculations.ALS
The one-center approximation is used for one- and two-electron spin-orbit integrals.FLS
The effective Fock-matrix approximation is used for the internal part too.AFLS|AMFI
The one-center approximation is used for one- and two-electron spin-orbit integrals, and the effective Fock-matrix approximation for the internal part.ECPLS
Effective core potentials are used for all atoms at which they are defined; contributions of all other atoms are neglected (see below).
In case that the effective Fock matrix is used for all contributions, and no spin-orbit integrals are pre-calculated and stored on disk (i.e., the LSINT
command is not given), the Fock matrices are evaluated in direct mode and no integrals are stored on disk. When this is combined with the one-center approximation (AMFI), the computing and I/O times are drastically reduced, and this makes spin-orbit calculations quite fast even for larger molecules.
Also, the treatment of ECP-type of spin-orbit interaction has been changed and now allows for treating both ECP and non-ECP atoms in one calculation. Thus, in molecules containing both heavy and light atoms, the heavy atoms can be described using ECPs and the light atoms using all-electron basis sets. If the operator type is LS
, ALS
, FLS
, or AFLS
, then for the atoms having an ECP spin-orbit operator defined in the basis input the ECP operator is used, while the full BP-operator is used for all other atoms (couplings are neglected). Both one-center and AMFI approximations can be used in this case. If, on the other hand, one specifies the operator type as ECPLS
, then the behavior is the same as in the previous versions, i.e., only the ECP contributions are considered and the contributions from all other atoms are neglected.
Calculation and diagonalization of the entire SO-matrix
HLSMAT
,type, record1, record2, record3, …
Computes the entire SO matrix and diagonalizes it using all states which are contained in the records record1, record2, record3, …. All records must have been generated using the SAVE
directive of the of the MULTI or MRCI programs. type may be either LS
for Breit-Pauli calculations, or ECP
for ECP-LS calculations. By default, the eigenvalues and dipole transition matrix elements between the ground and excited states are printed.
As with the TRANLS
card, the HLSMAT
is recognized only by the MRCI program and must be preceded by a CI
card (this is no longer necessary with molpro2021.3 or later). Also, the OCC
and CLOSED
cards must be the same for all states used in a HLSMAT
calculation.
Modifying the unperturbed energies
Often it may be sufficient to compute the spin-orbit matrix elements in a smaller basis or at a lower computational level than the energies. It is therefore possible to replace the energy eigenvalues by precomputed values, which are passed to the spin-orbit program by the MOLPRO
variable HLSDIAG
. The energy values in HLSDIAG
must be in exactly the same order as the states in the records given on the HLSMAT
card. Before any spin-orbit calculation, the variable HLSDIAG
must either be undefined or cleared (then the original energies are used), or must contain exactly the number of energies as the number of states treated in the subsequent spin-orbit calculation (use CLEAR,HLSDIAG
to clear any previous values in the variable). It is the user’s responsibility that the order of the energies in HLSDIAG
is correct!
See example in section SO calculation for the S-atom using the BP operator
Print Options for spin-orbit calculations
PRINT
,option$_1$=value$_1$, option$_2$=value$_2, \ldots$
where option can be
HLS
HLS=-1
only the SO energies and transition matrix elements between ground and excited states are printed (default).
HLS
$\ge 0$: The SO matrix is printed.
HLS
$\ge 1$: The property matrices are printed.
HLS
$\ge 2$: The individual matrix elements are printed (same as OPTION,MATEL
).
HLS
$\ge 3$: Debugging information is printed.
VLS
VLS=-1
: No print of eigenvectors (default).
VLS
$\ge 0$: The eigenvectors are printed.
Options for spin-orbit calculations
Some options can be set using the OPTION
directive (in any order)
OPTIONS
[,WIGNER
=value] [,HLSTRANS
=value] [,MATEL
=value]
where
WIGNER
This option determines whether the Wigner-Eckart theorem should be used when the SO matrix is determined.WIGNER=1
(default) uses the theorem,WIGNER=0
calculates each SO matrix element individually. This option is needed for test purposes only.HLSTRANS
This option determines whether a SO matrix calculation should be performed in the not spin-symmetry adapted basis set (HLSTRANS=0
), in the spin-symmetry adapted basis set (HLSTRANS=1
, default) or with both basis sets (HLSTRANS=2
). At present, symmetry adaption can only be performed for triplet states, where the following notation is used to indicate the symmetry adapted spin functions: $|S,M_S\rangle_+ = \frac{1}{\sqrt{2}} (|S,M_S\rangle + |S,-M_S\rangle)$, $|S,M_S\rangle_- = \frac{1}{\sqrt{2}} (|S,M_S\rangle - |S,-M_S\rangle)$. If only singlet and triplet states are considered, the spin-orbit matrix is blocked according to double-group symmetry and the eigenvalues for each each block are printed separately. In all other cases theHLSTRANS
option is ignored.MATEL
If the entire SO matrix is calculated usingHLSMAT
, the individual matrix elements are normally not shown. When the optionMATEL=1
is given, the individual matrix elements and the contributions of the internal and external configuration classes are printed.
Examples
SO calculation for the S-atom using the BP operator
- examples/s_so.inp
***,SO calculation for the S-atom geometry={s} basis={spd,s,vtz} !use uncontracted basis {rhf;occ,3,2,2,,2;wf,16,4,2} !rhf for 3P state {multi !casscf wf,16,4,2;wf,16,6,2;wf,16,7,2; !3P states wf,16,1,0;state,3;wf,16,4,0;wf,16,6,0;wf,16,7,0} !1D and 1S states {ci;wf,16,1,0;save,3010.1;state,3;noexc} !save casscf wavefunctions using mrci {ci;wf,16,4,0;save,3040.1;noexc} {ci;wf,16,6,0;save,3060.1;noexc} {ci;wf,16,7,0;save,3070.1;noexc} {ci;wf,16,4,2;save,3042.1;noexc} {ci;wf,16,6,2;save,3062.1;noexc} {ci;wf,16,7,2;save,3072.1;noexc} {ci;wf,16,1,0;save,4010.1;state,3} !mrci calculations for 1D, 1S states ed=energy(1) !save energy for 1D state in variable ed es=energy(3) !save energy for 1S state in variable es {ci;wf,16,4,2;save,4042.1} !mrci calculations for 3P states ep=energy !save energy for 3P state in variable ep {ci;wf,16,6,2;save,4062.1} !mrci calculations for 3P states {ci;wf,16,7,2;save,4072.1} !mrci calculations for 3P states text,only triplet states, casscf lsint !compute so integrals text,3P states, casscf {ci;hlsmat,ls,3042.1,3062.1,3072.1} !Only triplet states, casscf text,3P states, mrci {ci;hlsmat,ls,4042.1,4062.1,4072.1} !Only triplet states, mrci text,3P, 1D, 1S states, casscf {ci;hlsmat,ls,3010.1,3040.1,3060.1,3070.1,3042.1,3062.1,3072.1} !All states, casscf text,only triplet states, use mrci energies and casscf SO-matrix elements hlsdiag=[ed,ed,es,ed,ed,ed,ep,ep,ep] !set variable hlsdiag to mrci energies {ci;hlsmat,ls,3010.1,3040.1,3060.1,3070.1,3042.1,3062.1,3072.1}
SO calculation for the I-atom using ECPs
- examples/i_ecp.inp
***,I gprint,orbitals,civector,basis; gthresh,energy=1.d-8,coeff=1.d-8; geometry={I}; basis={ ! ! Iodine-ECP (Dirac-Fock) with SO-coupling ! ecp,I,46,4,3; 1; 2, 1.00000000, 0.00000000; ! lokal term = 0 2; 2, 3.50642001, 83.09814545; 2, 1.74736492, 5.06370919; ! s-terme 4; 2, 2.99860773, 1/3* 81.88444526; 2, 3.01690894, 2/3* 83.41280402; ! p-terms with weights 2, 1.59415934, 1/3* 2.32392477; 2, 1.19802939, 2/3* 2.72079843; 4; 2, 1.03813792, 2/5* 6.40131754; 2, 1.01158599, 3/5* 6.21328827; ! d-terms with weights 2, 2.04193864, 2/5* 19.11604172; 2, 1.99631017, 3/5* 19.08465909; 4; 2, 2.64971585,-3/7* 24.79106489; 2, 2.75335574,-4/7* 24.98147319; ! f-terms with weights 2, 0.49970082,-3/7* 0.27936581; 2, 0.79638982,-4/7* 0.70184261; 4; 2, 2.99860773,-2/3* 81.88444526; 2, 3.01690894, 2/3* 83.41280402; ! ECP-SO for p-terms 2, 1.59415934,-2/3* 2.32392477; 2, 1.19802939, 2/3* 2.72079843; 4; 2, 1.03813792,-2/5* 6.40131754; 2, 1.01158599, 2/5* 6.21328827; ! ECP-SO for d-terms 2, 2.04193864,-2/5* 19.11604172; 2, 1.99631017, 2/5* 19.08465909; 4; 2, 2.64971585, 2/7* 24.79106489; 2, 2.75335574,-2/7* 24.98147319; ! ECP-SO for f-terms 2, 0.49970082, 2/7* 0.27936581; 2, 0.79638982,-2/7* 0.70184261; ! ! Iodine-basis ! s,I,0.2027624,0.4080619,0.8212297,1.6527350,3.3261500; c,1.5,-0.4782372,-0.5811680,0.2617769,0.4444120,-0.1596560; s,I,0.05,0.1007509; p,I,0.2027624,0.4080619,0.8212297,1.6527350,3.3261500; c,1.5,0.4251859,0.2995618,0.0303167,-0.2064228,0.0450858; p,I,0.05,0.1007509,0.01; ! diffuse p-Funktion wegen evt. neg. Part.Ldg d,I,0.2,0.4; f,I,0.3; } {hf;occ,1,1,1,,1;wf,7,5,1} !scf for 2Pz {multi;occ,1,1,1,,1; !casscf with minmal active space wf,7,2,1;wf,7,3,1;wf,7,5,1} !average 2P states {ci;wf,7,2,1;noexc;save,5000.2} !save casscf vector for 2Px state {ci;wf,7,3,1;noexc;save,5100.2} !save casscf vector for 2Py state {ci;wf,7,5,1;noexc;save,5200.2} !save casscf vector for 2Pz state {ci;wf,7,2,1;save,6000.2} !mrci for 2Px state {ci;wf,7,3,1;save,6100.2} !mrci for 2Py state {ci;wf,7,5,1;save,6200.2} !mrci for 2Pz state {multi;occ,1,2,2,,2 !casscf with larger active space wf,7,2,1;wf,7,3,1;wf,7,5,1} !average 2P states {ci;wf,7,2,1;noexc;save,5010.2} {ci;wf,7,3,1;noexc;save,5110.2} {ci;wf,7,5,1;noexc;save,5210.2} {ci;wf,7,2,1;save,6010.2} {ci;wf,7,3,1;save,6110.2} {ci;wf,7,5,1;save,6210.2} text,casscf, occ,1,1,1,,1 {ci;hlsmat,ecp,5000.2,5100.2,5200.2} !do spin-orbit calculations text,casscf, occ,1,2,2,,2 {ci;hlsmat,ecp,5010.2,5110.2,5210.2} text,mrci, occ,1,1,1,,1 {ci;hlsmat,ecp,6000.2,6100.2,6200.2} text,mrci, occ,1,2,2,,2 {ci;hlsmat,ecp,6010.2,6110.2,6210.2}
- examples/i_ecp_cpp.inp
***,I gprint,orbitals,civector,basis; gthresh,energy=1.d-8,coeff=1.d-8; geometry={I}; basis={ ! ! Iodine-ECP (Dirac-Fock) with SO-coupling and cpp ! ecp,I,46,4,3; 1; 2, 1.00000000, 0.00000000; ! lokal term = 0 2; 2, 3.50642001, 83.09814545; 2, 1.74736492, 5.06370919; ! s-terme 4; 2, 2.99860773, 1/3* 81.88444526; 2, 3.01690894, 2/3* 83.41280402; ! p-terms with weights 2, 1.59415934, 1/3* 2.32392477; 2, 1.19802939, 2/3* 2.72079843; 4; 2, 1.03813792, 2/5* 6.40131754; 2, 1.01158599, 3/5* 6.21328827; ! d-terms with weights 2, 2.04193864, 2/5* 19.11604172; 2, 1.99631017, 3/5* 19.08465909; 4; 2, 2.64971585,-3/7* 24.79106489; 2, 2.75335574,-4/7* 24.98147319; ! f-terms with weights 2, 0.49970082,-3/7* 0.27936581; 2, 0.79638982,-4/7* 0.70184261; 4; 2, 2.99860773,-2/3* 81.88444526; 2, 3.01690894, 2/3* 83.41280402; ! ECP-SO for p-terms 2, 1.59415934,-2/3* 2.32392477; 2, 1.19802939, 2/3* 2.72079843; 4; 2, 1.03813792,-2/5* 6.40131754; 2, 1.01158599, 2/5* 6.21328827; ! ECP-SO for d-terms 2, 2.04193864,-2/5* 19.11604172; 2, 1.99631017, 2/5* 19.08465909; 4; 2, 2.64971585, 2/7* 24.79106489; 2, 2.75335574,-2/7* 24.98147319; ! ECP-SO for f-terms 2, 0.49970082, 2/7* 0.27936581; 2, 0.79638982,-2/7* 0.70184261; ! ! Iodine-basis ! s,I,0.2027624,0.4080619,0.8212297,1.6527350,3.3261500; c,1.5,-0.4782372,-0.5811680,0.2617769,0.4444120,-0.1596560; s,I,0.05,0.1007509; p,I,0.2027624,0.4080619,0.8212297,1.6527350,3.3261500; c,1.5,0.4251859,0.2995618,0.0303167,-0.2064228,0.0450858; p,I,0.05,0.1007509,0.01; ! diffuse p-Funktion wegen evt. neg. Part.Ldg d,I,0.2,0.4; f,I,0.3; } {cpp,init,1; ! core polarization potential I,2,1.028,,,1.23} ! Iod-Atom,form of cut-off function, static polarizability ! 1.23 = exponential-factor of cut-off function {hf;occ,1,1,1,,1;wf,7,5,1} !scf for 2Pz {multi;occ,1,1,1,,1; !casscf with minmal active space wf,7,2,1;wf,7,3,1;wf,7,5,1} !average 2P states {ci;wf,7,2,1;noexc;save,5000.2} !save casscf vector for 2Px state {ci;wf,7,3,1;noexc;save,5100.2} !save casscf vector for 2Py state {ci;wf,7,5,1;noexc;save,5200.2} !save casscf vector for 2Pz state {ci;wf,7,2,1;save,6000.2} !mrci for 2Px state {ci;wf,7,3,1;save,6100.2} !mrci for 2Py state {ci;wf,7,5,1;save,6200.2} !mrci for 2Pz state {multi;occ,1,2,2,,2 !casscf with larger active space wf,7,2,1;wf,7,3,1;wf,7,5,1} !average 2P states {ci;wf,7,2,1;noexc;save,5010.2} {ci;wf,7,3,1;noexc;save,5110.2} {ci;wf,7,5,1;noexc;save,5210.2} {ci;wf,7,2,1;save,6010.2} {ci;wf,7,3,1;save,6110.2} {ci;wf,7,5,1;save,6210.2} text,casscf, occ,1,1,1,,1 {ci;hlsmat,ecp,5000.2,5100.2,5200.2} !do spin-orbit calculations text,casscf, occ,1,2,2,,2 {ci;hlsmat,ecp,5010.2,5110.2,5210.2} text,mrci, occ,1,1,1,,1 {ci;hlsmat,ecp,6000.2,6100.2,6200.2} text,mrci, occ,1,2,2,,2 {ci;hlsmat,ecp,6010.2,6110.2,6210.2}