Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://classic.chem.msu.su/gran/gamess/JChemPhys_134_214113.pdf Äàòà èçìåíåíèÿ: Tue Jun 7 22:04:58 2011 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:08:33 2012 Êîäèðîâêà:

THE JOURNAL OF CHEMICAL PHYSICS 134, 214113 (2011)

Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory
Alexander A. Granovskya)
Firefly Project, Moscow 117593, Russia

(Received 23 October 2010; accepted 13 May 2011; published online 7 June 2011) The distinctive desirable features, both mathematically and physically meaningful, for all partially contracted multi-state multi-reference perturbation theories (MS-MR-PT) are explicitly formulated. The original approach to MS-MR-PT theory, called extended multi-configuration quasi-degenerate perturbation theory (XMCQDPT), having most, if not all, of the desirable properties is introduced. The new method is applied at the second order of perturbation theory (XMCQDPT2) to the 11 A ­ 21 A conical intersection in allene molecule, the avoided crossing in LiF molecule, and the 11 A1 to 21 A1 electronic transition in cis-1,3-butadiene. The new theory has several advantages compared to those of well-established approaches, such as second order multi-configuration quasi-degenerate perturbation theory and multi-state-second order complete active space perturbation theory. The analysis of the prevalent approaches to the MS-MR-PT theory performed within the framework of the XMCQDPT theory unveils the origin of their common inherent problems. We describe the efficient implementation strategy that makes XMCQDPT2 an especially useful general-purpose tool in the high-level modeling of small to large molecular systems. © 2011 American Institute of Physics. [doi:10.1063/1.3596699]
I. INTRODUCTION

Perturbation Theory (PT) has always been one of the most successful yet relatively inexpensive tools of quantum chemistry. In particular, the variant of PT suggested by MÜller and Plesset1 (MPn) turned out to be a very powerful tool already at the second order (MP2).2, 3 However, it is well known that simple single-reference constructs, which the MPn is, work very well for electronic states dominated by a single Slater determinant, but fail severely for multi-configuration wavefunctions. Several successful generalizations of MP or MP-like approaches to the case of multi-configuration reference states (MR-PT) have been suggested during the last two decades. Among them are the MR-MP2 approach by Hirao,4­7 second order complete active space perturbation theory (CASPT2) approach by Andersson, Malmqvist, and Roos,8­11 and second order n-electron valence state perturbation theory (NEVPT2) by Angeli et al.,12­15 the latter being a relatively recent development in this field. Overview of several other approaches to MR-PT can be found in Ref. 16. Nevertheless, the single-state approaches to MR-PT assume the "diagonalize then perturb" philosophy and require a zero-order wavefunction that is already well described by complete active space self-consistent field (CASSCF) or a similar procedure. In practice, this is more the exception than the rule; so the possibility for several reference states to mix within the MR-PT treatment is very important. The natural way to allow this mixing is to use the multi-state (MS) formulations of MR-PT approaches leading to various MS-MR-PT theories. Specifically, the MS formulation of MR-MP2 has long been known as the multi-configuration
a) Electronic mail: alex.granovsky@gmail.com.

quasi-degenerate perturbation theory (MCQDPT) approach developed by Nakano.17­23 The MS formulation of CASPT2 known as MS-CASPT2 was originally suggested by Finley et al.,24 while the MS version of NEVPT2 known as QDNEVPT2 was recently developed by Angeli et al.25 All of them are of the "diagonalize then perturb then diagonalize" type. In particular, this means that instead of perturbing the entire CASCI Hamiltonian, only several selected CI roots are perturbed and mixed, i.e., these methods employ the so-called partial contraction in the space of CI expansion coefficients. These approaches are more or less directly related to the formal theory of the effective Hamiltonians and quasi-degenerate perturbation theories (QDPT) that have been well established for a long time.26­37 Essentially, all three theories define approximations to the exact effective Hamiltonian acting within a model space,33, 37 i.e., a subspace spanned by selected CI roots. Energies of perturbed states are then obtained as eigenvalues of effective Hamiltonian while projections of perturbed states onto zero-order states are defined by the corresponding eigenvectors. It is worth to mention the SS-MRPT approach by Mukherjee and co-workers as an interesting state-specific alternative to described traditional multi-state scheme yet allowing relaxation of zero-order wavefunction within a model space.38­41 The most important distinctions in formulation of various MS-MR-PT schemes are the selection of zero-order Hamiltonian and of so called "perturbers," i.e., zero-order states allowed to interact with zero-order wavefunctions in PT treatment. In the case of MCQDPT2, "perturbers" are simply all possible individual CSFs or Slater determinants not belonging to CASCI space and obtained applying single and double excitations individually to every configuration state function (CSF) or determinant entering reference states. In
© 2011 American Institute of Physics

0021-9606/2011/134(21)/214113/14/$30.00

134, 214113-1

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-2

Alexander A. Granovsky

J. Chem. Phys. 134, 214113 (2011)

contrary, both MS-CASPT2 and QD-NEVPT2 are formulated using state-specific and more sophisticated "perturbers" defined as specific linear combinations of some selected classes of CSFs or determinants. Depending on the particular scheme in use, the latter approach is known as either internal8, 9, 42 or external43, 44 contraction. Both MR-MP2 and CASPT2 approaches apply oneparticle Fock-like zero-order Hamiltonian. In contract, NEVPT2 scheme adopts partially bi-electronic zero-order Hamiltonian initially suggested by Dyall.45 As will be discussed later, the situation is more complicated for the MS-MR-PT counterparts of all three approaches. The goal of the present paper is the further development in the field of MS-MR-PT theories. Namely, we present a new approach to MS-MR-PT and analyze some of the existing approaches within its framework. The rest of the paper is organized as follows. We start with formulation of some important properties that seem to be desirable for any partially contracted MS-MR-PTbased approach. Then we present a brief critical overview of MCQDPT and formulate the new approach to MS-MR-PT, which we call extended multi-configuration quasi-degenerate perturbation theory (XMCQDPT). We demonstrate that XMCQDPT has the desired properties for an improved, accurate and effective MS-MR-PT theory. We also present the efficient formulation of the working equations of XMCQDPT2 (XMCQDPT at second order). Next, we apply XMCQDPT2, as well as MCQDPT2 and different formulations of MSCASPT2, to model problems including a conical intersection in the allene molecule, an avoided crossing in LiF molecule, and a characterization of the 11 A1 to 21 A1 electronic transition in cis-1,3-butadiene. We demonstrate the superiority of the XMCQDPT2 methodology and the importance of the requirements we imposed on XMCQDPT. Discussion and final remarks conclude the paper.

II. PRELIMINARY CONSIDERATIONS

Let us formulate important properties that are desirable for any partially contracted MS-MR-PT-based approach and which we tried to achieve developing XMCQDPT. These properties can be considered as an attempt to bring various formulations of MS-MR-PT into better accordance with the formal theory of Effective Hamiltonians. In practical applications, various MR-PT and especially MS-MR-PT theories can be plagued by the so-called "intruder state problem,"46, 47 causing appearance of very small energy denominators in PT series and leading to spurious results of the entire PT calculation. The persistence of this problem varies with particular theory. For instance, it has been demonstrated that NEVPT2 is inherently less prone to this problem as compared with MR-MP2 and CASPT2.48, 49 However, the problem still probably persists in the QD version of NEVPT2 approach on the par with other MS-MR-PT schemes. For simplicity, we assume below that intruder states are either irrelevant, or that they were successfully eliminated using, for instance, intruder state avoidance (ISA) (Ref. 50) or similar techniques.

First, we assume that the effective Hamiltonian should be explicitly dependent on the dimension of the model space in any nontrivial order of PT. Evidently, this dependence should not be just the trivial one caused by the simple extension of the model space. It also should not be related only to the variations of orbital energies, orbitals, or both, caused by the averaging of one-particle density matrices over varying number of CI roots, and hence should not be related only to the implicit dependence of the effective Fock operator(s) and zero-order Hamiltonian on the dimension of the model space. Indeed, the exact effective Hamiltonian fits the explicit dependence of the true effective Hamiltonians on the model space extension. Hence, approximations to the exact effective Hamiltonian should approximate this explicit dependence on the model space size as well. One should realize that, in general, the restriction of an effective Hamiltonian to a subspace of the model space is not a new effective Hamiltonian associated with this smaller subspace. A proof of this statement based on the definition of effective Hamiltonians and variational principle can be found in Appendix SI of supplementary material.51 Second, the very important property is the convergence of energies as well as other observable properties with respect to model space extension at any fixed order of PT provided that the active space used in underlying MCSCF calculations is well-balanced and accurately designed so that there are no spurious zero-order states (i.e., solutions of non-linear MCSCF equations which do not describe real excited states52­54 ). This property is evidently mandatory for the low-lying electronic states that are typically of interest and is desirable for other states as well. The third property is that the effective Hamiltonian should be a function of the subspace spanned by the selected CI vectors, rather than a function of any particular choice of basis in this subspace. In particular, this means that the effective Hamiltonian should be the same regardless of whether the initial CI vectors, or any alternative vectors obtained as their arbitrary orthogonal transformation, were used for MS-MR-PT computations. Indeed, as many-body perturbation theory (MBPT) is a purely algebraic approach, it has well-defined algebraic invariance properties. Hence, it is natural to impose the same restriction on the MS-MR-PT. The very important consequence of this property is that the correctly formulated partially contracted theory should be exactly equivalent to the corresponding uncontracted theory provided that the dimension of the model space is equal to the dimension of the CI space (i.e., the overall number of CSFs or determinants spanning the CI space). Finally, the most important and physically evident requirement is that the computed energies must be uniquely defined, continuous and smooth functions of the molecular geometry and any other external parameters, with possible exceptions at the manifolds of their accidental degeneracy such as conical intersections. In modeling of many-electron molecular systems, the correct dependence of computed energies and other properties on the number of electrons in the correlation treatment plays crucial role. This requirement is related to such features of the underlying theories as exact or approximate sizeconsistency55 and separability.56­59 It is often believed that

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-3

