Quadratically Convergent Algorithm for Orbital Optimization in the Orbital- Optimized Coupled-cluster Doubles Method and in Orbital-optimized Second- Order Møller-plesset Perturbation Theory Quadratically Convergent Algorithm for Orbital Optimization in the Orbital-optimized Coupled-cluster Doubles 

Using a Lagrangian-based approach, we present a more elegant derivation of the equations necessary for the variational optimization of the molecular orbitals (MOs) for the coupled-cluster doubles (CCD) method and second-order Møller-Plesset perturbation theory (MP2). These orbital-optimized theories are referred to as OO-CCD and OO-MP2 (or simply " OD " and " OMP2 " for short), respectively. We also present an improved algorithm for orbital optimization in these methods. Explicit equations for response density matrices, the MO gradient, and the MO Hessian are reported both in spin-orbital and closed-shell spin-adapted forms. The Newton-Raphson algorithm is used for the optimization procedure using the MO gradient and Hessian. Further, orbital stability analyses are also carried out at correlated levels. The OD and OMP2 approaches are compared with the standard MP2, CCD, CCSD, and CCSD(T) methods. All these methods are applied to H 2 O, three diatomics, and the O + 4 molecule. Results demonstrate that the CCSD and OD methods give nearly identical results for H 2 O and diatomics; however, in symmetry-breaking problems as exemplified by O + 4 , the OD method provides better results for vibrational frequencies. The OD method has further advantages over CCSD: its analytic gradients are easier to compute since there is no need to solve the coupled-perturbed equations for the orbital response, the computation of one-electron properties are easier because there is no response contribution to the particle density matrices, the variational optimized orbitals can be readily extended to allow inactive orbitals, it avoids spurious second-order poles in its response function, and its transition dipole moments are gauge invariant. The OMP2 has these same advantages over canonical MP2, making it promising for excited state properties via linear response theory. The quadratically convergent orbital-optimization procedure converges quickly for OMP2, and provides molecular properties that are somewhat different than those of MP2 for most of the test cases considered (although they are similar for H 2 O). Bond lengths are somewhat longer, and vibrational frequencies somewhat smaller, for OMP2 compared to MP2. In the difficult case of O + 4 , results for several vibrational frequencies are significantly improved in going from MP2 to OMP2.


