A Jeziorski-Monkhorst fully uncontracted Multi-Reference perturbative treatment I: principles, second-order versions and tests on ground state potential energy curves

A Jeziorski-Monkhorst fully uncontracted Multi-Reference perturbative treatment I: principles, second-order versions and tests on ground state potential energy curves

Abstract

The present paper introduces a new multi-reference perturbation approach developed at second order, based on a Jeziorsky-Mokhorst 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 CAS-CI. Such a formulation accounts then for the coupling between the static and dynamic correlation effects. With our new definition of zeroth-order energies, these two approaches are strictly size-extensive 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 determinant-based perturbation theories. The accuracy of these methods has been investigated on ground-state 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 ground-state wave function of closed-shell molecules follows well-established paths. The perturbative expansions from the mean-field Hartree-Fock single determinant usually converge and may be used as basic tools, especially when adopting a mono-electronic zero-order Hamiltonian known as the Møller-Plesset Hamiltonian.(1) In this approach the wave function and the energy may be understood in terms of diagrams, which lead to the fundamental linked-cluster theorem.(2) The understanding of the size-consistency 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 zero-order wave function. These situations exhibit an intrinsic Multi-Reference (MR) character. A generalized linked-cluster theorem has been established by Brandow,(9) which gives a basis to the understanding of the size-consistency problem in this context, but the conditions for establishing this theorem are severe. They require a Complete Active Space (CAS) model space and a mono-electronic zero-order Hamiltonian. Consequently, the corresponding Quasi-Degenerate Perturbation Theory (QDPT) expansion cannot converge in most of the molecular MR situations(10); (11); (12). The research of theoretically satisfying (size-consistent) and numerically efficient MR treatments remains a very active field in Quantum Chemistry, as summarized in recent review articles concerning either perturbative(13) or Coupled-Cluster(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 zero-order 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 zero-order Hamiltonian. The simplest approach consists in using single determinants as outer-space 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 size-consistent and ii) it does not revise the model-space component of the wave function under the effect of its interaction with the outer-space. 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 state-specific(24) or multi-state(25) versions. Other methods which start from a CAS model space and use multi-determinantal outer-space vectors have been proposed later on and are now broadly used. The first one is the CASPT2 method,(26); (27) which employs a mono-electronic zero-order 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 size-consistent. The NEVPT2 method(28); (29); (30) also uses multi-determinantal perturbers (defined in two different ways in its partially and strongly contracted versions), it makes use of a more sophisticated bi-electronic Hamiltonian (the Dyall Hamiltonian(33)) to define the zero-order energies of these perturbers, it is parameter-free, intruder-state free, and size-consistent. 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). Multi-State 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 single-determinant perturbers, the origin of their size-inconsistency problem has been identified as due to the unbalance between the multi-determinant character of the zero-order 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 multi-reference linear coupled cluster equations in a stochastic framework of Full-CI Quantum Monte Carlo which also uses single Slater determinants as perturbers(39).

The present paper is composed as follows. In Section II, the here-proposed formalism is presented, whose main features are:

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

  2. the perturbers are single determinants (the method is externally non-contracted, according to the usual terminology);

  3. it is state-specific, and strictly separable when localized active MOs are used (see formal demonstration in the appendix);

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

  5. it is based on a Jeziorski-Monkhorst(40) (JM) expression of the wave operator and proceeds through reference-specific partitionings of the zero-order Hamiltonian, as it has been previously suggested in the so-called Multi-Partitionning(41); (42); (43) (MUPA) and also in the UGA-SSMPRT2.(38) Consequently, it does not define a unique zero-order energy to the outer-space determinants;

  6. it can be expressed either as a second-order energy correction or as a dressing of the CAS-CI matrix, which offers a full flexibility in the treatment of the feed-back effect of the post-CAS-CI correlation on the model space component of the wave function;

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

  8. thanks to our definition of the zeroth-order energies, all processes involving double excitations can be treated by using only the one- and two-boy density matrices, avoiding to loop on the perturbers;

  9. 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 multi-reference 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 JM-MRPT2 and JM-HeffPT2 methods. A numerical test of size-extensivity 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 JM-MRPT2 method.

Ii Working equations for the perturbation and effective Hamiltonian at second order

Figure 1: Example of interactions: the two determinants of the CAS interact through a bi-electronic operator involving the two active orbitals and , just as the two perturbers determinants generated by the same excitation operator on the two CAS determinants.

As demonstrated previously by one of the present authors and his collaborators,(34) the size-consistency problem in any multi-reference 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 CAS-CI 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 extra-diagonal 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 CAS-CI 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 bi-electronic 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 size-consistency 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 JM-MRPT2, 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 size-consistent 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 JM-HeffPT2 method.

ii.1 The JM-MRPT2 method

First-order perturbed wave function and second-order 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 CAS-CI 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 CAS-CI space. One should notice that such form is in principle not exact, as some changes of the coefficients within the CAS-CI space can formally occur when passing from the CAS-CI 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 CAS-CI 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 multi-reference methodologies will be investigated. The JM ansatz introduces the genealogy of the coefficients with respect to the Slater determinants within the CAS-CI 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 determinant-based 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 second-order correction to the energy as:

(8)

and the total second-order energy :

(9)

Definition of the energy denominators

The first-order wave function can be written explicitly in terms of the excitation operators :

(10)

However, one can notice that

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

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

  3. 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 first-order perturbed wave function directly in term of excitation operators applied on the each CAS-CI Slater determinant as:

(11)

where the is the part of the first-order 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øller-Plesset type operator over the doubly occupied and virtual orbitals. If one labels the active spin-orbitals, the spin-orbitals that are always occupied and the virtual spin-orbitals, the Dyall Hamiltonian can be written explicitly as:

(18)
(19)

where the and are defined as the spin-orbital energies associated with the density given by , and the effective active one-electron 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 inactive-virtual 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 inactive-virtual part of the excitation :

(23)

The inactive-virtual excitation energy is simply:

(24)

where and refer to, respectively, the inactive and virtual spin-orbitals 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 spin-adapted 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 reference-independent 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:

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

  2. any exchange terms involving two electrons of opposite spins (namely and )

ii.2 The JM-HeffPT2 method

An advantage of a determinant-based multi-reference 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 non-Hermitian operator :

(36)

and a dressed Hamiltonian as:

(37)

such that Eq. (35) becomes a non-symmetric linear eigenvalue equation within the CAS-CI space:

(38)

The second-order correction to the energy can be simply obtained as the expectation value of over the zeroth-order 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 CAS-CI wave function by treating the coupling that can exist between the correlation effects within and outside the CAS-CI space.

Iii Links with other multi-reference methods

iii.1 Strongly contracted NEVPT2

It is interesting to understand the similarities and differences between the present JM-MPRT2 and other strictly size-consistent MRPT2 methods, like the NEVPT2 and specially its strongly contracted variant (SC-NEVPT2). The first important similarity is that both the JM-MRPT2 and the NEVPT2 methods use the Dyall Hamiltonian. Then, the JM-MRPT2 uses perturbers that are individual Slater determinants, whereas the NEVPT2 uses linear combinations of Slater determinants. However, in the SC-NEVPT2, the contraction coefficients are closely related to the Hamiltonian matrix elements, just as in the JM-MRPT2 method. In order to better understand the differences between the SC-NEVPT2 and JM-MRPT2, let us take a practical example. Here, are inactive spin-orbitals, are active spin-orbitals and are virtual spin-orbitals. Considering a given semi-active double excitation , the first-order amplitude associated with in the JM-MRPT2 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 first-order perturbed wave function generated by the excitation :

(44)

which turns out to be:

(45)

In the SC-NEVPT2 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 first-order 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 JM-MRPT2 method. Then, the only difference between the SC-NEVPT2 and the JM-MRPT2 is the definition of the excitation energy occurring in Eqs. (45) and (47). In the SC-NEVPT2 method, the excitation energy is closely related to the excitation energy defined in JM-MRPT2:

(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 JM-MRPT2 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 SC-NEVPT2 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 JM-MRPT2 gives a different but rather crude excitation energy for each , and on the other hand the SC-NEVPT2 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 Multi-reference coupled cluster methods

The present formalism has also several links with other multi-reference 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 first-order corrected wave function recalls immediately a first-order Taylor expansion of the general JM-MRCC 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 reference-independent amplitudes defined in Eq. (33) one can define a unique double excitation operator as:

(54)

recalling thus the formalism of the internally-contracted-MRCC(50); (51); (52); (53) (ic-MRCC) which uses a unique excitation operator as in the single-reference coupled-cluster:

(55)

However, unlike the ic-MRCC formalism, the JM-MRPT2 equations do not suffer from the linear dependency problems.(51) In such a perspective, as the energy provided by the JM-MRPT2 equations is size-extensive, 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 Determinant-based multi-reference perturbation theories

The JM-MRPT2 presented here can be directly compared to the CIPSI method, just as the JM-HeffPT2 can be directly compared to the Shifted- method. Indeed, by using the following amplitudes:

(56)

in the equation of the second-order 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 size-consistency error of these methods comes from the unbalanced treatment between the variational energy of a multi-reference 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 JM-MRPT2 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 JM-HeffPT2 method to the Split-GAS(44) of Li Manni et al whose definition of the amplitude is:

(57)

In the Split-GAS framework, the correlation energy brought by the perturbers is included in the energy denominator, which introduces self consistent equations as in the Brillouin-Wigner perturbation theory.(45) However, the size-consistency 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 size-extensive MRPT2 methods, a clear advantage of the JM-MRPT2 is its simplicity. The NEVPT2 approach requires to handle the four-body density matrix and the CASPT2 needs to handle the three-body 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 JM-MRPT2 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 JM-MRPT2 is just a simple sum of contributions. This is in contrast with the UGA-SSMRPT2 equations which involve the handling of coupled amplitude equations.

iv.2 Removal of the determinant-based computational cost

The present formalisms are formally determinant-based 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 UGA-SSMRPT2 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 CAS-CI 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 second-order perturbation correction to the energy:

  1. the two-holes-two-particles excitation class (2h2p) which scales as

  2. the one-hole-two-particles excitation class (1h2p) which scales as

  3. the two-holes-one-particle excitation class (2h1p) which scales as

  4. the two-particles excitation class (2p) which scales as

  5. the two-holes excitation class (2h) which scales as

  6. the one-hole-one-particle excitation class (1h1p) which scales as

  7. the one-particle excitation class (1p) which scales as

  8. the one-hole 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 two-body density matrices of the zeroth-order wave function. To understand how, one can write the second-order correction to the energy as:

(59)

Consequently, as is necessarily of type:

(60)

one can reformulate the second-order correction to the energy in terms of many-body 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 internally-contracted formalisms. Of course, because of the restrictions in terms of holes and particles in the inactive and virtual orbitals, the handling of the four-body 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: