Transition Moments and Spin-Orbit Coupling
A review of various ways of computing spin-orbit coupling:
    D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon,
    Int.Rev.Phys.Chem. 22, 551-592(2003)
    GAMESS can compute transition moments and oscillator 
strengths for the radiative transitions between states 
written in terms of CI wavefunctions (GUGA only).  The 
transition moments are computed using both the "length 
form" and the "velocity form". In a.u., where h-bar=m=1, we 
start from
    [A,q] = -i dA/dp
For A=H, dH/dp=p, and p= -i d/dq,
    [H,q] = -i p = -d/dq.
For non-degenerate states,
    (Ea-Eb) = - 
This relates the dipole to the velocity form,
     = -1/(Ea-Eb) 
but the CI states will give different numbers for each 
side, since the states aren't exact eigenfunctions.  
Transition moment computation is OPERAT=DM in $TRANST.  For 
transition moments, the CI is necessarily performed on 
states of the same multiplicity.

    All other operators are various spin-orbit coupling
options.  There are two kinds of calculations possible,
which we will call SO-CI and SO-MCQDPT.  Note that there
is a hyphen in "spin-orbit CI" to avoid confusion with
"second order CI" in the sense of the SOCI keyword in $DRT
input.  For SO-CI, the initial states may be any CI wave-
function that the GUGA package can generate.  For SO-MCQDPT
the initial states for spin-orbit coupling are of CAS type,
and the operator mixing them corresponds to MCQDPT
generalised for spin-dependent operators (with certain

    GAMESS can compute the "microscopic Breit-Pauli
spin-orbit operator", which includes both a one and two
electron operator.  Additional information is given in a 
subsection below


    For transition moments, the states are generated by
CI calculations using the GUGA package.  These states are
the final states, and the results are just the transition
moments between these states.  The states are defined by
$DRTx input groups.

    For SO-CI, the energy of the CI states forms the
diagonal of a spin-orbit Hamiltonian, as in the state basis
the spin-free Hamiltonian is of course diagonal.  Addition
of the Pauli-Breit operator does not change the diagonal, 
but does add off-diagonal H-so elements.  For SO-MCQDPT,
the spin-free MCQDPT matrix elements are expanded into
matrices corresponding to all Ms values for a pair of
multiplicities.  These matrices are block-diagonal before
the addition of spin-orbit coupling terms, coupling Ms
values.  The diagonalization of this spin-orbit Hamiltonian
gives new energy levels, and spin-mixed final states.
Optionally, the transition dipoles between the final states
can be computed.  The input requirements are $DRTx or
$MCQDx groups which define the original pure spin states.

    We will call the initial states CAS-CI, since most of
the time they will be MCSCF states.  There may be cases
such as the Na example below where SCF orbitals are used,
or other cases where a FOCI or SOCI wavefunction might be
used for the initial states.  Please keep in mind that the
term does not imply the states must be MCSCF states, just
that they commonly are.

    In the above, x may vary from 1 to 64.  The reason for
allowing such a large range is to permit the use of Abelian
point group symmetry during the generation of the initial
states.  The best explanation will be an example, but the
number of these input groups depends on both the number of
orbital sets input, and how much symmetry is present.  The
next two subsections discuss these points.


    The orbitals for transition moments or for SO-CI can be
one common set of orbitals used by all CI states.  If one
set of orbitals is used, the transition moment or spin-
orbit coupling can be found for any type of GUGA CI wave-
function.  Alternatively, two sets of orbitals (obtained by
separate MCSCF orbital optimizations) can be used.  Two or
more separate CIs will be carried out.  The two MO sets
must share a common set of frozen core orbitals, and the
CI must be of the complete active space type.  These
restrictions are needed to leave the CI wavefunctions
invariant under the necessary rotation to corresponding
orbitals.  The non-orthogonal procedure implemented is a
GUGA driven equivalent to the method of Lengsfield, et al.
Note that the FOCI and SOCI methods described by these
workers are not available in GAMESS.

    If you would like to use separate orbitals during the
CI, they may be generated with the FCORE option in $MCSCF.
Typically you would optimize the ground state completely, 
then use these MCSCF orbitals in an optimization of the
excited state, under the constraint of FCORE=.TRUE.

    For SO-MCQDPT calculations, only one set of orbitals
may be input to describe all CAS-CI states.  Typically that
orbital set will be obtained by state-averaged MCSCF, see
WSTATE in $DET/$DRT, and also in the $MCQDx input.  Note
that although the RUNTYP=TRANSITN driver is tied to the
GUGA CI package, there is no reason the orbitals cannot be
obtained using the determinant CI package.  In fact, for
the case of spin-orbit coupling, you might want to utilize
the ability to state average over several spins, see PURES
in $DET.
    If there is no molecular symmetry present, transition
moment calculations will provide $DRT1 if there is one set
of orbitals, otherwise $DRT1 defines the CI based on $VEC1
and $DRT2 the CI based on $VEC2.  Also for the case of no
symmetry, a spin-orbit job should enter one $DRTx or $MCQDx
for every spin multiplicity, and all states of the same
multiplicity have to be generated from $VEC1 or $VEC2,
according to IVEX input.  


    The CAS-CI states are most efficiently generated using
symmetry, since states of different symmetry have zero
Hamiltonian matrix elements.  It is probably more efficient
to do four CI calculations in the group C2v on A1, A2, B1, 
and B2 symmetry, than one CI with a combined Hamiltonian
in C1 symmetry (unless the active space is very small), and
similar remarks apply to the SO-MCQDPT case.  In order to
avoid repeatedly saying $DRTx or $MCQDx, the following few
paragraphs say $DRTx only.

    Again supposing the group is C2v, and you are 
interested in singlet-triplet coupling.  After some 
preliminary CI calculations, you might know that the lowest 
8 states are two 1-a1, 1-b1, two 1-b2, one 3-a1, and two 3-
b2 states.   In this case your input would consist of five 
$DRTx, of  which you can give the three singlets in any 
order but these must preceed the two triplet input groups 
to follow the rule of increasing multiplicity.  Clearly it 
is not possible to write a formula for how many $DRTx there 
will be, this depends not only on the point group, but also 
the chemistry of the situation.

    If you are using two sets of orbitals, the generation
of the corresponding orbitals for the two sets will permute
the active orbitals in an unpredictable way.  Use STSYM to
define the desired state symmetry, rather than relying on
the orbital order.  It is easy and safer to be explicit
about the spatial orbital symmetry.

    The users are encouraged to specify full symmetry in
their $DATA input even though they may choose to set the
symmetry in $DRTx to C1.  The CI states will be labelled in
the group given in $DATA.  The use of non-Abelian symmetry
is limited by the absence of non-Abelian CI or MCQDPT.  In
this case the users can choose between setting full non-
Abelian symmetry in $DATA and C1 in $DRT or else an Abelian
subgroup in both $DATA and $DRT.  The latter choice appears
to be most efficient at present.  

    An example of SO-MCQDPT illustrating how the carbon 
atom of Kh symmetry (full rotation-reflection group) can be 
entered in D2h, Kh's highest Abelian group. The run time is 
considerably longer in C1 symmetry.

    As another example, consider an organic molecule with a
singly excited state, where that state might be coupled to
low or high spin, and where these two states might be close
enough to have a strong spin-orbit coupling.  If it happens
that the S1 and S0 states possess different symmetry, a
very reaasonable calculation would be to treat the S1 and 
T1 state with the same $VEC2 orbitals, leaving the ground
state described by $VEC1.  After doing an MCSCF on the S0
ground state for $VEC1, you could do a state-averaged MCSCF
for $VEC2 optimized for T1 and S1 simultaneously, using
PURES.  The spin orbit job would obtain its initial states
from three GUGA CI computations, S0 from $VEC1 and $DRT1,
S1 from $VEC2 and $DRT2, and T1 from $VEC2 and $DRT3.  Your
$TRANST would be NUMCI=3, IROOTS(1)=1,1,1, IVEX(1)=1,2,2.
Note that the second IROOTS value is 1 because S1 was
presumed to have a different symmetry than S0, so STSYM in
$DRT1 and $DRT2 will differ.  The calculation just outlined
cannot be done if S0 and S1 have the same spatial symmetry,
as IROOTS(1)=1,2,1 to obtain S1 during the second CI will
bring in an additional S0 state (one expressed in terms of
the $VEC2, at slightly higher energy).  This problem is the
origin of the statement several paragraphs above that a 
system with no symmetry will have one $DRTx for every spin
multiplicity included. 

    For transition moments, which do not diagonalize a 
matrix containing these duplicated states, it is OK to
proceed, provided you ignore all transition moments between
the same states obtained in the two different CIs.

spin orbit coupling

    Spin-orbit coupling calculations are always performed 
in a quasi-degenerate perturbative manner.  Typically the 
states close in energy are included into the spin-orbit 
coupling matrix. "Close" has a easily understandable 
meaning, since in the limit of small coupling the quasi-
degenerate treatment is reduced to a second order 
perturbative treatment, that is, the affect of a state upon 
the state of primary interest is given by the square of the 
spin-orbit coupling matrix element divided by the 
difference of the adiabatic energies.  This is useful to 
keep in mind when deciding how many CI states to include in 
the matrix.  The states that are included are treated in a 
fashion that is equivalent to infinite order perturbation 
theory (exact) whereas the states that are not included 
make no contribution.

    Spin-orbit runs can be done for even or odd numbers of 
electrons (any spin), for more than two different spin 
multiplicities at once, for general active spaces.  At 
times, when the spatial wavefunction is degenerate, a spin-
orbit run might involve only one spin multiplicity, e.g. a 
triplet-pi state in a linear molecule.  The most common 
case is two different spins, with non-zero spin orbit 
coupling possible only for delta-S=1: singlets spin-orbit 
couple with triplets, but not with quintets.  Use of three 
spins, such as S=0,1,2, will generate couplings between 
singlets and triplets, and between triplets and quintets, 
which together engender an indirect singlet/quintet mixing.

    As noted above, the treatment of spin-orbit involves 
first obtaining a handful of spin-pure states, whose 
energies form the diagonal of a model Hamiltonian.  The 
spin-orbit operator introduces off-diagonal couplings, and 
the resulting small Hamiltonian is diagonalized.  This 
generates spin-mixed states in the model space.  Since the 
model states are not fully relaxed (internally contracted), 
this is essentially a perturbative treatment:  certainly 
the spin-orbit effects have no influence on orbital 
optimizations or potential energy surfaces when treated in 
this manner, at the very last stage.

    The Breit-Pauli spin-orbit operator contains a one 
electron term arising from Pauli's reduction of the Dirac 
hydrogenic equation to a single-component form, and a two 
electron term due to Breit.  Computation of the full Breit-
Pauli operator is OPERAT=HSO2 (or HSOFF).  A close 
approximation to the latter is HSO2P (P=partial), which 
neglects all active-active two electron terms, which 
usually do not contribute very much to the total coupling, 
while saving substantial computer time.

    HSO1 completely omits the two electron terms, so is 
much faster than any of the two electron operators, but 
represents a potentially much greater loss of accuracy.  
HSO1's error can be remedied to some extent by regarding 
the nuclear charge in the one electron term as an 
adjustable parameter.  In addition, these effective charges 
are often used to compensate for missing nodes in valence 
orbitals of ECP runs, in which case the ZEFF are typically 
very far from the two nuclear charges.  ZEFTYP selects some 
built in values obtained by S.Koseki et al, but if you have 
some favorite parameters, they can be read in as the ZEFF 
input array.  Effective charges may be used for any OPERAT, 
but are most often used with HSO1.

    Theoretical considerations indicate that the Breit-
Pauli operator is not variationally stable, if it were to 
be used during the SCF iterations determining orbitals.  
However, a first order Douglas-Kroll type correction to the 
one electron part of the operator reduces its size by means 
of certain kinematic factors, removing this problem.  This 
can produce substantially better results, even if the 
operator is being treated pertubatively.  This correction 
(at first order) is automatically applied to the one 
electron part of the Breit-Pauli operator by any run that 
selects RELWFN=DK (at any order) in the scalar relativity 
treatment during the variational steps prior to the spin-
orbit perturbation.  See paper 32 below.

    Because the diagonalization of the model spin-orbit 
Hamiltonian leads to spin-mixed eigenvectors, approximate 
wavefunctions including spin-orbit coupling are generated.  
It is now possible to generate the density matrix for the 
spin-mixed states, so that property values for spin-mixed 
states can be found: see keyword ISTNO.  The natural 
orbitals of these spin-orbit density matrices turn out to 
be good approximations to the two spinors of the large 
components of full Dirac four component runs.  The total 
density of these spinors can be obtained for interpretation 
purposes.  Recognition that the spin-orbit coupling is a 
rotational operator (L dot S, where L = R cross P) when it 
acts on an orbital can lead to chemical interpretability of 
the spin orbit results.  See papers 40 and 41 below.

  It is also possible to obtain the dipole transition 
moments between the final spin-mixed wavefunctions, which 
of course do not any longer have a rigourous S quantum no.  
When the run is SO-MCQDPT, the transition moment are first 
computed only between CAS states, and then combined with 
the spin-mixed SO-MCQDPT coefficients.

technical matters:

     The only practical limitation on the computation of 
the Breit term is that HSO2FF is limited to 10 active 
orbitals on 32 bit machines, and to about 26 active 
orbitals on 64 bit machines.  The spin-orbit matrix 
elements vanish for |delta-S| > 1, but it is possible to 
include three or more spins in the computation.  Since 
singlets interact with triplets, and triplets interact with 
pentuplets, inclusion of S=0,1,2 simultaneously lets you 
pick up the indirect interaction between singlets and 
pentuplets that the intermediate triplets afford.

    The choice between HSO2 and HSO2FF is very often in
favor of the former.  HSO2 computes the matrix elements in 
CSF basis and then contracts them with CI coefficients, 
whereas HSO2FF uses a generalized density in AO basis 
computed for each pair of states, thus HSO2 is much more 
efficient in case of multiple states given in IROOTS. 
HSO2FF takes less memory for integral storage, thus it can 
be superior in case of small active spaces and large basis 
sets, in part because it does not store 2e SOC integrals on 
disk and secondly, it does not redundantly treat the same 
pair of determinants if they appear in different CSFs.  The 
numerical results with HSO2 and HSO2FF should be identical 
within machine and algorithmic accuracy.

    Various symmetries are used to avoid computing zero
spin-orbit matrix elements.  NOSYM in $TRANST allows some 
control over this: NOSYM=1 gives up point group symmetry
completely, while 2 turns off additional symmetries such 
as spin selection rules.  HSO1,2,2P compute all matrix
elements in a group (i.e. between two $DRTx groups with
fixed Ms(ket)-Ms(bra)) if at least one of them does not
vanish by symmetry, and HSO2P actually avoids computation
for each pair of states if forbidden by symmetry.  Setting
NOSYM=2 will cause HSO2FF to consider the elements between
two singlets, which are always calculated for HSO1,2,2P
when transition dipoles are requested as well.

    SYMTOL has a dramatic effect on the run speed.  This
cutoff is applied to CSF coefficcients, their products,
and these products times CSF orbital overlaps.  The value
permits a tradeoff of accuracy for run time, and since the
error in the spin-orbit properties approaches SYMTOL mainly
for SOCI functions, it may be useful to increase SYMTOL to
save time for CAS or FOCI functions.  Some experimenting
will tell you what you can get away with.  SYMTOL is also
used during CI state symmetry assignment, for NOIRR=-1
in $DRT.

   In case if you do not provide enough storage for the
form factors sorting then some extra disk space will be
used;  the extra disk space can be eliminated if you set
SAVDSK=.TRUE. (the amount of savings depends on the active
space and memory provided, it some cases it can decrease
the disk space up to one order of magnitude).  The form
factors are in binary format, and so can be transfered
between computers only if they have compatible binary
files.  There is a built-in check for consistency of a
restart file DAFL30 with the current run parameters.
input nitty-gritty
    The transition moment and spin-orbit coupling driver
is a rather restricted path through the GUGA CI part of
GAMESS.  Note that $GUESS is not read, instead the MOs will
be MOREAD in a $VEC1 and perhaps a $VEC2 group.  It is not
possible to reorder MOs.  For SO-CI,

2) $CIINP is not read.  The CI is hardwired to consist
   of CI DRT generation, integral transformation/sorting,
   Hamiltonian generation, and diagonalization.  This
   means $DRT1 (and maybe $DRT2,...), $TRANS, $CISORT,
   $GUGEM, and $GUGDIA input is read, and acted upon.
