Alternative definition of excitation amplitudes in Multi-Reference state-specific Coupled Cluster

TL;DR

This paper introduces a new reference-independent operator for state-specific Multi-Reference Coupled Cluster methods, improving the definition of excitation amplitudes and enhancing accuracy in molecular energy calculations.

Contribution

It proposes a novel operator to define excitation amplitudes in MR-CC, enabling better approximation and iterative refinement for more accurate molecular energies.

Findings

Achieved ~1 m$E_{ m h}$ non-parallelism error with FCI curves.

Tested on 6 molecules including triple-bond breaking states.

Results show similar performance for both variants of the method.

Abstract

A central difficulty of state-specific Multi-Reference Coupled Cluster (MR-CC) formalisms concerns the definition of the amplitudes of the single and double excitation operators appearing in the exponential wave operator. If the reference space is a complete active space (CAS) the number of these amplitudes is larger than the number of singly and doubly excited determinants on which one may project the eigenequation, and one must impose additional conditions. The present work first defines a state-specific reference-independent operator $\hat{\tilde{T}}^m$ which acting on the CAS component of the wave function $|\Psi_0^m \rangle$ maximizes the overlap between $(1+\hat{\tilde{T}}^m)|\Psi_0^m \rangle$ and the eigenvector of the CAS-SD CI matrix $|\Psi_{\rm CAS-SD}^m \rangle$. This operator may be used to generate approximate coefficients of the Triples and Quadruples, and a dressing of the CAS-SD CI matrix, according to the intermediate Hamiltonian formalism. The process may be iterated to convergence. As a refinement towards a strict Coupled Cluster formalism, one may exploit reference-independent amplitudes provided by $(1+\hat{\tilde{T}}^m)|\Psi_0^m \rangle$ to define a reference-dependent operator $\hat{T}^m$ by fitting the eigenvector of the (dressed) CAS-SD CI matrix. The two variants, which are internally uncontracted, give rather similar results. The new MR-CC version has been tested on the ground state potential energy curves of 6 molecules (up to triple-bond breaking) and a two excited states. The non-parallelism error with respect to the Full-CI curves is of the order of 1 m$E_{\rm h}$.

Tables

Table 1. Table 1: Non-parallelism errors (NPE) and maximum errors with respect to the Full-CI potential energy surface (m E h subscript 𝐸 h E_{\rm h} )

	CAS-SD		$\lambda$-MR-CCSD		Dressed CAS-SD		$\mu$-MR-CCSD
	NPE	Max Error	NPE	Max Error	NPE	Max Error	NPE	Max Error
C2H6	5.1	35.5	3.6	8.4	1.3	8.3	3.5	8.3
F2	3.8	19.8	3.1	4.0	0.9	4.2	1.6	3.9
C2H4 twist	1.5	27.7	0.3	6.7	0.7	7.2	0.5	6.7
BeH2	2.9	4.1	1.3	1.8	2.0	2.4	1.8	2.2
H2O	1.8	4.6	0.7	1.2	0.2	1.3	0.5	1.2
C2H4 stretch	2.8	22.1	1.6	5.3	1.7	6.2	1.6	5.2
N2	1.8	9.0	1.0	2.2	1.7	3.9	0.7	2.8
F2 ${}_{}{}^3\Sigma _u^{+}\left(m_s=1\right)$	2.5	18.6	1.3	3.3	1.3	3.5	1.2	3.3
F2 ${}_{}{}^3\Sigma _u^{+}\left(m_s=0\right)$	2.5	18.6	1.2	1.8	1.3	3.5	1.1	3.3
HF (ground state)	2.6	14.6	1.8	3.4	2.1	4.4	1.8	4.1
HF (excited state)	3.3	20.9	8.8	8.5	10.5	10.1	7.1	8.3
F2 (local)	3.8	19.8	1.2	3.2	1.5	3.1	1.0	4.6
N2 (local)	1.8	9.0	3.8	5.0	1.1	3.5	1.1	2.8

Equations

$$\Psi =\hat{\Omega }{\Phi }_0=\exp \left(\hat{T}{\Phi }_0\right)$$

$$|{\Psi }_{\mathrm{CAS}-\mathrm{SD}}^m⟩$$

$${\hat{T}}^m=\sum _I{\hat{T}}_I^m|I⟩⟨I|$$

$$|{\Psi }_{\mathrm{CAS}-\mathrm{SD}}⟩=\sum _I{C}_I|I⟩+\sum _i{c}_i|i⟩$$

$$\hat{\Omega }=\exp (\hat{T})$$

$$|{\Psi }_{\mathrm{CAS}-\mathrm{SD}}⟩=\sum _I{C}_I{\hat{T}}_I|I⟩$$

$${\hat{T}}_I=\sum _i{t}_{I\to i}{\hat{T}}_{I\to i}$$

$$C_i=\sum _I{t}_{I\to i}{C}_I$$

$$\frac{t_{I\to i}}{t_{J\to i}}=\frac{⟨i|\hat{H}|I⟩}{⟨i|\hat{H}|J⟩}.$$

$$t_{I\to i}={\lambda }_i⟨i|\hat{H}|I⟩$$

$${\lambda }_i=\frac{c_i}{⟨i|\hat{H}|{\Psi }_0⟩}$$

$$\hat{T}$$

$$=$$

$$\begin{array}{r}\underset{t}{\mathrm{argmin}}\Vert (1+\hat{T})|{\Psi }_0⟩-|{\Psi }_{\mathrm{CAS}-\mathrm{SD}}⟩\Vert \\ =\underset{t}{\mathrm{argmin}}\Vert \hat{T}|{\Psi }_0⟩-|{\Psi }_{\mathrm{SD}}⟩\Vert ,\end{array}$$

$$({\mathbf{A}}^{†}\mathbf{A})\mathbf{t}={\mathbf{A}}^{†}\mathbf{c}$$

$$\underset{t_{jk\to rs}}{\mathrm{argmin}}\Vert {\hat{T}}_{jk\to rs}|{\Psi }_0⟩t_{jk\to rs}-|{\Psi }_{\mathrm{SD}}⟩\Vert$$

$$\underset{t_{jk\to rs}}{\min }{\left(\sum _I{C}_I{t}_{jk\to rs}-\sum _i{c}_i⟨i|{\hat{T}}_{jk\to rs}|I⟩\right)}^2$$

$$t_{jk\to rs}=\frac{{\sum }_I{C}_I{c}_i}{{\sum }_I{C}_I^2}$$

$$t_{jk\to rs}=\frac{1}{{\sum }_I{C}_I^2}\left(\sum _I{C}_I^2\left(\frac{c_i}{C_I}\right)\right)$$

$$\left\{\begin{array}[]{l}{\bf t}_{0}={{\bf A}}^{\dagger}{\bf c}\\ {\bf t}_{n+1}={{\bf A}}^{\dagger}{\bf c}+\left({\bf I}-{{\bf A}}^{\dagger}{\bf A}\right){\bf t}_{n}\\ \end{array}\right.$$

$$\Psi =\hat{\Omega }{\Psi }_0$$

$$\hat{\Omega }=\exp (\hat{T})$$

$$c_{\alpha }=\frac{1}{2}⟨\alpha |{\hat{T}}^2|{\Psi }_0⟩.$$

$$\begin{array}{r}\left(⟨i|\hat{H}|i⟩-E\right)c_i+\sum _J⟨i|\hat{H}|J⟩C_J+\\ \sum _j⟨i|\hat{H}|j⟩c_j+\sum _{\alpha }⟨i|\hat{H}|\alpha ⟩c_{\alpha }=0.\end{array}$$

$$\langle i|\Delta|I\rangle=\sum_{\alpha}\langle i|\hat{H}|\alpha\rangle\left(\sum_{\left\{(p,q)\big{|}\hat{T}_{p}\hat{T}_{q}|I\rangle=|\alpha\rangle\right\}}t_{p}t_{q}\right)$$

$$\sum _I⟨i|\Delta |I⟩C_I=\sum _{\alpha }⟨i|\hat{H}|\alpha ⟩c_{\alpha }.$$

$$\tilde{c}_{i}=\langle i|\hat{T}|\Psi_{0}\rangle=\sum_{\left\{(I,l)\big{|}\hat{T}_{l}|I\rangle=|i\rangle\right\}}t_{l}\,C_{I}$$

$${\mu }_i=\frac{c_i}{{\tilde{c}}_i}$$

$$t_{I\to i}=t_{l,I}={\mu }_i{t}_l.$$

$$z=2.54-0.46x\text{(a.u.)}$$

Full text

I Introduction

The single-reference Coupled Cluster (CC) formalismCoester (1958); Coester and Kümmel (1960); Čiźěk (1966); Bartlett et al. (1990) is the standard technique in the study of the ground state of closed-shell molecules, i.e. those for which a mean-field treatment provides a reasonable zero-order single-determinant wave-function $\Phi_{0}$. This method incorporates the leading contributions to the correlation energy in a given basis set, it is based on the linked-cluster theoremGoldstone (1957) and is size-consistent since it is free from unlinked contributions. The method generates an approximate wave function under the action of a wave operator $\hat{\Omega}$ acting on the single-determinant reference $\Phi_{0}$, and assumes an exponential character to the wave operator

[TABLE]

The most popular version only introduces single and double excitation operators in $\hat{T}$, and is known as the Coupled Cluster Singles and Doubles (CCSD) approximation. It incorporates the fourth order correction of the quadruply excited determinants. The lacking fourth order contribution concerns the triply excited determinants, which may be added in a perturbative manner. The CC equations, obtained by projecting the eigenequation on each of the Singles and Doubles lead to coupled quadratic equations. In practice, guess values of the amplitudes of the $\hat{T}_{0\rightarrow i}$ operators appearing in the $\hat{T}$ operator may be taken as the coefficients of the Singles and Doubles $|i\rangle$ in the intermediate normalization of the SD CI vector. The solution of the CC equations may be obtained by treating the effect of the Triples and Quadruples as an iterative dressing of the SD CI matrix,Nebot-Gil et al. (1995) according to the Intermediate Effective Hamiltonian (IEH) theory.Malrieu, Durand, and Daudey (1985); Kirtman (1981) The field of application of this method, which is both intellectually satisfying and numerically efficient, is however limited to the systems and the situations where a single-determinant zero-order description is relevant. This is no longer the case when chemical bonds are broken, creating open shells, as occurs in most of the chemical reactions. The magnetic systems generally present several open shells and the low spin-multiplicity states are intrinsically of multiple-determinant character. Due to near degeneracies, most of the excited states are not only of multi-determinantal but of multi-configurational character. The conception of a multi-reference counterpart of the CCSD formalism is highly desirable, and has been the subject of intense research. The most comprehensive review has been given by R. Bartlett and his colleagues.Lyakh et al. (2012) For formal reasons and in particular to treat correctly the breaking of bonds, the reference space, or model space, is usually taken as a Complete Active Space (CAS), i.e. the Full-CI of a well-defined number of electrons (the active electrons) in a well-defined set of orbitals (the active MOs). The other MOs are called inactive. Let us label $|I\rangle,|J\rangle,\dots$ the reference determinants. The determinants $|i\rangle,|j\rangle,\dots$ which interact with the reference space are obtained under purely inactive or semi-active single and double excitations, they generate the CAS-SD CI space, the diagonalization of which provides a size-inconsistent energy $E_{\rm CAS-SD}^{m}$ and the corresponding eigenvector,

[TABLE]

with $\langle\Psi_{\rm CAS-SD}^{m}|\Psi_{\rm CAS-SD}^{m}\rangle=1$.

One strategy, which is not very aesthetic since it breaks the symmetry between degenerate reference determinants, but which has given rather satisfactory results, consists in selecting (eventually in an arbitrary manner) a specific single reference and in introducing in the wave operator the multiple excitations which generate the other references (the other determinants of the model space). The other strategies consider all the references on an equal footing, and are really multi-reference. Let us call $N$ the number of references, and $n$ the number of SD determinants. If the treatment pretended to provide $N$ eigenvectors simultaneously, one might define the $N\times n$ amplitudes sending from the references to the outer-space determinants, in a unique manner but this state-universal approach is not practicable when the model space is a CAS.

Most of the proposed formalisms are state-specific. In this case one faces the famous multi-parentage problem. This problem is recalled in section 2A. Sufficiency conditions have to be imposed.Mahapatra, Datta, and Mukherjee (1998) One solution was proposed by Mukherjee and coworkers, and has been widely tested.Das, Mukherjee, and Kállay (2010); Szabados (2011); Mahapatra, Datta, and Mukherjee (1999) Another one had been proposed earlier by one of us (JPM) and coworkers.Meller, Malrieu, and Caballol (1996) It consists, for a given outer-space determinant, in scaling the amplitudes of the various excitation operators $\hat{T}_{I\rightarrow i}$ on the interaction between the outer-space determinant and its parents. A recent work has implemented this second solution of the state-specific MR-CC problem and has tested its accuracy and robustness on a series of molecular benchmarks, comparing its results to the Full-CI (FCI) energies.Giner et al. (2016) The present work proposes an alternative process to define the amplitudes of the excitation operators. The state-specific MR-CC formalisms are usually based on the Jeziorski-MonkhorstJeziorski and Monkhorst (1981) splitting of the wave operator into a sum of operators acting individually on the various references

[TABLE]

We shall leave in a first time this assumption and define in section II.2 a reference-independent operator $\hat{T}$ which acting on the component of the desired state in the model space, $|\Psi_{0}^{m}\rangle$, provides a vector as close as possible to the CAS-SD eigenvector. This solution, defining reference-independent amplitudes of the excitations, may be exploited directly to generate approximate values of the coefficients of the triply and quadruply excited determinants, according to the exponential structure of the wave operator. From these coefficients one may dress the CAS-SD CI matrix, redefine amplitudes and iterate the process to convergence. This solution, presented in section 2C, is not an MR-CC technique, one may call it an exponential dressing of the CAS-SD CI matrix. Section 2D redefines reference-dependent excitation amplitudes from the reference-independent amplitudes by a fitting of the previous amplitudes on the coefficients of the Singles and Doubles of the (dressed) CAS-SD CI eigenvector. This represents an alternative solution to multi-parentage problem and opens the way to a strict MR-CC formalism. Section 3 presents a series of numerical tests on the bond breaking of single, double and triple bonds in ground states of molecules as well as a few tests on excited states. The results are compared to our previous proposal and with full Configuration Interaction (FCI) results.

II Formalisms

In this section, all the presented formalisms are state-specific. To simplify the notations we will consider that the state superscript $m$ is implicit for the wave functions ($\Psi^{m}\rightarrow\Psi$) and for the excitation operators ($\hat{T}^{m}\rightarrow\hat{T}$).

II.1 The multi-parentage problem in the Jeziorski-Monkhorst approach

Since one wants to produce a MR-CCSD method, one may start from a preliminary CAS-SD CI calculation which will help to fix guess values of the amplitudes of the excitation operators. Let us call $|I\rangle,|J\rangle,\dots$ the determinants of the CAS, i.e. the so-called reference vectors, and $|i\rangle,|j\rangle,\dots$ the Singles and Doubles which do not belong to the CAS and interact with them. The resulting approximate wave function of the targeted state $|\Psi\rangle$ is written

[TABLE]

Although this function is not size consistent one may note that the coefficients on the CAS determinants are no longer those of the CAS-CI : they incorporate the effect of the dynamical correlation on the composition of the CAS component of the wave function.

In CC formalisms the wave operator $\hat{\Omega}$ is assumed to take an exponential form

[TABLE]

and in our previous MR-CC formalismGiner et al. (2016) the Jeziorski-Monkhorst structure of the wave operator was adopted, introducing reference-specific wave operators acting specifically on each reference vector (Eq.3). One may exploit the knowledge of the CAS-SD CI eigenvector to determine guess operators $\hat{T}_{I}$ defined in such a manner that

[TABLE]

The $\hat{T}_{I}$ operators are a sum of single and double excitations $\hat{T}_{I\rightarrow i}$ possible on $|I\rangle$, multiplied by an amplitude $t_{I\rightarrow i}$

[TABLE]

In the single-reference CC the amplitudes of the excitation operators are obtained by projecting the eigenequation on the singly and doubly excited determinants, the number of unknowns is equal to the number of equations. This is no longer the case in the MR context : projecting the eigenequation on each on the singly or doubly excited vectors $|i\rangle$ is not sufficient to define the amplitudes $t_{I\rightarrow i}$ since for many classes of excitation an outer-space determinant interacts with several references, $|i\rangle=\hat{T}_{I\rightarrow i}|I\rangle=\hat{T}_{J\rightarrow i}|J\rangle$. The condition

[TABLE]

is not sufficient to define the amplitudes, even if one restricts the excitation operators to single and double excitations. Additional constraints have to be introduced to fix the amplitudes, and this is the famous multi-parentage problem. The number of amplitudes is larger than the number of outer-space determinants so that one cannot determine directly guess values of the amplitudes from Eq. 6. Different additional constraints have been proposed. One of them consists in scaling the amplitudes on the Hamiltonian interactions between the references and the outer space determinants,

[TABLE]

This constraint is expressed as

[TABLE]

where

[TABLE]

This solution has been recently implementedGiner et al. (2016) and shown to provide excellent agreements with Full-CI results on a series of molecular problems. It only presents minor stability problems in comparison with the present suggestion when the term ${\langle i|\hat{H}|\Psi_{0}\rangle}$ is small. From now on, we will refer to this method as the $\lambda$-MR-CCSD.

II.2 Introduction of reference-independent amplitudes

The present formalism will leave in the first step the Jeziorski-Monkhorst formulation of the wave operator and will consider the possibility to define a unique state-specific reference-independent operator $\hat{T}$, written as a sum of single and double excitation operators,

[TABLE]

where the indices $p$ and $q$ run on the virtual and active MOs and the indices $m$ and $n$ run on the inactive occupied and active MOs, excluding the possible occurrence of 4 active MOs.

This operator has the same form as the one introduced by the internally-contracted MR-CC (ic-MRCC) method by Evangelista and Gauss,Evangelista and Gauss (2011) and by Hanauer and Köhn,Hanauer and Köhn (2012) but it differs by both its determination and by the way we use it, as will appear later. The ic-MRCC method determines the amplitudes of the excitation by solving the projected Coupled Cluster equations, where the amplitudes appear as linear and quadratic terms. Hereafter we exploit the knowledge of the CAS-SD CI eigenvector to determine guess values of the reference-independent amplitudes. These excitation amplitudes will be used later on to estimate the coefficients of the Triples and Quadruples, and perform an iterative dressing of the CAS-SD CI matrix introducing the coupling between the Singles and Doubles with the Triples and Quadruples.

We propose a criterion to fix the amplitudes $t=\{t_{mn\rightarrow pq},t_{m\rightarrow p}\}$. Given the fact that we have at our disposal the CAS-SD wave function, a natural way to solve this overdetermined problem is to minimize the distance between the CAS-SD vector and the vector obtained by applying the $(1+\hat{T})$ operator on the CAS wave function

[TABLE]

$\hat{T}|\Psi_{0}\rangle$ being normalized such that $\|\hat{T}|\Psi_{0}\rangle\|=\||\Psi_{\rm SD}\rangle\|$.

To perform the minimization, we build the $N_{\rm SD}\times N_{\rm t}$ transformation matrix $A_{i,mn\rightarrow pq}=\langle i|\hat{T}_{mn\rightarrow pq}|\Psi_{0}\rangle$ which maps from the outer space of determinants $\{|i\rangle\}$ to the space of excited wave functions $\{\hat{T}_{mn\rightarrow pq}|\Psi_{0}\rangle\}$, and we search for the vector of amplitudes $\bf{t}$ which minimizes $\|{\bf A}.{\bf t}-{\bf c}\|$ by solving the normal equations

[TABLE]

Note that in the single-reference case, ${\bf A}$ is a permutation matrix and the CAS-SD wave function is exactly recovered.

The matrix ${\bf A}$ is usually so large that the use of standard singular value decomposition (SVD) routines to obtain the least squares solution is prohibitive.

Let us first consider the most numerous 2-hole-2-particle inactive double excitations $\hat{T}_{jk\rightarrow rs}$. Acting on a determinant $|I\rangle$ the operator creates a determinant $|i\rangle=\hat{T}_{jk\rightarrow rs}|I\rangle$ which can only be produced by this process. Therefore, the corresponding rows of ${\bf A}$ contain only one non-zero element located in the $jk\rightarrow rs$ column with value $A_{i,jk\rightarrow rs}=C_{I}$. The condition fixing the amplitude $t_{jk\rightarrow rs}$ is given by

[TABLE]

which is obtained by minimizing

[TABLE]

This condition turns out to be satisfied for

[TABLE]

One may notice that this is the weighted average of the ratios between the coefficients of the doubly excited determinants $|i\rangle$ and the coefficient of their unique reference generator,

[TABLE]

For all the remaining active excitations, A remains sparse since the maximum number of non-zero elements per column is equal to the number of reference determinants. Hence, we use Richardson’s iterative procedureRichardson (1911)

[TABLE]

which may be implemented very efficiently using sparse matrix products.

There are cases where multiple amplitudes applied to different references lead to same determinant : $\hat{T}_{jk\rightarrow rs}|I\rangle=\hat{T}_{lm\rightarrow tv}|J\rangle=|i\rangle$. If this determinant can be reached by no other process, there is an infinity of solutions for the amplitudes. The solution of Eq.(15) is ${\bf A}^{+}{\bf c}$ where ${\bf A}^{+}$ is the pseudo-inverse of ${\bf A}$, so the solution obtained minimizes the norm of the amplitude vector.Golub and Kahan (1965) In this way, the arbitrariness brought by the null space of ${\bf A}$ is minimized in the amplitude vector, and one obtains the most sensible solution.

II.3 Evaluation of the coefficients of Triples and Quadruples and iterative dressing of the CAS-SD CI matrix

The so-determined excitation operator $\hat{T}$ may be used to generate approximate values of the coefficient of the Triples and Quadruples as obtained by the action of $\frac{1}{2}\hat{T}^{2}$. Actually one may assume, in the spirit of the internally-contracted MR-CC methods, that the wave operator $\hat{\Omega}$ generating the correlated wave function $\Psi$ from $\Psi_{0}$,

[TABLE]

has an exponential structure,

[TABLE]

But this form will be simply used to estimate the coefficients of the triply and quadruply excited determinants $\{|\alpha\rangle\}$, leaving the internally-contracted structure of the outer-space. The coefficients of these determinants are estimated as

[TABLE]

In practice all the determinants $\{|\alpha\rangle\}$ are generated. For each $|\alpha\rangle$ one finds the reference determinants $\{|I\rangle\}_{\alpha}$ which differ by at most 4 orbital substitutions from $|\alpha\rangle$ (its grand-parents). One then identifies the set of complementary excitations \left\{(p,q)\,\big{|}\,\hat{T}_{p}\hat{T}_{q}|I\rangle=|\alpha\rangle\right\} as the products which generate $|\alpha\rangle$ from $\{|I\rangle\}_{\alpha}$. The genealogy of $|\alpha\rangle$ contributes to its associated coefficient $c_{\alpha}$ by the quantity $t_{p}t_{q}C_{I}$. Knowing $|\alpha\rangle$, one also knows the Singles and Doubles $\{|i\rangle\}_{\alpha}$ with which it interacts through the matrix elements $\langle i|\hat{H}|\alpha\rangle$, and in the eigenequation relative to $|i\rangle$

[TABLE]

One may replace the last sum by a dressing of the matrix elements between the determinant $|i\rangle$ and the references which are grand-parents of $|\alpha\rangle$,

[TABLE]

since

[TABLE]

The effect of the Triples and Quadruples is incorporated as a change of the columns of the CAS-SD CI matrix concerning the interaction between the references and the Singles and Doubles. This type of dressing was already employed in our previous MR-CC implementation.Giner et al. (2016) One will find in the same reference the practical procedure to make the dressed matrix Hermitian without any loss of information. Of course the whole process may be iterated. The diagonalization of the dressed CAS-SD CI matrix provides new values of the coefficients, not only of the Singles and Doubles which no longer suffer from the truncation, but also those of the references : the method is fully non-contracted. From the new wave function new amplitudes are obtained, a new dressing is defined and the process is repeated till convergence, which is usually rapidly obtained (3-4 iterations).

This formalism is not a strict MR-CC method since we exploit the CAS-SD CI function, and since this function slightly differs from the vector resulting from the action of $\hat{T}$ on the vector. Although the distance between these two vectors has been minimized they are not identical, $(1+\hat{T})|\Psi_{0}\rangle\neq|\tilde{\Psi}_{\rm CAS-SD}\rangle$.

Once the $\hat{T}$ operator has been obtained one might imagine a contracted exponential formalism calculating $\hat{T}^{2}|\Psi_{0}\rangle$ and the interaction between $\hat{T}|\Psi_{0}\rangle$ and $\hat{T}^{2}|\Psi_{0}\rangle$, but this calculation requires to return to the determinants. This formalism would remain internally contracted and would be less accurate than the procedure we propose. Actually in this version the deviations of the approximate reference-independent amplitudes from optimal ones, those which would generate the exact coefficients of the Singles and Doubles, only affects the evaluation of the coefficients of the Triples and Quadruples, and these deviations represent a minor source of error in the correction restoring the size extensivity. This reliability will be illustrated in the numerical tests.

II.4 State-specific MR-CC variant

In order to return to a MR-CC formalism, one may simply exploit the reference-independent amplitudes as an initial guess to define reference-dependent amplitudes. Currently the determinant $|i\rangle$ belonging to the Singles and Doubles has a coefficient $\tilde{c}_{i}$ in $\hat{T}|\Psi_{0}\rangle$

[TABLE]

which differs from the coefficient $c_{i}$ in $|\Psi_{\rm SD}\rangle$. One can define a parameter $\mu_{i}$, specific of the determinant $|i\rangle$,

[TABLE]

which multiplying $\tilde{c}_{i}$ will produce the exact coefficient $c_{i}$ of $|i\rangle$ in the (dressed) CAS-SD CI eigenvector. So the previous reference-independent amplitudes have now become reference-dependent. The excitation $\hat{T}_{l}$ which excites $|I\rangle$ to $|i\rangle$ ($|i\rangle=\hat{T}_{l}|I\rangle$) receives a reference-dependent amplitude

[TABLE]

The same excitation will receive a somewhat different amplitude when it acts on another reference $t_{l,J}\neq t_{l,I}$.

As the overlap between $(1+\hat{T})|\Psi_{0}\rangle$ and $|\Psi_{\rm CAS-SD}\rangle$ has been maximized the coefficients $\tilde{c}_{i}$ and $c_{i}$ are expected to be very close in particular if $c_{i}$ is large, and the parameter $\mu_{i}$ should be close to 1, at least for the determinants which contribute significantly to the wave function. In practice we observe this tendency, but the smallest coefficients are sacrificed during the maximization of the overlap and their $\mu_{i}$ can be very far from 1. This introduces some instabilities in the iterations, so we chose to limit the values of $\mu_{i}$ in the $[-2,2]$ range. The effect on the stability of the iterations is significant, and the effect on the energy is not noticeable.

This version returns to the Jeziorski-Monkhorst formalism as the wave operator again is a sum of reference-specific operators. The so-obtained amplitudes may be exploited to generate the coefficients of the Triples and Quadruples, and one may follow the same strategy as in our previous formalism, with an iterative column dressing of the interactions between the Singles and Doubles and the references. In a strict Coupled Cluster formalism, one should redefine the amplitudes of the double excitations by subtracting the products of the complementary single excitations they contain. For the sake of simplicity, as we did in our previous work, we did not proceed to this revision, the difference concerning only fifth-order perturbative corrections.

In what follows, we will refer to this method as $\mu$-MR-CCSD as it involves the $\mu_{i}$ (Eq.28).

III Numerical tests

In this section, we compare the here-proposed dressed CAS-SD and MR-CCSD to the MR-CCSD presented in ref Giner et al. (2016) on standard benchmark systems.Engels-Putzka and Hanrath (2009); Purvis and Bartlett (1982); Purvis et al. (1983); Laidig and Bartlett (1984); Mahapatra, Datta, and Mukherjee (1999); Kállay, Szalay, and Surján (2002); Hanrath (2005); Evangelista, Allen, and Schaefer (2006); Evangelista et al. (2010); Hanauer and Köhn (2011); Aoto and Köhn (2016); Hanauer and Köhn (2012); Das, Mukherjee, and Kállay (2010); Evangelista and Gauss (2011) To differentiate those two variants, we will label $\lambda$-MR-CCSD the variant of ref Giner et al. (2016) (Eq. 11) and the MR-CCSD of this work will be labeled $\mu$-MR-CCSD (Eq. 28).

The basis set used is Dunning’s cc-pVDZ,Dunning (1989) and the molecular orbitals were obtained using the CAS-SCF code present in GAMESS.Schmidt et al. (1993) All the following calculations were made using the Quantum Package,Scemama et al. (2016) an open-source library developed in our group. Full-CI energies were obtained using the CIPSI algorithm,Huron, Rancurel, and Malrieu (1973); Evangelisti, Daudey, and Malrieu (1983); Giner, Scemama, and Caffarel (2013) and the accuracy of the total energies is estimated to be of the order of $10^{-5}E_{\rm h}$. In all the calculations (Full-CI, CAS-SD and MR-CC), only the valence electrons are correlated (frozen core approximation).

III.1 Bond breaking

For all the applications we compare the dressed CAS-SD and $\mu$-MR-CCSD with the $\lambda$-MR-CCSD and the CAS-SD values. Results are also given using the reference-independent dressing of the CAS-SD CI matrix. All the applications are presented as energy differences with respect to the Full-CI energy estimated by a CIPSI calculation with a second-order perturbative correction. Figure 1 shows the difference of energy with respect to the Full-CI along the reaction coordinate. Table 1 summarizes the non-parallelism errors (NPE) and the maximum of the error obtained along the curve. The MR-CC treatment reduces the average and maximum error of the CAS-SD with respect to Full-CI by a factor close to 4. The correction is larger when the system involves an important number of inactive electrons (F2, C2H6) than when this number is small (BeH2, N2). One actually knows that the size-consistency error of the CAS-SD treatment increases with the number of inactive electrons, this error disappears in the MRCC treatment, which essentially misses some fourth-order connected effects of the Triples.

Single-bond breaking

We present here the single bond breaking of the $\sigma$ bonds of C2H6 and F2 molecules and of the $\pi$ bond of ethylene. The active spaces were chosen with two electrons in two MOs, the minimum wavefunctions to describe properly the dissociation of the molecules. In the case of ethane, the NPE of the CAS-SD is 5.1 m$E_{\rm h}$, and is reduced to 3.5 m$E_{\rm h}$ with the $\mu$-MR-CCSD. The curve of the dressed CAS-SD has the lowest NPE (1.3 m$E_{\rm h}$). The curves obtained by both MR-CCSD methods give equivalent results, with NPEs of 3.5 and 3.6 m$E_{\rm h}$.

In the case of F2 the NPE of the dressed CAS-SD is 0.9 m$E_{\rm h}$ and the NPE of the $\mu$-MR-CCSD is 1.6 m$E_{\rm h}$, both better than the NPE of the $\lambda$-MR-CCSD which has an NPE of 3.1 m$E_{\rm h}$. Also, one can remark here some numerical instabilities in the $\lambda$-MR-CCSD where the curve is not perfectly smooth.

In the next example, the $\pi$ bond of ethylene is broken by the rotation of the CH2 fragments. The CAS-SD has an NPE of 1.5 m$E_{\rm h}$, and using the dressed CAS-SD reduces the NPE to 0.7 m$E_{\rm h}$. The $\mu$-MR-CCSD gives an NPEs of 0.5 m$E_{\rm h}$, and the NPE obtained with the $\lambda$-MR-CCSD is slightly better with an NPE of 0.3 m$E_{\rm h}$.

Insertion of Be in H2

We present the results obtained by the insertion of a beryllium atom into the H2 molecule, which is a popular benchmark for MR-CC methods. The reference is still a CAS(2,2) for comparison with the literature, even though this choice of reference is not the most appropriate for a correct description of the reaction. The geometries are given by the relation

[TABLE]

where the beryllium atom is at the origin and the hydrogen atoms are at the coordinates $(x,0,\pm z)$. In this particular case, the $\mu$-MR-CCSD gives a NPE of 1.8 m$E_{\rm h}$ which is larger than the NPE of 1.3 m$E_{\rm h}$ obtained by the $\lambda$-MR-CCSD. This is due to only one point of the curve, the maximum which is higher by 0.4 m$E_{\rm h}$, all the other points being very close by less than 0.1 m$E_{\rm h}$. Here, the dressed CAS-SD and the $\mu$-MR-CCSD are equivalent.

Two bond breaking

For breaking two bonds we have used CAS(4,4) wave functions as the reference space. The first example is the simultaneous breaking of the two O—H bonds of the water molecule by stretching. Here, the CAS-SD exhibits a NPE of 1.8 m$E_{\rm h}$ which is significantly improved to 0.2 m$E_{\rm h}$ with the dressed CAS-SD. The $\mu$-MR-CCSD, with an NPE of 0.5 m$E_{\rm h}$, is slightly more parallel to the Full-CI curve than the $\lambda$-MR-CCSD which has an NPE of 0.7 m$E_{\rm h}$.

The second example is the double-bond breaking of ethylene by stretching. One should first clarify that the energy differences in the figure do not match those of the torsion along the bond because in the former example the reference was a CAS(2,2), and here it is a CAS(4,4). Dressing the CAS-SD reduces the NPE from 2.8 m$E_{\rm h}$ to 1.7 m$E_{\rm h}$. One can remark a discontinuity in the curve at large distances. The $\mu$-MR-CCSD and $\lambda$-MR-CCSD slightly improve the NPE to a value of 1.6 m$E_{\rm h}$, and both variants of the MR-CCSD are equivalent with smooth curves.

Triple-bond breaking

N2 is the typical benchmark for breaking a triple bond. Here, we have used a CAS(6,6) reference wave function. At the CAS-SD level, the NPE is 1.7 m$E_{\rm h}$, and the dressed CAS-SD doesn’t reduce the NPE. Here, it is necessary to use reference-dependant amplitudes to recover a low NPE : 1.0 m$E_{\rm h}$ with the $\lambda$-MR-CCSD, and 0.7 m$E_{\rm h}$ with the $\mu$-MR-CCSD.

III.2 Excited states

Triplet state of F2

We report here calculations on the triplet state ${}^{3}\Sigma_{u}^{+}$ of F2. The reference wave function was prepared in two different ways, both using restricted open-shell Hartree-Fock molecular orbitals. The first reference wave function labeled $m_{s}=1$ is a single open-shell determinant, and the second wave function is the triplet $m_{s}=0$, made of two determinants $1/\sqrt{2}(\alpha\beta-\beta\alpha)$.

To ensure that the CAS-SD is a strict eigenfunction of the $\hat{S^{2}}$ operator, we have included in $\Psi_{\rm SD}$ all the determinants with the same space part as the Singles and Doubles with respect to the CAS. These determinants are treated in the same way as Singles and Doubles and are treated variationally in the diagonalizations. Of course, those which are Triples or Quadruples with respect to $\Psi_{\rm ref}$ are excluded from the set of the $\{\alpha\}$ and have no effect in the dressing.

To reduce the computational cost, the Triples and Quadruples were not augmented with all the determinants with the same space part. The absence of some determinants gives rise to a slight deviation ($<10^{-6}$ a.u.) of $\langle\hat{S^{2}}\rangle$ from the desired eigenvalue, and it is expected to have some impact on the iterative dressing. It is worth checking the effect of this deviation from the exact spin multiplicity. The first test concerns the comparison of the $m_{s}=0$ and $m_{s}=1$ components of a triplet state.

In all the cases, the NPE of the CAS-SD (2.6 m$E_{\rm h}$) is improved to a value of 1.1–1.5 m$E_{\rm h}$. As expected the two variants of the MR-CCSD are strictly equivalent for $m_{s}=1$. Indeed, for both variants the usual single-reference amplitudes $c_{i}/c_{0}$ are recovered. The amplitudes of the $\lambda$-MR-CCSD lower the curve by 1 m$E_{\rm h}$ when going from $m_{s}=1$ to $m_{s}=0$. The dressed CAS-SD also gives a lower energy, but only by 0.5 m$E_{\rm h}$. This is due to the increased number of degrees of freedom in the fit of the amplitudes as no additional constraint is imposed to enforce the $m_{s}$-invariance. But when the reference-dependence is introduced via the $\mu_{i}$, it is imposed to recover the CAS-SD wave function which is $m_{s}$-invariant, and this step compensates the additional freedom gained in the fitting, and the $m_{s}=1$ and $m_{s}=0$ MR-CCSD curves differ by less than 0.1 m$E_{\rm h}$.

If one considers the error on the singlet-triplet gap with respect to the Full-CI reference, it appears clearly that the $\mu$-MR-CCSD gives the most accurate results, with errors lying between 0.1 m$E_{\rm h}$ and 1.3 m$E_{\rm h}$ along the curve.

Avoided crossing in HF and LiF

We have calculated the potential energy surfaces of the two lowest ${}^{1}\Sigma^{+}$ states of HF, using as reference wave function the CAS(2,2) with state-averaged CAS-SCF molecular orbitals in the aug-cc-pVDZ basis set. Figure 3 shows the NPEs of the ground and excited states. In the ground state, the NPE is 1.8 m$E_{\rm h}$ for both MR-CCSD variants, but the $\lambda$-MR-CCSD shows some numerical instabilities, as opposed to the $\mu$-MR-CCSD which gives a very smooth curve.

In the excited state, the situation is different : surprisingly the best NPE is obtained by the CAS-SD, but this may be due to the fact the molecule is particularly favorable to the CAS-SD : both H and H+ have no correlation energy, so the wave functions of the dissociated molecule can be expressed as a product. The two variants of the MR-CCSD agree at short and long distances, but they differ significantly between 2 and 3.0 Å, after the region of the avoided crossing. To understand these differences, we have plotted the two eigenvalues of the two state-specific Hamiltonians, one dressed for the ground state and one dressed for the excited state. It appears that between 2 and 3.0 Å, the lowest eigenvalue of the Hamiltonian dressed for the excited state is very badly described. The reason is that the fitting procedure for the amplitudes is a least-squares fit on the CAS-SD wave function of the state of interest, so the quality of the dressing for the determinants which have small coefficients on the state of interest but large coefficients on the other state will be very low. The $\lambda$-MR-CCSD has amplitudes which depend less on the wave function, so the quality is equivalent on both states, and the choice of these amplitudes is better suited for calculating excited states within the same symmetry.

In figure 5 we have represented the avoided crossing of LiF, also calculated with the aug-cc-pVDZ basis set. The physical situation is similar to HF, but the energy difference between the ground and the excited states is much smaller. A striking result is that the $\lambda$-MR-CCSD, although being state-specific, is able to reproduce very well the whole potential energy surfaces of both states. The position of the avoided crossing is very well reproduced by the three methods : the CAS-SD crosses at 6.3 Å, the Full-CI crosses at 6.8 Åand the dressed CAS-SD and the two MR-CCSD variants cross at 6.9 Å. The $\mu$-MR-CCSD and $\lambda$-MR-CCSD coincide in the short-range ($\leq 5~{}\AA$) and in the long range ($\geq 7.2~{}\AA$), but when the two states become very close in energy in the region of the crossing the dressed CAS-SD and the $\mu$-MR-CCSD are unable to give sensible values. This disappointing result motivates a future work on a multi-state $\mu$-MR-CCSD.

Sensitivity to the choice Molecular Orbitals

The $\mu$-MR-CCSD algorithm we propose is in the Jeziorski-Monkhorst framework, so it is not invariant with respect to the choice of molecular orbitals. In this section, we checked its sensitivity to the choice of the MO set by comparing results obtained with pseudo-canonical CAS-SCF orbitals and with localized MOs in the F2 and N2 molecules (figure 6).

In the F2 molecule, using localized MOs is a better choice than the pseudo-canonical MOs. The best NPE is obtained by the $\mu$-MR-CCSD method with a value of 1.1 m$E_{\rm h}$. In the case of N2, the situation is different : the NPE of the $\lambda$-MR-CCSD goes from 0.9 m$E_{\rm h}$ to 3.7 m$E_{\rm h}$, and the NPE of the $\mu$-MR-CCSD increases from 0.7 m$E_{\rm h}$ to 1.1 m$E_{\rm h}$. On the other hand, the dressed CAS-SD gives a better NPE with local orbitals, going from 1.7 m$E_{\rm h}$ to 1.1 m$E_{\rm h}$.

The fact that the $\mu$-MR-CCSD is less sensitive to the MO set than the $\lambda$-MR-CCSD can be understood. By changing the MO set, a single excitation rotates into a combination of single and double excitations. In the $\lambda$-MR-CCSD method, the amplitudes are calculated by taking into account the matrix elements of the Hamiltonian, which are of different nature depending on the degree of excitation, so the amplitudes are expected to change significantly. In the $\mu$-MR-CCSD variant, the amplitudes are adjusted in such a way that they fit the CAS-SD wave function, which is invariant by rotation of the MOs. Therefore, it is expected to be more robust with respect to the MO set.

IV Conclusions

We have proposed a method to determine reference-independent amplitudes by fitting the CAS-SD CI vector. These amplitudes may be used to perform a state-specific iterative dressing of the CAS-SD Hamiltonian in order to take into account the effect of the Triples and Quadruples in the spirit of the Coupled Cluster formalism. Alternatively, these amplitudes may be rescaled to reproduce the exact coefficients of the singles and doubles to introduce a reference-dependent character. In that case, the CAS-SD CI vector is recovered by the application of $(1+\hat{T})$ on the reference wave function, so we reach here the Jeziorski-Monkhorst Coupled Cluster formalism.

The CAS-SD dressed with reference-independent amplitudes gives excellent results for single-bond breaking (F2 and ethane) and the simultaneous breaking of the two O—H bonds of water, with a non-parallelism error lower than the milli-Hartree. When the active space becomes larger, it is necessary to go to the reference-dependent MR-CCSD introducing the $\mu$ factors in Eq. 28. In the case of ethylene and N2, this keeps the NPE to a value close to the milli-Hartree.

We have shown numerically that the here-proposed amplitudes are not very sensitive to the value of $m_{s}$ for open-shell systems, and to the choice of the molecular orbitals. This is clearly an improvement compared the amplitudes proposed earlierGiner et al. (2016). But we have also shown that the former amplitudes are a better choice when computing excited states of the same symmetry because the here-proposed amplitudes have a much more pronounced state-specific character which may be disadvantageous if the states are too close in energy. This problem can be cured by leaving the state-specific formalism for a multi-state formalismMalrieu (2013), and this will be the object of a future work.

Acknowledgments. This work has been made through generous computational support from CALMIP (Toulouse) under the allocation 2015-0510, and GENCI under the allocation x2015081738.