I. INTRODUCTION
The Brueckner orbitals are defined to be a set of orbitals for which the single-excitation coefficients are zero in the full configuration interaction (CI) wave function.In 1953, these orbitals were introduced by Brueckner 1 in a self-consistent method for the study of nuclear matter.In 1958, Nesbet 2 first introduced Brueckner orbitals into quantum chemistry.[5][6][7][8][9][10][11][12][13][14] In a 1981 study, Chiles and Dykstra 6 utilized Brueckner orbitals in coupled-cluster (CC) theory for the first time.They used coupled-cluster singles and doubles method 15-17 a) Electronic mail: ubozkaya@ccqc.uga.edu.
In 1991 studies, Kobayashi et al. presented analytic gradients with respect to nuclear coordinates for the BD (Ref.11) and BD(T) (Ref.12) methods.][20] These studies demonstrated that the differences in the total energies and physical properties between BD and CCSD or BD (T) and CCSD (T) are generally small.]13 For open-shell systems, spin-restricted Brueckner orbitals have been investigated by Crawford et al. 23 As discussed by Sherrill et al., 13 an alternative way to obtain approximate Brueckner orbitals is to exclude single excitation determinants from the wave function and find the orbitals that minimize the energy.This procedure yields what may be called "variational Brueckner double orbitals."It has been postulated that the variational Brueckner orbitals and projective Brueckner orbitals are the same in the Full CI limit, 13 even if for truncated wave functions this equality does not exist.However, as shown by Köhn and Olsen, 24 this equivalence is not precisely obtained in the Full CI limit, with differences of tens of microhartrees in total energies for projective Brueckner orbitals and optimized orbitals for small molecules.Nevertheless, for truncated cluster operators, there remain some advantages to variationally optimized orbitals, especially if one stops at low orders instead of attempting to converge to the full CI limit.Once the orbitals are optimized, the wave function will obey the Hellmann-Feynman theorem for orbital rotation parameters.Therefore, there is no need for orbital response terms in the evaluation of analytic gradients.In other words, it is unnecessary to solve the first order coupled-perturbed coupled-cluster doubles equations.Further, computation of one-electron properties is easier because there are no response contributions to the particle density matrices.Another advantage is that the variational optimized orbitals can be readily extended to allow inactive orbitals.[27] Additionally, orbital-optimized coupled-cluster avoids spurious second-order poles in its response function, and its transition dipole moments are gauge invariant. 28,29 riational orbitals for the CCD wave function were first considered by Purvis and Bartlett. 15In a 1983 study, Purvis et al. 30 investigated the C 2v insertion pathway for BeH 2 and concluded that variational orbitals could be useful in CC studies of reactions where the dominant CI coefficients change along the reaction path.In a 1987 study, Scuseria and Schaefer 9 optimized MOs for the CCD and CCSD wave functions.They rotated the orbitals until the orbital Z-vector 31 became zero.This idea is based on the fact that the variationally optimized orbitals make no orbital response contribution to the CC gradient at convergence. 32Scuseria and Schaefer observed that, for the CCD method, optimized orbitals were obtained without any difficulty in their test cases.However, for the CCSD method they observed an erratic behavior of the Zvector and they concluded that the problem was related to the presence of e T1 in the CCSD model which already describes orbital relaxation effects.Scuseria and Schaefer further noted that the CCD energies obtained from the optimized orbitals were very close to the CCSD energies from the Hartree-Fock orbitals for their test cases.The difference in convergence properties of CCD and CCSD orbital optimization procedures was also discussed by Jankowski et al. 33 In 1998, Sherrill et al. 13 presented energies and analytic gradients for a CCD model using variationally optimized orbitals.Givens rotations were used for the unitary transformation of MOs, 34 and they computed amplitude response via the Z-vector approach of Handy and Schaefer, 31 which was first applied to CC theory by Adamowicz et al. 35 The orbital optimization equations are n this study, we will present a quadratically convergent algorithm for optimization of the MOs for the CCD method.First and second derivatives of the coupled-cluster doubles -functional 37,38 (CCD-) with respect to orbital rotation parameters are presented.A Newton-Raphson algorithm is used for optimization via the MO gradient and Hessian.Further, MO stability analyses are carried out using the MO Hessian at correlated levels (note that, by considering only the orbital-orbital portion of the Hessian, we neglect any orbital/amplitude coupling).Differences between our approach and the previous approach, 13 can be summarized as follows: (1) we minimize the CCD-functional with respect to orbital rotation parameters, while Ref. 13 considered amplitude derivative contributions via the Z-vector approach; (2) we use the exact MO Hessian, while Ref. 13 used an approximate diagonal MO Hessian; (3) in order to achieve quadratic convergence we use a full Newton-Raphson step, while Ref. 13 employed the DIIS technique; (4) we used an exponential transformation, while Ref. 13  used Givens rotations for rotating the orbitals.
Although in this initial work we focus on ground electronic states, orbital optimization of correlated wavefunctions can also be very useful for excited state studies.Equationof-motion coupled-cluster (EOM-CC) wavefunctions for doublet states are improved when using an OD reference, and similarly, orbital optimization can improve open-shell CIS(D) computations. 39e also consider orbital optimization in second-order Møller-Plesset perturbation theory (OMP2).Like OD compared to CCD, OMP2 has a number of theoretical advantages over standard MP2.OMP2 allows the computation of response properties (such as excitation energies, transition moments, and frequency-dependent polarizabilites) to be computed, which is not the case for canonical MP2.
Thus, response computations based on OMP2 may become an attractive alternative to methods such as secondorder approximate coupled cluster singles and doubles model (CC2) (Ref.40) and second-order polarization propagator approach (SOPPA). 41Again, we have minimized the MP2 -functional 38,42 (MP2-) with respect to orbital rotation parameters.The main difference of our new approach for MP2 from the earlier OO-MP2 approaches is that in previous studies [43][44][45][46] the Hylleraas functional (J 2 ) was minimized with respect to orbital rotations where J 2 is the Hylleraas functional, (1) is the first order correction to wave function, Ĥ0 is the unperturbed Hamiltonian, E (0) and E (1) are zeroth and first order energy corrections, and V is the perturbation operator.The Hylleraas functional J 2 provides an upper limit for the second-order energy, hence the minimization of J 2 yields a variational approach for finding the first-order wavefunction and second-order energy. 47The stationary points of the the Hylleraas functional J 2 coincide with the stationary points of the MP2-functional. 48More generally, it is possible to construct a Hylleraas functional J 2n , the minimum of which corresponds to the wavefunction (n) and the energy E (2n) .However, this method does not provide a variational approach for the energy corrections of odd order.For example, for both second-and third-order energy corrections one only needs to the first-order wavefunction (1) .By minimizing J 2 one can find the (1) that minimizes E (2) but not E (3) .In other words, the Hylleraas functional J 2 is an upper bound for E (2) , but it is not an upper bound for E (3) .Hence, the complete minimization for MP3 is not possible with that approach.Thus, the Hylleraas functional based approaches are not straightforward to apply higher orders of perturbation.However, our method is quite general and can be easily extended to higher orders of many-body perturbation theory.We will pursue such extensions in the near future.
We apply our algorithms to the H 2 O, C 2 , N 2 , F 2 , and O + 4 molecules and compare results for total energies and selected molecular properties with those from the standard MP2, CCD, CCSD, and CCSD(T) methods.Our C++ codes for OD and OMP2 are written by present authors for this research and interfaced with a preliminary version of PSI4 package. 49Our new CCD lambda (CCD-) and response density matrices codes are tested against the ACES II program. 50In order to validate the MO rotation formalism we wrote a separate quadratically convergent self-consistent field (QC-SCF) program and tested against the standard SCF code of PSI3 program, CSCF. 51The overall OD code is tested against QCHEM OD code. 52We have also reproduced the results of Sherrill et al. 13 and Krylov et al. 25 for optimized doubles (OD) and valence optimized doubles (VOD) using the DZP basis set.Further, we have verified our MO Hessian code comparing eigenvalues with previously published results at the SCF level. 53,54 THEORETICAL APPROACH A. Quadratically convergent orbital-optimized coupled-cluster doubles method (OD)