3) The density matrices are not generated, and so no
   properties (other than the transition moment or the
   spin-orbit coupling) are computed.
4) There is no restart capability provided, except for
   saving some form-factor information.
5) $DRT1, $DRT2, $DRT3, ... must go from lowest to highest

6) IROOTS will determine the number of CI states in each
   CI for which the properties are calculated.  Use
   NSTATE to specify the number of CI states for the
   CI Hamiltonian diagonalization.  Sometimes the CI
   convergence is assisted by requesting more roots
   to be found in the diagonalization than you want to
   include in the property calculation.

For SO-MCQDPT, the steps are


2) the number of roots in each MCQDPT is controlled by
   $TRANST's IROOTS, and each such calculation is
   defined by $MCQD1, $MCQD2, ... input.  These must go
   from lowest multiplicity to highest.


The review already mentioned:
"Spin-orbit coupling in molecules: chemistry beyond the
 adiabatic approximation".
D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon,
Int.Rev.Phys.Chem. 22, 551-592(2003)

Reference for separate active orbital optimization:
 1. B.H.Lengsfield, III,  J.A.Jafri,  D.H.Phillips,
    C.W.Bauschlicher, Jr.  J.Chem.Phys. 74,6849-6856(1981)

References for transition moments:
2a. H.C.Longuet-Higgins
    Proc.Roy.Soc.(London) A235, 537-543(1956)
