How to do RHF, ROHF, UHF, and GVB calculations

general considerations
    These four SCF wavefunctions are all based on Fock 
operator techniques, even though some GVB runs use more 
than one determinant.  Thus all of these have an intrinsic 
N**4 time dependence, because they are all driven by 
integrals in the AO basis.  This similarity makes it 
convenient to discuss them all together.  In this section 
we will use the term HF to refer generically to any of 
these four wavefunctions, including the multi-determinate 
GVB-PP functions.  $SCF is the main input group for all 
these HF wavefunctions.

    As will be discussed below, in GAMESS the term ROHF 
refers to high spin open shell SCF only, but other open 
shell coupling cases are possible using the GVB code.
    Analytic gradients are implemented for every possible 
HF type calculation possible in GAMESS, and therefore 
numerical hessians are available for each.
    Analytic hessian calculation is implemented for RHF, 
ROHF, and any GVB case with NPAIR=0 or NPAIR=1.  Analytic 
hessians are more accurate, and much more quickly computed 
than numerical hessians, but require additional disk 
storage to perform an integral transformation, and also 
more physical memory.
    The second order Moller-Plesset energy correction (MP2) 
is implemented for RHF, UHF, ROHF, and MCSCF wavefunctions.  
Analytic gradients may be obtained for MP2 with RHF, UHF, 
or ROHF reference wavefunctions, and MP2 level properties 
are therefore available for these, see MP2PRP in $MP2.  All 
other cases give properties for the SCF function.
    Direct SCF is implemented for every possible HF type 
calculation.  The direct SCF method may not be used with 
DEM convergence.  Direct SCF may be used during energy, 
gradient, numerical or analytic hessian, CI or MP2 energy 
correction, or localized orbital computations.
direct SCF
    Normally, HF calculations proceed by evaluating a large 
number of two electron repulsion integrals, and storing 
these on a disk.  This integral file is read in once during 
each HF iteration to form the appropriate Fock operators.  
In a direct HF, the integrals are not stored on disk, but 
are instead reevaluated during each HF iteration.  Since 
the direct approach *always* requires more CPU time, the 
default for DIRSCF in $SCF is .FALSE.
    Even though direct SCF is slower, there are at least 
two reasons why you may want to consider using it.  The 
first is that it may not be possible to store all of the 
integrals on the disk drives attached to your computer. 
Second, what you are really interested in is reducing the 
wall clock time to obtain your answer, not the CPU time. 
Workstations, particularly nodes with multiple CPUs and 
only one disk subsystem, may have modest hardware I/O 
capabilities.  Other environments such as a mainframe 
shared by many users may also have very poor CPU/wall clock 
performance for I/O bound jobs such as conventional HF.
    You can estimate the disk storage requirements for 
conventional HF using a P or PK file by the following 
          nint = 1/sigma * 1/8 * N**4
          Mbytes = nint * x / 1024**2
Here N is the total number of basis functions in your run, 
which you can learn from an EXETYP=CHECK run.  The 1/8 
accounts for permutational symmetry within the integrals.  
Sigma accounts for the point group symmetry, and is 
difficult to estimate accurately.  Sigma cannot be smaller 
than 1, in no symmetry (C1) calculations.  For benzene, 
sigma would be almost six, since you generate 6 C's and 6 
H's by entering only 1 of each in $DATA.  For water sigma 
is not much larger than one, since most of the basis set is 
on the unique oxygen, and the C2v symmetry applies only to 
the H atoms.  The factor x is 12 bytes per integral for 
basis sets smaller than 255, and 16 otherwise. Finally, 
since integrals that are very close to zero need not be 
stored on disk, the actual power dependence is not as bad 
as N**4, and in fact in the limit of very large molecules 
can be as low as N**2.  Thus plugging in sigma=1 should 
give you an upper bound to the actual disk space needed.  
If the estimate exceeds your available disk storage, your 
only recourse is direct HF.
    What are the economics of direct HF?  Naively, if we 
assume the run takes 10 iterations to converge, we must 
spend 10 times more CPU time computing the integrals on 
each iteration.  However, we do not have to waste any CPU 
time reading blocks of integrals from disk, or in unpacking 
their indices.  We also do not have to waste any wall clock 
time waiting for a relatively slow mechanical device such 
as a disk to give us our data.
    There are some less obvious savings too, as first noted 
by Almlof.  First, since the density matrix is known while 
we are computing integrals, we can use the Schwarz 
inequality to avoid doing some of the integrals.  In a 
conventional SCF this inequality is used to avoid doing 
small integrals.  In a direct SCF it can be used to avoid 
doing integrals whose contribution to the Fock matrix is 
small (density times integral=small).  Secondly, we can 
form the Fock matrix by calculating only its change since 
the previous iteration.  The contributions to the change in 
the Fock matrix are equal to the change in the density 
times the integrals.  Since the change in the density goes 
to zero as the run converges, we can use the Schwarz 
screening to avoid more and more integrals as the 
calculation progresses.  The input option FDIFF in $SCF 
selects formation of the Fock operator by computing only 
its change from iteration to iteration.  The FDIFF option 
is not implemented for GVB since there are too many density 
matrices from the previous iteration to store, but is the 
default for direct RHF, ROHF, and UHF.
    So, in our hypothetical 10 iteration case, we do not 