XMCQDPT: The new approach to MS-MR-PT

J. Chem. Phys. 134, 214113 (2011)

the exact size-consistency is very desirable for any high-level computational approach designed to describe large systems. However, we do not consider small deviations from the exact size-consistency as a serious deficiency of MS-MR-PT theories, provided that these deviations are more or less constant regardless of the model space dimension and do not cause any significant unphysical effects even for large systems. It seems that the exact or approximate core-separability59 is the most important requirement, while small deviations from the strict separability are quite acceptable. In particular, the possibility of applying the correctly formulated MS-MR-PT-based methods to the calculation of vertical electronic excitation energies should not be affected by these deviations.
III. OVERVIEW OF THE MCQDPT APPROACH

into a zero-order Hamiltonian H0 and a perturbation V:
0 H = HMCQDPT + V .

(1)

The model Fock-like operator F : F = Con s t +
pq

f pq a + aq pq , p

(2a)

0 used to define the zero-order Hamiltonian HMCQDPT is diagonal:

^ F = Con s t +
p

p a + a p . p

(2b)

Because XMCQDPT is closely related to the MCQDPT (Multi-Configuration Quasi-Degenerate Perturbation Theory) approach suggested by Nakano17 as a multi-state generalization of the MR-MP2 theory by Hirao,4 it is helpful to review the basics of MCQDPT. Additional details can be found in the original papers on this theory.17­23 MCQDPT is the multi-state multi-reference Van Vlecktype perturbation theory of partially contracted type that employs isometric (i.e., unitary through any PT order) normalization.60­62 More precisely, a model space is spanned by several CI vectors (partial contraction) that are typically obtained as a result of the state-averaged CASSCF procedure (the P subspace), while the secondary, first-order interacting space is formed by the individual CSFs or determinants and thus is not contracted (the S subspace). It is also helpful to define the O subspace as a subspace of the CASCI vectors that are complementary to P, with PO forming the entire CASCI subspace R, and OS forming the subspace Q, as shown in Fig. 1. The derivation of master equations of MCQDPT can be found in Ref. 17. Traditionally, the total Hamiltonian is split

Here, indices p and q run over set of all MOs, is the spin variable, Const is an arbitrary fixed constant, and p are the or0 ^ bital energies. As stated in Ref. 17, HMCQDPT is identical to F . The zero-order energy of an arbitrary CSF B is given by ^ E 0 = B | F | B = Const + B
i

n i ( B )i ,

(3)

where index i runs over all orbitals occupied in CSF B and n i are the corresponding occupation numbers. The zero-order energy of any CI vector with CI expansion coefficients C B 63 is then (the so-called barycentric formula )
0 0 ^ E = H = | F | = Const + B

C

2 B

E0 . B

(4)

The resolvent operator for S subspace RS , as defined in Ref. 17, is I | ( R S A )| J = I | S ( R S A ) P | J = 0,
1 E 0 -E J
0 I

I | A |J ,

I S, otherwise

JP

(5)

and similarly for the Q subspace. The expressions for the effective Hamiltonian in the lowest orders of PT are
(0) 0 | Hef f | = E ,

(6)

(0 | Hef +1) | = E f

MCSCF

,

(7)

(2) | Hef f | =

1 | V ( R S V ) | + H .c., 2

(8)

where and denote CASCI eigenvectors belonging to model space. However, there is a very important point to be mentioned, 0 namely, the exact form of HMCQDPT is not actually the one given by Eq. (2b), but is rather as follows. 0 The PP and OO blocks of HMCQDPT are diagonal and have the following form: H
0 PP
MCQDPT

=


^ | | F | |,

(9a)

FIG. 1. Illustration on various subspaces and partitioning used throughout the paper.

H

0 OO

MCQDPT

=


|

^ | F |

|,

(9b)

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-4

Alexander A. Granovsky

J. Chem. Phys. 134, 214113 (2011)
0 However, HMCQDPT does not have any well-defined algebraic transformation properties with respect to transformations of basis within the model space. Moreover, it is evident 0 that for any non-trivial model space HMCQDPT is actually a many-particle operator; and hence, the perturbation is also artificially made to be a many-particle operator rather than the two-particle one as would naturally expect. These two points apply equally to MS-CASPT2 and QD-NEVPT2 theories. More precisely, MCQDPT results in perturbation being a many-particle operator in the R subspace, MS-CASPT2 results in perturbation being a many-particle operator in both R and S subspaces, while QD-NEVPT2 results in perturbation being a many-particle operator in the S subspace. It is worth noting that both MS-CASPT2 and QD-NEVPT2 apply internal contraction and a multi-partitioning scheme by Zaitsevskii and Malrieu74 which introduce an additional source of non-invariance into these theories (as was pointed out by anonymous Referee, the redefinition of S subspace as suggested in Ref. 75 could be in principle used in construction of invariant form of MS-CASPT2). As will be shown below, these facts altogether are the source of various difficulties the above mentioned methods encounter in practical applications even when being applied to small molecules like ozone76, 77 or allene.78 In particular, none of these methods have the desired properties discussed in Sec. II. For instance, let us assume one holds fixed the number of CI roots and values of their corresponding weights in computing the averaged one-particle density matrix that is ^ used to define the model Fock operator F . In this case, the effective second-order MCQDPT2 Hamiltonians of progressively decreasing dimension form the trivial sequence with the n â n effective operator being exactly equal to the (n+1) â (n+1) operator with the row and column, corresponding to the (n+1)th extra CI root, being removed. Evidently, there is no any explicit dependence of the effective Hamiltonian on the dimension of model space in this approach, as the n â n effective operator is simply enclosed into (n+1) â (n+1) one.

where summation runs over CI vectors in the model space for PP block (Eq. (9a)) and over CI vectors in the complement of the model space in CASCI space for OO block (Eq. (9b)). The SS block is indeed as stated H
0 SS
MCQDPT

^ = FSS .

(10)

0 All other blocks of HMCQDPT are zero; in particular, PO, OP, SP, and PS blocks are zero. This is evident from the expression for the resolvent used in MCQDPT (Eq. (5)), as well as from the explicit expressions for matrix elements of the second-order effective Hamiltonian (Eqs. (38) and (39) of Ref. 17). Moreover, if we limit ourselves to the second order of PT obeying Eqs. (8) and (5), it is sufficient to define only PP (Eq. (9a)), SS (Eq. (10)), OP, PO, SP, and PS (ze0 ros) blocks of HMCQDPT ; while the OO, SO, and OS blocks can be arbitrary. It is worth noting that essentially the same form of H0 was used in the context of MS-MR-PT in the pioneering work by Spiegelmann and Malrieu.64 Similar forms of the "restriction" of a model operator onto model space (Eq. (9)) are used in the MS-CASPT2 (Eqs. 20­23 of Ref. 24) and QD-NEVPT2 (Eqs. 15 and 16 of Ref. 25) theories as well. Upon selecting the interacting subspace and obtaining a formal expression for PT series, one needs a way to uniquely define H0 to complete the construction of the particular MSMR-PT theory. In the case of MCQDPT this is equivalent to ^ the definition of F , that is, of orbitals and orbital energies. MCQDPT typically employs semi-canonical Fock orbitals, which diagonalize blocks of the effective Fock operator yet leaving the solutions of MCSCF equations invariant. The corresponding eigenvalues are used as the orbital energies. The effective Fock operator itself is obtained using closedshell Hartree-Fock-like expression with the first order density matrix taken as the weighted average of the density matrices for selected subset of CI roots spanning the model space. The interesting modification of this scheme known as IVOMCQDPT265­67 replaces semi-canonical MCSCF orbitals by improved virtual orbitals (IVO) and MCSCF procedure by CASCI in the basis of these orbitals (IVO-CASCI).68, 69 The particular choice of interacting space and zero-order Hamiltonian in MCQDPT theory results in a H0 that is fully diagonal in P, O, and S subspaces. The attractive consequence of this feature is a very simple expression for the resolvent (Eq. (5)) that eliminates any need to solve system(s) of linear equations. Instead, the second-order effective Hamiltonian is expressed as the direct sum of different contributions corresponding to each particular class of diagrams, with energy denominators of a very simple structure (Eqs. (38) and (39) of Ref. 17). This opens a way to a very efficient implementation of MCQDPT2 as was shown in70, 71 and is implemented within the FI R E FL Y (formerly, the PC GAMESS) quantum chemistry package,72 which is partially based on the GAMESS (US) (Ref. 73) source code. In particular, the so-called resolvent fitting approach is computationally very efficient.71 The current implementation allows MCQDPT2 calculations of systems with active spaces up to several millions of CSFs and with the overall number of molecular orbitals up to 2000­3000 to be routinely performed on a standalone single-CPU workstation or desktop computer.

IV. THE XMCQDPT THEORY

The difference between MCQDPT and XMCQDPT (Extended Multi-Configuration Quasi-Degenerate Perturbation Theory) is the choice of the H0 operator. More precisely, XMCQDPT employs H0 that is much closer to the model ^ operator F given by the Eq. (2) above. Let us analyze first the explicit expressions for various blocks of H0 used within XMCQDPT framework. Unless explicitly stated otherwise, ^ throughout this section we assume that the model operator F is diagonal and is defined by the Eq. (2b). As the SS block is formed by the individual CSFs, the SS part of H0 is diagonal and is given by Eq. (3). Similarly, the RS and SR blocks are zero. However, the RR part has a nontrivial structure, namely, for two arbitrary normalized vectors and belonging to the P subspace, the expression for matrix elements of H0 is
0 ^ H = | H 0 | = | F | = B

C

B

+

CB E0 . B



(11)

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-5

XMCQDPT: The new approach to MS-MR-PT

J. Chem. Phys. 134, 214113 (2011)

Hence, unlike MCQDPT, the PP block of H0 is evidently nondiagonal in the basis of CASCI eigenvectors. To uniquely define H0 , we now need to define OO, PO, and OP blocks. There is freedom how this can be done. The simplest form of the theory corresponds to the OO block defined by Eq. (9b), with PO and OP blocks being exactly zero (XMCQDPT-I approach). Alternatively, to be more consistent, we can use Eq. (11) to define the OO block of H0 , once again with fully decoupled P and O spaces (XMCQDPT-II). Evidently, at the second order of PT, XMCQDPT-I and XMCQDPT-II result in the same expression for the second-order effective Hamiltonian (Eqs. (12)­ (15) below). Overall, the XMCQDPT-II approach seems to be the preferred formulation of XMCQDPT. The present work is focused on XMCQDPT2 which is the second order formulation of XMCQDPT-I and II. Alternative variants of theory, based on ideas similar to GVVPT2 theory by Khait, Song, and Hoffmann79 and RMS approach by Staroverov and Davidson,80 will be reported elsewhere. As the PP block of H0 is non-diagonal, the formal working equations of XMCQDPT are more complicated than those of MCQDPT, even at the second order of PT. In the case of XMCQDPT2, the only important non-diagonality is that of the PP block of H0 . It can be easily shown that XMCQDPT2 is defined by the following equations (adopting Van Vleckstyle notations):
(0) (0) Hef f = H pp ,

(12)

(1) Hef f = V pp ,

(13)

1 + H ps v 1 + v 1 Hsp , 2 where v 1 is defined by the system of linear equations:
(2) Hef f = 0 0 Hsp + Hss v 1 = v 1 H pp .

(14)

(15)

is also satisfied taking into account the structure of Eqs. (14)­(15) and the non-diagonal nature of the PP block of H0 in the basis of CI vectors forming the model space. As will be shown in Sec. VII, XMCQDPT2 results in the dependence of Heff on the dimension of a model space that approximates the dependence inherent to exact effective Hamiltonians. The approximate size-extensivity, the rate of convergence of the computed XMCQDPT2 energies with respect to increase of the dimension of model space and the continuity of the potential energy surfaces will be discussed in Secs. V­ VII. Considered formally, Eq. (15) is the large set of independent small (nân, where n is the dimension of model space) linear equations defining amplitudes of various excitations. This looks like a serious deficiency for practical applications of XMCQDPT2 for the modeling of large systems. However, one can completely avoid the necessity to solve these equations explicitly. By taking into account the invariance of the theory, one can diagonalize H0 pp and use the rotated CI vectors (the so-called intermediate basis) throughout the entire calculations rather than the original set of CASCI eigenvectors. The eigenvalues of H0 pp define the zero-order energies of intermediate states, with the resolvent given by the same formal expression (Eq. (5)) as for MCQDPT2 approach. This opens a way to a very efficient implementation of XMCQDPT at second order. In our practical implementation, MCQDPT2 and XMCQDPT2 share virtually the same code with only a few extra lines of code being specific to XMCQDPT2. The computational requirements are identical for both methods. While the described formulation of XMCQDPT2 is computationally efficient, it is not invariant with respect to unitary transformations within subspaces of inactive orbitals (labeled by indices i, j in Eq. (18) below), active orbitals (indices k, l) and external orbitals (indices a, b). Nonetheless, the fully invariant generalization of our approach is straightforward and ¯ is based on the model operator F with the following block structure:

Equations (5)-(7) above can be considered as the special case of Eqs. (12)­(15) corresponding to the specific choice of H0 pp within the MCQDPT approach (Eq. (9a)). Equations (14) and (15) defining H(2) eff obey the invariance property that was formulated in Sec. I. Namely, if Eq. (15) is satisfied for some v 1 , then
0 0 Hsp U + Hss (v 1 U ) = (v 1 U ) U + H pp U

¯ F = Con s t +
ij

f¯ij ai+ a

j

+
kl

+ f¯kl ak al + ab 0 pp

+ f¯ab aa ab .

For instance, the invariant definition of H

(0) and Hef f is

(18)

(16)
(0) Hef f 0 ¯ = H = | H 0 | = | F | .

also holds for any unitary U. Applying U to Eq. (14), one gets
(2) U + Hef f U =

(19)

(17) Thus, any unitary transformation of the basis set in the model space transforms H(2) eff according to the usual law of transformation of an arbitrary linear operator. Since the invariance property is satisfied, the XMCQDPT2 theory is evidently equivalent to the fully uncontracted approach provided that the dimension of the model space is equal to the dimension of the CI space. The "explicit dependence on the model space extension" property

1 (U + H ps )(v 1 U ) + U + v 2

+ 1

( Hsp U ) .

Equation (11) above is the special case of Eq. (19) corresponding to the CASSCF reference states and to the use of semi-canonical MOs. The generalized formulation of XMCQDPT2 has the same invariance properties as the underlying MCSCF procedure of general type; however, it is very inefficient computationally due to necessity in solution of huge system of linear equation(s). Provided that the diag¯ onal ansatz to the active-active block of F is reasonable, one can use non unitary-invariant formulation of XMCQDPT2 as an approximation to its fully invariant generalization.

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-6

Alexander A. Granovsky

J. Chem. Phys. 134, 214113 (2011)

V. EXAMPLE 1: A CONICAL INTERSECTION IN THE ALLENE MOLECULE

As has been mentioned above, MCQDPT, MS-CASPT, and QD-NEVPT2 theories assume the PP block of H0 that is defined using projectors onto individual CI vectors and does not have well-defined transformation properties with respect to the unitary transformation of the model space. This results in the non-invariant perturbation theories and introduces an artificial many-particle nature into the perturbation. In contrast, XMCQDPT applies H0 resulting in the invariant theory with balanced perturbation being approximation to the true two-particle operator on either the P subspace (XMCQDPT-I) or in the entire R subspace (XMCQDPT-II) and the true twoparticle operator in the S space. A discussion on the importance of balanced perturbation in PT treatment can be found in Ref. 13 in the context of NEVPT approach. It is natural to expect that non-invariance should lead to significant erratic behavior of the non-invariant theories in the regions of the molecular geometries where the corresponding CI vectors undergo rapid change or are not uniquely defined at all. More precisely, all non-invariant theories are not defined at the geometries corresponding to conical intersections at the underlying MCSCF level, and behave badly in the vicinities of the avoided crossings. In Example 1, we demonstrate this poor behavior for MCQDPT2 and two variants of MS-CASPT2 (SS-SRCASPT2 and MS-MR-CASPT2) in the vicinity of the Minimum Energy Conical Intersection (MECI) point of 11 A and 21 A states of the distorted allene molecule. The geometry at this point is graphically displayed in Fig. 2. The MECI was located at the State-Averaged Complete Active Space SCF (SA-CASSCF)(4,4) level of theory using the GAMESS (US)-style variation of the Dunning-Hays basis set augmented by a single polarization spherical d-shell on each carbon atom. The Cs symmetry was enforced on the allene system. The (4,4) active space was formed by three A and one A orbitals that results in 12 CSFs in the A subspace. The MECI geometry, orbitals (depicted in Fig. S1) and basis set are provided in the supplementary material.51

FIG. 2. The 11 A ­ 21 A MECI geometry of allene molecule shown from thee different perspectives.