The CCD-energy functional and amplitude equations
In this section, for the MO indexing we will use conventional notation: i, j, k, l, m, n for occupied orbitals; a, b, c, d, e, f for virtual orbitals; and p, q, r, s, t, u, v, w for general spin orbitals.The Hamiltonian operator can be written using second-quantization formalism [55][56][57] as Ĥ = p,q h pq p † q + 1 4 p,q,r,s pq||rs p † q † ŝ r, where h pq is the one-electron Hamiltonian matrix element, pq||rs is the antisymmetrized two-electron integral, and p † and q are creation and annihilation operators.The coupled-cluster doubles (CCD) energy and amplitude equations can be summarized as follows: 56,58 where e T2 is the usual cluster double excitation operator, |0 is the reference determinant, t ab ij is the CCD t 2 -amplitude, ab ij | is doubly excited Slater determinant, E ccd is the CCD energy and E scf is the reference, self-consistent field (SCF), energy.
In order to obtain a variational energy functional ( E ccd ) it is convenient to introduce the following functional (CCD-): 38 where ˆ 2 is the double de-excitation operator and λ ab ij is a particular de-excitation amplitude.

The parametrization of OD wave function
The orbital variations may be expressed by means of an exponential unitary operator 48,[62][63][64][65] | = e K | = e K e T2 |0 , (13)   where K is the orbital rotation operator The effect of the orbital rotations on the MO coefficients can be written as where C (0) is the initial MO coefficient matrix and C(κ) is the new MO coefficient matrix as a function of κ.Now, let us define a variational energy functional (Lagrangian) as a function of κ, and first and second derivatives of the energy with respect to the κ parameter at κ = 0 Then the energy can be expanded up to second-order as follows: where w is the MO gradient vector, κ is the MO rotation vector, and A is the MO Hessian matrix.Therefore, minimizing the energy with respect to κ yields This final equation corresponds to the usual Newton-Raphson (NR) step.

Response (relaxed) density matrices
66 Therefore, it is appropriate to first introduce response density matrices for the CCD-functional.However, it should be noted that in some articles the definition of the response density matrices might be slightly different from ours.We have defined response density matrices as follows: where P± (pq) is defined as with P(pq) acts to permute the indices p and q.Using these definitions the one-and two-particle density matrices (PDMs) have the same permutational symmetries as the one-electron and antisymmetrized two-electron integrals, respectively.Furthermore, we can decompose the PDMs into reference and correlation contributions as follows: The energy of the CCD-functional can be expressed in terms of PDMs as follows: where E ccd is the correlation energy for the CCD-functional and f pq is the matrix element of the Fock operator.Correlation contributions for unique non-zero blocks of response one-and two-particle density matrices and some useful intermediates are presented in Appendix A.

Generalized-Fock and orbital (MO) gradient
Before writing the explicit equation for the MO gradient let us first define a generalized-Fock matrix (MO Lagrangian) as where F tu is a generalized-Fock matrix element.Then, the MO gradient can be written in terms of the elements of the generalized-Fock as follows: The MO gradient is defined by the asymmetricity of the generalized-Fock matrix and at convergence the generalized-Fock matrix will be symmetric.

Orbital (MO) Hessian
We can express the MO Hessian in terms of the generalized-Fock and response density matrices as follows: Intentionally, we keep our equations quite general in Eq. (37); 67 however, we need only to consider Hessian block(s) corresponding to non-redundant rotations (see below).One can show that for the CCD wave function with all orbitals active, only occupied-virtual (O-V) rotations are non-redundant. 48More generally, there may be inactive occupied (INO) and virtual (INV) orbitals (with frozen occupied and virtual approximations).In such a case, active occupiedinactive occupied (ACO-INO) and active virtual-inactive virtual (ACV-INV) rotations also should be considered in addition to all occupied-all virtual (O-V) rotations. 13nother important point is that, once the MO gradient vector is constructed, it can be seen that some rotations al-ready have zero gradients due to symmetry constraints.Further, some rotations may have numerically very small gradients, which are can be disregarded if the corresponding MO gradients are lower than a predefined cutoff value such as 10 −10 .Therefore, there is no need to construct Hessian or κ vector elements for those rotations.The Hessian and κ vector are only required for the independent-pairs (IDP) for which gradients are computationally non-zero.Hence, the number of computationally significant IDPs may be smaller than the number of non-redundant pairs (NRPs).We observed that the number of IDPs is significantly smaller than the number of total O-V type rotations (NRPs) in the case where all orbitals are active for test cases considered in this research.Constructing the MO Hessian explicitly only for the computationally significant independent-pairs, one can achieve both lower storage and quadratic convergence.
The computational cost for construction of the full MO Hesssian scales as n 6 , hence, one may consider solving Eq. ( 25) using iterative procedures with one index transformed quantities which scales as n 5 . 48However, constructing the MO Hessian only for computationally significant IDPs one can decrease the scaling to n 2 idp .n 2 , which is ∼ n 4 − n 5 for test cases considered in this research.Moreover, as OD and OMP2 iterations proceed, the number of numerically significant non-zero IDPs decreases, and after several iterations we observed that the number of IDPs is significantly reduced.Thus, it appears that our current algorithm is even more effcient than one-index transformation procedure.Further, in case of all orbitals active, one only needs to construct VOVO block of the MO Hessian, in case of frozen-core approximation one additionally needs VOOO and OOOO blocks.Hence, for the most cases where the CC methods are applicable one can store the MO Hessian in core memory since it has only two virtual index.Nevertheless, for large molecules, in order to avoid explicit construction and inversion of the MO Hessian matrix we also implemented the preconditioned conjugate gradient algorithm 68 for iterative solution of Eq. ( 25).

