- Email: [email protected]
Full effect of triples in a valence universal multi-reference coupled cluster calculation
30 July 2002
Chemical Physics Letters 361 (2002) 298–306 www.elsevier.com/locate/cplett
Full effect of triples in a valence universal multi-reference coupled cluster calculation Sudip Chattopadhyay a, Ashis Mitra b, Debasis Jana a, Pradipta Ghosh a, Dhiman Sinha c,d,* a
Department of Physical Chemistry, Indian Association for the Cultivation of Science, Calcutta 700 032, India b Department of Chemistry, R.K. Mahavidyalaya, Kailasahar, Tripura 799 277, India c Department of Pure Chemistry, University College of Science, 92 A P C Road, Calcutta 700 009, India d Department of Chemistry, Tripura University, Suryamaninagar, Tripura 799 130, India Received 2 April 2002; in final form 10 June 2002 This article is dedicated to Professor Mihir Chowdhury on the happy occassion of his 65th birthday.
Abstract Full three-electron coupled cluster (CC) operators (triples) along with all singles and doubles are included in a valence universal multi-reference CC calculation to estimate the change in the correlation effect for ionization potential (IP) calculations. We study the vertical IPs of two typical systems, N2 and CO, which go in the opposite direction, i.e., the results were worsened on inclusion of triples when working with relatively smaller basis. However, upon reaching the basis set saturation the expected trend in the results is observed and they are close within 0.1 eV to experimental values almost in every IP. Ó 2002 Published by Elsevier Science B.V.
1. Introduction The single-reference (SR) coupled-cluster (CC) method in varying degrees of truncation is now well established as a very successful tool for treating electron correlation of closed-shell systems to a high degree of accuracy [1,2]. This is size-extensive at any level of truncation of the rank of the cluster operator. However, the SRCC theory, truncated at the singles and doubles (SD) excited cluster level is not adequate
*
Corresponding author. E-mail address: dhimansinha@rediffmail.com (D. Sinha).
for the states with quasi-degeneracy emerging from multi-reference (MR) character. There are several ways to overcome this difficulty. One of them being an explicit accounting of three- and/ or four-body cluster amplitudes [3,4] within the SRCC framework. Another alternative scheme estimates the higher-body cluster amplitudes by augmenting the SR-CCSD equations by inclusion of these cluster amplitudes from a complete active space function [5]. The numerical implementation along this line is computationally demanding. Rather, conceptually the most appropriate one would be to treat quasi-degeneracy stemming from non-dynamical correlation effects by starting from an active space of quasi-
0009-2614/02/$ - see front matter Ó 2002 Published by Elsevier Science B.V. PII: S 0 0 0 9 - 2 6 1 4 ( 0 2 ) 0 0 9 7 4 - 0
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
degenerate determinants, while the remaining dynamical correlations can be included via a cluster expansion of the wavefunction involving virtual determinants. This approach is thus a pure MR generalization of SRCC method. The early multi-reference CC (MRCC) formulations were all based on effective hamiltonian spanned by the reference determinants, constituting a complete model space. There are two broad categories of such MRCC methods, viz. the valence universal (VU), or the Fock space [6–8] and the state universal (SU), or the Hilbert space [9,10] ones. While the SU-MRCC method is more suited for the studies of potential energy surfaces, the VU-MRCC approach is convenient to provide energy differences of spectroscopic interest (IPs, EAs and EEs) directly relative to some common reference state. There exists another suite of methods of direct calculation of energy differences starting from the CC or CC like representation of the ground state. These are SAC-CI [11,12], CC-based linear response theory (LRT) [13], and the CC-based equation-of-motion (EOM) [14] methods. While all these three methods are equivalent, they are generally different from the corresponding VU-MRCC methods. However, it has been shown that, for the principal IPs and EAs the VU-MRCC [7,8] is equivalent to SAC-CI [11,12], LRT [13] and EOM [14] IP and EA CC methods. In this Letter we focus on the IP problem via the VUMRCC method. The VU-MRCC, despite its formal rigor, is seriously beset by the problem of intruder states [15]. However, implementation of incomplete model spaces [16,17] might seem to be a possible way out, but it is not so in practice. This is because, during the generation of energy differences of spectroscopic interest, it is appropriate to provide at least a qualitatively correct description of the high-lying states perturbed by the intruder states as they form a part of the same manifold of excited states, rather than ignoring them completely from the description. The concept of intermediate hamiltonian [18] offers a real possibility of avoiding the determinantal effects of intruders while preserving most of the advantages of the effective hamiltonian framework.
299
Mukherjee and co-workers [13,19] had earlier proposed a technique to improve the convergence behaviour of the VU-MRCC equations [7,8] in the presence of intruders using the idea of eigenvalueindependent partitioning (henceforth to be referred as EIP-MRCC). The EIP-MRCC does not differ from the traditional VU-MRCC method as far as the basic theoretical constructs, viz. the definition of the valence rank, the hierarchical decoupling of the various valence rank problems, etc., are concerned. It nonetheless provides an alternative avenue for solving the VU-MRCC equations having certain computationally strategic advantages. The EIP technique converts the nonlinear VU-MRCC equations for any model space into a set of non-hermitian eigenvalue equations. This opens up the interesting possibility of arriving at the desired solutions through the use of stable root-search and root-homing procedures used in eigenproblems [20,21] to determine the desired low-lying roots of spectroscopic interest one by one iteratively, preserving the norm of the corresponding eigenvector. Thus the method is less prone to divergence than the corresponding nonlinear equations of traditional VU-MRCC. It is noteworthy that the intermediate hamiltonian method via EIP in the context of CC framework have also been developed by Mukherjee and coworkers [22]. The recently formulated state-specific (SS) method [23] is the most effective strategy to tackle the intruder state problem, while retaining most of the important traits of a good many-body model. Pertinent to mention is that the SS methods based on the Fock space MRCC approach still remains an uncharted area of electronic structure theory. As noted, unlike the IP calculation employing the SRCC methods, which require separate calculations for each state, the VU-MRCC approach provides results for many states at once which are numerically at par with other methods. It is worthwhile to point out one of the major distinguishing features of the VU-MRCC approach from the SR-based strategies. This is in regard to the excitation operators in the two formulations – certain double excitations in the SR method only appear as triple excitations in the VU-MRCC. Thus it would be interesting to explore the possi-
300
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
bility of developing VU-MRCC schemes including higher-body cluster operators (triples and quadruples) and compare their performance with other simple VU-MRCC methods and the SR calculations. Haque and Kaldor [24], Mukherjee and co-workers [25] and Bartlett and co-workers [26] have included the lowest order effects of triples in their VU-MRCC formulations with some demonstrative results. Also, some methods like the complete fourth order Møller–Plesset perturbation theory (MP4) calculation [27], and the SAC-CI (general-R) method [28], CC based-LRT [29] and EOM [30] have taken into account the triples/ higher excitations extensively. In the present Letter, we report the effects of the inclusion of triples to our EIP-MRCC approach and examine the performance of this by some pilot IP calculations. In this application for the excited states, we have considered the full set of determinants generated by the doubles and the triples, below a certain threshold energy. However, in the ground state, which is a closed-shell, the contribution of triples are considered through the term VN T2 =D (see [3] for a detailed discussion) only, where VN and T2 are, respectively, the two-body part of the hamiltonian and two-body closed-shell CC operators, both in normal order, and D is the accompanying orbital energy difference. An important feature of the inclusion of triples is that the dimension of the determinantal set increases tremendously, by several tens of thousands or even much more, as against the set of only a few hundred determinants appearing while using no triples at all. We have observed that a careless truncation in the set of determinants may tamper the results significantly. In this Letter we have calculated all the vertical IPs of the outer valence region of N2 and CO using four different basis sets each, augmented gradually towards the saturation limit. Our aim is to see, to what extent the full onset of triples in this method enhances the CC result towards experiment. Interestingly, relatively smaller basis sets gave apparently poorer IPs for both the systems while including triples in the calculation. However, upon achieving a basis set saturation, we could obtain consistently better outer valence IPs for the two systems. A similar trend has also been observed by
Raghavachari [31] in regard to the effect of basis set and electron correlation effects on electron affinities ðEAÞ of some first row atoms. They found that with relatively smaller basis, second order Møller–Plesset perturbation (MP2) results for EA were sometimes, e.g., for O and F, much better than the complete MP4 one, even though the latter included higher excitations extensively. But the MP4 results turned out to be consistent for all the atoms they had taken when using the bigger basis sets. We conclude this section by putting forth the organization of this Letter. In the section to follow, we present a succinct discussion on the EIP-MRCC theory. In Section 3, we discuss the algorithmic consideration that leads to the numerical implementation of our method. In this section itself we provide the vertical IPs of two interesting systems (N2 and CO) as our illustrative application to judge the efficacy of our EIPMRCC method. We conclude the Letter by summarizing remarks in Section 4.
2. Theory We briefly describe here the mathematical framework of our EIP-MRCC method. One starts from the VU Bloch equation [7,8] H XP ðk;lÞ ¼ XHeff P ðk;lÞ ; 0 6 k 6 m; 0 6 l 6 n;
ð1Þ
where H is the atomic or molecular hamiltonian in normal order and P ðk;lÞ is the projector for (k holes and l particles), or simply (k; l), model space [7,8]. We shall explain the effective hamiltonian Heff shortly. The wave operator, X is expressed as a VU normal ordered exponential ansatz [7,8]: X ¼ exp S ½m;n ( !) m X n X ðk;lÞ ¼ exp S ; ð2Þ k¼0
l¼0
S ðk;lÞ being the (k; l) valence component [7,8,19] of the cluster operator S ½m;n , where [m; n] means the collection of all (k; l) valence components upto the highest (m; n) valence component. S ðk;lÞ acts upon P ðk;lÞ space determinants to generate the virtual Q
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
space and R space determinants (we sub-divide the entire virtual space into two sub-spaces, viz. Q and R for reasons cited below). The VU-MRCC method calculates ionization potentials as direct difference energies instead of calculating the state energies Egr ðN Þ and EðN 1Þ individually to find out their difference, where N stands for the total number of electrons in the closed-shell system. In the approximate level calculation, that is, when we apply truncation scheme, the methods that calculate difference energies directly do a balanced cancellation of the common correlation terms as compared to the methods that calculate state energies separately. In the VU-MRCC framework, the IP calculation is called a [1,0] valence problem, and thus it includes the CC operators for (0,0) and (1,0) valence components. Consequently, X ¼ exp S ð0;0Þ þ S ð1;0Þ ð3Þ ¼ Xð0;0Þ exp S ð1;0Þ ; where Xð0;0Þ ¼ exp S ð0;0Þ , is the closed-shell wave operator, S ð0;0Þ stands for the cluster operators of the closed-shell CC calculation, and S ð1;0Þ is the cluster operator for the single ionization with simultaneous excitations of that system. The standard choices of S ð0;0Þ and S ð1;0Þ are ð0;0Þ
S ð0;0Þ ¼ S1
ð0;0Þ
þ S2
ð1;0Þ
and S ð1;0Þ ¼ S2
;
ð4Þ
where the maximum electron rank at both the ð0;0Þ closed-shell and open-shell situations is two, S1 ð0;0Þ and S2 stand for singles and doubles in the closed-shell situation. If all the hole orbitals are kept active, which has been the case in this applið1;0Þ cation, no singles S1 (Fig. 1a) should figure out in the calculation. Then, at the present level of ð1;0Þ approximation, only the doubles S2 (Fig. 1b) act upon the P -subspace that contains 1h determinants only to generate the 2h 1p determinants of the Q-subspace.
Following the procedures given in [19], one gets, on using Eq. (4), the EIP-MRCCSD matrix equation (Method I) (S and D stand for singles and doubles, respectively):
XPP YPP YPQ YQP
¼
YQQ XPP XQP
ð1;0Þ
ð1;0Þ
Fig. 1. (a) Singles, S1 , and (b) doubles, S2 erators for open-shell.
of the CC op-
MM
ðND þMÞM
ð1;0Þ
ð1;0Þ
Heff XPP ¼ XPP EPP ; where
ð5Þ
ð1;0Þ Heff
ð6Þ ð0;0Þ
¼ P ð1;0Þ H X
P ð1;0Þ , and
Rð2Þ XPP ¼ XQP or 1 Rð2Þ ¼ XQP XPP ;
ð7Þ ð1;0Þ S2
ð2Þ
where R represents the matrix of the amplið1;0Þ tudes. EPP are actually the IPs as difference energies. M and ND are the dimensions of the P ð1;0Þ subspace and the Q-subspace, respectively. To incorporate the effect of triples in our EIPMRCC theory, we choose the following scheme: (a) ð0;0Þ the three-body zero-valence operator, S3 is only approximately computed, as in the CCSDT1-a truncation [3], this corresponds to a non-iterative ð0;0Þ inclusion of S3 ; (b) the effect of the one-valence ð1;0Þ three-body operator, S3 , is, however, included ð0;0Þ exactly. Thus, S and S ð1;0Þ in our scheme are ð0;0Þ
þ S2
ð1;0Þ
þ S3
S ð0;0Þ ¼ S1 ð1;0Þ S3
(b)
XQP ðND þMÞM h i ð1;0Þ EPP :
ðND þMÞðND þMÞ
It is interesting to note that Eq. (5) is identical in form to CC-based LRT [13,32] equations for IP calculations using 1h and 2h 1p determinants. The terms in Eq. (5) are explained as follows. YPP , YPQ , YQP and YQQ are derived from the conð0;0Þ tracted composite ðH X Þ obtained from the left hand side of Eq. (1), after doing some manipulations guided by Wick’s algebra on it. The matrix X relates to the following equations [19],
¼ S2
(a)
301
ð0;0Þ ð1;0Þ
ð0;0Þ
þ S3 :
and S ð1;0Þ ð8Þ
shown in Fig. 2 corresponds to double shakeup excitations (Fig. 2a) and also correlation relaxation (Fig. 2b). In (Fig. 2b), the lines with double arrow indicate spectator scattering. We calculate the full effect of these, which is the key
302
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
(a)
(b) ð1;0Þ
Fig. 2. The various types of triples, S3 for open-shell.
of the CC operators
issue of the present work. In the figures the lines having single arrow pointing from left to right indicate the inactive hole lines, those from right to left are the inactive particle lines. The lines having double arrow from left to right are active hole lines. By similar induction as that of EIP-MRCCSD scheme, using Eq. (8) we get the corresponding EIP-MRCCSDT matrix equation as (Method II): 2 3 2 3 YPP YPQ YPR XPP 4 YQP YQQ YQR 5 4 XQP 5 YRP YRQ YRR ðNDT þMÞðNDT þMÞ XRP ðNDT þMÞM 2 3 XPP h i ð1;0Þ ¼ 4 XQP 5 EPP ; ð9Þ MM XRP ðNDT þMÞM
Fig. 3. Representative diagrams for Y -matrix elements involving P - and Q-subspaces. The circle stands for the one-body part and the square stands for the two-body and the three-body ð0;0Þ parts of H X .
where T represents that triples were taken in full for the excited state, but for the ground state they were considered at the CCSDT1-a [3] level of approximation only. Like Eq. (7), we get a similar relation 1 Rð3Þ ¼ XRP XPP ;
ð10Þ ð1;0Þ
where Rð3Þ represents the matrix of S3 amplitudes. The R-subspace contains the (3h 2p) determinants. NDT is the sum of the dimensions of the Q- and R-subspaces together. The first few roots of the Y -matrices in Eqs. (5) and (9) are determined by Rettrup iterative algorithm [20]. The difference in the so obtained two similar roots show the change due to the inclusion of triples. An interesting point about Eqs. (5) and (9) is that we start ð1;0Þ ð1;0Þ with X which explicitly contains S2 and S3 , but we need not deal with them finally, although the information of all the orbitals on the left of S ð1;0Þ vertex is certainly necessary for locating the Y -matrix elements in the computer programme.
Fig. 4. Y -matrix elements involving R-subspace. The square in ð0;0Þ YRP stands for the three-body part of HX , and that in ð0;0Þ YRQ and YRR stand for the two- and three-body parts of H X .
This is described by the disconnected lines in Figs. 3 and 4. Of course, if required, one can at once determine S ð1;0Þ amplitudes through Eqs. (7) and (10).
3. Results and discussions In this Letter we present two separate sets of calculations done for IPs with each basis. In one, comprising the Method-I, the full CCSD calcula-
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
tion was followed by the calculation of the transformed Hamiltonian, H ¼ fH Xð0;0Þ g, by using the corresponding CC amplitudes, and finally all outer valence IPs were determined by direct diagonalization or by iterative diagonalization using the Rettrup algorithm [20], depending upon the size of the IP matrix. In the other, comprising the Method-II, first a CCSDT1-a [3] calculation was done to be followed by the calculation of H using the new CC amplitudes, and finally again the outer valence IPs were determined but now only by using the Rettrup algorithm as the size of the matrix is very large. In either of the methods, the maximum electron excitation rank of the CC correlation operators both at the closed-shell and the open-shell levels were kept the same, which was for estimating the electron correlations at the two levels on the same footing. For example, in Method-I, the maximum electron rank at both the closed- and open-shells was two, whereas in Method-II, the rank was three. The CCSD and the CCSDT1-a correlation energies are given in Table 1. While constructing H in Method-II, the closedshell CC amplitudes for singles and doubles obtained after CCSDT1-a calculation were used. That is, no closed-shell CC amplitude related to triples was employed. However, the three-body H ð0;0Þ operators through VN S2 ð¼ VN T2 Þ were included. The convergence with the use of Rettrup algorithm was pretty fast and en route to convergence it took only 10–12 iterations provided that the starting iterate was properly chosen. Furthermore, by making the Rettrup FORTRAN code independent of input-output (I/O) operations, the wall-clock time is drastically reduced. Instead of reading the
303
stored Y -matrix elements from the hard disk at each iteration, they are calculated afresh, one row at a time, in core memory and simultaneously utilized. Although such a strategy increases the total CPU time, this is, however, very much less than the wall-clock time required while following the I/O dependent Rettrup code. All the virtual IPs of the outer valence region of N2 are compared with the experiment in Table 2, and the corresponding results for CO are presented in Table 3. We have performed the calculations employing four different basis sets (I, II, III and IV) for each molecule which are described below. The first basis set does not have any polarization function, whereas the latter ones include d, d and f , and d, f and g polarization functions, respectively. The d, f and g functions in the bases comprise of 6, 10 and 15 Cartesian gaussians, respectively. Basis (on each atom) I II III
IV
TZ (10s,6p)/[5s,3p] cc-pVDZ (9s,4p,1d)/ [3s,2p,1d] cc-pVTZ (10s,5p,2d,1f)/ [4s,3p,2d,1f] cc-pVQZ (12s,6p,3d,2f,1g)/ [5s,4p,3d,2f,1g]
All the basis sets were collected from: http:// www.emsl.pnl.gov:2080/form/basisform.html In Tables 2 and 3, the variation of the Hartree– Fock (HF) energies and the Koopmans IPs with
Table 1 Hartree–Fock (HF) energies and correlation energies (in a.u.) System
Basis
EHF
Corr. energy (CCSD)
Corr. energy (CCSDT-1a)
N2
I II III IV
)108.900582 )108.954759 )108.984150 )108.991433
)0.251158 )0.319201 )0.405526 )0.430544
)0.262531 )0.332081 )0.425578 )0.452686
CO
I II III IV
)112.707913 )112.749783 )112.781003 )112.789142
)0.231969 )0.303745 )0.388982 )0.432384
)0.245393 )0.316097 )0.407992 )0.453468
re for N2 2.074 a.u., for CO 2.13 a.u.
304
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
Table 2 Vertical ionization potential (in eV) of N2 States
Basis
Koopmans
EIP-MRCCSD
EIP-MRCCSDT
3r1 g
I II III IV
17.15 17.06 17.23 17.25
14.90 15.22 15.62 15.74
14.77 15.07 15.45 15.56
I II III IV
16.97 16.56 16.68 16.70
17.33 16.94 17.22 17.33
16.95 16.69 16.96 17.07
I II III IV
21.09 21.08 21.18 21.19
18.00 18.51 18.87 18.97
17.75 18.28 18.60 18.71
1p1 u
3r1 g
a
Experimenta 15.6
17.0
18.8
Ref. [33].
Table 3 Vertical ionization potential (in eV) of CO States
Basis
Koopmans
EIP-MRCCSD
EIP-MRCCSDT
Experiment
5r1
I II III IV
15.20 14.94 15.07 15.09
13.56 13.83 14.16 14.26
13.18 13.51 13.84 13.94
14.00a 14.01b
I II III IV
17.65 17.28 17.41 17.43
16.94 16.76 17.07 17.19
16.84 16.66 16.95 17.07
I II III IV
21.76 21.84 21.88 21.89
19.13 19.50 19.78 19.91
18.90 19.27 19.59 19.65
1p1
4r1
a b
16.90a ;b
19.60a 19.70b
Ref. [34]. Ref. [35].
the variation in the basis sets indicate the gradual progress to the saturation limit. Only the very high lying virtual orbitals were kept frozen in the post HF calculations. Inclusion of triples lowers the correlated energy of the ground state by 0.60 eV for N2 and by 0.57 eV for CO in case of basis IV. The lowering is considerable and unless a proper consideration of the triples in the ground state correlation is made effective, the IPs would have been underestimated. Since the VU-MRCC method, from the viewpoint of neither the state energies nor the differ-
ence energies, is variational, when employing the truncations in the CC operators may overestimate or underestimate the result. The variations of IPs in Tables 2 and 3 with various basis sets clearly suggest that the EIP-MRCCSD calculations overestimate all outervalence IPs (as we have calculated only them), of both molecules as compared to the EIP-MRCCSDT method. One might argue over the basis of deciding whether a result has been overestimated or underestimated. Seemingly though, the IPs are better with EIP-MRCCSD than with EIP-MRCCSDT when the basis set is way
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
off from the saturation limit. But this is only to be regarded as fortuitous if we examine the entire table for IPs. We put forth the reason as follows. At the exact level, there is no truncation which in principle means that the starting atomic orbital (AO) basis set and the Slater determinantal set are both infinite dimensional and hence are complete sets. And only then, the usual expansion theorem can be implemented for having a proper wavefunction by the linear combination of all the possible Slater determinants subject to spatial and spin symmetry requirements. But for practical reasons, we are compelled to impose truncation in the AO basis set as well as the set of Slater determinants. However, the basis set saturation limit attained by the saturation in the total HF ground state energies as well as the individual orbital energy on augmenting the AO basis gradually from a starting one assures us that an infinite dimensional AO basis is not essential. Rather a finite dimensional, but sufficiently large, AO basis can serve reliably. To our satisfaction, the matter with AO basis is now well settled on the basis of extensive state-of-the-art computational experiments. But in case of determinants, building a truncation scheme is more difficult as every theoretical model has its own built-up and hence its own advantages and limitations. Certainly, a truncation scheme that obliterates all determinants corresponding to a certain kind of higher excitation, such as triples, may effect the result significantly. To our knowledge, the only calculation of IPs of N2 and CO with systematic inclusion of higher excitations was due to Nakatsuji and his group [28]. They showed how the inclusion of n-tuple excitations (n > 2) subject to satisfying an energy threshold gave IPs for N2 and CO that are much nearer to the full CI ones than the calculations using excitations for (n 6 2) only. They did not compare their results with experiment because of the restrictive nature of the basis and the active space of MOs for the construction of Slater determinants they took. The same feature was also reflected in our case. Thus, the EIP-MRCCSDT calculation gave IPs all (except for one case, 1p1 of CO) being within 0.1 eV to experiment as compared to those due to EIPMRCCSD calculation which are off the experiments by from a low of 0.14 eV to a high of 0.33
305
eV. Even for 1p1 of CO, which is not within 0.1 eV from experiment, the EIP-MRCCSDT value is yet improved than the EIP-MRCCSD by 0.12 eV. Lastly, we emphasized that, although here we have referred to only the outer valence IPs, similar kind of observation is expected for the main IPs in the inner valence region too.
4. Summarising remarks The traditional VU-MRCC method [7,8] in vogue suffers seriously from the intruder state problem. This creeps in the effective hamiltonian based strategies since all the roots obtained by diagonalising the hamiltonian matrix are constrained to represent the eigenvalues of the hamiltonian, in spite of the fact that some of the associated functions are poorly represented by the functions spanning the model space. Our EIPMRCC approach opens up the possibility to improve the convergence behaviour of the traditional VU-MRCC equations even in the presence of intruders. In this Letter, we have presented the EIP-MRCC method by including full effect of three-body cluster operators along with the singles and doubles to calculate the vertical IPs of two typically interesting molecular systems, viz. N2 and CO, in an attempt to demonstrate the effect of higher-body cluster operator corrections for IP calculations using different basis sets. These systems, when studied in relatively smaller bases, show worsening of IP values upon the inclusion of triples. However, on reaching the limit of basis set saturation, the usual trend is observed, as is reflected by our results. We observe that our calculated IPs differ from the experimental values by less than 0.1 eV. Our inclusion of full three-electron cluster operators thus, in the language of CCSDT1-a [3], leads to a balanced description for vertical IP calculations.
Acknowledgements The authors gratefully acknowledge some very suggestive discussions with Professor Debashis Mukherjee. A part of the computation has been
306
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
carried out using the ICOSER (DST) facility of the level III computer center, IACS. We thank Professor S. Ray of the computer center for help and advice. Assistance from Mrs. Sarbani Saha is sincerely acknowledged. A substantial part of the work has been done in the Department of Chemistry, Tripura University, under the CSIR Project (01(1498)/97/EMR-II), New Delhi, India. P.G. acknowledges the CSIR for his fellowship. References [1] J. Cizek, Adv. Chem. Phys. 14 (1969) 35. [2] R.J. Bartlett, in: D.R. Yarkony (Ed.), Modern Electronic Structure Theory, World Scientific, Singapore, 1995. [3] Y.S. Lee, S.A. Kucharskii, R.J. Bartlett, J. Chem. Phys. 81 (1984) 5906. [4] P. Piecuch, S.A. Kucharski, R.J. Bartlett, J. Chem. Phys. 110 (1999) 6103, and references therein. [5] X. Li, J. Paldus, J. Chem. Phys. 107 (1997) 6257. [6] D. Mukherjee, R.K. Moitra, A. Mukhopadhyay, Mol. Phys. 33 (1977) 955. [7] I. Lindgren, Int. J. Quantum Chem. S12 (1978) 33. [8] D. Mukherjee, S. Pal, Adv. Quantum Chem. 20 (1989) 561, and references therein. [9] B. Jeziorski, H.J. Monkhorst, Phys. Rev. A 24 (1981) 1668. [10] J. Paldus, in: S. Wilson, G.H.F. Diercksen (Eds.), Methods in Computational Molecular Physics NATO, vol. 293, Plenum, New York, 1992, p. 99. [11] H. Nakatsuji, K. Hirao, J. Chem. Phys. 68 (1978) 2953, 4279. [12] K. Hirao, J. Chem. Phys. 79 (1983) 5000. [13] D. Mukhopadhyay, S. Mukhopadhyay, R. Chaudhuri, D. Mukherjee, Theor. Chim. Acta 80 (1991) 441.
[14] M. Nooijen, R.J. Bartlett, J. Chem. Phys. 102 (1995) 3629. [15] T.H. Schucan, H.A. Weidenm€ uller, Ann. Phys. 73 (1972) 108. [16] D. Mukhopadhyay, D. Mukherjee, Chem. Phys. Lett. 163 (1989) 3827. [17] L. Meissner, R.J. Bartlett, J. Chem. Phys. 92 (1990) 561. [18] J.P. Malrieu, P. Durand, J.P. Daudey, J. Phys. A (Math. Gen.) 18 (1985) 809. [19] D. Sinha, S.K. Mukhopadhyay, R. Chaudhuri, D. Mukherjee, Chem. Phys. Lett. 154 (1989) 544. [20] S. Rettrup, J. Comput. Phys. 45 (1982) 100. [21] K. Hirao, H. Nakatsuji, J. Comput. Phys. 45 (1982) 246. [22] B. Datta, R. Chaudhuri, D. Mukherjee, J. Mol. Struct. (Theochem) 361 (1996) 21. [23] U.S. Mahapatra, B. Datta, D. Mukherjee, J. Chem. Phys. 110 (1999) 6171, and references therein. [24] A. Haque, U. Kaldor, Chem. Phys. Lett. 120 (1985) 261. [25] N. Vaval, S. Pal, D. Mukherjee, Theor. Chem. Acc. 99 (1998) 100. [26] D.E. Bernholdt, R.J. Bartlett, Adv. Quantum Chem. 34 (1999) 271. [27] J.A. Pople, P. von R. Schleyer, J. Kaneti, G.W. Spitznagel, Chem. Phys. Lett. 145 (1988) 359. [28] M. Ehara, H. Nakatsuji, Chem. Phys. Lett. 282 (1998) 347. [29] H. Koch, O. Christiansen, P. Jørgensen, A.S. de Meras, T. Helgaker, J. Chem. Phys. 106 (1997) 1808. [30] J.D. Watts, R.J. Bartlett, Chem. Phys. Lett. 258 (1996) 581. [31] K. Raghavachari, J. Chem. Phys. 82 (1985) 4142. [32] S. Ghosh, D. Mukherjee, S.N. Bhattacharyya, Mol. Phys. 43 (1981) 173. [33] K.H. Tan, E.C. Brion, Ph.E. Vanderleeuw, M.J. Vander Wiel, Chem. Phys. 29 (1978) 299. [34] U. Gellius, E. Basilier, S. Svensson, T. Bergmark, K. Siegbahn, J. Electron. Spectrosc. 2 (1974) 405. [35] S. Dey, A.T. Dixon, K.R. Lassey, I.E. McCarthy, P.J.O. Tuebner, E. Weigold, P.S. Bagus, E.K. Viinikka, Phys. Rev. A 15 (1977) 102.
Chemical Physics Letters 361 (2002) 298–306 www.elsevier.com/locate/cplett
Full effect of triples in a valence universal multi-reference coupled cluster calculation Sudip Chattopadhyay a, Ashis Mitra b, Debasis Jana a, Pradipta Ghosh a, Dhiman Sinha c,d,* a
Department of Physical Chemistry, Indian Association for the Cultivation of Science, Calcutta 700 032, India b Department of Chemistry, R.K. Mahavidyalaya, Kailasahar, Tripura 799 277, India c Department of Pure Chemistry, University College of Science, 92 A P C Road, Calcutta 700 009, India d Department of Chemistry, Tripura University, Suryamaninagar, Tripura 799 130, India Received 2 April 2002; in final form 10 June 2002 This article is dedicated to Professor Mihir Chowdhury on the happy occassion of his 65th birthday.
Abstract Full three-electron coupled cluster (CC) operators (triples) along with all singles and doubles are included in a valence universal multi-reference CC calculation to estimate the change in the correlation effect for ionization potential (IP) calculations. We study the vertical IPs of two typical systems, N2 and CO, which go in the opposite direction, i.e., the results were worsened on inclusion of triples when working with relatively smaller basis. However, upon reaching the basis set saturation the expected trend in the results is observed and they are close within 0.1 eV to experimental values almost in every IP. Ó 2002 Published by Elsevier Science B.V.
1. Introduction The single-reference (SR) coupled-cluster (CC) method in varying degrees of truncation is now well established as a very successful tool for treating electron correlation of closed-shell systems to a high degree of accuracy [1,2]. This is size-extensive at any level of truncation of the rank of the cluster operator. However, the SRCC theory, truncated at the singles and doubles (SD) excited cluster level is not adequate
*
Corresponding author. E-mail address: dhimansinha@rediffmail.com (D. Sinha).
for the states with quasi-degeneracy emerging from multi-reference (MR) character. There are several ways to overcome this difficulty. One of them being an explicit accounting of three- and/ or four-body cluster amplitudes [3,4] within the SRCC framework. Another alternative scheme estimates the higher-body cluster amplitudes by augmenting the SR-CCSD equations by inclusion of these cluster amplitudes from a complete active space function [5]. The numerical implementation along this line is computationally demanding. Rather, conceptually the most appropriate one would be to treat quasi-degeneracy stemming from non-dynamical correlation effects by starting from an active space of quasi-
0009-2614/02/$ - see front matter Ó 2002 Published by Elsevier Science B.V. PII: S 0 0 0 9 - 2 6 1 4 ( 0 2 ) 0 0 9 7 4 - 0
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
degenerate determinants, while the remaining dynamical correlations can be included via a cluster expansion of the wavefunction involving virtual determinants. This approach is thus a pure MR generalization of SRCC method. The early multi-reference CC (MRCC) formulations were all based on effective hamiltonian spanned by the reference determinants, constituting a complete model space. There are two broad categories of such MRCC methods, viz. the valence universal (VU), or the Fock space [6–8] and the state universal (SU), or the Hilbert space [9,10] ones. While the SU-MRCC method is more suited for the studies of potential energy surfaces, the VU-MRCC approach is convenient to provide energy differences of spectroscopic interest (IPs, EAs and EEs) directly relative to some common reference state. There exists another suite of methods of direct calculation of energy differences starting from the CC or CC like representation of the ground state. These are SAC-CI [11,12], CC-based linear response theory (LRT) [13], and the CC-based equation-of-motion (EOM) [14] methods. While all these three methods are equivalent, they are generally different from the corresponding VU-MRCC methods. However, it has been shown that, for the principal IPs and EAs the VU-MRCC [7,8] is equivalent to SAC-CI [11,12], LRT [13] and EOM [14] IP and EA CC methods. In this Letter we focus on the IP problem via the VUMRCC method. The VU-MRCC, despite its formal rigor, is seriously beset by the problem of intruder states [15]. However, implementation of incomplete model spaces [16,17] might seem to be a possible way out, but it is not so in practice. This is because, during the generation of energy differences of spectroscopic interest, it is appropriate to provide at least a qualitatively correct description of the high-lying states perturbed by the intruder states as they form a part of the same manifold of excited states, rather than ignoring them completely from the description. The concept of intermediate hamiltonian [18] offers a real possibility of avoiding the determinantal effects of intruders while preserving most of the advantages of the effective hamiltonian framework.
299
Mukherjee and co-workers [13,19] had earlier proposed a technique to improve the convergence behaviour of the VU-MRCC equations [7,8] in the presence of intruders using the idea of eigenvalueindependent partitioning (henceforth to be referred as EIP-MRCC). The EIP-MRCC does not differ from the traditional VU-MRCC method as far as the basic theoretical constructs, viz. the definition of the valence rank, the hierarchical decoupling of the various valence rank problems, etc., are concerned. It nonetheless provides an alternative avenue for solving the VU-MRCC equations having certain computationally strategic advantages. The EIP technique converts the nonlinear VU-MRCC equations for any model space into a set of non-hermitian eigenvalue equations. This opens up the interesting possibility of arriving at the desired solutions through the use of stable root-search and root-homing procedures used in eigenproblems [20,21] to determine the desired low-lying roots of spectroscopic interest one by one iteratively, preserving the norm of the corresponding eigenvector. Thus the method is less prone to divergence than the corresponding nonlinear equations of traditional VU-MRCC. It is noteworthy that the intermediate hamiltonian method via EIP in the context of CC framework have also been developed by Mukherjee and coworkers [22]. The recently formulated state-specific (SS) method [23] is the most effective strategy to tackle the intruder state problem, while retaining most of the important traits of a good many-body model. Pertinent to mention is that the SS methods based on the Fock space MRCC approach still remains an uncharted area of electronic structure theory. As noted, unlike the IP calculation employing the SRCC methods, which require separate calculations for each state, the VU-MRCC approach provides results for many states at once which are numerically at par with other methods. It is worthwhile to point out one of the major distinguishing features of the VU-MRCC approach from the SR-based strategies. This is in regard to the excitation operators in the two formulations – certain double excitations in the SR method only appear as triple excitations in the VU-MRCC. Thus it would be interesting to explore the possi-
300
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
bility of developing VU-MRCC schemes including higher-body cluster operators (triples and quadruples) and compare their performance with other simple VU-MRCC methods and the SR calculations. Haque and Kaldor [24], Mukherjee and co-workers [25] and Bartlett and co-workers [26] have included the lowest order effects of triples in their VU-MRCC formulations with some demonstrative results. Also, some methods like the complete fourth order Møller–Plesset perturbation theory (MP4) calculation [27], and the SAC-CI (general-R) method [28], CC based-LRT [29] and EOM [30] have taken into account the triples/ higher excitations extensively. In the present Letter, we report the effects of the inclusion of triples to our EIP-MRCC approach and examine the performance of this by some pilot IP calculations. In this application for the excited states, we have considered the full set of determinants generated by the doubles and the triples, below a certain threshold energy. However, in the ground state, which is a closed-shell, the contribution of triples are considered through the term VN T2 =D (see [3] for a detailed discussion) only, where VN and T2 are, respectively, the two-body part of the hamiltonian and two-body closed-shell CC operators, both in normal order, and D is the accompanying orbital energy difference. An important feature of the inclusion of triples is that the dimension of the determinantal set increases tremendously, by several tens of thousands or even much more, as against the set of only a few hundred determinants appearing while using no triples at all. We have observed that a careless truncation in the set of determinants may tamper the results significantly. In this Letter we have calculated all the vertical IPs of the outer valence region of N2 and CO using four different basis sets each, augmented gradually towards the saturation limit. Our aim is to see, to what extent the full onset of triples in this method enhances the CC result towards experiment. Interestingly, relatively smaller basis sets gave apparently poorer IPs for both the systems while including triples in the calculation. However, upon achieving a basis set saturation, we could obtain consistently better outer valence IPs for the two systems. A similar trend has also been observed by
Raghavachari [31] in regard to the effect of basis set and electron correlation effects on electron affinities ðEAÞ of some first row atoms. They found that with relatively smaller basis, second order Møller–Plesset perturbation (MP2) results for EA were sometimes, e.g., for O and F, much better than the complete MP4 one, even though the latter included higher excitations extensively. But the MP4 results turned out to be consistent for all the atoms they had taken when using the bigger basis sets. We conclude this section by putting forth the organization of this Letter. In the section to follow, we present a succinct discussion on the EIP-MRCC theory. In Section 3, we discuss the algorithmic consideration that leads to the numerical implementation of our method. In this section itself we provide the vertical IPs of two interesting systems (N2 and CO) as our illustrative application to judge the efficacy of our EIPMRCC method. We conclude the Letter by summarizing remarks in Section 4.
2. Theory We briefly describe here the mathematical framework of our EIP-MRCC method. One starts from the VU Bloch equation [7,8] H XP ðk;lÞ ¼ XHeff P ðk;lÞ ; 0 6 k 6 m; 0 6 l 6 n;
ð1Þ
where H is the atomic or molecular hamiltonian in normal order and P ðk;lÞ is the projector for (k holes and l particles), or simply (k; l), model space [7,8]. We shall explain the effective hamiltonian Heff shortly. The wave operator, X is expressed as a VU normal ordered exponential ansatz [7,8]: X ¼ exp S ½m;n ( !) m X n X ðk;lÞ ¼ exp S ; ð2Þ k¼0
l¼0
S ðk;lÞ being the (k; l) valence component [7,8,19] of the cluster operator S ½m;n , where [m; n] means the collection of all (k; l) valence components upto the highest (m; n) valence component. S ðk;lÞ acts upon P ðk;lÞ space determinants to generate the virtual Q
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
space and R space determinants (we sub-divide the entire virtual space into two sub-spaces, viz. Q and R for reasons cited below). The VU-MRCC method calculates ionization potentials as direct difference energies instead of calculating the state energies Egr ðN Þ and EðN 1Þ individually to find out their difference, where N stands for the total number of electrons in the closed-shell system. In the approximate level calculation, that is, when we apply truncation scheme, the methods that calculate difference energies directly do a balanced cancellation of the common correlation terms as compared to the methods that calculate state energies separately. In the VU-MRCC framework, the IP calculation is called a [1,0] valence problem, and thus it includes the CC operators for (0,0) and (1,0) valence components. Consequently, X ¼ exp S ð0;0Þ þ S ð1;0Þ ð3Þ ¼ Xð0;0Þ exp S ð1;0Þ ; where Xð0;0Þ ¼ exp S ð0;0Þ , is the closed-shell wave operator, S ð0;0Þ stands for the cluster operators of the closed-shell CC calculation, and S ð1;0Þ is the cluster operator for the single ionization with simultaneous excitations of that system. The standard choices of S ð0;0Þ and S ð1;0Þ are ð0;0Þ
S ð0;0Þ ¼ S1
ð0;0Þ
þ S2
ð1;0Þ
and S ð1;0Þ ¼ S2
;
ð4Þ
where the maximum electron rank at both the ð0;0Þ closed-shell and open-shell situations is two, S1 ð0;0Þ and S2 stand for singles and doubles in the closed-shell situation. If all the hole orbitals are kept active, which has been the case in this applið1;0Þ cation, no singles S1 (Fig. 1a) should figure out in the calculation. Then, at the present level of ð1;0Þ approximation, only the doubles S2 (Fig. 1b) act upon the P -subspace that contains 1h determinants only to generate the 2h 1p determinants of the Q-subspace.
Following the procedures given in [19], one gets, on using Eq. (4), the EIP-MRCCSD matrix equation (Method I) (S and D stand for singles and doubles, respectively):
XPP YPP YPQ YQP
¼
YQQ XPP XQP
ð1;0Þ
ð1;0Þ
Fig. 1. (a) Singles, S1 , and (b) doubles, S2 erators for open-shell.
of the CC op-
MM
ðND þMÞM
ð1;0Þ
ð1;0Þ
Heff XPP ¼ XPP EPP ; where
ð5Þ
ð1;0Þ Heff
ð6Þ ð0;0Þ
¼ P ð1;0Þ H X
P ð1;0Þ , and
Rð2Þ XPP ¼ XQP or 1 Rð2Þ ¼ XQP XPP ;
ð7Þ ð1;0Þ S2
ð2Þ
where R represents the matrix of the amplið1;0Þ tudes. EPP are actually the IPs as difference energies. M and ND are the dimensions of the P ð1;0Þ subspace and the Q-subspace, respectively. To incorporate the effect of triples in our EIPMRCC theory, we choose the following scheme: (a) ð0;0Þ the three-body zero-valence operator, S3 is only approximately computed, as in the CCSDT1-a truncation [3], this corresponds to a non-iterative ð0;0Þ inclusion of S3 ; (b) the effect of the one-valence ð1;0Þ three-body operator, S3 , is, however, included ð0;0Þ exactly. Thus, S and S ð1;0Þ in our scheme are ð0;0Þ
þ S2
ð1;0Þ
þ S3
S ð0;0Þ ¼ S1 ð1;0Þ S3
(b)
XQP ðND þMÞM h i ð1;0Þ EPP :
ðND þMÞðND þMÞ
It is interesting to note that Eq. (5) is identical in form to CC-based LRT [13,32] equations for IP calculations using 1h and 2h 1p determinants. The terms in Eq. (5) are explained as follows. YPP , YPQ , YQP and YQQ are derived from the conð0;0Þ tracted composite ðH X Þ obtained from the left hand side of Eq. (1), after doing some manipulations guided by Wick’s algebra on it. The matrix X relates to the following equations [19],
¼ S2
(a)
301
ð0;0Þ ð1;0Þ
ð0;0Þ
þ S3 :
and S ð1;0Þ ð8Þ
shown in Fig. 2 corresponds to double shakeup excitations (Fig. 2a) and also correlation relaxation (Fig. 2b). In (Fig. 2b), the lines with double arrow indicate spectator scattering. We calculate the full effect of these, which is the key
302
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
(a)
(b) ð1;0Þ
Fig. 2. The various types of triples, S3 for open-shell.
of the CC operators
issue of the present work. In the figures the lines having single arrow pointing from left to right indicate the inactive hole lines, those from right to left are the inactive particle lines. The lines having double arrow from left to right are active hole lines. By similar induction as that of EIP-MRCCSD scheme, using Eq. (8) we get the corresponding EIP-MRCCSDT matrix equation as (Method II): 2 3 2 3 YPP YPQ YPR XPP 4 YQP YQQ YQR 5 4 XQP 5 YRP YRQ YRR ðNDT þMÞðNDT þMÞ XRP ðNDT þMÞM 2 3 XPP h i ð1;0Þ ¼ 4 XQP 5 EPP ; ð9Þ MM XRP ðNDT þMÞM
Fig. 3. Representative diagrams for Y -matrix elements involving P - and Q-subspaces. The circle stands for the one-body part and the square stands for the two-body and the three-body ð0;0Þ parts of H X .
where T represents that triples were taken in full for the excited state, but for the ground state they were considered at the CCSDT1-a [3] level of approximation only. Like Eq. (7), we get a similar relation 1 Rð3Þ ¼ XRP XPP ;
ð10Þ ð1;0Þ
where Rð3Þ represents the matrix of S3 amplitudes. The R-subspace contains the (3h 2p) determinants. NDT is the sum of the dimensions of the Q- and R-subspaces together. The first few roots of the Y -matrices in Eqs. (5) and (9) are determined by Rettrup iterative algorithm [20]. The difference in the so obtained two similar roots show the change due to the inclusion of triples. An interesting point about Eqs. (5) and (9) is that we start ð1;0Þ ð1;0Þ with X which explicitly contains S2 and S3 , but we need not deal with them finally, although the information of all the orbitals on the left of S ð1;0Þ vertex is certainly necessary for locating the Y -matrix elements in the computer programme.
Fig. 4. Y -matrix elements involving R-subspace. The square in ð0;0Þ YRP stands for the three-body part of HX , and that in ð0;0Þ YRQ and YRR stand for the two- and three-body parts of H X .
This is described by the disconnected lines in Figs. 3 and 4. Of course, if required, one can at once determine S ð1;0Þ amplitudes through Eqs. (7) and (10).
3. Results and discussions In this Letter we present two separate sets of calculations done for IPs with each basis. In one, comprising the Method-I, the full CCSD calcula-
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
tion was followed by the calculation of the transformed Hamiltonian, H ¼ fH Xð0;0Þ g, by using the corresponding CC amplitudes, and finally all outer valence IPs were determined by direct diagonalization or by iterative diagonalization using the Rettrup algorithm [20], depending upon the size of the IP matrix. In the other, comprising the Method-II, first a CCSDT1-a [3] calculation was done to be followed by the calculation of H using the new CC amplitudes, and finally again the outer valence IPs were determined but now only by using the Rettrup algorithm as the size of the matrix is very large. In either of the methods, the maximum electron excitation rank of the CC correlation operators both at the closed-shell and the open-shell levels were kept the same, which was for estimating the electron correlations at the two levels on the same footing. For example, in Method-I, the maximum electron rank at both the closed- and open-shells was two, whereas in Method-II, the rank was three. The CCSD and the CCSDT1-a correlation energies are given in Table 1. While constructing H in Method-II, the closedshell CC amplitudes for singles and doubles obtained after CCSDT1-a calculation were used. That is, no closed-shell CC amplitude related to triples was employed. However, the three-body H ð0;0Þ operators through VN S2 ð¼ VN T2 Þ were included. The convergence with the use of Rettrup algorithm was pretty fast and en route to convergence it took only 10–12 iterations provided that the starting iterate was properly chosen. Furthermore, by making the Rettrup FORTRAN code independent of input-output (I/O) operations, the wall-clock time is drastically reduced. Instead of reading the
303
stored Y -matrix elements from the hard disk at each iteration, they are calculated afresh, one row at a time, in core memory and simultaneously utilized. Although such a strategy increases the total CPU time, this is, however, very much less than the wall-clock time required while following the I/O dependent Rettrup code. All the virtual IPs of the outer valence region of N2 are compared with the experiment in Table 2, and the corresponding results for CO are presented in Table 3. We have performed the calculations employing four different basis sets (I, II, III and IV) for each molecule which are described below. The first basis set does not have any polarization function, whereas the latter ones include d, d and f , and d, f and g polarization functions, respectively. The d, f and g functions in the bases comprise of 6, 10 and 15 Cartesian gaussians, respectively. Basis (on each atom) I II III
IV
TZ (10s,6p)/[5s,3p] cc-pVDZ (9s,4p,1d)/ [3s,2p,1d] cc-pVTZ (10s,5p,2d,1f)/ [4s,3p,2d,1f] cc-pVQZ (12s,6p,3d,2f,1g)/ [5s,4p,3d,2f,1g]
All the basis sets were collected from: http:// www.emsl.pnl.gov:2080/form/basisform.html In Tables 2 and 3, the variation of the Hartree– Fock (HF) energies and the Koopmans IPs with
Table 1 Hartree–Fock (HF) energies and correlation energies (in a.u.) System
Basis
EHF
Corr. energy (CCSD)
Corr. energy (CCSDT-1a)
N2
I II III IV
)108.900582 )108.954759 )108.984150 )108.991433
)0.251158 )0.319201 )0.405526 )0.430544
)0.262531 )0.332081 )0.425578 )0.452686
CO
I II III IV
)112.707913 )112.749783 )112.781003 )112.789142
)0.231969 )0.303745 )0.388982 )0.432384
)0.245393 )0.316097 )0.407992 )0.453468
re for N2 2.074 a.u., for CO 2.13 a.u.
304
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
Table 2 Vertical ionization potential (in eV) of N2 States
Basis
Koopmans
EIP-MRCCSD
EIP-MRCCSDT
3r1 g
I II III IV
17.15 17.06 17.23 17.25
14.90 15.22 15.62 15.74
14.77 15.07 15.45 15.56
I II III IV
16.97 16.56 16.68 16.70
17.33 16.94 17.22 17.33
16.95 16.69 16.96 17.07
I II III IV
21.09 21.08 21.18 21.19
18.00 18.51 18.87 18.97
17.75 18.28 18.60 18.71
1p1 u
3r1 g
a
Experimenta 15.6
17.0
18.8
Ref. [33].
Table 3 Vertical ionization potential (in eV) of CO States
Basis
Koopmans
EIP-MRCCSD
EIP-MRCCSDT
Experiment
5r1
I II III IV
15.20 14.94 15.07 15.09
13.56 13.83 14.16 14.26
13.18 13.51 13.84 13.94
14.00a 14.01b
I II III IV
17.65 17.28 17.41 17.43
16.94 16.76 17.07 17.19
16.84 16.66 16.95 17.07
I II III IV
21.76 21.84 21.88 21.89
19.13 19.50 19.78 19.91
18.90 19.27 19.59 19.65
1p1
4r1
a b
16.90a ;b
19.60a 19.70b
Ref. [34]. Ref. [35].
the variation in the basis sets indicate the gradual progress to the saturation limit. Only the very high lying virtual orbitals were kept frozen in the post HF calculations. Inclusion of triples lowers the correlated energy of the ground state by 0.60 eV for N2 and by 0.57 eV for CO in case of basis IV. The lowering is considerable and unless a proper consideration of the triples in the ground state correlation is made effective, the IPs would have been underestimated. Since the VU-MRCC method, from the viewpoint of neither the state energies nor the differ-
ence energies, is variational, when employing the truncations in the CC operators may overestimate or underestimate the result. The variations of IPs in Tables 2 and 3 with various basis sets clearly suggest that the EIP-MRCCSD calculations overestimate all outervalence IPs (as we have calculated only them), of both molecules as compared to the EIP-MRCCSDT method. One might argue over the basis of deciding whether a result has been overestimated or underestimated. Seemingly though, the IPs are better with EIP-MRCCSD than with EIP-MRCCSDT when the basis set is way
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
off from the saturation limit. But this is only to be regarded as fortuitous if we examine the entire table for IPs. We put forth the reason as follows. At the exact level, there is no truncation which in principle means that the starting atomic orbital (AO) basis set and the Slater determinantal set are both infinite dimensional and hence are complete sets. And only then, the usual expansion theorem can be implemented for having a proper wavefunction by the linear combination of all the possible Slater determinants subject to spatial and spin symmetry requirements. But for practical reasons, we are compelled to impose truncation in the AO basis set as well as the set of Slater determinants. However, the basis set saturation limit attained by the saturation in the total HF ground state energies as well as the individual orbital energy on augmenting the AO basis gradually from a starting one assures us that an infinite dimensional AO basis is not essential. Rather a finite dimensional, but sufficiently large, AO basis can serve reliably. To our satisfaction, the matter with AO basis is now well settled on the basis of extensive state-of-the-art computational experiments. But in case of determinants, building a truncation scheme is more difficult as every theoretical model has its own built-up and hence its own advantages and limitations. Certainly, a truncation scheme that obliterates all determinants corresponding to a certain kind of higher excitation, such as triples, may effect the result significantly. To our knowledge, the only calculation of IPs of N2 and CO with systematic inclusion of higher excitations was due to Nakatsuji and his group [28]. They showed how the inclusion of n-tuple excitations (n > 2) subject to satisfying an energy threshold gave IPs for N2 and CO that are much nearer to the full CI ones than the calculations using excitations for (n 6 2) only. They did not compare their results with experiment because of the restrictive nature of the basis and the active space of MOs for the construction of Slater determinants they took. The same feature was also reflected in our case. Thus, the EIP-MRCCSDT calculation gave IPs all (except for one case, 1p1 of CO) being within 0.1 eV to experiment as compared to those due to EIPMRCCSD calculation which are off the experiments by from a low of 0.14 eV to a high of 0.33
305
eV. Even for 1p1 of CO, which is not within 0.1 eV from experiment, the EIP-MRCCSDT value is yet improved than the EIP-MRCCSD by 0.12 eV. Lastly, we emphasized that, although here we have referred to only the outer valence IPs, similar kind of observation is expected for the main IPs in the inner valence region too.
4. Summarising remarks The traditional VU-MRCC method [7,8] in vogue suffers seriously from the intruder state problem. This creeps in the effective hamiltonian based strategies since all the roots obtained by diagonalising the hamiltonian matrix are constrained to represent the eigenvalues of the hamiltonian, in spite of the fact that some of the associated functions are poorly represented by the functions spanning the model space. Our EIPMRCC approach opens up the possibility to improve the convergence behaviour of the traditional VU-MRCC equations even in the presence of intruders. In this Letter, we have presented the EIP-MRCC method by including full effect of three-body cluster operators along with the singles and doubles to calculate the vertical IPs of two typically interesting molecular systems, viz. N2 and CO, in an attempt to demonstrate the effect of higher-body cluster operator corrections for IP calculations using different basis sets. These systems, when studied in relatively smaller bases, show worsening of IP values upon the inclusion of triples. However, on reaching the limit of basis set saturation, the usual trend is observed, as is reflected by our results. We observe that our calculated IPs differ from the experimental values by less than 0.1 eV. Our inclusion of full three-electron cluster operators thus, in the language of CCSDT1-a [3], leads to a balanced description for vertical IP calculations.
Acknowledgements The authors gratefully acknowledge some very suggestive discussions with Professor Debashis Mukherjee. A part of the computation has been
306
S. Chattopadhyay et al. / Chemical Physics Letters 361 (2002) 298–306
carried out using the ICOSER (DST) facility of the level III computer center, IACS. We thank Professor S. Ray of the computer center for help and advice. Assistance from Mrs. Sarbani Saha is sincerely acknowledged. A substantial part of the work has been done in the Department of Chemistry, Tripura University, under the CSIR Project (01(1498)/97/EMR-II), New Delhi, India. P.G. acknowledges the CSIR for his fellowship. References [1] J. Cizek, Adv. Chem. Phys. 14 (1969) 35. [2] R.J. Bartlett, in: D.R. Yarkony (Ed.), Modern Electronic Structure Theory, World Scientific, Singapore, 1995. [3] Y.S. Lee, S.A. Kucharskii, R.J. Bartlett, J. Chem. Phys. 81 (1984) 5906. [4] P. Piecuch, S.A. Kucharski, R.J. Bartlett, J. Chem. Phys. 110 (1999) 6103, and references therein. [5] X. Li, J. Paldus, J. Chem. Phys. 107 (1997) 6257. [6] D. Mukherjee, R.K. Moitra, A. Mukhopadhyay, Mol. Phys. 33 (1977) 955. [7] I. Lindgren, Int. J. Quantum Chem. S12 (1978) 33. [8] D. Mukherjee, S. Pal, Adv. Quantum Chem. 20 (1989) 561, and references therein. [9] B. Jeziorski, H.J. Monkhorst, Phys. Rev. A 24 (1981) 1668. [10] J. Paldus, in: S. Wilson, G.H.F. Diercksen (Eds.), Methods in Computational Molecular Physics NATO, vol. 293, Plenum, New York, 1992, p. 99. [11] H. Nakatsuji, K. Hirao, J. Chem. Phys. 68 (1978) 2953, 4279. [12] K. Hirao, J. Chem. Phys. 79 (1983) 5000. [13] D. Mukhopadhyay, S. Mukhopadhyay, R. Chaudhuri, D. Mukherjee, Theor. Chim. Acta 80 (1991) 441.
[14] M. Nooijen, R.J. Bartlett, J. Chem. Phys. 102 (1995) 3629. [15] T.H. Schucan, H.A. Weidenm€ uller, Ann. Phys. 73 (1972) 108. [16] D. Mukhopadhyay, D. Mukherjee, Chem. Phys. Lett. 163 (1989) 3827. [17] L. Meissner, R.J. Bartlett, J. Chem. Phys. 92 (1990) 561. [18] J.P. Malrieu, P. Durand, J.P. Daudey, J. Phys. A (Math. Gen.) 18 (1985) 809. [19] D. Sinha, S.K. Mukhopadhyay, R. Chaudhuri, D. Mukherjee, Chem. Phys. Lett. 154 (1989) 544. [20] S. Rettrup, J. Comput. Phys. 45 (1982) 100. [21] K. Hirao, H. Nakatsuji, J. Comput. Phys. 45 (1982) 246. [22] B. Datta, R. Chaudhuri, D. Mukherjee, J. Mol. Struct. (Theochem) 361 (1996) 21. [23] U.S. Mahapatra, B. Datta, D. Mukherjee, J. Chem. Phys. 110 (1999) 6171, and references therein. [24] A. Haque, U. Kaldor, Chem. Phys. Lett. 120 (1985) 261. [25] N. Vaval, S. Pal, D. Mukherjee, Theor. Chem. Acc. 99 (1998) 100. [26] D.E. Bernholdt, R.J. Bartlett, Adv. Quantum Chem. 34 (1999) 271. [27] J.A. Pople, P. von R. Schleyer, J. Kaneti, G.W. Spitznagel, Chem. Phys. Lett. 145 (1988) 359. [28] M. Ehara, H. Nakatsuji, Chem. Phys. Lett. 282 (1998) 347. [29] H. Koch, O. Christiansen, P. Jørgensen, A.S. de Meras, T. Helgaker, J. Chem. Phys. 106 (1997) 1808. [30] J.D. Watts, R.J. Bartlett, Chem. Phys. Lett. 258 (1996) 581. [31] K. Raghavachari, J. Chem. Phys. 82 (1985) 4142. [32] S. Ghosh, D. Mukherjee, S.N. Bhattacharyya, Mol. Phys. 43 (1981) 173. [33] K.H. Tan, E.C. Brion, Ph.E. Vanderleeuw, M.J. Vander Wiel, Chem. Phys. 29 (1978) 299. [34] U. Gellius, E. Basilier, S. Svensson, T. Bergmark, K. Siegbahn, J. Electron. Spectrosc. 2 (1974) 405. [35] S. Dey, A.T. Dixon, K.R. Lassey, I.E. McCarthy, P.J.O. Tuebner, E. Weigold, P.S. Bagus, E.K. Viinikka, Phys. Rev. A 15 (1977) 102.