Upon locating the MECI point, two-dimension nonrelaxed surface scans in its vicinity were performed for SACASSCF(4,4) as well as for MCQDPT2, SS-SR-CASPT2 (the originally proposed MS-CASPT224 ), MS-MR-CASPT2 (the variant of MS-CASPT2 procedure81, 82 that is specific to MOLPRO (Ref. 83) package) and XMCQDPT2 approaches. The energies of both crossing states were collected for the effective Hamiltonians of dimensions 2 â 2, 6 â 6, and 12 â 12. Three chemical core orbitals were excluded from the correlation treatment. Convergence was not achieved in 6 â 6 and 12 â 12 MS-CASPT2 calculations, even when applying large values (up to 1.0) of the energy denominators shift. The default ISA shift of 0.02 was used in MCQDPT2 and XMCQDPT2 calculations. Nonetheless, it should be noted that it did not have any serious impact on the results of scans as compared to the control set of calculations with zero ISA shift. The geometry scan variables were defined as follows. The first variable was the value of C-C-C bend angle. The second variable corresponded to the simultaneous and equal change of two C-C-C-H torsion angles for two hydrogen atoms located on the same carbon atom. Upon changing the geometry at one of the terminal C atoms, the Cs symmetry was enforced on the whole molecule. Effectively, the second variable describes pyramidalization of the terminal C atoms. The change of each variable covered the range from -10 degrees to +10 degrees, with CASSCF MECI point geometry being the coordinate origin. The step along each variable was set to 0.25 degrees. The scan grid was formed by 81 points along each variable, 6561 points overall. The MOLPRO package was used for SS-SR-CASPT2 and MS-MR-CASPT2 calculations with its default settings controlling the particular choice of H0 within the CASPT2 approach. The FI R E FL Y package was used for SA-CASSCF, MCQDPT2, and XMCQDPT2 calculations. To avoid any misleading effects caused by the use of various convergence thresholds, cutoffs, by numerical noise, and so on, the extra high accuracy was requested throughout all stages of all calculations. The calculated absolute values of the energy difference of two states were then visualized by computing twodimensional color-mapped surface plots of 64 iso-surfaces each. The final results are depicted in Figs. 3­10. As can be seen, the results of the MCQDPT2 and MS-CASPT2 calculations are very similar. In particular, surfaces obtained in the 2â2 MCQDPT2 (Fig. 4) and the MS-MR-CASPT2 (Fig. 6) calculations are virtually the same while the structure of surface obtained using 2â2 SS-SR-CASPT2 approach is slightly different (Fig. 7). Overall, results are exactly as expected. There is a singularity at the CASSCF MECI point, and a rather large area of the entire PES is clearly an artifact of the results from these approaches lacking any physical meaning (Figs. 4­7). The artificial distortion of PES tends to increase upon the extension of the model space (Fig. 5). In contrast, SA-CASSCF (Fig. 3) and XMCQDPT2 (Figs. 8­10) calculations provide smooth and continuous surfaces (naturally, the SA-CASSCF surface has the cusp at the coordinate origin). Moreover, the results obtained in 6 â 6 XMCQDPT2 calculations (Fig. 9) are virtually the same

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-7

XMCQDPT: The new approach to MS-MR-PT
SA-CASSCF(4,4), NState=2

J. Chem. Phys. 134, 214113 (2011)
MCQDPT2, Heff: 12x12, ISA shift=0.02
0

10 8

10 8 6

0 0.005000 0.01000 0.01500 0.02000 0.02500

0.006875

6 4 2
0.01375

4
0.02063

0 -2

Variable 2

2 0 -2 -4

Variable 2

0.02750

0.03437

-4 -6 -8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10
0.04125

-6
0.04813

0.03000 0.03500 0.04000

-8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

0.05500

Variable 1

Variable 1

FIG. 3. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. SA-CASSCF(4,4).

FIG. 5. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. MCQDPT2 with 12 states.
MS-MR-CASPT2, Heff: 2x2
10 8 6
0.01750 0 0.008750

as those for the complete CASCI space of 12 states corresponding to the fully uncontracted limit of XMCQDPT2 (Fig. 10). It should be noted that the existence of some artifacts on the potential energy surfaces computed within either the MCQDPT2 or MS-CASPT2 approaches has already been mentioned in the literature.76, 77 However, to the best of our knowledge, the present work is the first one that gives the correct theoretical explanation of these peculiarities.

4

Variable 2

2 0 -2 -4 -6

0.02625 0.03500 0.04375 0.05250 0.06125 0.07000

VI. EXAMPLE 2: THE LIF BENCHMARK

-8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

The avoided crossing of neutral and ionic 1 + states of LiF molecule near its dissociation limit is usually considered as one of the standard benchmark tests for any quasidegenerate multi-reference theory. This section summarizes our results for this molecule.

Variable 1

FIG. 6. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. Two-state MS-MR-CASPT2.
SS-SR-CASPT2, Heff: 2x2

MCQDPT2, Heff: 2x2, ISA shift=0.02
10 8 6 4 2 0.01750 0.02625 0.03500 0.04375 0.05250 0.06125 0 0.008750

10 8

0 0.008750

6 4
0.01750 0.02625 0.03500 0.04375

Variable 2

2 0 -2 -4 -6 -8 -10

Variable 2

0 -2 -4 -6 -8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

0.05250 0.06125 0.07000

0.07000

-10

-8

-6

-4

-2

0

2

4

6

8

10

Variable 1

Variable 1

FIG. 4. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. Two-state MCQDPT2.

FIG. 7. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. Two-state SS-SR-CASPT2.

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-8

Alexander A. Granovsky
XMCQDPT2, Heff: 2x2

J. Chem. Phys. 134, 214113 (2011)

10 8

0 0.008750

6 4
0.01750 0.02625 0.03500 0.04375

Variable 2

2 0 -2 -4 -6 -8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

0.05250 0.06125 0.07000

Variable 1

FIG. 8. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. Two-state XMCQDPT2.

Similar to Angeli et al.,25 we adopt the same Li(9s5p)/[4s2p], F(9s6p1d)/[4s3p1d] basis set employed by Bauschlicher and Langhoff84 in their full CI (FCI) study on LiF. However, there are several important differences in our computational methodology as well as in the choice of reference calculations. The zero-order wavefunctions were obtained using the same SA-CASSCF(6,6) procedure with the same choice of active space and with averaging over two 1 + states with equal weights. However, unlike Angeli et al.,25 we do not consider reported FCI results84 as a proper reference for comparison. Indeed, the FCI wavefunction includes excitations that do not belong to the first-order interacting space. Next, because there are no more than eight electrons in the correlation treatment, it is natural to expect that the overall effect due to missed unlinked quadruple excitations in multireference singles and doubles CI (MRSDCI) is small. Therefore, it is more consistent to consider the second and third orders of MS-MR-PT theories based on the CAS reference space as approximations to MRSDCI rather than to the FCI

and hence to use MRSDCI energy curves for comparison purposes. Taking this into account, we did not try to mimic the FCI calculations reported in Ref. 56 as closely as possible, e.g., unlike Angeli et al.,25 we did not freeze three sigma orbitals forming SA-CASSCF(6,6) core space at SA-CASSCF(2,2) level, but rather allowed them to completely relax. Moreover, as the true chemical core is formed by only two of these orbitals, the third one was included both in PT and MRSDCI calculations. This seems to be an especially important point since otherwise there would not be any double-occupied inactive orbitals in the PT treatment, and thus calculations would lack several important classes of excitations. All calculations were based on canonicalized SA-CASSCF(6,6) orbitals. MRSDCI calculations were performed using occupation restricted multiple active spaces (ORMAS) (Refs. 85 and 86) code that is the part of the GAMESS ( US ) package. SS-SR-CASPT2 and MS-MRCASPT2 were performed using MOLPRO, while the FI R E FL Y package was used for MCQDPT2 and XMCQDPT2 calculations. The model space was formed by the two states of interest. The results of our calculations are presented on Figs. 11 and 12. The SS-SR-CASPT2 potential energy curves (PECs) and energy splitting curve (ESC) are evidently in the worst agreement with MRSDCI data, while the MS-MR-CASPT2 PECs and ESC are a significant improvement. The energy splitting curves are very similar for both MS-MR-CASPT2 and MCQDPT2 and are reasonably close to that of MRSDCI. Nonetheless, there are significant artifacts on the SS-SRCASPT2, MS-MR-CASPT2 and especially MCQDPT2 PECs and ESCs in the vicinity of the avoided crossing of underlying CASSCF at internuclear distance of 8.5 Bohr. XMCQDPT2 gives artifacts-free PECs and ESC, which are in the best agreement with the MRSDCI PECs and ESC. The XMCQDPT2 value of energy splitting at the avoided crossing point is also best when compared with the MCQDPT2 and MS-MR-CASPT2 results, although the

XMCQDPT2, Heff: 6x6
10 8 6
0.01750 0 0.008750

XMCQDPT2, Heff: 12x12 (fully uncontracted limit)
10 8
0.008750 0

6 4
0.02625 0.03500 0.04375 0.05250 0.01750 0.02625 0.03500 0.04375 0.05250

4

Variable 2

Variable 2

2 0 -2 -4 -6 -8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

2 0 -2 -4 -6

0.06125 0.07000

-8 -10 -10 -8 -6 -4 -2 0 2 4 6 8 10

0.06125 0.07000

Variable 1

Variable 1

FIG. 9. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. Six-state XMCQDPT2.

FIG. 10. Color mapped iso-surface plot of energy separation between of 21 A and 11 A states of allene molecule as a function of scan variables. XMCQDPT2 with 12 states (fully uncontracted limit).

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-9
-106.94 -106.95 -106.96 -106.97

XMCQDPT: The new approach to MS-MR-PT

J. Chem. Phys. 134, 214113 (2011)

-106.98 -106.99 -107.00 -107.01 -107.02 -107.03 7 8 9 10 11 12 13 14
MRSDCI MRSDCI XMCQDPT2 XMCQDPT2 MCQDPT2 MCQDPT2 MS-MR-CASPT2 MS-MR-CASPT2 SS-SR-CASPT2 SS-SR-CASPT2

ing SA-CASSCF procedure degrades. This can be a serious issue for non-invariant MS-MR-PT approaches defined using projectors, which can result in incorrect description of interstate mixing and PESs. In contrast, situation is much better for invariant theories such as XMCQDPT2. Indeed, the final correlated states depend only on the subspace spanned by zero-order wavefunctions (i.e., on the model space) rather than on the zero-order wavefunctions themselves. Therefore, if the number of physically important CSFs shared by zeroorder SA-CASSCF states is less or equal to the dimension of model space, one can a priori expect better and more balanced description of the interstate mixing and PESs by invariant MS-MR-PT theories as compared with non-invariant formulations.
VII. EXAMPLE 3: 11 A1 AND 21 A1 STATES OF CIS-1,3-BUTADIENE

Energy, Hartree

Internuclear distance, Bohr

FIG. 11. Potential energy curves (PECs) of the LiF molecule near the avoided crossing area computed using different MS-MR-PT approaches. MRSDCI PECs are given as reference. PECs were intentionally left nonshifted to common dissociation limit.

position of the minimum is slightly shifted (by ca. 0.25 Bohr as compared with MRSDCI). The absolute values of off-diagonal elements of half-sum symmetrized effective Hamiltonians at the internuclear separation of 11 Bohr are 0.00341 (SS-SR-CASPT2), 0.00317 (MS-MR-CASPT2), 0.00329 (MCQDPT2) and 0.00240 hartree (XMCQDPT2). Following the referee's suggestion, computed spectroscopic constants of LiF molecule are given in Table SI of supplementary materials.51 In multireference electronic structure calculations for mixed electronic states, the relative contributions of different reference configurations tend to differ substantially between the reference SA-CASSCF wavefunctions and the final correlated wavefunction. For instance, in the example above the position of avoided crossing point between neutral and ionic PECs of 1 + states of LiF molecule differs by ca. 3 Bohr between SA-CASSCF and MRSDCI calculations. Obviously, in the region between SA-CASSCF and correlated crossing points, the quality of zero-order wavefunctions obtained us-

1.5

States energy splitting, eV

MRSDCI XMCQDPT2 MCQDPT2 SS-SR-CASPT2 MS-MR-CASPT2

1.0

0.5

0.0 7 8 9 10 11 12 13 14

Internuclear distance, Bohr

FIG. 12. Interstate separation in the LiF molecule near avoided crossing area computed using different MS-MR-PT approaches. MRSDCI curve is given as reference.

As has been already mentioned, we consider convergence of state energies as a function of the dimension of the model space as a very important requirement for any correctly formulated MR-MS-PT. The approximate size-consistency, especially exact or approximate core-separability, is another desirable property. Here we address these properties in the context of the XMCQDPT2 approach. MCQDPT2 is neither size-consistent nor core-separable. In contrast, different variants of CASPT2 as well as NEVPT2 are exactly core-separable; moreover, NEVPT2 is exactly size consistent and strictly separable. While analysis of their MS counterparts is more complicated, both theories are at least exactly core-separable. In the special case of a one-dimensional model space, XMCQDPT2 is completely equivalent to the MR-MP2 approach.4 Due to the structure of its energy denominators, MR-MP2 is neither exactly size-consistent, nor coreseparable. This can be easily verified by performing test computations on simple model systems like the Be2 molecule at large internuclear distances. However, the deviations from exact size-consistency are typically small. Keeping these properties of MR-MP2 in mind, it is logical to expect that XMCQDPT2 should not be in general exactly size-consistent or core-separable as well. This is indeed the case. However, we mention as an important fact that the fully uncontracted limit of XMCQDPT is an exactly sizeconsistent and separable theory through the fourth order of Van Vleck PT expansion.40 Assuming that the convergence of energies with increase of the dimension of model space is reasonably fast, XMCQDPT2 should be approximately sizeconsistent and core-separable for low-lying states even for the quite modest dimensions of effective Hamiltonians, with sizeinconsistency errors rapidly decreasing upon the extension of the model space. These properties of XMCQDPT2 are investigated below in comparison with MCQDPT2 using the interesting testcase of 1 1 A1 and 21 A1 states of cis-1,3-butadiene. As the 21 A1 state is of strongly multireference character, we consider this system as a good test for any MS-MR-PT. The geometry of cis-1,3-butadiene was optimized at the MP2(frozen core)/cc-pVTZ in C2v symmetry group. The equilibrium geometry results are available as part of the supplementary material.51 Located equilibrium structure was

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-10

Alexander A. Granovsky
-155.46 -155.47 -155.48 -155.49 -155.50 -155.51 -155.52 -155.69 -155.70 -155.71 -155.72 -155.73 -155.74 2 4
MCQDPT2/cc-pV5Z(-h) MCQDPT2'/cc-pV5Z(-h) XMCQDPT2/cc-pV5Z(-h) XMCQDPT2'/cc-pV5Z(-h)

J. Chem. Phys. 134, 214113 (2011)

used throughout all subsequent calculations. The PT2 calculations on energies of the 11 A1 and 21 A1 states of cis-1,3butadiene were performed using standard basis set cc-pVDZ, cc-pVTZ, cc-pVQZ, as well as cc-pV5Z(-h) (h shell was removed because current versions of FI R E FL Y do not support h shells). Four chemical core orbitals were excluded from the PT2 treatment. All PT2 calculations employed recommended ISA energy denominators shift of 0.02. As the 21 A1 state is well known to be of valence type,87, 88 the minimal -valence active space (depicted in Fig. S1)51 was used for SA-CASSCF procedure, and the Rydberg-type orbitals were not added to the basis sets. Thus, only 12 CSFs form the 1 A1 subspace of CASSCF(4,4) Hamiltonian. The two lowest roots of 1 A1 symmetry were state-averaged with equal weights both in SA-CASSCF and in all PT2 calculations regardless of the dimension of the model space. As all four active orbitals have significantly non-trivial occupation numbers, we consider this averaging to be a well suited for the definition of zero order states, orbitals and orbital energies aimed to describe 11 A1 and 21 A1 states within the MR-MSPT approach. Good description of higher-energy states would require averaging over larger number of states; however, this was not the goal of our study. Yet, when entering effective Hamiltonians, higher-energy states allow efficient relaxation of 11 A1 and 21 A1 states within CASCI space and improve their overall description. Indeed, CASSCF usually tends to converge to solutions that are too delocalized in the space of CI coefficients, so additional states in the model space are needed to eliminate this deficiency at XMCQDPT2 level. The results of our calculations are shown in Figs. 13 and 14 as a function of the dimension of effective Hamiltonians. The dependence on the basis set is demonstrated in Figs. S2­S12. Here, XMCQDPT2 and MCQDPT2 acronyms are used for modified XMCQDPT2 and MCQDPT2 calculations that apply the MP2-like expression for double excitations from double occupied inactive orbitals to external MOs. This eliminates the dominant contribution violating exact coreseparability and makes both theories approximately coreconsistent. However, as can be easily seen, this modification is not really necessary for XMCQDPT2 as it results in virtually the same state energies as the exact XMCQDPT2, surprisingly even for very small dimensions of the model spaces. This allows us to conclude that XMCQDPT2 is almost exactly core-consistent, at least for this specific model system, and we expect this property to hold for other systems as well. On the other hand, this modification clearly has a very significant impact on the computed MCQDPT2 energies and seriously affects the 11 A1 ­ 21 A1 excitation energy as well. In the case of XMCQDPT2, we can see in Figs. 13 and 14 and S2­S10, the energies of both states are smooth and almost monotonic functions of the dimension of model space n, saturating at n = 7. This behavior is very encouraging and is common for all basis sets we used. Moreover, energies of both states become even smoother for the larger cc-pVQZ and ccpV5Z(-h) basis sets. Analysis of the corresponding eigenvectors of effective Hamiltonians shows that the seventh CASCI state is the last one that significantly contributes to the projection of the perturbed 11 A1 and 21 A1 states onto the model

Energy, Hartree

6

8

10

12

Dimension of a model space

