Domain-based Local Pair Natural Orbital version of Mukherjee's state

Article Cite This: J. Chem. Theory Comput. 2018, 14, 1370−1382

pubs.acs.org/JCTC

Domain-Based Local Pair Natural Orbital Version of Mukherjee’s State-Specific Coupled Cluster Method Jiri Brabec,† Jakub Lang,†,§ Masaaki Saitow,‡ Jiří Pittner,† Frank Neese,‡ and Ondřej Demel*,† †

J. Heyrovský Institute of Physical Chemistry, v.v.i., Academy of Sciences of the Czech Republic, Dolejškova 3, 18223 Prague 8, Czech Republic ‡ Max Planck Institute for Chemical Energy Conversion, 45470 Mülheim an der Ruhr, Germany § Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles University in Prague, Hlavova 2030, 12840 Prague 2, Czech Republic ABSTRACT: This article reports development of a local variant of Mukherjee’s state-specific multireference coupled cluster method based on the domain-based pair natural orbital approach (DLPNO-MkCC). The current implementation is restricted to connected single and double excitations and model space with up to biexcited references. The performance of the DLPNOMkCCSD was tested on calculations of tetramethyleneethane. The results show that above 99.9% of the correlation energy was recovered, with respect to the conventional MkCC method. To demonstrate the applicability of the method to large systems, singlet−triplet gaps of triangulene and bis(1-(2,6-diisopropylphenyl)-3,3,5,5-tetramethylpyrrolidine-2-ylidene)beryllium complex were studied. For the last system (105 atoms), we were able to perform a calculation in cc-pVTZ with 2158 basis functions on a single CPU in less than 9 days.

1. INTRODUCTION The coupled cluster method1 is one of the most successful approaches for the description of dynamical correlation. Among its advantages are compactness of the exponential parametrization of the wave function, size-extensivity, invariance to orbital rotation within orbital subspaces, and a systematic hierarchy of approximations converging to the full configuration interaction limit. However, the traditional single-reference CC approach also has several disadvantages. One of them is relatively poor performance at low truncation of the cluster operator for systems where nondynamical correlation is important. However, the multireference generalization of coupled cluster methods is not unique, and many types of approaches have been developed over years. These methods can be divided into several classes: Hilbert space methods,2−7 which are based on CAS-like reference and use a distinct cluster operator for each reference; Fock space and internally contracted methods8−17 that also use a CAS-like reference but the cluster operator is common for all references; and methods that keep the framework of single-reference methods.18−37 In Hilbert space approaches, the wave operator is assumed in the form of a Jeziorski-Monkhorst ansatz.2 For state universal methods, all of the Heff eigenvalues are physical, whereas for state-specific methods,5−7 only one has a physical meaning. Among the state-specific Hilbert space approaches, the Mukherjee’s variant (MkCC)5,38 is most commonly used due to its rigorous size-extensivity. Other advantages include © 2018 American Chemical Society

insensivity to intruder problem, which plagues the state universal approaches, and a simplicity of structure of the working equations, which enables reusing much of singlereference code in its implementation. This method was developed in many groups, including Mukherjee, Schaefer, Gauss, Kállay, and Pittner.39−43 Nowadays, codes with iterative and noniterative triple excitations,41,44 higher excitations,42 analytic gradient,45 F12 variant,46 and uncoupled approximation,47,48 as well as massively parallel codes49 are available. Another drawback of the canonical CC method is the high computational cost due to the unfavorable scaling (O(n6) for CCSD, O(n7) for CCSD(T)). This represents a severe limitation for applicability to large molecules. One possibility how to overcome the high cost of calculation is massive parallelization. CC codes using distributed memory tensor networks, e.g., Super Instruction Assembly Language,50,51 Tensor Contraction Engine,49,52−54 TiledArray,55 and Cyclops Tensor Framework,56 have been developed. However, this does not remove the steep scaling, and the methods are still applicable only to medium-sized molecules, even on the most powerful machines available today. Furthermore, only a fraction of computational chemistry groups have an unlimited access to clusters with tens or hundreds of thousands of cores, which are needed for these calculations. Received: November 28, 2017 Published: January 18, 2018 1370

DOI: 10.1021/acs.jctc.7b01184 J. Chem. Theory Comput. 2018, 14, 1370−1382

Article

Journal of Chemical Theory and Computation

In order to achieve near-linear scaling, the above-mentioned scheme was improved by Riplinger et al., yielding the domainbased local PNO coupled cluster method (DLPNO-CCSD).81 The main differences compared to LPNO-CCSD are (i) construction of the PNOs as linear combinations of PAOs, (ii) setting up the connected singles in the PNO basis, (iii) prescreening of the contributing pair energies using dipole approximation, and (iv) efficient integral transformation based on the PAO basis. In this way, the bottlenecks of the LPNOCCSD method are removed and the resulting method can be applied to systems with over 500 atoms and 10 000 bases of basis set functions,81−83 where the cost of correlation treatment is comparable or even smaller than the cost of the preceding Hartree−Fock calculation. At the same time, the robustness of the PNO approaches is preserved. Nowadays, the PNO-based approaches are investigated by several groups, including Werner and Hättig,83−91 and are widely applied for studies of real chemical interest.92−101 In this article, we report development of the domain-based local PNO variant of Mukherjee’s state-specific multireference method (MkCC). The MkCC was selected because of the simple structure of working equations closely resembling the single-reference CCSD, as well as size-extensivity, which is crucial for large molecules. The DLPNO-MkCC method has been implemented in the ORCA package. To test the performance of the DLPNO-MkCCSD method, benchmark calculations of isomers of naphthynes, the rotational barrier of the tetramethyleneethane, [Be(MeL)2] complex, where MeL = 1-(2,6-diisopropylphenyl)-3,3,5,5-tetramethylpyrrolidine-2-ylidene, and triangulene were carried out. The TME and naphthyne systems were used also in the previous LPNOMkCCSD study102 and are thus suitable to test the accuracy and efficiency of the DLPNO-MkCCSD method.

A well-established way to reduce the scaling are the local approaches introduced by Pulay,57,58 which are based on the short-range (R−6) character of electron correlation. In the basis of localized orbitals, the matrix of HN operator is sparse, and nonzero elements are only encountered when all of the indices correspond to spatially closed orbitals. For occupied orbitals, standard Foster−Boys and Pipek Mezey procedures59,60 are commonly used. For virtual space, there are several possibilities. In the original local coupled cluster method of Werner and Schütz,61−63 the projected atomic orbitals (PAOs) are used. These are a priori local and can be easily used to form domains corresponding to local occupied orbitals. Pairs of occupied orbitals are classified according to the real space distance of the respective centers. The strong pairs are treated at the full level of theory (CC), weak and distant pairs only approximately (MP2), while very distant pairs are neglected altogether. In this way, linear scaling methods were developed.64,65 The remaining disadvantage is the sensitivity of the results to the values of cutoff parameters. Many other local approaches have been developed based on various variants of divide-and-conquer strategy: fragment molecular orbital method,66 the divide-and-conquer method,67 divide−expand−consolidate method,68,69 the incremental scheme,70 and the cluster-in-molecule approach,71 among many others. Here, the system is divided into fragments. Series of separate calculations are performed for individual fragments or their pairs, and therefore, the methods are embarrassingly parallel. The drawbacks of these method are redundancies originating from overlapping fragments and the potentially limiting high cost of individual fragment calculations for large fragments. An alternative approach is to use a virtual space spanned by pair natural orbitals (PNOs).72,73 The PNOs were introduced in the 1960s by Edmiston and Krause and were shown to provide a compact parametrization of virtual space, yielding an accelerated convergence of the CI expansion. The CISD and CEPA code of Meyer based on PNOs74,75 provided state-ofthe-art description by that time for a very modest computational cost. The PNO approach was also used for the multireference MR CISD method by Taylor76 and MR CEPA by Fink and Staemmler.77 The full potential of PNOs could only be realized in combination with local correlation approaches, modern hardware, and modern integral transformation technologies, especially the density fitting (DF) or resolution of the identity (RI) methodology.78 In late 2000s, the PNO technique was used in the LPNO CEPA and LPNO-CCSD.79,80 Within this scheme, each pair of localized occupied orbitals is assigned a separate set of PNOs. These orbitals provide the most compact description of the virtual space of the corresponding pair and are as local (or nonlocal) as the physical situation requires. In this approach, only a limited number of cutoff parameters are used: TCutPairs for contributing electron pairs as the minimum MP2 pair energy for the pair to be kept; TCutPNO as the minimum occupation number of a PNO to be kept; and TCutMKN for expansion of the auxiliary basis. None of these cutoff parameters is based on real space distance. The advantages of this approach include compactness of the LPNO-based wave function, high accuracy measured as the percentage of correlation energy relative to calculation with canonical orbitals, smooth dependence of the remaining error on the cutoff parameters, and black box character of the method.

