How to do MCSCF (and CI) calculations

    Multi-configuration self consistent field (MCSCF) 
wavefunctions are the most general SCF possible.  MCSCF 
allows for a natural description of chemical processes 
involving the separation of electrons (bond breaking, 
electronic excitation, etc), which are often not well 
represented using the single configuration SCF methods.

    MCSCF wavefunctions, as the name implies, contain more 
than one configuration, each of which is multiplied by a 
configuration interaction (CI) coefficient, optimized to 
determine its weight.  In addition, the shapes of the 
orbitals used to form each of the configurations are 
optimized, just as in a simpler SCF, to self consistency.

    Typically every different chemical problem requires 
that an MCSCF wavefunction be designed to treat it, on a 
case by case basis, by choosing an "active space".  For 
example, one may be interested in describing the reactivity 
of a particular functional group, instead of elsewhere in 
the molecule.  So, the active electrons and active orbitals 
will be those that are "active" on that functional group.  
Orbitals elsewhere in the molecule just remain doubly 
occupied, as for RHF.  This means some attention must be 
paid to orbitals in order to obtain the desired results.

    Procedures for the selection of configurations (which 
amounts to choosing the number of active electrons and 
active orbitals), for the two mathematical optimizations 
just mentioned, ways to interpret the resulting MCSCF 
wavefunction, and treatments for dynamical electron 
correlation of MCSCF wavefunctions are the focus of a 
review article:
    "The Construction and Interpretation of MCSCF
         M.W.Schmidt and M.S.Gordon,  
         Annu.Rev.Phys.Chem. 49,233-266(1998)
One section of this article is devoted to the problem of 
designing the correct active space to treat your problem.  
Additional reading is listed at the end of this section.  A 
much newer and very valuable review article is
    "Multiconfiguration self-consistent field and
       multireference configuration interaction
              methods and applications"
 P.G.Szalay, T.Muller, G.Gidofalvi, H.Lischka, R.Shepard
         Chem.Rev. 112, 108-181(2012)

    These pages describe a powerful and mature MCSCF 
program, allowing computation of the MCSCF energy, nuclear 
gradient, and nuclear hessian for pure states.  Practical 
procedures for generating starting orbitals are available.  
Localized orbital analysis of the final active orbitals is 
provided.  State-averaged energies and their analytic 
gradients can be obtained.  Minimum energy crossing points 
and conical intersections between surfaces may be found, 
see RUNTYP=MEX and CONINT.  Non-adiabatic coupling matrix 
elements (NACME) can be computed, see RUNTYP=NACME.  
Diabatic potential energy surfaces can be obtained, see 
DIABAT in $MCSCF.  If desired, spin-orbit couplings or 
transition dipole moments can be found, see elsewhere in 
this chapter.  Efficient perturbative treatments of the 
dynamical correlation energy for all electrons, whether in 
active and filled orbitals, are provided.  Of course, 
parallel computation has been enabled.

    The most efficient technique implemented in GAMESS for 
finding the dynamic correlation energy of MCSCF is second 
order perturbation theory, in the variant known as MCQDPT 
(known as MRMP for one state).  MCQDPT is discussed in a 
different section of this chapter.

    The use of CI, probably in the form of second order CI, 
will be described below, en passant, during discussion of 
the input defining the configurations for MCSCF.  Selection 
of a CI following any type of SCF (except UHF) is made with 
CITYP in the $CONTRL group, and masterminded by $CIINP.

MCSCF implementation

    With the exception of the QUAD converger, the MCSCF 
program is of the type termed "unfolded two-step" by Roos.  
This means the orbital and CI coefficient optimizations are 
separated.  The latter are obtained in a conventional CI 
diagonalization, while the former are optimized by a 
separate orbital improvement step.

    Each MCSCF iteration (except for the JACOBI and QUAD 
convergers) consists of the following steps:
1) transformation of AO integrals to the current MO basis,
2) generation of the Hamiltonian matrix and optimization
   of the CI coefficients by a Davidson diagonalization,
3) generation of the first and second order density matrix, 
4) improvement of the molecular orbitals.
During the first iteration at the first geometry, you will 
receive verbose output from each of these steps, but each 
subsequent iteration produce only a single summary line.

    The CI problem in steps two and three has four options 
for the many electron basis, namely ALDET, ORMAS, or GENCI 
using determinants, or GUGA using CSFs.  This choice is 
made with the keyword CISTEP in $MCSCF.  Much more will be 
said below about the differences between determinants and 
CSFs.  The word "configuration" will be used throughout 
this section to refer to either determinants or CSFs, when 
a generic term is needed for the many-electron functions.  
Most people use CSF and configuration interchangeably, so 
please note the distinction made here.

    The orbital update in step four has five options, 
namely FOCAS, SOSCF, FULLNR, JACOBI, and QUAD, listed here 
in roughly the order of their increasing mathematical 
sophistication, convergence characteristics, and of course, 
their computer resource requirements.  Again, these are 
chosen by keywords in the $MCSCF group.  More will be said 
just below about the relative merits of these.

    Depending on the converger chosen, the program will 
select the appropriate kind of integral transformation at 
step one. There's seldom need to try to fine tune this, but 
note that the $TRANS group allows you to choose an AO 
integral direct transformation, with the DIRTRF flag.

    The type of CI and the type of orbital converger are to 
some extent "mix and match".  This is particularly true for 
the two full CI programs, GUGA or ALDET, where either 
produces exactly the same CI density matrices.  Here is a 
chart of the ways to combine CI and orbital optimizers:
               parallel run's
     orbital   transformation   CI computation via CISTEP=
    converger      memory       GUGA   ALDET  GENCI  ORMAS
    ---------  --------------   ----   -----  -----  -----
     FOCAS       replicated      ok     ok    silly  silly
     SOSCF       replicated      ok     ok     ok     ok
     FULLNR      distributed     ok     ok     ok     ok
     QUAD          serial        ok     xx     xx     xx
     JACOBI        serial        ok     ok     ok     ok
