* * * Other open shell SCF cases (GVB) * * *
Genuine GVB-PP 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.
It is possible to do other open shell cases with the GVB code, which can handle
the following cases:
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 ROHF module can do better.
Note
that all of
these
cases
are
really
ROHF,
since
the
default
for
NPAIR in $SCF
is 0.
Many open shell states with degenerate open shells (for example, in diatomic
molecules) can be treated as well. There is a sample of this in the 'Input
Examples' section of this manual.
If you would like to do any cases other than those shown above, you must derive
the coupling coefficients ALPHA and BETA, and input them with the
occupancies F
in the $SCF group.
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
derived from the values in Appendix A of the book "A General SCF Theory" by
Ramon Carbo and
Josep M. Riera, Springer-Verlag (1978). The basic rule is
(1) F(i) = 1/2 * omega(i)
(2) ALPHA(i) = alpha(i)
(3) BETA(i) = - beta(i),
where omega, alpha, and beta are the names used by Ramon in his Tables.
The variable NSETO should give the number of open shells, and NO should give the
degeneracy
of each open shell. Thus the 5-S state of carbon would have NSETO=2, and
NO(1)=1,3.
Some specific examples, for the lowest term in each of the atomic P N
configurations are
! 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.333333333333333
ALPHA(1)= 2.0 0.66666666666667 0.16666666666667
BETA(1)= -1.0 -0.33333333333333 -0.16666666666667 $END
! p**3 4-S state
$CONTRL SCFTYP=ROHF MULT=4 $END
! 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
! 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
Be sure to give all the digits, as these are part of a double precision energy
formula.
Coupling constants for d N configurations are
! 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
! 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
The source for these values is R.Poirier, R.Kari, and I.G.Csizmadia's book
Handbook of Gaussian
Basis Sets, Elsevier, Amsterdam, 1985.
Note that GAMESS can do a proper calculation on the ground terms for the d 2 , d
3 , d 7 , and d 8
configurations only by means of state averaged MCSCF. For d 8 , use
$CONTRL SCFTYP=MCSCF MULT=3 $END
$DRT GROUP=C1 FORS=.TRUE. NMCC=xx NDOC=3 NALP=2 $END
$GUGDIA NSTATE=10 $END
$GUGDM2 WSTATE(1)=1,1,1,1,1,1,1,0,0,0 $END
Open shell cases such as s**1,d**n are probably most easily tackled with the
state-averaged
MCSCF program.
* * * 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
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 all 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
(if
desired,
with
possible
open
shells)
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 1 y 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.