2. THEORY This section is divided into three parts. First, we review basics of MkCC theory. Next, we describe how to obtain PNOs for subsequent MkCC calculation. Finally, we use these PNOs to formulate the DLPNO-MkCC method. 2.1. Basics of MkCC Method. The reference function is assumed to be a linear combination of reference configurations spanning a model space M

|ΨωP ⟩ =

∑ Cμω|Φμ⟩ (1)

μ=1

Furthermore, let us define the wave operator Ω̂ω as |Ψω⟩ = Ω̂ω|ΨωP ⟩

(2)

and the effective Hamiltonian Heff as eff ̂ ̂ Ω̂ωP ̂ Ĥ = PH

(3)

where Ψω stands for the exact wave function, P for the projection operator on the model space, and index ω for the studied state. eff Ĥ |ΨωP ⟩ = ,ω|ΨωP ⟩

(4)

The matrix representation within the model space reads as

∑ Hμνeff Cνω = ,ωCμω ν

1371

(5) DOI: 10.1021/acs.jctc.7b01184 J. Chem. Theory Comput. 2018, 14, 1370−1382

Article


the system size, the number of ip and pq pairs is thus negligible compared to ij pairs, and therefore, the prescreening for ip and pq pairs is neither required nor performed. The next step is evaluation of perturbation energies of the “surviving” pairs and subsequent construction of PNOs. We have devised two schemes, one based on MP2 and one based on NEVPT2. In the former scheme, both the pair energies and the PNOs are constructed at the MP2 level. For each pair of occupied/ active spatial orbitals surviving the prescreening step we calculate first-order amplitudes

The wave operator is taken in the form of the Jeziorski− Monkhorst ansatz M

Ω̂ω =

∑ eT̂(μ) Φμ⟩⟨Φμ (6)

μ=1

where T(μ) is a cluster operator assigned to the μth reference configuration. Assuming that a complete model space (CMS) is employed, the amplitudes corresponding to internal excitations are zero due to intermediate normalization. The remaining cluster amplitudes are obtained by solving the cluster equations ⟨Φ(μ)ϑ |e +

−T (μ)

He

tia′ j′′b ′(μ)

T (μ)

|Φ(μ)⟩Cμω

∑ Hμνeff Cνω⟨Φ(μ)ϑ |e−T(μ)eT(ν)|Φ(μ)⟩ = 0 ν≠μ

K ia′ j′′b ′

=

ϵa ′(μ) + ϵb ′(μ) − ϵi ′(μ) − ϵj ′(μ)

(9)

where indices a′,b′ correspond to PAOs and are restricted to the domain of the i′j′ pair, K ia′ j′′b ′ is the exchange integral and ϵa ′(μ) is the diagonal element of the Fock matrix of the μth reference. Then we use these amplitudes to calculate the unrestricted MP2 energy. The UMP2 energy is summed over spin and averaged over M′ references where i′, j′ are both occupied. If the average pair energy is greater than or equal to the TCutPairs cutoff parameter (default value 1.0 × 10−4), the pair is kept, otherwise it is discarded. No pq pairs are neglected in this step because the number of these pairs is negligible and all of these pairs are assumed to be important. For the surviving pairs, the RHF MP2 densities for the αβ spin case are constructed for each reference

(7)

where index ϑ corresponds to an excitation outside of the model space. Here, the first term is identical to the single-reference CC residual calculated with a set of cluster amplitudes corresponding to the μ-th reference configuration, and only the second term, referred to as the coupling term, connects cluster amplitudes of different references. The matrix elements of the effective Hamiltonian are given as eff eff Hμν = ⟨Φ(μ)|Ĥ Φ(ν)⟩ ̂ = ⟨Φ(μ) Ĥ Φ(μ)⟩δμν + ⟨Φ(μ)|Ĥ N(ν)eT(ν)|Φ(ν)⟩

†

(8)

Di ′ j ′(μ) = Tĩ ′ j ′(μ)Ti ′ j ′(μ) + Tĩ ′ j ′(μ)Ti†′ j ′(μ)

2.2. Prescreening and PNO Generation for DLPNOMkCCSD. In this section, we describe how to generate a common set of PNOs for all reference configurations used in the subsequent DLPNO-MkCCSD calculation. In this paragraph, indices i, j, ... stand for inactive (occupied) orbitals, p, q, ... for active, a, b, ... for external orbitals, i′, j′, ... for orbitals occupied in a given reference, and a′, b′, ... for orbitals unoccupied in a given reference. The DLPNO-MkCC calculation starts with a CAS SCF calculation. In the next step, the inactive orbitals are localized using the standard Foster−Boys procedure. The orbital domains of PAOs are set up in the same way as those in the single-reference DLPNOCCSD.81 The PAOs are canonicalized with respect to the CAS SCF Fock matrix. In the next step, the prescreening of pair energies is performed. At this step, orbital pairs involving active orbitals, diagonal pairs, as well as pairs with differential overlap greater than the TCutDO cutoff parameter (default 1.0 × 10−2) are always included. For the remaining pairs, the pair energy is estimated at the strongly contracted NEVPT2 level103−105 using the dipole approximation. The contributions Vij (i.e., class including ij → pq excitations, the active indices are omitted in the name of the class), Vija (i.e., class including ij → ap excitations), and Vijab (i.e., class including ij → ab excitations) are all included in the pair energy. Details can be found in Ref 86, e.g., in eq 24. Pairs with estimated pair energy greater than or equal to TCutPre (default value 1.0 × 10−6 Hartree) are kept for further calculation, while the remaining pairs (i.e., offdiagonal ij pairs with low pair energy and differential overlap) are neglected, and only their pair energies are collected. As was established earlier, the dipole approximation works sufficiently well for those pairs with small pair energy contributions. Here, the number of active orbitals is assumed to be independent of

(10)

where Tĩ ′ j ′(μ) =

1 (4Ti ′ j ′(μ) − 2Tj ′ i ′(μ)) 1 + δi ′ j ′

(11)

The densities are then averaged over M′ references Di ′ j ′ =

1 M′

M′

∑ Di′ j ′(μ) μ

(12)

where both i′ and j′ are occupied. Subsequent diagonalization of the Di′j′ matrix yields the PNOs Di ′ j ′dix′ j ′ = nix′ j ′dix′ j ′

(13)

that are used for all spin cases of the pair. Only those PNOs with occupation numbers nix′ j ′ greater or equal to TCutPNO (default value 3.33 × 10−7) are kept. The energy correction for neglected PNOs is calculated for each reference separately as the difference between unrestricted MP2 energy in the PAO basis and the truncated PNO basis. This correction is added to pair energies of pairs neglected in previous steps, for each reference. The PNOs for connected singles are constructed as PNOs of diagonal i = j or p = q pairs with a reduced value of the cutoff parameter of TCutPNOSingles = 0.03 × TCutPNO, analogous to eq 13. In the alternative PNO construction scheme, the semicanonical partially contracted N-electron valence perturbation theory (PC-NEVPT2) energies are calculated for the “surviving” pairs. For the ij pairs, the energy contributions from Vijab, Vija, and Vij classes are added together. Detailed expressions are presented in Ref 86. Similarly, the Viab, Via, and 1372


Article


for the connected doubles. The first terms (direct terms) are calculated using a slightly modified single-reference open-shell DLPNO-CCSD code83 in ORCA. The changes are needed to properly handle different occupation of active orbitals for the references. Because in the SR DLPNO-CCSD code all α-spin singly occupied molecular orbitals (SOMOs) are occupied and all β-spin SOMOs are unoccupied, it is efficient to generate amplitudes and intermediates within these specific orbital windows, which is different from the orbital layout used here (see above). Further changes are needed to zero out the amplitudes and residuals corresponding to excitations from unoccupied/to occupied spin orbitals for the respective reference. This technique was used previously in the LPNOMkCCSD method.102 The second terms in eqs 17 and 18, the coupling terms, can be evaluated as