"xx" means QUAD converger is coded only for CISTEP=GUGA.
"silly" means that this converger ignores active-active
     rotations, so these runs are likely to be divergent,
     or perhaps converge to a false solution.
"serial" means this can only run sequentially at present.

    The next two sections provide more information on the 
two mathematical optimizations, first how the orbital shape 
is refined, and then the determinantion of CI coefficients.

orbital updates
    There are presently five orbital improvement options, 
namely FOCAS, SOSCF, FULLNR, JACOBI, and QUAD.  All but the 
JACOBI update run in parallel.  Each converger is discussed 
briefly below, in order of increasing robustness.  The most 
commonly used convergers are SOSCF and FULLNR.

   The input to control the orbital update step is the 
$MCSCF group, where you can pick the convergence procedure.  
Most of the input in this group is rather specialized, but 
note in particular MAXIT and ACURCY, which control the 
convergence behavior.

    FOCAS is a first order, complete active space MCSCF 
optimization procedure.  It is based on a novel approach 
due to Meier and Staemmler, using very fast but numerous 
microiterations to improve the convergence of what is 
intrinsically a first order method.  Since FOCAS requires 
only one virtual orbital in the integral transformation to 
compute the Lagrangian (whose asymmetry is the orbital 
gradient, and must fall below ACURCY at convergence), the 
total MCSCF job may take less time than a second order 
method, even though it may require many more iterations to 
converge.  The use of microiterations is crucial to FOCAS' 
ability to converge.  It is important to take a great deal 
of care choosing the starting orbitals.

    SOSCF is a method built upon the FOCAS code, which 
seeks to combine the speed of FOCAS with approximate second 
order convergence properties.  Thus SOSCF is an approximate 
Newton-Raphson, based on a diagonal guess at the orbital 
hessian, and in fact has much in common with the SOSCF 
option in $SCF.  Its time requirements per iteration are 
like FOCAS, with a convergence rate better than FOCAS but 
not as good as true second order.  Storage of only the 
diagonal of the orbital hessian allows the SOSCF method to 
be used with much larger basis sets than exact second order 
methods.  Because SOSCF usually requires the least CPU 
time, disk space, and memory needs, it is the default.  
Good convergence by the SOSCF method requires that you 
prepare starting orbitals carefully, and read in all MOs in 
$VEC, as providing canonicalized virtual orbitals increases 
the diagonal dominance of the orbital hessian.  Parallel 
computations are possible with SOSCF, but only to a modest 
number of nodes.

    FULLNR means a full Newton-Raphson orbital improvement 
step is taken, using the exact orbital hessian.  FULLNR is 
a robust convergence method, and normally takes the fewest 
iterations to converge.  Computing the exact orbital 
hessian requires two virtual orbital indices be included in 
the integral transformation, making this step quite time 
consuming, and of course memory for storage of the orbital 
hessian must be available.  Because both the transformation 
and orbital improvement steps of FULLNR are time consuming, 
FULLNR is not the default.  You may want to try FULLNR when 
convergence is difficult, assuming you have already tried 
preparing good starting orbitals by the hints below.

    The serial FULLNR code uses the augmented hessian 
matrix approach to solve the Newton-Raphson equations.  
There are two suboptions for computation of the orbital 
hessian: DM2 is faster, but takes more memory than TEI.  
The parallel implementation of FULLNR avoids explicit 
storage of the orbital hessian, by recomputing the product 
of the hessian times orbital rotation vector during the 
subiterations solving the Newton-Raphson problem.  The 
partial integral transformation used to set up the FULLNR 
converger has been changed to use distributed memory, and 
will scale like the MP2 energy/gradient programs, to many 
nodes.  Parallel FULLNR requires large MEMORY only for the 
CI step (if the active space is big), but always requires a 
large MEMDDI.  The parallel FULLNR program is essentially 
diskless, apart from storage of converged CI vectors.

    The JACOBI method uses a series of 2 by 2 orbital 
rotations by an angle predicted to lower the energy.  This 
should essentially ensure convergence after sweeping 
through all possible orbital pairs enough times.  The 
procedure was created to converge selected (general) 
determinant MCSCF functions, but of course it can be used 
will full lists as well in difficult cases.  The JACOBI 
calculation will consist of a full four index 
transformation over all MOs before the iterations begin.  
Each iteration consists of
 1. a small 4 index transformation over active orbitals 
 2. optimization of the CI vector
 3. generation of the 1e- and 2e- density matrices
 4. sweeps over Jacobi rotations, using MO integrals in
    memory to generate each rotation, with a subsequent
    update after each pair is rotated.
 5. when sufficient energy lowering has been achieved,
    begin a new iteration.
This procedure never generates the orbital Lagrangian! 
Unfortunately this means that at present it is not possible 
to compute nuclear gradients. Due to lack of a Lagrangian, 
ACURCY is of course irrelevant, so the convergence test is 

    QUAD uses a fully quadratic, or second order approach 
and is thus the most powerful MCSCF converger.  The QUAD 
code is programmed only for CISTEP=GUGA.  QUAD runs begin 
with normal unfolded FULLNR iterations, until the orbitals 
approach convergence sufficiently.  QUAD then begins the 
simultaneous optimization of CI coefficients and orbitals, 
and convergence should be obtained in 3-4 additional MCSCF 
iterations.  The QUAD method requires building the full 
electronic hessian, including orbital/orbital, orbital/CI, 
and CI/CI blocks, which is a rather big matrix.  In 
principle, this is the most robust method available, but it 
is limited to perhaps 50-100 CSFs only, because it is a 
memory hog.  QUAD may be helpful in converging excited 
electronic states, but note that you may not use state 
averaging with QUAD.  In practice, QUAD has not received 
very much use compared to the unfolded convergers.

CI coefficient optimization

    Determinants or configuration state functions (CSFs) 
