Second Order Perturbation Theory
The perturbation theory techniques available in GAMESS
expand to the second order energy correction only, but
permit use of nearly any zeroth order SCF wavefunction.
Since MP2 theory for systems well described by the chosen
zeroth order reference recovers about 80-85% of the
dynamical correlation energy (assuming the use of large
basis sets), MP2 is often a computationally effective
theory. For higher accuracy, you can instead choose the
more time consuming coupled cluster theory. When using
MPLEVL=2, it is important to ensure that your system is
well described at zeroth order by your choice of SCFTYP.
The input for second order pertubation calculations
based on SCFTYP=RHF, UHF, or ROHF is found in $MP2, while
for SCFTYP=MCSCF, see $MRMP.
By default, frozen core MP2 calculations are performed.
RHF and UHF reference MP2
These methods are well defined, due to the uniqueness of
the Fock matrix definitions. These methods are also well
understood, so there is little need to say more, except to
point out an overview article on RHF or UHF MP2 gradients:
C.M.Aikens, S.P.Webb, R.L.Bell, G.D.Fletcher,
M.W.Schmidt, M.S.Gordon
Theoret.Chem.Acc. 110, 233-253(2003)
The distributed memory parallel MP2 gradient program is
described in
G.D.Fletcher, M.W.Schmidt, M.S.Gordon
Adv.Chem.Phys. 110, 267-294(1999)
and that for UMP2 in
C.M.Aikens, M.S.Gordon
J.Phys.Chem.A 108, 3103-3110(2004)
One point which may not be commonly appreciated is that
the density matrix for the first order wavefunction for the
RHF and UHF case, which is generated during gradient runs
or if properties are requested in the $MP2 group, is of the
type known as "response density", which differs from the
more usual "expectation value density". The eigenvalues of
the response density matrix (which are the occupation
numbers of the MP2 natural orbitals) can therefore be
greater than 2 for frozen core orbitals, or even negative
values for the highest 'virtual' orbitals. The sum is of
course exactly the total number of electrons. We have seen
values outside the range 0-2 in several cases when the
single configuration HF wavefunction is not an appropriate
description of the system, and thus these occupancies may
serve as a guide to the wisdom of using a HF reference:
M.S.Gordon, M.W.Schmidt, G.M.Chaban, K.R.Glaesemann,
W.J.Stevens, C.Gonzalez J.Chem.Phys. 110,4199-4207(1999)
high spin ROHF reference MP2
There are a number of open shell perturbation theories
described in the literature. It is important to note that
these methods give different results for the second order
energy correction, reflecting ambiguities in the selection
of the zeroth order Hamiltonian and in defining the ROHF
Fock matrices. See
K.R.Glaesemann, M.W.Schmidt
J.Phys.Chem.A 114, 8772-8777(2010)
for a figure showing 4 different ROHF-based perturbation
theory potentials, which are highly parallel, but have
different total energy values.
Two of the perturbation theories mentioned below, RMP
and ZAPT, are available in GAMESS using SCFTYP=ROHF (see
OSPT in $MP2). Nuclear gradients can be obtained for ZAPT.
The OPT1 results can be generated using MPLEVL=2 with
SCFTYP=MCSCF, using an active space where every orbital is
singly occupied with the highest MULT possible (which is
single-determinant).
One theory is known as RMP, which it should be pointed
out, is entirely equivalent to the ROHF-MBPT2 method and to
the CUMP2 method. The perturbation theory is as UHF-like
as possible, and can be chosen in GAMESS by selection of
OSPT=RMP. The second order energy is defined by
1. P.J.Knowles, J.S.Andrews, R.D.Amos, N.C.Handy,
J.A.Pople Chem.Phys.Lett. 186, 130-136(1991)
2. W.J.Lauderdale, J.F.Stanton, J.Gauss, J.D.Watts,
R.J.Bartlett Chem.Phys.Lett. 187, 21-28(1991).
The submission dates are in inverse order of publication
dates, and -both- papers should be cited when using this
method. Here we will refer to the method as RMP in keeping
with much of the literature. The RMP method diagonalizes
the alpha and beta Fock matrices separately, so their
occupied-occupied and virtual-virtual blocks are
canonicalized. This generates two distinct orbital sets,
whose double excitation contributions are processed by the
usual UHF MP2 program, but an additional energy term from
single excitations is required. The recent CUHF method
converges a UHF calculation directly to these semi-
canonical orbitals, rather than generating them after ROHF
converges. The CUMP2 perturbation theory is thus identical
to methods described above:
3. G.E.Scuseria, T.Tsuchimochi
J.Chem.Phys. 134, 064101/1-14(2011)
CUMP2's input is SCFTYP=UHF, MPLEVL=2, CUHF=.TRUE. rather
than SCFTYP=ROHF, MPLEVL=2, OSPT=RMP.
RMP's use of different orbitals for different spins adds
to the CPU time required for integral transformations, of
course. just like UMP2. RMP is invariant under all of the
orbital transformations for which the ROHF itself is
invariant. Unlike UMP2, the second order RMP energy does
not suffer from spin contamination, since the reference
ROHF wave-function has no spin contamination. The RMP
wavefunction, however, is spin contaminated at 1st and
higher order, and therefore the 3rd and higher order RMP
energies are spin contaminated. Other workers have
extended the RMP theory to gradients and hessians at second
order, and to fourth order in the energy,
3. W.J.Lauderdale, J.F.Stanton, J.Gauss, J.D.Watts,
R.J.Bartlett J.Chem.Phys. 97, 6606-6620(1992)
4. J.Gauss, J.F.Stanton, R.J.Bartlett
J.Chem.Phys. 97, 7825-7828(1992)
5. D.J.Tozer, J.S.Andrews, R.D.Amos, N.C.Handy
Chem.Phys.Lett. 199, 229-236(1992)
6. D.J.Tozer, N.C.Handy, R.D.Amos, J.A.Pople, R.H.Nobes,
Y.Xie, H.F.Schaefer Mol.Phys. 79, 777-793(1993)
We deliberately omit references to the ROMP precursor of
the RMP formalism. RMP gradients are not available.
The Z-averaged perturbation theory (ZAPT) formalism for
ROHF perturbation theory is the preferred implementation of
open shell spin-restricted perturbation theory (OSPT=ZAPT
in $MP2). The ZAPT theory has only a single set of
orbitals in the MO transformation, and therefore runs in a
time similar to the RHF perturbation code. The second
order energy is free of spin-contamination, but some spin-
contamination enters into the first order wavefunction (and
hence properties). This should be much less contamination
than for OSPT=RMP. For these reasons, OSPT=ZAPT is the
default open shell method.
References for ZAPT are
7. T.J.Lee, D.Jayatilaka Chem.Phys.Lett. 201, 1-10(1993)
8. T.J.Lee, A.P.Rendell, K.G.Dyall, D.Jayatilaka
J.Chem.Phys. 100, 7400-7409(1994)
The formulae for the seven terms in the ZAPT energy are
clearly summarized in the paper
9. I.M.B.Nielsen, E.T.Seidl
J.Comput.Chem. 16, 1301-1313(1995)
The ZAPT gradient equations are found in
10. G.D.Fletcher, M.S.Gordon, R.L.Bell
Theoret.Chem.Acc. 107, 57-70(2002)
11. C.M.Aikens, G.D.Fletcher, M.W.Schmidt, M.S.Gordon
J.Chem.Phys. 124, 014107/1-14(2006)
We would like to thank Tim Lee for his gracious assistance
in the implementation of the ZAPT energy.
There are a number of other open shell theories, with
names such as HC, OPT1, OPT2, and IOPT. The literature for
these is
12. I.Hubac, P.Carsky Phys.Rev.A 22, 2392-2399(1980)
13. C.Murray, E.R.Davidson
Chem.Phys.Lett. 187,451-454(1991)
14. C.Murray, E.R.Davidson
Int.J.Quantum Chem. 43, 755-768(1992)
15. P.M.Kozlowski, E.R.Davidson
Chem.Phys.Lett. 226, 440-446(1994)
16. C.W.Murray, N.C.Handy
J.Chem.Phys. 97, 6509-6516(1992)
17. T.D.Crawford, H.F.Schaefer, T.J.Lee
J.Chem.Phys. 105, 1060-1069(1996)
The latter two of these give comparisons of the various
high spin methods, and the numerical results in ref. 17 are
the basis for the conventional wisdom that restricted open
shell theory is better convergent with order of the
perturbation level than unrestricted theory. Paper 8 has
some numerical comparisons of spin-restricted theories as
well. We are aware of one paper on low-spin coupled open
shell SCF perturbation theory
18. J.S.Andrews, C.W.Murray, N.C.Handy
Chem.Phys.Lett. 201, 458-464(1993)
but this is not implemented in GAMESS. See the MCSCF
reference perturbation code for this case.
GVB based MP2
This is not implemented in GAMESS. Note that the MCSCF
perturbation program discussed below should be able to
develop the perturbation corrections to open shell
singlets, by using a $DRT input such as
NMCC=N/2-1 NDOC=0 NAOS=1 NBOS=1 NVAL=0
which generates a single CSF if the two open shells have
different symmetry, or for a one pair GVB function
NMCC=N/2-1 NDOC=1 NVAL=1
which generates a 3 CSF function entirely equivalent to
the two configuration TCSCF, a.k.a GVB-PP(1). For the
record, note that if we attempt a triplet state with the
MCSCF program,
NMCC=N/2-1 NDOC=0 NALP=2 NVAL=0
we get a result equivalent to the OPT1 open shell method
described above, not the RMP or ZAPT result. It is
possible to generate the orbitals with a simpler SCF
computation than the MCSCF $DRT examples just given, and
read them into the MCSCF based MP2 program described below,
by RDVECS=.TRUE..
MCSCF reference perturbation theory
Just as for the open shell case, there are several ways
to define a multireference perturbation theory. The most
noteworthy are the CASPT2 method of Roos' group, the MRMP2
method of Hirao, the closely related MCQDPT2 method of
Nakano, and the MROPTn methods of Davidson. Although the
total energies of each method are different, energy
differences should be rather similar. In particular, the
MRMP/MCQDPT method implemented in GAMESS gives results for
the singlet-triplet splitting of methylene in close
agreement to CASPT2, MRMP2(Fav), and MROPT1, and differs by
2 Kcal/mole from MRMP2(Fhs), and the MROPT2 to MROPT4
methods.
The MCQDPT method implemented in GAMESS is a multistate
perturbation theory due to Nakano. If applied to 1 state,
it is the same as the MRMP model of Hirao. When applied to
more than one state, it is of the philosophy "perturb
first, diagonalize second". This means that perturbations
are made to both the diagonal and off-diagonal elements to
give an effective Hamiltonian, whose dimension equals the
number of states being treated. The effective Hamiltonian
is diagonalized to give the second order state energies.
Diagonalization after inclusion of the off-diagonal
perturbation ensures that avoided crossings of states of
the same symmetry are treated correctly. Such an avoided
crossing is found in the LiF molecule, as shown in the
first of the two papers on the MCQDPT method:
H.Nakano, J.Chem.Phys. 99, 7983-7992(1993)
H.Nakano, Chem.Phys.Lett. 207, 372-378(1993)
The closely related single state "diagonalize, then
perturb" MRMP model is discussed by
K.Hirao, Chem.Phys.Lett. 190, 374-380(1992)
K.Hirao, Chem.Phys.Lett. 196, 397-403(1992)
K.Hirao, Int.J.Quant.Chem. S26, 517-526(1992)
K.Hirao, Chem.Phys.Lett. 201, 59-66(1993)
Computation of reference weights and energy contributions
is illustrated by
H.Nakano, K.Nakayama, K.Hirao, M.Dupuis
J.Chem.Phys. 106, 4912-4917(1997)
T.Hashimoto, H.Nakano, K.Hirao
J.Mol.Struct.(THEOCHEM) 451, 25-33(1998)
Single state MCQDPT computations are very similar to MRMP
computations. A beginning set of references to the other
multireference methods used includes:
P.M.Kozlowski, E.R.Davidson
J.Chem.Phys. 100, 3672-3682(1994)
K.G.Dyall J.Chem.Phys. 102, 4909-4918(1995)
B.O.Roos, K.Andersson, M.K.Fulscher, P.-A.Malmqvist,
L.Serrano-Andres, K.Pierloot, M.Merchan
Adv.Chem.Phys. 93, 219-331(1996).
and a review article is available comparing these methods,
E.R.Davidson, A.A.Jarzecki in "Recent Advances in Multi-
reference Methods" K.Hirao, Ed. World Scientific, 1999,
pp 31-63.
The CSF (GUGA-based) MRMP/MCQDPT code was written by
Haruyuki Nakano, and was interfaced to GAMESS by him in the
summer of 1996. This program makes extensive use of disk
files during its specialized transformations and the
perturbation steps. Its efficiency is improved if you can
add extra physical memory to reduce the number of file
reads. In practice we have used this program up to about
12 active orbitals, and with very large disks, to about 500
AOs. In 2005, Joe Ivanic programmed a determinant based
MRMP/MCQDPT program. This uses the normal integral
transformation routines already present in GAMESS, and
direct CI technology to avoid disk I/O. The determinant
program is able to handle larger active spaces than the CSF
program, and has already been used for cases with 16
electrons in 16 orbitals, and basis sets up to 500 AOs.
When proper care is taken with numerical cutoffs, such
as CI vector convergence and the generator cutoff in the
CSF code, both programs produce identical results. Both
are enabled for parallel execution. The more mature CSF
program has several interesting options not found in the
determinant program: perturbative treatment of spin-orbit
coupling, energy denominators which are a band-aid for the
horribly named "intruder states", and the ability to find
the weight of the MCSCF reference in the 1st order
wavefunction. Neither program produces a density matrix
for property evaluation, nor are analytic gradients
programmed.
Second order perturbation corrections are also available
for ORMAS-type references, due to a program by Luke Roskop.
This is accessed by MRPT=DETMRPT, defining the reference by
CISTEP=ORMAS in $MCSCF, with a $ORMAS group. In its
present form, this program is limited to cases with equal
numbers of alpha and beta spins (i.e. singlets, triplets,
... with SZ=0, but not doublets, quartets...). See
L.Roskop, M.S.Gordon J.Chem.Phys. 135, 044101/1-11(2012)
Finally, note the existence of the MRMP=GMCPT program by
Haruyuki Nakano, for other non-CAS references.
- - - - - - - - - - - -
We end with an input example to illustrate open shell
reference and multi-reference pertubation computations on
the ground state of NH2 radical:
! 2nd order perturbation test on NH2, following
! T.J.Lee, A.P.Rendell, K.G.Dyall, D.Jayatilaka
! J.Chem.Phys. 100, 7400-7409(1994), Table III.
! State is 2-B-1, 69 AOs, either 1 or 49 CSFs.
!
! For 1 CSF reference,
! E(ROHF) = -55.5836109825
! E(ZAPT) = -55.7763947115
! [E(ZAPT) = -55.7763947289 at lit's ZAPT geom]
! E(RMP) = -55.7772299958
! E(OPT1) = -55.7830422945
! [E(OPT1) = -55.7830437413 at lit's OPT1 geom]
!
! For 49 CSF full valence MCSCF reference,
! CSFs: E(MRMP2) = -55.7857440268
! dets: E(MRMP2) = -55.7857440267
!
$contrl scftyp=mcscf mplevl=2 runtyp=energy mult=2 $end
$system mwords=1 memddi=1 $end
$guess guess=moread norb=69 $end
$mcscf fullnr=.true. $end
!
! Next set of lines carry out a MRMP computation,
! after a preliminary MCSCF orbital optimization.
!
! using determinants
$det stsym=B1 ncore=1 nact=6 nels=7 $end
! using CSFs, for the very same calculation.
--- $mcscf cistep=guga $end
--- $drt group=c2v stsym=B1 fors=.t.
--- nmcc=1 ndoc=3 nalp=1 nval=2 $end
--- $mrmp mrpt=mcqdpt $end
--- $mcqdpt stsym=B1 nmofzc=1 nmodoc=0 nmoact=6 $end
! Next lines carry out a single reference OPT1.
--- $det stsym=B1 ncore=4 nact=1 nels=1 $end
--- $mrmp mrpt=mcqdpt rdvecs=.true. $end
--- $mcqdpt nmofzc=1 nmodoc=3 nmoact=1 stsym=B1 $end
! Next lines are single reference RMP and/or ZAPT
--- $contrl scftyp=rohf $end
--- $mp2 ospt=rmp $end
$data
2-B-1 state...TZ2Pf basis, RMP geom. of Lee, et al.
Cnv 2
Nitrogen 7.0
S 6
1 13520.0 0.000760
2 1999.0 0.006076
3 440.0 0.032847
4 120.9 0.132396
5 38.47 0.393261
6 13.46 0.546339
S 2
1 13.46 0.252036
2 4.993 0.779385
S 1 ; 1 1.569 1.0
S 1 ; 1 0.5800 1.0
S 1 ; 1 0.1923 1.0
P 3
1 35.91 0.040319
2 8.480 0.243602
3 2.706 0.805968
P 1 ; 1 0.9921 1.0
P 1 ; 1 0.3727 1.0
P 1 ; 1 0.1346 1.0
D 1 ; 1 1.654 1.0
D 1 ; 1 0.469 1.0
F 1 ; 1 1.093 1.0
Hydrogen 1.0 0.0 0.7993787 0.6359684
S 3 ! note that this is unscaled
1 33.64 0.025374
2 5.058 0.189684
3 1.147 0.852933
S 1 ; 1 0.3211 1.0
S 1 ; 1 0.1013 1.0
P 1 ; 1 1.407 1.0
P 1 ; 1 0.388 1.0
D 1 ; 1 1.057 1.0
$end
OPT1 geom: H 1.0 0.0 0.7998834 0.6369401
RMP geom: H 1.0 0.0 0.7993787 0.6359684
ZAPT geom: H 1.0 0.0 0.7994114 0.6357666
E(ROHF)= -55.5836109825, E(NUC)= 7.5835449477, 9 ITERS
$VEC
...omitted...
$END
created on 7/7/2017