Vi contributions are added for ip pairs, as are Vab and Va added for pq terms. Among ij and ip pairs, only those with PCNEVPT2 pair energies greater than or equal to TCutPairs are included for the subsequent coupled cluster step. Again, no pq pairs are neglected. For the surviving pairs, the PC-NEVPT2 pair densities are constructed. For ij pairs, this includes contributions from classes Vijab and Vija, for ip pairs contributions from classes Viab and Via, and for pq pairs Vab and Va contributions. The densities are diagonalized in order to obtain the PNOs. Only those PNOs with occupation number greater than or equal to TCutPNO are kept. The energy contribution of the remaining PNOs is calculated as the difference between the pair energy calculated in the PAO basis and that for the truncated PNO basis. It was expected that the latter scheme will be preferable due to more satisfactory treatment of the multireference nature of the studied systems. However, as detailed below, this expectation was not confirmed by benchmark calculations. 2.3. DLPNO-MkCCSD Working Equations. In the next step, we make use of the fact that we have a common set of PNOs for all of the references. This means that we can construct a common set of molecular integrals for all references in order to save computational effort for the rather costly integral transformation as well as to cut disk storage requirements. For the same reason, we want to have the same PNOs for all spin cases (α−α, α−β, and β−β). However, the active orbitals are occupied for some references and unoccupied for others. It is therefore convenient to include the active orbitals to both occupied and virtual orbital windows. The virtual orbital space is spanned by PNOs, and these are in principle different for each pair of occupied orbitals. To keep the scheme simple, we keep the active orbitals unchanged and include them to the set of PNOs of all pairs for occupied orbitals. Therefore, we also adapt the notation in the subsequent text: indices i, j, ... stand for both internal and active orbitals; a, b, ... for active and external in the MO basis; and ai̅ j, b̅ij, ... for both active and external in the PNO basis of the ij pair. The cluster operators are truncated at the connected singles and doubles level and are constructed in the PNO basis set as described above. T (μ) = T1(μ) + T2(μ) T1(μ) =

a b̅

3. RESULTS AND DISCUSSION In this section, we present results of benchmark calculations on tetramethyleneethane (TME), isomers of naphthynes, the [Be(MeL)2] complex, where MeL = 1-(2,6-diisopropylphenyl)3,3,5,5-tetramethylpyrrolidine-2-ylidene, and triangulene. In all calculations, we utilized CMS(2,2) orbitals obtained from CASSCF calculations. 3.1. Accuracy of the DLPNO-MkCCSD Approach with Respect to Canonical MkCCSD. In order to compare the performance of our new method with the canonical version of MkCCSD, we performed benchmark calculations on two selected systems. The size of those systems cannot be too small because the number of the pairs and the PNOs would not be appropriate for benchmarking purposes. On the other hand, it cannot be too large; otherwise, we could not perform canonical MkCCSD calculations. We have thus chosen TME and naphthyne isomers, which have strong multireference character requiring CMS(2,2) model space and where the DLPNOMkCCSD results can be compared with recently published LPNO-MkCCSD results (see ref 102). 3.1.1. Rotational Barrier of TME. First, we calculated the rotational barrier of TME using the geometries from ref 44. For

i

(15)

ai̅

∑ ∑ t(μ)ija̅ b ̅ aa†̅ ab†̅ ajai ij ij

ij

ij

ij

aij̅ bij̅

(16)

The DLPNO-MkCCSD cluster eq 7 becomes ⟨Φ(μ)iaii̅ |e−T(μ)H eT(μ)|Φ(μ)⟩Cμα +

∑ Hμνeff Cνα⟨Φ(μ)ia̅ i |e−T(μ)eT(ν)|Φ(μ)⟩ = 0 i

(17)

ν≠μ

for the connected singles and a b̅

⟨Φ(μ)ijij̅ ij |e−T(μ)H eT(μ)|Φ(μ)⟩Cμα +

∑ Hμνeff Cνα⟨Φ(μ)ija̅ b ̅ |e−T(μ)eT(ν)|Φ(μ)⟩ = 0 ij ij

ν≠μ

a b̅

⟨(Φ(μ))ijij̅ ij |e−T(μ)eT(ν)|Φ(μ)⟩ = Δtij ij̅ ij(ν /μ , μ) 1 a b̅ + P(ij)P(ab)Δti ij̅ (ν /μ , μ)Δt j ij(ν /μ , ν) (20) 2 a... a... for the connected doubles. Here, Δti... (ν/μ,μ) = ti... (ν/μ) − a... a... ta... i... (μ) and ti... (ν/μ) are equal to ti... (ν) if all i... are occupied and all a... unoccupied for both the ν-th and μ-th references and zero otherwise. In contrast with the previous LPNO-MkCCSD variant, the singles couplings are formulated in the singles PNO basis. Because the PNOs are identical for all of the references, this causes no complications. The doubles couplings in eq 20 are equivalent to their counterparts in LPNO-MkCCSD. The diagonal elements of the effective Hamiltonian are calculated as the CC energies corresponding to the respective references plus the correction for neglected pairs and truncated PNOs, denoted as EREST. If we restrict ourselves to active spaces with mutually up to biexcited references, the off-diagonal Heff μν elements can be collected as the residuals corresponding to the internal excitations.

i

i

(19)

for the connected singles and

(14)

∑ ∑ t(μ)ia̅ aa†̅ ai

T2(μ) =

⟨(Φ(μ))iai̅ |e−T(μ)eT(ν)|Φ(μ)⟩ = Δtiai̅ (ν /μ , μ)

(18) 1373


Article


Figure 1. (a) Rotational barrier of TME calculated by the DLPNO-MkCCSD method with different sources of PNOs and (b) percentage of correlation energy recovered for TME by DLPNO-MkCCSD with different sources of PNOs.

this system, a model space consisting of two references HOMO2 LUMO0 and HOMO0 LUMO2 was employed; the remaining two references do not contribute because of symmetry consideration. The standard correlation-consistent cc-pVTZ basis set106 in spherical representation was used throughout the calculations. The aim of these calculations was to answer several basic questions: (a) what percentage of correlation energy is recovered, and how does that depend on cutoff parameters, (b) what is the accuracy for relative energies, and (c) how does it all depend on the method used for the PNO construction. As a first step, we calculated the rotational barriers as well as percentages of correlation energy recovered with the default cutoff settings with both NEVPT2 and MP2 PNOs. The results are shown in Figure 1a,b. In both cases, approximately 99.9% of the canonical correlation energy was recovered. However, for relative energies, the two curves show significant differences especially in the region of 30−50° where the multireference character is more complicated (see ref 107). Clearly, DLPNOMkCCSD based on MP2 PNOs gives a better description, which can be also seen from the values of the nonparallelity error (NPE): 0.17 kcal/mol for MP2 PNOs vs 0.66 kcal/mol for NEVPT2 PNOs. To understand the origin of these differences, we checked individual pair energies at the MP2 and NEVPT2 levels and the corresponding number of PNOs. As a reference geometry, we took the angle 40°. The most interesting difference concerned the energy of the pq pairs, where the NEVPT2 pair energies were almost 2 orders of magnitude smaller than their MP2 counterparts (0.2 vs 14.8 mHartree). Correspondingly, the number of PNOs obtained for these pairs was at the NEVPT2 level smaller by a factor of 4. To illustrate these findings, we have plotted diagrams showing ip and pq pair energies for NEVPT2, MP2, MkCCSD at first iteration, and converged MkCCSD. We also plotted NEVPT2 pair energies summed over the active index in order to obtain an adequate quantity comparable with MP2 and MkCCSD results (MP2 or MkCCSD pair energies are reference-specific; for the example of TME in CMS(2,2) containing only two closed-shell references, only one of the ip pairs contributes to each reference for a given inactive orbital i, while for NEVPT2 all pair energies are nonzero and these are not reference-specific). The results are presented in Figure 2. In all cases, MP2 ip pair energies are consistently closer to their MkCCSD counterparts, slightly overestimating MkCCSD. NEVPT2 energies are smaller than MP2 energies in all cases,

Figure 2. Pair energies of ip (on the left) and pq (on the right) pairs for TME (α = 40°) estimated at NEVPT2, MP2, and MkCCSD levels, sorted by their magnitude. NEVPT2 (sum) stands for NEVPT2 pair = ∑p ENEVPT2 . energies summed over the active label: ENEVPT2(sum) i ip Default thresholds were used.

