A JeziorskiMonkhorst fully uncontracted MultiReference perturbative treatment I: principles, secondorder versions and tests on ground state potential energy curves
Abstract
The present paper introduces a new multireference perturbation approach developed at second order, based on a JeziorskyMokhorst expansion using individual Slater determinants as perturbers. Thanks to this choice of perturbers, an effective Hamiltonian may be built, allowing for the dressing of the Hamiltonian matrix within the reference space, assumed here to be a CASCI. Such a formulation accounts then for the coupling between the static and dynamic correlation effects. With our new definition of zerothorder energies, these two approaches are strictly sizeextensive provided that local orbitals are used, as numerically illustrated here and formally demonstrated in the appendix. Also, the present formalism allows for the factorization of all double excitation operators, just as in internally contracted approaches, strongly reducing the computational cost of these two approaches with respect to other determinantbased perturbation theories. The accuracy of these methods has been investigated on groundstate potential curves up to full dissociation limits for a set of six molecules involving single, double and triple bond breaking. The spectroscopic constants obtained with the present methods are found to be in very good agreement with the full configuration interaction (FCI) results. As the present formalism does not use any parameter or numerically unstable operation, the curves obtained with the two methods are smooth all along the dissociation path.
I Introduction
The research of the groundstate wave function of closedshell molecules follows wellestablished paths. The perturbative expansions from the meanfield HartreeFock single determinant usually converge and may be used as basic tools, especially when adopting a monoelectronic zeroorder Hamiltonian known as the MøllerPlesset Hamiltonian.(1) In this approach the wave function and the energy may be understood in terms of diagrams, which lead to the fundamental linkedcluster theorem.(2) The understanding of the sizeconsistency problem led to the suggestion of the Coupled Cluster approximation,(3); (4); (5); (6); (7) which is now considered as the standard and most efficient tool in the study of such systems in their ground state, especially in its CCSD(T) version where linked corrections by triple excitations are added perturbatively.(8) The situation is less evident when considering excited states, chemical reactions and molecular dissociations, since it then becomes impossible to find a relevant single determinant zeroorder wave function. These situations exhibit an intrinsic MultiReference (MR) character. A generalized linkedcluster theorem has been established by Brandow,(9) which gives a basis to the understanding of the sizeconsistency problem in this context, but the conditions for establishing this theorem are severe. They require a Complete Active Space (CAS) model space and a monoelectronic zeroorder Hamiltonian. Consequently, the corresponding QuasiDegenerate Perturbation Theory (QDPT) expansion cannot converge in most of the molecular MR situations(10); (11); (12). The research of theoretically satisfying (sizeconsistent) and numerically efficient MR treatments remains a very active field in Quantum Chemistry, as summarized in recent review articles concerning either perturbative(13) or CoupledCluster(14) methods.
The present work concentrates on the search of a new MR perturbative approach at second order (MRPT2). Of course, pragmatic proposals have been rapidly formulated, consisting first in the identification of a reference model space, defined on the set of single determinants having large components in the desired eigenstates of the problem. Diagonalizing the Hamiltonian in this reference space delivers a zeroorder wave function. Then one must define the vectors of the outer space to be used in the development and, in a perturbative context, choose a zeroorder Hamiltonian. The simplest approach consists in using single determinants as outerspace eigenvectors, and this has been used in the CIPSI method(15); (16) which is iterative, increasing the model space from the selection of the perturbing determinant of largest coefficients and their addition to the model space. From a practical point of view this method is very efficient and is now employed to reach near exact Full Configuration Interaction (FCI) energies on small molecules,(17); (18) and also as trial wave function in the context of quantum Monte Carlo.(19); (20); (21); (22) But the method suffers two main defects: i) it is not sizeconsistent and ii) it does not revise the modelspace component of the wave function under the effect of its interaction with the outerspace. This last defect is avoided if one expresses the effect of the perturbation as a change of the matrix elements of the model space CI matrix, according to the Intermediate Effective Hamiltonian (IEH) theory,(23) as done in the statespecific(24) or multistate(25) versions. Other methods which start from a CAS model space and use multideterminantal outerspace vectors have been proposed later on and are now broadly used. The first one is the CASPT2 method,(26); (27) which employs a monoelectronic zeroorder Hamiltonian. The method suffers from intruder state problems, to be cured in a pragmatic manner through the introduction of some parameters, and is not strictly sizeconsistent. The NEVPT2 method(28); (29); (30) also uses multideterminantal perturbers (defined in two different ways in its partially and strongly contracted versions), it makes use of a more sophisticated bielectronic Hamiltonian (the Dyall Hamiltonian(33)) to define the zeroorder energies of these perturbers, it is parameterfree, intruderstate free, and sizeconsistent. Both methods are implemented in several popular codes, and use a contracted description of the model space component of the desired eigenfunction (fixed by the diagonalization of the Hamiltonian in the model space). MultiState versions exist to give some flexibility to the model space component, in particular around weakly avoided crossings, but this flexibility is very limited(31); (32). If one returns to methods using singledeterminant perturbers, the origin of their sizeinconsistency problem has been identified as due to the unbalance between the multideterminant character of the zeroorder wave function and the single determinant character of the perturbers.(34) It is in principle possible to find size consistent formulations but they require rather complex formulations,(35); (36); (37); (38) and face some risk of numerical instabilities since they involve divisions by possibly small coefficients, the amplitudes of which may be small. Finally, one should mention a very recent approach based on the rewriting of the multireference linear coupled cluster equations in a stochastic framework of FullCI Quantum Monte Carlo which also uses single Slater determinants as perturbers(39).
The present paper is composed as follows. In Section II, the hereproposed formalism is presented, whose main features are:

it considers a CAS model space (to achieve the strict separability requirement), usually obtained from a preliminary CASSCF calculation;

the perturbers are single determinants (the method is externally noncontracted, according to the usual terminology);

it is statespecific, and strictly separable when localized active MOs are used (see formal demonstration in the appendix);

it makes use of the Dyall Hamiltonian to define the excitation energies appearing in the energy denominators;

it is based on a JeziorskiMonkhorst(40) (JM) expression of the wave operator and proceeds through referencespecific partitionings of the zeroorder Hamiltonian, as it has been previously suggested in the socalled MultiPartitionning(41); (42); (43) (MUPA) and also in the UGASSMPRT2.(38) Consequently, it does not define a unique zeroorder energy to the outerspace determinants;

it can be expressed either as a secondorder energy correction or as a dressing of the CASCI matrix, which offers a full flexibility in the treatment of the feedback effect of the postCASCI correlation on the model space component of the wave function;

the contributions of the various classes of excitations are easily identified (as in the CASPT2 and NEVPT2 methods);

thanks to our definition of the zerothorder energies, all processes involving double excitations can be treated by using only the one and twoboy density matrices, avoiding to loop on the perturbers;

given a set of molecular orbitals, it is parameter free and does not contain any threshold to avoid numerical instabilities.
After having presented the working equations of the present formalism, section III proposes a comparison with other existing MR approaches, such as some special cases of multireference coupled cluster (MRCC) and MRPT2. Then, Section V presents the numerical results for the ground state potential energy curves of six molecules involving single, double and triple bond breaking with both the JMMRPT2 and JMHeffPT2 methods. A numerical test of sizeextensivity is provided, together with the investigation of the dependency of the results on the locality of the active orbitals. Finally, Section VI summarizes the main results and presents its tentative developments. The reader can find in the Section Appendix: proof of separability a mathematical proof of strong separability of the JMMRPT2 method.
Ii Working equations for the perturbation and effective Hamiltonian at second order
As demonstrated previously by one of the present authors and his collaborators,(34) the sizeconsistency problem in any multireference perturbative expansion using single Slater determinants as perturbers comes from the unbalanced zeroth order energies that occur in the denominators. More precisely, if the zeroth order wave function is a CASCI eigenvector, its energy is stabilized by all the interactions within the active space, whereas a perturber treated as a single Slater determinant does not benefit from these extradiagonal Hamiltonian matrix elements. However, if one considers the set of perturber determinants created by the application of a given excitation operator on all Slater determinants of the CASCI wave function, most of the interactions found within the active space will occur within this set of perturbers (see Figure 1 for a pictorial example). Therefore, the use of linear combinations of Slater determinants as perturbers together with a bielectronic zeroth order operator, as it is the case in the NEVPT2 framework which uses the Dyall zeroth order operator, leads to balanced energy differences and removes the sizeconsistency problem.
On the basis of such considerations, the present work proposes an approach that uses single Slater determinants as perturbers and takes benefits of a new definition of energy denominators as expectation values of the Dyall zeroth order Hamiltonian over a specific class of linear combinations of Slater determinants. We first expose the definition of this perturbation theory, namely the JMMRPT2, which is strictly separable provided that local orbitals are used.
A large benefit from this new definition is that one may go beyond the sole calculation of the energy and improve the reference wave function by taking into account, in a strictly sizeconsistent way, the correlation effects brought by the perturbers on the reference space. In a second step, we reformulate the approach as a dressing of the Hamiltonian matrix within the set of Slater determinants belonging to the reference wave function, which is diagonalized. This approach will be referred to as the JMHeffPT2 method.
ii.1 The JMMRPT2 method
Firstorder perturbed wave function and secondorder energy
The formalism presented here is state specific and is not therefore restricted to ground state calculations. Nevertheless, for the sake of clarity and compactness, we omit the index referring explicitly to a specific eigenstate.
The zeroth order wave function is assumed to be a CASCI eigenvector expanded on the set of reference determinants :
(1) 
Such a wave function has a variational energy :
(2) 
Starting from a normalized (i.e. ), we assume that the exact wave function can be expressed as:
(3) 
where here the are all possible Slater determinants not belonging to the CASCI space. One should notice that such form is in principle not exact, as some changes of the coefficients within the CASCI space can formally occur when passing from the CASCI eigenvector to the FCI one, but such an approximated form for the exact wave function is the basis of many MRPT2 approaches like NEVPT2, CASPT2 or CIPSI.
As in any projection technique, the exact energy can be obtained by projecting the Schrödinger equation on :
(4)  
and one only needs to compute the coefficients of the that interact with , which consist in all individual Slater determinants being singly or doubly excited with respect to any Slater determinant in the CASCI space. From now on we implicitly refer to as any single Slater determinant belonging to such a space.
The coefficients are then written according to the JM ansatz(40), whose general expression for wave function is not explicitly needed here, and will be therefore given in the section III.2 when the comparison of the present method with other multireference methodologies will be investigated. The JM ansatz introduces the genealogy of the coefficients with respect to the Slater determinants within the CASCI space:
(5) 
where the quantity is the excitation amplitude related to the excitation process that leads from to :
(6) 
Here, we restrict to be a single or double excitation operator. Within this JM formulation of , a very general first order approximation of the amplitudes can be expressed as:
(7) 
where the excitation energy depends explicitly of the couple . Such a definition is different from other determinantbased MRPT2 like the CIPSI or shifted where the excitation energy does not depend on the couple but only on the . With this definition of , one can write the secondorder correction to the energy as:
(8)  
and the total secondorder energy :
(9)  
Definition of the energy denominators
The firstorder wave function can be written explicitly in terms of the excitation operators :
(10)  
However, one can notice that

