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
wavefunctions"
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
on ENGTOL.
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
CSFs
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
CSFs.
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
$DET STSYM=B1 NSTATE=3 NCORE=xx NELS=8 NACT=6 $END
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
maxe(1)=10,x
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
mine(1)=2,1,0
maxe(1)=4,5,2
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
maxe(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
maxe(1)=4,1
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
maxe(1)=6,6,2
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
maxe(1)=8,6,2
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.
CSF CI
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
NFZC=1 NDOC=2 NAOS=2 NBOS=2 NALP=1 NVAL=3
which gives the reference CSF
FZC,DOC,DOC,AOS,BOS,AOS,BOS,ALP,VAL,VAL,VAL
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
MCC,MCC,MCC,DOC,DOC,DOC,VAL,VAL
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
FZC,FZC,FZC,DOC,DOC,DOC,VAL,VAL,EXT,EXT,...
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
FZC FZC FZC DOC DOC DOC VAL VAL ... VAL
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
here.
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
orthogonalization).
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
correctly.
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
http://code.google.com/p/wxmacmolplt
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
wavefunctions"
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
Wavefunctions?"
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
algorithm"
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.Ivanic
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