Updating MO coefficients
In each iteration, we need to update the MO-coefficients and transform one-and two-electron integrals from atomic orbital (AO) basis to MO basis where C (0) is the old MO-coefficients matrix, C is the new MO-coefficients matrix, and U is the MO-rotation matrix.In order to use Eq. ( 38) we can expand the exponential term up to the second order However, this approximate U matrix will not be unitary.Therefore, at each step we orthogonalize the U matrix with a modified Gram-Schmidt procedure. 68

MO optimization procedure
Our quadratically convergent orbital-optimized CCD wave function is defined by a set of orbital rotation parameters κ, and double excitation, and de-excitation amplitudes t 2 and λ 2 .In order to optimize the orbitals we need to construct the MO gradient via the generalized-Fock matrix and the MO Hessian.For computational efficiency, we optimize t 2 , λ 2 , and κ simultaneously, similar to the process of Sherrill et al. 13 We can summarize the MO optimization procedures as follows: (1) Perform the SCF procedure and store necessary parameters.
(2) Transform the AO integrals to the MO space.

B. The closed-shell (RHF reference) spin adapted equations
In this section, we present the spin-adapted equations for the closed-shell (RHF reference) case.The closed-shell spinadapted equations for the PDMs, generalized-Fock matrix, and MO Hessian are presented.One can find closed-shell spin adapted equations for t 2 amplitudes in the papers of Scuseria et al. 16,17 and for λ 2 amplitudes (as Z-vector) in the paper of Scheiner et al. 66

Spin-adapted equations for response density matrices
The spin integration for PDMs has been carried out diagrammatically following the prescription of Bartlett and Musial. 70For the closed-shell case, we have defined the spinadapted PDMs as follows: P QRS = 2 p α q α r α s α + p α q β r α s β , (43)   where uppercase indices indicate the spatial orbitals.For the one-particle density matrix (OPDM) we will use the same definition as in Eq. ( 26).However, the two-particle density matrix (TPDM) given in Eq. ( 27) loses some permutational symmetry properties with spin-adaptation.Therefore, in order to achieve the same (8-fold) permutational symmetry as the two-electron integrals in the spin-adapted formalism we have re-symmetrized TPDM.Our new definition can be written as Hereafter, we will drop the tilde notation and all two-particle density matrix expressions refer to the fully symmetrized TPDM, Eq. ( 44).Reference contributions for non-zero blocks of the PDMs are Spin-adapted equations for unique non-zero blocks of one-particle density matrix are while spin-adapted equations for unique non-zero blocks of two-particle density matrix are The intermediates G, V , V , and V are provided in Appendix B.

Spin-adapted equations for the generalized-Fock and MO gradient
For the closed-shell, case we have defined the spinadapted generalized-Fock and MO gradient as follows: Spin-adapted equations for the generalized-Fock and MO gradient are where h T P and P Q|RT are spin-free core Hamiltonian matrix element and two-electron integrals (h T P = h t α p α and P Q|RT = p α q α |r α t α ).

Spin-adapted equations for the MO Hessian
For the closed-shell case, we have defined the spinadapted MO Hessian as follows: The explicit spin-adapted equation for the MO Hessian is Now, Eq. ( 25) can be rewritten as We have kept the equations for the MO hessian, Eq. ( 60), and the generalized-Fock matrix, Eq. ( 57), quite general.However, for computational efficiency we have decomposed these equations in terms of unique non-zero blocks of the PDMs.Further, we have constructed only necessary block(s) of the MO Hessian.

C. Quadratically convergent orbital-optimized second-order Møller-Plesset perturbation theory (OMP2)
Because the MP2 method may be regarded as an approximate CCD method, we may use a similar approach for orbital optimization, although the detailed equations will of course differ.The MP2 energy and amplitude equations can be summarized as follows: ŴN = 1 4 p,q,r,s pq||rs { p † q † ŝ r}, where fN and ŴN are one-and two-electron components of the normal ordered Hamiltonian operator and subscript c means only connected diagrams are included.The MP2 amplitude equation can be written more explicitly as follows: 42,71,72 t ab (1) ij t ab (1) im f mj The {t ab (1) ij } are the first-order amplitudes (expansion coefficients for first-order correction to the wave function).In the regular RHF based MP2 method the last two terms of Eq. ( 66) may be set to zero because of the use of the canonical orbitals.However, since during the OMP2 iterations the Fock matrix will not be diagonal anymore, we add these offdiagonal Fock contributions to the amplitude equation.Alternatively, one may semicanonicalize the active Fock matrix at each iteration to get rid of those off-diagonal Fock contributions.
The MP2-correlation energy functional can be defined as 73 Correlation contributions for the unique non-zero blocks of PDMs are Second-order OPDM equations are similar to those for the CCD case.We can obtain the second-order G (2) ij and G (2) ab intermediates by replacing both t 2 -and λ 2 -amplitudes in Eq. (A4) and Eq.(A5) with first-order t 2 -amplitudes.In the usual canonical SCF case, diagonal parts of fN appear in zeroth-order and ŴN in first-order of perturbative energy corrections. 71,74 herefore, to obtain an overall second-order energy correction, (OO) and (VV) blocks of OPDM appear in second-order and the TPDM appears in first-order.Thus, we can write the correlation energy of the MP2-functional in terms of PDMs as follows: (2) pq f pq + p,q,r,s corr(1) pqrs pq||rs .(73)   We did not consider single excitations, although the Fock matrix will not be diagonal during OMP2 iterations, by analogy to the OD method.