the excitation operators do not explicitly depend on as they are general single or double excitation operators, just as in the Hamiltonian for instance;

a given excitation operator contributes to the coefficients of several ( = = );

the application of all the single and double excitation operators on each generates the entire set of as the reference is a CAS.
Therefore one can rewrite the firstorder perturbed wave function directly in term of excitation operators applied on the each CASCI Slater determinant as:
(11) 
where the is the part of the firstorder wave function associated with the excitation process :
(12) 
In order to fully define our perturbation theory and intermediate Hamiltonian theory, one needs to select an expression for the energy denominators occurring in the definition of . We propose to take a quantity that does not depend explicitly on the reference determinant but only depends on the excitation process :
(13) 
Consequently, in the expression of (see Eq. (12)), the energy denominator can be factorized:
(14)  
where is simply:
(15) 
Also, one can notice that, as and differ by a simple constant factor, they have the same normalized expectation values:
(16) 
Then, the excitation energy is simply taken as the difference of the normalized expectation values of the Dyall Hamiltonian over and :
(17)  
This ensures the strong separability when localized orbitals are used, as will be illustrated numerically in the section V.
The Dyall Hamiltonian is nothing but the exact Hamiltonian over the active orbitals, and a MøllerPlesset type operator over the doubly occupied and virtual orbitals. If one labels the active spinorbitals, the spinorbitals that are always occupied and the virtual spinorbitals, the Dyall Hamiltonian can be written explicitly as:
(18) 
(19) 
where the and are defined as the spinorbital energies associated with the density given by , and the effective active oneelectron operator . With a proper choice of the constant in Eq. (19),
(20) 
one has:
(21) 
Because the Dyall Hamiltonian acts differently on the active and inactivevirtual orbitals, the excitation energy is the sum of an excitation energy associated with the inactive and virtual orbitals and of an excitation energy associated with the active orbitals:
(22) 
Also, it is useful to differentiate the active part from the inactivevirtual part of the excitation :
(23) 
The inactivevirtual excitation energy is simply:
(24) 
where and refer to, respectively, the inactive and virtual spinorbitals involved in the excitation operator . Conversely, the active excitation energy has a more complex expression, namely:
(25) 
Practical consequences: the difference between single and double excitation operators
From Eq. (25), one must differentiate the class of the pure single excitation operators from the pure double excitation operators. For the sake of clarity, we define the spinadapted bielectronic integrals as:
(26) 
where is the spin variable of the spin orbitals and . If one considers a given double excitation involving four different spin orbitals , , and :
(27) 
one can notice that the Hamiltonian matrix elements associated with this double excitation only depend, up to a phase factor, on the four indices involved in the . Indeed, if is possible on both and , one has:
(28)  
and as
(29)  
it becomes:
(30) 
Therefore, as the hamiltonian matrix elements of type can be factorized both in the numerator and the dominator of the expression of the active part of the excitation energy (see Eq. (25)). Finally, the expression of the active part of the excitation energy for a given double excitation is simply:
(31)  
As a consequence, the amplitudes and associated with the same excitation for different parents and are also equal:
(32)  
and one can define a unique excitation operator which does not depend on the reference determinant on which it acts. The explicit form of the referenceindependent excitation operator is
(33) 
In the case where is a pure single excitation operator, the term may strongly depend on and Eq. (25) cannot be simplified.
Precaution for spin symmetry
As the formalism proposed here deals with Slater determinants, it cannot formally ensure to provide spin eigenfunctions. In order to ensure the invariance of the energy with the value of a given spin multiplicity, we introduced a slightly modified version of the Dyall Hamiltonian which does not consider:

any exchange terms in the Hamiltonian matrix elements when active orbitals are involved

any exchange terms involving two electrons of opposite spins (namely and )
ii.2 The JMHeffPT2 method
An advantage of a determinantbased multireference perturbation theory is that it can be easily written as a dressing of the Hamiltonian matrix within the reference space. Starting from the Schrödinger equation projected on a given reference determinant one has:
(34) 
Using the expression for the first order coefficients , it becomes:
(35) 
Therefore, one can define an nonHermitian operator :
(36) 
and a dressed Hamiltonian as:
(37) 
such that Eq. (35) becomes a nonsymmetric linear eigenvalue equation within the CASCI space:
(38) 
The secondorder correction to the energy can be simply obtained as the expectation value of over the zerothorder wave function:
(39)  
Finally, one can define a Hermitian operator :
(40) 
and a corresponding eigenpair (, ) verifying:
(41) 
The diagonalization of such a Hamiltonian allows then to improve the CASCI wave function by treating the coupling that can exist between the correlation effects within and outside the CASCI space.
Iii Links with other multireference methods
iii.1 Strongly contracted NEVPT2
It is interesting to understand the similarities and differences between the present JMMPRT2 and other strictly sizeconsistent MRPT2 methods, like the NEVPT2 and specially its strongly contracted variant (SCNEVPT2). The first important similarity is that both the JMMRPT2 and the NEVPT2 methods use the Dyall Hamiltonian. Then, the JMMRPT2 uses perturbers that are individual Slater determinants, whereas the NEVPT2 uses linear combinations of Slater determinants. However, in the SCNEVPT2, the contraction coefficients are closely related to the Hamiltonian matrix elements, just as in the JMMRPT2 method. In order to better understand the differences between the SCNEVPT2 and JMMRPT2, let us take a practical example. Here, are inactive spinorbitals, are active spinorbitals and are virtual spinorbitals. Considering a given semiactive double excitation , the firstorder amplitude associated with in the JMMRPT2 formalism is given by:
(42) 
where the active part of the excitation energy directly comes from Eq. (31):
(43) 
Note that such a quantity can be thought as an approximation of the electron affinity of the molecule, as it is the change in energy when one introduces “brutally” an electron in spin orbital without relaxing the wave function. Consequently, as it has been emphasized in section III.2, one can consider the part of the firstorder perturbed wave function generated by the excitation :
(44) 
which turns out to be:
(45)  
In the SCNEVPT2 framework, one does not consider explicitly a given but has to consider a unique excitation which is a linear combination of all possible for all active spin orbitals , with proper contraction coefficients. To be more precise, the firstorder perturbed wave function associated with is:
(46) 
where the excitation energy associated with is unique for all the excitation operators , and can be thought as an average excitation energy over all . Consequently, one can express the part of that comes from the as:
(47) 
which we can compare to Eq. (45) in the case of the JMMRPT2 method. Then, the only difference between the SCNEVPT2 and the JMMRPT2 is the definition of the excitation energy occurring in Eqs. (45) and (47). In the SCNEVPT2 method, the excitation energy is closely related to the excitation energy defined in JMMRPT2:
(48)  
where the quantity is the same for all active orbitals and defined as:
(49) 
Under this perspective, one sees that the quantity is related to defined in Eq. (43):