may be used to form the many electron basis set.  It is 
necessary to explain these in a bit of detail so that you 
can understand the advantages of each.

   A determinant is a simple object: a product of spin 
orbitals with a given Sz quantum number, that is, the 
number of alpha spins and number of beta spins are a 
constant.  Matrix elements involving determinants are 
correspondingly simple, but unfortunately determinants are 
not necessarily eigenfunctions of the S**2 operator.

    To expand on this point, consider the four familiar 2e- 
functions which satisfy the Pauli principle.  Here u, v are 
space orbitals, and a, b are the alpha and beta spin 
functions.  As you know, the singlet and triplets are:
       S1 = (uv + vu)/sqrt(2) * (ab - ba)/sqrt(2)
       T1 = (uv - vu)/sqrt(2) *  aa
       T2 = (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2)
       T3 = (uv - vu)/sqrt(2) *  bb
It is a simple matter to multiply out S1 and T2, and to 
expand the two determinants which have Sz=0,
       D1 = |ua vb|          D2 = |va ub|
This reveals that
       S1 = (D1+D2)/sqrt(2)   or   D1 = (S1 + T2)/sqrt(2)
       T2 = (D1-D2)/sqrt(2)        D2 = (S1 - T2)/sqrt(2)
Thus, one must take a linear combination of determinants in 
order to have a wavefunction with the desired total spin. 
There are two important points to note:
  a) A two by two Hamiltonian matrix over D1 and D2 has
     eigenfunctions with -different- spins, S=0 and S=1.
  b) use of all determinants with Sz=0 does allow for the
     construction of spin adapted states.  D1+D2, or D1-D2,
     are -not- spin contaminated.
By itself, a determinant such as D1 is said to be "spin 
contaminated", being a fifty-fifty admixture of singlet and 
triplet.  (It is curious that calculations with just one 
such determinant are often called "singlet UHF", when this 
is half triplet!).  Of course, some determinants are spin 
adapted all by themselves, for example the spin adapted 
functions T1 and T3 above are single determinants, as are 
the closed shells
       S2 = (uu) * (ab - ba)/sqrt(2).
       S3 = (vv) * (ab - ba)/sqrt(2).
It is possible to perform a triplet calculation, with no 
singlet states present, by choosing determinants with Sz=1 
such as T1, since then no state with Sz=0 exists in the 
determinant basis set (as is required when S=0).  To 
summarize, the eigenfunctions of a Hamiltonian formed by 
determinants with any particular Sz will be spin states 
with S=Sz, S=Sz+1, S=Sz+2, ... but will not contain any S 
values smaller than Sz.

    CSFs are an antisymmetrized combination of a space 
orbital product, and a spin adapted linear combination of 
simple alpha-beta products.  Namely, the following CSF
       C1 = A (uv) * (ab-ba)/sqrt(2)
which has a singlet spin function is identical to S1 above 
if you write out what the antisymmetrizer A does, and the 
       C2 = A (uv) * aa
       C3 = A (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2)
       C4 = A (uv) * bb
equal T1-T3.  Since the three triplet CSFs have the same 
energy, GAMESS works with the simpler form C2.  Singlet and 
triplet computations using CSFs are done in separate runs, 
because when spin-orbit coupling is not considered, the 
Hamiltonian is block diagonal in a CSF basis.  Technical 
information about the CSFs is that they use Yamanouchi-
Kotani spin couplings, and matrix elements are obtained 
using a GUGA, or graphical unitary group approach.  
    Determinant and CSF are both primarily used for MCSCF 
wavefunctions, but can be used in CI (see CITYP in 
$CONTRL).  Other comparisons between the determinant and 
CSF implementations, as they exist in GAMESS today, are
                             determinants      CSFs
    parallel execution           yes            yes
    direct CI                    yes             no
    use Abelian group symmetry   yes            yes
    state average mixed spins    yes             no
    first order density          yes            yes
    state averaged densities     yes            yes
    analytic nuclear hessian     yes             no
    can form CI Lagrangian        no            yes
In nearly every circumstance the determinant CI will run 
faster than GUGA, so it is the default.  Here are timings 
for N electrons in N orbitals, no symmetry used:
    N in N    ALDET      GENCI   --- GUGA ---
       8          1         1        1      0
      10          8        38       19     33
      12        228      3122      534   2209
      14       7985        --    15377 130855
The reason there are two numbers under GUGA is that the 
first is for writing the loops (Hamiltonian data) to disk, 
and the second is for the actual diagonalization.  Note 
that the GUGA Hamiltonian time alone is greater than the 
entire ALDET computation!  The ALDET program does not store 
anything on disk, and so runs at CPU/wall ratios of 1.  The 
quality of the initial guess of the CI eigenvector in the 
various determinant codes is much better than in the GUGA 
CSF code, so the chances of converging to an incorrect 
excited state root is much less.  Finally, determinant 
input is generally easier.  No wonder determinants are now 
the default configurations!

    Two of the determinant CI programs, namely ALDET or 
ORMAS (but not GENCI) have been changed to use replicated 
memory parallelism, with modest scaling.  The GUGA program 
is parallel in its solving, but not its Hamiltonian 
generation, and as already noted, its H formation takes 
more time than the entire determinant CI (translation: use 
CISTEP=ALDET, not CISTEP=GUGA).  The GENCI program will run 
as a serial bottleneck (no speedup) in parallel runs.

    The next two sections describe in detail the input for 
specification of the configurations, either determinants or 

determinant CI
    Three determinant CI codes are provided for MCSCF, one 
for full CI spaces (ALDET), another named the Occupation 
Restricted Multiple Active Spaces (ORMAS), and finally 
there is a program for arbitrary (selected) determinant 
lists (GENCI).  For straight CI, but not MCSCF, there is a 
fourth program, the full second order CI (CITYP=FSOCI), 
whose purpose is MR-CISD.

    ALDET is a full CI within the chosen active space.  It 
is possible to go up to 16 electrons in 16 orbitals, if 
your computer has a lot of memory.  ALDET is the only 
CISTEP for which analytic nuclear hessian is possible, and 
it is also the most scalable CI code (using replicated 
memory to store CI vectors).  A sample input for ALDET is
Keywords in this group actually relate to all determinant 
programs, and are described below.

    The $DET input group is basic to all determinant CI 
codes.  Keywords GROUP and STSYM specify the desired 
spatial symmetry of the determinants.  Most runs need give 
only the orbital and electron counts:  NCORE, NACT, and 
NELS.  The number of electrons is 2*NCORE+NELS, and will be 
checked against the charge implied by ICHARG.  The MULT 
given in $CONTRL is used to determine the desired Sz value, 
by extracting S from MULT=2S+1, then by default Sz=S.  If 
you wish to include lower spin multiplicities, which will 
increase the CPU time of the run, but will let you know 
what the energies of such states are, just input a smaller 
value for SZ.  The states whose orbitals will be MCSCF 
optimized will be those having the requested MULT value, 
unless you choose otherwise with the PURES flag.

    The remaining parameters in the $DET group give extra 
control over the diagonalization process.  Most are not 
given in normal circumstances, except NSTATE, which you may 
need to adjust to produce enough roots of the desired MULT 
value.  The only important keyword which has not been 
discussed is the WSTATE array, giving the weights for each 
state in forming the first and second order density matrix 
elements, which drive the orbital update methods.  Note 
that analytic gradients are available only when the WSTATE 
array is a unit vector, corresponding to a pure state, such 
as WSTATE(1)=0,1,0 which permits gradients of the first 
excited state to be computed.  When used for state averaged 
MCSCF, WSTATE is normalized to a unit sum, thus input of 
WSTATE(1)=1,1,1 really means a weight of 0.33333...  for 
the each of the states being averaged.

    ORMAS (Occupation Restricted Multiple Active Space) is 
a program designed to limit the size of the full CI 
problem, and may be useful when the number of active 
orbitals is 10 or higher.  By dividing your total active 
space into multiple subspaces, and specifying a range of 
electrons to occupy each subspace, most of the full CI's 
effect can be included.  ORMAS generates a full CI within 
each orbital subspace, taking the product of each small 
full CI to generate the determinant list.

    Here are some ideas on how to use ORMAS, which is a 
very flexible CI program:

a) single reference, arbitrary excition level CI-X, from
a closed shell reference:
       $det   ncore=y nact=z nels=10     (y+z=entire basis)
       $ormas nspace=2 mstart(1)=y+1,y+6 mine(1)=10-x,0
      This excites the 5 doubly occupied orbitals, to the
      desired excitation level of X.

      An open shell example of CI-SD from ...22111 might be
       $contrl mult=4
       $det    ncore=y nact=z nels=7     (y+z=entire basis)
       $ormas  nspace=3 mstart(1)=y+1,y+3,y+6
      No more than 2e- are allowed to be promoted from the
      doubly occupied or singly occupied spaces, and no
      more than 2 are allowed to enter the singly occupied
      or empty spaces.

   b) simple product of active spaces
      For example, consider furan, with two active
      subspaces.  Keeping the 5 true core and the 4 CH
      bonds in the core space, the sigma subspace might
      contains 5 ring sigma, one oxygen lone pair, and 5
      ring sigma antibonds, with a total of 12 e-.  The pi
      active space contains 5 pi orbitals and 6 e-:
       $det    ncore=9 nact=16 nels=18
       $ormas  nspace=2 mstart(1)=10,21 mine(1)=12,6
      Having the minimum and maximum electron counts the
      same is what makes this the simple product of two
      separate active spaces.  In other words, this is
      similar to the QCAS procedure of Nakano and Hirao,
      but ORMAS limits only the total electron counts, 
      not separately the numbers of alpha and beta e-,
      in other words all spin couplings are used.

   c) flexible occupancy between active subspaces
      Imagine that you are interested in excited states of
      formaldehyde, some of which will have Rydberg
      character, dominated by single excitations into 
      diffuse orbitals.  H2CO's valence states arise from 3
      orbitals, the CO pi and pi* and one oxygen lone pair.
      Placing the O sp lone pair and three sigma bonds into
      the filled space, and centering diffuse s,p,d shells
      on the carbon:
       $det    ncore=6 nact=12 nels=4
       $ormas  nspace=2 mstart(1)=7,10 mine(1)=3,0
      This is a 4e-, 3 orbital n,pi,pi* space to describe
      valence states, and excites one electron into the 9
      diffuse orbitals to describe Rydberg states.  It is
      many fewer determinants than a 4e- in 12 orbital FCI.

   d) RAS-like CI
      The previous example is reminiscent of Roos' RAS-SCF.
      In fact ORMAS can do RAS-SCF, which is three spaces:
      the lowest space is allowed to excite only a few
      electrons, a middle space that is the rest, and a top
      space into which only a few electrons can be excited.
      Suppose there are 10 e-, 10 orbitals, that the bottom
      and top spaces involve 3 orbitals, and that a "few"
      means specifically 2 e-:
       $det    ncore=20 nact=10 nels=10 $end
       $ormas  nspace=3 mstart(1)=21,24,28 mine(1)=4,2,0
      However, ORMAS can use more than 3 orbital subspaces.

   e) first or second order CI.
      Consider C2H4, with a 4 orbital active space of CC
      sigma, pi, pi*, and sigma*.  In order to correlate
      the four valence CH orbitals by double excitations,
      an MCSCF based on $DET, followed by SOCI based on
      $CIDET and $ORMAS, is:
       $contrl scftyp=mcscf cityp=ormas
       $mcscf  cistep=aldet
       $det    ncore=6 nact=4 nels=4
       $cidet  ncore=2 nact=y nels=12  (y=rest of basis)
       $ormas  nspace=3 mstart(1)=3,7,11 mine(1)=6,2,0
      which permits singles and doubles out of the CH and
      CC spaces, into the CC and external spaces.

    ORMAS is a full CI (or several full CI's) within each 
orbital subspace.  However, ORMAS does not generate all 
excitation levels between spaces (just those implied by the 
minimum and maximum electron counts you give).  This means 
ORMAS MCSCF runs must optimize active-active rotations 
between the subspaces, and therefore you should expect 
better convergence from FULLNR than SOSCF.

    ORMAS is sure to require orbital reordering.  For the 
furan example just mentioned, there is no reason to expect 
that the RHF occupied orbitals will not have the filled 
sigma and pi orbitals intermingled.  You must use the 
NORDER and IORDER keywords in $GUESS to carefully partition 
starting orbitals into sigma and pi subspaces.

    The selected (general) determinant list is used if 
CISTEP=GENCI, and the list is controlled by two input 
groups.  The first is $GEN, which is identical to $DET 
except for inclusion of an additional keyword GLIST=INPUT. 
This reads the determinants (as space products) from an 
additional input group $GCILST.  Completely arbitrary 
choices for the space products may be made, but peculiar 
lists may lead to poor MCSCF convergence.  The FOCAS 
converger should not be used, as it assumes full CI spaces.

    If you are doing straight CI calculations, the required 
input for each determinant CITYP is:
      ALDET needs $CIDET
      ORMAS needs $CIDET and $ORMAS
      GENCI needs $CIDET and $CIGEN and probably $GCILST
      FSOCI needs $CIDET and $SODET
In other words, $CIDET replaces $DET, and $CIGEN replaces  
$GEN, but the keywords in the groups mean the same thing. 
The reason for different names is to allow CI calculations 
to follow MCSCF in the same run, without clashing input 
group names.


    The GUGA-based CSF package was originally a set of 
different programs, so its input is spread over several 
input groups.  The CSFs are specified by a $CIDRT group in 
the case of CITYP=GUGA, and by a $DRT group for MCSCF 
wavefunctions.  Thus it is possible to perform an MCSCF 
defined by a $DRT input (or perhaps using $DET during the 
MCSCF), and follow this with a second order CI defined by a 
$CIDRT group, in the same run.

    The remaining input groups used by the GUGA CSFs are 
$CISORT, $GUGEM, $GUGDIA, and $GUGDM2 for MCSCF runs, with 
the latter two being the most important, and in the case of 
CI computations, $GUGDM and possibly $LAGRAN groups are 
relevant.  Perhaps the most interesting variables outside 
the $DRT/$CIDRT group are NSTATE in $GUGDIA to include 
excited states in the CI computation, IROOT in $GUGDM to 
select the CI state for properties, and WSTATE in $GUGDM2 
to control which state's orbitals are optimized, and 
possible state-averaging.
    The $DRT and $CIDRT groups are almost the same, with 
the only difference being orbitals restricted to double 
occupancy are called MCC in $DRT, and FZC in $CIDET. 
Therefore the rest of this section refers only to "$DRT".

    The CSFs are specified by giving a reference CSF, 
together with a maximum degree of electron excitation from 
that single CSF.  The MOs in the reference CSF are filled 
in the order MCC or FZC first, followed by DOC, AOS, BOS, 
ALP, VAL, and EXT (the Aufbau principle).  AOS, BOS, and 
ALP are singly occupied MOs.  ALP means a high spin alpha 
coupling, while AOS/BOS are an alpha/beta coupling to an 
open shell singlet.  This requires the value NAOS=NBOS, and 
their MOs alternate.  An example is
which gives the reference CSF
This is a doublet state with five unpaired electrons.  VAL 
orbitals are unoccupied only in the reference CSF, they 
will become occupied as the other CSFs are generated.  This 
is done by giving an excitation level, either explicitly by 
the IEXCIT variable, or implicitly by the FORS, FOCI, or 
SOCI flags.  One of these four keywords must be chosen, and 
during MCSCF runs, this is usually FORS.

    Consider another simpler example, for an MCSCF run,
      NMCC=3 NDOC=3 NVAL=2
which gives the reference CSF
having six electrons in five active orbitals.  MCSCF 
calculations are usually of the Full Optimized Reaction 
Space (FORS) type.  Some workers refer to FORS as CASSCF, 
complete active space SCF.  These are the same, but the 
keyword is spelled FORS in GAMESS.  In the present 
instance, choosing FORS=.TRUE. gives an excitation level of 
4, as the 6 valence electrons have only 4 holes available 
for excitation.  MCSCF runs typically have only a small 
number of VAL orbitals.  It is common to summarize this 
example as "six electrons in five orbitals".

    The next example is a first or second order multi-
reference CI wavefunction, where
      NFZC=3 NDOC=3 NVAL=2 NEXT=-1
leads to the reference CSF
FOCI or SOCI is chosen by selecting the appropriate flag, 
the correct excitation level is automatically generated. 
Note that the -1 for NEXT causes all remaining MOs to be 
included in the external orbital space.  One way of viewing 
FOCI and SOCI wavefunctions is as all singles, or all 
singles and doubles, from the entire MCSCF wavefunction as 
a reference.  An equivalent way of saying this is that all 
CSFs with N electrons (in this case N=6) distributed in the 
valence orbitals in all ways (that is the FORS MCSCF 
wavefunction) make up the reference wavefunction. To this, 
FOCI adds all CSFs with N-1 electrons in active and 1 
electron in external orbitals.  SOCI adds all CSFs with N-2 
electrons in active orbitals and 2 in external orbitals.  
SOCI is often prohibitively large, but is also a very 
accurate wavefunction.  SOCI can also be performed with 
determinants, as CITYP=FSOCI, or CITYP=ORMAS.  The latter 
may be the most efficient way to generate SOCI energies.  
For larger molecules, where SOCI is impractical, the most 
effective way to recover dynamic correlation energy is the 
multireference perturbation method.

    Sometimes people use the CI package for ordinary single 
reference CI calculations, such as
        NFZC=3 NDOC=5 NVAL=34
which means the reference RHF wavefunction is
and in this case NVAL is a large number conveying the total 
number of -virtual- orbitals into which electrons are 
excited.  The excitation level would be given as IEXCIT=2, 
perhaps, to perform a SD-CI.  All excitations smaller than 
the value of IEXCIT are automatically included in the CI.  
Note that NVAL's spelling was chosen to make the most sense 
for MCSCF calculations, and so it is a bit of a misnomer 

     Before going on, there is a quirk related to single 
reference CI that should be mentioned.  Whenever the single 
reference contains unpaired electrons, such as
       NFZC=3 NDOC=4 NALP=2 NVAL=33
some "extra" CSFs will be generated.  The reference here 
can be abbreviated 
    2222 11 000 000 000 000 000 000 000 000 000 000 000
Supposing IEXCIT=2, the following CSF
    2200 22 000 011 000 000 000 000 000 000 000 000 000
will be generated and used in the CI.  Most people would 
prefer to think of this as a quadruple excitation from the 
reference, but acting solely on the reasoning that no more 
than two electrons went into previously vacant NVAL 
orbitals, the GUGA CSF package decides it is a double. So, 
an open shell SD-CI calculation with GAMESS will not give 
the same result as other programs, although the result for 
any such calculation with these "extras" is correctly 
computed.  Note that if you also select the INTACT option, 
the extra space products are eliminated, but that some of 
the spin couplings for the truly IEXCIT'd space products 
are also eliminated.  Note that this kind of problem does 
not arise if you use ORMAS!

    As was discussed above, the CSFs are automatically 
spin-symmetry adapted, with S implicit in the reference 
CSF.  The spin quantum number you appear to be requesting 
in $DRT (basically, S = NALP/2) will be checked against the 
value of MULT in $CONTRL.  The total number of electrons, 
2*NMCC(or NFZC) + 2*NDOC + NAOS + NBOS + NALP will be 
checked against the input given for ICHARG.

    The CSF package is also able to exploit spatial 
symmetry, which like the spin and charge, is implicitly 
determined by the choice of the reference CSF.  The keyword 
GROUP in $DRT governs the use of spatial symmetry.
    The CSF program works with Abelian point groups, which 
are D2h and any of its subgroups.  However, $DRT allows the 
input of some (but not all) higher point groups.  For non-
Abelian groups, the program automatically assigns the 
orbitals to an irrep in the highest possible Abelian 
subgroup.  For the other non-Abelian groups, you must at 
present select GROUP=C1.  Note that when you are computing 
a Hessian matrix, many of the displaced geometries are 
asymmetric, hence you must choose C1 in $DRT (however, be 
sure to use the highest symmetry possible in $DATA!).

    The symmetry of the reference CSF given in your $DRT is 
one way to determine the symmetry of the CSFs which are 
generated.  As an example, consider a molecule with Cs 
symmetry, and these two reference CSFs
      ...MCC...DOC DOC VAL VAL
      ...MCC...DOC AOS BOS VAL
Suppose that the 2nd and 3rd active MOs have symmetries a' 
and a".  Both of these generate singlet wavefunctions, with 
4 electrons in 4 active orbitals, but the former constructs 
1-A' CSFs, while the latter generates 1-A" CSFs.  However, 
if the 2nd and 3rd orbitals have the same symmetry type, an 
identical list of CSFs is generated.  The alternative is to 
enter the spatial symmetry with the STSYM keyword.

    In cases with high point group symmetry, it may be 
possible to generate correct state degeneracies only by 
using no symmetry (GROUP=C1) when generating CSFs.  As an 
example, consider the 2-pi ground state of NO.  If you use 
GROUP=C4V, which will be mapped into its highest Abelian 
subgroup C2v, the two components of the pi state will be 
seen as belonging to different irreps, B1 and B2. The only 
way to ensure that both sets of CSFs are generated is to 
enforce no symmetry at all, so that CSFs for both 
components of the pi level are generated.  This permits 
state averaging (WSTATE(1)=0.5,0.5) to preserve cylindrical 
symmetry.  It is however perfectly feasible to use C4v or 
D4h symmetry in $DRT when treating sigma states.

     The use of spatial symmetry decreases the number of 
CSFs, and thus the size of the Hamiltonian that must be 
computed.  In molecules with high symmetry, this may lead 
to faster run times with the GUGA CSF code, compared to the 
determinant code.

starting orbitals
    The first step is to partition the orbital space into 
core, active, and external sets, in a manner which is 
sensible for your chemical problem.  This is a bit of an 
art, and the user is referred to the references quoted at 
the end of this section.  Having decided what MCSCF to 
perform, you now must consider the more pedantic problem of 
what orbitals to begin the MCSCF calculation with.
    You should always start an MCSCF run with orbitals from 
some other run, by means of GUESS=MOREAD.  Do not expect to 
be able to use HUCKEL!  At the start of a MCSCF problem, 
use orbitals from some appropriate converged SCF run.  A 
realistic example of an MCSCF calculation is GAMESS 
examples 8 and 9.  Once you get an MCSCF to converge, you 
can and should use these MCSCF MOs at other nearby 
geometries (MOREAD will apply an appropriate Schmidt 

    Starting from SCF orbitals can take a little bit of 
care.  Most of the time (but not always) the orbitals you 
want to correlate will be the highest occupied orbitals in 
the SCF.  Fairly often, however, the correlating orbitals 
you wish to use will not be the lowest unoccupied virtuals 
of the SCF.  You will soon become familiar with NORDER=1 in 
$GUESS, as reordering is needed in 50% or more cases.

    The occupied and especially the virtual canonical SCF 
MOs are often spread out over regions of the molecule other 
than "where the action is".  Orbitals which remedy this can 
generated by two additional options at almost no CPU cost.

    The best way to improve upon the SCF canonical virtual 
orbitals as starting MOs is to generate valence virtual 
orbitals (VVOs), after any RHF, ROHF, or GVB calculation.  
These are constructed by projection of internally stored 
atomic core and valence orbitals onto the SCF external 
orbitals, so by construction, the resulting VVOs are 
valence in character.  See VVOS in $SCF.

   An alternative choice, usually not as good as VVOs, are 
the modified virtual orbitals.  MVOs are obtained by 
diagonalizing the Fock operator of a very positive ion, 
within the virtual orbital space only. As implemented in 
GAMESS, MVOs can be obtained at the end of any RHF, ROHF, 
or GVB run by setting MVOQ in $SCF nonzero, at the cost of 
a single SCF cycle.  Typically, we use MVOQ=+6.  Generating 
MVOs does not change any of the occupied SCF orbitals of 
the original neutral, but gives more valence-like LUMOs.

    Another way to improve SCF starting orbitals is by a 
partial localization of the occupied orbitals.  Typically 
MCSCF active orbitals are concentrated in the part of the 
molecule where bonds are breaking, etc.  Canonical SCF MOs 
are normally more spread out.  By choosing LOCAL=BOYS along 
with SYMLOC=.TRUE. in $LOCAL, you can get orbitals which 
are localized, but still retain orbital symmetry to help 
speed the MCSCF along.  In groups with an inversion center, 
a SYMLOC Boys localization does not change the orbitals, 
but you can instead use LOCAL=POP.  LOCAL=RUEDNBRG may also 
be used, but requires more machine resources, and the other 
two localizations are normally good enough for starting 
orbital purposes.  Localization tends to order the orbitals 
fairly randomly, so be prepared to inspect them, and then 
reorder them appropriately.

    In case VVOS are generated in the same run as the 
localization, the localization is also applied within the 
valence virtual space.  The effect is to localize occupied 
orbitals into lone pairs and bond pairs, and in the VVOs 
space, to find localized antibond pairs.  When you choose 
the bonds and antibonds of chemical interest as the 
starting orbitals for the active space, convergence to the 
desired MCSCF wavefunction should follow.

    If you take the time to design your active space 
sensibly, select appropriate starting orbitals from the 
occupied and VVO unoccupied spaces, possibly by localizing 
these two subspaces, and carefully inspect your converged 
results, you will be able to carry out MCSCF computations 
    Convergence of MCSCF is by no means guaranteed.  Poor 
convergence can invariably be traced back to either a poor 
initial selection of orbitals, or poor design of the active 
space.  The best advice is, before you even start: 
    "Look at the orbitals."
    "Then, look at the orbitals again".
Later, if you have any trouble:
    "Look at the orbitals some more".
Few people are able to see the orbital shapes in the LCAO 
matrix in a log file, and so need a visualization program.  
In particular, you should download a copy of MacMolPlt from 
This runs on all popular desktop operating systems (Apple, 
Linux, and Windows), making it easy to see your final MCSCF 
orbital shapes.

    Even if you don't have any trouble, look at the 
orbitals to see if they converged to what you expected, and 
have reasonable occupation numbers.  It is particularly 
useful to check the oriented localized MCSCF orbitals (see 
the discussion of this in the section on localized orbitals 
in this section for more information).  MCSCF is by no 
means the sort of "black box" that RHF is these days, so 
please look very carefully at your final results.

miscellaneous hints
    It is very helpful to execute a EXETYP=CHECK run before 
doing any MCSCF or CI run.  The CHECK run will tell you the 
total number of configurations and check the charge and 
multiplicity and electronic state symmetry, based on your 
input.  The CHECK run also lets the program feel out the 
memory that will be required to actually do the run. Thus 
the CHECK run can potentially prevent costly mistakes, or 
tell you when a calculation is prohibitively large.

    A very common MCSCF wavefunction has 2 electrons in 2 
active MOs.  This is the simplest possible wavefunction 
describing a singlet diradical.  While this function can be 
obtained in an MCSCF run (using NACT=2 NELS=2 or NDOC=1 
NVAL=1), it can be obtained much faster by use of the GVB 
code, with one GVB pair.  This GVB-PP(1) wavefunction is 
also known in the literature as two configuration SCF, or 
TCSCF.  The two configurations of this GVB are equivalent 
to the three configurations used in this MCSCF, as orbital 
optimization in natural form (configurations 20 and 02) 
causes the coefficient of the 11 configuration to vanish.

    If you are using a large active space (say, 12 or more 
orbitals), the main bottleneck in the MCSCF calculation is 
the formation and diagonalization of the Hamiltonian, not 
the integral transformation and orbital updates.  Of 
course, since determinants are much faster than CSFs, and 
do not use large disk files, you should use determinants 
for large active spaces.  In this case, you would be wise 
to switch to FULLNR, which will minimize the total number 
of iterations, and thus the number of CI calculations.  
Note that by selecting ITERMX=5 in $DET or $GEN, you can 
avoid fully converging the CI during each MCSCF iteration, 
saving a bit of time.  Since each iteration's CI 
calculation starts with the previous iteration's result, 
the CI vectors will become fully converged during the MCSCF 
cycles.  The total run time may decrease, although a few 
additional MCSCF iterations may be required.  For small 
active spaces, where the CI step takes trivial time, you 
should use a bigger ITERMX to ensure fully converged CI 
states are generated every iteration.
    If you choose to use ORMAS, a general determinant CI, 
or if you select an CSF excitation level IEXCIT smaller 
than that needed to generate the FORS space, you must use 
the SOSCF, JACOBI, or FULLNR method as these can optimize 
active-active rotations.  Be sure to set FORS=.FALSE. in 
$MCSCF when for non-full CI cases, or else very poor 
convergence will result.  Actually, the convergence for 
incomplete active spaces is likely to be poorer than for 
full active spaces, anyway.

    A good way to check the active space is to localize the 
orbitals, to see if they resemble the atomic orbitals which 
you imagined formed the bonds, antibonds, and lone pairs in 
the active space.  The ORIENT keyword in $LOCAL will print 
a density matrix analysis, showing active electron bonding 
and antibonding patterns (see reference 18/19 below).

                       - - - - -

    The MCSCF technology in GAMESS is the result of some 
considerable programming effort:  The FOCAS, serial FULLNR, 
and QUAD convergers were adapted from Michel Dupuis' HONDO 
program.  The SOSCF converger was written by Galina Chaban, 
the parallel FULLNR converger is due to Graham Fletcher, 
and the JACOBI converger is due to Joe Ivanic.  The GUGA CI 
programs were written by Bernie Brooks and others, while 
all determinant CI codes (ALDET, GENCI, ORMAS, and FSOCI) 
stem from Joe Ivanic.  Analytic nuclear Hessians were 
programmed by Tim Dudley.  The CSF-based multireference 
pertubation program was written by Haruyuki Nakano, with a 
determinant implementation provided by Joe Ivanic.  Shiro 
Koseki and Dmitri Fedorov are responsible for the spin-
orbit coupling and transition moment codes.  The expertise 
of Klaus Ruedenberg in MCSCF wavefunctions has been the 
inspiration for many of these developments!
MCSCF references
    There are several review articles about MCSCF listed 
below.  Of these, the first two are a nice overview of the 
subject, the final 3 are more technical.
  1.  "The Construction and Interpretation of MCSCF
      M.W.Schmidt and M.S.Gordon,  
         Ann.Rev.Phys.Chem. 49,233-266(1998)
 2a. "The Multiconfiguration SCF Method"
      B.O.Roos, in "Methods in Computational Molecular
        Physics", edited by G.H.F.Diercksen and S.Wilson
        D.Reidel Publishing, Dordrecht, Netherlands,
        1983, pp 161-187.
 2b. "The Multiconfiguration SCF Method"
      B.O.Roos, in "Lecture Notes in Quantum Chemistry",
        edited by B.O.Roos, Lecture Notes in Chemistry v58, 
        Springer-Verlag, Berlin, 1994, pp 177-254.
  3. "Optimization and Characterization of a MCSCF State"
     J.Olsen, D.L.Yeager, P.Jorgensen
        Adv.Chem.Phys. 54, 1-176(1983).
  4. "Matrix Formulated Direct MCSCF and Multiconfiguration
       Reference CI Methods"
     H.-J.Werner,  Adv.Chem.Phys.  69, 1-62(1987).
  5. "The MCSCF Method"
     R.Shepard,  Adv.Chem.Phys.  69, 63-200(1987).
    There is an entire section on the choice of active 
spaces in Reference 1.  As this is a matter of great 
importance, here are two alternate presentations of the 
design of active spaces:
  6. "The CASSCF Method and its Application in Electronic
       Structure Calculations"
     B.O.Roos, in "Advances in Chemical Physics", vol.69,
        edited by K.P.Lawley, Wiley Interscience, New York,
        1987, pp 339-445.
  7. "Are Atoms Intrinsic to Molecular Electronic
     K.Ruedenberg, M.W.Schmidt, M.M.Gilbert, S.T.Elbert
       Chem.Phys. 71, 41-49, 51-64, 65-78 (1982).
    Two papers germane to the FOCAS implementation are
  8. "An Efficient first-order CASSCF method based on
        the renormalized Fock-operator technique."
     U.Meier, V.Staemmler  Theor.Chim.Acta 76, 95-111(1989)
  9. "Modern tools for including electron correlation in
        electronic structure studies"
     M.Dupuis, S.Chen, A.Marquez, in "Relativistic and
        Electron Correlation Effects in Molecules and
        Solids", edited by G.L.Malli, Plenum, NY 1994

    The paper germane to the the SOSCF converger is

 10. "Approximate second order method for orbital
      optimization of SCF and MCSCF wavefunctions"
     G.Chaban, M.W.Schmidt, M.S.Gordon
     Theor.Chem.Acc. 97: 88-95(1997)
    Two papers germane to the FULLNR converger, and two 
discussing implementation details are
 11. "General second order MCSCF theory: A Density Matrix
        Directed Algorithm"
     B.H.Lengsfield, III, J.Chem.Phys. 73,382-390(1980).
 12. "The use of the Augmented Matrix in MCSCF Theory"
     D.R.Yarkony, Chem.Phys.Lett. 77,634-635(1981).
 13. M.Dupuis, P.Mougenot, J.D.Watts, in "Modern Techniques
     in Theoretical Chemistry", E.Clementi, editor, ESCOM,
     Leiden, 1989, chapter 7.
 14. "A parallel multi-configuration self-consistent field
     G.D.Fletcher, Mol.Phys. 105, 2971-2976(2007)

    The paper describing the JACOBI converger is

 15. "A MCSCF method for ground and excited states based on
      full optimizatons of successive Jacobi rotations"
     J.Ivanic, K.Ruedenberg
     J.Comput.Chem. 24, 1250-1262(2003)

    For determinant CI codes, see

 16. "Identification of deadwood in configuration spaces
      through general direct configuration interaction"
     J.Ivanic, K.Ruedenberg
     Theoret.Chem.Acc. 106, 339-351(2001)
 17. "Direct configuration interaction and multi-
      configurational self-consistent-field method for
      multiple active spaces with variable occupancies.
      Part I.Method  Part II.Applications"
     J.Chem.Phys.  119, 9364-9376 and 9377-9385(2003)

    For CSFs, see

 18. "GUGA approach to the electron correlation problem"
     B.R.Brooks, H.F.Schaefer
       J.Chem.Phys.  70, 5092-5106(1979)

    Orientation of localized MCSCF active orbitals for 
bonding analysis:

 19. J.Ivanic, G.M.Atchity, K.Ruedenberg
       Theoret.Chem.Acc. 120, 281-294(2008)
 20. J.Ivanic, K.Ruedenberg
       Theoret.Chem.Acc. 120, 295-305(2008)

created on 7/7/2017