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
formulae:
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
factor.
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
difficult.
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
appropriate.
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
problems.
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
TROUBLE WITH DIIS IS WHEN THE ENERGY RISES TO ITS FINAL
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
events.
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
cleanup.
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,
0.666666666666666667
beta(1)=0.666666666666666667,
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
processes.
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:
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO(1)=1 $END
two open shells, triplet:
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=2 NO(1)=1,1 $END
two open shells, singlet:
$CONTRL SCFTYP=GVB MULT=1 $END
$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)
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
$CONTRL SCFTYP=GVB MULT=2 $END
$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
$CONTRL SCFTYP=GVB MULT=3 $END
$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
$CONTRL SCFTYP=ROHF MULT=4 $END
which is equivalent to
$CONTRL SCFTYP=GVB MULT=4 $END
$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
$CONTRL SCFTYP=GVB MULT=3 $END
$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
$CONTRL SCFTYP=GVB MULT=2 $END
$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
$CONTRL SCFTYP=GVB MULT=2 $END
$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
$CONTRL SCFTYP=GVB MULT=3 $END
$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
$CONTRL SCFTYP=GVB MULT=4 $END
$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
$CONTRL SCFTYP=GVB MULT=5 $END
$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
$CONTRL SCFTYP=ROHF MULT=6 $END
! d**6 5-D state
$CONTRL SCFTYP=GVB MULT=5 $END
$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
$CONTRL SCFTYP=GVB MULT=4 $END
$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
$CONTRL SCFTYP=GVB MULT=3 $END
$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
$CONTRL SCFTYP=GVB MULT=2 $END
$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:
http://userpage.fu-berlin.de/~dandrae/
openshell/openls/openls.html
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
program
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 = --------------------------------
N(N**2-4F**2)
and
4F(N-1) + N(N-4) + 4S(S+1)
B-oo = --------------------------
(N**2-4F**2)
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
converges.
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
i
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:
$CONTRL SCFTYP=GVB NOSYM=1 $END
$SCF NCO=3 NSETO=2 NO(1)=1,1 $END
for the x**1y**1 state, or for the x**2-y**2 state,
$CONTRL SCFTYP=GVB NOSYM=1 $END
$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