III. RESULTS AND DISCUSSION
Results from the OD and OMP2 methods have been obtained for H 2 O, three diatomics (N 2 , F 2 , and C 2 ), and O + 4 for comparison with those from the standard MP2, CCD, CCSD, and CCSD(T) methods.The geometry and harmonic vibrational frequencies of the H 2 O molecule were determined using numerical differentiation of the total energies.For diatomic molecules, bond distances, vibrational frequencies, and spectroscopic constants (B e , D e , and α e ) are obtained using the DIATOMIC (Ref.75) program.For each method five single point energies, uniformly distributed about the equilibrium bond length (±0.005Å, ±0.010 Å), are provided to the DIATOMIC program.The CCSD, CCSD(T), and MP2 energies are obtained from PSI3 program, 51 while the CCD, OD, and OMP2 energies are obtained from our new codes.Pople's 6-311G(d) and 6-311G(d,p) basis sets, [76][77][78] Dunning's correlation-consistent polarized core-valence doubleand triple-ζ basis sets (cc-pCVDZ and cc-pCVTZ), 79,80 and Huzinaga-Dunning DZP and TZ2P basis sets [81][82][83][84] were used without the frozen core approximation for correlated procedures.For single point energy computations, 10 −10 was used as a convergence criterion for the total energies.
Table I presents comparison of total number of iterations for quadratically convergent (QC) and DIIS algorithms in the OD method.Total number of iterations are very similar for both algorithms.However, it should be noted that our new algorithm is more reliable.It is well-known that in problematic cases, such as the molecules which have flat orbital spaces, DIIS may fail to converge or may locate a undesired stationary point when initial guess orbitals are not enough close to the optimized orbitals. 85Especially in case of active space approximations, such as the valence active space, the initial MO Hessian may include negative eigenvalue(s).In such a case the DIIS algorithm may converge to a saddle-point rather than a minimum.However, our algorithm along with the augmented Hessian method guarantees convergence to a minimum if there is any. 69

A. H 2 O
Table II presents the total energies, equilibrium O-H bond distances (r OH ), H-O-H bond angles (θ H OH ), and harmonic vibrational frequencies (ω 1 , ω 2 , ω 3 ) for the H 2 O ( X 1 A 1 ) molecule with the six correlated methods.The total energies of CCSD and OD remain nearly identical due to the very similar equilibrium geometries.The energy difference between OD and CCSD is roughly 0.1-0.2millihartree.The OD and OMP2 methods provide total energies that are 1-2 milihartrees lower than those of CCD and MP2.In general as the basis set size increases the energy difference between orbital-optimized and unoptimized methods increases slightly.
The computed O-H bond distance and H-O-H bond angle are in reasonable agreement with the experimental values 86 of r OH = 0.9580 Å and θ H OH = 104.5 • especially with triple-ζ quality basis sets.The equilibrium O-H bond distance and H-O-H bond angle of the OD method are almost the same with the CCSD, and slightly longer (by 0.0007-0.0009Å), than those for the CCD.Similarly, the OMP2 method provides roughly 0.001-0.002Å longer bond lengths than MP2.
The predicted harmonic vibrational frequencies are in satisfactory agreement with the experimental harmonic values 87 of 3832 (ω 1 ), 1648 (ω 2 ), and 3943 (ω 1 ) cm −1 .The OD method provides almost identical harmonic vibrational frequencies with the CCSD method ( ω ≤ 2 cm −1 ), while yielding significantly lower frequencies than CCD.The MP2 and OMP2 methods provide slightly different frequencies; the difference between predicted ω 2 values is less than 4 cm −1 , while ω 1 and ω 3 values differ by as much as 33 cm −1 .Another salient feature appearing in Table II is that the orbital optimized methods yield slightly lower vibrational frequencies than corresponding unoptimized methods, since orbital optimization leads to longer bond distances.

B. N 2
Table III presents the total energies, equilibrium bond distances (r e ), harmonic (ω e ) and fundamental (anharmonic) (ν e ) vibrational frequencies, rotational constants (B e ), centrifugal distortion constants (D e ), and vibration-rotation coupling constants (α e ) for the N 2 molecule (X 1 + g ) with the six correlated methods.The OD and CCSD total energies are again close to each other.However, the energy differences, 1-2 milihartrees, are 10 times higher than in the H 2 O molecule, 0.1-0.2milihartrees.Further, the OMP2 method provides energies 5-8 milihartrees lower than MP2, while the OD method provides energies lower than CCD by 3-4 milihartrees.Thus for the N 2 molecule, the MP2 method is more sensitive to the orbitals than CCD, and the Hartree-Fock (HF) orbitals are far from being optimal choice for the MP2 method.Since the MP2 energies are already close to the CCSD energies, by optimizing the orbitals the OMP2 energies fall below the CCSD energies by 3-8 milihartrees.However, the OMP2 energies are still higher than those for CCSD(T) by 7-12 milihartrees.All the predicted bond distances are in satisfactory agreement with the experimental value 88 of 1.0977 Å.The OD method provides slightly shorter bond lengths than the CCSD method, whereas OMP2 provides longer bond lengths than MP2.The difference between OD and CCSD is r e = 0.0005-0.0007Å, while it is r e = 0.0068-0.0074Å between MP2 and OMP2.Further, the MP2, OMP2, and CCSD(T) methods predict longer bond distances than experiment with all basis sets considered in this research.The OMP2 method in particular overestimates bond distances with the double-ζ basis sets.Thus, our results indicate that in order to obtain an accurate bond distance for the N 2 molecule using the OMP2 method one needs to use large basis sets (at least triple-ζ quality basis sets).
For the harmonic vibrational frequency of the N 2 molecule, the MP2 and OMP2 methods provide significantly lower frequencies than other methods ( ω e ∼ 250 and ∼ 350 cm −1 , respectively), since these methods yield longer bond distances.The predictions of other methods are in reasonable agreement with the experiment. 88A similar feature is also observed for the fundamental frequencies.For spectroscopic constants, again the OD and CCSD methods yield very close results.Centrifugal distortion and vibration-rotation coupling constants are virtually the same for these methods.However, rotational constants differ by as much as 0.026 cm −1 because of the slightly different bond distances.The MP2 and OMP2 methods again provide significantly different results compared to the other four methods.The rotational constant is underestimated with these methods due to the overestimated bond distances.Another prominent feature of Table III is that the CCSD(T) method gives better results for spectroscopic constants than other five methods comparing to experiment.

C. F 2
Table IV presents total energies and spectroscopic constants for the F 2 (X 1 + g ) molecule with the six correlated methods.The OMP2 method provides lower total energies than the CCSD method by 3-7 milihartrees with triple-ζ basis sets.However, OMP2 energies are still higher than those of CCSD(T) by 10-14 milihartrees.
Except in the cc-pCVTZ basis, both OMP2 and CCSD(T) methods significantly overestimate the bond length.However, with the cc-pCVTZ basis these methods provide better results than the other methods.These results again suggest that one should use the OMP2 method with large basis sets (at least triple-ζ quality basis sets) in order to obtain reliable bond distances.Further, since OMP2 and CCSD(T) methods overestimate the bond length, they also give lower vibrational frequencies compared to other methods.The OMP2 and CCSD(T) methods underestimate the rotational constant due to overestimated bond distance except for the cc-pCVTZ basis.The difference between the OD and CCSD methods is B e ≤ 0.0015 cm −1 , while for OMP2 and MP2 it is B e ≤ 0.0318 cm −1 .For centrifugal distortion constant results, all methods are in reasonable agreement with the experimental value 88 of 3.3 × 10 −6 cm −1 .However, all methods underestimate the vibration-rotation coupling constant.MP2 and CCD methods provide smaller constants than other methods.

D. C 2
Table V presents the total energies and spectroscopic constants for the C 2 molecule (X 1 + g ) with the six correlated methods.The most remarkable feature of Table V is that the OMP2 total energies are even lower than those for CCSD(T) except for the cc-pCVTZ basis.Since MP2 energies are already lower than those for CCSD, by optimizing orbitals the OMP2 energies fall below the CCSD(T) by 2-3 milihartrees.
For the CCSD method T 1 -diagnostics 89-91 are 0.03, which are higher than the reference value of 0.02.Further, the percentage of reference determinants in the OD and OMP2 wavefunctions are approximately 89% and 88%, respectively.These results indicate a slight multireference character for the C 2 molecule.It is well-known that the C 2 molecule has lowlying excited states which strongly interact with the ground state. 92Hence, low total energies of OMP2 may be attributed to multireference character of the C 2 molecule, since the perturbation theory assumes that the exact wavefunction is dominated by a single determinant.Another interesting point is that with each basis set, the energy of the LUMO is negative.Further analyses with natural orbitals (NO) show that occupation numbers are 0.2424 and 0.2562 for the LUMO at OD/cc-pCVTZ and OMP2/cc-pCVTZ levels, respectively.The NO analysis show that the occupied character of the LUMO is more than for a regular virtual orbital.TABLE III.Predicted total energies (in hartrees), geometries (bond distances in angstroms), harmonic and anharmonic vibrational frequencies (in cm −1 ), and spectroscopic constants (B e and α e in cm −1 , D e in 10 6 × cm −1 ) for the N 2 (X Since the standard MP2 method already overestimates the bond distance, by optimizing the MOs the OMP2 method provides even longer bond distances.However, as the basis set size increases the OMP2 method gives a more reasonable bond distance.As observed for N 2 and F 2 , one should use the OMP2 method with large basis sets (larger basis sets than triple-ζ quality basis sets) in order to obtain a reliable bond distance.The predicted harmonic vibrational frequencies are in reasonable agreement with the experimental harmonic value 88 of 1828 cm −1 except for the OMP2 method which provides values around 1650 cm −1 .The low vibrational frequencies of the OMP2 method are arising from the overestimated bond distances.The difference between predicted frequencies with the OD and CCSD methods is ω e ≤ 13 cm −1 , while it is ω e ≤ 219 cm −1 with OMP2 and MP2.A similar situation is also observed for the fundamental frequencies.For spectroscopic constants of the C 2 molecule results of all methods except for OMP2 are in general agreement with experiment. 88The OMP2 method significantly underestimates the rotational constant due to the overestimated bond distances.However, for the vibration-rotation coupling constant the OMP2 method provides lower errors than other methods compared to experiment, differing by 0.0013 cm −1 or less from experiment.

E. O + 4
The O + 4 molecule is an important species in atmospheric ion chemistry. 93In a 1994 study, Lindh and Barnes 94 demonstrated that the standard wavefunctions, even including the complete active space SCF (CASSCF) method, suffer from an artifactual symmetry breaking problem.97] Spatial symmetry problems in reference (SCF) wavefunctions arise when the MO Hessian has a negative eigenvalue corresponding to a rotation which mixes the orbitals of different irreducible representations.The energy is a stationary point with respect to symmetry breaking rotations due to the symmetry constraints; however, that stationary point is not necessarily minimum.If the energy is a saddle-point with respect to symmetry breaking rotations, it is possible to find a lower energy symmetry broken solution.Such instabilities indicate the presence of different electronic states with close energies.Symmetry breaking problems manifest themselves when second-and higher-order molecular properties, such as vibrational frequencies, are computed. 48Symmetry breaking in the wavefunction leads to anomalous behavior in force constants, hence yielding spurious vibrational frequencies. 97,98 errill et al. 13 investigated the performance of the OD method for the O + 4 molecular ion.They showed that the OD method yields reliable vibrational frequencies as opposed to the standard CCSD method and discussed the symmetrybreaking problem in O + 4 in detail.Hence, we will not repeat that discussion.However, we are interested to compare the performance of the MP2 and OMP2 methods for this challenging system.Table VI presents total energies, geometries, and harmonic vibrational frequencies, for rectangular O + 4 ( 4 B 1g ) using the 6-31G(d) basis set.R oo denotes the bond length in each O 2 subunit, while R cc denotes the distance between the two parallel O 2 units.For the UHF MP2 and UHF CCD methods, geometry optimizations and frequency computations were carried out via analytic derivatives with the CFOUR program, 99 while UHF OMP2 com-putations were performed with our new codes via numerical derivatives.
Comparing to other methods both the OMP2 and MP2 methods, especially OMP2, yield a longer bond length for R oo , while yielding a shorter distance for R cc .The OMP2 method yields a longer bond length by 0.057-0.172Å for R oo , while it yields a shorter distance for R cc by 0.015-0.061Å than other methods.The MP2 method yield an anomalous vibrational frequency of 3540 cm −1 for ω 6 (b 3u ), which differs by a remarkable 2200 cm −1 from the experimental value of 1320 cm −1 . 100The OMP2 method gives an underestimated vibrational frequency of 769 cm −1 for ω 6 , which differs by 550 cm −1 from experiment.The low vibrational frequency of 769 cm −1 may be arising from the longer R oo value comparing to other methods.Although ω 6 of OMP2 is not in a good agreement with experiment, it is still much better than MP2.Further, the MP2 method gives a spuriously high vibrational frequency for ω 5 (b 2u ) due to the symmetry breaking.The OMP2 method yields a frequency of 144 cm −1 , which is in agreement with the OD and BD methods.
For ω 2 (a g ) and ω 4 (a u ) OMP2 predictions differ by only 1 and TABLE V. Predicted total energies (in hartrees), geometries (bond distances in Å), harmonic and anharmonic vibrational frequencies (in cm −1 ), and spectroscopic constants (B e and α e in cm −1 , D e in 10 6 × cm −1 ) for the C 2 (X

IV. CONCLUSIONS
We have presented a quadratically convergent algorithm for variational optimization of the molecular orbitals for the CCD and MP2 methods.We have reported explicit spinorbital and closed-shell spin-adapted equations for the response density matrices, generalized-Fock matrix (thus the MO gradient), and the MO Hessian explicitly.We have compared our OD and OMP2 approaches with the standard MP2, CCD, CCSD, and CCSD(T) methods.
For variational optimization of the molecular orbitals for OMP2 method, we have considered MP2 as an approximate CCD and we have used the same algorithm.Hence, we have constructed MP2-functional and minimized it with respect to orbital-rotation parameters.Comparing to earlier attempts to obtain variationally optimized MOs for MP2, [43][44][45][46] we have minimized MP2-functional instead of the Hylleraas functional.It should be pointed out that our method is more general than the Hylleraas functional based approaches and can be easily extended to higher orders of many-body perturbation theory.
When we compare the scaling of each method, CCD and CCSD scale as n 6 , while MP2 scales n 5 (due to integral transformation).However, scaling of OD is 3n 6 (t 2 -amplitude n 6 , λ 2 -amplitude n 6 , response PDMs n 6 , generalized-Fock n 5 , and MO hessian ∼ n 4 − n 5 ), whereas scaling of OMP2 is n 5 .Thus, for large molecules, a single point OD computation will be approximately 3-times more expensive than a CCSD computation.Similarly, the OMP2 method will be several times more expensive than the standard MP2 method.However, it should be noted that the OMP2 method generally converges in fewer iterations (6-8 iterations) than CCD, CCSD, and OD.For a geometry optimization procedure the computational times for CCSD and OD, or MP2 and OMP2, will be comparable.
For the molecules considered in this research we observed that the CCSD and OD methods give very close results.However, the OMP2 and MP2 methods provide the significantly different total energies.The OMP2 method yields 1-2 milihartrees lower energies than MP2 for the H 2 O molecule.For diatomics (N 2 , F 2 , and C 2 ) the energy differences between the OMP2 and MP2 methods are larger, 6-33 milihartrees.Further, for diatomics we have observed that OMP2 provides lower energies than CCSD.For the C 2 molecule, perhaps due to its small HOMO-LUMO gap, the OMP2 method yields even lower energies than CCSD (T).For each molecule considered in this research, the MO stability analyses have been carried out at OMP2/cc-pCVTZ and OD/cc-pCVTZ levels.The results of the MO stability analyses show that the OD and OMP2 methods properly converge to minima in orbital rotation space.
For geometrical parameters, the OMP2 method generally yields longer bond lengths than the MP2 and CCSD methods by 0.0207-0.0433Å.For the H 2 O and F 2 molecules the OMP2 method provides reasonable results for bond distances and frequencies, even better than CCSD, when compared to experiment.However, for the N 2 and C 2 molecules the OMP2 method overestimates bond lengths and underestimates vibrational frequencies.For the test cases considered, if the standard MP2 total energy is higher than that of CCD, then OMP2 yields reasonable predictions for bond distances and vibrational frequencies.In such cases, OMP2 predictions for bond lengths and frequencies are close to those from CCSD (T).However, if the MP2 total energy is lower than the CCD energy, then OMP2 considerably overestimates bond lengths and underestimates vibrational frequencies.Finally, the OMP2 method provides better vibrational frequencies than MP2 method in case of symmetry breaking problems (for O + 4 ).

ACKNOWLEDGMENTS
This research was supported by the U.S. National Science Foundation, Grants CHE-1054286 (to HFS) and CHE-1011360 (to CDS).U.B. would like to thank the Scientific and Technological Research Council of Turkey (TÜB İTAK) for supporting his scientific studies at the University of Georgia and Middle East Technical University.We would like to thank Dr. Andrew Simmonett for many helpful discussions.

V
MN = 2V i α j β m α n β − V i α j β n α m β , ABCD = 2V a α b β c α d β − V a α b β d α c β , B = 2V i α a β j α b β + V i α a β j β b α , (B10) V I AJ B = 2V i α a β j β b α + V i α a β j α b β , Sherrill et al. obtained new θ by scaling the current orbital rotation angles by a crude approximation to the diagonal elements of the orbital Hessian.They used a steepest-descent algorithm and accelerated convergence by employing Pulay's direct inversion of the iterative subspace (DIIS) method.
Check the convergence for the energy, root-meansquare (rms) MO gradient, maximum MO gradient, rms κ vector, maximum κ vector, rms t 2 , and rms λ 2 .(19) If convergence criteria are not satisfied for the seven parameters in Step 18 go back to Step 7.

TABLE I .
Comparison of total number of iterations for quadratically convergent (QC) and DIIS algorithms in the OD method.
1+ g ) molecule with the six correlated methods.

