WIXLIB

Rij ( ja kb) aikab ,(9)jkabwia ( jk bc ) aikac bbj .(10)jkbcEq 6 defines the second-order correction to the CIS excitation energy, CIS(D) , leading atotal excitation energy that is CIS CIS(D) .B. Resolution-of-the-identity in CIS(D) Theory.Let us first introduce theauxiliary basis as a resolution-of-the-identity (RI) approximation in the CIS(D) theory.The RI approximation describes all electron-repulsion integrals (ERIs) in eqs 7 – 10 as( pq rs ) RI ( pq P)( P Q ) 1(Q rs )PQ BRpq(11)BrsR ,Rwith the B matrix defined asRB pq ( pq P )( P R ) 1/2 .(12)PFrom a computational point of view, it is advantageous to define three other relatedquantities:VaiR bib BabR ,(13)bOaiR baj BijR ,(14)jDaiR VaiR OaiR .(15)With these definitions, it is easy to show that eq 7 can be transformed intouijab DaiP BbjP DbiP BajP BaiP DbjP BbiP DajP .P9(16)

Eqs 8 – 10 can be obtained by combining the ERI (ia jb) RI and the amplitude aijab .An efficient algorithm for this resolution of the identity formulation of CIS(D)theory is presented in Figure 1, where one can easily see that its cost scales with the fifthpower of system size. In addition, the disk transfer cost is fourth order, with the size ofthe storage space requirement scaling as third order. In this algorithm, a batching schemeis introduced to minimize the cost for disk input-output (I/O), especially for DaiR (Line 4).The I/O cost decreases with a larger batch size, and the maximum batch size can be easilycalculated from the size of available memory and disk space. One important point is thatfor a calculation of S excited states, the total cost grows as 2O 2V 2 XS .The algorithmpossesses three additional fifth order steps related to the computation of (ia jb) RI , Rab ,and Rij , but these do not depend on S.It is interesting to note the possibility of a minor modification of the abovealgorithm. Based on the formal similarity to RI-MP2 gradient theory,31,32 eqs 8 and 9 canbe rearranged asPRab BbkP ak,(17)kPRij BajP aiP ,(18)aPwhere, following RI-MP2 gradient theory, the three-center two-particle density matrix is: aiP aijab ( jb Q )( P Q ) 1/2 .(19)jbQThis leads to an alternative working expression for wia as well:wia BbjP aiP bbj ( jc kb) aikacbbj .jbPjkbc10(20)

The potential benefit of using aiP will be reduced disk IO cost. While eqs 8 and 9require fourth order disk IO related to the storage of ERIs and the a-amplitudes (Line 17and 18 in Figure 1), the use of aiP would only require third order disk access for thecalculations of Rab and Rij . However, the CPU cost of calculating eq 19 will be O2V2X.In fact, this is larger than the combined cost of eqs 8 and 9 (O2V2N). This additional CPUtime will become more important as the system size grows.Therefore, it is moredesirable to generate Rab and Rij based on the RI-approximated integrals without using aiP . However, the use of aiP will be crucial in the efficient implementation of the SOSCIS(D) theory as will be shown later.C. The SCS-CIS(D) Method. We define SCS-CIS(D) theory in an analogy tothe manner in which Grimme first proposed the corresponding ground state SCS-MP2method, by scaling the same-spin and opposite-spin components of the energy.Inaddition to this split spin component treatment, an empirical damping factor 0 1 forthe CIS excitation energy is introduced for the direct term as ab V U1 0 ijab V U1 0 11ijab,(21a)U 0 ij ab4 ijab a b i j 4 ijab a b i j ijSS2U 0 OS2 ab V U1 0 ijijab a b i j . abij (21b)We will detail the role of in a later section. With the obvious components of theindirect term as in SCS-MP2, the spin-component scaling modification of eq 3 willbecome:11

E SCS-CIS(D) CIS V cUOSU 2OS cUSSU 2SS 0 CIS V cTOST2OS cTSST2SS U1 0 .(22)With the independent scaling of U2-term and the use of the damping factor ( ), thisequation differs significantly from a previous suggestion for defining SCS-CIS(D),27which left the first term of eq 3 unmodified, but replaced the second term as we havedone in the above.When the ground state correlation energy contribution is separated as in eq 6, theSCS correlation correction to the CIS energy can be written as SCS-CIS(D) cUOS wUOS cTOS wTOS cUSS wUSS cTSS wTSS ,(23)with the obvious definition for each of the terms. In practice, the opposite-spin (OS) andthe same-spin (SS) component splitting of the U operator can be performed without anyadditional computational cost during the first summation in eq 6. In contrast, the splittingof the T operator requires separate evaluations of OS and SS contributions to Rab , Rij ,and wia . In closed-shell systems, this is attained at an additional cost of O 2V 2 N (eqs 8and 9). However, compared to the leading cost of 2O 2V 2 XS , this additional cost of SCSCIS(D) is negligible especially when excitation energies to multiple states are calculatedat the same time.D. The SOS-CIS(D) Theory. The opposite spin part of the CIS(D) correctioncan be extracted from eq 22 asE OS-CIS(D) CIS V U 2OS 0 CIS V T2OSU1 0 .12(24)

From the symmetry of uijab with respect to the permutation of indices, it is easy to showthat the first term in this equation becomes CIS V U 0 OS2ijab u abij2 a b i j ,(25)where we again use the empirical damping factor 0 1 , and haveuijab DaiP BbjP BaiP DbjP .(26)PLikewise, the second term in eq 24 can be expressed as CIS V T2OSU1 0 E OS-MP2 bia bib Rab bia bib Rabiabiab b b Rij b b Rij bia wia bia wia ,c ci jijccicjijcia(27)iawithRab ( jc kb) a cajk ,(28)jkcRij ( ja kb ) aikab ,(29)kabwia ( jb kc ) aacikjkbcbbj ( jk bc ) aikac bbj ,(30)jkbcThe beta spin intermediates, Rab , Rij , and wia are defined analogously. Also, E OS-MP2denotes the opposite spin component of MP2 correction.24This opposite-spin formalism can be transformed into a fourth order algorithm24through the use of a Laplace transform with discrete numerical quadratures26 x 1 dt exp( xt ) t t exp( xt )0in conjunction with the RI approximation. Firstly, eq 25 can be transformed as13(31)

wI t DaiP BbjP BaiP DbjP DaiQ BbjQ BaiQ DbjQ et ( a b i j ) tPQ ijab(32) t DaiP BbjP BaiP DbjP D aiQ B bjQ B aiQ D bjQ e t ,tPQ ijabwith B aiR BaiR e( i a ) t and D aiR DaiR e( i a ) t . Let us denote terms from eq 27 involving Rijand Rab as wII :wII bia bib Rab bia bib Rab bicb cj Rij bic bcj Rij .iabiabijc(33)ijcIn addition, it is easy to show that the last terms in eq 27 that involve wia can beexpressed as wIII t ( f P f P ) f Q BckP B ckQ BcjPV jkP B ckQ f Q tjkc PQ kc PQ ( f P f P ) f Q B B B V B f Q ,kc PQjkc PQ PckQckPcjPjkQck(34) with the definitions: VijR a bia BajR , f R ai bia BaiR , f R ai bia B aiR and their analogsin the -spin space.A fourth order algorithm can be implemented by carefully rearranging the orderof summations in various terms. When X, Y, and Z are defined as X PQ BaiP B aiQ ,(35)ai YPQ DaiP B aiQ D aiP BaiQ ,ai(36)ai Z PQ DaiP D aiQ ,(37)aitogether with their obvious -spin analogs, the first OS-CIS(D) correction term becomes wI t e t X PQZ PQ YQPYPQ YPQYQP Z PQX PQ .tPQ14(38)

Also, in analogy to the RI-CIS(D) case, when aiP is introduced as aiP t B aiQ X PQ,t(39)Qit is trivial to show that eqs 28 and 29 are equivalent toRab aiP BbiP ,(40)iPRij BajP aiP .(41)aPFinally, when G and H matrices are defined asGck BcjPV jkP ,(42)jPH ck B ckP f P ,(43)jPthe last term of the OS-CIS(D) correction becomes Gck H ckwIII t ( f P f P ) f Q X PQtkc PQ ( f P f P ) f Q X PQ Gck H ck .PQkc (44)By collecting the above expressions, the scaled opposite spin CIS(D) excitation energiesare obtained as SOS-CIS(D) CIS V cUU 2OS 0 CIS V cT T2OSU1 0 cU wI cT ( wII wIII ),(45)with two empirical scaling parameters, cU and cT , where the latter is already fixed fromthe ground state SOS-MP2 energy, and the former is to be determined by comparingagainst either higher accuracy calculations or experiments.15

Of the various working expressions listed in the above, only eqs 36 and 37 need tobe evaluated for each excited state and each Laplace quadrature point. Accordingly, thecomputational cost of this method will be dominated by the evaluation of these twoequations, requiring a total of 2OVX2ST operations, with S and T denoting the numbers ofexcited states and quadrature points, respectively. The resulting overall algorithm isshown in Figure 2. Comparing against the ground state SOS-MP2 method, we concludethat the cost per state (for S not too small) will be approximately twice the cost of thecorresponding ground state SOS-MP2 calculation. Also, it should be noted that theLaplace transform in our algorithm does not require any aggressive integral screeningscheme,33 which is practically required for an efficient treatment of the same-spincomponent calculation.III. Optimizations of ParametersAs shown in the previous section, the proposed methods require optimization ofvarious parameters. The most straightforward way will be to use experimental data in thedetermination of these parameters. In this work, the extensive set of organic moleculesadopted by Grimme and co-workers27 have been used again. This set only includesvalence transitions with * and n * characters. To make the training set morecomplete, we have added well-characterized Rydberg transitions of water and ammonia.(See Table 1 for the complete list of the transitions.)In the calculations of organic molecules, molecular geometries were obtained atthe HF/6-311G(d,p) level for the ground states and the CIS/6-311G(d,p) level for theexcited states. These levels of theory are roughly comparable in quality to DFT methods,16

although there is a systematic tendency to make bondlengths slightly too short, and thusvibrational frequencies slightly too high. To obtain 0-0 transition energies, correctionsfor zero-point energies must be computed for both the ground and the excited states.Frequencies obtained from analytic Hessians at the above levels of theory have been usedfor this purpose after scaling with a factor of 0.9. In the correlated excitation energycalculations at the optimized geometries, the aug-cc-pVTZ basis34 was employed togetherwith its corresponding auxiliary basis set.35The CIS and HF components of thecalculation were performed without the RI approximation.In the case of Rydberg states, 0-0 transitions may not be experimentallyobservable36 because of potentially large Franck-Condon shifts. Accordingly, we haveused vertical excitation energies37,38 for these transitions.In the excitation energycalculations, we have used 6-311(2 ,2 )G(d,p) basis together with the auxiliary basis ofaug-cc-pVTZ.34 Even though this auxiliary basis was not specifically optimized for thePople-style basis, the RI approximation error with the basis was always found to besmaller than 0.001 eV, similar to the report for the ground state energy calculations.31 Allcalculations were performed with a development version of Q-Chem 3.0.39A. Performance of RI-CIS(D). Because the present SCS- and SOS-CIS(D)theories are based on CIS(D), it is natural to look to the performance of this method to getinsight for possible improvements toward SCS- and SOS-CIS(D). Figure 3 presents theerrors (against the experimental values) of RI-CIS(D) for the molecules in the training set.Firstly, one can clearly see that there exists a systematic overestimation in the valencetransitions: the method tends to give larger transition energies than experiment as17

represented by its mean signed error (MSE) of 0.19 eV for these transitions. In addition,but more notably, the method tends to severely underestimate the transition energies ofthe Rydberg transitions. This underestimation is indeed a generic problem of the method.As the energy denominator in eq 5 becomes smaller, the magnitude of the direct termbecomes larger.Because the correlation correction from the direct term is alwaysnegative, an over-correction caused by small denominator (or large CIS) leads to thistendency of underestimation in the total transition energies. This effect predominantlyappears for Rydberg transitions because qualitatively they involve the lowest lying (mostdiffuse) virtual orbitals, and therefore the smallest energy denominators which are mostsensitive to When the new SCS- and SOS-CIS(D) methods were applied with 1(no damping), a similar defect was observed for these Rydberg transitions in the trainingset.This defect of an unbalanced description of valence and Rydberg excited stateswill be removed as one introduces higher correction terms of the coupled cluster theory.(Recall that CIS(D) may be expressed as a low order truncated solution of linear responsecoupled cluster theory.) However, such an approach is not a realistic option for ourpresent work, where the design of an efficient algorithm is under pursuit. Instead, werecall that the problem is mostly remedied when the excitation energy is iterativelycalculated in the quasi-degenerate variant CIS(D0) theory,29 where is omitted in thecalculation of the excitation amplitudes bijab (eq 5).The more balanced behavior ofCIS(D0) theory is one motivation for the introduction of the empirical damping factor .A closely related motivation is the above discussion of the difference in the important18

virtual orbitals between Rydberg excited states (low-lying) and valence states (higherlying antibonding orbitals).We stress that the use of the damping factor is likely to most improve the CIS(D)method when it is combined with spin component scaling. Because the effect of dampingwill be to decrease the direct correlation correction (in other words, the U2-term willbecome less negative), it will tend to degrade the performance of CIS(D) for valencetransitions. We aim to compensate this potential problem through the use of scalingparameters.In theory, we can test the behavior of both SCS- and SOS-CIS(D) as a function ofthe value. Because different values affect every individual component of the directterm in a different manner (Line 13 in Figure 1), such a test will require tremendouscomputational effort in the SCS-CIS(D) case. In the SOS-CIS(D) case, however, onlyone set of calculations can be used to obtain excitation energies at all different (Line 9in Figure 2). For this practical reason, we will only use SOS-CIS(D) to obtain theoptimal value of the damping factor.B. Numerical Quadratures for Laplace Transform.For the SOS-CIS(D)method, we need to define the quadrature scheme used to evaluate the Laplace transform.Here, we employed the same scheme previously reported with SOS-MP2 theory,24 theground state counterpart of SOS-CIS(D). Specifically, 10 numerical quadrature pointswere obtained by minimizing the integrated error 2xmaxxmin 1 dx t t exp( xt ) x 2(46)19

according to Wilson and Almlöf40 with xmin 0.01 and xmax 400 a.u.Thex a b i j CIS values for all the molecules tested in this work actually fell inthis range for all possible values ( 0 1 ). The contribution to the excitation energyafter the seventh quadrature point was found to be negligible (less than 0.001 eV) in alltest results as was found previously for the ground state case.24 Because the contributionfrom the seventh quadrature point appeared to be considerably smaller than the overalluncertainty level of SOS-CIS(D) (discussed later), one might consider a reduction in thenumber of quadrature points to improve the efficiency. To preserve the consistency withthe ground state description, however, we did not try this in the present work.In fact, the above scheme will not be the most efficient strategy for the numericalintegration of the Laplace transform. The best accuracy with the least number of points isexpected if the points and weights are actually determined for the given system,40potentially with two separate quadrature schemes for the direct term and the indirect term.In addition, different quadrature schemes33,41 may further reduce the computational cost.We do not consider such possibilities in this work for the following reasons. First,system-specific optimal quadrature points will surely depend on the energy eigenvaluesof canonical molecular orbitals and potentially on CIS excitation energies when is nonzero.This dependency will introduce an undesirable complication when analyticgradients of the ground and excited state surfaces are considered.42 Also, using differentquadratures for the direct and the indirect terms will result in different B matrices in theterms, increasing the associated computational cost. For the time being, our presentapproach based on simple least-square Gaussian quadrature achieves sufficient accuracyand efficiency ( 0.001 eV error with only seven quadrature points).20

C. Determination of Damping Factor. Now that the quadrature scheme isdefined, we can determine the optimal damping factor as follows. For any given value,SOS-CIS(D) method has two adjustable empirical parameters. Of the two, the parameterrelated to the indirect term with the T2 operator will be transferred from the counterpartground state theory (SOS-MP2) for consistency (namely, cT 1.3).24 This leaves onlythe parameter related to the U2 opera

of second order Møller-Plesset perturbation theory (MP2) can be systematically improved by separate scaling of same-spin and opposite-spin contributions to the correlation energies. 23. This spin-component scaling scheme of the MP2 excitation amplitudes was later also applied to CIS(D), though only on the so-called indirect term (see next .