Modification for spin-adapted version of configuration interaction singles with perturbative doubles

Chemical Physics Letters 443 (2007) 389–397 www.elsevier.com/locate/cplett

Modification for spin-adapted version of configuration interaction singles with perturbative doubles Yuji Mochizuki

a,b,c,*

, Kiyoshi Tanaka

c,d

a

Department of Chemistry, Faculty of Science, Rikkyo University, 3-34-1 Nishi-ikebukuro, Toshima-ku, Tokyo 171-8501, Japan b CREST Project, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan c Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8904, Japan d Advancesoft, Center for Collaborative Research, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8904, Japan Received 13 April 2007; in final form 12 June 2007 Available online 17 June 2007

Abstract We propose a couple of modifications for the spin-adapted version of configuration interaction singles with perturbative doubles, originated by Head-Gordon et al. as CIS(D). A partial renormalization technique is first introduced for the T2 doubles operator of second-order Møller–Plesset perturbation (MP2) theory. Second, an extra correction from the T1 singles operator is taken into account. The combined modifications can be abbreviated as PR-CIS(Ds) (partially renormalized variant of configuration interaction singles with perturbative doubles and extra singles). A series of calculations with a benchmark set of molecules are performed for both singlet and triplet excited states. Through these tests, it is observed that favorable improvements are obtainable by the proposed modifications for low-lying valence states. Ó 2007 Elsevier B.V. All rights reserved.

1. Introduction Calculations of the configuration interaction singles (CIS) [1] provide semi-quantitative estimates for the lowlying states, which are predominantly characterized by single excitations from the Hartree–Fock (HF) ground state, at the low computational cost of iterative N4 (where N is the number of basis functions as a symbolic index of molecular size) like the HF calculation. For the most of closedshell molecules, a typical overestimation in CIS excitation energies has been known to be 1–2 eV relative to experimental values. Foresman et al. [1] first proposed a correction recipe to introduce the correlation effects for a certain CIS state in analogy with the second-order Møller–Plesset (MP2) perturbation theory for the ground state.

* Corresponding author. Address: Department of Chemistry, Faculty of Science, Rikkyo University, 3-34-1 Nishi-ikebukuro, Toshima-ku, Tokyo 171-8501, Japan. Fax: +81 3 3985 2407. E-mail address: [email protected] (Y. Mochizuki).

0009-2614/$ - see front matter Ó 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2007.06.059

This CIS-MP2 method consists of doubles and triples contributions, where the cost of computations is N5 for the former and N6 for the latter. Head-Gordon et al. [2] pointed out that the triples contribution in CIS-MP2 [1] violates the size-consistency. They thus devised an alternative recipe by factorizing the triples operator into the product of the MP2 doubles operator and the CIS operator. Through the factorization, the resulting CIS(D) (meaning the CIS with perturbative doubles correction) method is made size-consistent, and the scaling of computational cost is favorably reduced to be N5 as just in the case of MP2 calculation. Ref. [2] actually showed that the numerical results by CIS(D) outperform those by CIS-MP2 for singlet states of benchmark molecules. Head-Gordon et al. [2,3] discussed that the CIS(D) can be considered as a non-iterative second-order approximation to the excited states methods based on the coupled-cluster singles and doubles (CCSD) [4–6]. In Ref. [3], they proposed a partial third-order extension as CIS(D3) having the computational scaling of N6. A complete third-order extension and a partial fourth-order extension were achieved by Hirata [7] as CIS(3) and

390

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

CIS(4)p, respectively. For the issue of near-degeneracy of states, the iterative variants of CIS(D) were developed by Head-Gordon et al. as I-CIS(D) [3,8] and also as more extensive CIS(Dn) (n = 0, 1, 1) [9], where the matrix space of the correlation-dressed Hamiltonian is maintained to be singles-only in the cases of n = 0, 1: this feature differs from the second-order approximate CCSD called as CC2 [10]. Grimme and Izgorodina [11] reported rather an empirical modification in which the MP2 amplitudes are altered by the spin-component scaling (SCS) with respect to ‘parallel’ and ‘anti-parallel’ terms [12]. They applied this SCSCIS(D) method to a variety of molecules for both singlet and triplet excited states, but the observed trend of improvements was not uniform unfortunately [11]. Meissner [13] highlighted a potential risk possessed in the CIS(D) [2] and its higher-order extensions [3,7] in the case of inter-unit charge transfer states. For this issue, he claimed the need of ‘perturb-then-diagonalize’ type corrections. Recently, Fan and Hirata [14] proposed some new versions of doubles correction named CIS(2) or DCIS(2), and the performance of these recipes has been in testing. From a view point of practical applicability to large molecular systems, the non-iterative N5 scaling of original CIS(D) calculations [2] has still been attractive, whereas more extensive approaches [3,7–9,13,14] have been proposed. In the methodology of fragment molecular orbital (FMO) [15], we developed a parallelized integral-direct CIS(D) module and equipped it into a local version of ABINIT-MP program [16] by adopting the multilayer framework (MLFMO-CIS(D) [17]). It is notable that the implementation of this CIS(D) module was based not on the spin-orbital formulation of Ref. [2] but on our original formulation with the spin-adaptation for a separate treatment of spin singlet and triplet states [17]. As Ref. [18], we reported a large-scale MLFMO-CIS(D)/6-31G* study of red fluorescent protein, DsRed. The computational chromophore system of DsRed involved the space consisting of 106 occupied MOs and 537 virtual MOs (where N is 684) for the largest case, implying that total number of spinadapted doubles incorporated implicitly was as many as 1.6 billion. The excitation and emission energies were then estimated to be 2.28 eV and 2.21 eV, respectively, being in reasonable agreement with the corresponding experimental values of 2.22 eV and 2.13 eV [18]. Very recently [19], the MLFMO-CIS(D) method [17] has also been used in a simulation study for hydrated formaldehyde in the combination of FMO-based molecular dynamics (MD) to sample the fluctuating configurations of droplet. Although the CIS(D) method [2] could be useful in calculating the excitation energies of large systems efficiently [17–19], care should be taken for its inherent limitation as an approximation to the CCSD excited states methods [3–6]. If some modification to improve the reliability or performance of CIS(D) corrections with keeping the computational feature of N5 is available, it could become a valuable tool to be used concomitantly.

In this Letter, we report a couple of modification recipes for the CIS(D) method [2], based on our spin-adapted formulation [17]. The first way employs the partially renormalized (PR) MP2 approach which was proposed by Dykstra and Davidson [20] by introducing the topological factor devised for the coupled pair functional (CPF) [21]. This spin-adapted PR scaling is conceptually different from the SCS of empirical type [11,12,14]. The second way concerns the inclusion of extra singles operator associated with the MP2 doubles operator in a second-order wavefunction, where the product of this singles operator and the CIS operator is actually evaluated for the additional correction to CIS(D) energy. These two modifications with N5 cost could be consistently combined when needed, leading to PR-CIS(Ds) by abbreviating ‘partially renormalized variant of CIS with perturbative doubles and extra singles’. Remaining of this Letter is configured as follows. In Section 2, the details of two modifications are described in the respective subsections. Test calculations with a benchmark set of small molecules, e.g., ethylene, are presented for both singlet and triplet states in Section 3. Furthermore, the cytosine molecule is taken up from a potential need for MLFMO-CIS(D) calculations [17–19] on DNA systems. 2. Spin-adapted modification 2.1. Partial renormalization of doubles amplitudes For later convenience, we first summarize the essential equations of CIS(D) [2] using mainly our spin-adapted notation [17] in the paragraph below. The indices of ijk and abc specify the occupied and virtual MOs, respectively. The explicit labeling of singlet or triplet states is omitted in this Section. A certain CIS state with excitation energy x is described as the linear combination of singly excited configuration state functions (CSFs) from the HF ground state W0 of closed-shell type, X ~ba U e a: e 1 W0 ¼ WCIS ¼ U ð1Þ i i ia

The tilde clarifies the spin-adaptation, and ~bai corresponds e 1 . The correspondto the CIS amplitudes associated with U ing CIS(D) correction is given by U 2i U 1i xCISðDÞ ¼ Eh e þ EheT 2 e ;

ð2Þ

U 2i e 2 W0 i; ¼ hWCIS jV j U Eh e

ð3Þ

U 1i e 1 W 0 ic : EheT 2 e ¼ hWCIS jV j Te 2 U

ð4Þ

e 2 yields the relaxation energy of The doubles operator U excited state, which should be a negative quantity, through the Møller–Plesset perturbation potential V as Eq. (3). This doubles correction of CIS(D) [2] is the same as that appeared in the CIS-MP2 method [1]. The spin e 2 is clearly common with U e 1 . In contrast, Te 2 state of U

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

is the MP2 doubles operator written with five types of CSFs [17] as X X X e aa e ab e aa ~ ~ ~ Te 2 W0 ¼ aaa aab aaa ii U ii þ ii U ii þ ij U ij ia

þ

i>j;a

i;a>b

X

ab ðSSÞ e ab ½~ aij ðSSÞ U ij i>j;a>b

þ

ab ðTTÞ e ab ~ aij ðTTÞ U ; ij

ð5Þ

U 1i of Eq. (4) is Now, the PR modification of EheT 2 e described. We follow the fundamental symbol and notation of Refs. [20,21] for consistency. A symbolic MP2 correlation function for an electron pair ij with spin-coupling p may be defined as X e P ¼ðij;pÞ :¼ e ab ~aab U ð9Þ ij;p U ij;p ; ab

EMP2 ¼ hW0 jV j Te 2 W0 i;

ð6Þ

Te 2 W0 :¼

X

eP ¼ U

P

where the fourth summation has two spin-coupling paths of singlet–singlet (denoted as superscripts of ‘(SS)’) and triplet–triplet (‘(TT)’). The expression of MP2 doubles amplitudes in Eq. (5) is found in Table 1 of Ref. [17]. The subscript ‘c’ in the right-hand side of Eq. (4) means the adoption of the connected part of operator contractions e 1 [2], and the disconnected part correbetween Te 2 and U MP2 sponds just to E . This factorization is a kernel of the CIS(D) method [2] which satisfies the size-consistency unlike the CIS-MP2 with non-factorized triples [1]. In fact, U 1i EheT 2 e can be considered as the differential MP2 correlation energy from the ground state, and this quantity may be positive since a hole is created in the occupied space by the single excitation. The gross CIS(D) correction of Eq. (2) can thus be determined by the balance of cancellaU 2i U 1i tion with negative Eh e . The final form of EheT 2 e evaluation is given by X U 1i ~ EheT 2 e ¼ ð7Þ bai~vai : ia

The crucial matrix ~vai is constructed through a series of integral contractions with the CIS amplitudes and also with the MP2 doubles amplitudes. The actual expression of ~vai in our spin-adapted formulation [17] is too lengthy to show here, in contrast to a simple tensorial form in the spin-orbital notation of Ref. [2], 1X a cb b ac vai ¼ hjkkbci½bbi aca ð8Þ jk þ bj aik þ 2bj aik ; 2 jkbc where caution is taken for the spin-orbital index of matrices (without tilde) and Dirac’s anti-symmetrized integral. It is noted that Grimme and Izgorodina [11] made the SCS modification [12] for the MP2 amplitudes of this equation.

391

XX ij;p

e ab ~aab ij;p U ij;p :

ð10Þ

ab

The MP2 energy can then be rewritten as the sum of pairwise correlation energies: X X ePi ¼ EMP2 ¼ hW0 jV j U eP : ð11Þ P

P

Dykstra and Davidson [20] introduced the PR technique for the MP2 energy (see Scheme I of Ref. [20]): EPR-MP2 ¼

X hW0 jV j U ePi NP

P

NP ¼ 1 þ

X

¼

X eP ; NP P

ð12Þ

e QjU e Q i; T PQ h U

ð13Þ

Q

1 T PQ ¼ ðdik þ dil þ djk þ djl Þ for P ¼ ðij; pÞ and Q ¼ ðkl; qÞ; 4 ð14Þ where NP and TPQ are the PR factor and the topological factor, respectively, which were devised by Ahlrichs et al. [21] for the CPF method as a size-consistent extension of CI singles and doubles (CISD). By setting all the TPQ in Eq. (13) to zero, the usual MP2 energy is recovered. The topological factor [21] was designed to equip both the size-consistency and orbital-invariance from the concept of type-1 coupled electron pair approximation (CEPA-1) in which excitations higher than doubles are effectively taken into account. Ref. [20] showed that the PR-modified MP2 method works better than does the naive MP2 in comparison to the CCD method as the reference, even when some near-degeneracy takes place. The application of this PR-MP2 to the CIS(D) method can be straightforward. Namely, the MP2 doubles amplitudes in Te 2 are scaled with NP, and the origin of differential correlation energy turns to EPR-MP2. In the present PR-CIS(D), we do U 2i not intend to modify Eh e since the topological factor

Table 1 Excitation energies and decomposition of CIS(Ds) corrections (in a.u.) for the lowest singlet and triplet states of carbon monoxidea ~

~

~

State

CIS

PRb

CIS(Ds)

E hU 2 i c

E hT 2 U 1 i

E hT 1 U 1 i

Singlet

0.328753

No Yes

0.324163 0.322389

0.092407

0.091331 0.089474

0.003514 0.003432

Triplet

0.211411

No Yes

0.234406 0.232769

0.057704

0.085060 0.083324

0.004360 0.004261

a

The bond length of C–O was optimized by the MP2/6-31G* procedure with the GAUSSIAN program [26,27] under no frozen-core restriction. The MP2 correlation energy (EMP2) was 0.293691 a.u. The excited state calculations were performed also with the 6-31G* basis set [27]. The states of interest were 1,3 P of the rp* transition with twofold degeneracy. b The partial renormalization [20] yielded 0.287691 a.u. as EPR-MP2. ~ c EhU 2 i was kept unchanged in the PR modification of CIS(D) correction.

392

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

[21] of Eq. (14) is defined for the ground state. It should be noted here that Piecuch and Kowalski et al. [22,23] utilized a similar way of renormalization to develop a variety of CCSD-based calculations which are applicable to systems having substantial near-degeneracies: for example, the completely renormalized variant of equation-of-motion CCSD with perturbative triples (CR-EOMCCSD(T) [23]) was usable even for the ozone molecule. Finally, the computational processing of PR-CIS(D) is addressed although the description of integral-direct implementation with parallelism [17] is out of purpose in this Letter. The NP list is prepared at an obvious cost of N5 just as the usual MP2 calculation. This preparation should be done before the construction of ~vai starts since the completed NP is used as the scaling factor of Te 2 .

The second way of modification utilizes the MP2 singles operator associated with the doubles operator as below. With the spin-adapted notation [17], an approximate second-order wavefunction of the ground state [3,24] may be written as ð15Þ

where the costly doubles–doubles operator is neglected as clarified with the dash symbol. The MP2 singles operator Te 1 is then given by X e a; ~ Te 1 W0 ¼ ð16Þ aai U i ia

The actual processing of ~aai construction would be performable with an integral-direct parallelism [15,17] when a full U 1i implementation is pursued. In a similar way to EheT 2 e , e the energy correction from T 1 is evaluated as ð21Þ

where the corresponding disconnected part is a zero quantity for the ground state by the Brillouin theorem. The new term of Eq. (21) may present an extra relaxation energy induced by the MP2 correlation for excited state, and the quantity can be positive or negative depending on the state of interest. The modified CIS(D) correction with extra singles contribution, CIS(Ds), is thus given by e ee ee xCISðDsÞ ¼ Eh U 2 i þ Eh T 2 U 1 i þ Eh T 1 U 1 i :

ð22Þ

By introducing ~xai just like ~vai , Eq. (21) is rewritten in a product form X ee ~ba~xa : Eh T 1 U 1 i ¼ ð23Þ i i ia

e a jV j Te 2 W0 i hU ~ai ¼  i a ; Dai

ð17Þ

Dai ¼ ea  ei :

ð18Þ

Eq. (18) is the orbital energy difference due to single excitation. The construction of ~ aai can be performed at the cost of 5 N , as addressed later. It is notable that ~ aai commits a sort of relaxation induced by the MP2 correlation although its energy correction for the ground state is zero because of Brillouin’s theorem, hW0 jV j Te 1 W0 i ¼ 0:

ðk>iÞ ðbc 2 j>k;b ) ðc>aÞ 1 X abðSSÞ caðSSÞ ½ðjc; ikÞ þ ðkc; ijÞ~ajk þðkb; ijÞ~ajk þ pffiffiffi : ð20Þ 2 j>k;c

ee e 1 W0 i ; Eh T 1 U 1 i ¼ hWCIS jV j Te 1 U c

2.2. Inclusion of singles contribution

WMP20 ¼ ð1 þ Te 2 þ Te 1 ÞW0 ;

ab

lists of (ia, bc) and (ia, jk) [24]. For example, the ~ aij ðSSÞ part 5 is evaluated by a series of N accumulations as ( ðjc

ð19Þ

In Ref. [2], Head-Gordon et al. omitted the MP2 singles in making the CIS(D) method from the CCSD excited states equations, by expecting that the amount of energy correction could be small. This contribution was later incorporated in a part of the CIS(D3) extension in which higher interactions such as doubles–doubles were rather focused on [3]. Unfortunately, the solo contribution from the MP2 singles has not been illuminated, to author’s knowledge, in spite of a potentially attractive feature with N5 cost to evaluate its additional correction for CIS(D) energy [2,17]. We therefore highlight this contribution. The matrix of ~ aai amplitudes is constructed by the contraction between the MP2 doubles and the MO integral

After some algebra, the final expression of ~xai is obtained as 1 X ~xai ¼ pffiffiffi ½ðjc; abÞ þ 2ðjb; acÞð1  SÞ~bbj ~aci 2 jbc 1 X þ pffiffiffi ½2ðjc; abÞ  ðjb; acÞ~bbi ~acj 2 jbc 1 X  pffiffiffi ½ðkb; ijÞ þ 2ðjb; ikÞð1  SÞ~bbj ~aak 2 jkb 1 X  pffiffiffi ½2ðkb; ijÞ  ðjb; ikÞ~baj ~abk ; ð24Þ 2 jkb where the singlet and triplet cases take S = 0 and S = 1, respectively. The ~xai matrix would be constructed with the so-called Fock-like contraction technique of N4 scaling [1,9,15,17] once the ~aai list is available. The PR modification by scaling the MP2 doubles amplitudes of Te 2 , described in the previous Subsection, is applicable in a consistent manner, leading to PR-CIS(Ds). 3. Results and discussion The proposed modifications were implemented as a prototype code on a single workstation, in order to test the numerical performance by using a benchmark set of mole-

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

cules. This prototype PR-CIS(Ds) code was interfaced with the GAMESS program [25]. The atomic orbital (AO) integral lists and the AO-MO coefficients of the HF ground state were imported from GAMESS via files, where C1 symmetry was imposed throughout. Prior to proceeding the benchmark molecules [2], we first examined whether the modifications do not degrade the desirable properties of size-consistency or orbitalinvariance equipped in the original CIS(D) method [2]. Non-interacting carbon monoxide molecules, (CO)n (n = 1, 2, 3), were used for this purpose. The C–O length was optimized at the MP2/6-31G* procedure with the GAUSSIAN program [26,27] under no frozen-core restriction. The lowest rp* states of singlet and triplet were calculated with the same setting of basis set and correlation space. Table 1 lists the energy results for carbon monoxide ‘monomer’. The decomposition of CIS(Ds) correction indicates U 1i that the amount of EheT 1 e (about 0.1 eV) is much smaller he U 2i he T 2e than those of E and E U 1 i , as expected from the perturbation order. For the singlet state, the cancellation U 2i U 1i between Eh e and EheT 2 e is coincidentally large, and then he T 1e U 1i E dominates the gross correction. The triplet case does not show such a large cancellation, and the extra sin-

393

gles correction is supplementary. Actually, each value in this table has twofold degeneracy within six decimal places because of the p symmetry in a linear molecule. In ‘dimer’ or ‘trimer’, the MOs could delocalize over the CO units and also the px,y-components could mix each other under C1 symmetry. Even for these conditions, the CIS(Ds) energy showed fourfold degeneracy and sixfold degeneracy for ‘dimer’ and ‘trimer’, respectively, regardless of the PR scaling. These facts confirm that the size-consistency and orbital-invariance are maintained in the present modifications. Foresman et al. applied the CIS and CIS-MP2 calculations [1] to ethylene [28] and butadiene [29]. Head-Gordon et al. [2,3,8,9] and Hirata [7] employed these molecules to check how their respective methods of excited states work in comparison with the CCSD-based reference values. We followed this option of testing. Table 2 shows the calculated excitation energies of ethylene and butadiene, where details on the calculation condition are found in captions. The states shown here are selected according to the adequacy of the CIS zeroth-order description indicated in Refs. [7–9]. Since the triplet CCSD values have not been available, the experimental data compiled in Refs. [28,29] are used as reference: the CIS-MP2 values are omitted.

Table 2 Calculated excitation energies (in eV) of selected lower singlet and triplet excited states of ethylenea and butadieneb PRc

No d

Yes

State

CIS

CIS(D)

CIS(Ds)

CIS(D)

CIS(Ds)

CCSDe

Ethylene Singlet 1B3u (V) 1B1g (V) 1B1u (V) 1B2g (V) 2Ag (R)

7.13 7.71 7.74 7.86 8.09

7.21 7.84 8.04 7.86 8.18

7.12 7.76 8.00 7.76 8.08

7.16 7.80 7.99 7.81 8.13

7.07 7.71 7.95 7.71 8.04

7.31 7.96 8.14 7.99 8.34

Triplet 1B1u (V) 1B3u (V) 1B1g (V) 1B2g (V) 1Ag (R)

3.57 6.91 7.63 7.75 7.77

4.57 7.11 7.81 7.79 8.09

4.52 7.03 7.72 7.70 8.01

4.51 7.06 7.76 7.74 8.04

4.46 6.98 7.67 7.65 7.96

Butadiene Singlet 1Bg (R) 1Bu (V) 1Au (R)

6.11 6.21 6.45

6.11 6.29 6.44

5.98 6.24 6.30

6.05 6.22 6.38

5.92 6.17 6.25

Triplet 1Bu (V) 1Ag (V) 1Bg (R)

2.62 4.33 6.00

3.51 5.29 6.08

3.46 5.26 5.95

3.44 5.22 6.02

3.39 5.19 5.89

Experimentalf

7.11 7.80 7.60 8.01 8.29 4.36 6.98 7.79 8.15

6.20 6.42 6.53

6.22 5.91

3.20 4.90 6.1

a The calculations of ethylene were performed with the 6-311(2+,2+)G** basis set at the MP2/6-31G* geometry [2,26–28]. Our prototype code reproduced the CIS excitation energies of Ref. [28] and also the MP2 correlation energy and singlet CIS(D) energies of Ref. [2]. b The calculations of butadiene were performed with the 6-311(2+)G* basis set at the MP2/6-31G* geometry [2,26,27,29]. The numerical data of Refs. [2,29] were reproduced by our prototype code, as in the case of ethylene. Frozen-core restriction was imposed for butadiene. c Flag of the partial renormalization. The modified MP2 energies for ethylene and butadiene were 0.327004 a.u. and 0.527592 a.u., respectively. d Label and character (‘V’ for valence type and ‘R’ for Rydberg type) of states are taken from Refs. [2,28,29]. e The CCSD-based reference values for singlet states are taken from Ref. [2]. f Experimental data were taken from the compilation of Refs. [28,29].

394

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

For the valence (V) type states, the modification as CIS(Ds) yields a preferable direction of corrections toward the experimental data for both molecules in most cases of singlet and triplet. The PR modification of Te 2 are also effective. The correspondence by the present modifications to experimental values is seemingly better than that by the reference CCSD for singlet. It should, however, be kept in mind that this situation might be accidental due to cancellation of errors, since various higher terms existing in the excited state CCSD equations [4–6,22] have been neglected in the approximation of CIS(D) [2] and also CIS(Ds). Such intrinsic limitations are made apparent for the Rydberg (R) type excitations. Namely, the excitation energies with CIS(D) corrections are already too low relative to the CCSD values and the experimental data [2]. The modifications do not play a positive role. The results of formaldehyde and acetaldehyde [2,30] are listed in Table 3. For the Rydberg states, the poor performance of CIS(D) becomes more prominent for these carbonyl compounds than for ethylene and butadiene in U 1i Table 2. The EheT 1 e correction has a right direction, but its amount is too small. The PR modification shows rather negative results. Essentially, the discrepancy is attributed to U 2i an overshoot of Eh e lowering [2,17]. This implies that the orbital relaxation in the ‘ion-like core’ of Rydberg state is hard to be satisfactorily treated at the second-order level,

in analogy with the propagator calculations for ionization potential [24]. Head-Gordon et al. [3] and Hirata [7] certainly obtained much improved values of 7.09 eV by CIS(D3) and 7.39 eV by CIS(4)p, respectively, for the singlet B2 state of formaldehyde as a representative case. In contrast, the CIS(D) values [2] themselves are well for the lower valence states of singlet and triplet, and additional corrections are obtained by the CIS(Ds) or PR-CIS(Ds) in some cases. The results of four benchmark molecules in Tables 2 and 3 [2,28–30] lead to an observation that the present modifications would nonetheless provide some improvements for low-lying valence excited states though they are not always positive. We thus recommend the concomitant use of the original CIS(D) values and the modified values by the contribution from extra singles and also by the PR scaling. In such a way, the validity of CIS(D) calculations might be checked with a range of corrected energies. Gran˜a et al. [31] reported a spectral study by the CIS(D) method [2] for the nitric acid and chlorine nitrate, which are important in atmospheric chemistry. We took up these molecules as inorganic examples having an inherent ionicity which might be overestimated at the HF level of the ground state. In Ref. [31], the lowest triplet CIS(D) energies with 6-31G* [27] were compared to the results of quadratic CISD (QCISD) and its perturbative triples

Table 3 Calculated excitation energies (in eV) of selected lower singlet and triplet excited states of formaldehydea and acetaldehydeb PRc

No d

Yes

State

CIS

CIS(D)

CIS(Ds)

CIS(D)

CIS(Ds)

CCSDe

Experimentalf

Formaldehyde Singlet 1A2 (V) 1B2 (R) 2B2 (R)

4.48 8.63 9.36

3.98 6.44 7.26

4.05 6.50 7.31

3.91 6.41 7.23

3.98 6.47 7.28

3.95 7.06 7.89

4.07 7.11 7.97

Triplet 1A2 (V) 1A1 (V) 1B2 (R)

3.67 4.64 8.28

3.49 6.04 6.47

3.55 6.03 6.52

3.43 5.96 6.44

3.48 5.95 6.49

Acetaldehyde Singlet 1A00 (V) 2A 0 (R) 3A 0 (R)

4.89 8.51 9.22

4.28 6.13 7.04

4.36 6.16 7.06

4.21 6.10 7.00

4.30 6.13 7.02

4.15 4.99 8.24

3.86 6.21 6.17

3.93 6.22 6.19

3.80 6.14 6.14

3.86 6.15 6.16

Triplet 1A00 (V) 1A 0 (V) 2A 0 (R) a

3.50 5.86 6.83

4.26 6.78 7.49

4.28 6.82 7.46 3.97 5.99 6.81

The calculations of formaldehyde were performed with the 6-311(2+,2+)G** basis set at the MP2/6-31G* geometry [2,26,27,30]. The numerical data of Refs. [2,30] were reproduced by our prototype code for singlet states. The triplet CIS energies were not presented in Ref. [30]. b The calculations of acetaldehyde were performed with the 6-311(2+)G* basis set at the MP2/6-31G* geometry [2,26,27,30]. The numerical data of Refs. [2,30] were reproduced by our prototype code for singlet states. c Flag of the partial renormalization. The modified MP2 energies for formaldehyde and acetaldehyde were 0.380253 a.u. and 0.525308 a.u., respectively. d Label and character (‘V’ for valence type and ‘R’ for Rydberg type) of states are taken from Refs. [2,30]. e The CCSD-based reference values for singlet states are taken from Ref. [2]. f Experimental data were taken from the compilation of Ref. [30].

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

395

Table 4 Calculated excitation energies (in eV) of lowest singlet and triplet states of nitric acid and chlorine nitratea PRb

No c

Yes

State

CIS

CIS(D)

CIS(Ds)

CIS(D)

CIS(Ds)

Ref. value

Nitric acid Singlet 1A00 (np*) 2A 0 (pp*)

5.76 8.20

4.89 7.39

5.06 7.66

4.79 7.28

4.95 7.55

4.99d 7.49d

Triplet 1A 0 (pp*) 1A00 (np*)

3.60 5.12

4.73 4.59

4.82 4.73

4.62 4.50

4.70 4.64

4.33 (5.04)e 4.46 (4.96)e

Chlorine nitrate Singlet 1A00 (nr*) 2A 0 (nr*)

4.83 5.70

4.53 5.43

4.50 5.43

4.48 5.38

4.44 5.38

4.46d 5.33d

Triplet 1A00 (nr*) 1A 0 (nr*)

3.46 4.25

3.50 4.41

3.49 4.42

3.44 4.36

3.43 4.36

3.41 (3.43)e 4.52 (4.79)e

a

The calculations of nitric acid and chlorine nitrate were performed with the 6-31G* basis set [27] at the QCISD/6-31G* geometry [26,27,31]. Frozencore restriction was imposed. The CIS and CIS(D) values of Ref. [31] were reproduced by our prototype code for both singlet and triplet states. b Flag of the partial renormalization. c Label and transition type of states are taken from Ref. [31]. d The reference values for singlet states were taken from the CCSD results in Ref. [31]. e The QCISD values [31], whose reference determinant is of UHF type, are adopted for triplet states (see text). The more correlated QCISD(T) values are also cited in parenthesis.

correction (QCISD(T)) [32] based on the unrestricted HF (UHF) determinant. We adopted such a strategy and used the same condition for excited states calculations [31]. Table 4 presents the lowest excitation energies of singlet and triple states. In comparison with the reference values of CCSD for singlet and of QCISD and QCISDT(T) for triplet, it is found that the proposed CIS(Ds) yields certain improvements (more than 0.1 eV for singlet) over the unmodified CIS(D) for the nitric acid, demonstrating an effect of MP2-induced relaxation via Te 1 . The PR modification works too. Similar improvements by PR-CIS(Ds) are found for the chlorine nitrate. The results in Table 4 imply a potential that the modified CIS(D) recipes are useful in estimating lower excitation energies of inorganic molecules as long as the issue of near-degeneracy [3,8,9] is not severe. Manifestly, the excited states of base molecules in DNA have attracted great interest in photobiochemistry from both experimental and theoretical sides. Valiev and Kowalski [33] recently performed an interesting study in which a model DNA system was simulated by a classical MD scheme and the sampled configurations of a cytosine moiety were then subjected to the equation-of-motion type CCSD calculations (EOMCCSD and CR-EOMCCSD(T) [22,23]) with Dunning’s cc-pVDZ basis set [34] under the presence of environmental point charges. They finally obtained the lowest pp* and np* energies of singlet as the statistically averaged values, and these values were compared to the respective energies of free cytosine molecule [33]. We have been interested in modeling the DNA systems by the MLFMO-CIS(D) methodology [17,18] fur-

ther including the configuration fluctuation with MD [19]. We thus employed the cytosine molecule as a pilot test. The lowest triplet states were calculated as were the singlet states, since the singlet–triplet energy gap of cytosine has been of interest [35]. The basis sets of cc-pVDZ [34] and 6-31G* [27] were used in a series of CIS(D) calculations. The reference values were calculated by using the EOMCCSD module [22,23] available in GAMESS [25]. Table 5 summarizes these calculated values, where the direct comparison to experimental data for cytosine is not the present purpose. For the singlet states, the modification of CIS(Ds) leads to a better correspondence to the EOMCCSD reference value. Ref. [31] discussed the importance of doubles contributions in CCSD expansions and then claimed the advantage of CR-EOMCCSD(T) [22,23] over EOMCCSD. However, the shift in excitation energies as the environmental effect (‘free’ versus ‘in DNA’) was well evaluated for both pp* and np* states at the EOMCCSD level [33]. This would be encouraging by considering the fact that the present recipes of CIS(D) modification reproduce the reference values reasonably, as shown in Table 5. For the triplet pp* state, the PRCIS(Ds) procedure yields the excitation energies of 4.02 eV with cc-pVDZ and 4.09 eV with 6-31G*. Abouaf et al. [35] estimated an energy range of 3.75–3.97 eV by the UHF-based CCSD calculations for the lowest triplet excitation of vertical process. The range reported in Ref. [35] is comparable to corresponding our estimates. As a whole, the present attempt on cytosine could be considered as a touchstone for future MLFMO-CIS(D) calculations [17–19] on DNA systems.

396

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397

Table 5 Calculated excitation energies (in eV) of lowest singlet and triplet states of cytosinea PRb

No

Yes

State

CIS

CIS(D)

CIS(Ds)

CIS(D)

CIS(Ds)

CCSDc

cc-pVDZ Singlet 2A 0 (pp*) 1A00 (np*)

6.11 6.90

4.99 5.30

5.05 5.35

4.90 5.23

4.96 5.28

5.05 5.49

Triplet 1A 0 (pp*) 1A00 (np*)

3.64 5.12

4.11 4.87

4.10 4.87

4.03 4.78

4.02 4.79

6-31G* Singlet 2A 0 (pp*) 1A00 (np*)

6.20 6.92

5.02 5.34

5.10 5.40

4.93 5.27

5.01 5.33

Triplet 1A 0 (pp*) 1A00 (np*)

3.70 5.13

4.18 4.87

4.17 4.90

4.10 4.79

4.09 4.82

5.09 5.53

a

Two types of basis sets, cc-pVDZ [34] and 6-31G* [27], were used for excited state calculations. The molecular geometries were optimized with the respective basis sets at the MP2 level by using the GAUSSIAN program [26]. Frozen-core restriction was imposed throughout. b Flag of the partial renormalization. c The EOMCCSD module [22,23] available in GAMESS [25] was used to obtain the reference values for singlet states. In Ref. [33], the EOMCCSD values with cc-pVDZ basis [34] were reported as 5.02 eV for 2A 0 state and 5.44 eV for 1A00 state, at a density functional theory-optimized geometry (see this original paper for details).

4. Summary

Acknowledgements

In this Letter, we proposed a couple of modifications for the spin-adapted version [17] of CIS(D) method [2], based on the partial renormalization for the Te 2 operator of MP2 doubles [20] and also on the inclusion of the contribution from the Te 1 operator of extra singles [24]. The combined modifications could be termed PR-CIS(Ds) which maintains the computational cost of N5 scaling. With a prototype code, a series of test calculations were performed for both singlet and triplet states on a benchmark set of molecules: ethylene [28], butadiene [29], formaldehyde and acetaldehyde [30]. The nitric acid and chlorine nitrate [31] were employed as inorganic examples. Finally, the cytosine molecule [33,35] was calculated from our interest in modeling DNA systems. These calculated excitation energies were compared with the available CCSD reference values and experimental data. Comparisons showed that some improvements could be obtained through the proposed modifications for low-lying valence states, although negative cases were also observed. In other words, the calculations of CIS(Ds) or PR-CIS(Ds) would provide supplementary values by which the validity of CIS(D) is monitored with a range of corrected energies. A full implementation with parallelism in the ABINIT-MP program [15,16] has been under way at the enhancement of MLFMO-CIS(D) ability [17–19]. Further calculations with a wider variety of molecules would be interesting to investigate the performance of the modified CIS(D) recipes in more detail.

YM would thank Prof. Shigenori Tanaka, Dr. Tatsuya Nakano, Dr. Takeshi Ishikawa and Dr. Kaori Fukuzawa for their interest and encouragement. YM would also be grateful to the anonymous referees for comments to improve the manuscript. The work reported here was supported primarily by the ‘Revolutionary Simulation Software for 21st Century’ (RSS21) project operated by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) and partially by the CREST project operated by the Japan Science and Technology Agency (JST). YM would also be grateful for the S.F.R. aid by Rikkyo University. References [1] J.B. Foresman, M. Head-Gordon, J.A. Pople, M.J. Frisch, J. Phys. Chem. 96 (1992) 135. [2] M. Head-Gordon, R.J. Rico, M. Oumi, T.J. Lee, Chem. Phys. Lett. 219 (1994) 21. [3] M. Head-Gordon, T.J. Lee, in: R.J. Bartlett (Ed.), Recent Advances in Coupled Cluster Methods, World Scientific, Singapore, 1997. [4] H. Nakatsuji, Chem. Phys. Lett. 59 (1978) 362. [5] H. Koch, P. Jørgensen, J. Chem. Phys. 93 (1990) 3333. [6] J.F. Stanton, R.J. Bartlett, J. Chem. Phys. 98 (1993) 7029. [7] S. Hirata, J. Chem. Phys. 122 (2005) 094105. [8] M. Oumi, D. Maurice, T.J. Lee, M. Head-Gordon, Chem. Phys. Lett. 279 (1997) 151. [9] M. Head-Gordon, M. Oumi, D. Maurice, Mol. Phys. 96 (1999) 593. [10] O. Christiansen, H. Koch, P. Jørgensen, Chem. Phys. Lett. 243 (1995) 409.

Y. Mochizuki, K. Tanaka / Chemical Physics Letters 443 (2007) 389–397 [11] [12] [13] [14] [15]

[16] [17] [18] [19] [20] [21] [22] [23]

S. Grimme, E.I. Izgorodina, Chem. Phys. 305 (2004) 223. S. Grimme, J. Chem. Phys. 118 (2003) 9095. L. Meissner, Mol. Phys. 104 (2006) 2073. P.-D. Fan, S. Hirata, to be published. E.B. Starikov, J.P. Lewis, S. Tanaka, Modern Methods for Theoretical Physical Chemistry of Biopolymers, Elsevier, Amsterdam, 2006 (FMO documents written by Fedorov et al. and Nakano et al. therein). T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, K. Kitaura, Chem. Phys. Lett. 351 (2002) 475. Y. Mochizuki et al., Theor. Chem. Acc. 117 (2007) 541. Y. Mochizuki, T. Nakano, S. Amari, T. Ishikawa, K. Tanaka, M. Sakurai, S. Tanaka, Chem. Phys. Lett. 433 (2007) 360. Y. Mochizuki, Y. Komeiji, T. Ishikawa, T. Nakano, H. Yamataka, Chem. Phys. Lett. 437 (2007) 66. C.E. Dykstra, E.R. Davidson, Int. J. Quant. Chem. 78 (2000) 226. R. Ahlrichs, P. Scharf, C. Ehrhardt, J. Chem. Phys. 82 (1985) 890. P. Piecuch et al., Theor. Chem. Acc. 112 (2004) 349, references therein. K. Kowalski, P. Piecuch, J. Chem. Phys. 120 (2004) 1715.

397

¨ hrn, Propagators in Quantum Chemistry, second [24] J. Linderberg, Y. O edn., Wiley, Hoboken, NJ, 2004. [25] M.W. Schmidt et al., J. Comp. Chem. 14 (1993) 1347, GAMESS, (http://www.msg.ameslab.gov/GAMESS/). [26] M.J. Frisch et al., GAUSSIAN (http://www.gaussian.com). [27] J.B. Foresman, A. Frisch, Exploring Chemistry with Electronic Structure Methods, second edn., Gaussian Inc., Pittsburgh, 1996. [28] K.B. Wiberg, C.M. Hadad, J.B. Foresman, W.A. Chupka, J. Phys. Chem. 96 (1992) 10756. [29] K.B. Wiberg, C.M. Hadad, B. Ellison, J.B. Foresman, J. Phys. Chem. 97 (1993) 13586. [30] C.M. Hadad, J.B. Foresman, K.B. Wiberg, J. Phys. Chem. 97 (1993) 4293. [31] A.M. Gran˜a, T.J. Lee, M. Head-Gordon, J. Phys. Chem. 99 (1995) 3493. [32] J.A. Pople, M. Head-Gordon, K. Raghavachari, J. Chem. Phys. 87 (1987) 5968. [33] M. Valiev, K. Kowalski, J. Chem. Phys. 125 (2006) 21101. [34] T.H. Dunning, J. Chem. Phys. 90 (1989) 1007. [35] R. Abouaf, J. Pommier, H. Dunet, P. Quan, P.-C. Nam, M.T. Nguyen, J. Chem. Phys. 121 (2004) 11668.