FIG. 13. Computed PT2/cc-pV5Z(-h) energies of 11 A1 and 21 A1 states of cis-1,3-butadiene as the functions of the dimension of a model space. Lower panel: 11 A1 state. Upper panel: 21 A1 state.

space. At the same time, there are large changes in the state energies in the low model space dimension regions of both MCQDPT2 and MCQDPT2 energies, and there is no saturation but rather a constant steady decrease of energies, up to n = 12. The latter can be easily explained taking into account the variational principle and the fact that with the stateaveraging procedure we used, MCQDPT2 and MCQDPT2 lead to the (n+1) â (n+1) effective Hamiltonian being simply the augmented version of the n â n one. Given these results, we conclude that effective Hamiltonians of increasing dimension computed within the XMCQDPT2 framework tend to approximate the dependence on the dimension of a model space inherent to exact effective Hamiltonians. Examination of the effective Hamiltonians computed within the MCQDPT2 approach reveals multiple offdiagonal elements of unrealistically large magnitude. Use of MCQDPT2 instead of MCQDPT2 decreases their magnitude, sometimes by a factor of 1.5. However, on average they are still severely overestimated as compared to the
6.7 6.6

Vertical excitation energy, eV

6.5 6.4 6.3 6.2 6.1 6.0 5.9 5.8 5.7 2 4 6

MCQDPT2/cc-pV5Z(-h) MCQDPT2'/cc-pV5Z(-h) XMCQDPT2/cc-pV5Z(-h) XMCQDPT2'/cc-pV5Z(-h)

8

10

12

Dimension of a model space

FIG. 14. Computed PT2/cc-pV5Z(-h) 11 A1 to 21 A1 vertical excitation energies of cis-1,3-butadiene as the functions of the dimension of a model space.

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-11

XMCQDPT: The new approach to MS-MR-PT

J. Chem. Phys. 134, 214113 (2011)

XMCQDPT2 results. These overestimated off-diagonal elements cause observed changes and oscillations. The overestimation of the off-diagonal elements within the MCQDPT2 and MS-CASPT2 approaches has already been documented in the literature by several authors.89­97 We postpone the discussion on the nature of this overestimation to Sec. VIII. However, we note in passing that the absolute energies of individual states (or, more generally, PES segments) computed within the MCQDPT2 approach seem to have little, if any, physical sense, at least for large effective Hamiltonians. Let us now examine the behavior of the 11 A1 ­ 21 A1 excitation energies. Again, XMCQDPT2 results in excitation energies being almost smooth and monotonic functions of the model space dimension showing only minor variations upon reaching n = 7. The excitation energy curves become smoother as the quality of the basis set improves. The excitation energies obtained using cc-pVDZ basis set are significantly overestimated being compared with our results for other basis sets; so we do not recommend use of this basis set in any high-level XMCQDPT2 calculations. On the other hand, results for cc-pVTZ, cc-pVQZ, and cc-pV5Z(-h) basis sets are all similar. In particular cc-pVQZ and cc-pV5Z(-h) results are virtually the same while cc-pVTZ somewhat overestimates excitation energy. As expected, MCQDPT2 gives oscillating estimates of excitation energy up to n = 7. At larger n there are only small to modest oscillations near some average value. Nevertheless, the curves are not as smooth compared with those of XMCQDPT2 curves and do not seem to exhibit a tendency to become smoother as the quality of basis set improves. Once again, there is a considerable difference between the MCQDPT2 and MCQDPT2 curves. MCQDPT2 tends to attenuate oscillations, although is not able to remove them completely in the initial region. It also shows much better agreement with converged and basis-set saturated XMCQDPT2 results. The converged and nearly basis-set saturated XMCQDPT2 estimate of the 11 A1 ­ 21 A1 excitation energy in cis-1,3-butadiene is ca. 6.12­6.15 eV (Fig. S12), while the unperturbed SA-CASSCF/cc-pV5Z(-h) value is 6.535 eV. Since the exact experimental excitation energy is not known,87, 88, 98­100 the XMCQDPT2 estimate can be compared with 6.04 eV CASPT2-based result by Serrano-AndrÈs et al.88 and 6.13 eV CI4 + Q result by Cave and Davidson.87
VIII. DISCUSSION AND CONCLUDING REMARKS

action increases, MCQDPT2 performance starts to degrade. For instance, one can apriori expect large zero-order interactions in situations when two or more CI roots share a common set of leading CSFs, resulting in questionable applicability of MCQDPT2 as well as of all other non-invariant theories formulated using projectors onto the model space. The neglect of off-diagonal terms introduces inaccuracy into both diagonal and off-diagonal elements of effective Hamiltonians. One recalls that in the basis of non-rotated CI vectors, the working equations of XMCQDPT2 (Eq. (15)) are the set of linear equations of the form: H with solution v = H ¯
0 pp 0 pp

-D



¯ v = b , ¯
-1

(20)

-D



¯ b ,

(21)

where D is a diagonal matrix and index runs over CSFs belonging to S subspace. For an arbitrary square matrix A, let us denote the matrix obtained from A by zeroing all its off-diagonal elements as [A]. Similarly, {A} denotes the matrix A with all diagonal elements zeroed. The equations for MCQDPT2 are then v ¯
,MCQDPT

=

H

0 pp

-D



-1

¯ b .

(22)

0 Considering off-diagonal elements of matrix H pp - D as perturbation and taking into account standard arguments of perturbation theory, one gets the following first-order approx0 imation to the exact inverse of H pp - D : 0 H pp - D -1

=

H

0 pp

-D H
0 pp



1+ H
-1 -1

0 pp

-D
0 pp



-1 -1 -1

H H H

0 pp

-1

= 1+ H
0 pp

-D

H
0 pp

0 pp

-D H



-1

-D

-

H

-D



0 pp

0 pp

-D



. (23)

-1

Using Eq. (23) instead of the exact inverse, one gets v ¯ â H
0 pp

-D
0 pp



-1

¯ b -
-1

H

0 pp

-D



-1

H H

0 pp 0 pp

H

-D H



¯ b = v ¯
-1

,MCQDPT

-

-D



-1

âH

0 pp

0 pp

-D

¯ b .

(24)

As has been mentioned previously, MCQDPT2 and MS-CASPT2 tend to overestimate off-diagonal elements of effective Hamiltonians.89­97 This is true, although to a lesser extent, for MCQDPT2 as well. This can be explained as follows. One can consider MCQDPT2 as the approximation to XMCQDPT2 which completely neglects the off-diagonal elements of H0 in the CASCI space. The quality of this approximation depends on the magnitude of the neglected off-diagonal elements. If the zero-order interaction of states forming the model space is small, these two theories should be in the close agreement. However, as soon as zero-order inter-

The second term of Eq. (24) is therefore the first-order correction to the MCQDPT2 amplitudes within XMCQDPT2 approach. The structure of this correction taken together with Eqs. (14) and (15) allows one to conclude that large offdiagonal elements h 0 of H0 pp result in large systematic errors in H(2) eff, elements of H(2) eff , MCQDPT . Similar arguments can be applied to other MS-MR-PT theories which neglect offdiagonal terms, i.e., to MS-CASPT2. In the case of MCQDPT2 , the leading MP2-like contribution to the correlation energy is treated differently so it results only in the equal shift of all diagonal elements of H(2) eff . Hence, this reduces the accumulated errors in offdiagonal elements of H(2) eff, MCQDPT . However, this does not eliminate errors entirely. For example, for the all-trans protonated Shiff-base form of retinal molecule (540 MOs: 25 frozen core orbitals, 63 doubly occupied inactive, 12 active

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-12

Alexander A. Granovsky

J. Chem. Phys. 134, 214113 (2011)

orbitals, 440 external orbitals), the two-state MCQDPT2/SACASSCF(12,12) treatment at PBE0/cc-pVDZ geometry51 results in the off-diagonal element of the effective Hamiltonian equal to 0.0833 hartree, with the difference of diagonal values of the effective Hamiltonian equal to 0.0763 hartree. The computed S0-S1 vertical excitation energy is 0.1833 hartree (249 nm) and is unphysical. MCQDPT2 gives an off-diagonal value of 0.0294 hartree that is still severely overestimated, with a diagonal difference of 0.0744 hartree, resulting in a S0S1 vertical excitation energy of 0.0948 hartree (481 nm). XMCQDPT2 gives an off-diagonal value of 0.0027 hartree, with a diagonal difference of 0.0767 hartree, and a S0-S1 vertical excitation energy of 0.07687 hartree (593 nm). Note that for MCQDPT2 and MCQDPT2 (but not for XMCQDPT2), the diagonal values of Heff are approximations to the state-specific MR-MP2 energies. The unperturbed SA-CASSCF value of S0-S1 vertical excitation energy is 0.09935 hartree (458 nm). Experimentally, the S0-S1 band in the gas phase is broad, ranging from 450 nm to 650 nm with an essentially flat top extending from 530 nm to 610 nm.101 In the case of XMCQDPT2, there is reasonable agreement with the experimental data that can be further improved using larger basis sets, extended model spaces, MP2-optimized geometry, and averaging over thermally available vibrational states.102 Remarkably, XMCQDPT2 computations on retinal molecule, including all required integral transformation and CI steps, took as fast as ca. 3250 seconds of wall clock time running on fourcore single CPU Intel Core i7 2600K-based desktop system. The systematic overestimation of off-diagonals within MCQDPT2 has been documented previously (2005­ 2006).89­91 The development of XMCQDPT was completed in 2008,96 and initially led to development of the so-called augmented effective Hamiltonian approach to MCQDPT2 (aug-MCQDPT2, or, more precisely, aug-MCQDPT2 ) by Bochenkova.89­95 This procedure can be considered as the attempt to apply a variation of the stabilization graph technique to the MCQDPT2 effective Hamiltonians of sequentially enlarged dimensions aimed to minimize effects caused by overestimated off-diagonal elements of the effective Hamiltonian. Given our results, and the accompanied analysis, we believe that at least some, if not most, of the previously reported problematic MCQDPT2 and MS-CASPT2 cases that were attributed to intruders or interaction with discretized continuum, are probably related to the intrinsic property of these methods to overestimate off-diagonal coupling. As to real intruders, our preliminary experience shows that XMCQDPT2 is somewhat less prone to intruders than MCQDPT2 and MS-CASPT2. The sensitivity of XMCQDPT2 to intruders will be examined in details in another study. Over the last decades, minimum energy conical intersections (MECIs) and associated dynamics have been attracting constantly growing attention of the Scientific Community. At present, MECIs are typically located using the SA-CASSCF procedure. It appears that for large molecular systems, both MCQDPT2 and MS-CASPT2 are not suitable for locating conical intersections103, 104 because of the uncertainty in the off-diagonal elements of corresponding effective Hamiltonians. However, the possibility of locating MECI points using