in the JMMRPT2 method, the quantity explicitly refers to the “brutal” addition of en electron in orbital , whatever the inactive orbitals or virtual orbitals involved in ;

the quantity involved in the SCNEVPT2 is an average electronic affinity over all possible excitation processes within the active space, but keeping a trace of the inactive and virtual excitation processes involved in thanks to the interaction .
Consequently, the quantity contains also the interactions between the various . To summarize, on one hand, the JMMRPT2 gives a different but rather crude excitation energy for each , and on the other hand the SCNEVPT2 has a unique and sophisticated excitation energy for all . Of course, one can extend this comparison to all the other classes of double excitations.
iii.2 Multireference coupled cluster methods
The present formalism has also several links with other multireference methods. First of all, as it uses a JM genealogical definition for the coefficients (see Eqs. (5) and (7)), the wave function corrected at first order can be written as:
(50)  
By introducing the excitation operator acting only on as:
(51) 
the expression of in Eq. (50) becomes:
(52) 
Such a parameterization for the firstorder corrected wave function recalls immediately a firstorder Taylor expansion of the general JMMRCC ansatz:
(53) 
Also, within the present formalism, the class of the double excitations can be factorized as shown in the previous section. Therefore, using the referenceindependent amplitudes defined in Eq. (33) one can define a unique double excitation operator as:
(54) 
recalling thus the formalism of the internallycontractedMRCC(50); (51); (52); (53) (icMRCC) which uses a unique excitation operator as in the singlereference coupledcluster:
(55)  
However, unlike the icMRCC formalism, the JMMRPT2 equations do not suffer from the linear dependency problems.(51) In such a perspective, as the energy provided by the JMMRPT2 equations is sizeextensive, it can be seen as a linearized coupled cluster version using a hybrid parameterization of the wave function: internally contracted ansatz for the double excitation operators and JM ansatz for the single excitation operators.
iii.3 Determinantbased multireference perturbation theories
The JMMRPT2 presented here can be directly compared to the CIPSI method, just as the JMHeffPT2 can be directly compared to the Shifted method. Indeed, by using the following amplitudes:
(56) 
in the equation of the secondorder correction on the energy (see Eq. (8)), one obtains the CIPSI energy, and by introducing in the definition of the dressed Hamiltonian , one obtains the Shifted energy. As mentioned previously, it has been shown that the sizeconsistency error of these methods comes from the unbalanced treatment between the variational energy of a multireference wave function such as and the variational energy of the single Slater determinant . Such an error is not present within the definitions of the excitation energies in the JMMRPT2 method as the latter introduces expectation values of the Hamiltonian over linear combinations of perturber Slater determinants.
In a similar context, one can compare the JMHeffPT2 method to the SplitGAS(44) of Li Manni et al whose definition of the amplitude is:
(57) 
In the SplitGAS framework, the correlation energy brought by the perturbers is included in the energy denominator, which introduces self consistent equations as in the BrillouinWigner perturbation theory.(45) However, the sizeconsistency error in such a method is even more severe than in the Shifted as the excitation energies are much larger due to the presence of the total correlation energy .
Iv Computational cost
iv.1 Mathematical complexity and memory requirements
Compared to other sizeextensive MRPT2 methods, a clear advantage of the JMMRPT2 is its simplicity. The NEVPT2 approach requires to handle the fourbody density matrix and the CASPT2 needs to handle the threebody density matrix. Both of these computationally intensive phases can be skipped in our formalism as one only needs to compute expectation values whose number are relatively small compared to NEVPT2 and CASPT2. The most involved quantity to be computed is
(58) 
for all pairs where is an inactive orbital and is a virtual orbital. These quantities need to be only computed once since they can all fit in memory. Each is, from the computational point of view, equivalent to an expectation value over the CASSCF wave function. As all the are independent, the computation of these quantities can be trivially parallelized. Regarding the memory footprint of the JMMRPT2 method, it scales as ( being the number of active orbitals) for the storage of the and quantities.
Regarding the complexity of the equations for the amplitudes, it is clear that once computed the active part of the denominator, the JMMRPT2 is just a simple sum of contributions. This is in contrast with the UGASSMRPT2 equations which involve the handling of coupled amplitude equations.
iv.2 Removal of the determinantbased computational cost
The present formalisms are formally determinantbased methods, which implies that the computational cost should be proportional to the number of perturbers that one has to generate to compute the corrections to the energy or the dressing of the Hamiltonian matrix, just as in the CIPSI, Shifted or UGASSMRPT2 methods. To understand the main computational costs, one can divide the excitation classes according to the difference dedicated CI (DDCI) framework,(46) which classifies the Slater determinants in terms of numbers of holes in the doubly occupied orbitals and particles in the virtual orbitals. If is the number of Slater determinants of the CASCI zeroth order wave function, , and respectively the number of doubly occupied, active and virtual orbitals, one can then classify each excitation class according to the number of perturbers needed to compute their contribution to the secondorder perturbation correction to the energy:

the twoholestwoparticles excitation class (2h2p) which scales as

the oneholetwoparticles excitation class (1h2p) which scales as

the twoholesoneparticle excitation class (2h1p) which scales as

the twoparticles excitation class (2p) which scales as

the twoholes excitation class (2h) which scales as

the oneholeoneparticle excitation class (1h1p) which scales as

the oneparticle excitation class (1p) which scales as

the onehole excitation class (1h) which scales as
Nevertheless, our formalism presents several mathematical simplifications that allow one to basically remove any browsing over the Slater determinants , and once more there is a difference between the single and double excitations processes.
iv.3 Factorization of the most numerous double excitation processes
As the five most computationally demanding excitation classes involve only double excitation operators in their equations, their contribution can be formalized directly thanks to the one and twobody density matrices of the zerothorder wave function. To understand how, one can write the secondorder correction to the energy as:
(59)  
Consequently, as is necessarily of type:
(60) 
one can reformulate the secondorder correction to the energy in terms of manybody density matrices:
(61)  
Such a formulation avoids completely to run over Slater determinants, and consequently kills the prefactor in involved in each of the excitation classes, just as in the internallycontracted formalisms. Of course, because of the restrictions in terms of holes and particles in the inactive and virtual orbitals, the handling of the fourbody density matrix never occurs in our formalism. We report here the explicit equations for the energetic corrections of the five most numerous double excitation classes: