A Comparison of Three Approaches to the Reduced-Scaling Coupled

Article pubs.acs.org/JCTC

A Comparison of Three Approaches to the Reduced-Scaling Coupled Cluster Treatment of Non-Resonant Molecular Response Properties Harley R. McAlexander and T. Daniel Crawford* Department of Chemistry, Virginia Tech, Blacksburg, Virginia 24061, United States ABSTRACT: We have investigated the performance of the reduced-scaling coupled cluster method based on projected atomic orbitals (PAOs), pair natural orbitals (PNOs), and orbital specific virtuals (OSVs) for the prediction of linear response properties. These methods introduce different degrees of controllable sparsity in the ground-state and perturbed coupled cluster wave functions, leading to localization errors in properties such as dynamic polarizabilities and specific optical rotations. Using a series of chiral test compounds, we find that the inherent costs associated with computing response properties are significantly greater than those for determining the ground-state energy. As the dimensionality of the molecular system increasesfrom (pseudo)linear structures, such as fluoroalkanes, to cagelike structures, such as β-pinenethe crossover point between canonical-orbital and localized-orbital algorithms increases substantially. Furthermore, both the OSV and PNO methods provide greater reduction in cost (as measured by the size of the double-excitation space) than do PAOs, and PNOs provide the greatest level of sparsity for the systems examined here. Single-excitation truncation induces much larger errors than corresponding doubles truncation due to the fact that the first-order contribution to the one-electron perturbed wave function appears in the singles amplitudes. Both the PNO and OSV methods perform reasonably well for frequency-dependent polarizabilities provided appropriate thresholds are used for the occupation-number and weak-pair cutoffs on which each method depends. Specific rotations, however, are very sensitive to wave function truncation, to the extent that aggressive thresholds can yield the incorrect sign of the rotation, due to the delicate balance of positive and negative wave function contributions to the mixed electric-/magnetic-field response.

1. INTRODUCTION Theoretical chemistry is an essential tool of scientific investigation, standing as “a full partner with experiment” in discovery.1−4 Unfortunately, the most robust and reliable theoretical methods also come with a steep computational price, and accurate methods such as coupled cluster theory rapidly become impractical for use for large chemical systems. 5−11 The coupled cluster singles and doubles (CCSD) method, for example, scales as 6(N 6) (iterative), where N represents the size of the molecular system, whereas the most popular perturbational correction for connected triples, (T), scales as 6(N7) (noniterative). In the effort to reduce these high-degree polynomial scalings, a number of techniques have been developed to introduce greater sparsity in the coupled cluster wave function, including methods based on orbital localization, molecular fragments, orbital partitioning, singular value decomposition, and other techniques12−43 with varying degrees of success. Among these, one of the most widely used is the projectedatomic-orbital (PAO) approach pioneered by Pulay and Sæbø,12,44−46 which has been greatly advanced within coupled cluster theory by Werner, Schütz, and co-workers.14,16−20,47−49 For properties related to ground-state energies, e.g., reaction enthalpies and other thermodynamic constants for which the response or derivative of the wave function is not needed,14,16−20,45−48 such methods have proven to be © 2015 American Chemical Society

invaluable for treating large systems with little loss of accuracy compared to canonical approaches. Even for properties such as excitation energies, much progress has been made, and a number of reduced-scaling schemes have been introduced.50−56 However, the application of localization methods to higherorder response properties, e.g., (hyper)polarizabilities, magnetizabilities, and chiroptical response, has been much more limited.57−63 One of the first examinations of the performance of local coupled cluster methods for field-dependent response properties was reported by Korona and Werner58 where they considered dipole moments and static polarizabilities with the local coupled cluster singles and doubles (LCCSD) approach.58 Using the PAO recipe, they found that such electric-fielddependent properties were much more sensitive to orbital relaxation effects compared to untruncated calculations, but nonetheless, the average error for their test set remained small. They also showed that they could further reduce such errors by extending the virtual space for significant pairs beyond the standard domains that would otherwise have been sufficient for ground-state energy calculations. Russ and Crawford also applied the PAO approach to both static and dynamic dipole polarizabilities using coupled cluster theory.59,61 In their first work, they introduced a new orbitalReceived: September 18, 2015 Published: December 1, 2015 209

DOI: 10.1021/acs.jctc.5b00898 J. Chem. Theory Comput. 2016, 12, 209−222

Article

Journal of Chemical Theory and Computation

Herein, we report the first application of the PNO and OSV schemes to the calculation of coupled cluster nonresonant response properties,70−72 particularly dynamic polarizabilities and specific optical rotations. We choose as our benchmark systems a set of pseudolinear, two-dimensional, and threedimensional molecular structures. We consider the impact of various truncation schemes, including those limited only to the double excitation components of the wave function vs truncation of both singles and doubles, as well as exclusion of so-called weak pairs. The relative performance of each method is compared to both conventional coupled cluster methods and the PAO approach.

domain selection scheme based on the first-order response of the ground-state molecular orbitals to the external electric field as determined via the coupled-perturbed Hartree−Fock (CPHF) equations.59 (A similar approach was later utilized by Kats, Korona, and Schütz in the calculation of first-order properties and excitation energies using local CC methods.60) This approach naturally led to larger orbital domains than those required for ground-state energies in agreement with the findings of Korona and Werner. They later extended this work to dynamic response properties, including frequency-dependent polarizabilities and mixed electric-/magnetic-field-dependent properties such as optical rotations of chiral compounds.61 They found that such methods exhibited localization errors of only a few percent compared to the canonical approach for (pseudo)linear molecular structures. However, two- and threedimensional structures such as the cagelike compound (1R,4R)norbornenone required significantly larger domains, and thus much greater computational demands, for similarly accurate optical rotations.64 This work was later extended by McAlexander et al. to the optimized-orbital coupled cluster (OO−CC) method62 in an effort to develop a reduced-scaling coupled cluster model that was also gauge invariant for magnetic-field-dependent properties. One of the drawbacks of conventional (truncated) coupled cluster methods when applied to chiroptical properties, such as optical rotations, is that the computed property depends on the choice of length vs velocity representation of operators such as the electric dipole operator.65 As shown by Pedersen and Koch,66,67 variational optimization of the underlying molecular orbitals, in conjunction with large one-electron basis sets and/ or gauge-including atomic orbitals (GIAOs),68,69 can overcome this inadequacy. McAlexander and co-workers applied the PAO-LCC approach to the OO−CC formalism to determine if such optimization impacted the scaling of the method.62 To maintain a strongly sparse wave function structure, they relocalized the molecular orbitals in each iteration of the orbital optimization. They found that, in some instances, oscillations in the domain sizes caused the orbital optimization to fail, but in general, the accuracy of the local orbital-optimized coupled-cluster doubles (LOO−CCD) calculations were similar in accuracy to corresponding LCCSD results. Most recently, Friedrich et al. applied the incremental scheme to CCSD frequency-dependent dipole polarizabilities. Compared to the PAO approach, the incremental scheme13,29,30 employs a fragmentation-like, N-body structure to reduce the computational cost. At the second-order level, they achieved results with an average error of 0.6% for a 47molecule test set. Furthermore, similar to Crawford and Russ,61 they noted the increased error for complex, three-dimensional systems where the number of important interactions is larger than for simple (pseudo)linear structures. Although the PAO approach has been the most heavily applied reduced-scaling coupled cluster method, two other schemes have seen significant development in the last several years: the local pair natural orbital (PNO) approach resurrected by Neese and co-workers31,32 and the orbital-specific virtual (OSV) approach by Chan and co-workers.37,38 Given the poorer scaling exhibited by the PAO-based LCC approach as noted above for the response properties of three-dimensional structures, a natural question arises as to the efficacy of these newly emerging local correlation approaches for such properties.

2. THEORY 2.1. Response Theory. Coupled cluster linear response (CCLR) theory70,71 has proven to be an effective and useful tool for calculating various properties of small molecules. The linear response function may be derived as the second derivative of a time-averaged Lagrangian (quasi-energy) with respect to external perturbations A⃗ and B⃗ , leading to the expression ⎡ 1 ±ω ⟨⟨A⃗ ; B⃗ ⟩⟩ω = Ĉ P[̂ A( −ω), B( +ω)]⎢⟨ϕ0|(1 + Λ̂) ⎣ 2 ⎤ 1 B A B [A̅ , X̂ω ]|ϕ0⟩ + ⟨ϕ0|(1 + Λ̂)[[H̅ , X̂ −ω], X̂ω ]|ϕ0⟩⎥ ⎦ 2 (1)

