Higher Energy Derivatives in Hilbert Space Multi-Reference Coupled Cluster Theory : A Constrained Variational Approach

: In this paper, we present formulation based on constrained variational approach to compute higher energy derivatives upto third order in Hilbert Space Multi-Reference Coupled Cluster (HSMRCC) Theory. This is done through the use of a functional with Lagrange multipliers corresponding to HSMRCC method, as done by Helgaker, Jorgensen and Szalay. We derive explicit expressions upto third order energy derivatives. Using (2 n + 1) and (2 n + 2) rules, the cancellation of higher order derivatives of functional parameters that are not necessary according to these rules, is explicitly demonstated. Simpliﬁed expressions are presented. We discuss several aspects of the functional used and its potential implications.


INTRODUCTION
In last two decades, coupled-cluster (CC) method [1][2][3] has emerged as the most efficient and accurate method for calculation of electronic structure and spectra where inclusion of electron correlation is an important factor.For systems dominated by single determinant, the dynamical electron correlation is the most important part of correlation.The single-reference coupled-cluster (SRCC) method compactly sums up several infinite order perturbation terms through the solution of a non-linear equation leading to an accurate description of dynamical correlation.By its very method of formulation, the much required property of size-extensivity property is automatically included in CC method.
Another reason for the emergence of SRCC as state-of-the-art method is development of efficient analytic energy derivative techniques for the calculation of molecular electronic properties [4,5].The linear response theory is the appropriate framework for this [4,6].In this theory, molecular energy gradients, hessians, dipole-moment, polarizability and other frequency dependent properties appear as energy response quantities of appropriate order in perturbation.Stationary electronic theories i.e., theories where some or all parameters are determined from stationarity of energy functional, enjoy an advantage for calculating energy derivatives.Such theories obey the Hellmann-Feynmann theorem and its generalization (2n + 1) rule, which drastically simplify the evaluation of energy derivatives [7].There have been several attempts to formulate a stationary theory based on coupled-cluster ansatz [11][12][13][14].All such attempts have resulted in theories much more complicated theories.As a result, the standard CC method is non-stationary.
As a consequence of non-stationary nature of the standard CC method, the calculation of first order energy response requires that first order response of cluster amplitudes be calculated for every mode of perturbation.However, Bartlett and co-workers [5,8] have shown that the requirement of first derivative of cluster amplitudes for each mode of perturbation can be replaced by a single perturbation-independent quantity, known as Z-vector.There have been further developments in obtaining compact expressions for first and second derivatives including orbital response [5, ?, 10].By constructing an unconstrained Lagrangian functional, Helgaker, Jorgensen and coworkers have incorporated all such developments in a single formulation which is easily applicable to higher orders [15][16][17][18] Description of molecular excited states and potential energy surfaces can not be easily done with single-reference theories such as SRCC.These are the cases which are marked by nearly equal domination of a number of determinants referred to as quasi-degenerate.This is essentially known as non-dynamical electron correlation, which can be efficiently introduced by using a multi-determinantal description of zeroth order situation [3,19,20].There have been several approaches developed in last two decades.Dominant among them is the effective Hamiltonian method [21,22], where an effective Hamiltonian is diagonalized within a suitably chosen quasidegenerate model space to approximately reproduce a part of the spectrum associated with the exact Hamiltonian.Pursuing coupled-cluster ideas have led to two different ansatz for wave operator to introduce dynamical correlation.First is the Fock-Space (FS) multi-reference coupledcluster (MRCC) method [19].It is efficient for studying states associated with a model space of specified number of of valence electrons, while at the same time considering all lower valence model spaces with fewer number of electrons in them.This makes it suitable for studying spectroscopic quantities such as excitation energies.The second method, Hilbert-Space (HSMRCC) [23], is useful for study of several interacting states with a fixed number of electrons and is suitable for potential energy surfaces.To overcome the intruder state problems associated with complete model spaces usually used in effective Hamiltonian theories several different directions are being pursued in recent years.Use of incomplete model space [19,24] intermediate Hamiltonian [21,25] and state-specific multi-reference approaches [26,27] are major ones.
It may be mentioned that the one-valence FSMRCC is equivalent to the Equation of Motion CC (EOM-CC) method.In the context of EOM-CC fully relaxed analytic response has been developed and implemented by Gauss and Stanton [28,29].However, FSMRCC method for one valence upward, and HSMRCC in particular, have no parallel in EOM-CC or any other SR based theories.Thus, it is important to formulate the linear response of effective Hamiltonian MRCC methods [30].Analytic linear response approach for a general effective Hamiltonian based on coupled-cluster method was originally formulated by Pal [31].It has been implemented to the FSMRCC [32].This has been recently used to calculate dipole-moments of open-shell molecules [33].Szalay has outlined an approach [34] based on undetermined Lagrange multipliers method used by Hegaker, Jorgensen and coworkers in SRCC [15][16][17], to outline the calculation of first-derivatives of MRCC methods, and analyze the cost of MRCC derivatives as compared to SRCC derivatives.The present authors have introduced Z-vector method for effective Hamiltonian MRCC methods and derived expressions for the first energy derivative for HSMRCC theory [35].
The current work focuses on deriving expressions for higher-order energy derivatives (specifically, upto third-derivative) for HSMRCC theory through appropriate higher-order generalization of Z-vector method.For this we propose an undetermined Lagrange multiplier functional which is similar to the one used by Szalay [34].The equations for model space coefficients, cluster amplitudes, Lagrange multipliers and their first derivatives are obtained.In section II, SRCC and MRCC linear response, Z-vector method and its relation to the method based on contrained Lagrange multipliers are briefly summarized.Z-vecror method for HSMRCC is also discussed, outlining its relation to the functional proposed by Szalay.In section III, using a functional with undetermined Lagrange multipliers corresponding to HSMRCC method leading to Z-vector equations, expressions upto third derivatives for a specific state are derived.Equations for the Lagrange multipliers and its derivatives are presented.Section IV summarizes the results with some general comments on the results of section III.

SRCC AND MRCC LINEAR RESPONSE THEORY
In this section, key developments in SRCC and MRCC linear response theory are discussed.The relation of Z-vector to constrained variational method is examined in detail.The advantage of the latter method is highlighted for use in computing higher-order derivatives.

Single Reference Coupled Cluster Method
The SRCC method has been thoroughly discussed in literature.It is obtained by using exponential ansatz e T acting on a determinant in Schrodinger equation where T is the cluster operator with usual second quantized definition.
In the above expression i, j, . . .refer to orbitals occupied in |Φ , and a, b, . . .refer to orbitals unoccupied in |Φ and a a and a i are orbital creation and annihilation operators respectively.Although T is N -body operator, it is usually truncated.Truncation of T to one and two-body parts is quite common and referred to as CCSD, although upto four-body truncations have been reported in literature.By substituting ansatz into equation and premultiplication with e −T followed by projection with hole-particle determinants |Φ ab... ij... (henceforth denoted as |Φ q with q denoting collective hole-particle excitations respect to |Φ ), we get the CC equations.
The equations for cluster amplitudes are non-linear in nature and can be solved by iteration.

Analytic Linear Response for SRCC method.
A computational framework to analytically calculate various properties with SRCC method was first outlined by Monkhrost using linear response approach [36].In this framework, various molecular properties are defined to be various order derivatives of energy with respect to a small external perturbation strength parameter g.By expanding Energy E, Hamiltonian H and the cluster amplitudes T as a Taylor series around zero perturbation strength, and collecting terms of same order in g, expressions for various derivatives of CC Energy can be obtained.For first order, neglecting orbital response, the equations are, AT (1) = B (g) where, AT (1) = Φ q |e −T [H, T (1) ]e T |Φ , ∀q (8) In above equations, superscript (1) on T and H refers to first derivative of respective quantities with respect to external perturbation strength parameter g.Terms within {} are used to denote column vectors and superscript t denotes transpose.Similar equations can be derived for higher derivatives of T and E. Inclusion of orbital response (which is important for geometric derivatives) was first done by Jorgensen and Simons [37] who derivided detailed expressions upto second derivatives using the second-quantization framework.

Z-vector method and other developments
As outlined, SRCC, being non-stationary theory, does not have the advantage of (2n + 1) rule.Thus, T (1) is required for calculating first derivative of energy E (1) .This has to be done for every mode of perturbation which is a disadvantage.A step towards eliminating this disadvantage was taken by Bartlett and coworkers [8,5,9].They make use of a technique based on Dalgarno's interchange theorem [38] used by Handy and Schaefer for configuration interaction (CI) energy derivatives [39].The technique, known as Z-vector technique, derives algebraic expressions for reducing the number of linear equations for orbital rotations and cluster amplitudes.By inverting equation (7) and substituting in ( 6), E (1) can be reorganized as, Equation ( 13) is a perturbation independent equation providing the Z t vector.The advantage of such reorganization is, unlike earlier equation (7), equation ( 13) needs to solved only once.This Z-vector method is in some sense equivalent to (2n + 1) type of rule for non-stationary methods.Further simplifications have been carried out by introducing effective CC density matrices [5,9] much akin to CI derivative developments.The technique of Rice and coworkers [40] has also been applied to reduce the number of AO to MO transformations.First applications of CC analytic derivatives have been reported by Bartlett and coworkers [41] and Scheiner and coworkers [42].