spend as much as 10 times more time in integral evaluation.  
Additionally, the run as a whole will not slow down by 
whatever factor the integral time is increased.  A direct 
run spends no additional time summing integrals into the 
Fock operators, and no additional time in the Fock 
diagonalizations.  So, generally speaking, a RHF run with 
10-15 iterations will slow down by a factor of 2-4 times 
when run in direct mode.  The energy gradient time is 
unchanged by direct HF, and this is a large time compared 
to HF energy, so geometry optimizations will be slowed down 
even less.  This is really the converse of Amdahl's law:  
if you slow down only one portion of a program by a large 
amount, the entire program slows down by a much smaller 
    To make this concrete, here are some times for GAMESS 
for a job which is a RHF energy for a SbC4O2NH4.  These 
timings were obtained an extremely long time ago, on a 
DECstation 3100 under Ultrix 3.1, which was running only 
these tests, so that the wall clock times are meaningful.  
This system is typical of Unix workstations in that it uses 
SCSI disks, and the operating system is not terribly good 
at disk I/O.  By default GAMESS stores the integrals on 
disk in the form of a P supermatrix, because this will save 
time later in the SCF cycles.  By choosing NOPK=1 in 
$INTGRL, an ordinary integral file can be used, which 
typically contains many fewer integrals, but takes more CPU 
time in the SCF.  Because the DECstation is not terribly 
good at I/O, the wall clock time for the ordinary integral 
file is actually less than when the supermatrix is used, 
even though the CPU time is longer.  The run takes 13 
iterations to converge, the times are in seconds.
                           P supermatrix   ordinary file
   # nonzero integrals      8,244,129       6,125,653
   # blocks skipped            55,841          55,841
   CPU time (ints)              709              636
   CPU time (SCF)              1289             1472
   CPU time (job total)        2123             2233
   wall time (job total)       3468             3200
When the same calculation is run in direct mode (integrals 
are processed like an ordinary integral disk file when 
running direct),
      iteration 1:         FDIFF=.TRUE.   FDIFF=.FALSE.
   # nonzero integrals       6,117,416      6,117,416
   # blocks skipped             60,208         60,208
      iteration 13:
   # nonzero integrals       3,709,733      6,122,912
   # blocks skipped            105,278         59,415
   CPU time (job total)         6719            7851
   wall time (job total)        6764            7886
    Note that elimination of the disk I/O dramatically 
increases the CPU/wall efficiency.  Here's the bottom line 
on direct HF:
      best direct CPU / best disk CPU = 6719/2123 = 3.2
      best direct wall/ best disk wall= 6764/3200 = 2.1
Direct SCF is slower than conventional disk SCF, but not 
outrageously so!  From the data in the tables, we can see 
that the best direct method spends about 6719-1472 = 5247 
seconds doing integrals.  This is an increase of about 
5247/636 = 8.2 in the time spent doing integrals, in a run 
that does 13 iterations (13 times evaluating integrals). 
8.2 is less than 13 because the run avoids all CPU charges 
related to I/O, and makes efficient use of the Schwarz 
inequality to avoid doing many of the integrals in its 
final iterations.

convergence accelerators

    Generally speaking, the simpler the HF function, the 
better its convergence.  In our experience, the majority of 
RHF and ROHF runs converge readily from GUESS=HUCKEL.  UHF 
often takes considerably more iterations than either of 
these, due to the extremely common occurrence of heavy spin 
contamination.  GVB runs typically require GUESS=MOREAD, 
although the Huckel guess usually works for NPAIR=0.  GVB 
cases with NPAIR greater than one are particularly 

    Unfortunately, not all HF runs converge readily.  The 
best way to improve your convergence is to provide better 
starting orbitals!  In many cases, this means to MOREAD 
orbitals from some simpler HF case.  For example, if you 
want to do a doublet ROHF, and the HUCKEL guess does not 
seem to converge, do this:  Do an RHF on the +1 cation. RHF 
is typically more stable than ROHF, UHF, or GVB, and 
cations are usually readily convergent.  Then MOREAD the 
cation's orbitals into the neutral calculation which you 
wanted to do at first.

    GUESS=HUCKEL does not always start with the correct 
electronic configuration.  It may be useful to use PRTMO in 
$GUESS during a CHECK run to examine the starting orbitals, 
and then reorder them with NORDER if that seems 

    Of course, by default GAMESS uses the convergence 
procedures which are usually most effective.  Still, there 
are cases which are difficult, so the $SCF group permits 
you to select several alternative methods for improving 
convergence.  Briefly, these are

    EXTRAP.  This extrapolates the three previous Fock 
matrices, in an attempt to jump ahead a bit faster.  This 
is the most powerful of the old-fashioned accelerators, and 
normally should be used at the beginning of any SCF run.  
When an extrapolation occurs, the counter at the left of 
the SCF printout is set to zero.

    DAMP.  This damps the oscillations between several 
successive Fock matrices.  It may help when the energy is 
seen to oscillate wildly.  Thinking about which orbitals 
should be occupied initially may be an even better way to 
avoid oscillatory behaviour.

    SHIFT.  Level shifting moves the diagonal elements of 
the virtual part of the Fock matrix up, in an attempt to 
uncouple the unoccupied orbitals from the occupied ones.  
At convergence, this has no effect on the orbitals, just 
their orbital energies, but will produce different (and 
hopefully better) orbitals during the iterations.

    RSTRCT.  This limits mixing of the occupied orbitals 
