JCP.128.094105

Formulation and implementation of direct algorithm for the symmetry- adapted cluster and symmetry-adapted cluster–configuration interaction method Ryoichi Fukuda1,2,a� and Hiroshi Nakatsuji1,2,b� 1Quantum Chemistry Research Institute,c� Kyodai Katsura Venture Plaza, North Building 106, 36-1 Goryo Oohara, Nishikyo-ku, Kyoto 615-8245, Japan 2Department of Synthetic Chemistry and Biological Chemistry, Graduate School of Engineering, Kyoto University, Kyoto-Daigaku-Kaysura, Nishikyo-ku, Kyoto 615-8510, Japan �Received 26 July 2007; accepted 14 December 2007; published online 4 March 2008� We present a new computational algorithm, called direct algorithm, for the symmetry-adapted cluster �SAC� and SAC–configuration interaction �SAC-CI� methodology for the ground, excited, ionized, and electron-attached states. The perturbation-selection technique and the molecular orbital index based direct sigma-vector algorithm were combined efficiently with the use of the sparse nature of the matrices involved. The formal computational cost was reduced to O�N2�M� for a system with N-active orbitals and M-selected excitation operators. The new direct SAC-CI program has been applied to several small molecules and free-base porphin and has been shown to be more efficient than the conventional nondirect SAC-CI program for almost all cases. Particularly, the acceleration was significant for large dimensional computations. The direct SAC-CI algorithm has achieved an improvement in both accuracy and efficiency. It would open a new possibility in the SAC/SAC-CI methodology for studying various kinds of ground, excited, and ionized states of molecules. © 2008 American Institute of Physics. �DOI: 10.1063/1.2832867� I. INTRODUCTION The symmetry-adapted cluster �SAC�1 and SAC– configuration interaction �SAC-CI�2 methodology proposed and first coded by one of the authors in 1978 is an electron correlation methodology for ground, excited, ionized, and electron-attached states of molecules. It has been success- fully applied to diverse chemistry, physics, and biology in- volving various kinds of electronic and vibrational states �for recent reviews, see Ref. 3�. The methodology is based on the cluster expansion formalism combined with the variational principle.2 Later, theoretically identical methodologies, such as coupled-cluster linear response theory4 and equation-of- motion coupled cluster,5 have been reported. These methods established a highly accurate way to study excited electronic structures of molecules. The earlier version of the SAC/SAC-CI program6 was published in 1985. A great deal of effort has been made to the theory and program development, and a lot of new ideas have been implemented. The SAC-CI general-R method,7 multireference version of the SAC/SAC-CI method,8 expo- nentially generated wave function idea,9 extension up to sep- tet spin multiplicity,10 and the analytical energy gradient method11 were the representative fruits of these efforts. Re- cently, the giant SAC/SAC-CI theory for giant molecular crystals has been proposed.12 The SAC/SAC-CI method has been implemented in the GAUSSIAN03 software program package and has been widely used in universities and industries.13 So far, the SAC/SAC-CI code has been written with the excitation-operator driven algorithm1,2 with the use of the perturbation-selection technique.14,15 The integral driven al- gorithm was introduced to the SAC/SAC-CI method by Hirao,16 but his algorithm was not combined with the perturbation-selection technique. Molecular orbital �MO� integral driven direct algorithm for electron correlation methods was first proposed by Roos for singles and doubles �SD� CI method.17 In such a proce- dure, the iteration vector �usually termed as sigma vector� is directly constructed from MO integrals, without an explicit construction of a Hamiltonian matrix. At the same time, when the SAC/SAC-CI SD code was established,1,2 the coupled-cluster doubles �CCD� method only for ground state was formulated by Pople et al. in integral driven form.18 An efficient computation algorithm of CCSD was designed by Scuseria et al.,19 who introduced intermediate arrays for ef- ficiency. These MO integral driven approaches include all SD or D excitation operators. We consider that this feature conflicts with the policy of our SAC/SAC-CI program ex- plained below. A policy of our SAC/SAC-CI program is that we discard minor unimportant terms by introducing some selection pro- cedures with thresholds. Accordingly, depending on the de- sired accuracy, we introduce appropriate thresholds of selec- tions, and less important but time-consuming terms are neglected.20 The perturbation selection of the linked opera- tors is a particularly important technique.14,15 By virtue of this technique, we can use a single theory and a single pro- gram for various research subjects;3 from fine theoretical a�Electronic mail: fukuda@qcri.or.jp. b�Electronic mail: h.nakatsuji@qcri.or.jp. c�URL: http://www.qcri.or.jp/. THE JOURNAL OF CHEMICAL PHYSICS 128, 094105 �2008� 0021-9606/2008/128�9�/094105/14/$23.00 © 2008 American Institute of Physics128, 094105-1 Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp spectroscopy of small molecules21–23 to photobiology of big biological molecules.24,25 This feature enables us to investi- gate various chemical phenomena from an equal viewpoint. The perturbation-selection technique is particularly impor- tant for excited states. The atomic orbital �AO� basis func- tions are usually optimized for the ground state. We need a flexible AO basis to describe various kinds of excited states. Diffuse functions and/or sizable Rydberg functions are often required. Therefore, the number of active MOs �N� for excited-state calculations is much larger than that for ground state. Order N6 algorithms �the optimal scale of CCSD method� will soon face the limit of computer. The limitations of O�N6� algorithms would be much more severe for excited states than for ground states. A problem in the algorithm of the conventional SAC-CI program lies in the calculations of unlinked terms. The Hamiltonian matrix elements between selected excitation op- erators are evaluated within the loops driven by the labels of excitation operators. Consequently, the unlinked integral part becomes an O�M3� step, where M is the number of selected excitation operators. This step is so time consuming that the additional approximation shall be introduced for practical calculations, and the approximation was to cut less important excitation operators than a given threshold off the loops. This cutoff approximation is, however, problematic because the unlinked terms are essentially important in the SAC/SAC-CI theory. An efficient use of the core memory is also very impor- tant. We use a projective reduction formula26,27 to evaluate the Hamiltonian matrix elements between selected excitation operators. The MO integrals are randomly requested from the projective reduction routine in the loops of the excitation operator labels. Therefore, we have to load all the MO elec- tron repulsion integrals �ERIs� on a core memory. If the core memory is not enough to load them on, we have to use a multipass algorithm for separated MO ERIs, but this multi- pass algorithm is inefficient. In this paper, a solution for these problems is provided by combining the direct algorithm with the perturbation- selection technique. We have developed the direct SAC and SAC-CI algorithm running within the selected excitation op- erators. The MO based direct algorithm that utilizes the sparseness of the Hamiltonian matrices provides an efficient algorithm that works within the selected excitation operators. This also results in an efficient minimal memory require- ment. The theoretical framework and the algorithm of the di- rect SAC/SAC-CI method are described in the next section. The implementation of the new direct SAC-CI algorithm is given in Sec. III, and applications to several test cases are reported in Sec. IV. The present article focuses on single- point calculations. The analytical energy gradient of the SAC/SAC-CI method has been given in GAUSSIAN 03. We have also adapted the direct algorithm for the analytical en- ergy gradient calculations, and it will be reported in the subsequent paper. II. THEORY A. SAC/SAC-CI method The details of the SAC/SAC-CI methodology have been reviewed in several articles.1–3 Here, we summarize the points pertinent to the present study. The SAC expansion for a totally symmetric singlet ground state is written as �SAC = exp�S��0� , �1� where �0� is a closed-shell Hartree-Fock �HF� single determi- nant and S = � I cISI † . �2� SI † is a symmetry-adapted excitation operator. The SAC co- efficient cI is calculated by solving the nonvariational equa- tions �0��H − ESAC���SAC� = 0 �3� and �0�SL�H − ESAC���SAC� = 0, �4� where H is the Hamiltonian and ESAC is the SAC energy. The SAC theory provides not only the SAC wave func- tion �SAC, but also a set of functions that constitute the basis for the excited states. Namely, the set of functions �K = PSK † ��SAC� �5� provides an adequate basis for expanding the excited states. Here, P is a projection operator P=1− ��SAC���SAC�. There- fore, we expand our excited states by linear combinations of �K as �SAC-CI �p� = � K dK �p��K, �6� which is the SAC-CI wave function for the excited states. Here, the index p labels the excited state. The SAC-CI wave function can also be defined for the excited state having dif- ferent symmetries and for the ionized and electron-attached states. Including these states, Eq. �5� is generalized as �K = PRK † ��SAC� , �7� where RK † denotes a set of excitation, ionization, and electron-attachment operators. The SAC-CI coefficients and energies are calculated by solving the nonvariational equa- tion �0�RL�H − ESAC-CI �p� ���SAC-CI �p� � = 0. �8� In the SAC/SAC-CI program, we can use the perturbation-selection method. We select only important ex- citation operators by the second-order perturbation theory and reduce the labors in the SAC and SAC-CI calculations.14,15,20 In the singles and doubles method �SAC/ SAC-CI SD-R�, we include all singles and apply a perturba- tion selection to doubles. In the ground state SAC calcula- tion, the doubly excitation operator SK † is included if 094105-2 R. Fukuda and H. Nakatsuji J. Chem. Phys. 128, 094105 �2008� Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp �Es�SK † ��� �g, �9� where Es�SK † � is the second-order energy contribution of SK † to the ground state, Es�SK † � = �0�SKH�0��0�HSK † �0� �0�SKHSK † �0� − �0�H�0� . �10� For excited, ionized, and electron-attached states, the doubly excitation operators which satisfy �Ep�RK † ��� �e �11� with Ep�RK † � = �0�RKH��Ref �p� ���Ref �p� �HRK † �0� ERef �p� − �0�RKHRK † �0� �12� are included in the SAC-CI calculations. Here, �Ref �p� is the reference wave function which has the energy ERef �p� . The rec- ommended usage of the perturbation-selection technique is found in the SAC-CI Guide.20 B. Conventional nondirect SAC/SAC-CI algorithm The conventional version of the SAC/SAC-CI program was formulated to be driven with the excitation operator la- bels. For example, the SAC-CI secular equation is written as � K �HLK + ULK�dK �p� = ESAC-CI �p� � K dK �p�SLK. �13� Here, HLK, ULK, and SLK denote linked, unlinked, and over- lap matrix elements, respectively, HLK = �0�RLHRK † �0� , �14� ULK = � I �0�RLHRK † SI †�0�cI, �15� and SLK = �0�RLRK † �0� . �16� For simplicity, we dropped higher-order unlinked terms in Eqs. �13� and �15� and the projection operator that will ap- pear in Eq. �16�. The details of this simplification were de- scribed in Ref. 20. The generalized eigenvalue problem of Eq. �13� is solved iteratively with the modified Davidson’s procedure,28 where the sigma vector, �L �m� = � K �HLK + ULK�bK �m� , �17� is constructed in each iteration step from the basis vector b and the matrix elements. In the conventional code, the matrix elements are evaluated with the projective reduction formula and are stored on disk. The loops for the integral evaluation �Eqs. �14�–�16�� and the sigma-vector construction �Eq. �17�� are driven through the excitation operator labels I, L, and K. We call this algorithm conventional “nondirect” algorithm, in contrast to the “direct” algorithm introduced in this paper. The important features of the conventional nondirect al- gorithm are as follows: �1� The perturbation-selection technique is easily intro- duced into the program. The perturbation selection re- duces the number of the matrix elements to be evalu- ated and the range of the summation in Eqs. �13� and �17�, thus dramatically reducing the total computation time. �2� The extensions to include general excitation operators are straightforward. This feature was important for the generalizations of the SAC-CI code, like the SAC-CI general-R method7 and SAC-CI for high-spin multiplicities.10 These extensions were done by ex- panding the nature of the excitation operators and could be easily implemented in the conventional code by ex- panding the loops of the operators labels I, L, and K. The matrix elements are calculated using the projective reduction formula, PROJR.27 The PROJR program evalu- ates the matrix elements by comparing the left and right configuration state functions. The algorithm repre- sented by Eqs. �14�–�17� does not depend on the spe- cific form of the R operators. Because of these attractive features, the conventional SAC-CI program has been extended to cover a wide range of chemistry and accuracy3 and has been successfully applied from fine spectroscopy21–23 of rather small molecules to biospectroscopy24 and photobiology involving moderately large molecules.25 However, the conventional algorithm has following demerits in comparison with the direct algorithm introduced in this paper. �1� The calculation is time consuming, particularly for the unlinked terms. The unlinked terms U in Eq. �15� have three indices of excitation operators, which result in a triply nested loop structure. However, many of the terms, �RL�H�RK † SI †�, are identically zero because of Slater’s rule, but the calculation of these unlinked terms can be done only at the innermost part of the loop. �2� A whole MO ERI has to be loaded on core memory. This is because the MO ERI is randomly accessed from the PROJR subroutine. There is no regulation in the ac- cess of the MO ERI. In constructing the direct SAC/SAC-CI code, we want to overcome these demerits, keeping the merits of the perturba- tion selection. C. Direct SAC/SAC-CI algorithm In the direct SAC method, the excitation operators S† and the coefficients c are defined by the MO labels instead of the excitation operator labels. For singles �S1� and doubles �S2�, they are S† = S1 + S2 = � i � a ci aSi a† + 1 2�ij �ab cij abSij ab† . �18� Hereafter, we use indices i , j ,k , . . . for occupied MOs, a ,b ,c , . . . for virtual MOs, and p ,q ,r ,s , . . . for general MO’s. Inserting Eq. �18� into Eqs. �3� and �4�, we obtain the expressions of the SAC energy, the SAC equations for S1 and S2. In the direct SAC algorithm, we first write down these working equations explicitly as summarized in Table X of the Appendix, where fqp= �p�f �q� and vqspr= �pq �rs�= �pr �qs� 094105-3 Direct algorithm for SAC-CI J. Chem. Phys. 128, 094105 �2008� Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp denote the Fock matrix element and the MO ERI, respec- tively. We introduced the permutation operator Pijab�¯�ijab = �¯�ijab + �¯� jiba �19� for convenience. Note that the expressions in Table X are not unique. We will introduce intermediate arrays to make the number of operations minimal.19 The expressions in Table X were designed to be effective in the perturbation-selection technique. The SAC-CI secular equation �Eq. �8�� is written in the nonsymmetric eigenvalue problems Hd�p� = ESAC−CI �p� Sd�p� �20� and H*d¯ �p� = ESAC−CI �p� S*d¯ �p�, �21� where d and d¯ denote the right-hand and left-hand eigenvec- tors, respectively. The superscript � denotes a Hermitian transpose. In the generalized Davidson procedure,28 we use m basis vectors bi and b¯ i, which are collected in matrices B and B¯ as B = �b1,b2, . . . ,bm�, B¯ = �b¯1,b¯2, . . . ,b¯m� . �22� Here, m is some small number used in generalized David- son’s procedure.28 The basis vectors bi and b¯ i satisfy the biorthonormal relation B¯ *SB = 1 , �23� where 1 is an m�m unit matrix. Then, we form a small Davidson matrix as H˜ = B¯ *HB = B¯ *� = �¯ *B . �24� The elements of �=HB and �¯ =H*B¯ are the so-called sigma vectors for the right-hand and left-hand basis vectors, respec- tively. Explicitly, they are �i = Hbi, �¯i = H*b¯ i. �25� To impose a biorthonormal relation to the basis, we define tau vectors as �i = Sbi, �¯i = S*b¯ i. �26� In each iteration, the small matrix C¯ *H˜ C = D �27� is diagonalized. After convergence, we obtain the SAC-CI vectors by the following expression: d = BC, d¯ = B¯ C¯ . �28� In the SAC-CI SD-R method, the sigma and tau vectors have three blocks: R0 �HF�, R1 �singles�, and R2 �doubles� as � = ��0�1 �2 � = � �0�ia �ij ab� . �29� The MO-indexed representations of the working equation are obtained from Eq. �8� by using the MO-indexed R operators. For singlet excited states, the R operators are R† = R0 + R1 + R2 = d0 + � i � a di aRi a† + 1 2�ij �ab dij abRij ab† . �30� The SAC-CI coefficients are obtained from the basis b as follows: dij ab = � m bij ab�m�C�m�. �31� The resultant working equations for the singlet excited states are summarized in Table XI of the Appendix. For the left- hand projection, we use the following notation for the Her- mite conjugation: ��¯ij ab�* = �¯ab ij , �b¯ij ab�* = b¯ab ij . �32� To simplify our representation, upper bars were dropped for the left-hand projections in Table XI. For triplet states, the MO-indexed R operator is R† = R1 + R2 = � i � a di aRi a† + � ij � a�b dij abRij ab† . �33� The working equations for the triplet states are summarized in Table XII of the Appendix. The equations for cation and anion doublet states are obtained from the triplet equations: Cation doublet states are obtained by replacing one of the unoccupied MO indices to infinitely separated orbital, e.g., Ri �† . Similarly, the electron-attached �anion doublet� states are written as an electron transfer from an infinitely sepa- rated orbital to one of the unoccupied orbitals like �→a. Thus, the R operators for the cation doublet and anion dou- blet are written as R† = � i di �Ri �† + � ij � b dij �bRij �b† �34� and R† = � a d� a R� a† + � i � a�b d�j abR�j ab† , �35� respectively. With the replacements to infinitely separated orbitals in Eqs. �34� and �35�, the working equations for the cation doublet and anion doublet are obtained. Note that the matrix elements including indices �, such as f i�, vij�a, etc., are zero. III. IMPLEMENTATION OF THE DIRECT SAC/SAC-CI ALGORITHM In the conventional nondirect program, the loops are driven with the labels, I, of excitation operators. The MO indices corresponding to the excitation operator are taken from the predefined labels of the excitation operators. This 094105-4 R. Fukuda and H. Nakatsuji J. Chem. Phys. 128, 094105 �2008� Author complimentary copy. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp correspondence may be written as SI †→Sijab† and is kept in the LABEL array, which is a 4�M matrix, where M is the number of the selected operators and each operator has four MO indices. Because th