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