Method of Undetermined Lagrange Multipliers
The conceptually simple procedure of eliminating T (1) in favor of a Z-vector is somewhat cumbersome for higher orders and has been pursued by Salter and coworkers [10].However, Helgakar, Jorgensen and coworkers [17], have pursued an attractive alternative formulation of CC derivative which automatically incorporates Z-vector technique to all orders.Such a formulation proceeds by constructing a Lagrange functional with undetermined multipliers Λ = {λ q , ∀q} corresponding to CC equations (5) as follows.
The functional may be viewed as approximation to the full Extended Coupled Cluster (ECC) functional of Arponen [13].The full ECC functional been used by Pal and coworkers to derive the first and higher order energy derivatives using a stationary approach [43][44][45].Optimization of functional ( 14) leads to same equations as cluster amplitude equations (5) and Z-vector equations ( 13) with Λ taking the role of Z-vector.In this formulation, it is transparant to derive expressions for higher order derivatives.While T (1) obey (2n + 1) rule, the undetermined multipliers λ obey (2n + 2) rule [16].

Multi-Reference Coupled Cluster (MRCC) Theories.
Extension of coupled-cluster method to quasi-degenerate situations demanding multi-reference description has been non-trivial.Different approaches were followed by different researchers.
Among them, the approach of constructing an effective Hamiltonian in a complete model space has emerged as the standard one [3,21,22].We briefly summarize the approach in the following.
A strongly interacting M -dimensional complete model space is considered by distributing a given number of electrons among a conveniently chosen set of valence orbitals.This model space, P 0 , is assumed to approximate the quasi-degenerate target space P , spanned by M exact states of exact Hamiltonian H. Hence, it contains the zeroth-order approximations for all the exact states under consideration.Dynamical correlation is brought in via wave-operator Ω mapping model and target space as, The wave-operator, Ω, through an appropriate parametrization, is constructed by solving generalized Bloch equation.
Bloch equation is simpliy a restatement of the schrodinger equation for all targeted states.The exact energy and zeroth order approximations of target states within model space, |Ψ (0) i , are obtained by diagonalizing the effective Hamiltonian.
Two kinds of parametrizations for the wave-operator Ω have been widely discussed in literature.The first one involves the use of a common reference vaccum |Φ for defining a single secondquantized cluster operator T connecting model spaces with different number of valence electrons from zero to a specified number used to build the model space.This leads to Fock-Space MRCC (FSMRCC) [19].T is so constructed that the corresponding wave-operator, Ω = e T , is able to yield exact states for different model spaces connected by T .FSMRCC has been successful in accurate computations of spectroscopic quantities such as ionization potential, electron affinity, excitation energies [46,47].
The second formulation does not use a single reference vaccum, but uses as many vacua as the number of states in the manifold, with one cluster operator for each vaccum.This is referred to as Hilbert-Space MRCC (HSMRCC) [23].Unlike FSMRCC, HSMRCC does not require the use of model spaces with different number of valence electrons.The wave-operator is written as, where T µ is the cluster operator with respect to vaccum |Φ µ similar in structure to SRCC cluster operator.T µ does not contain excitations leading to states within model space.Using the above ansatz in equations ( 16) and ( 17) and projecting with determinants, |Φ q (µ) , hole-particle excited with respect to |Φ µ (q denotes collective orbital indices involved in the excitation), we get the following HSMRCC equations. µ Ciν and C µi refer to left and right eigenvectors of the effective Hamiltonian H ef f .HSMRCC has been pursued by various researchers and has been found to be useful for description of potential energy surfaces [48][49][50].Spin adaptated formulations have also been pursued and successfully applied to various chemical systems [51].
There have been several other MRCC approaches such as Intermediate Hamiltonian approach [21,25], use of incomplete model spaces [19,24] and state-specific MRCC [26,27].However, the current work focuses on complete model space effective Hamiltonian based MRCC approaches with special emphasis on HSMRCC theory.

Analytic Linear Response for MRCC Theories
Application of linear response to obtain expressions for MRCC derivatives was initiated by Pal, who, following Monkhrost's approach outlined the formulation for the effective Hamiltonian based MRCC theory [31].Specific expressions were obtained for one hole, one particle and hole-particle sectors of FSMRCC theory [32].The structure of the original nonvariational formulation did not include prescription of Z-vector and thus were not applicable for higher energy derivatives and even gradients.Applications of the formalism to compute first order molecular properties have however been carried out in recent years [33].The formalism has also been extended to time-dependent perturbation case [52].
On the other hand, Szalay was the first to extend Helgaker and Jorgensen's method by proposing a Lagrangian functional for a chosen state in the manifold which yields MRCC cluster amplitude and effective Hamiltonian diagonalization equations [34].His work mainly concerned with esti-mating the relative cost of MRCC first order response calculations as compared to SRCC response equations.
First efforts towards directly introducing Handy-Schefer Z-vector type technique for the effective Hamiltonian, H ef f , were carried out by Pal and coworkers [53].Working on FSMRCC response equations for one hole model space, they concluded that only highest sector amplitude derivatives can be eliminated from H ef f (1) under degenerate diagonal H ef f .Further work on HSM-RCC theory has shown that it is not possible to eliminate cluster amplitude derivatives in effective Hamiltonian first derivative with a single Z-vector for all states, although M 2 number of independent Z-vectors would be sufficient.Recently, however, the present authors have showed that it is possible to eliminate the cluster amplitude derivatives for a chosen state energy derivative.
They have derived detailed equations for HSMRCC Z-vector for a chosen state and developed efficient expressions for HSMRCC energy first derivative by extending the idea of effective CC density matrices [35].

HSMRCC HIGHER ENERGY DERIVATIVES
Extension of algebraic Z-vector method is difficult for higher order energy derivatives, although it has been pursued in SRCC context by Salter and coworkers [10] to obtain second derivative expressions.It is advantageous to go over to constrained variation approach [17].The approach generalizes the Z-vector method to all orders retaining the simplicity of HSMRCC method.In this section, we derive expressions for HSMRCC energy derivatives upto third-order for a chosen state.For this, we propose a functional with Lagrange multipliers corresponding using HSMRCC equations, ( 19)-( 23), as constraints.The functional leads to equations for Lagrange multipliers corresponding to cluster amplitude equation (19) to be the same as Z-vector equations derived recently [35].Using (2n + 1) rule for cluster amplitudes and model space coefficients, (2n + 2) rule for Lagrange multipliers, the expressions for second and third order are written in a simple form.

The Lagrangian functional and obtaining linear response
Before giving explicit expression for the Lagrangian functional which corresponds to HSMRCC method, we assume the following abbreviations.We denote the set of cluster amplitudes for all vacua by T , the set of all model space coefficients for the i-th state by C i and the set of all Lagrange multpliers corresponding to HSMRCC cluster amplitude equations (19) by Λ.We further collectively refer the set of quantities T ,C i ,Λ and the energy of i-th state E i (the Lagrange multiplier corresponding to biorthogonal conditions (23)), by Θ.The following equations summarize the above abreviations.
Now, the Lagrangian functional which corresponds to HSMRCC equations ( 19)-( 23) is constructed as, where HSMRCC equation for T , (19), is denoted by E q (µ) = 0. Making J stationary with respect to Θ generates sufficient number of equations to determine these parameters.The above functional is similar to one used by Szalay [34] in context with HSMRCC first derivatives.However, the present functional does not introduce C i dependence in second term.As a consequence, J leads to Z-vector equations derived by eliminating the first order response of cluster amplitudes, T (1) from energy first derivative expression.Although current formulation results in slightly complicated expressions for effective CC density matrices, it will be advantageous as discussed in next section.
To derive expressions for various order derivatives, linear response theory is used.A small perturbation H (1) with strength parameter g is introduced into the Hamiltonian H as, H(g) = H + gH (1)  (29) Since Θ is determined at all strenghts of perturbation g, the functional J becomes a function of g, denoted by J (g, Θ).Response equations can be obtained in two different approaches.First approach, followed by Helgaker, Jorgensen and coworkers [17] and Bartlett and coworkers [10] is to expand the functional J and the stationary equations obtained for Θ as a Taylor series in strength parameter g.Terms of same order in g in Θ equations are collected and equated to obtain hierarchical equations for various response quantities Θ (n) .Second approach, proposed by Pal [12,54] in the context with stationary CC theories, is to derive response equations of any required order as a stationary equations.In this approach, J (g, Θ) is expanded as, J (g, Θ) = J (0) + gJ (1) It should be noted that J (n) is a functional of quantities Θ (m) m = 0, n .All response equations upto a required order n can be derived by making the functionals J (k) k = 0, n stationary with respect to Θ (m) m = 0, n (m ≤ k) .This leads to the following equations.
This includes the n = 0 case corresponding to the unperturbed (zeroth order) equations as well.Henceforth, we drop the superscipt for zeroth-order quantities.It has been shown that there is a large amount of redundancy in the above equations.Hence it is sufficient to solve This leads to hierarchical equations for Θ (m) ∀m = 0, n .In this work, we follow the second approach in deriving response equations.According to (2n + 1) rule associated with stationary methods, we need Θ and Θ (1) to obtain expressions upto thrid order energy derivatives, and hence equations (32) need to be solved upto n = 1.It should be noted that both approaches are entirely equivalent and lead to identical equations with a given functional.

Response response equations upto first order
Detailed expressions for J (n) upto n = 3 are given in Appendix A, equations (A-4)-(A-6).To derive response equations for Θ upto first order, we use expression for J along with equation (A-1).Zeroth-Order quantities, Θ, can be obtained by When applied for different parameters Θ, the above equation gives HSMRCC equations for cluster amplitudes (19) and the eigenvalue equations ( 20)-( 21).In addition it gives equations for Λ as, where subscript τ i (µ) indicates differentiation with respect to specific cluster amplitude t i (µ) and τ i (µ) is the hole-particle excitation operator associated with t i (µ).The quantity [E q (η)] τ i (µ) has been referred to CC Jacobian by Jorgensen and coworkers [17] and in this case it is the HSMRCC Jacobian.Diagrammatic representation for the above equation can be easily obtained as outlined in [35].The above equation is the same as the Z-vector equations derived recently by the present authors [35] using elimination technique of Handy and Schaefer [39].
The first order quantities, Θ (1) are necessary for calculating higher order energy derivatives and can be obtained by, The above equations lead to following equations for first order quantities Θ (1) .It should be noted that these equations depend on zeroth-order quantities Θ.
All terms used here are defined in Appendix A. These equations should be solved in the same order.Equation ( 40) is the equation for Λ (1) .It not only depends on first derivatives of T and C i , but also on zeroth order quantity Λ.The equations ( 36) and ( 40) reveal the same structure pointed out by Jorgensen and coworkers [17] i.e., T (1) and Λ (1) are related by the same HSMRCC Jacobian.The only difference between both equations is in the inhomogenous part.

Simplified expressions for Energy derivatives
Energy derivative of n-th order, E (n) i , is just the value of functional J (n) , denoted by J (n) opt , when the stationary values of Θ (m) m = 0, n are substituted in it.Hence, J (n) opt can be considered as the required energy derivative and treat E i as another Lagrange multiplier.It has been shown that C (n) i and T (n) obey (2n + 1) rule and Lagrange multipliers Λ (n) and E (n) i obey (2n + 2) rule.However, the expressions obtained by simple substitution as above do not take advantage of the above rules.Hence the expressions must be simplified by explicit application of these rules.This eliminates any unnecessary higher order derivatives present in these expressions.Elimination is carried out by referring to appropriate response equations, including the zeroth-order response equations.
For J (n) opt n = 1, 3 the quantities which need to be eliminated are given below, by the application of (2n + 1) and (2n + 2) rules.
i , T (3) , T (2) , E i , Λ (3) , Λ (2)   In Appendix B, expressions for J (n) opt n = 1, 3 obtained from (A-4)-(A-6) are rearranged to indicate explicit cancellation of higher derivatives of C i , Λ and E i .Terms indicated in {} can be eliminated by application of response equations of appropriate order.This is discussed in following subsections.Elimination of higher derivatives of T from the remaining expressions will be demonstrated seperately.
Elimination of T (1) which is present in the surviving terms of (B-1), requires application of zeroth order Λ equations.For demonstate this, H νµ ef f (1) and E (1) q (η) are expanded to seperate out terms containing T (1) .All such terms cancel precicely from Λ equations, (34).Hence, simplified expression for J (1) opt is given by, λ q (η) E (1)  q (η) where subscript T (n) indicates retention of terms containing T (m) m = 0, n .This expression has been further simplified by making use of state-dependent effective CC density matrices [35].