TABLE IV .
Predicted total energies (in hartrees), geometries (bond distances in angstroms), harmonic and anharmonic vibrational frequencies (in cm −1 ), and spectroscopic constants (B e and α e in cm −1 , D e in 10 6 × cm −1 ) for the F 2 (X 1 + g ) molecule with the six correlated methods.

TABLE VI .
Predicted total energies (in hartrees), geometries (bond distances in Å), and harmonic vibrational frequencies (in cm −1 ), for rectangular O + 4 ( 4 B 1g ) using the 6-31G(d) basis set.−1 from those of the OD and BD methods, while for ω 3 (b 1g ) the OMP2 result differs by 40 and 10 cm −1 from those of OD and BD, respectively.For ω 3 , the MP2 prediction of 609 cm −1 differs by 237 and 267 cm −1 from those of OD and BD.For ω 1 , both OMP2 and MP2 methods yield lower vibrational frequencies by 1108 and 562 cm −1 than OD, respectively, while they yield lower frequencies by 1113 and 567 cm −1 than BD.Overall, then, both MP2 and OMP2 provide significantly different vibrational frequencies than the presumably more reliable OD and BD frequencies.Although OMP2 underestimates ω 1 more severely than MP2 does, it provides significantly improved predictions for ω 3 , ω 5 , and ω 6 .
a This paper.b Sherrill et al. (Ref.13).c Barnes and Lindh (Ref.4).d Jacox and Thompson (Ref.100).e Spurious vibrational frequency due to symmetry breaking in the reference wave function.Downloaded 29 May 2013 to 130.207.50.154.This article is copyrighted as indicated in the abstract.Reuse of AIP content is subject to the terms at: http://jcp.aip.org/about/rights_and_permissions 2 cm