with the empty ones, especially the flipping of the HOMO 
and LUMO to produce undesired electronic configurations or 
states.  This should be used with caution, as it makes it 
very easy to converge on incorrectly occupied electronic 
configurations, especially if DIIS is also used.  If you 
use this, be sure to check your final orbital energies to 
see if they are sensible.  A lower energy for an unoccupied 
orbital than for one of the occupied ones is a sure sign of 

    DIIS.  Direct Inversion in the Iterative Subspace is a 
modern method, due to Pulay, using stored error and Fock 
matrices from a large number of previous iterations to 
interpolate an improved Fock matrix.  This method was 
developed to improve the convergence at the final stages of 
the SCF process, but turns out to be quite powerful at 
forcing convergence in the initial stages of SCF as well. 
By giving ETHRSH as 10.0 in $SCF, you can practically 
guarantee that DIIS will be in effect from the first 
iteration.  The default is set up to do a few iterations 
with conventional methods (extrapolation) before engaging 
DIIS.  This is because DIIS can sometimes converge to 
solutions of the SCF equations that do not have the lowest 
possible energy.  For example, the 3-A-2 small angle state 
of SiLi2 (see M.S.Gordon and M.W.Schmidt, Chem.Phys.Lett., 
132, 294-8(1986)) will readily converge with DIIS to a 
solution with a reasonable S**2, and an energy about 25 
milliHartree above the correct answer.  A SURE SIGN OF 
VALUE.  However, if you obtain orbitals at one point on a 
PES without DIIS, the subsequent use of DIIS with MOREAD 
will probably not introduce any problems.   Because DIIS is 
quite powerful, EXTRAP, DAMP, and SHIFT are all turned off 
once DIIS begins to work.  DEM and RSTRCT will still be in 
use, however.

    SOSCF.  Approximate second-order (quasi-Newton) SCF 
orbital optimization.  SOSCF will converge about as well as 
DIIS at the initial geometry, and slightly better at 
subsequent geometries.  There's a bit less work solving the 
SCF equations, too.   The method kicks in after the orbital 
gradient falls below SOGTOL.  Some systems, particularly 
transition metals with ECP basis sets, may have Huckel 
orbitals for which the gradient is much larger than SOGTOL. 
In this case it is probably better to use DIIS instead, 
with a large ETHRSH, rather than increasing SOGTOL, since 
you may well be outside the quadratic convergence region. 
SOSCF does not exhibit true second order convergence since 
it uses an approximation to the inverse hessian.  SOSCF 
will work for MOPAC runs, but is slower in this case. SOSCF 
will work for UHF, but its convergence may be better than 
DIIS.  SOSCF will work for non-Abelian cases, but may 
encounter problems if the open shell is degenerate.

    It should be clear that SOSCF and DIIS are the two 
work-horse convergers, with DAMP (and possibly SHIFT) 
useful in cases where the initial guess is such that these 
two are not engaged immediately.  If you compute many 
different types of molecules, you will find cases where 
SOSCF works but DIIS does not, but also cases where DIIS 
works but SOSCF does not (although often both will work).  
If you do not obtain convergence with one of these, try the 
other one!  If you still have problems, attempt to get 
better starting orbitals.

    DEM.  Direct energy minimization should be your last 
recourse.  It explores the "line" between the current 
orbitals and those generated by a conventional change in 
the orbitals, looking for the minimum energy on that line. 
DEM should always lower the energy on every iteration, but 
is very time consuming, since each of the points considered 
on the line search requires evaluation of a Fock operator.  
DEM will be skipped once the density change falls below 
DEMCUT, as the other methods should then be able to affect 
final convergence.   While DEM is working, RSTRCT is held 
to be true, regardless of the input choice for RSTRCT.  
Because of this, it behooves you to be sure that the 
initial guess is occupying the desired orbitals.  DEM is 
available only for RHF.  The implementation in GAMESS 
resembles that of R.Seeger and J.A.Pople, J.Chem.Phys. 65, 
265-271(1976).  Simultaneous use of DEM and DIIS resembles 
the ADEM-DIOS method of H.Sellers, Chem.Phys.Lett. 180, 
461-465(1991).  DEM does not work with direct SCF.

high spin open shell SCF (ROHF)

    Open shell SCF calculations are performed in GAMESS by 
both the ROHF code and the GVB code.  Note that when the 
GVB code is executed with no pairs, the run is NOT a true 
GVB run, and should be referred to in publications and 
discussion as a ROHF calculation.  Low spin couplings are 
possible with the GVB program.

    The ROHF module in GAMESS can handle any number of open 
shell electrons, provided these have a high spin coupling.  
For example: one open shell, doublet:
                $CONTRL SCFTYP=ROHF MULT=2 $END
             two open shells, triplet:
                $CONTRL SCFTYP=ROHF MULT=3 $END
             m open shells, high spin:
                $CONTRL SCFTYP=ROHF MULT=m+1 $END

    John Montgomery (who was then at United Technologies) 
is responsible for the ROHF implementation in GAMESS.  The 
following discussion is due to him, dating from 1988 when 
his method of forming a combined Fock operator was included 
in GAMESS.  Other choices (Euler and two "canonical" sets) 
were added to the table in 2009/2010.

    The Fock matrix in the MO basis has the form
                   closed       open        virtual
        closed      F2      |     Fb     | (Fa+Fb)/2
        open        Fb      |     F1     |    Fa
        virtual   (Fa+Fb)/2 |     Fa     |    F0
where Fa and Fb are the usual alpha and beta Fock matrices 
any UHF program produces.  All ROHF methods agree on these, 
as they are the variational conditions that separate the 
doubly occupied, alpha occupied, and empty orbital spaces.  
The diagonal blocks can be written
               F2 = Acc*Fa + Bcc*Fb
               F1 = Aoo*Fa + Boo*Fb
               F0 = Avv*Fa + Bvv*Fb
Some choices for the canonicalization coefficients to 
define the diagonal blocks are
                          Acc  Bcc  Aoo Boo  Avv  Bvv
 Guest and Saunders       1/2  1/2  1/2 1/2  1/2  1/2
 Roothaan single matrix  -1/2  3/2  1/2 1/2  3/2 -1/2
 Davidson/1988            1/2  1/2   1   0    1    0
 Binkley, Pople, Dobosh   1/2  1/2   1   0    0    1
 McWeeny and Diercksen    1/3  2/3  1/3 1/3  2/3  1/3
 Faegri and Manne         1/2  1/2   1   0   1/2  1/2
 GVB program/Euler        1/2  1/2  1/2  0   1/2  1/2
 "canonical 1"            0    1    1    0    1    0
 "canonical 2"     (2S+1)/2S -1/2S  0    1  -1/2S (2S+1)/2S
See below for how these last two rows connect to ionization 

    The 1988 GAMESS ROHF program using a now deleted 
Davidson-type ROHF produced final orbitals matching the 
line "1988" above.  This differs from the choices made in 
Davidson's own MELD program.  The MELD program itself has 
always done a "cleanup" of the virtual space, after 
convergence, using Avv=Bvv=1/2, producing orbitals which 
are the same as Faegri/Manne.  If MELD's MP2 option is 
chosen, the occupied space is also altered after 
convergence, using Aoo=Boo=1/2, which is the Guest/Saunders 
line above.  Thus the term "Davidson orbitals" is used here 
to refer to the behavior of the now-deleted 1988 ROHF code 
in GAMESS, which didn't have either type of final orbital 

    The choice of the diagonal blocks is arbitrary, as ROHF 
is converged when the off diagonal blocks go to zero.  The 
exact choice for these blocks can however have an effect on 
the convergence rate.  This choice also affects the MO 
coefficients, and orbital energies, as the different 
choices produce different canonical orbitals within the 
three subspaces.  All methods, however, will give identical 
total wavefunctions, and hence identical properties such as 
nuclear gradients and hessians.  Some of the perturbation 
theories for open shell cases are defined in terms of a 
particular canonicalization, if so, GAMESS automatically 
canonicalizes after convergence so the desired orbitals and 
energies are given to the perturbation theory codes.
    The default coupling case in GAMESS is the Roothaan 
single matrix set.  Note that pre-1988 versions of GAMESS 
produced "Davidson/1988" orbitals.  If you would like to 
fool around with any of these other canonicalizations, the 
Acc, Aoo, Avv and Bcc, Boo, Bvv parameters can be input as 
the first three elements of ALPHA and BETA in $SCF.  For 
example, the McWeeny/Diercksen canonicalization, in double 
precision, is obtained by
 $scf    couple=.true. alpha(1)=0.333333333333333333,
                                0.333333333333333333 $end

    Here is some idea of the range of eigenvalues that 
result from using the various canonicalization schemes in 
the table above.  The system is 6-31G nitrogen atom, all 
runs give E= -54.3820511123 (matching all 10 decimals):
                       1s      2s      2p     "3p"    "3s"
Roothaan           -15.5514 -0.5306 -0.1774 +0.7666 +0.8704
McWeeny,Diercksen  -15.6214 -0.8745 -0.1183 +0.8984 +0.9684
Davidson           -15.6355 -0.9432 -0.5657 +0.8457 +0.9292
Guest,Saunders     -15.6355 -0.9432 -0.1774 +0.9248 +0.9880
Binkley,Pople,Dob. -15.6355 -0.9432 -0.5657 +1.0039 +1.0469
Faegri,Manne       -15.6355 -0.9432 -0.5657 +0.9248 +0.9880
GVB (aka Euler)    -15.6355 -0.9432 -0.2828 +0.9248 +0.9880
"canonical 1"      -15.5933 -0.7370 -0.5657 +0.8457 +0.9292
"canonical 2"      -15.7065 -1.2863 +0.2109 +1.0567 +1.0861

    Since the ionization potential (IP) for a 2p electron 
in Nitrogen is 0.53 Hartree, it is clear that most of the 
orbital energies above do not approximately predict this 
IP.  Recent work by Boris Plakhutin and co-workers (see the 
3 papers in the ROHF references above) leads to two sets of 
orbitals and eigenvalues, for the prediction of both IP and 
electron affinities (EA), for various ionization events, 
starting from a high spin ROHF state:
        canonical 1 (to produce high spin final states)
   A1 = remove  beta e- from filled space,
   B1 = remove alpha e- from open shell,
   C1 = attach alpha e- to virtual space.
        canonical 2 (to produce  low spin final states)
   A2 = remove alpha e- from filled space, 
   B2 = attach  beta e- to filled space,
   C2 = attach  beta e- to virtual space.
Correct handling of spin requires the value of the spin S 
of the initial state in the 2nd canonicalization set.  Two 
different ROHF runs are necessary to get all six EA and IP 

    Additional discussion about ROHF orbital energies may 
be found in
   K.R.Glaesemann, M.W.Schmidt
   J.Phys.Chem.A 114, 8772-8777(2010)   [available on-line]
along with a demonstration of the non-uniqueness of ROHF-
based perturbation theories.

    High-spin ROHF results are automatically obtained if a 