where ϕ0 is the reference wave function, Λ̂ is the de-excitation operator associated with the left-hand wave function amplitudes, Ĉ is a symmetrizer with respect to both interchange of the sign of the field frequency and taking the complex conjugate of the expression, and P̂ symmetrizes with respect to the perturbation operators. Quantities with an overbar have been similarity transformed with the ground state coupled cluster exponential operator, e.g., for the electronic Hamiltonian, ̂

̂

(2) H̅ = e−T HeT ω To obtain the perturbed amplitudes, X̂ A , the linear systems of equations must be solved: ω ⟨ϕ|(H̅ − ω)XÂ |ϕ0⟩ = −⟨ϕ|A̅ |ϕ0⟩

(3)

In the case of frequency-dependent polarizabilities, where the only perturbation is the electric dipole operator, six sets of perturbed amplitudes (one for each Cartesian component of μ at both positive and negative values of ω) need be computed. For optical rotations, on the other hand, we must obtain perturbed wave function parameters for all Cartesian components of both the electric and magnetic dipole operators, μ and m, respectively, and both positive and negative field frequencies for a total of 12 sets of perturbed amplitudes. The linear response function then takes on the form ⎧ ±ω ⎡ ω ⟨⟨μ; m⟩⟩ω = Im⎨Ĉ P(̂ μ−ω , mω)⎢⟨ϕ0|(1 + Λ̂)[μ ̅ , X̂ m ] ⎣ ⎩ |ϕ0⟩ +

⎤⎫ 1 −ω ω ⟨ϕ0|(1 + Λ̂)[[H̅ , Xμ̂ ], X̂ m ]|ϕ0⟩⎥⎬ ⎦⎭ 2

(4)

Given the number of unperturbed and perturbed amplitudes that must be determined to obtain such response functions, local correlation approaches, such as PAOs, PNOs, and OSVs, have the potential to provide useful predictions of molecular 210


Article

Journal of Chemical Theory and Computation ⎡ fi (C)́ = min⎢⎣

properties for large chemical systems. However, the question remains as to whether such methods will maintain the inherent complexity of field-perturbed wave functions necessary for robust and accurate results while still providing sufficient sparsity to improve the computational costs of the methods. To answer these questions, we test three local correlation techniques applied to coupled cluster calculations of frequency-dependent polarizabilities and optical rotations. 2.2. Local Correlation Approaches. The rationale behind the reduced-scaling coupled cluster methods considered here is the localization of the occupied and virtual molecular orbital spaces, which leads to increased sparsity in the correlated wave function. For the occupied space, numerous methods are available and widely found to be effective, including Pipek− Mezey,73 Boys,74 natural bond orbitals,75 or others. However, the virtual space, which is often viewed as a byproduct of the self-consistent field procedure, is much more problematic, and the choice of virtual MO localization is the key distinguishing feature among the PAO, PNO, and OSV methods examined here. Although the PAO approach defines a single new virtual space that is available to all occupied orbitals/pairs, the OSV and PNO methods provide a unique virtual orbital space for each occupied orbital or occupied orbital pair, respectively. In the absence of truncation, the PAO space is the size of the oneelectron basis set, whereas the OSV space is Nocc × v and the PNO space is Npairs × v, where Nocc is the number of active occupied orbitals and Npairs is the number of unique pairs of active occupied orbitals. In practice, however, the virtual space is truncated (and linear dependencies removed) to achieve substantial computational savings. The orbital conventions used in this paper are as follows: μ, ν, λ,... pertain to atomic orbitals, i, j, k,... for localized, occupied orbitals, a, b, c,... for canonical virtual molecular orbitals, μ̅, ν,̅ λ̅,... for projected atomic orbitals, and finally ai̅ j, b̅ij, ci̅ j,... for virtual orbitals in the PNO or OSV basis where the subscript denotes the pair/orbital-specific basis to which the orbital belongs. Additionally, quantities with an overbar are in the PNO/OSV basis, whereas quantities with a hat are in the pair/ orbital-specific semicanonical basis. 2.2.1. Projected Atomic Orbitals. In the PAO approach, the occupied orbitals are subjected to a conventional localization scheme, such as the Boys74 or Pipek−Mezey73 procedures, after the initial Hartree−Fock calculation. The virtual space is then transformed to a localized basis by projecting away the contributions from the occupied space, such that ⎛ |χμ ⟩ = ⎜⎜1 − ̅ ⎝

i

⎠

μ ∈ [i] ν

(6)