2b. F.Weinhold, J.Chem.Phys. 54,1874-1881(1970)
 3. C.W.Bauschlicher, S.R.Langhoff
    Theoret.Chim.Acta 79:93-103(1991)
 4. "Intermediate Quantum Mechanics, 3rd Ed." Hans A.
    Bethe, Roman Jackiw   Benjamin/Cummings Publishing,
    Menlo Park, CA (1986), chapters 10 and 11.
 5. S.Koseki, M.S.Gordon
    J.Mol.Spectrosc. 123, 392-404(1987)
References for HSO1 spin-orbit coupling, and Zeff values:
 6. S.Koseki, M.W.Schmidt, M.S.Gordon  
    J.Phys.Chem.  96, 10768-10772 (1992)
 7. S.Koseki, M.S.Gordon, M.W.Schmidt, N.Matsunaga
    J.Phys.Chem.  99, 12764-12772 (1995)
 8. N.Matsunaga, S.Koseki, M.S.Gordon
    J.Chem.Phys.  104, 7988-7996 (1996)
 9. S.Koseki, M.W.Schmidt, M.S.Gordon
    J.Phys.Chem.A  102, 10430-10435 (1998)
10. S.Koseki, D.G.Fedorov, M.W.Schmidt, M.S.Gordon
    J.Phys.Chem.A  105, 8262-8268 (2001)