UHF calculation is constrained so that the occupied beta 
orbitals lie entirely within the occupied alpha orbital 
space.  Such a constraint is called CUHF, and while it will 
appear to have different alpha/beta orbitals, and 
alpha/beta eigenvalues, its spin expectation value will be 
=2S+1, and its energy will be exactly that of ROHF.  
The CUHF method is given by
   G.E.Scuseria, T.Tsuchimochi
   J.Chem.Phys. 134, 064101/1-14(2011)
CUHF runs as a UHF calculation, solving two distinct SCF 
equations in the alpha and beta space, so its convergence 
behavior and its final eigenvalues will differ from that 
for any of the ROHF canonicalizations.  See the section on 
perturbation theory in this chapter regarding CUHF's 
perturbation theory, CUMP2.

other open shell SCF cases (GVB)
    Genuine GVB-PP "perfect pairing" runs will be discussed 
later in this section.  First, we will consider how to do 
open shell SCF with the GVB part of the program, with 
NPAIR=0.  This is an extremely powerful open shell SCF 
program, capable of handling low spin couplings and/or 
degenerate open shells, with the energy formula

   E = 2 sum F-i*h-ii + sum sum ALPHA-ij*J-ij + BETA-ij*Kij
          i              i   j

Here F-i (meaning subscript i) is a fractional occupancy, 
and is 1.0 for the filled shell that is usually present 
below the open shell(s).  A d**2 shell has F-i=0.2.  The 
h,J,K are the usual one electron, coulomb, and exchange 
operators.  For a few common cases, the F, ALPHA, BETA are 
internally stored in the program,

one open shell, doublet:
     $SCF    NCO=xx NSETO=1 NO(1)=1 $END
two open shells, triplet:
     $SCF    NCO=xx NSETO=2 NO(1)=1,1 $END
two open shells, singlet:
     $SCF    NCO=xx NSETO=2 NO(1)=1,1 $END
    Note that the first two cases duplicate runs which the 
high spin ROHF module can readily do.  The last case is 
also an ROHF calculation, albeit for a low-spin coupling.  
One should note that the open shell singlet case is a 
variational calculation IF AND ONLY IF the two singly 
degenerate orbitals have different space symmetry.

    A publication should refer to any run that has NPAIR=0 
as being ROHF type, rather than GVB, since there are no 
perfect pairs in use.  Note, however, that the GVB program 
can, if you wish, have both open shells and genuine valence 
bond pairs at the same time!  See the following section 
about the use of GVB's perfect pairs.  The use of open 
shells and pairs is illustrated in one of the GAMESS 
standard example inputs.

    If you would like to do any cases other than those 
shown above, including many cases with orbital degeneracy, 
you must derive the coupling coefficients ALPHA and BETA, 
and input them with the occupancies F in the $SCF group.  
Fortunately, many cases can be looked up!

    Mariusz Klobukowski of the University of Alberta has 
shown how to obtain coupling coefficients for the GVB open 
shell program for many such open shell states.  These can 
be obtained from Appendix A of the book "A General SCF 
Theory" by Ramon Carbo and Josep M. Riera, Springer-Verlag 
(1978).  The rule is
       (1)      F(i) = 1/2 * omega(i)
       (2)  ALPHA(i) =       alpha(i)
       (3)   BETA(i) =      - beta(i),
where omega, alpha, beta are symbols used in these Tables.

    The keyword NSETO gives the number of open shells, and 
keyword NO gives the degeneracy of each open shell.  Thus 
the 5-S state of carbon (from configuration s1,p3) would 
enter NSETO=2 and NO(1)=1,3.  NCO is an easy keyword: that 
is the number of filled orbitals.  The total number of 
electrons in an open shell GVB run is
        NE = 2*NCO + sum 2*F(i)*NO(i)
and this value will be checked against your ICHARG input.

                        - - - -

    Specific input for all terms found in the atomic p**N 
configurations follow.  Be sure to enter 14 digits, as 
these values are part of a double precision energy formula!

    Values for the excited terms in the p**N configurations 
were extracted from C.F.Jackels and E.R.Davidson Int. J. 
Quantum Chem. 8, 707-714(1974), which explains the concept 
of averaging over equivalent determinants to enforce a 
symmetric density matrix, which preserves radial symmetry 
in the atomic orbitals.  The ALPHA and BETA values can also 
be extracted from Roothaan's 1960 A and B coefficients by 
the prescription detailed below.

!   p**1   2-P state
 $SCF    NCO=xx   NSETO=1  NO=3   COUPLE=.TRUE.
      F(1)=  1.0  0.16666666666667
  ALPHA(1)=  2.0  0.33333333333333  0.00000000000000
   BETA(1)= -1.0 -0.16666666666667 -0.00000000000000  $END

!   p**2   3-P state
 $SCF    NCO=xx   NSETO=1  NO=3   COUPLE=.TRUE.
      F(1)=  1.0  0.33333333333333
  ALPHA(1)=  2.0  0.66666666666667  0.16666666666667
   BETA(1)= -1.0 -0.33333333333333 -0.16666666666667  $END

For the 1-D excited state, change the open-shell parameters 
to ALPHA(3)=0.1 and BETA(3)= +0.03333333333333
For the 1-S excited state, change the open-shell parameters 
to ALPHA(3)=0.0 and BETA(3)= +0.33333333333333

!   p**3   4-S state
which is equivalent to
 $SCF    NCO=xx   NSETO=1  NO=3   COUPLE=.TRUE.
      F(1)=  1.0  0.50000000000000
  ALPHA(1)=  2.0  1.00000000000000  0.50000000000000
   BETA(1)= -1.0 -0.50000000000000 -0.50000000000000  $END

For 2-D, use ALPHA(3)= 0.4              and BETA(3)= -0.2
for 2-P, use ALPHA(3)= 0.33333333333333 and BETA(3)=  0.0

!   p**4   3-P state
 $SCF    NCO=xx   NSETO=1  NO=3   COUPLE=.TRUE.
      F(1)=  1.0  0.66666666666667
  ALPHA(1)=  2.0  1.33333333333333  0.83333333333333
   BETA(1)= -1.0 -0.66666666666667 -0.50000000000000  $END

For 1-D, use ALPHA(3)= 0.76666666666667 and BETA(3)= -0.3
for 1-S, use ALPHA(3)= 0.66666666666667 and BETA(3)=  0.0

!   p**5   2-P state
 $SCF    NCO=xx   NSETO=1  NO=3   COUPLE=.TRUE.
      F(1)=  1.0  0.83333333333333
  ALPHA(1)=  2.0  1.66666666666667  1.33333333333333
   BETA(1)= -1.0 -0.83333333333333 -0.66666666666667  $END

                        - - - -

Coupling constants for the highest spin state(s) in d**N 
configurations are taken from "Handbook of Gaussian Basis 
Sets", R.Poirier, R.Kari, I.G.Csizmadia, Elsevier, 
Amsterdam, 1985.

!     d**1   2-D state
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.1
         ALPHA(1)= 2.0, 0.20, 0.00
          BETA(1)=-1.0,-0.10, 0.00  $END

!     d**2   average of 3-F and 3-P states
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.2
         ALPHA(1)= 2.0, 0.40, 0.05
          BETA(1)=-1.0,-0.20,-0.05  $END
Note: "average" means a degeneracy weighted combination of 
determinants with Sz=S: here, 2 electrons in 5 orbitals 
coupled as a triplet is ten determinants with Sz=1.  The 
energy from these parameters is E=[7xE(3-F)+3xE(3-P)]/10.

!     d**3   average of 4-F and 4-P states
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.3
         ALPHA(1)= 2.0, 0.60, 0.15
          BETA(1)=-1.0,-0.30,-0.15  $END

!     d**4   5-D state
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.4
         ALPHA(1)= 2.0, 0.80, 0.30
          BETA(1)=-1.0,-0.40,-0.30 $END

!     d**5   6-S state

!     d**6   5-D state
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.6
         ALPHA(1)= 2.0, 1.20, 0.70
          BETA(1)=-1.0,-0.60,-0.50 $END

!     d**7   average of 4-F and 4-P states
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.7
         ALPHA(1)= 2.0, 1.40, 0.95
          BETA(1)=-1.0,-0.70,-0.55  $END

!     d**8   average of 3-F and 3-P states
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.8
         ALPHA(1)= 2.0, 1.60, 1.25
          beta(1)=-1.0,-0.80,-0.65  $end

!     d**9   2-D state
 $SCF    NCO=xx NSETO=1 NO=5 COUPLE=.TRUE.  F(1)=1.0,0.9
         ALPHA(1)= 2.0, 1.80, 1.60
          BETA(1)=-1.0,-0.90,-0.80 $END

    Note that GAMESS can do a proper calculation on the 
ground term for the d**2, d**3, d**7, d**8 configurations 
only by means of state averaged MCSCF.  For d**8, use
 $contrl scftyp=mcscf mult=3 $end
 $det    group=c1 ncore=xx nact=5 nels=8
         nstate=10 wstate(1)=1,1,1,1,1,1,1,0,0,0 $end
to correctly average the 7 lowest roots (3-F) with no 
weight given to the highest three roots (3-P).  Although 
this is done with the SCFTYP=MCSCF program, this is still a 
SCF calculation since only the orbitals, but not any CI 
coefficients, are optimized.

    Dirk Andrae in Berlin has provided a great many 
examples of using the MCSCF program to do terms found in 
the p**n, d**n, and even f**n configurations:
Type that web reference on just one line, of course.

    Open shell cases such as s**1,d**n are probably most 
easily tackled with the state-averaged MCSCF program.  The 
ORMAS CI code may be convenient in fixing the number of 
electrons found in each open shell.

    If you are a true afficionado of atomic calculations, 
including open shell f configurations, look for Russ 
Pitzer's version of the famous Roothaan/Bagus ATOM-SCF 
    R.M.Pitzer, Comput.Phys.Commun. 170, 239-264(2005)

                         - - - -

    The literature contains an alternate energy formula, 
involving the so-called Roothaan A and B parameters:

  E = 2 sum f-i*h-ii + sum sum f-i*f-j[2Aij*Jij - Bij*Kij]
         i              i   j

Comparing this to the GVB program's energy formula above, 
we immediately see that its fractional occupancy f is the 
same as GVB's F, and that
   ALPHA-ij = +2 * f-i * f-j * A-ij
    BETA-ij = -1 * f-i * f-j * B-ij
Let c=the closed shell orbitals, and o=a single open shell.  
Tables such as Roothaan's 1960 paper give only the A-oo and 
B-oo values, since A-cc = B-cc = A-co = B-co = 1.0.  It is 
easy to convert such A,B values to all ALPHA-cc,-co,-oo and 
BETA-cc,-co,-oo values.  Degenerate open shells in linear 
molecules follow immediately from Roothaan's 1960 table:

!       2-Pi from configuration pi**1:
 $contrl scftyp=gvb mult=2 $end
 $scf    nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.25
         alpha(1)= 2.0, 0.50, 0.0
          beta(1)=-1.0,-0.25, 0.0  $end

!       3-Sigma-minus from configuration pi**2:
 $contrl scftyp=gvb mult=3 $end
 $scf    nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50
         alpha(1)= 2.00, 1.00, 0.50
          beta(1)=-1.00,-0.50,-0.50 $end

!       1-Delta from configuration pi**2:
 $contrl scftyp=gvb mult=1 $end
 $scf    nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50
         alpha(1)= 2.00, 1.00, 0.00
          beta(1)=-1.00,-0.50,+0.50 $end

!       1-Sigma-plus from configuration pi**2:
 $contrl scftyp=gvb mult=1 $end
 $scf    nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.50
         alpha(1)= 2.00, 1.00, 0.25 
          beta(1)=-1.00,-0.50, 0.00 $end

!       2-Pi from configuration pi**3:
 $contrl scftyp=gvb mult=2 $end
 $scf    nco=xx nseto=1 no=2 couple=.true. f(1)=1.0,0.75
         alpha(1)= 2.00, 1.50, 1.00
          beta(1)=-1.00,-0.75,-0.50 $end

The same inputs apply to delta**N configurations: only the 
names of the energy terms (states) change!

    It is possible to do every term arising from the 
s**1,p**N configurations.  Roothaan-style A,B parameters 
are as follows:
     Acc = Bcc = 1.0 (as always)
     Aco = Bco = 1.0, for o=s and for o=p (as always)
     Ass = Bss = 0.0, since there is no other s e- to repel
     Asp = 1.0
     Bsp = -3K101 + 1 
     App and Bpp: take from the parent term in config p**N
and then apply the rule turning A,B into ALPHA,BETA, for 
Fc=1, Fs=0.5, Fp=N/6.  K101 should be taken from Roothaan 
and Bagus in Methods Comput. Phys. 2, 49(1963).

                       - - - -

    It may be useful to obtain a "configuration average", 
especially whenever a specific term energy cannot be 
computed.  Boris Plakhutin in Novosibirsk has kindly 
provided a recipe for the Roothaan A,B values that give the 
average energy of configuration O**N, meaning N electrons 
in some open shell O.  Let D be the degeneracy of the open 
shell O:
                 atoms: O=p, d, f... has D= 3, 5, 7...
  non-linear molecules: O=e, t, g, h has D= 2, 3, 4, 5.
The fractional occupation number F-o = N/(2*D).  The open 
shell Roothaan couplings are A-oo = B-oo = (N-1)/(N-F).

For example, the case d**2 has F-o=1/5 so A-oo=B-oo=5/9.  
Using the formulae to generate all ALPHA and BETA causes 
the GVB program to optimize the following spin-and-space 
weighted average of all terms appearing in d**2:
   E-avg = [21xE(3-F) + 9xE(3-P) +
             9xE(1-G) + 5xE(1-D) + E(1-S)]/45
by this input
 $contrl scftyp=gvb mult=3 $end
 $scf    nco=xx nseto=1 no=5 couple=.true.
             f(1)= 1.0, 0.2
         alpha(1)= 2.00, 0.40, 0.044444444444444444
          beta(1)=-1.00,-0.20,-0.022222222222222222 $end
Note that here, as in all other runs with the GVB module, 
the spin in $CONTRL is irrelevant: the energy that is 
optimized is determined by F, ALPHA, and BETA.  In this 
last case, these data implicitly include both triplet and 
singlet determinants in the energy averaging.

    Boris Plakhutin also provided a formula in case you 
want to average over a particular spin S within an open 
shell O**N with some degeneracy D implying a fractional 
occupancy of F.  The result is
              N**2(N-1) + FxN(N-4) + 4FxS(S+1)
       A-oo = --------------------------------
              4F(N-1) + N(N-4) + 4S(S+1)
       B-oo = --------------------------
For the d**2 configuration, the space-and-spin degeneracy 
weighted average of the 3-F and 3-P terms has Aoo=5/8 and 
Boo=10/8, while the average of the 1-G, 1-D, and 1-S terms 
has Aoo=5/12 and Boo=-10/12.  These optimize orbitals for 
the average
      E-triplet = [21xE(3-F) + 9xE(3-P)]/30
which is not quite the same as the space degeneracy average 
shown above, and for the average
      E-singlet = [9xE(1-G) + 5xE(1-D) + E(1-S)]/15
The literature reference for the overall and the spin-
specific configurational averages is
   B.N.Plakhutin, J. Mathematical Chem. 22, 203-233(1997)

true GVB perfect pairing runs
    True GVB runs are obtained by choosing NPAIR nonzero. 
If you wish to have some open shell electrons in addition 
to the geminal pairs, you may add the pairs to the end of 
any of the GVB coupling cases shown above.  The GVB module 
assumes that you have reordered your MOs into the order: 
NCO double occupied orbitals, NSETO sets of open shell 
orbitals, and NPAIR sets of geminals (with NORDER=1 in the 
$GUESS group).
    Each geminal consists of two orbitals and contains two 
singlet coupled electrons (perfect pairing).  The first MO 
of a geminal is probably heavily occupied (such as a 
bonding MO u), and the second is probably weakly occupied 
(such as an antibonding, correlating orbital v).  If you 
have more than one pair, you must be careful that the 
initial MOs are ordered u1, v1, u2, v2..., which is -NOT- 
the same order that RHF starting orbitals will be found in.  
Use NORDER=1 to get the correct order.
    These pair wavefunctions are actually a limited form of 
MCSCF.  GVB runs are much faster than MCSCF runs, because 
the natural orbital u,v form of the wavefunction permits a 
Fock operator based optimization.  However, convergence of 
the GVB run is by no means assured.  The same care in 
selecting the correlating orbitals that you would apply to 
an MCSCF run must also be used for GVB runs.  In 
particular, look at the orbital expansions when choosing 
the starting orbitals, and check them again after the run 
    GVB runs will be carried out entirely in orthonormal 
natural u,v form, with strong orthogonality enforced on the 
geminals.  Orthogonal orbitals will pervade your thinking 
in both initial orbital selection, and the entire orbital 
optimization phase (the CICOEF values give the weights of 
the u,v orbitals in each geminal).  However, once the 
calculation is converged, the program will generate and 
print the nonorthogonal, generalized valence bond orbitals.  
These GVB orbitals are an entirely equivalent way of 
presenting the wavefunction, but are generated only after 
the fact.
    Convergence of true GVB runs is by no means as certain 
as convergence of RHF, UHF, ROHF, or GVB with NPAIR=0. You 
can assist convergence by doing a preliminary RHF or ROHF 
calculation, and use these orbitals for GUESS=MOREAD. Few, 
if any, GVB runs with NPAIR non-zero will converge without 
using GUESS=MOREAD.  Generation of MVOs during the 
preliminary SCF can also be advantageous.  In fact, all the 
advice outlined for MCSCF computations below is germane, 
for GVB-PP is a type of MCSCF computation.
    The total number of electrons in the GVB wavefunction
is given by the following formula:

        NE = 2*NCO + sum 2*F(i)*NO(i) + 2*NPAIR

The charge is obtained by subtracting the total number of 
protons given in $DATA.  The multiplicity is implicit in 
the choice of alpha and beta constants.  Note that ICHARG 
and MULT must be given correctly in $CONTRL anyway, as the 
number of electrons from this formula is double checked 
against the ICHARG value.
the special case of TCSCF
    The wavefunction with NSETO=0 and NPAIR=1 is called 
GVB-PP(1) by Goddard, two configuration SCF (TCSCF) by 
Schaefer or Davidson, and CAS-SCF with two electrons in two 
orbitals by others.  Note that this is just semantics, as 
these are identical.  This is a very important type of 
wavefunction, as TCSCF is the minimum acceptable treatment 
for singlet biradicals.  The TCSCF wavefunction can be 
obtained with SCFTYP=MCSCF, but it is usually much faster 
to use the Fock based SCFTYP=GVB.  Because of its 
importance, the TCSCF function (together with open shells, 
if desired) permits analytic hessian computation.

a caution about symmetry
    Caution!  Some exotic calculations with the GVB program 
do not permit the use of symmetry.  The symmetry algorithm 
in GAMESS was "derived assuming that the electronic charge 
density transforms according to the completely symmetric 
representation of the point group", Dupuis/King, JCP, 68, 
3998(1978).   This may not be true for certain open shell 
cases, and in fact during GVB runs, it may not be true for 
closed shell singlet cases!  

    First, consider the following correct input for the 
singlet-delta state of NH:
 $SCF    NCO=3 NSETO=2 NO(1)=1,1 $END
for the x**1y**1 state, or for the x**2-y**2 state,
 $SCF    NCO=3 NPAIR=1 CICOEF(1)=0.707,-0.707 $END
Neither gives correct results, unless you enter NOSYM=1. 
The electronic term symbol is degenerate, a good tip off 
that symmetry cannot be used.  However, some degenerate 
states can still use symmetry, because they use coupling 
constants averaged over all degenerate states within a 
single term, as is done in EXAM15 and EXAM16.  Here the 
"state averaged SCF" leads to a charge density which is 
symmetric, and these runs can exploit symmetry.
    Secondly, since GVB runs exploit symmetry for each of 
the "shells", or type of orbitals, some calculations on 
totally symmetric states may not be able to use symmetry. 
An example is CO or N2, using a three pair GVB to treat the 
sigma and pi bonds.  Individual configurations such as 
(sigma)**2,(pi-x)**2,(pi-y*)**2 do not have symmetric 
charge densities since neither the pi nor pi* level is 
completely filled.  Correct answers for the sigma-plus 
ground states result only if you input NOSYM=1.
   Problems of the type mentioned should not arise if the 
point group is Abelian, but will be fairly common in linear 
molecules.  Since GAMESS cannot detect that the GVB 
electronic state is not totally symmetric (or averaged to 
at least have a totally symmetric density), it is left up 
to you to decide when to input NOSYM=1.  If you have any 
question about the use of symmetry, try it both ways.  If 
you get the same energy, both ways, it remains valid to use 
symmetry to speed up your run.

   And beware!  Brain dead computations, such as RHF on 
singlet O2, which actually is a half filled degenerate 
shell, violate the symmetry assumptions, and also violate 
nature.  Use of partially filled degenerate shells always 
leads to very wild oscillations in the RHF energy, which is 
how the program tries to tell you to think first, and 
compute second.  Configurations such as pi**2, e**1, or 
f2u**4 can be treated, but require GVB wavefunctions and F, 
ALPHA, BETA values from the sources mentioned.

created on 7/7/2017