where χ́i, with the corresponding MO coefficients Ć iμ, is the approximation to χi constructed from the current set of domain-included AOs. If the chosen subset of atoms (and thus their associated set of AOs) does not produce a value of f i within a user-defined cutoff (0.02 is a typical value), then additional atoms are added to the domain, normally in a ranked order based, for example, on Mulliken orbital populations. Because each PAO is mapped to an original AO, the PAOs associated with these atoms then comprise the domain of that particular occupied orbital

χμ ∈ A ∈ χi

(7)

̅

Alternate domain schemes have been designed to take into account magnetic and/or electric perturbations, e.g., through the coupled-perturbed Hartree−Fock (CPHF) equations,59,61 where the HF property tensor contributions are decomposed on an atom-by-atom basis in a manner similar to the original domain scheme described above. For electric perturbations, for example, i

αxy =

i

vir

∑ ∑ Uaixμaiy = occ

a

AO

∑ ∑ U ρxiμρyi occ

ρ

(8)

where αxy is a component of the dipole polarizability, Uxai is a solution to the CPHF equations for a dipole perturbation, μyai denote electric dipole integrals, and the final equality is obtained by back-transformation of the virtual orbitals to the AO basis. Partitioning this AO summation over atomic contributions yields, for a single localized occupied orbital, AO

αxyiA =

∑ U ρxiμρyi

(9)

ρ∈A

and analogous expressions may be derived for other perturbations. In order to avoid the cancellation of terms arising from positive and negative contributions to the sum, Russ and Crawford proposed59,61 a completeness check analogous to that of Boughton and Pulay, viz. AO

εxyi =

AO

∑ |U ρxiμρyi | − ∑ ρ

⎞

∑ |χi ⟩⟨χi |⎟⎟|χμ ⟩

∫ (χi − χí )2 dτ ⎤⎥⎦ = 1 − ∑ ∑ Cμ́ i SμνCνi

ρ ∈ [i]

|U ρxiμρyi |

(10)

2.2.2. Pair-Natural Orbitals and Orbital-Specific Virtuals. Following localization of the occupied orbitals via a standard approach, the virtual space in the PNO and OSV approaches is defined based on natural orbitals derived from a low level model of electron correlation, normally second-order Møller− Plesset (MP2) perturbation theory. Given the amplitudes of the first-order wave function

(5)

where |χμ̅⟩ is a PAO, |χμ⟩ is an atomic orbital (basis function), and |χi⟩ is an occupied MO. Thus, although the PAOs form a nonorthogonal set, they remain orthogonal to the occupied space. The reduction of the size of the correlated wave function in the PAO approach arises from the construction of occupied orbital domains such that the virtual space for a given occupied orbital domain (denoted [i], for example) is a certain subset of the full list of PAOs. The selection of the appropriate subset is typically based on an atom-by-atom scheme in which all the (nucleus-centered) basis functions associated with a given atom are used to compute the value of a completeness function76 for the occupied orbital, e.g.,

ij Tab =

⟨ab|ij⟩ fii + f jj − ϵa − ϵb

(11)

where ⟨ab|ij⟩ is a two-electron integral in Dirac’s notation, f ii and f jj are diagonal elements of the occupied block of the Fock matrix, and ϵa and ϵb are (canonical) virtual orbital energies. A pair specific (each ij) or orbital specific (each ii) virtual density may be constructed, such that 211


Article

Journal of Chemical Theory and Computation Dij =

2 (TijT̃ ij † + Tij †T̃ ij) 1 + δij

carried out a series of calculations of dipole polarizabilities and specific optical rotations for a set of hydrogen helices, a sequence of fluoroalkanes, and the paradigmatic chiral species (S)-1-phenylethanol, (1R,4R)-norbornenone, and (1R,5R)-βpinene. For selected systems, we also calculated results for the PAO approach for comparison. For the latter three molecules, we employed the CC2 method,78 whereas for the hydrogen helices and fluoroalkanes, we used CCSD. For the polarizability and optical rotation calculations, the aug-cc-pVDZ basis set79,80 was used. Given that the development of production-level implementations of local correlation methods requires a considerable time investment, we have chosen to carry out simulations of the effects of the PAO, PNO, and OSV truncations on molecular response properties using our PSI481 canonical-MO coupled cluster code as a foundation. This approach, which has been used in previous investigations of PAO-based codes,14,45,50,61,62 first requires transformation of the coupled cluster amplitude residuals Rij from the canonical MO basis to the localized basis in each iteration

(12)

where

T̃ ij = 2Tij − Tij †

(13)

Diagonalization of the density yields the PNOs or OSVs for a given pair or orbital, Qij, with associated occupation numbers nij, such that

Dij Qij = Qijn ij

(14)

Orbitals with occupation numbers below a user-selected threshold are neglected to introduce truncation, such that ij Tab ≈

∑

ij ij Q aa Taij b ̅ Q bb ̅

aijbij ∈ [ij]

ij̅

ij̅ ij

ij

(15)

where ai̅ j ∈ [ij], such that > the cutoff. In addition, for the OSV approach, the nondiagonal transformation matrix is formed by a direct sum of the transformation vectors for pairs ii and jj. nijaij̅

Qij ≡ Qii ⊕ Q jj

R̅ ij = Qij †RijQij

(16)

Next, the residuals are transformed to a semicanonical pair- or orbital-specific basis where the virtual block of the Fock matrix is diagonal

Unlike the PAO and OSV approaches, PNOs are orthogonal among themselves for a given pair, but both OSVs and PNOs are nonorthogonal between pairs, yielding an overlap

Sij ,kl = Qij †Qkl

(18)

ij R̂ = Lij †R̅ ijLij

(17)

(19)

Diagonalizing the virtual Fock block in the local basis provides a pair-/orbital-specific semicanonical basis

The computational costs of solving the CC linear-response equations within the PNO and OSV are similar to those for the ground-state amplitude equations. A significant difference, however, arises in the CC perturbed wave function equations, eq 3, which are normally formulated in terms of the similaritytransformed Hamiltonian, H̅ . When density-fitting is used to represent the two-electron part of the bare electronic Hamiltonian, as is typically done to reduce the costs of the PNO and OSV methods, then H̅ must be (re)constructed in each iteration of the solution of eq 3 (see ref 77 for a related discussion in the context of the equation-of-motion CC method for excited states). 2.2.3. Pair Truncations. To approach linear scaling in practical calculations, truncating only the virtual space is not sufficient, and additional approximations must be introduced. One of the most common techniques is to categorize the various occupied pairs such that only a subset are treated in the full manner, whereas the other set or sets are neglected or treated in a less costly way. In PAO calculations, for example, the pair rankings are often straightforwardly based on the number of common atoms between the domains of any two occupied orbitals.14 For instance, occupied pairs with no shared atoms are typically labeled “weak” and approximated at a lower level, such as MP2. In addition, the pair hierarchy may be defined using distance cutoffs. For the OSV and PNO approaches, the virtual space for each occupied pair is distinct, and an alternative pair truncation scheme must be employed. The typical approach for truncating pairs in these methods is based on the approximate MP2 pair energies,31 such that when the pair energy falls below a threshold, it may be dropped from the calculation or treated in a more approximate manner.

ij

F̅ ijV Lij = SijLijF̂ V

(20)

F̂ijV

ij

where S is the overlap matrix in the given pair space and is a diagonal matrix with entries ϵ̂aij, ϵ̂ijb, ... Next, the energy denominators are applied to obtain the amplitude increments, where the occupied energies are the diagonal elements of the occupied block of the Fock matrix, and the virtual energies are in the specific pair/orbital semicanonical basis, such that ij Δ̂ab =

ij

R̂ab fii + f jj − ϵ̂ija − ϵ̂bij

(21)

The update in the semicanonical basis may then be backtransformed to the canonical virtual basis ij

Δij = QijLijΔ̂ Lij †Qij †

(22)

and the new amplitudes computed in the usual manner Tijnew = Tijold + Δij

(23)

In an analogous fashion, the local filter technique may be applied to singles using the diagonal pairs for the corresponding transformations i

Δi = QiiLiiΔ̂ Lii †Qii †

(24)

and tinew = tiold + Δi

(25)

This same filtration technique is applied to the Λ amplitudes as well as to each set of perturbed wave functions, XAω. Although no computational savings are yet achieved with this simulation approach, it allows for more rapid testing of the local correlation techniques to determine if further implementation

3. COMPUTATIONAL DETAILS To test the performance of the OSV and PNO local correlation methods with coupled cluster linear response theory, we have 212


Article


computed using eq 3 with a static electric-field perturbation oriented along the A (long) axis of the molecule. The occupied orbital space was localized according to the Pipek−Mezey recipe,73 and the virtual space is that defined above for each local correlation method but without truncation so that the level of sparsity introduced by each approach can be compared. For all three methods, we observe that the perturbed amplitudes are less sparse than their unperturbed counterparts, which is a result of the fact that the perturbed wave function corresponds to the derivative of the ground-state wave function with respect to the external field, as opposed to resulting from the application of a finite field.64 For the PAO and OSV bases, the peak positions for the perturbed amplitude distributions are shifted to roughly an order of magnitude more positive than their unperturbed counterparts. For instance, the peaks for the OSV curves are around 10−5 and 10−4 for the unperturbed and perturbed amplitudes, respectively. However, despite the similar degree of sparsity, the OSV distribution is much narrower than the PAO distribution and also exhibits an asymmetric shape, presumably due to the requisite elimination of linear dependencies in the construction of the pairwise virtual space (eq 16). For the PNO distributions, the peak shift between unperturbed and perturbed amplitudes is larger than for the OSV or PAO basis with the perturbed amplitudes approximately 2 orders of magnitude less sparse than the unperturbed amplitudes (versus 1 order of magnitude for OSVs or PNOs). Furthermore, both PNO distributions are shifted to the left relative to their OSV and PAO analogues, indicating a greater degree sparsity in the former. These distributions indicate that, for all three methods, the cost of computed response properties that depend on the first-order perturbed wave function will necessarily be higher than for properties that require only the electron density. However, the PNO basis could potentially provide a more efficient means for determination of such properties than the PAO or OSV approaches, provided that the inherent sparsity of the wave function can be effectively utilized. 4.2. (P)-(H2)n Helices. For the initial performance testing, we selected a simple and straightforward benchmark, (P)-(H2)n helices. The structure of each molecule was chosen to have a H−H bond distance of 0.75 Å and a helical angle of 60° with a distance between H atoms on the helical axis of 1.5 Å. These hypothetical systems have been used previously in probing local correlation techniques,61 but only the PAO approach has been considered. The weak, nonbonding interactions between hydrogen dimers result in the chains’ strongly localized electronic structure, which should thus be a “best case” system for local correlation techniques. In addition, the presence of a stereogenic axis and low molecular weight combine to yield relatively large specific rotations, whose sensitivity to truncation of the space should be easily determined. Figure 2 reports the dynamic polarizability and specific rotation (both at 589 nm) as well as the T̂ 2 ratio for the smallest [(H2)2] and largest [(H2)7] of the helices examined here, computed at the CCSD/aug-cc-pVDZ level of theory. For the polarizability, the OSV and PNO values consistently approach the canonical value from below and do not oscillate or cross it as we tighten the truncation threshold. That is, as the size of the virtual space increases, new contributions to the polarizability are always positive. This can be understood simply from examination of the leading terms of eqs 1 and 3 using the electric dipole operator for both perturbations A⃗ and B⃗ ; the largest contribution arises from the term in eq 1 that is

is worthwhile. In addition, the simulation code may be used for verification of subsequent production-level programs. A similar approach has been used in previous investigations.14,45,50,61,62 We should note that identical truncations and/or domainselection schemes must be applied both to the unperturbed and perturbed wave functions for the response functions we obtain to match the corresponding energy derivatives obtained, for example, using finite-difference (finite-field) techniques. Furthermore, to assess the extent of truncation in the local correlation calculations, we compare the T̂ 2 ratios between methods, defined simply as the ratio of the number of doubles amplitudes in the truncated and canonical bases. For the OSV and PNO approaches, several occupation thresholds were tested. We used both the Boughton−Pulay domains76 and coupled-perturbed Hartree−Fock-based domains61 for the PAO results. All property values are calculated using the modified-velocity gauge82 representation of the electric dipole operator. Additionally, three different types of truncations were examined. We explored the effects of truncating only the doubles (D) amplitudes, doubles and singles (D+S), and doubles with neglecting weak pairs (D +WP), where occupied pairs with an MP2 pair energy less than 10−4 Eh were dropped. For the sparsity tests (vide infra), the 631+G* split-valence basis set was selected.83 Occupied MOs corresponding to core orbitals were kept frozen in all calculations.

4. RESULTS AND DISCUSSION 4.1. Sparsity: Simple Alkanes. In exploring and comparing the local correlation techniques, we first examine the relative sparsity of the coupled cluster wave function when represented in the virtual-MO space defined by each method. For this test, we have focused on a series of n-alkanes in a manner similar to that reported earlier by one of the present authors for the PAO approach.64 Figure 1 presents the distribution of unperturbed and perturbed CCSD/6-31+G* double-excitation (T̂ 2) amplitudes for the largest alkane considered here, n-octane. The perturbed amplitudes were

Figure 1. Statistical distribution of CCSD/6-31+G* unperturbed and perturbed double-excitation (T̂ 2) amplitudes for n-octane. The perturbation corresponds to a static electric-field oriented along the long axis of the carbon backbone, as determined using eq 3. 213


Article


Figure 2. Sodium D-line dynamic polarizabilities (in au), specific rotations [in deg dm−1 (g/mL)−1], and T̂ 2 ratios for (H2)2 and (H2)7 helices as a function of −log (occupation threshold) for the OSV (red) and PNO (blue) methods. The canonical (untruncated) values are shown in yellow.

the negative of the linear response function in eq 1, each new contribution is necessarily positive as the size of the virtual space grows. Optical rotations, on the other hand, converge in an oscillatory manner toward the canonical result as we approach the full virtual space. For (H2)2, the OSV and PNO curves cross

linear in the perturbed amplitudes, the leading contribution to which is the electric dipole operator between the reference wave function, |ϕ0⟩, and singly excited determinants. Thus, the largest contribution to the polarizability arises from the square of the dipole integrals, (μ⃗ ia)2, divided by (negative) MO energy differences. Recognizing that the polarizability corresponds to 214


Article

215

0.0180 0.00594 0.00264 −49.1 −41.0 −32.6 116.6 186.6 257.4 0.218 0.105 0.0628 −55.4 −46.1 −37.0 120.2 192.7 266.4

T̂ 2 ratio 10

[α]D α T̂ 2 ratio

10

[α]D α

0.107 0.0497 0.0284 118.6 190.0 262.4 121.9 195.8 271.0 1-fluoropropane 1-fluoropentane 1-fluoroheptane

−60.8 −50.2 −40.5

−53.6 −45.2 −36.5

T̂ 2 ratio 10

[α]D α [α]D

canonical

α molecule

Occupation number thresholds of 10−6 and 10−7 were used for the OSV and PNO approaches, and Boughton−Pulay (BP) domains with a cutoff of 0.02 were used for the PAO approach with only double-excitation amplitude truncation. a

0.182 0.0901 0.0534 −44.3 −39.8 −32.4 120.0 192.3 265.7 0.0559 0.0198 0.00910 −55.3 −44.5 −35.7 119.1 190.8 263.7

[α]D

BP 10

α

[α]D

T̂ 2 ratio

α

PAO

−7

PNO −6 −7

OSV −6

Table 1. Dynamic Polarizabilities (au) and Specific Rotations [deg dm−1 (g/mL)−1] computed at 589 nm at the CCSD/aug-cc-pVDZ Level of Theory for a Series of (M)Fluoroalkanesa

the canonical result several times, but the curves for the longer helix (H2)7 only cross twice. This behavior may be understood using the same analysis as above for the polarizabilities; however, in the case of the mixed electric-dipole/magneticdipole polarizability required for the specific rotation, the necessary products of electric- and magnetic-dipole integrals, m⃗ ia and μ⃗ ia, may produce both positive and negative contributions. As a result, fortuitous error cancellation arises for these ideal systems at looser thresholds. For (H2)2, the rotations are more accurate at 10−6 than the 10−8 threshold, and a similar scenario is seen for (H2)7 for the 10−4 and 10−5 cutoffs. The localization errors in the polarizabilities are relatively small, around 1%, for both OSV and PNO methods, particularly for a truncation threshold of 10−6 and smaller (tighter than is typically required for energy-only calculations). Corresponding errors for optical rotations are larger and less well-behaved due to the oscillatory behavior described above. For (H2)7, for example, a threshold of 10−6 yields a roughly 10% error for the OSV approach and a 20% error for PNOs. Tightening the threshold to 10−7 decreases these errors to 5 and 0.5%, respectively. However, further constraining the threshold to 10−8 decreases the OSV error to 2% but increases the PNO error to 3%. Of course, the T̂ 2 ratio for each method (which is one measure of the cost of a production-level implementation of the model) increases as the threshold is tightened, though the rise is slower for the PNO method than for OSVs with the former requiring approximately 5% of the total wave function at a cutoff of 10−7 vs the latter, which requires ∼20% of the wave function at the same threshold. 4.3. (M)-1-Fluoroalkanes. The linear, chainlike fluoroalkanes serve as a next step in testing local correlation methods. Similar to the hydrogen dimer helices, these molecules have a pseudo-one-dimensional structure but are larger with a more complex electronic structure with stronger bonding interactions (the B3LYP/6-31G**optimized geometry of each molecule was used for all computations in this work). As noted in previous studies,61,62 although these systems have zero optical rotation due to conformational averaging, we select a single, optically active conformer with a negative F−C1-C2−C3 dihedral angle and thus having (M) axial chirality. Unlike earlier works featuring this test set, we limit our investigation to only the aug-cc-pVDZ basis set. Although basis set effects are considerable in response calculations, diffuse functions have been shown to be essential for accurate results.84,85 Also, for fluoropentane, the 6-31+G* and aug-cc-pVDZ basis sets yielded similar localization errors in a previous study.62 Table 1 reports dynamic polarizabilities and optical rotations at 589 nm for the three fluoroalkanes. The polarizability increases with the size of the molecule, as expected for a sizeextensive property. However, unlike the hydrogen-molecule helices, the specific rotation decreases as the chain length grows in this case because the extension of the carbon backbone does not increase the optically active region of the molecule. For the polarizabilities, the relative error for each approach increases only slightly with increasing chain length: 3% for the OSV with a 10−6 cutoff, less than 2% for OSV with 10−7 cutoff, 5% for PNO (10−6), less than 3% for PNO (10−7), and less than 2% for PAO. Although the largest error arises for fluoroheptane with the PNO method with an occupation threshold of 10−6, the fraction of local T̂ 2 amplitudes relative to canonical amplitudes is quite small in this case at around 0.003 (i.e., only 0.3% of the doubles amplitudes are required in the local

T̂ 2 ratio



Article

216

258.7 303.7 320.0 324.1 0.0288 0.160 0.467 0.593 −158.0 −182.4 −192.7 −195.4 304.1 316.2 319.6 320.2 0.0442 0.236 0.728 0.974

[α]D

307.0 320.5 324.4 325.0 10−6 10−8 10−10 10−12

−128.5 −150.1 −160.9 −164.3

α T2 ratio D

[α]D α cutoff

For comparison, the canonical polarizability is 325.1 au, and the canonical specific rotation is −165.2 deg dm−1 (g/mL)−1. Truncation of the localized wave function is limited to doubles only (D), doubles truncation plus neglect of weak pairs (D+WP), and truncation of both doubles and singles (D+S) (NB: T̂ 2 ratios for the D and D+S truncations are by definition identical). a

0.00550 0.0465 0.181 0.415 −130.4 −175.1 −189.7 −194.0 295.4 312.3 318.3 319.9 0.00574 0.0535 0.248 0.644 −101.2 −142.2 −156.7 −162.7 297.3 316.0 322.8 324.7 −72.8 −107.8 −127.9 −164.5

D

252.8 301.1 319.1 323.9

D+S α T2 ratio D+WP

[α]D [α]D

D+WP

T2 ratio

α

D+S

[α]D

α

[α]D

T2 ratio

α

PNO OSV

Table 2. Dynamic Polarizabilities (au) and Specific Rotations [deg dm−1 (g/mL)−1] Computed at 589 nm at the CC2/aug-cc-pVDZ Level of theory for (S)-1-Phenylethanol as a Function of the Occupation-Number Cutoff in the OSV and PNO Approachesa

PNO approach). The T̂ 2 ratio for all three methods naturally decreases with increasing chain length with the PNO approach remaining below 1% even at a cutoff of 10−7. In addition, the PAO approach yields polarizability errors roughly equivalent to the OSV method with an occupation cutoff of 10−7. However, we note that the use of expanded PAO domains selected through the coupled-perturbed Hartree−Fock (CPHF) method, as described in refs 61 and 62, yields lower errors (not shown) but forces the use of the entire virtual space. Although the local correlation methods correctly reproduce the (negative) sign of the specific rotation of all three fluoroheptanes, the relative localization errors are significantly larger than for polarizabilities. The OSV method with a threshold of 10−7 yields the smallest errors of ∼9% (and the largest T̂ 2 ratios), whereas the PAO approach gives the largest errors of up to nearly 20−30%. The OSV and PNO approaches produce mostly consistent relative errors, though the most effective balance in terms of errors (10%) relative to cost (T̂ 2 ratio of

A Comparison of Three Approaches to the Reduced-Scaling Coupled

Recommend Documents