and NEVPT2(sum) energies are smaller in all cases except for the four largest contributions, where they are much larger than the ones from other methods. For pq pairs, the difference between NEVPT2 and MP2 or MkCCSD is dramatically larger. Overall, the MP2 approach gives a much better description of pairs involving active orbitals, which is crucial for proper construction of ip and pq PNOs. The same discrepancy in pq pair energies was observed for other systems that we studied: naphthynes as well as the beryllium complex. Furthermore, significant differences were observed in the energy of the neglected pairs and PNOs (EREST) for the ip pairs. For NEVPT2 PNOs, the EREST contribution was 10 times larger than that for MP2 PNOs, which indicates a worse description in the PNO basis for given pairs. The same was not observed for the pq pairs, where the EREST is approximately 0.03 mHartee. However, that was only due to the very small NEVPT2 pair energy. When the MP2 energy is calculated for the pq pairs using NEVPT2 PNOs, EREST rises by 2 orders of magnitude, which leads to serious overestimation of the correlation energy. Again, this shows that, surprisingly, the NEVPT2 PNOs are not well suited for Hilbert space MRCC methods. To verify this further, we tried to selectively tighten the TCutPNO cutoff parameter by a factor 0.01 for pq pairs and 0.1 for ip pairs. This has led to a decrease in NPE from 0.66 to approximately 0.33 kcal/mol. Furthermore, we performed a series of calculations with mixed PNOs: (a) MP2 PNOs for pq 1374


Article


Figure 3. (a) Rotational barrier of TME with respect to TCutPNO and (b) corresponding NPE (in kcal/mol) as a function of TCutPNO.

Figure 4. Percentage of canonical correlation energy recovered for three angles (α = 0°, 40°, and 90°) with respect to (a) TCutPNO and (b) TCutPairs thresholds.

10−5 and TCutPNO = 1.0 × 10−7, the NPE drops below 0.1 kcal/mol. Furthermore, we have studied how the percentage of correlation energy changes with the cutoffs used in DLPNOMkCCSD. The results are shown in Figure 4a for TCutPNO, Figure 4b for TCutPairs, and Figure 5 for TCutDO. In all cases, the dependency on cutoff parameter is smooth. The error is very small or negligible with TCutPairs smaller than 1 × 10−4 Hartree (default value) because the number of cut pairs is very small. As TCutPairs increases, the percentage of recovered correlation energy (as well as the error) steeply increases to more than 100%, which is caused by overestimation of the pair

pairs, NEVPT2 for ip and ij pairs and (b) MP2 PNOs for pq and ip pairs, NEVPT2 for ij pairs. These curves are shown in Figure 1a. The corresponding NPEs of 0.30 or 0.35 kcal/mol, respectively, lie between the NPEs for NEVPT2 and MP2 PNOs. It should be also pointed out that part of the relative success of calculations with MP2 PNOs may be connected with the relatively small active space employed. In the CMS(2,2) space, only one reference has nonzero contributions for a given pq pair because at least one of the orbitals p, q is unoccupied for the remaining references. Therefore, the PNOs for these pairs are optimized for the only contributing reference. When only two references of the CMS(2,2) space are used, the same is true also for the ip pairs. On the basis of the results we obtained in these benchmark calculations, we employed MP2 PNOs in all subsequent calculations. In Figure 3a, we compare the rotational barrier with respect to the TCutPNO threshold, where all PNOs were constructed from the MP2 amplitudes. These results show very similar curves for all of the values of TCutPNO. For TCutPNO smaller than 1.0 × 10−6, the relative DLPNO-MkCCSD energy with respect to α = 90° is below the canonical MkCCSD energy at each point. The NPE associated with the TCutPNO parameter is shown in Figure 3b. For the default threshold, the NPE is only 0.17 kcal/mol; for lower values of TCutPNO, it decreases to 0.12 kcal/mol. However, similar improvement can be obtained by tightening TCutPairs to 5.0 × 10−5. This can be of practical importance because the computational cost for the two sets of calculations is very similar. With TCutPairs = 5.0 ×

Figure 5. Percentage of correlation energy recovered for three TME molecules for α = 0°, 40°, and 90°, with respect to the TCutDO threshold. 1375


Article


default TCutPNO threshold, approximately 99.90−99.93% of the correlation energy is recovered for all four isomers. This represents a significant improvement over LPNO-MkCCSD, which recovered less than 99.8% (ref 102). The minimum is between 1.0 × 10−6 and 1.0 × 10−5; however, DLPNOMkCCSD still recovers more than 99.8% of the correlation energy. The smallest value 99.82% is for 1,8-naphthyne at TCutPNO = 1.0 × 10−5, while 2,6-naphthyne has a minimum of 99.87% at TCutPNO = 1.0 × 10−6. The minimum for LPNOMkCCSD was approximately 0.2% lower.102 Next, the isomerization energies with respect to 2,6naphthyne were calculated. Canonical MkCCSD yields values of 0.57 kcal/mol for 2,7-naphthyne, 2.18 kcal/mol for 1,6naphthyne, and 5.03 kcal/mol for 1,7-naphthyne. For DLPNOMkCCSD, the errors of isomerization energies with respect to canonical MkCCSD are listed in Figure 7a,b. At the default value of TCutPNO, the deviations for 2,7-naphthyne are 0.06 kcal/mol and those for 1,8- and 1,7-naphthyne are approximately 0.26 kcal/mol. With TCutPairs = 5 × 10−5 and default TCutPNO, i.e., without significant increase of computational cost, the error is reduced to approximately 0.17 kcal/mol for both 1,7- and 1,8-naphthyne. Importantly, when TCutPNO is reduced to 1.0 × 10−7 and TCutPairs to 5 × 10−5, the error is reduced to only 0.06 kcal/mol or less for all three isomers. For 1,7-naphthyne, this represents an improvement over the LPNO-MkCCSD result102 by approximately 0.1 kcal/mol from 0.36 to 0.26 kcal/mol. For the other two isomers, the effect is smaller; for 1,8-naphthyne, the error is increased by 0.05 kcal/ mol, and for 2,7-naphthyne, it is reduced by 0.02 kcal/mol. However, the tightening of TCutPNO at the LPNO-MkCCSD level had only a small effect. An explanation for the improvement for tight values of cutoffs is likely due to the inclusion of the terms in cluster equations that are neglected at the LPNO level. These consist mainly of three-external dressings of four-external doubles contributions and several terms in dressing of two-external doubles contribution. It can be expected that these terms are larger in the multireference case, where the T1 amplitudes are in general larger due to nonzero Fock elements. Figure 8a,b shows the change of the percentage of recovered correlation energy for 1,8-naphthyne as a function of TCutPairs, TCutDO, and TCutMKN thresholds. One can clearly see that for TCutPairs the energy remains constant up to 5.0 × 10−5. For larger thresholds, the energy increases monotonically as the number of neglected pairs grows significantly, and their pair energy is

energy of the neglected pairs at the perturbation theory level. Next, let us investigate the TCutPNO curve. For very large values of TCutPNO, only a few PNOs per pair survive and the MP2 used for the neglected PNOs significantly overestimates the correlation energy. The curve has then a minimum near TCutPNO = 1.0 × 10−6, and then the percentage of recovered correlation grows ideally to 100%. The obtained curve is highly satisfactory because it is smooth; over 99.9% of the correlation energy is recovered along all of the curve, and for very small TCutPNO, it approaches 100% of the correlation energy. We also calculated the energy with respect to the TCutDO (Figure 5) threshold. The energy is almost constant with TCutDO up to 0.01, which is also the default value. With looser values of TCutDO, the amount of correlation energy TCutPNO quickly decreases as important atoms are cut out of the domains. 3.1.2. Naphthyne Isomers. As the second benchmark, we have investigated four isomers, (1,7), (1,8), (2,6), and 2,7naphthyne, which exhibit multireference character. All calculations employed the cc-pVTZ basis set. For 1,7-naphthyne, also the two open-shell references of the CMS(2,2) space were included, while only the two closed-shell references are needed for the remaining three isomers due to the symmetry considerations. The geometries for single-point calculations were taken from ref 49. The percentage of correlation energies recovered by the DLPNO-MkCCSD method are shown in Figure 6. With the

Figure 6. Percentage of correlation energy recovered for naphthyne isomers by DLPNO-MkCCSD as a function of TCutPNO.

Figure 7. Relative ΔΔE computed as ΔEx − ΔE2,6, where ΔEx = EDLPNO‑MkCCSD − EMkCCSD calculated with (a) TCutPairs = 1.0 × 10−4 and (b) x x −5 TCutPairs = 5.0 × 10 . 1376


Article


Figure 8. Percentage of correlation energy recovered for the 1,8-naphthyne isomer by DLPNO-MkCCSD with respect to (a) TCutPairs and (b) TCutDO and TCutMKN thresholds.

Å), the differences are approximately 20−25 μHartree, and in the dissociation limit, they disappear because the weights are the same in both sets of calculations. Therefore, we conclude that the errors originating from the weighting method are significantly smaller than other errors encountered in DLPNO approaches. 3.2. Application to Larger Systems. To illustrate the effectiveness of the DLPNO-MkCCSD approach, we have performed calculations on the [Be(MeL)2] complex, where MeL = 1-(2,6-diisopropylphenyl)-3,3,5,5-tetramethylpyrrolidine-2ylidene, and on triangulene diradical (see Figure 10). The

overestimated at the perturbative level. At the default value of 1.0 × 10−4, the error is only 0.01%. TCutDO and TCutMKN have an opposite effect because for larger values, the cutoffs of the orbital domains (TCutDO) and fitting domains (TCutMKN) start to lose essential components, which leads to substantial loss of the correlation energy contributions. Figure 8b shows that for up to default values of TCutDO = 1.0 × 10−2 and TCutDO = 1.0 × 10−3 the percentage is almost constant, and for larger thresholds, the correlation energy significantly decreases. 3.1.3. Potential Curve for F2. In the process of generation of PNOs, the pair density matrices are averaged over contributing references. In principle, this could potentially lead to significant errors when the references have considerably different weights. In order to investigate this, we have calculated two curves obtained by the DLPNO-MkCCSD method in the cc-pVQZ basis, one where all contributing references have the same weight and one where weights from preceding CAS SCF calculations are used. Here, the CAS SCF weights of the references change drastically along the curve. The difference between the two curves is shown in Figure 9. As is clear from the graph, the differences between the two curves are very small. The largest difference of approximately 60 μHartree is observed at 1.6 Å. Near the equilibrium (1.42

Figure 10. Molecules (a) [Be(MeL)2] and (b) triangulene.

[Be(MeL)2] complex was recently prepared as one of the first compounds containing Be(0) as a metal center, where obtained results indicate that the stability is given by a strong threecenter, two-electron bond C−Be−C (ref 108). Triangular nonKekule structures as triangulene have a high-spin state as the ground state, in this particular case the triplet state. The first excited singlet state has an open-shell character; the published S−T energy gap is 0.6−0.8 eV at MR-CISD level (ref 109) or 0.16 eV at the spin-polarized DFT level (ref 110). The closedshell singlet state has been predicted to be about 0.35 eV above the ground triplet state at the spin-polarized DFT level (ref 110). In our calculations, we focused at closed-shell singlet states. For a proper description of multireference character of singlet states, we employed CMS(2,2) and utilized the ccpVTZ basis set. The geometry of [Be(MeL)2] was taken from ref 108, and for triangulene, we used the geometry from ref 109. The triangulene has 828 basis functions, while the [Be(MeL)2] complex was computed using 2158 basis functions. As the auxiliary basis set for RI projection, we utilized ccpVQZ/C, giving 3564 or 9086 basis functions, respectively. The results obtained for default thresholds are shown in Table 1. The number of RHF pairs for triangulene that were

Figure 9. Difference between the DLPNO-MkCCSD energy of F2 calculated in the cc-pVQZ basis with (a) all contributing references having the same weight in the construction of PNOs and (b) reference weights obtained from preceding CAS SCF calculation. 1377


Article


We also performed a timing test on a single node (Table 2), using a 1 CPU core and 128 GB of memory. As expected, the

Table 1. Total Number of Pairs (Pairs), Number of Surviving Pairs, Total Number of Amplitudes That Would Appear in Canonical Calculation, Number of Optimized PNO Amplitudes in DLPNO-MkCCSD, Average Number of PNOs per Pair, Energy of the Lowest Singlet State, Corresponding Correlation Energy, Weak-Pairs Correction, and Singlet−Triplet Gap Energya system number of atoms number of basis functions number of auxiliary basis functions total number of RHF/UHF pairs number of RHF/UHF pairs kept total number of PNOs number of canonical amplitudes number of PNO amplitudes average number of PNOs per (RHF) pair total energy [au] correlation energy [au] right eigenvector of Heff weak-pair energy correction [au] S−T gap [eV]

triangulene

Table 2. Computational Timings for [Be(MeL)2]a time in s Fock matrix formation organizing maps RI-PNO integral transformation initial guess state vector update sigma-vector construction ⟨0|H|D⟩ ⟨0|H|S⟩ ⟨D|H|D⟩(0-ext) ⟨D|H|D⟩(2-ext Fock) ⟨D|H|D⟩(2-ext) ⟨D|H|D⟩(4-ext) ⟨D|H|D⟩(4-ext-corr) CCSD doubles correction ⟨S|H|S⟩ ⟨S|H|D⟩(1-ext) ⟨D|H|S⟩(1-ext) ⟨S|H|D⟩(3-ext) CCSD singles correction Fock-dressing singles amplitudes (ik|jl)-dressing (ij|ab),(ia|jb)-dressing pair energies

[Be(MeL)2]

34 828 3564 2601/5151 1505/2959 59604 1807549275 4506028 (0.25%) 39

105 2158 9086 13924/27730 3229/6322 99512 28474384260 5958516 (0.02%) 30

−843.767587 −3.493519 (−0.712,−0.703) −0.024 0.78

−1682.370676 −7.608054 (−0.937,0.350) −0.080 −0.42

The S−T gap is defined as Esinglet − Etriplet. The triplet state energy was calculated using the single-reference open-shell DLPNO-CCSD code (see ref 83 for the method and implementation). a

included in the coupled cluster treatment was 1505, i.e., 58% of the original 2601 pairs. Among the neglected pairs, 46 were excluded in the prescreening step. For [Be(MeL)2], the percentage of pairs kept was reduced to 23%. Also, the number of pairs dropped in prescreening has risen to 1536, i.e., 11% of the original pairs. In terms of the number of cluster amplitudes, for triangulene and the DLPNO-MkCCSD method, 0.25% of the canonical amplitudes is needed. For larger [Be(MeL)2], the number is only 0.02%. The dropped pairs energy contribution represented approximately 70% of the weak pairs energy, while the remaining 30% corresponds to neglected PNOs. The right eigenvectors of the effective Hamiltonian are (−0.712,−0.703) for triangulene and (−0.937,0.350) for [Be(MeL)2], respectively, which confirms strong multireference character of both systems. The singlet−triplet (S−T) gap obtained for triangulene is 0.78 eV, which is larger than the recently published 0.35 eV in ref 110. This is in line with the situation for the lower openshell singlet state, where the MR CISD value was also considerably higher than the spin-polarized DFT one (see above). A different situation was observed for [Be(MeL)2], where the singlet state is approximately 0.42 eV below the triplet state. In order to quantify the influence of the TCutPairs threshold (as discussed in the previous section) on relative energy, we performed also [Be(MeL)2] calculations with TCutPairs = 5 × 10−5. The total DLPNO-MkCCSD energy increased for about 6.4 mHartree to −1682.364315 Hartree, whereas the singlereference triplet energy converged to −1682.349709 Hartree, which is for 5.6 mHartree higher than that for the default cutoff. This results in a small change in the S−T gap from −0.42 to −0.40 eV. On the other hand, the total number of UHF pairs treated at a strongly correlated level increased by 19% from 6322 to 7499, which negatively affects timings.

64540.3 6.9 195480.1 77744.7 73.1 421561.8 0.9 0.8 55771.4 296.2 145998.3 6913.1 64900.7 574.1 2163.4 4022.2 138.0 661.5 0.4 142704.3 587.3 15950.7 2680.3 20.9

time in % 8.5 0.0 25.7 10.2 0.0 55.5 0.0 0.0 13.2 0.1 34.6 1.6 15.4 0.1 0.5 1.0 0.0 0.2 0.0 33.9 0.1 3.8 0.6 0.0

of of of of of of of of of of of of of of of of of of

sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma sigma

a

The calculation was carried out with default thresholds, using a 1 CPU core and 128 GB of memory.

most expensive part of the calculation is an evaluation of sigma vectors followed by RI-PNO integral transformation. Prescreening and PNO construction took 10% of the total execution time, and 8.5% time was spent on Fock matrix formation. The total computation time was less than 9 days. As far as sigma construction, the most time was spent on linear two-external doubles contributions ⟨D|H|D⟩(2-ext), dressing of the Fock matrix, and three-external dressings of four-external doubles contributions ⟨D|H|D⟩(4-ext-corr), which is consistent with the single-reference DLPNO-CCSD code.

4. CONCLUSIONS We have introduced a new DLPNO-MkCCSD approach, which is a Hilbert space multireference method applicable to large systems (more than 100 atoms and thousands of basis functions). Even using only one CPU core, we can perform multireference CCSD calculation with 105 atoms and 2158 basis functions in less than 9 days, which is not feasible in the canonical orbitals representation within any reasonable time. On two benchmark systems, TME and isomers of naphthynes, we have demonstrated the accuracy of the DLPNO-MkCCSD with respect to the canonical MkCCSD method. It is shown that by using the default set of thresholds one could obtain more than 99.9% of the correlation energy. Also, the nonparallelity error tested on the TME system is 0.17 kcal/mol, and errors for isomerization energies of naphthynes are smaller than 0.3 kcal/mol. We also studied two larger systems, the triangulene diradical and [Be(MeL)2] complex. The obtained singlet−triplet gap 0.78 eV for triangulene is larger than the recently published DFT 1378


Article

Journal of Chemical Theory and Computation energy, while for [Be(MeL)2] we found the singlet state as the ground state with a S−T energy gap of −0.42 eV. Compared to the previous LPNO-MkCCSD method, we have shown that DLPNO-MkCCSD is more accurate in terms of percentage of the correlation energy recovered and is applicable to larger systems due to the almost linearly scaling computational cost. Furthermore, it is possible to develop a perturbative triples correction within the DLPNO-MkCCSD scheme, which has a computational cost comparable with that of the corresponding singles and doubles method, unlike in the LPNO-MkCCSD framework. This is particularly important because it is known that the triples correction is highly beneficial for the accuracy of the MkCC method, and it will be the next step of our further development.

■

(12) Hanauer, M.; Kö hn, A. Pilot applications of internally contracted multireference coupled cluster theory, and how to choose the cluster operator properly. J. Chem. Phys. 2011, 134, 204111. (13) Hanauer, M.; Kö hn, A. Pilot applications of internally contracted multireference coupled cluster theory, and how to choose the cluster operator properly. J. Chem. Phys. 2011, 134, 204111. (14) Nooijen, M. State Selective Equation of Motion Coupled Cluster Theory: Some Preliminary Results. Int. J. Mol. Sci. 2002, 3, 656−675. (15) Datta, D.; Nooijen, M. Multireference equation-of-motion coupled cluster theory. J. Chem. Phys. 2012, 137, 204107. (16) Demel, O.; Datta, D.; Nooijen, M. Additional global internal contraction in variations of multireference equation of motion coupled cluster theory. J. Chem. Phys. 2013, 138, 134108. (17) Nooijen, M.; Demel, O.; Datta, D.; Kong, L.; Shamasundar, K. R.; Lotrich, V.; Huntington, L. M.; Neese, F. Communication: Multireference equation of motion coupled cluster: A transform and diagonalize approach to electronic structure. J. Chem. Phys. 2014, 140, 081102. (18) Oliphant, N.; Adamowicz, L. Multireference coupled-cluster method using a single-reference formalism. J. Chem. Phys. 1991, 94, 1229−1235. (19) Piecuch, P.; Oliphant, N.; Adamowicz, L. A state-selective multireference coupled-cluster theory employing the single-reference formalism. J. Chem. Phys. 1993, 99, 1875−1900. (20) Piecuch, P.; Adamowicz, L. State-selective multireference coupled-cluster theory employing the single-reference formalism: Implementation and application to the H8 model system. J. Chem. Phys. 1994, 100, 5792−5809. (21) Adamowicz, L.; Piecuch, P.; Ghose, K. B. The state-selective coupled cluster method for quasi-degenerate electronic states. Mol. Phys. 1998, 94, 225. (22) Piecuch, P. Active-space coupled-cluster methods. Mol. Phys. 2010, 108, 2987−3015. (23) Paldus, J.; Č ížek, J.; Takahashi, M. Approximate account of the connected quadruply excited clusters in the coupled-pair manyelectron theory. Phys. Rev. A: At., Mol., Opt. Phys. 1984, 30, 2193− 2209. (24) Piecuch, P.; Toboła, R.; Paldus, J. Approximate account of connected quadruply excited clusters in single-reference coupledcluster theory via cluster analysis of the projected unrestricted HartreeFock wave function. Phys. Rev. A: At., Mol., Opt. Phys. 1996, 54, 1210− 1241. (25) Li, X.; Paldus, J. Reduced multireference CCSD method: An effective approach to quasidegenerate states. J. Chem. Phys. 1997, 107, 6257. (26) Paldus, J.; Li, X. Externally corrected coupled-cluster approaches: Energy versus Amplitude corrected CCSD. Collect. Czech. Chem. Commun. 2003, 68, 554−586. (27) Kowalski, K.; Piecuch, P. The method of moments of coupledcluster equations and the renormalized CCSD[T], CCSD(T), CCSD(TQ), and CCSDT(Q) approaches. J. Chem. Phys. 2000, 113, 18. (28) Kinoshita, T.; Hino, O.; Bartlett, R. J. Coupled-cluster method tailored by configuration interaction. J. Chem. Phys. 2005, 123, 074106. (29) Veis, L.; Antalík, A.; Brabec, J.; Neese, F.; Legeza, Ö .; Pittner, J. Coupled Cluster Method with Single and Double Excitations Tailored by Matrix Product State Wave Functions. J. Phys. Chem. Lett. 2016, 7, 4072−4078. PMID: 27682626. (30) Piecuch, P.; Kowalski, K.; Pimienta, I. S. O.; McGuire, M. J. Recent advances in electronic structure theory: Method of moments of coupled-cluster equations and renormalized coupled-cluster approaches. Int. Rev. Phys. Chem. 2002, 21, 527. (31) Piecuch, P.; Kowalski, K.; Pimienta, I. S. O. Method of Moments of Coupled-Cluster Equations: Externally Corrected Approaches Employing Configuration Interaction Wave Functions. Int. J. Mol. Sci. 2002, 3, 475−497. (32) Kowalski, K.; Piecuch, P. Extension of the method of moments of coupled-cluster equations to excited states: The triples and

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jiri Brabec: 0000-0002-7764-9890 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Czech Science Foundation (Project 15-00058Y), DAAD (DAAD-16-17), and Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations (Project IT4Innovations National Supercomputing Center LM2015070).

■

REFERENCES

(1) Č ížek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256−4266. (2) Jeziorski, B.; Monkhorst, H. J. Coupled-cluster method for multideterminantal reference states. Phys. Rev. A: At., Mol., Opt. Phys. 1981, 24, 1668−1681. (3) Kucharski, S. A.; Bartlett, R. J. Hilbert space multireference coupled-cluster methods. I. The single and double excitation model. J. Chem. Phys. 1991, 95, 8227−8238. (4) Balková, A.; Kucharski, S. A.; Bartlett, R. J. The multi-reference Hilbert space coupled-cluster study of the Li2 molecule. Application in a complete model space. Chem. Phys. Lett. 1991, 182, 511−518. (5) Mahapatra, U. S.; Datta, B.; Mukherjee, D. A size-consistent statespecific multireference coupled cluster theory: Formal developments and molecular applications. J. Chem. Phys. 1999, 110, 6171−6188. (6) Másǐ k, J.; Hubač, I. Multireference Brillouin-Wigner coupledcluster theory. Single-root approach. Adv. Quantum Chem. 1998, 31, 75−104. (7) Hanrath, M. An exponential multireference wave-function Ansatz. J. Chem. Phys. 2005, 123, 084102. (8) Banerjee, A.; Simons, J. The coupled-cluster method with a multiconfiguration reference state. Int. J. Quantum Chem. 1981, 19, 207−216. (9) Hoffmann, M. R.; Simons, J. Analytical energy gradients for a unitary coupled-cluster theory. Chem. Phys. Lett. 1987, 142, 451−454. (10) Mukherjee, D. Normal ordering and a Wick-like reduction theorem for fermions with respect to a multi-determinantal reference state. Chem. Phys. Lett. 1997, 274, 561−566. (11) Evangelista, F. A.; Gauss, J. An orbital-invariant internally contracted multireference coupled cluster approach. J. Chem. Phys. 2011, 134, 114102. 1379


Article

Journal of Chemical Theory and Computation quadruples corrections to the equation-of-motion coupled-cluster singles and doubles energies. J. Chem. Phys. 2002, 116, 7411. (33) Włoch, M.; Lodriguito, M. D.; Piecuch, P.; Gour, J. R. Two new classes of non-iterative coupled-cluster methods derived from the method of moments of coupled-cluster equations. Mol. Phys. 2006, 104, 2991. (34) Lodriguito, M. D.; Kowalski, K.; Wloch, M.; Piecuch, P. Noniterative coupled-cluster methods employing multi-reference perturbation theory wave functions. J. Mol. Struct.: THEOCHEM 2006, 771, 89−104. (35) Piecuch, P.; Kowalski, K.; Pimienta, I. S. O.; Fan, P.-D.; Lodriguito, M.; McGuire, M. J.; Kucharski, S. A.; Kuś, T.; Musiał, M. Method of moments of coupled-cluster equations: a new formalism for designing accurate electronic structure methods for ground and excited states. Theor. Chem. Acc. 2004, 112, 349. (36) Kowalski, K.; Piecuch, P. New type of noniterative energy corrections for excited electronic states: extension of the method of moments of coupled-cluster equations to the equation-of-motion coupled-cluster formalism. J. Chem. Phys. 2001, 115, 2966−2978. (37) Piecuch, P.; Kowalski, K. In Computational Chemistry: Reviews of Current Trends; Leszczyński, J., Ed.; World Scientific: Singapore, 2011; Vol. 5; pp 1−104. (38) Mahapatra, U. S.; Datta, B.; Mukherjee, D. Development of a size-consistent state-specific multireference perturbation theory with relaxed model-space coefficients. Chem. Phys. Lett. 1999, 299, 42−50. (39) Evangelista, F. A.; Allen, W. D.; Schaefer, H. F., III High-order excitations in state-universal and state-specific multireference coupledcluster theories: Model systems. J. Chem. Phys. 2006, 125, 154113. (40) Evangelista, F. A.; Allen, W. D.; Schaefer, H. F., III Coupling term derivation and general implementation of state-specific multireference coupled cluster theories. J. Chem. Phys. 2007, 127, 024102. (41) Evangelista, F. A.; Simmonett, A. C.; Allen, W. D.; Schaefer, H. F., III; Gauss, J. Triple excitations in state-specific multireference coupled cluster theory. J. Chem. Phys. 2008, 128, 124104. (42) Das, S.; Mukherjee, D.; Kállay, M. Full implementation and benchmark studies of Mukherjee’s state-specific multireference coupled cluster ansatz. J. Chem. Phys. 2010, 132, 074103. (43) Bhaskaran-Nair, K.; Demel, O.; Pittner, J. Multireference StateSpecific Mukherjee’s coupled cluster method with non-iterative triexcitations. J. Chem. Phys. 2008, 129, 184105. (44) Bhaskaran-Nair, K.; Demel, O.; Šmydke, J.; Pittner, J. Multireference state-specific Mukherjee’s coupled cluster method with noniterative triexcitations using uncoupled approximation. J. Chem. Phys. 2011, 134, 154106. (45) Prochnow, E.; Evangelista, F. A.; Schaefer, H. F., III; Allen, W. D.; Gauss, J. Analytic gradients for the state-specific multireference coupled cluster singles and doubles model. J. Chem. Phys. 2009, 131, 064109. (46) Demel, O.; Kedzuch, S.; Svana, M.; Ten-no, S.; Pittner, J.; Noga, J. An explicitly correlated Mukherjee’s state specific coupled cluster method: development and pilot applications. Phys. Chem. Chem. Phys. 2012, 14, 4753−4762. (47) Das, S.; Bera, N.; Ghosh, S.; Mukherjee, D. An externallycorrected size-extensive single-root MRCC formalism: Its kinship with the rigorously size-extensive state-specific MRCC theory. J. Mol. Struct.: THEOCHEM 2006, 771, 79−87. Modelling Structure and Reactivity: the seventh triennial conference of the World Association of Theoritical and Computational Chemists (WATOC 2005). (48) Demel, O.; Bhaskaran-Nair, K.; Pittner, J. Uncoupled multireference state-specific Mukherjee’s coupled cluster method with triexcitations. J. Chem. Phys. 2010, 133, 134106. (49) Brabec, J.; Bhaskaran-Nair, K.; Kowalski, K.; Pittner, J.; van Dam, H. J. J. Towards large-scale calculations with State-Specific Multireference Coupled Cluster methods: Studies on dodecane, naphthynes, and polycarbenes. Chem. Phys. Lett. 2012, 542, 128−133. (50) Lotrich, V.; Flocke, N.; Ponton, M.; Yau, A. D.; Perera, A.; Deumens, E.; Bartlett, R. J. Parallel implementation of electronic structure energy, gradient, and Hessian calculations. J. Chem. Phys. 2008, 128, 194104.

(51) Deumens, E.; Lotrich, V. F.; Perera, A.; Ponton, M. J.; Sanders, B. A.; Bartlett, R. J. Software design of ACES III with the super instruction architecture. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 895−901. (52) Kowalski, K.; Krishnamoorthy, S.; Olson, R. M.; Tipparaju, V.; Aprà, E. Scalable Implementations of Accurate Excited-state Coupled Cluster Theories: Application of High-level Methods to Porphyrinbased Systems. Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis; New York, 2011; pp 72:1−72:10. (53) Aprà, E.; Klemm, M.; Kowalski, K. Efficient Implementation of Many-body Quantum Chemical Methods on the Intel Xeon Phi Coprocessor. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis; Piscataway, NJ, 2014; pp 674−684. (54) Brabec, J.; Krishnamoorthy, S.; van Dam, H.; Kowalski, K.; Pittner, J. Massively parallel implementation of the Brillouin-Wigner MRCCSD method. Chem. Phys. Lett. 2011, 514, 347−351. (55) See https://github.com/ValeevGroup/tiledarray for a massivelyparallel block-sparse tensor network written in C++. (Accessed May 13, 2017). (56) Solomonik, E.; Matthews, D.; Hammond, J. R.; Stanton, J. F.; Demmel, J. A massively parallel tensor contraction framework for coupled-cluster computations. Journal of Parallel and Distributed Computing 2014, 74, 3176−3190. Domain-Specific Languages and High-Level Frameworks for High-Performance Computing (57) Pulay, P. Localizability of dynamic electron correlation. Chem. Phys. Lett. 1983, 100, 151−154. (58) Sæbø, S.; Pulay, P. Local configuration interaction: An efficient approach for larger molecules. Chem. Phys. Lett. 1985, 113, 13−18. (59) Boys, S. F. Construction of Some Molecular Orbitals to Be Approximately Invariant for Changes from One Molecule to Another. Rev. Mod. Phys. 1960, 32, 296−299. (60) Pipek, J.; Mezey, P. G. A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys. 1989, 90, 4916−4926. (61) Hampel, C.; Werner, H.-J. Local treatment of electron correlation in coupled cluster theory. J. Chem. Phys. 1996, 104, 6286−6297. (62) Schütz, M.; Werner, H.-J. Low-order scaling local electron correlation methods. IV. Linear scaling local coupled-cluster (LCCSD). J. Chem. Phys. 2001, 114, 661−681. (63) Schutz, M. A new, fast, semi-direct implementation of linear scaling local coupled cluster theory. Phys. Chem. Chem. Phys. 2002, 4, 3941−3947. (64) Werner, H.-J.; Schütz, M. An efficient local coupled cluster method for accurate thermochemistry of large systems. J. Chem. Phys. 2011, 135, 144116. (65) Yang, J.; Chan, G. K.-L.; Manby, F. R.; Schütz, M.; Werner, H.-J. The orbital-specific-virtual local coupled cluster singles and doubles method. J. Chem. Phys. 2012, 136, 144105. (66) Fedorov, D. G.; Kitaura, K. Coupled-cluster theory based upon the fragment molecular-orbital method. J. Chem. Phys. 2005, 123, 134103. (67) Kobayashi, M.; Nakai, H. Extension of linear-scaling divide-andconquer-based correlation method to coupled cluster theory with singles and doubles excitations. J. Chem. Phys. 2008, 129, 044103. (68) Kristensen, K.; Ziółkowski, M.; Jansík, B.; Kjærgaard, T.; Jørgensen, P. A Locality Analysis of the Divide−Expand−Consolidate Coupled Cluster Amplitude Equations. J. Chem. Theory Comput. 2011, 7, 1677−1694. PMID: 26596432. (69) Høyvik, I.-M.; Kristensen, K.; Jansik, B.; Jørgensen, P. The divide-expand-consolidate family of coupled cluster methods: Numerical illustrations using second order Møller-Plesset perturbation theory. J. Chem. Phys. 2012, 136, 014105. (70) Stoll, H. The correlation energy of crystalline silicon. Chem. Phys. Lett. 1992, 191, 548−552. 1380


Article


Cluster with Pair Natural Orbitals (PNO-LCCSD). J. Chem. Theory Comput. 2017, 13, 3650−3675. PMID: 28661673. (92) Antony, J.; Grimme, S.; Liakos, D. G.; Neese, F. Protein−Ligand Interaction Energies with Dispersion Corrected Density Functional Theory and High-Level Wave Function Based Methods. J. Phys. Chem. A 2011, 115, 11210−11220. PMID: 21842894. (93) Anoop, A.; Thiel, W.; Neese, F. A Local Pair Natural Orbital Coupled Cluster Study of Rh Catalyzed Asymmetric Olefin Hydrogenation. J. Chem. Theory Comput. 2010, 6, 3137−3144. PMID: 26616776. (94) Liakos, D. G.; Neese, F. Interplay of Correlation and Relativistic Effects in Correlated Calculations on Transition-Metal Complexes: The (Cu2O2)2+ Core Revisited. J. Chem. Theory Comput. 2011, 7, 1511−1523. PMID: 26610142. (95) Zade, S. S.; Zamoshchik, N.; Reddy, A. R.; Fridman-Marueli, G.; Sheberla, D.; Bendikov, M. Products and Mechanism of Acene Dimerization. A Computational Study. J. Am. Chem. Soc. 2011, 133, 10803−10816. PMID: 21710966. (96) Kubas, A.; Brase, S.; Fink, K. Theoretical Approach Towards the Understanding of Asymmetric Additions of Dialkylzinc to Enals and Iminals Catalysed by [2.2]Paracyclophane-Based N,O-Ligands. Chem. Eur. J. 2012, 18, 8377−8385. (97) Ashtari, M.; Cann, N. Proline-based chiral stationary phases: A molecular dynamics study of the interfacial structure. Journal of Chromatography A 2011, 1218, 6331−6347. (98) Zhang, J.; Dolg, M. Dispersion Interaction Stabilizes Sterically Hindered Double Fullerenes. Chem. - Eur. J. 2014, 20, 13909−13912. (99) Minenkov, Y.; Chermak, E.; Cavallo, L. Accuracy of DLPNOCCSD(T) Method for Noncovalent Bond Dissociation Enthalpies from Coinage Metal Cation Complexes. J. Chem. Theory Comput. 2015, 11, 4664−4676. PMID: 26574257. (100) Sparta, M.; Neese, F. Chemical applications carried out by local pair natural orbital based coupled-cluster methods. Chem. Soc. Rev. 2014, 43, 5032−5041. (101) Liakos, D. G.; Sparta, M.; Kesharwani, M. K.; Martin, J. M. L.; Neese, F. Exploring the Accuracy Limits of Local Pair Natural Orbital Coupled-Cluster Theory. J. Chem. Theory Comput. 2015, 11, 1525− 1539. (102) Demel, O.; Pittner, J.; Neese, F. A Local Pair Natural OrbitalBased Multireference Mukherjee’s Coupled Cluster Method. J. Chem. Theory Comput. 2015, 11, 3104−3114. PMID: 26575747. (103) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. Introduction of n-electron valence states for multireference perturbation theory. J. Chem. Phys. 2001, 114, 10252−10264. (104) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. J. Chem. Phys. 2002, 117, 9138−9153. (105) Angeli, C.; Pastore, M.; Cimiraglia, R. New perspectives in multireference perturbation theory: the n-electron valence state approach. Theor. Chem. Acc. 2007, 117, 743−754. (106) Dunning, T. H. The Cl + Bk extrapolation method. Application to hydrogen fluoride11Work performed under the auspices of the Office of Basic Energy Sciences of the U.S. Department of Energy. Chem. Phys. 1979, 42, 249−258. (107) Pozun, Z. D.; Su, X.; Jordan, K. D. Establishing the Ground State of the Disjoint Diradical Tetramethyleneethane with Quantum Monte Carlo. J. Am. Chem. Soc. 2013, 135, 13862−13869. PMID: 23947763. (108) Arrowsmith, M.; Braunschweig, H.; Celik, M. A.; Dellermann, T.; Dewhurst, R. D.; Ewing, W. C.; Hammond, K.; Kramer, T.; Krummenacher, I.; Mies, J.; Radacki, K.; Schuster, J. K. Neutral zerovalent s-block complexes with strong multiple bonding. Nat. Chem. 2016, 8, 890−894. (109) Das, A.; Müller, T.; Plasser, F.; Lischka, H. Polyradical Character of Triangular Non-Kekulé Structures, Zethrenes, pQuinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study. J. Phys. Chem. A 2016, 120, 1625−1636. PMID: 26859789.

(71) Li, S.; Ma, J.; Jiang, Y. Linear scaling local correlation approach for solving the coupled cluster equations of large systems. J. Comput. Chem. 2002, 23, 237−244. (72) Edmiston, C.; Krauss, M. Configuration-Interaction Calculation of H3 and H2. J. Chem. Phys. 1965, 42, 1119−1120. (73) Ahlrichs, R.; Kutzelnigg, W. Ab initio calculations on small hydrides including electron correlation. Theoretica chimica acta 1968, 10, 377−387. (74) Meyer, W. Ionization energies of water from PNO-CI calculations. Int. J. Quantum Chem. 1971, 5, 341−348. (75) Meyer, W. PNO-CI and CEPA studies of electron correlation effects. Theoretica chimica acta 1974, 35, 277−292. (76) Taylor, P. R. A rapidly convergent CI expansion based on several reference configurations, using optimized correlating orbitals. J. Chem. Phys. 1981, 74, 1256−1270. (77) Fink, R.; Staemmler, V. A multi-configuration reference CEPA method based on pair natural orbitals. Theoretica chimica acta 1993, 87, 129−145. (78) Vahtras, O.; Almlöf, J.; Feyereisen, M. Integral approximations for LCAO-SCF calculations. Chem. Phys. Lett. 1993, 213, 514−518. (79) Neese, F.; Wennmohs, F.; Hansen, A. Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method. J. Chem. Phys. 2009, 130, 114108. (80) Neese, F.; Hansen, A.; Liakos, D. G. Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis. J. Chem. Phys. 2009, 131, 064103. (81) Riplinger, C.; Neese, F. An efficient and near linear scaling pair natural orbital based local coupled cluster method. J. Chem. Phys. 2013, 138, 034106. (82) Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. J. Chem. Phys. 2013, 139, 134101. (83) Saitow, M.; Becker, U.; Riplinger, C.; Valeev, E. F.; Neese, F. A new near-linear scaling, efficient and accurate, open-shell domainbased local pair natural orbital coupled cluster singles and doubles theory. J. Chem. Phys. 2017, 146, 164105. (84) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. J. Chem. Phys. 2015, 143, 034108. (85) Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse mapsA systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. (86) Guo, Y.; Sivalingam, K.; Valeev, E. F.; Neese, F. SparseMapsA systematic infrastructure for reduced-scaling electronic structure methods. III. Linear-scaling multireference domain-based pair natural orbital N-electron valence perturbation theory. J. Chem. Phys. 2016, 144, 094111. (87) Werner, H.-J.; Knizia, G.; Krause, C.; Schwilk, M.; Dornbach, M. Scalable Electron Correlation Methods I.: PNO-LMP2 with Linear Scaling in the Molecular Size and Near-Inverse-Linear Scaling in the Number of Processors. J. Chem. Theory Comput. 2015, 11, 484−507. PMID: 26580908. (88) Menezes, F.; Kats, D.; Werner, H.-J. Local complete active space second-order perturbation theory using pair natural orbitals (PNOCASPT2). J. Chem. Phys. 2016, 145, 124115. (89) Helmich, B.; Hättig, C. A pair natural orbital implementation of the coupled cluster model CC2 for excitation energies. J. Chem. Phys. 2013, 139, 084114. (90) Schmitz, G.; Helmich, B.; Hättig, C. A scaling PNO−MP2 method using a hybrid OSV−PNO approach with an iterative direct generation of OSVs. Mol. Phys. 2013, 111, 2463−2476. (91) Schwilk, M.; Ma, Q.; Köppl, C.; Werner, H.-J. Scalable Electron Correlation Methods. 3. Efficient and Accurate Parallel Local Coupled 1381


Article

Journal of Chemical Theory and Computation (110) Pavliček, N.; Mistry, A.; Majzik, Z.; Moll, N.; Meyer, G.; Fox, D. J.; Gross, L. Synthesis and characterization of triangulene. Nat. Nanotechnol. 2017, 12, 308−311.

1382


Domain-based Local Pair Natural Orbital version of Mukherjee's state

Recommend Documents