Correlation correction to configuration interaction singles from coupled cluster perturbation theory
Abstract
A new state specific correlation correction to configuration interaction singles (CIS) excitation energies is preseted using coupled cluster perturbation theory (CCPT). General expressions for CIS-CCPT are derived and expanded explicitly to first order in the wavefunction and second order in the energy. By virtue of the nature of CCPT this method is a priori size extensive and incorporates infinite order effects into the wavefunction. This results in a balanced singles space excited state theory that at second order is an improvement over the ubiquitous CIS(D) method and comparable in quality to equation of motion coupled cluster (EOM-CC). A modest test set composed of the first four excited states from nine small organic molecules was used to quantify the accuracy and consistency of the CIS-CCPT2 excitation energies and density of states. We find that CIS-CCPT2 has a standard deviation error of 0.18 eV for excitation energies and 0.14 eV for density of states compared to EOM-CC, a factor of two better than CIS(D) with a significant reduction in the maximum deviation as well.
I Introduction
The study of electronically excited states of molecules is an important point of intersection between active experimental investigations and theoretical methodological progress. With the ubiquitous use of lasers in optical spectroscopy the precision obtainable with modern experiments is difficult to replicate with modern theoretical methods. Despite this, theory can play an important role in optical spectroscopy by aiding with excited state assignments, predicting a priori the density of states (DOS) that should be expected in a certain spectrum range and even as a sanity check for the various complications that can arrise in an experimental apparatus.
While progress continues in perfecting time dependent density functional theory (TD-DFT) for use in excited statesDreuw and Head-Gordon (2005) the defacto standards remain either multireference configuration interactionSzalay et al. (2012) (MRCI) or equation of motion coupled clusterBartlett and Musiał (2007) (EOM-CC) theory. Typically truncated to only include single and double excitations, EOM-CCSDStanton and Bartlett (1993) is qualitatively consistent and approaches quantitative accuracy in many cases. Inclusion of triple (either iterativelyWatts and Bartlett (1994) or perturbativelyWatson Jr. et al. (2013)) excitations in the EOM framework generally brings the method into quantitative agreement with experiment. However EOM based methods come at a strong computational cost, with EOM-CCSD having a formal computational scaling of with significant storage and I/O costs which become demanding for larger molecules.
Approximations to the excited state problem have been under active research for many years now, the simplest and oldest being configuration interaction singles (CIS). Using CIS, the excited state is approximated as a linear combination of single excitations including the ground state Hartree-Fock wavefunction. As expected of a theory that is essentially a mean-field solution to excited states (lacking all correlation energy as such) the accuracy of CIS is akin to Hartree-Fock solutions for the ground state. Early attempts to add correlation to the CIS wavefunction using second order perturbation theory gave rise to CIS-MP2Foresman et al. (1992) and the more widely used size extensivity corrected CIS(D)Head-Gordon et al. (1994) method, with various improvements in accuracyRhee and Head-Gordon (); Liu et al. (2013); Liu and Subotnik (2014) and extensions to higher orders.Hirata (2005)
In this article, we approach the CIS correlation problem using coupled cluster perturbation theoryBartlett et al. (2010) (CCPT) with the demands that the resulting method should be a priori size-extensive, should contain no cluster opperators of higher particle rank than doubles and include essential infinite order contributions. We start in Sec. II by outlining the general form of single determinant CCPT, followed by deriving our modification to include the CIS wavefunction into a zero’th order ansatze in Sec. III. The resulting CIS-CCPT energy is explicitly expanded to second order in Sec. III.1, with exploratory numerical calculations presented in Sec. V.
Ii Coupled cluster perturbation theory
Single determinate CCPTBartlett et al. (2010) starts with the exponential wavefunction ansatz
(1) |
where the cluster operator,
(2) |
is defined in terms of the excitation operators (here limited to single and double excitations, with higher order operators defined similarly)
(3) | ||||
(4) |
Throughout this work the indices are reserved for particle states while refer to hole states. With this ansatz the similarity transformed Schrödinger equation can be written as
(6) |
where is the reference state projection operator, , and is the usual normal ordered Hamiltonian
(7) |
Here is the Fock one-electron matrix, is the antisymmetrized two-electron integral matrix and denotes normal ordering of the enclosed operators. The excitation operator amplitudes are obtained by solving the system of projected equations
(8) |
with defined to be the operator spanning the complimentary space to ,
(9) |
For convenience we also can refer to the complimentary space by the particle excitation rank as for the single excitation space, for double excitations and so forth.
Partitioning the Hamiltonian into a zero’th order and perturbation component,
(10) |
where the exact contribution to and is left undefined, the correlation energy (the zero’th order energy and wavefunction is defined to be the corresponding mean field solution) can be expressed to order as
(11) |
where
(12) |
is the similarity transformed perturbation Hamiltonian reduced to only connected (denoted by the subscript) terms by application of the Baker-Campbell-Hausdorff (BCH) identity. The ’th order exponential wavefunction is defined here to include all possible products of ’th () order cluster operators (including the ’th order contribution)
(13) |
The amplitudes of which can be determined by evaluating Eq. 8 order by order using
(14) |
Hamiltonian Partitioning | |
---|---|
MBPT | Particle Rank |
With the CCPT framework so defined, the remaining choice is how to partition the Hamiltonian. With future considerations in mind, we write the normal ordered Hamiltonian (Eq. 7) in terms of the one () and two () particle operators, separated by the associated particle excitation rank:
(15) |
Using the standard MBPT choice of Hamiltonian partitioning (see Table 1) naturally reproduces the usual MBPT energy and wavefunction,Shavitt and Bartlett (2009) however we are free to choose alternative partitionings should we desire. Motivated by the performance of LinearCCD (coupled cluster doubles truncated to remain linear in ) Bartlett et al.Bartlett et al. (2010) partition the Hamiltonian in terms of particle excitation rank. The effect of this partitioning choice is that any given order in perturbation theory infinite order contributions are included in the cluster operators, an important property found in coupled cluster based wavefunctions.
Iii CIS-CCPT(n)
Excitation Energies | ||||
CIS | CIS(D) | CIS-CCPT2 | EOM-CCSD | |
Acetone | ||||
Benzene | ||||
Formaldehyde | ||||
Furan | ||||
^{1}^{1}1CIS overlap %. | ||||
Naphthalene | ||||
Norbornadiene | ||||
Propanamide | ||||
Pyrazine | ||||
MD | ||||
MAD | ||||
RMSD | ||||
Relative Density of States Statistics | ||||
MD | ||||
MAD | ||||
RMSD |
To include the CIS wavefunction in the CCPT framework we redefine our zero’th order wavefunction to be
(16) |
where the subscript denotes that only terms from the exponential up to linear order (as indicated in Eq. 16) are retained and
(17) |
is a single excitation operator with amplitudes obtained using the configuration interaction singles (CIS) methodForesman, Head-Gordon, and Pople (1992) for the ’th state. Using Eq. 16 the projection operator becomes
(18) |
With this definition of the starting wavefunction the ’th order reference energy (given by ) becomes shifted by the excitation energy
(19) |
We continue to use the ground state projection operator defined as , but note that must be altered to remain orthogonal to . Substitution of Eq. 16 into Eq. 11 gives
(20) |
Subtracting the ground state correlation energy, (computed using Eq. 11), from Eq. 20 and expanding the denominator to cancel unlinked terms we obtain the excitation correlation energy
(21) |
where denotes that only linked terms are retained.
The associated amplitude projected equations with direct substitution of Eq. 16 into Eq. 14 are
(22) |
However this definition does not ensure that all the contributing operators are connected. We can strictly enforce connectivity by performing a double similarity transformation of Eq. 14 by premultiplying with before projecting against the complimentary space (an operation we are free to perform so long as is orthogonal to ), obtaining
(23) |
where Eq. 23 is obtained by repeated use of the BCH identity.
iii.1 Second-order energy
Symmetry | CIS^{2}^{2}2This work. | CIS(D)^{1}^{1}footnotemark: 1 | CIS(3)^{3}^{3}3Ref. Hirata, 2005. | CIS(4)^{2}^{2}footnotemark: 2 | MBPT(2)^{4}^{4}4Ref. Gwaltney, Nooijen, and Bartlett,1996 | CIS-CCPT2^{1}^{1}footnotemark: 1 | EOM-CCSD^{1}^{1}footnotemark: 1 | Exp.^{5}^{5}5Adopted from Ref. Head-Gordon et al.,1994 |
---|---|---|---|---|---|---|---|---|
^{6}^{6}6CIS overlap %. | ||||||||
The first order wave function defined by Eq. 23 is explicitly given as
(24) |
(25) |
with no higher excitations included at this order due to the connected nature of Eq. 23. The second order energy, obtained from Eq. 21 using the amplitudes defined in Eqs. 24 and 25, can be explicitly expressed as
(26) |
Assuming the use of canonical Hartree-Fock starting orbitals () the amplitude equations reduce to
(27) |
and by defining
(28) |
along with
(29) |
the total CIS-CCPT2 excitation energy is compactly expressed as
(30) |
For explicit spin-orbital equations corresponding to Eqs. 27-29 we refer to Appendix A.
The CIS-CCPT2 equations above introduce infinite order correlation to the CIS excitation energy root by root with the computational cost of an iterative scaling. This computational cost can be contrasted on one hand with EOM-CCSD, where the scaling limit is also with an additional cost factor of due to having both left and right hand vectors.^{7}^{7}7Additional storage and computer I/O scaling costs further penalize EOM-CCSD over CIS-CCPT2. On the other hand the prevailing CIS based second order perturbation theory, CIS(D), scales as but is not a priori size extensive. This failing is well known to plague the CIS(n) perturbation methods,Head-Gordon et al. (1994); Hirata (2005) requiring an ad hoc correction to restore extensivity.
Iv Electronic Structure Calculations
Symmetry | CIS^{8}^{8}8This work. | CIS(D)^{1}^{1}footnotemark: 1 | CIS(3)^{9}^{9}9Ref. Hirata, 2005. | CIS(4)^{2}^{2}footnotemark: 2 | MBPT(2)^{10}^{10}10Ref. Gwaltney, Nooijen, and Bartlett,1996 | CIS-CCPT2^{1}^{1}footnotemark: 1 | EOM-CCSD^{1}^{1}footnotemark: 1 | Exp.^{11}^{11}11Adopted from Ref. Head-Gordon et al.,1994 |
---|---|---|---|---|---|---|---|---|
^{12}^{12}12CIS overlap %. | ||||||||
To be consistent with the existing computational excited state literatureHirata (2005); Schreiber et al. (2008); Silva-Junior et al. (2010) referenced in this work, molecular equilibrium geometries were obtained using the MP2/6-31G* level of theory. All correlation calculations employ the frozen-core approximation, where the orbitals belonging to first row heavy atoms are excluded from the correlation space. Electronic structure calculations were performed on the University of Florida HiPerGator supercomputer using the massively parallel ACES IIILotrich et al. (2008) quantum chemistry program.
The use of the TZVP basis set from Schäfer et al.Schäfer, Horn, and Ahlrichs (1992) was considered here for retaining consistency with existing calculations. However because of the dependency of our CIS-CCPT2 method on the quality of the CIS wavefunction, and the tendency of the CIS method to over characterize states with Rydberg character necessitating the inclusion of diffuse functions in the basis set, we opted instead to use the significantly more expensive aug-cc-pVTZ correlation consistent basis set.Kendall, Dunning Jr, and Harrison (1992) Investigative calculations comparing the TZVP and aug-cc-pVTZ at the EOM-CCSD level of theory show good agreement for the first excitation energies of the molecules considered in Table 2. However, the excited state ordering and density of states for the higher excitation spectrum varied greatly between the two basis sets for a number of the molecules under consideration. This is easily attributable to the better characterization of the excited states through the inclusion of more angular momentum in the basis set ( contracted primitive set in the aug-cc-pVTZ basis while the TZVP basis set is limited to contracted primitives) in addition to the added diffuse functions. However this ordering of states and DOS difference demands great care when comparing with other published results computed with the smaller TZVP basis set.
V Numerical Results and Discussion
To assess the effectiveness of the CIS-CCPT2 method described here we consider two numerical metrics: absolute magnitude and ordering of the lowest four excitation energies and the energy spacing (namely the DOS) relative to the first excited state. These metrics were examined using a subset selection of small molecules from the TBE excited state test set.Schreiber et al. (2008); Silva-Junior et al. (2010) Additionally we examined in detail two small molecules (Acetaldehyde and Formaldehyde) known to be problematic for CIS based perturbation theory.Hirata (2005) All performance statistics (MD for maximum deviation, MAD for mean absolute deviation and RMSD for root mean square deviation) are computed relative to the corresponding EOM-CCSD root. Excited state results presented here are sorted by the CIS-CCPT2 energy ordering unless stated otherwise, with the corresponding EOM-CCSD energy assigned by symmetry.
The CIS-CCPT2 predicted excited state spectrum for a selection of small molecules (a subset of the TBE-1Schreiber et al. (2008) test set) is presented in Table 2 as well as the corresponding CIS and CIS(D) excitation energy for that same root. We find excellent agreement between the CIS-CCPT2 and EOM-CCSD predicted excitation energies for most of the molecules examined with an excellent mean deviation of only 0.16 eV. This can be compared with CIS(D), the prevalent second-order CIS based perturbation theory, which has a similarly respectable 0.23 eV average deviation. Further comparison shows however that the CIS-CCPT2 method is much more consistent than CIS(D) with a standard deviation of 0.18 eV compared to 0.33 eV as well as comparing the maximum deviation, 0.31 eV from CIS-CCPT2 and 0.8 eV for CIS(D).
Additionally we find that the CIS-CCPT2 predicted state ordering agrees exactly with the EOM-CCSD for most of the molecules listed in Table 2 with the exception of Furan and Propanamide. This is contrasted with the hit or miss predicted state ordering from CIS(D). In case of Propanamide the EOM-CCSD state ordering of is only narrowly determined by an energy difference of eV between the and states. Due to this small energy difference it is difficult to claim that the EOM-CCSD ordering is indeed definitive.^{13}^{13}13Exploratory calculations using EOM-CCSD(T) on these propanamide states have state orderings that agree with those predicted by CIS-CCPT2. To properly explain the state ordering failure in Furan it is necessary to examine the overlap of the CIS wavefunction with that of the associated EOM-CCSD wavefunction, given by
(31) |
where and are the left and right hand EOM-CCSD singles eigenvectors for the ’th root. For simplicity we approximate this overlap by assuming that . The closer is to unity the more accurate the CIS wavefunction is as zero’th-order wavefunction. For most states investigated here , however for the problematic case of the Furan state this overlap is poor (), indicative of a stronger doubles contribution to the final wavefunction.
We further examine the impact of a poor by computing the lowest six roots of Formaldehyde and Acetaldehyde, both molecules known to have states with significant Rydberg characterGwaltney and Bartlett (1995) and poor CIS overlap which is problematic for CIS based excited state perturbation theory.Hirata (2005) It is evident from Table 3 that for Formaldehyde the state, with a %, is poorly described by either the CIS(n) methods shown or CIS-CCPT2. This can be seen again for the miss-ordering of the state of Acetaldehyde, where the estimated CIS overlap is quite low ().
In addition to the absolute energy spectrum we also examine the excited DOS relative to the first excited state, the statistics of which are listed at the bottom of Table 2. An average deviation of 0.10 eV and a standard deviation of 0.14 eV for CIS-CCPT2 is quite promising (contrasted with the 0.19 eV average and 0.29 standard deviation for CIS(D)), this is further illustrated in Fig. 1 where the relative energy spectrum for the small molecules considered here are plotted. As can be seen, CIS-CCPT2 does well in matching the features exhibited in the EOM-CCSD excited state spectrum, suggesting that this method can be accurate as EOM-CCSD for predicting relative energy spacings. This is of great usefulness in computing solvent related shifts, a topic of significant interest to us.
Vi Conclusion
Coupled cluster perturbation theory has been used to derive a correlation correction to the configuration interaction singles reference wavefunction. The resulting CIS-CCPT method is a singles space multi-determinant perturbation theory with a well defined wavefunction that includes infinite-order effects and is a priori size extensive. The latter being an important consideration as the usual CIS based perturbation treatments are only size extensive by ad hoc omission of certain terms.Head-Gordon et al. (1994); Hirata (2005) The first order wavefunction and second order energy for CIS-CCPT has been explicitly presented in Section III.1 (see Eqs. 27-29) and implemented in the ACES IIILotrich et al. (2008) quantum chemistry package.
The quality of this CIS-CCPT2 method is demonstrated by examining the excited state spectrum of a number of small molecules (see Table 2). The resulting average (standard) deviation of 0.16 (0.18) eV as compared to EOM-CCSD is a marked improvement over CIS(D) with an average (standard) deviation of 0.23 (0.33) eV. Additionally, the density of states relative to the first excited state also compares well with EOM-CCSD, with an average (standard) deviation of 0.10 (0.14) eV. To illustrate this last comparison we plot in Fig. 1 the relative excited state spectrum of the molecules considered here, where CIS-CCPT2 can be seen to trace the features seen in the EOM-CCSD excited spectrum. It can also be noted that the CIS-CCPT2 and CIS(D) DOS generally bracket the corresponding EOM-CCSD spectrum, a convenient trend for identifying problematic states. To further investigate the behavior of CIS-CCPT2 for difficult cases we compare the first six excited states of Formaldehyde and Acetaldehyde (known to be problematic for CIS based perturbation theoryHirata (2005)) to both EOM-CCSD and experimental results. We find that CIS-CCPT2 has excellent agreement with the comparison values except in the expected situation where the CIS wavefunction overlap with the EOM-CCSD singles vector (Eq. 31) is poor.
Vii Acknowledgments
The authors would like to acknowledge funding support from the United States Air Force Office of Scientific Research.
Appendix A Spin-Orbital Equations for CIS-CCPT2
Explicit spin-orbital equations that define the CIS-CCPT2 amplitude and energy expressions are presented below. The amplitude Eq. 27 is (assuming Einstein notation)
(32) |
where
(33) |
are antisymmetric two electron integrals,
(34) |
are the Hartree-Fock orbital energies and
(35) |
is the antisymmetrized permutation operator. From the CIS-CCPT2 energy expression (Eq. 30) the the contribution (Eq. 28) is
(36) |
while (Eq. 29) is given by
(37) |
Appendix B Connection to CIS(n) perturbation theory
A satisfying connection to CIS(n) perturbation theory can be drawn from the CIS-CCPT2 equations with some judicious truncation and approximation. We start by approximating the equation from Eq. 25 by the zero’th iteration approximation
(38) |
using the denominator shifted resolvent (infinite order methods like CIS-CCPT(n) are of course invariant to denominator shifts)
(39) |
With this approximation Eq. 28 is reduced to
(40) |
which is generally small and so is neglected, while Eq. 29 becomes
(41) |
To finish the ad hoc link to the second order CIS(D) equations it is necessary to realize that in the CIS(n) derivation does not start with the ansatz but instead the fundamentally different starting wavefunction. This means that in the CIS(n) expansion only two of the four terms^{14}^{14}14These two surviving terms are generally of the same magnitude and opposite sign. in Eq. 41 can occur. Simplification of Eq. 41 directly gives the CIS(D) energy equation:
(42) |
References
- Dreuw and Head-Gordon (2005) A. Dreuw and M. Head-Gordon, Chem. Rev. 105, 4009 (2005).
- Szalay et al. (2012) P. G. P. Szalay, T. T. Müller, Gidofalvi, G. G. Lischka, and R. R. H. H.; Shepard, Chem. Rev. 112, 108 (2012).
- Bartlett and Musiał (2007) R. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
- Stanton and Bartlett (1993) J. F. Stanton and R. J. Bartlett, J. Chem. Phys. 98, 7029 (1993).
- Watts and Bartlett (1994) J. D. Watts and R. J. Bartlett, J. Chem. Phys. 101, 3073 (1994).
- Watson Jr. et al. (2013) T. J. Watson Jr., V. F. Lotrich, P. G. Szalay, A. Perera, and R. J. Bartlett, J. Phys. Chem. A 117, 2569 (2013).
- Foresman et al. (1992) J. B. Foresman, M. Head-Gordon, J. A. Pople, and M. J. Frisch, J. Phys. Chem. 96, 135 (1992).
- Head-Gordon et al. (1994) M. Head-Gordon, R. J. Rico, M. Oumi, and T. J. Lee, Chem. Phys. Lett. 219, 21 (1994).
- (9) Y. M. Rhee and M. Head-Gordon, .
- Liu et al. (2013) X. Liu, Q. Ou, E. Alguire, and J. E. Subotnik, J Chem. Phys. 138, 221105 (2013).
- Liu and Subotnik (2014) X. Liu and J. E. Subotnik, J. Chem. Theory Comput. 10, 1004 (2014).
- Hirata (2005) S. Hirata, J. Chem. Phys. 122, 094105 (2005).
- Bartlett et al. (2010) R. J. Bartlett, M. Musial, V. F. Lotrich, and T. Kus, in Recent Progress in Coupled-Cluster Methods, Vol. 11, edited by P. Carsky, J. Paldus, and J. Pittner (Springer, Dordrecht, 2010) Chap. 1, pp. 1–34.
- Shavitt and Bartlett (2009) I. Shavitt and R. J. Bartlett, Many-Body Methods in Chemistry and Physics (Cambridge, New York, 2009).
- Foresman, Head-Gordon, and Pople (1992) J. Foresman, M. Head-Gordon, and J. Pople, J. Phys. Chem. 96, 135 (1992).
- Gwaltney, Nooijen, and Bartlett (1996) S. R. Gwaltney, M. Nooijen, and R. J. Bartlett, Chem. Phys. Lett. 248, 189 (1996).
- (17) Additional storage and computer I/O scaling costs further penalize EOM-CCSD over CIS-CCPT2.
- Schreiber et al. (2008) M. Schreiber, M. R. Silva-Junior, S. P. A. Sauer, and W. Thiel, J Chem. Phys. 128, 134110 (2008).
- Silva-Junior et al. (2010) M. R. M. Silva-Junior, M. M. Schreiber, S. P. A. S. Sauer, and W. W. Thiel, J Chem. Phys. 133, 174318 (2010).
- Lotrich et al. (2008) V. Lotrich, N. Flocke, M. Ponton, A. Yau, A. Perera, E. Deumens, and R. Bartlett, J. Chem. Phys. 128, 2722 (2008).
- Schäfer, Horn, and Ahlrichs (1992) A. Schäfer, H. Horn, and R. Ahlrichs, J. Chem. Phys. 97, 2571 (1992).
- Kendall, Dunning Jr, and Harrison (1992) R. Kendall, T. Dunning Jr, and R. Harrison, J. Chem. Phys. 96, 6796 (1992).
- (23) Exploratory calculations using EOM-CCSD(T) on these propanamide states have state orderings that agree with those predicted by CIS-CCPT2.
- Gwaltney and Bartlett (1995) S. R. Gwaltney and R. J. Bartlett, Chem. Phys. Lett. 241, 26 (1995).
- (25) These two surviving terms are generally of the same magnitude and opposite sign.