high-level correlated methods appears attractive. With this in mind, two further developments seem to be very important. The first is the formulation and development of the efficient analytical gradient procedure for XMCQDPT2. The second is the "XMCQDPT-like" reformulation of the existing MS-CASPT2 and QD-NEVPT2 approaches. To summarize, we formulated and implemented the new approach to MS-MR-PT, called XMCQDPT, which is a multistate generalization of MR-MP2 based on the non-diagonal, invariant choice of zero-order Hamiltonian inside the model space. This approach possesses several unique properties as discussed above. It is an interesting, useful, and attractive tool for high level modeling of the electronic properties of small to large molecular systems. After its development and implementation in 2008, this method was made routinely available for FI R E FL Y users, and has already been successfully applied by some of them to the modeling of challenging systems.105­111
ACKNOWLEDGMENTS

Author is very grateful to Dr. Anastasia V. Bochenkova, Dr. Vladimir I. Pupyshev, and Dr. Andrei V. Zaitsevskii for numerous stimulating discussions and to Dr. James W. Kress and Dr. Bruce S. Greer for reviewing draft of manuscript. Multiple helpful comments made by anonymous Referees are appreciated. This work was partially supported by the RFBR Grant (project No 11­03-01214-a).
1 2 3 4 5 6 7 8 9 10 11

12 13 14 15 16

17 18 19 20 21 22

C. MÜller and M. S. Plesset, Phys. Rev. 46, 618 (1934). J. A. Pople, J. S. Binkley, and R. Seeger, Int. J. Quantum Chem., Quantum Chem. Symp. S10, 1 (1976). M. J. Frisch, M. Head-Gordon, and J. A. Pople, Chem. Phys. Lett. 166, 275 (1990). K. Hirao, Chem. Phys. Lett. 190, 374 (1992). K. Hirao, Chem. Phys. Lett. 196, 397 (1992). K. Hirao, Int. J. Quantum Chem. S 26, 517 (1992). K. Hirao, Chem. Phys. Lett. 201, 59 (1993). K. Andersson, P-å. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, J. Phys. Chem. 94, 5483 (1990). K. Andersson, P.-å. Malmqvist, and B. O. Roos, J. Chem. Phys. 96, 1218 (1992). K. Andersson and B. O. Roos, Int. J. Quantum Chem. 45, 591 (1993). B. O. Roos, K. Andersson, M. P. FÝlscher, P.-å. Malmqvist, L. SerranoAndrÈs, K. Pierloot, and M. MerchÀn, in Advances in Chemical Physics: New Methods in Computational Quantum Mechanics, Edited by I. Prigogine and S. A. Rice (Wiley, New York, 1996), vol. XCIII, p. 219. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, J. Chem. Phys. 114, 10252 (2001). C. Angeli, R. Cimiraglia, and J.-P. Malrieu, J. Chem. Phys. 117, 9138 (2002). C. Angeli, R. Cimiraglia, and J.-P. Malrieu, Chem. Phys. Lett. 350, 297 (2001). C. Angeli, S. Borini, and R. Cimiraglia, Theor. Chem. Acc. 111, 352 (2004). R. K. Chaudhuri, K. F. Freed, G. Hose, P. Piecuch, K. Kowalski, M. Wloch, S. Chattopadhyay, D. Mukherjee, Z. Rolik, à. Szabados, G. TÑth, and P. R. SurjÀn, J. Chem. Phys. 122, 134105 (2005). H. Nakano, J. Chem. Phys. 99, 7983 (1993). H. Nakano, Chem. Phys. Lett. 207, 372 (1993). H. Nakano, K. Nakayama, K. Hirao, and M. Dupuis, J. Chem. Phys. 106, 4912 (1997). T. Hashimoto, H. Nakano, and K. Hirao, J. Mol. Struct.: (THEOCHEM) 451, 25 (1998). H. Nakano, R. Uchiyama, and K. Hirao, J. Comput. Chem. 23, 1166 (2002). M. Miyajima, Y. Watanabe, and H. Nakano, J. Chem. Phys. 124, 044101 (2006).

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-13
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

XMCQDPT: The new approach to MS-MR-PT

J. Chem. Phys. 134, 214113 (2011) pdf) for presentation on resolvent-fitting approach to MCQDPT2 and XMCQDPT2 theories. Accessed October 20, 2010. A. A. Granovsky, FI R E FL Y Quantum Chemistry Package, version 7.1.G, http://classic.chem.msu.su/gran/firefly/index.html. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem. 14, 1347 (1993). A. Zaitsevskii and J.-P. Malrieu, Chem. Phys. Lett. 223, 597 (1995). P. J. Knowles and H.-J. Werner, Theor. Chim. Acta 84, 95 (1992). L. De Vico, L. Pegado, J. Heimdal, P. SÆderhjelm, and B. O. Roos, Chem. Phys. Lett. 461, 136 (2008). See I. Ryabinkin, http://classic.chem.msu.su/gran/gamess/mcqdpt-vsxmcqdpt-in-O3.pdf (mirrored at http://phys148-1.chem.msu.ru/gran/ gamess/mcqdpt-vs-xmcqdpt-in-O3.pdf) for presentation on performance of MCQDPT2 and XMCQDPT2 in modeling of excited states of ozone molecule. Accessed October 20, 2010. See A. A. Granovsky, http://classic.chem.msu.su/gran/gamess/conic-qdpt. pdf (mirrored at http://phys148-1.chem.msu.ru/gran/gamess/conic-qdpt. pdf) for presentation discussing performance of various MS-MR approaches in modeling of conical intersection in allene molecule. Accessed October 20, 2010. Y. G. Khait, J. Song, and M. R. Hoffmann, J. Chem. Phys. 117, 4133 (2002). V. N. Staroverov and E. R. Davidson, Chem. Phys. Lett. 296 435 (1998). H.-J. Werner, Mol. Phys. 89, 645 (1996). P. Celani and H.-J. Werner, J. Chem. Phys. 112, 5546 (2000). MOLPRO , version 2006.1, a package of ab initio programs, H.-J. Werner, P. J. Knowles, R. Lindh, et al., see http://www.molpro.net. C. W. Bauschlicher and S. R. Langhoff, J. Chem. Phys. 89, 4246 (1988). J. Ivanic, J. Chem. Phys. 119, 9364 (2003). J. Ivanic, J. Chem. Phys. 119, 9377 (2003). R. Cave and E. Davidson, J. Phys. Chem. 91, 4481 (1987). L. Serrano-AndrÈs, B. O. Roos, and M. MerchÀn, Theor. Chim. Acta 87, 387 (1994). A. V. Nemukhin, A. V. Bochenkova, K. B. Bravaya, and A. A. Granovsky, Proc. SPIE 6449, 64490N (2007). K. Bravaya, A. Bochenkova, A. Granovsky, and A. Nemukhin, J. Am. Chem. Soc. 129(43), 13035 (2007). L. Kessel, I. B. Nielsen, A. V. Bochenkova, K. B. Bravaya, and L. H. Andersen, J. Phys. Chem. A 111(42), 10537 (2007). K. B. Bravaya, A. V. Bochenkova, A. A. Granovsky, A. P. Savitsky, and A. V. Nemukhin, J. Phys. Chem. A 112(37), 8804 (2008). K. Bravaya, A. Bochenkova, A. Granovsky, A. Nemukhin, Russ. J. Phys. Chem. B 2, 671 (2008). T. Rocha-Rinza, O. Christiansen, J. Rajput, A. Gopalan, D. B. Rahbek, L. H. Andersen, A. V. Bochenkova, A. A. Granovsky, K. B. Bravaya, A. V. Nemukhin, K. L. Christiansen, and M. B. Nielsen, J. Phys. Chem. A 113(34), 9442 (2009). M. G. Khrenova, A. V. Bochenkova, and A. V. Nemukhin, Proteins: Structure, Function, and Bioinformatics 78, 614 (2010). See A. A. Granovsky, http://classic.chem.msu.su/gran/gamess/xmcqdpt. pdf (mirrored at http://phys148-1.chem.msu.ru/gran/gamess/xmcqdpt. pdf) for presentation on the initial implementation of XMCQDPT2 in FI R E FL Y package. Accessed October 20, 2010. L. Serrano-AndrÈs, M. MerchÀn, and R. Lindh, J. Chem. Phys. 122, 104107 (2005). O. A Mosher, W. M. Flicker, and A. Kuppermann, J. Chem. Phys. 59, 6502 (1973). J. P. Doering and R. McDiarmid, J. Chem. Phys. 75, 2477, (1981). R. R. Chadwick, D. P. Gerrity, and B. S. Hudson, Chem. Phys. Lett. 115, 24 (1985). J. Rajput, D. B. Rahbek, L. H. Andersen, A. Hirshfeld, M. Sheves, P. AltoÕ, G. Orlandi, and M. Garavelli, Angew. Chem. 122, 1834 (2010). A. V. Bochenkova and A. A. Granovsky, "Anomalous gas-phase photoabsorption profile of Retinal: a high level ab initio study," Phys. Rev. Lett. (to be published). B. G. Levine, J. D. Coe, and Todd J. MartÌnez, J. Phys. Chem. B 112, 405 (2008). Toshifumi Mori and Shigeki Kato, Chem. Phys. Lett. 476, 97 (2009). I. V. Polyakov, B. L. Grigorenko, E. M. Epifanovsky, A. I. Krylov, and A. V. Nemukhin, J. Chem. Theory Comput. 6, 2377 (2010).

52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70

71

R. Ebisuzaki, Y. Watanabe, and H. Nakano, Chem. Phys. Lett. 442, 164 (2007). J. Finley, P.-å. Malmqvist, B. O. Roos, and L. Serrano-AndrÈs, Chem. Phys. Lett. 288, 299 (1998). C. Angeli, S. Borini, M. Cestari, and R. Cimiraglia, J. Chem. Phys. 121, 4043 (2004) C. Bloch, Nucl. Phys. 6, 329 (1958). B. H. Brandow, Rev. Mod. Phys. 39, 771 (1967). B. H. Brandow, Int. J. Quant. Chem. 15, 207 (1979); I. Lindgren, J. Phys. B 7, 2441 (1974). I. Lindgren, Int. J. Quant. Chem. 14(S12), 33 (1978). G. Hose and U. Kaldor, J. Phys. B 12, 3827 (1981). J.-P. Daudey and J.-P. Malrieu, in Current Aspects of Quantum Chemistry 1981, edited by R. CarbÑ (Elsevier, New York, 1982), p. 35. Ph. Durand and J.-P. Malrieu, Adv. Chem. Phys. 67, 321 (1987). K. F. Freed and M. G. Sheppard, J. Chem. Phys. 75, 4507 (1981). K. F. Freed and M. G. Sheppard, J. Chem. Phys. 75, 4525 (1981). V. Hurtubise and K. F. Freed, Adv. Chem. Phys. 83, 465 (1993). J.-P. Malrieu, J.-L. Heully, and A. Zaitsevskii, Theor. Chim. Acta 90, 167 (1995). U. S. Mahapatra, B. Datta, and D. Mukherjee, Chem Phys Lett 299, 42 (1999). U. S. Mahapatra, B. Datta, and D. Mukherjee, J. Phys. Chem. A 103, 1822 (1999). P. Ghosh, S. Chattopadhyay, D. Jana, and D. Mukherjee, Int. J. Mol. Sci. 3, 733 (2002) U. S. Mahapatra, S. Chattopadhyay, and R. K. Chaudhuri , J. Chem. Phys. 129, 024108 (2008). H.-J. Werner and P. J. Knowles, J. Chem. Phys. 89, 5803 (1988). P. E. M. Siegbahn, Chem. Phys. 25, 197 (1977). H. Zhang, J.-P. Malrieu, P. Reinhardt, and J. Ma, J. Chem. Phys. 132, 034108 (2010). K. G. Dyall, J. Chem. Phys. 102, 4909 (1995). T. H. Schucan and H. WeidenmÝller, Ann. Phys. (N.Y.) 73, 108 (1972). T. H. Schucan and H. WeidenmÝller, Ann. Phys. (N.Y.), 76, 483 (1973). C. Camacho, R. Cimiraglia, and H. A. Witek, Phys. Chem. Chem. Phys. 12, 5058 (2010). C. Camacho, H. A. Witek, and R. Cimiraglia, J. Chem. Phys. 132, 244306 (2010). H. A. Witek, Y.-K. Choe, J. P. Finley, and K. Hirao, J. Comput. Chem. 23, 957 (2002). See supplementary material at http://dx.doi.org/10.1063/1.3596699 for Appendix SI, Table SI, Figures S1-S12 and data on basis sets, equilibrium geometries, and MOs. M. Lewin, Arch. Ration. Mech. Anal. 171, 83 (2004). è. CancÕs, H. Galicher, and M. Lewin, J. Comp. Phys. 212, 73 (2006). M. Lewin, J. Math. Chem. 44, 967 (2008). J. A. Pople, J. S. Binkley, and R. Seeger, Int. J. Quantum Chem., Symp. 10, 1 (1976). K. A. Brueckner, Phys. Rev. 100, 36 (1955). A. Hugenholtz, Physica 23, 481 (1957). R. J. Bartlett, Ann. Rev. Phys. Chem. 32, 359 (1981). M. Hanrath, Chem. Phys. 356, 31 (2009). B. Kirtman, J. Chem. Phys. 49, 3890 (1968). B. Kirtman, J. Chem. Phys. 75, 798 (1981). I. Shavitt and L. R. Redmon, J. Chem. Phys. 73, 5711 (1980). B. Huron, J. P. Malrieu, and P. Rancurel, J. Chem. Phys. 58, 5745 (1973). F. Spiegelmann and J.-P. Malrieu, J. Phys. B 17 1235 (1984). R. K. Chaudhuri, K. F. Freed, S. Chattopadhyay, and U. S. Mahapatra, J. Chem. Phys. 128, 144304 (2008). S. Chattopadhyay, R. K. Chaudhuri, and U. S. Mahapatra, J. Chem. Phys. 129, 244108 (2008). R. K. Chaudhuri and K. F. Freed, J. Chem. Phys. 129, 054308 (2008). D. M. Potts, C. M. Taylor, R. K. Chaudhuri, and K. F. Freed, J. Chem. Phys. 114, 2592 (2001). S. Chattopadhyay, R. K. Chaudhuri, and K. F. Freed, Phys. Chem. Chem. Phys. 13, 7514 (2011). See A. A. Granovsky, http://classic.chem.msu.su/gran/gamess/qdpt2.pdf (mirrored at http://phys148-1.chem.msu.ru/gran/gamess/qdpt2.pdf) for presentation on efficient strategy for summation of PT series in MCQDPT2 and XMCQDPT2 implementation. Accessed October 20, 2010. SeeA.A.Granovsky, http://classic.chem.msu.su/gran/gamess/table_qdpt2. pdf (mirrored at http://phys148-1.chem.msu.ru/gran/gamess/table_qdpt2.

72 73

74 75 76 77

78

79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94

95 96

97 98 99 100 101 102

103 104 105

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp


214113-14
106 107

Alexander A. Granovsky
109

J. Chem. Phys. 134, 214113 (2011) A. Y. Freidzon, A. A. Bagatur'yants, E. N. Ushakov, S. P. Gromov, and M. V. Alfimov, Int. J. Quantum Chem. 111, 2649 (2011). A. Udvarhelyi and T. Domratcheva, Photochem. Photobiol. 87, 554 (2011). A. Ya. Freidzon, A. V. Scherbinin, A. A. Bagaturyants, and M. V. Alfimov, J. Phys. Chem. A 115, 4565 (2011).

108

E. Epifanovsky, I. Polyakov, B. Grigorenko, A. Nemukhin, and A. I. Krylov, J. Chem. Phys. 132, 115104 (2010). C. Bruno, F. Paolucci, M. Marcaccio, R. Benassi, C. Fontanesi, A. Mucci, F. Parenti, L. Preti, L. Schenetti, and D. Vanossi, J. Phys. Chem. B 114(26), 8585 (2010). P. M. Kozlowski, T. Kamachi, M. Kumar, T. Nakayama, and K. Yoshizawa, J. Phys. Chem. B 114, 5928 (2010).

110 111

Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp