Symmetry-adapted formulation of the G-particle-hole hypervirial equation method

Highly accurate 2-body reduced density matrices of atoms and molecules have been directly determined without calculation of their wave functions with the use of the G-particle-hole hypervirial (GHV) equation method (Alcoba et al. in Int. J. Quantum Chem. 109:3178, 2009). Very recently, the computational efficiency of the GHV method has been significantly enhanced through the use of sum factorization and matrix-matrix multiplication (Alcoba et al. in Int. J. Quantum Chem 111:937, 2011). In this paper, a detailed analysis of the matrix contractions involved in GHV calculations is carried out. The analysis leads to a convenient strategy for exploiting point group symmetry, by which the computational efficiency of the GHV method is further improved. Implementation of the symmetry-adapted formulation of the method is reported. Computer timings and hardware requirements are illustrated for several representative chemical systems. Finally, the method is applied to the well-known challenging calculation of the torsional potential in ethylene.

The sparsity of these matrices due to symmetry can be exploited by carrying out a detailed analysis of the matrix contractions involved in the GHV method. This analysis, which is carried out here only for Abelian groups, leads to a symmetry-adapted formulation of the GHV algorithm which generates significant computational savings in both floating-point operations and memory storage.
The paper is organized as follows: In the next section the notation, definitions and general theoretical background of the GHV methodology are given. Section 3 is devoted to report the symmetry-adapted formulation of the GHV. The efficiency of this formulation and the results obtained in a set of applications of the method are reported and analyzed in Sect. 4. Finally, the conclusions of this work are given in the last section.

General notation
In what follows we will consider a many-electron system having a fixed and welldefined number of particles, N . We will also consider that the one-electron space is spanned by a finite basis set of 2K orthonormal spin-orbitals. Under these conditions the many-body Hamiltonian may be written in second quantization language [46] as: H = 1 2 p,q,r,s 0 H rs pq a p † a q † a s a r (1) where the operators a p † and a s refer to the fermion creator and annihilator operators respectively, and 0 H is a second-order matrix which collects the one-and two-electron integrals over the basis set. In this formalism a p-RDM, p D, corresponding to an N-electron state may be defined as: where pˆ is a p-electron density operator. The anticommuting rules satisfied by the fermion operators, together with the resolution of the identity operator, render possible the decomposition of a p-RDM element into a sum of terms involving lower-order RDM elements and additional terms describing q-body correlation effects [3,27,28,30,[40][41][42][43][50][51][52][53][54][55][56][57][58][59]. Let us consider the decomposition of the 2-RDM which provides the simplest example. Thus, it can easily be shown that The last term of this equation defines the elements of the second-order correlation matrix (CM), 2 C, or equivalently, of the G-particle-hole matrix, 2 G, which have the form: As can be seen, both matrices have the same elements but these elements occupy different positions in both matrices; hence, they are different matrices with different properties [56,59]. The corresponding operator definitions are therefore: These operators, which are very different from the 2-electron density operator, are at the center of the G-particle-hole hypervirial equation methodology.

The GHV methodology
By applying a matrix-contracting mapping involving the G-particle-hole operator to the matrix representation of the hypervirial of the N-electron density operator, one obtains the G-particle-hole hypervirial (GHV) equation [40][41][42][43][44][45], whose compact form is: As shown explicitly in ref. [40], using the second-quantized definition of the Hamiltonian in Eq. (6) and rearranging the creation and annihilation operators, we can express the GHV equation in terms of only third-and lower-order matrices: where (3;2,1) C is a third-order CM, which is interrelated to the 3-RDM as follows [60]: When is not a Hamiltonian eigenstate, the residual of the GHV equation, r.h.s of Eq. (7), does not generally vanish. Instead of looking for the wave-function that solves this equation, in the GHV method one looks directly for the corresponding CMs, or equivalently, and for computational convenience, for the RDMs. This can be achieved by solving the following system of differential equations for the 2-RDM in a transformation parameter λ [42]: where the operatorŜ, which depends on the residual of the GHV equation, is chosen at each λ to minimize the energy along its gradient Equation (9) only determines the changes in the 2-RDM with λ, and hence, without a knowledge of the 3-RDM, the differential equations are indeterminate. To remove this indeterminacy, the 3-RDM elements are approximated with a modified version of Nakatsuji-Yasuda's algorithm [19,42,44]. With this approximation algorithm, the system of differential equations can be propagated in λ until either (i) the least-squares error of the GHV equation, or (ii) the least-squares error of its contraction into the 1-electron space [40] ceases to decrease [42,44]. Computationally, this is achieved by integrating the differential equations with an adaptive variable step method due to Fehlberg [42,44,61].

Symmetry-adaptation of the GHV method
Evaluations of the GHV method have very recently been implemented using sum factorization and matrix-matrix multiplication at computational costs of K 6 in floating-point operations and K 4 in memory storage [42]. In this section it is shown that all intermediate arrays resulting from sum factorization can be defined by covariant equations. Therefore, these arrays retain symmetry properties of the ordinary density and electron integral matrices. The sparsity of all these (ordinary and intermediate) matrices due to symmetry can be exploited to further enhance the efficiency of the GHV method.

Analysis of the matrix contractions involved in the GHV equation method
As previously mentioned, the present version of our computational code implements Eq. (7) in the form which explicitly depends on the RDMs. When replacing in Eq. (7) the approximation algorithm for the 3-RDM elements, the structure of the terms showing the most time-consuming operations, which are of the type r,s, p is transformed into a sum of terms which are themselves formed by a triple or quadruple product of matrix elements [42] r,s, p (12) whereÂ is the antisymmetrizer operator, 1 D ( * ) and 1D( * ) are the 1-RDM and the first-order hole reduced density matrix (1-HRDM) corresponding to a Hartree-Fock reference calculation, and 2 is the second-order cumulant matrix, which is defined as: In order to render the construction time of the r.h.s. of this equation proportional to K 6 , each nest of loops appearing in the equation is sum factorized [42,[62][63][64]. As an example, let us just consider in detail the sum factorization of one of the elementary products occurring in this expression: This procedure allows one to perform the operation efficiently in two stages starting from the inner brackets outwards, with the auxiliary matrices defined by: A detailed analysis of the mathematical operations involved in calculation of the auxiliary matrices resulting from sum factorization reveals that all these intermediate matrices are defined by covariant equations. For instance, matrices 2 X and 2 R can be expressed in terms of elementary tensorial operations as follows: (1,3,2,4)→(4,1,3,2) con con ri ml where The covariance of these equations is very important since it implies that all the auxiliary matrices involved in GHV method retain symmetry properties of the ordinary density and electron integral matrices which can be exploited to further enhance the efficiency of the methodology. This is done in the following paragraphs.

Symmetry-adapted formulation of the GHV equation method
It is well-known that the operations in the symmetry group of a molecule, group F, leave the coefficients of the second-order matrix 0 H unchanged, and therefore 0 H is an invariant (2,2)-tensor for the group F [65]. Analogously, if the N -electron state is a non degenerate one (or more generally belongs to a 1-dimensional representation of F), then the p-RDMs are invariant ( p, p)-tensors for F [65]. Therefore, when the spin-orbitals are symmetry-adapted and ordered according to their irreducible representations, the 1-RDM is block diagonal [65]. Moreover, the 2-RDM and 0 H are sparse and when F is Abelian they are also block diagonal, after an appropriate reordering. As the auxiliary arrays resulting from sum factorization in the GHV methodology are defined in terms of covariant equations, they are also invariant tensors for the group F and has the same block structure as the density and electron integral matrices. The structure of the symmetry forbidden coefficients in all these first-and secondorder matrices is easier to analyze when the group F is Abelian, and hence only this kind of groups will be considered hereafter. When the studied electronic system has non-Abelian symmetry group, an Abelian subgroup will be considered.
In order to show how the sparsity of all the ordinary and intermediate matrices due to symmetry can be exploited to further enhance the efficiency of the GHV calculations, let us just discuss in detail how to apply the block structure of the tensors to perform the evaluation of the GHV operations for each of the auxiliary operations resulting from sum factorization of Eq. (14). The auxiliary matrix 2 X ri ml defined in Eq. (17) is a (2, 2)-tensor for the group F, and therefore the non null blocks are associated to irreducible representations π r , π i , π m , π l of F such that π r ⊗ π i = π m ⊗ π l . Hence, one could avoid the evaluation of the symmetry forbidden elements, and calculate the remaining elements as follows: 2 X ri ml = π s ,π p π r ⊗π s =π p ⊗π m π i ⊗π p =π s ⊗π l s∈π s , p∈π p 0 H rs In a similar way, the non null blocks of elements 2 R i j ml in Eq. (18) are associated to irreducible representations π i , π j , π m , π l of the group F such that π i ⊗π j = π m ⊗π l , and for each of these blocks one evaluates At this point, it is important to note that these constrained nested loops may be efficiently implemented with the help of auxiliary matrices that signal when the direct product of the irreducible representations of orbitals contains the totally symmetric representation of the molecular point group. Using these matrices it is possible both to avoid matrix reordering as well as to skip whole symmetry forbidden blocks of spin-orbitals in one single step.
The remaining matrix operations involved in the GHV methodology, which have been described in [42], can be analysed and evaluated in a similar way. Hence, if ordinary density and electron integral matrices entering in the GHV equation are invariant tensors for the group F, then at each iteration of the GHV methodology both the residual of the GHV equation and the intermediate auxiliary matrices used to evaluate it are invariant tensors for F. Therefore, it is possible to exploit their block structure to improve the efficiency of the GHV computations and reduce the memory requirements. In the next Section the computational advantages of this symmetry-adapted formulation of the GHV (sa-GHV) approach will be discussed and illustrated.

Efficiency of the sa-GHV equation method
While the exact number of matrix operations involved in the sa-GHV method is difficult to obtain in general, it is possible to develop a rough estimation. Thus, let us consider that the group F has f irreducible representations. If we assume that the partitioning of molecular spin-orbitals according to irreducible representation is strictly regular, then a straightforward calculation shows that (2, 2)-tensors have f blocks of size K 2 / f × K 2 / f , so they have K 4 / f non-null coefficients, and the operations in Eq. (22), which is one of the terms showing the most time-consuming operations, have a time proportional to f × K 2 / f 3 = K 6 / f 2 . These estimates show that the computational costs of the GHV method can be reduced by as much as a factor of f in storage and f 2 in floating-point operations. It must be noted that there are other kind of operations of lower complexity involved in the sa-GHV method, and therefore these relations are, as mentioned, rough estimations of the efficiency. In order to compare these estimates with those actually achieved, we have carried out a number of calculations on a series of selected small to medium sized systems in their ground states at equilibrium experimental geometries [66]. These systems, which have different point groups, have been chosen in order to explore the efficiency of the sa-GHV method and to compare the improvements implemented by the algorithm in each group. The PSI3 quantum chemistry package [67] has been used to calculate the one-and two-electron integrals, the orthonormal molecular orbitals and the initial values, at the HF level of approximation, of all the matrices required for initiating the iterative GHV process. In order to fairly assess the performance improvement due to symmetry, two sets of calculations have been carried out using the same algorithms. Thus, in one set of calculations we have assumed the largest Abelian subgroup of the point group describing the full symmetry of the system determined by PSI3 package and in another set the group assumed corresponds to C 1 symmetry. Consequently, the gains due to symmetry do not arise from spurious reasons but rather directly reflect the savings inherent in the symmetry-adapted method.
Tables 1 and 2 report the statistics pertaining to the computational cost and hardware requirements of GHV calculations for a number of molecules using three different basis sets. Instead of giving absolute values, which are strongly dependent on hardware facilities, the tables document the ratios of the computer time and memory requirements between the calculations performed in the largest Abelian subgroup of the point group describing the full symmetry of the system determined by PSI3 and those performed in C 1 symmetry. As can be appreciated from the documented data presented in Tables 1 and 2, the improvement increases not only with the order of the group but also with the size of the basis set considered, because the orbitals forming the larger basis sets are usually more evenly distributed among the representations of the group. The results show that computational efficiency ranges from 1.5 to 24.2 in floating-points operations rates and from 1.88 to 7.30 in memory allocation. These computed factors of reduction due to symmetry are indeed close to the theoretical estimates in most of the cases. This asymptotic limit is only actually achieved when the symmetry blocking of the orbitals is optimum. However, values of ∼ 0.3 f 2 in computer times and ∼ 0.9 f in memory are achieved in cases where the partitioning of molecular orbitals according to irreducible representation is far from regular, such as in 6-31G(d) acetylene which has 10, 1, 3, 3, 1, 10, 3 and 3 orbitals of a g , b 1g , b 2g , b 3g , a u , b 1u , b 2u , and b 3u symmetries, respectively.

Application of the sa-GHV equation method to the torsional potential in ethylene
Finally, as an application of the newly developed program, we have studied the challenging torsional potential in ethylene. This calculation, which is computationally very costly without the symmetry adaptation, is accurately carried out as will now be shown. While at its equilibrium (D 2h ) geometry ethylene is a well-behaved closed-shell molecule whose π -valence ground state can be described accurately by single-reference methods, it becomes a diradical at the barrier, when the π -bond is completely broken. Thus, at the twisted (C 2v at 90 • , and D 2 otherwise) geometry the ethylene's ground state is two-configurational. The tortional potential has been calculated by freezing all degrees of freedom except the torsional angle, and using experimental geometrical parameters [66]. Total energies and unoptimized barrier height calculated by the GHV method using a minimal (STO-3G) basis set are presented in Table 3. For comparison, we also report usual ab initio spin-restricted single-reference HF, MP2, CISD, and CCSD model results, and FCI ones, which is why a minimal basis set has been considered. Corresponding potential energy curves are shown in Fig. 1. As can be appreciated, HF, MP2 and CISD methods lead to severe unbalanced treatment of (π ) 2 and (π * ) 2 configurations, which results in unphysical shapes of the PES, i.e., a sharp cusp at 90 • and large errors in barrier heights. In contrast, the CCSD and GHV models perform well, producing very similar cuspless PESs which closely follow the FCI one, the GHV barrier height being slightly overestimated due to a degraded treatment of both the untwisted ethylene at 0 • and twisted ethylene at 90 • . Predictions of total energies and unoptimized barrier height calculated using a non-minimal (6-31G) basis set are also reported in Table 3. The results show that in this case HF, MP2, CISD, CCSD and GHV models lead to unphysical shapes of the PES when compared with the best available variational one (CISDTQQ within a frozen core approximation where inner core spin-orbitals have been assumed to be fully occupied in every configuration).  Fig. 1 Ethylene torsion, minimal basis set. In contrast to HF, MP2 and CISD, CCSD and GHV curves do not exhibit an unphysical cusp and closely follow the reference FCI curve In spite of this, CCSD and GHV still represents a definite quantitative advantage over both the HF, MP2 and CISD results.

Conclusions
In this paper, we have outlined a powerful scheme for including point group symmetry in GHV calculations. The method provides a means for exploiting sparsity in the density, integral and sum factorization intermediate matrices due to symmetry and is amenable to an efficient computational implementation. Highly symmetric molecules containing a few tens of atoms no longer represent a formidable computational obstacle since the cpu and memory requirements of calculations using this approach are not limited by the total number of spin-orbitals forming the basis set but rather by the maximum number of spin-orbitals belonging to the irreducible representations of the point group describing the full symmetry of the system. Finally, let us remark that the reported strategy for exploiting symmetry within the GHV method may also greatly accelerate other approaches for computing directly the 2-RDM such as the solution of the second-order contracted Schrödinger equation and related equations.