Simplified expression for J (3) opt
Expression for Energy third derivative J (3) opt is given in (B-3) along with groupings necessary for elimination.Equation (40) for Λ (1) is the only additional equation to be solved in addition to T (1) and C (1) i equations.C (3) i and C (2) i can be eliminated through C i and C (1) i equations respectively.Similary, Λ (3) , Λ (2) , E (3) i , E (2) i can be trivially eliminated by applying appripriate response equations.

DISCUSSION
Results of the previous section demonstrate the advantage of constrained variational approach for non-stationary theories over elimination technique.In this formulation, it is particularly easy to derive expressions for higher derivatives and apply (2n + 1) and (2n + 2) rules to simplify them.
The nature of cancellation of terms in J (n) opt containing higher order derivatives of Θ, is also clear.Derivatives of Lagrange numtipliers Λ and E i , which are not required by the (2n + 2) rule, cancel from lower order response equations of cluster amplitudes T and biorthogonality equations.Cancellation of higher order derivatives of C i that are not required by the (2n + 1) rule obeyed by it, occurs again by lower order response equation of its conjugate model space coefficient.Since the Lagrange multiplier E i appears in these equations, some derivatives of E i go into these cancellations.It should be noticed that while Lagrange multipliers Λ and Λ (1) appear in J (3) opt , the same does not happen for E i and E (1) i .In J (3) , only E (1) i appears, while E i goes into cancellations involving higher order derivatives of C i .Cancellation of higher derivatives of T not required by (2n+1) rule, happen through lower order Λ response equations.Finally, for even order derivatives, as exemplified in J (2) opt , further simplifications are possible through the use of response equations for C i , resulting in eliminating the terms where both C (1) i and T (1) are present.In the context of SRCC linear response theory, it has been known that the constrained variational approach is identical to Z-vector method introduced much earlier, with the Lagrange multipliers Λ becoming identical to the Z-vector [45].However, this is not the case in MRCC, where several possibilities of functional open up.The functional proposed in this work leads to the same equations as the elimination method followed in Z-vector approach.There are several other possible functionals worth pursuing.One of them has been proposed by Szalay [34].In his functional, a factor of Ciµ C µi is introduced in the second term of (28).This leads to expressions for effective CC density matrices which are closer in structure to SRCC counterparts [9].On the other hand, such a functional leads to Λ equations whose homogenous part depends on C i .As a result, T (n) and Λ (n) do not share the same HSMRCC Jacobian, as it happens in SRCC.On the other hand, the functional in this work leads to same HSMRCC Jacobian for T (n) and Λ (n) .Although this is a clear advantage, other possibile functionals should be explored.However, it should be noted that all such functional are equivalent in the sense that they provide same energies and derivatives differing only in the form of expressions.
Another important aspect of constrained variational approach is the intrinsic state-dependency of the functional.It should be noted that the state-dependency of Z-vector is not special to MRCC methods, but it has been observed in Equation of Motion CC method [28,29].The present authors have investigated whether T (1) can be eliminated from H ef f (1) in favor of a single Z-vector.It has been concluded to be not possible for a general case, because of matrix nature of H ef f (1) .Hence, it is necessary to become state-selective and eliminate T (1) from the energy deriative of a specific state.
It has been recognized that constrained variation approach for SRCC is related to bivariational CC methods investigated by Arponen [13].The Single-Reference Extended CC (SR-ECC) formulated by Arponen is attractive in many respects and has been investigated for its use in calculating molecular properties [43][44][45].Apart from being stationary, terminating and size-extensive nature of functional, with all their advantages, it is also known to sum much larger class of perturbation theory diagrams as compared to SRCC.Hence it is expected to be highly accurate and appropriate for use to calculate nonlinear molecular electronic properties.Generalization of SR-ECC to multi-reference situations (MR-ECC) has not been reported in literature.The nature of constrained variation approach for HSMRCC clearly indicates that MR-ECC has to be state-selective theory, much akin to the decontracted state-selective MRCC approach proposed by Mukherjee and coworkers [27].The non-uniqueness of the constrained variational functional indicates that many different formulations of MR-ECC may be possible.This is similar to the possibilty of having varieties of state-selective non-stationary MRCC theories.This line of research is worth pursuing.