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 

   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 

   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 

   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
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, 

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 

   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 

   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

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


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

created on 7/7/2017