Spin-orbit-coupling

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.

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.

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.

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.

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.

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,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.

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 the HLSTRANS option is ignored.
• MATEL If the entire SO matrix is calculated using HLSMAT, the individual matrix elements are normally not shown. When the option MATEL=1 is given, the individual matrix elements and the contributions of the internal and external configuration classes are printed.

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}

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}