11. S.Koseki, T.Matsushita, M.S.Gordon
    J.Phys.Chem.A  110, 2560-2570(2006)

References for full Breit-Pauli spin-orbit coupling:
20. T.R.Furlani, H.F.King
    J.Chem.Phys.  82, 5577-5583 (1985)
21. H.F.King, T.R.Furlani
    J.Comput.Chem.  9, 771-778 (1988)
22. D.G.Fedorov, M.S.Gordon
    J.Chem.Phys. 112, 5611-5623 (2000)
Paper 22 contains information on the HSO2P partial two 
electron operator method.

Symmetry in spin-orbit coupling:
23. D.G.Fedorov, M.S.Gordon
    ACS Symposium Series 828, pp 1-22(2002)

Reference for SO-MCQDPT:
25. D.G.Fedorov, J.P.Finley  Phys.Rev.A 64, 042502 (2001)

Reference for Spin-Orbit with Model Core Potentials:
30. D.G.Fedorov, M.Klobukowski
    Chem.Phys.Lett. 360, 223-228(2002)
31. T.Zeng, D.G.Fedorov, M.Klobukowski 
    J.Chem.Phys. 131, 124109/1-17(2009)
32. T.Zeng, D.G.Fedorov, M.Klobukowski
    J.Chem.Phys. 132, 074102/1-15(2010)
The last two of these also discuss the 1st order Douglas-
Kroll transformation of the 1e- part of the spin orbit 

Reference for properties, interpretations, and Spin-Orbit 
Natural Spinors:
40. T. Zeng, D. G. Fedorov, M. W. Schmidt, M. Klobukowski
    J. Chem. Phys. 134, 214107/1-9(2011)
41. T. Zeng, D. G. Fedorov, M. W. Schmidt, M. Klobukowski
    J. Chem. Theor. Comput. 7, 2864-2875(2011)
with an application being
42. T.Zeng, D.G.Federov, M.W.Schmidt, M.Klobukowski
    J.Chem.Theo.Comput. 8, 3061-3071(2012).

Recent applications (see also 32,40,41):
50. S.P.Webb, M.S.Gordon
    J.Chem.Phys. 109, 919-927(1998)
51. D.G.Fedorov, M.Evans, Y.Song, M.S.Gordon, C.Y.Ng
    J.Chem.Phys. 111, 6413-6421 (1999)
52. D.G.Fedorov, M.S.Gordon, Y.Song, C.Y.Ng
    J.Chem.Phys. 115, 7393-7400 (2001)
53. B.J.Duke  J.Comput.Chem. 22, 1552-1556 (2001)
54. C.M.Aikens, M.S.Gordon
    J.Phys.Chem.A 107, 104-114(2003)
55. D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon
    K.Hirao and Y.Ishikawa (eds.) Recent Advances in
    Relativistic Molecular Theory, Vol. 5,
    (World Scientific, Singapore), 2004, pp 107-136.

                          * * *

    Special thanks to Bob Cave and Dave Feller for their 
assistance in performing check spin-orbit coupling runs 
with the MELDF programs.  Special thanks to Tom Furlani for 
contributing his 2e- spin-orbit code and answering many 
questions about its interface.  Special thanks to Haruyuki 
Nakano for explaining the spin functions used in the MCQDPT 

    We end with 2 examples.  Note that you must know what
you are doing with term symbols, J quantum numbers, point
group symmetry, and so on in order to make skillful use of
this part of the program.  Seeing your final degeneracies
turn out like a text book says it should is beautiful!
!  Compute the splitting of the famous sodium D line.
!  Joseph von Fraunhofer (Denkschriften der Koeniglichen
!  Akademie der Wissenschf. zu Muenchen, 5, 193(1814-1815))
!  observed the sun through good prisms, finding 700 lines,
!  and named the brightest ones A, B, C... just in order.
!  He was able to resolve the D line into two lines, which
!  occur at 5895.940 and 5889.973 Angstroms.  It would take
!  a century to understand the D line is Na's 3s <-> 3p
!  transition, and that spin-orbit coupling is what splits
!  the D line into two.  Charlotte Moore's Atomic Energy
!  Levels, volume 1, gives the experimental 2-P interval
!  as 17.1963, since the three relevent levels are at 
!  2-S-1/2= 0.0, 2-P-1/2= 16,956.183, 2-P-3/2= 16,973.379.

1. generate ground state 2-S orbitals by conventional ROHF.
   the energy of the ground state is -161.8413919816
--- $contrl scftyp=rohf mult=2 $end
--- $system kdiag=3 memory=300000 $end
--- $guess  guess=huckel $end
2. generate excited state 2-P orbitals, using a state-
averaged SCF wavefunction to ensure radial degeneracy of 
the 3p shell is preserved.  The open shell SCF energy is
-161.7682895801. The computation is both spin and space 
restricted open shell SCF on the 2-P Russell-Saunders term.  
Starting orbitals are reordered orbitals from step 1.
--- $contrl scftyp=gvb mult=2 $end
--- $system kdiag=3 memory=300000 $end
--- $guess  guess=moread norb=13
---         norder=1 iorder(6)=7,8,9,6 $end
--- $scf    nco=5 nseto=1 no(1)=3 rstrct=.t. couple=.true.
---             f(1)=  1.0  0.16666666666667
---         alpha(1)=  2.0  0.33333333333333  0.0
---          beta(1)= -1.0 -0.16666666666667  0.0 $end

3. compute spin-orbit coupling in the 2-P term.  The use of 
C1 symmetry in $DRT1 ensures that all three spatial CSFs 
are kept in the CI function.  In the preliminary CI, the 
spin function is just the alpha spin doublet, and all three 
roots should be degenerate, and furthermore equal to the 
GVB energy at step 2.  The spin-orbit coupling code uses 
both doublet spin functions with each of the three spatial 
wavefunctions, so the spin-orbit Hamiltonian is a 6x6 
matrix.  The two lowest roots of the full 6x6 spin-orbit 
Hamiltonian are the doubly degenerate 2-P-1/2 level, while 
the other four roots are the degenerate 2-P-3/2 level.
 $contrl scftyp=none cityp=guga runtyp=transitn mult=2 $end
 $system memory=2000000 $end
 $basis  gbasis=n31 ngauss=6 $end
 $gugdia nstate=3 $end
 $transt operat=hso1 numvec=1 numci=1 nfzc=5 nocc=8 
         iroots=3 zeff=10.04 $end
 $drt1   group=c1 fors=.true. nfzc=5 nalp=1 nval=2 $end

Na atom...2-P excited state...6-31G basis
Dnh 2

Na 11.0
--- GVB ORBITALS --- GENERATED AT  7:46:08 CST 30-MAY-1996
Na atom...2-P excited state
E(GVB)=     -161.7682895801, E(NUC)=     .0000000000,    5 
 1  1 9.97912679E-01 8.83038094E-03 0.00000000E+00...
      ... orbitals from step 2 go here ...
13  3-1.10674398E+00 0.00000000E+00 0.00000000E+00

   As an example of both SO-MCQDPT, and the use of as much
symmetry as possible, consider carbon.  The CAS-CI uses
an active space of 2s,2p,3s,3p orbitals, and the spin-orbit
job includes all terms from the lowest configuration,
2s2,2p2.  These terms are 3-P, 1-D, and 1-S.  If you look
at table 58 in Herzberg's book on electronic spectra, you
will be able to see how the Kh spatial irreps P, D, S are
partitioned into the D2h irreps input below.

!   C SO-MRMP on all levels in the s**2,p**2 configuration.
!  levels        CAS         and     MCQDPT
!   1           .0000                 .0000 cm-1      3-P-0
!   2-4       12.6879-12.8469       13.2721-13.2722   3-P-1
!   5-9       37.8469-37.8470       39.5638-39.5639   3-P-2
!  10-14   12169.1275            10251.7910           1-D-2
!  15      19264.4221            21111.5130           1-S-0
!   The active space consists of (2s,2p,3s,3p) with 4 e-.
!   D2h symmetry speeds up the calculation considerably,
!   on the same computer D2h = 78 and C1 = 424 seconds.
 $contrl scftyp=none cityp=none mplevl=2
         runtyp=transitn $end
 $system memory=5000000 $end
!            below is input to run in C1 subgroup
--- $transt operat=hso2 numvec=-2 numci=2 nfzc=1 nocc=9
---         iroots(1)=6,3 parmp=3
---         ivex(1)=1,1 $end
--- $mrmp mrpt=mcqdpt rdvecs=.t. $end
--- $MCQD1  nosym=1 nstate=6 mult=1 iforb=3 
---         nmofzc=1 nmodoc=0 nmoact=8
---     wstate(1)=1,1,1,1,1,1 thrcon=1e-8 thrgen=1e-10 $END
--- $MCQD2  nosym=1 nstate=3 mult=3 iforb=3 
---         nmofzc=1 nmodoc=0 nmoact=8
---         wstate(1)=1,1,1 thrcon=1e-8 thrgen=1e-10 $END
!            below is input to run in D2h subgroup
 $transt operat=hso2 numvec=-7 numci=7 nfzc=1 nocc=9
         iroots(1)=3,1,1,1, 1,1,1   parmp=3
         ivex(1)=1,1,1,1,1,1,1 $end
 $mrmp mrpt=mcqdpt rdvecs=.t. $end
 $MCQD1  nosym=-1 mult=1 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=Ag wstate(1)=1,1,1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD2  nosym=-1 mult=1 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B1g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD3  nosym=-1 mult=1 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B2g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD4  nosym=-1 mult=1 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B3g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD5  nosym=-1 mult=3 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B1g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD6  nosym=-1 mult=3 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B2g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
 $MCQD7  nosym=-1 mult=3 iforb=3 
         nmofzc=1 nmodoc=0 nmoact=8 
     stsym=B3g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
!     input  to prepare the 3-P ground state orbitals
!     great care is taken to create symmetry equivalent p's
--- $contrl scftyp=mcscf cityp=none mplevl=0
---         runtyp=energy mult=3 $end
--- $guess  guess=moread norb=55 purify=.t. $end 
--- $mcscf  cistep=guga fullnr=.t. $end
--- $drt    group=c1 fors=.true.
---         nmcc=1 ndoc=1 nalp=2 nval=5 $end
--- $gugdia nstate=9 maxdia=1000 $end
--- $gugdm2 wstate(1)=1,1,1 $end
C...aug-cc-pvtz (10s,5p,2d,1f) -> [4s,3p,2d,1f] 
Dnh 2

C 6.0
 S   8
  1        8236.000000         0.5310000000E-03
  2        1235.000000         0.4108000000E-02
  3        280.8000000         0.2108700000E-01
  4        79.27000000         0.8185300000E-01
  5        25.59000000         0.2348170000
  6        8.997000000         0.4344010000
  7        3.319000000         0.3461290000
  8       0.3643000000        -0.8983000000E-02
 S   8
  1        8236.000000        -0.1130000000E-03
  2        1235.000000        -0.8780000000E-03
  3        280.8000000        -0.4540000000E-02
  4        79.27000000        -0.1813300000E-01
  5        25.59000000        -0.5576000000E-01
  6        8.997000000        -0.1268950000
  7        3.319000000        -0.1703520000
  8       0.3643000000         0.5986840000
 S   1
  1       0.9059000000          1.000000000
 S   1
  1       0.1285000000          1.000000000
 P   3
  1        18.71000000         0.1403100000E-01
  2        4.133000000         0.8686600000E-01
  3        1.200000000         0.2902160000
 P   1
  1       0.3827000000          1.000000000
 P   1
  1       0.1209000000          1.000000000
 D   1
  1        1.097000000          1.000000000
 D   1
  1       0.3180000000          1.000000000
 F   1
  1       0.7610000000          1.000000000
 S   1
  1       0.440200000E-01      1.00000000
 P   1
  1       0.356900000E-01      1.00000000
 D   1
  1       0.100000000          1.00000000
 F   1
  1       0.268000000          1.00000000

E(MCSCF)=      -37.7282408589, 11 ITERS
 1  1 9.75511467E-01 ...snipped...

created on 7/7/2017