Externally corrected CCSD with renormalized perturbative triples (R-ecCCSD(T)) and density matrix renormalization group and selected configuration interaction external sources
Seunghoon Lee, Huanchen Zhai, Sandeep Sharma, Cyrus J. Umrigar, Garnet Kin-Lic Chan
EExternally corrected CCSD with renormalizedperturbative triples (R-ecCCSD(T)) and densitymatrix renormalization group and selectedconfiguration interaction external sources
Seunghoon Lee, ∗ , † Huanchen Zhai, † Sandeep Sharma, ∗ , ‡ C. J. Umrigar, ∗ , ¶ and Garnet Kin-Lic Chan ∗ , † † Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena,California 91125, USA ‡ Department of Chemistry, The University of Colorado at Boulder, Boulder, Colorado 80302,USA ¶ Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853, USA
E-mail: [email protected]; [email protected];[email protected]; [email protected]
Abstract
We investigate the renormalized perturbative triples correction together with the externallycorrected coupled-cluster singles and doubles (ecCCSD) method. We take the density matrixrenormalization group (DMRG) and heatbath CI (HCI) as external sources for the ecCCSDequations. The accuracy is assessed for the potential energy surfaces of H O, N , and F .We find that the triples correction significantly improves on ecCCSD and we do not seeany instability of the renormalized triples with respect to dissociation. We explore how to a r X i v : . [ phy s i c s . c h e m - ph ] F e b alance the cost of computing the external source amplitudes with respect to the accuracy ofthe subsequent CC calculation. In this context, we find that very approximate wavefunctions(and their large amplitudes) serve as an efficient and accurate external source. Finally, wecharacterize the domain of correlation treatable using the externally corrected method andrenormalized triples combination studied in this work via a well-known wavefunction diagnostic. Electronic structure methods that can efficiently handle both static and dynamic correlation remainan important area of investigation. Because there is a wide spectrum of strongly correlated problems,ranging from mildly "quasi"-degenerate scenarios (e.g. in the electronic structure of diradicals )to extensively near-degenerate problems (e.g. in the electronic structure of multi-centre transitionmetal clusters ) a variety of theoretical strategies have been proposed.For highly-degenerate problems, it is common to combine dynamic correlation methods withan explicit treatment of a multi-reference state. The representation of the multi-reference statecan range from an exact complete active space (CAS) representation (for small numbers oforbitals), to density matrix renormalization group (DMRG) (e.g. when almost all orbitals aredegenerate), to selected configuration interaction and Monte Carlo approximations forintermediate cases. On top of these, various flavours of perturbation theory, configurationinteraction, and exponential approximations have been explored. However, the combinationof dynamic correlation with multi-reference representations is not straightforward and usuallyleads to added conceptual, implementation, and computational complexity.For quasi-degenerate problems, an alternative strategy can be used, which incorporates a limitedamount of static correlation on top of an existing SR method. This has been particularly popularin conjunction with SR coupled cluster methods.
Some examples include variants of tailoredcoupled cluster and externally corrected coupled cluster methods.
The computational costof such SR static correlation methods is often lower than that of true MR dynamic correlationmethods and thus large active spaces become affordable. However, the approximations are limited2o problems with only a modest amount of degeneracy.In this work, we will focus on quasi-degenerate problems, and in particular, we will investigatethe externally corrected coupled cluster method.
This extracts static correlation from a MRmethod by using the MR wavefunction as an "external" source of higher order coupled clusteramplitudes. For example, in the ecCCSD approximation, the T and T amplitudes are extractedfrom the external source and a new set of T and T amplitudes are computed in their presence.Should the T and T amplitudes be exact, then the T and T amplitudes and the energy will beexact. A different, spiritually related, approximation is tailored CCSD, which has been ofrenewed interest of late. Here, instead of higher order cluster amplitudes, the large (activespace) T and T amplitudes are fixed from an external source.The ecCCSD method has a long history and external sources, ranging from unrestricted Hartree-Fock to CASSCF and CASCI and, most recently, full configuration interaction quantum MonteCarlo, have been used. One of the more successful applications of ecCCSD is the reducedmultireference (RMR) CCSD method, which uses a MRCISD wave function as the externalsource. RMR CCSD(T), which incorporates some of the residual dynamic correlation throughperturbative triples, has also been studied. Despite the promising performance of RMRCCSD(T)in several studies, it suffers from two main limitations. First, conventional MRCISD can onlybe applied for modest sizes of active spaces (typically, up to about orbitals as limited by theexact CAS treatment). Second, the (T) correction, although not divergent like its single-referencecounterpart, still overcorrects the dynamical correlation in the bond-stretched region. In this work, we make two modifications to ecCCSD to overcome and ameliorate the abovelimitations. First, we utilize variational DMRG and HCI wave functions as external sources forecCCSD. This allows for the use of larger quasidegenerate active spaces, of the size typicallytreated by DMRG and HCI. Second, we explore the renormalized perturbative triples correction.This has been shown in the single reference setting to ameliorate the overcorrection of standardperturbative triples, without affecting the computational scaling.We describe the use of DMRG and HCI as external sources for ecCCSD in Section 2.1. It3s possible to create a near-exact method by using a near-exact external source. Since near-exactexternal sources can be computed by DMRG and HCI for the small systems we employ as testcases in this paper, the critical question is not simply the accuracy of the method, but the balancebetween cost and accuracy. To this end, we explore a variety of approximate treatments of theexternal DMRG and HCI sources as discussed in Sec. 2.2. The various triples approximations forecCCSD are discussed in Section 2.3. Computational details are provided in Section 3 and theaccuracy of the renormalized perturbative triples correction is assessed for three potential energysurfaces (PESs) in Sections 4.1-4.4. We characterize the range of quasi-degenerate correlationscaptured in this work in Section 4.6. Finally we discuss the limitations of this method in Section4.7 and summarize our findings in Section 5.
In this work, we use a reference configuration | Φ (cid:105) with the same occupancy as the Hartree-Fock(HF) determinant. It is useful to define a projection operator onto the space of k -tuply excitedconfigurations relative to the reference; we denote this Q k . The external source is used to providean important subset of the triply and quadruply excited configurations; the projector onto thissubset is denoted Q eck , and its complement is Q ck . Thus Q k = Q ec k + Q c k , k = 3 , , (1) In ecCCSD, the coupled cluster operator T is given by T = T + T + Q ec3 T + Q ec4 T . (2)4ere, T n , n = 1 , , , , are the n -fold cluster operators. In this work, we extract Q ec k T k , k = 3 , from the DMRG or the HCI variational wave function. For the DMRG wave function, the triplyand quadruply excited configurations D p that define Q ec k , k = 3 , are chosen to be those where themagnitude of the CI coefficient c p is above a threshold, i.e. | c p | > s √ ω, (3)where s is an arbitrary scaling factor and ω is the largest discarded weight of the density matrix atthe maximum bond dimension in two-dot DMRG sweeps (carried out without noise). An efficientalgorithm to convert a matrix-product state to CI coefficients above a given threshold is described inAppendix A. For the HCI variational wave function, D p is included in the projectors Q ec k , k = 3 , using the heatbath algorithm with a threshold (cid:15) , i.e., it is included if |(cid:104) D p | H | D q (cid:105) c q | > (cid:15), (4)for at least one determinant D q which is already in the variational space. The extracted CIcoefficients are then converted into cluster amplitudes.The T and T amplitudes are obtained by solving the ecCCSD equations using fixed Q ec3 T and Q ec4 T , Q + Q )( H N e T + T + Q ec3 T + Q ec4 T ) C | Φ (cid:105) , (5)where H N is the Hamiltonian in normal-ordered form, and the subscript C denotes the connectedpart of the corresponding operator expression. With the relaxed T and T , the ecCCSD correlationenergy of the ground state is obtained as ∆ E ecCCSD0 = (cid:104) Φ | ( H N e T + T + Q ec3 T + Q ec4 T ) C | Φ (cid:105) . (6)5 .2 Approximations in the external source While DMRG and HCI can, in small molecules, be a source of nearly exact T and T amplitudeseven in the full orbital space, this provides no computational advantage as such calculations aremore expensive than the subsequent coupled cluster calculation. Consequently, it is importantto balance the cost of the external source calculation and that of the subsequent coupled clustercalculation by making approximations in the external source. This introduces the additional complicationthat one must ensure that errors introduced into the external source do not lead to unacceptableerrors in the final coupled cluster calculation.In this work, we consider six different types of approximate external sources with differentsizes of active spaces and different values of parameters summarized in Table 1. Type I usesCASSCF-like external sources in minimal active spaces. Since the minimal active space of F , i.e.,two electrons in two orbitals (2e,2o), does not contain T and T , we perform minimal active spaceecCC calculations for only H O and N (with (4e,4o) and (6e,6o), respectively). These provideamplitudes that are very close to the exact CASSCF amplitudes. However, these amplitudes lackthe relaxation that comes from allowing excitations within a larger space of orbitals, and of coursethe amplitudes outside the minimal active space are completely absent. Type III uses CASCI-likeexternal sources (the orbitals are not optimized to save computer time) in larger active spaces((8e,18o), (10e,16o), and (14e,16o) for H O, N , and F , respectively).In either case, one can potentially introduce bad external amplitudes if the effect of relaxationon the amplitude upon going to a larger space is large relative to the size of the amplitude (e.g.changes its sign). Thus, we study also Type II and Type IV external sources which are similar toType I and Type II external sources respectively except that they employ an additional thresholdto screen out all except the largest T and T amplitudes. The absolute values of the T and T elements at the most stretched geometry of each molecule are sorted in a single large vector. Thenorm of the vector is computed. Only the largest elements of the vector are retained such that theresulting norm is more than 80% of the norm of the full vector. Along PESs of each molecule, weused the same set of elements of T and T (but with the appropriate values for each geometry) as6he external sources, to maintain the smoothness of the PESs.As discussed in Section 4.4, the type-III and type-IV sources improve upon the PESs obtainedfrom the type-I and type-II sources, but the DMRG calculations to obtain the sources incur ahigher computational cost than the subsequent CC calculations. To reduce the cost, we havealso tried type-V sources, which employ loosely converged DMRG wave functions with smallbond dimensions of M = 25 , , and . Finally, type-VI sources employ large thresholds (cid:15) = 0 . , and . to obtain loosely converged HCI wave functions in the full orbital spaces.This combination has the advantage that it can be considered a black-box method wherein a singleparameter (cid:15) controls the tradeoff between accuracy and cost.Table 1: Six types of approximate external sources. Detailed explanation to be found in the maintext. Type Method Active space Cut-off parametersWave function Cluster amplitudeI DMRG Minimal a M = 2000 100 % with s = 0 . HCI (cid:15) = 10 − %II DMRG Minimal a M = 2000 80 % with s = 0 . HCI (cid:15) = 10 − %III DMRG Larger b M = 2000 100 % with s = 0 . IV DMRG Larger b M = 2000 80 % with s = 0 . V DMRG Larger b M = 25 , ,
100 80 % for H O and N , % for F with s = 0 . VI HCI Full c (cid:15) = 0 . , .
003 100 % a (4e , O , (6e , b (8e , O , (10e , , (14e , c (10e , O , (14e , , (14e , .3 Perturbative triples corrections Expressions for the standard, renormalized, and completely renormalized perturbative triples correctionscan be written down in analogy with their single reference definitions.
We first define thecompletely renormalized (CR)-ecCCSD(T) correction. We use the state | Ψ ecCCSD(T) (cid:105) defined as | Ψ ecCCSD(T) (cid:105) = (1 + T + T + Q ec3 T + Q c3 T [2]3 + Q c3 Z ) | Φ (cid:105) , (7) T [2]3 | Φ (cid:105) = R (3)0 ( V N T ) C | Φ (cid:105) , (8) Z | Φ (cid:105) = R (3)0 V N T | Φ (cid:105) . (9)where V N is the two-body part of the Hamiltonian in normal-ordered form and R (3)0 denotes thethree-body component of the reduced resolvent operator in many-body perturbation theory, givenby differences of orbital energies in the denominator. The resulting formula for the CR-ecCCSD(T)energy correction is δ CR − ecCCSD(T)0 = (cid:104) Ψ ecCCSD(T) | Q c3 ( H N e T ) C | Φ (cid:105)(cid:104) Ψ ecCCSD(T) | e T | Φ (cid:105) . (10)The energy corrections for renormalized (R)-ecCCSD(T) and ecCCSD(T) can be obtained bytaking the lowest-order estimates of the correction and by assuming the denominator to be one,i.e., δ R − ecCCSD(T)0 = (cid:104) Ψ ecCCSD(T) | Q c3 ( V N T ) C | Φ (cid:105)(cid:104) Ψ ecCCSD(T) | e T | Φ (cid:105) , (11) δ ecCCSD(T)0 = (cid:104) Ψ ecCCSD(T) | Q c3 ( V N T ) C | Φ (cid:105) . (12)Unlike perturbative triples without external correction (as in CCSD(T)), the approximate perturbativeexpression T [2]3 is only evaluated for determinants omitted in the external source. Thus it is notexpected to diverge as long as the external source includes all degeneracies. Nonetheless, it canstill overestimate the triples correlation. The role of the denominator in the "renormalized" triples8pproximations is to rescale this correction, which can be expected to reduce the overestimation. All CC calculations were performed using cc-pVTZ basis sets.
The ecCC calculations wereperformed with the six different types of external sources summarized in Table 1 and discussed inSection 2.2. The CASSCF-like external sources (type I and II) used natural orbitals in the minimalactive spaces and CASSCF orbitals for the core and external spaces, while the CASCI-like externalsources (type III-VI) used HF orbitals.All CC calculations were carried out using a local version of PySCF interfaced with StackBLOCK for DMRG and Arrow for HCI. We used Dice to get accurate SHCI PESs.
The dissociation PESs of H O and N , shown in Fig. 1, were obtained by the ecCC methods usingthe type-I external sources in Table 1 (i.e. near exact wave functions in the minimal active spacesand all amplitudes of T and T ). As expected, the ecCC curves with tightly converged HCI andDMRG external sources are almost identical. Thus, we show the PESs of ecCC using only oneof these external sources (colored solid lines) and the PESs of CC (colored dotted lines) in Fig.1. These are compared against accurate PESs represented as black lines obtained by SHCI in thefull space. The accurate energies from SHCI are given in the Supporting Information. The meanabsolute errors (MAE) and the non-parallelity errors (NPE) of the PESs are listed in Table 2.The CCSD curves (blue dotted lines) have an unphysical dip due to an inadequate treatmentof static correlation. This problem is completely eliminated within ecCCSD (blue solid line).However, there are significant errors with respect to the black curve at large distances, giving a(MAE,NPE) of ( . , . ) m E H for H O and ( . , . ) m E H for N .9e see that the (T) correction captures much of the missing dynamic correlation, such thatthe green ecCCSD(T) curves reach an accuracy of 1 mE H around the equilibrium geometry.However, when the bond is stretched, the (T) correction overestimates the dynamic correlation,leading to another unphysical dip in the PESs. Although this overestimation can easily be reducedby increasing the size of the active space for the small systems treated here, this would be veryexpensive for large systems. Alternatively, we can use the renormalized triples formula to damp the(T) correction in Eq. (11). R-ecCCSD(T) (the solid red curves) completely removes the unphysicaldips in the PESs. These attain (MAE,NPE) = ( . , . ) and ( . , . ) m E H for H O and N ,respectively. The difference between the type-I and type-II external sources is that the former use all theamplitudes of T and T while the latter use only a small number of the largest amplitudes, whichcontribute about 80% of the total weight of T and T (see Table 1 and Section 2.2). For theminimal active space of H O, i.e. (4e,4o), the external source contains only one non-zero T amplitude, and all elements of T are zero to within numerical noise. Thus, the PESs of ecCC using80% of the external amplitudes (type-II) are almost identical to those using 100% of the externalamplitudes (type-I). On the other hand, for the minimal active space of N , i.e. (6e,6o), the externalsource contains several large T and T amplitudes. At the stretched geometry corresponding toa bond length of a , 80% of the total T and T amplitude weight is recovered by the ninelargest elements of T . Figure 2 shows that the type-II external sources, although using feweramplitudes, improve the PESs of ecCCSD and R-ecCCSD(T) relative to the type-I sources. ThePES of R-ecCCSD(T) displays a (MAE,NPE) = ( . , . ) m E H .One concern with partial use of the amplitudes in the minimal active space is the possibilityof divergence in (T). Assuming a HF occupation of the natural orbitals, the recomputed orbitalenergies (diagonal parts of the Fock matrix) of the holes and those of the particles are such thatthe system retains a sizable gap. N at the . a bond length, for example, has a gap of around10igure 1: PESs of H O and N in the top and bottom panels, respectively, obtained with the CCmethods and the ecCC methods using the type-I external sources. The type-I sources correspondto near exact wave functions from the minimal active space, and use all the amplitudes of T and T . The blue, green, and red solid lines are the ecCCSD, ecCCSD(T), and R-ecCCSD(T) PESs,respectively. These are to be compared with the PESs of CC represented as dotted lines. The blacklines are accurate PESs obtained by SHCI. Bond length ( Å ) -76.35-76.25-76.15-76.05-75.95-75.85-75.75 T o t a l e n e r g y ( E H ) Type-I (4e,4o)ecCCSDecCCSD(T)R-ecCCSD(T) accurateCCSDCCSD(T)R-CCSD(T)
Bond length ( a ) -109.41-109.31-109.21-109.11-109.01-108.91-108.81-108.71 T o t a l e n e r g y ( E H ) Type-I (6e,6o)ecCCSDecCCSD(T)R-ecCCSD(T) accurateCCSDCCSD(T)R-CCSD(T) H ON . Hartrees between the three occupied orbitals and the three unoccupied orbitals. This preventsthe divergence of the (T) and the renormalized (T) corrections at the bond lengths we considered(although (T) largely overestimates dynamic correlation).Figure 2: PESs of N obtained by the ecCC methods using all the external amplitudes (type-I) andsome of the largest external amplitudes (type-II) from the minimal active space (solid and dottedlines, respectively). Other descriptions are the same as in Fig.1 Bond length ( a ) -109.41-109.31-109.21-109.11-109.01-108.91-108.81-108.71 T o t a l e n e r g y ( E H ) Type-I (6e,6o)ecCCSDecCCSD(T)R-ecCCSD(T)accurate Type-II (6e,6o)ecCCSDecCCSD(T)R-ecCCSD(T) N In addition, we investigated PESs of R-ecCCSD(T) using larger active spaces with all amplitudes(type-III) and 80% of the amplitudes (type-IV). We used active spaces of (8e,18o), (10e,16o), and(14e,16o) for H O, N , and F , respectively, and obtained near exact wave functions in the spaces.The PESs of ecCC using the type-III and type-IV external sources are shown as colored solid anddotted lines, respectively, in Figure 3. The MAE and NPE of the PESs are given in Table 2.For H O and N , the resulting PESs of R-ecCCSD(T) with all amplitudes (red solid lines inthe top and middle panels) achieve (MAE,NPE) = ( . , . ) and ( . , . ) m E H , respectively.12hese are minor improvements compared to the results using the minimal active space externalsources. However, the PESs obtained using only 80% of the external amplitudes (red dotted lines)are much closer to the black, accurate PES, and reach a (MAE,NPE) = ( . , . ) and ( . , . )m E H . Similarly to in the minimal active space, 80% of the amplitude weight in the larger activespace corresponds to only one element of T for H O and the nine largest elements of T for N .When we can find such large elements, the partial usage of the amplitudes clearly has advantagesin both accuracy and efficiency. It furthermore points to the importance of accounting for the errorfrom missing relaxation in the external source amplitudes.Unlike in H O and N , there are no particularly large elements in the T and T amplitudesof F . To truncate the amplitudes to 80% of their weight, we used the approximately 400 largestelements of T and T summing to 80% of the total T and T weights at the bond length of . Å. The red solid and dotted lines in the bottom panel of Figure 3 show the PESs of R-ecCCSD(T)using all and 80% of the amplitudes, respectively. These two PESs are very close to the accurateblack curve and reach a (MAE,NPE) = ( . , . ) and ( . , . ) m E H for 100% and 80% of theamplitudes, respectively. Although both are accurate, the partial use of the external source in thiscase (where there is not a small number of large elements) leads to slightly worse accuracy. In the previous section, we showed that the use of larger active spaces significantly improves theecCC PESs. However, obtaining tightly converged DMRG wave functions in large active spacesrequires more CPU time than the subsequent CC calculation. For example, optimizing an externalwave function with 16 orbitals requires around a few minutes of CPU time for one DMRG sweepwith M = 2000 . Although a few minutes is not prohibitively large in many applications, it is largecompared to the subsequent ecCC calculation which only takes tens of seconds, at least for thesmall molecules considered here. In addition, the fact that in some cases only a partial use of theamplitudes led to better results in the last section suggests that it is not a good use of computationaltime to tightly converge the external source. 13igure 3: PESs of H O, N , and F in the top, middle, and bottom panels, respectively, for theecCC methods using the larger than minimal active spaces and the near exact external sources.The solid and dotted lines correspond to PESs obtained using all external amplitudes of T and T (type-III) and only the largest amplitudes of T (type-IV), respectively. The descriptions areotherwise the same as those in Fig.2. Bond length ( Å ) -76.35-76.25-76.15-76.05-75.95-75.85-75.75 T o t a l e n e r g y ( E H ) Type-III (8e,18o)ecCCSDecCCSD(T)R-ecCCSD(T)accurate Type-IV (8e,18o)ecCCSDecCCSD(T)R-ecCCSD(T)
Bond length ( a ) -109.41-109.31-109.21-109.11-109.01-108.91-108.81-108.71 T o t a l e n e r g y ( E H ) Type-III (10e,16o)ecCCSDecCCSD(T)R-ecCCSD(T)accurate Type-IV (10e,16o)ecCCSDecCCSD(T)R-ecCCSD(T)
Bond length ( Å ) -199.30-199.28-199.26-199.24-199.22-199.20-199.18-199.16-199.14 T o t a l e n e r g y ( E H ) Type-III (14e,16o)ecCCSDecCCSD(T)R-ecCCSD(T)accurate Type-IV (14e,16o)ecCCSDecCCSD(T)R-ecCCSD(T) F H ON E H ) of H O, N , and F PESs in a range of geometries R ∈ [0 . , . Å, R ∈ [1 . , . a , and R ∈ [1 . , . Å, respectively. ecCC methods using the near exactexternal sources with all elements of T and T amplitudes ("100%") or only the largest elementsof the T amplitudes ("80%") from the external source.Method T , T H O N F Active MAE NPE Active MAE NPE Active MAE NPECCSD 15.9 30.4 32.4 102.3 34.8 57.7ecCCSD 100% 4e,4o 24.0 65.9 6e,6o 64.1 149.9 2e,2o a a80% 24.0 65.9 57.5 114.5100% 8e,18o 17.6 49.9 10e,16o 48.4 113.0 14e,16o 25.2 41.480% 15.6 39.4 41.4 74.0 27.0 40.9R-CCSD(T) 6.8 43.0 32.1 162.1 6.2 10.6R-ecCCSD(T) 100% 4e,4o 3.9 23.1 6e,6o 12.8 41.4 2e,2o a a80% 3.9 23.1 5.4 13.9100% 8e,18o 2.9 18.0 10e,16o 12.1 44.3 14e,16o 0.9 1.480% 1.1 8.1 4.6 8.4 1.7 1.6 a no T and T in the minimal active space of F We thus now consider loosely converged DMRG sources (type-V) in the larger active spaces.We re-computed PESs of R-ecCCSD(T) shown in Figure 3 using a DMRG source with bonddimensions 25, 50, and 100. We used the truncation to 80% amplitude weight for H O and N ,while we used all the external amplitudes for F . Table 3 shows the corresponding MAE and NPEin the same range of geometries in Table 2.For all cases, when we reduced the bond dimension to M = 100 , the MAE increased by . ∼ . m E H , and the NPE increased by . ∼ . m E H , compared to using M = 2000 .However for M = 100 , one DMRG sweep with 16 orbitals took only a few seconds of CPUtime at the 10 a bond length of N , giving a better computational balance between the DMRGcalculation and the subsequent ecCC calculations. When we further reduced the bond dimensionto M = 25 , the MAE increased by . m E H and the NPE increased by . and . m E H for H Oand N . In the case of F , which does not have a small number of large T and T elements, theMAE and NPE increased more, by . and . m E H .15 .5 PESs with type-VI external sources Type VI sources use the full orbital space and employ the single HCI parameter, (cid:15) , to select thelarge T and T amplitudes. They do not require a choice for the CAS space, and we did not furtherthreshold the T and T amplitudes. However, for systems where the molecule has a low-spinground state, but the dissociated fragments have high-spin ground states, the HCI wave functionsat stretched geometries are strongly spin-contaminated when (cid:15) is large, even when time-reversalsymmetry is employed to reduce the spin contamination. (This type of spin-contamination isavoided in the small bond dimension DMRG calculations via the full use of spin symmetry).Although this has little effect on the HCI energy, it makes the amplitudes unsuitable for ecCCcalculations.In the second part of Table 3 the MAE and NPE from type-VI sources using (cid:15) = 0 . and . are shown. These are evaluated over a smaller range of geometries than those used in the first partof Table 3 for the reason mentioned above. Since N is a singlet that dissociates into atoms thatare quartets, it has particularly large errors. The errors improve upon going from (cid:15) = 0 . to (cid:15) = 0 . , in particular the PEC has an unphysical dip at (cid:15) = 0 . which disappears at (cid:15) = 0 . ,but the errors are still substantial. Similarly to how using all the amplitudes could lead to largererrors than only using some of the amplitudes in the previous sections, the errors incurred from thetype VI sources here emphasize that the quality of the external amplitudes cannot be judged solelyfrom the energy obtained by the external method. In this section, we present the errors of R-ecCCSD(T) for the systems in this work and analyzethem using the well known CC error diagnostic D , defined by the matrix 2-norm of the T amplitudes. The magnitude of D can be used to distinguish between the SR and MR characterof the different geometries on the PESs. Organic molecules are sometimes considered to have MRcharacter when D is larger than . . Figure 4 shows the absolute errors on a log scale versusthe D diagnostic. Each symbol represents a geometry on the PESs of one of the molecules, H O,16able 3: MAE and NPE (m E H ) of R-ecCCSD(T) PESs for H O, N , and F using the approximateDMRG and HCI sources with a small bond dimension ( M ) or large threshold ( (cid:15) ). Detailedexplanation to be found in the main text. M H O N F Active MAE a NPE a Active MAE a NPE a Active MAE a NPE a
25 (8e,18o) 3.8 8.7 (10e,16o) 5.5 11.1 (14e,16o) 4.3 6.850 3.4 8.2 5.2 9.6 3.5 3.4100 3.0 8.4 4.9 9.4 1.4 2.12000 2.9 8.1 4.6 8.4 0.9 1.4 (cid:15) H O N F Active MAE b NPE b Active MAE b NPE b Active MAE b NPE b . (10e,58o) 1.6 7.1 (14e,60o) 8.1 61.5 (14e,58o) 3.6 4.7 . a R ∈ [0 . , . Å, [1 . , . a , and [1 . , . Å for H O, N , and F PESs, respectively. b R ∈ [0 . , . Å, [1 . , . a , and [1 . , . Å for H O, N , and F PESs, respectively.Figure 4: Absolute errors of R-CCSD(T) and ecCC methods on a log scale plotted against the D diagnostic. Detailed explanation in the main text. D | E rr o r | ( E H ) R-CCSD(T)ecCCSDecCCSD(T)R-ecCCSD(T) , and F . Square symbols denote R-CCSD(T) and circle symbols denote the various externallycorrected theories, using the type-V and the type-VI external sources. For D ranging from to . , the absolute errors of R-ecCCSD(T) (red circles) are less than . E H and most of the errorsare smaller than those of R-CCSD(T) (black squares) and ecCCSD (blue circles). The absoluteerrors of R-ecCCSD(T) are less than those of ecCCSD(T) (green circles) in the MR region where D > . , while they are mostly greater than those of ecCCSD(T) in the range . < D < . .Overall, it is clear to see that R-ecCCSD(T) offers the most balanced treatment of errors across awide range of SR and MR character. However, for very weakly correlated systems, the original(T) correction is slightly more accurate than the renormalized (T) correction. Although the work above shows that it is possible to obtain quantitative accuracy across the fullpotential energy surface, at a reasonable computational cost, using the combination of externalsources and the renormalized (T) correction, SR ecCC approaches have a fundamental limitationwhen the CI coefficient of the reference configuration (the HF configuration in this work) in theexternal source is exactly or numerically zero. This leads to overflow due to the assumption ofintermediate normalization. The smallest reference coefficient value we encountered in this workwas . at the stretched a geometry of N . Although this coefficient is not numerically zero, itis difficult to converge the ecCCSD energy. We extracted initial amplitudes of T and T from theexternal source and then iteratively converged the ecCCSD energy using a damping parameter of . to update the amplitudes. (We did not use the direct inversion in the iterative subspace (DIIS)algorithm). At this geometry, the energy could be converged monotonically up to a threshold of − Hartrees, although the norm of the amplitudes could not be converged, and we simply usedamplitudes from the last iteration in the set of monotonically decreasing energies.18
Conclusion
In this work, we explored externally corrected coupled cluster with a renormalized triples correction(R-ecCCSD(T)) using DMRG and HCI and external sources. The critical question is how to bestbalance the accuracy and cost of computing the external source with the cost of the overall method.To this end, we considered multiple types of external sources: "exact" external sources, where theDMRG and HCI wavefunctions were tightly converged, and "approximate" external sources wherethey were loosely converged. We also considered both "full" usage of the T and T amplitudes,and "partial" usage, wherein we retained only the largest elements.For all systems considered here, we found that R-ecCCSD(T) can significantly improve on thedescription of the potential energy surface given either by the external source alone or CC alone.For example, the unphysical dips in the PES of H O and N in the bond-stretched region can becompletely eliminated. The use of approximate external sources, possibly with truncation to onlythe large T and T amplitudes, appears to be a practical way to balance the cost of the externalcalculation and the coupled cluster calculation in small molecules. Using the D diagnostic tocharacterize the different points on the potential energy surfaces, we find that R-ecCCSD(T) givesabsolute errors of less than 15 m E H in the range of D from . to . . In fact, the errors ofR-ecCCSD(T) are less than those of ecCCSD and R-CCSD(T) in almost all cases, except when D is very small, where the renormalized (T) correction appears to be slightly less accurate than thesimple (T) correction.There are several interesting questions remaining which lie beyond what we have consideredin this work. For example, while R-ecCCSD(T) appears quite stable up to large values of the D diagnostic, what is the largest amount of multi-reference character which can be handled? Here, thedifficulty in solving the CC equations, and the divergence of the amplitudes reflecting the problemsof intermediate normalization, cannot be ignored. In addition, in the realm of quasidegenerateproblems, we can ask whether other non-iterative corrections such as the "completely renormalized"triples and quadruples corrections, corresponding to CR-ecCCSD(T) and CR-ecCCSD(T,Q),would further improve on the present R-ecCCSD(T) method. Finally, the current approximation,19ith its modest computational requirements on the external source, is applicable to the same scaleof systems that can be handled by single reference coupled cluster methods. Thus understandingthe performance of this method in larger correlated systems is of interest. Note : As this work was being prepared for submission, we were made aware of a related recentsubmission to the arxiv, that also discusses selected configuration interaction as an externalsource and perturbative triples corrections in the context of externally corrected coupled clustermethods.
Acknowledgement
Work by SL, HZ and GKC was supported by the US National Science Foundation via AwardCHE-1655333. GKC is a Simons Investigator. SS was supported by the US National ScienceFoundation grant CHE-1800584 and the Sloan research fellowship. CJU was supported in part bythe AFOSR under grant FA9550-18-1-0095.
Supporting Information Available eferences (1) Abe, M. Diradicals. Chem. Rev. , , 7011–7088.(2) Nakano, M. Electronic Structure of Open-Shell Singlet Molecules: Diradical CharacterViewpoint. Top. Curr. Chem. , .(3) Kurashige, Y.; Chan, G. K.-L.; Yanai, T. Entangled quantum electronic wavefunctions ofthe Mn 4 CaO 5 cluster in photosystem II. Nature chemistry , , 660–666.(4) Sharma, S.; Sivalingam, K.; Neese, F.; Chan, G. K.-L. Low-energy spectrum of iron–sulfurclusters directly from many-particle quantum mechanics. Nat. Chem. , , 927–933.(5) Li, Z.; Guo, S.; Sun, Q.; Chan, G. K.-L. Electronic landscape of the P-cluster of nitrogenaseas revealed through many-electron quantum wavefunction simulations. Nat. Chem. , , 1026–1033.(6) Das, G.; Wahl, A. C. New Techniques for the Computation of MulticonfigurationSelf-Consistent Field (MCSCF) Wavefunctions. The Journal of Chemical Physics , ,1769–1775.(7) Werner, H.-J.; Meyer, W. A quadratically convergent multiconfiguration–self-consistentfield method with simultaneous optimization of orbitals and CI coefficients. The Journalof Chemical Physics , , 2342–2356.(8) Werner, H.-J.; Knowles, P. J. A second order multiconfiguration SCF procedure withoptimum convergence. The Journal of chemical physics , , 5053–5063.(9) White, S. R.; Noack, R. M. Real-space quantum renormalization groups. Physical reviewletters , , 3487.(10) White, S. R. Density matrix formulation for quantum renormalization groups. Physicalreview letters , , 2863. 2111) White, S. R. Density-matrix algorithms for quantum renormalization groups. PhysicalReview B , , 10345.(12) White, S. R.; Martin, R. L. Ab initio quantum chemistry using the density matrixrenormalization group. The Journal of chemical physics , , 4127–4130.(13) Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial costalgorithm: A study of the density matrix renormalization group. The Journal of chemicalphysics , , 4462–4476.(14) Legeza, Ö.; Röder, J.; Hess, B. Controlling the accuracy of the density-matrixrenormalization-group method: The dynamical block state selection approach. PhysicalReview B , , 125114.(15) Chan, G. K.-L. An algorithm for large scale density matrix renormalization groupcalculations. The Journal of chemical physics , , 3172–3178.(16) Legeza, Ö.; Noack, R.; Sólyom, J.; Tincani, L. Computational Many-Particle Physics ;Springer, 2008; Vol. 739.(17) Chan, G.; Zgid, D. The Density Matrix Renormalization Group in Quantum Chemistry,Annu. Rep.
Comput. Chem , , 149–162.(18) Marti, K. H.; Reiher, M. The density matrix renormalization group algorithm in quantumchemistry. Zeitschrift für Physikalische Chemie , , 583–599.(19) Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. Annual review of physical chemistry , , 465–481.(20) Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithmsfor quantum chemistry. The Journal of chemical physics , , 124121.(21) Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantumchemistry. The European Physical Journal D , , 1–20.2222) Szalay, S.; Pfeffer, M.; Murg, V.; Barcza, G.; Verstraete, F.; Schneider, R.; Legeza, Ö.Tensor product methods and entanglement optimization for ab initio quantum chemistry. International Journal of Quantum Chemistry , , 1342.(23) Yanai, T.; Kurashige, Y.; Mizukami, W.; Chalupsk`y, J.; Lan, T. N.; Saitow, M. Densitymatrix renormalization group for ab initio Calculations and associated dynamic correlationmethods: A review of theory and applications. International Journal of Quantum Chemistry , , 283–299.(24) Ivanic, J.; Ruedenberg, K. Identification of deadwood in configuration spaces throughgeneral direct configuration interaction. Theoretical Chemistry Accounts , ,339–351.(25) Huron, B.; Malrieu, J.; Rancurel, P. Iterative perturbation calculations of ground andexcited state energies from multiconfigurational zeroth-order wavefunctions. The Journalof Chemical Physics , , 5745–5759.(26) Garniron, Y.; Scemama, A.; Giner, E.; Caffarel, M.; Loos, P.-F. Selected configurationinteraction dressed by perturbation. J. Chem. Phys. , , 064103.(27) Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. AMountaineering Strategy to Excited States: Highly Accurate Reference Energies andBenchmarks. J. Chem. Theory Comput. , , 4360–4379.(28) Garniron, Y.; Scemama, A.; Loos, P.-F.; Caffarel, M. Hybrid stochastic-deterministiccalculation of the second-order perturbative contribution of multireference perturbationtheory. J. Chem. Phys. , .(29) Buenker, R. J.; Peyerimhoff, S. D. Individualized configuration selection in CI calculationswith subsequent energy extrapolation. Theoretica chimica acta , , 33–58.2330) Evangelisti, S.; Daudey, J.-P.; Malrieu, J.-P. Convergence of an improved CIPSI algorithm. Chemical Physics , , 91–102.(31) Harrison, R. J. Approximating full configuration interaction with selected configurationinteraction and perturbation theory. The Journal of chemical physics , , 5021–5031.(32) Steiner, M.; Wenzel, W.; Wilson, K.; Wilkins, J. The efficient treatment of higher excitationsin CI calculations: A comparison of exact and approximate results. Chemical Physics Letters , , 263 – 268.(33) Neese, F. A spectroscopy oriented configuration interaction procedure. The Journal ofChemical Physics , , 9428–9443.(34) Abrams, M. L.; Sherrill, C. D. Important configurations in configuration interaction andcoupled-cluster wave functions. Chemical Physics Letters , , 121 – 124.(35) Bytautas, L.; Ruedenberg, K. A priori identification of configurational deadwood. ChemicalPhysics , , 64 – 75, Moving Frontiers in Quantum Chemistry:.(36) Evangelista, F. A. A driven similarity renormalization group approach to quantummany-body problems. The Journal of chemical physics , , 054109.(37) Knowles, P. J. Compressive sampling in configuration interaction wavefunctions. MolecularPhysics , , 1655–1660.(38) Schriber, J. B.; Evangelista, F. A. Communication: An adaptive configuration interactionapproach for strongly correlated electrons with tunable accuracy. The Journal of ChemicalPhysics , , 161106.(39) Tubman, N. M.; Lee, J.; Takeshita, T. Y.; Head-Gordon, M.; Whaley, K. B. A deterministicalternative to the full configuration interaction quantum Monte Carlo method. The Journalof chemical physics , , 044112. 2440) Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward full CI. Journal of chemical theory andcomputation , , 1169–1178.(41) Caffarel, M.; Applencourt, T.; Giner, E.; Scemama, A. Recent Progress in Quantum MonteCarlo ; ACS Publications, 2016; pp 15–46.(42) Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-Bath Configuration Interaction: AnEfficient Selected Configuration Interaction Algorithm Inspired by Heat-Bath Sampling.
J.Chem. Theory Comput. , , 3674–3680.(43) Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. SemistochasticHeat-Bath Configuration Interaction Method: Selected Configuration Interaction withSemistochastic Perturbation Theory. Journal of Chemical Theory and Computation , , 1595–1604.(44) Li, J.; Otten, M.; Holmes, A. A.; Sharma, S.; Umrigar, C. J. Fast Semistochastic Heat-BathConfiguration Interaction. J. Chem. Phys. , , 214110.(45) Booth, G. H.; Thom, A. J.; Alavi, A. Fermion Monte Carlo without fixed nodes: A gameof life, death, and annihilation in Slater determinant space. The Journal of chemical physics , , 054106.(46) Cleland, D.; Booth, G. H.; Alavi, A. Communications: Survival of the fittest: Acceleratingconvergence in full configuration-interaction quantum Monte Carlo. 2010.(47) Ghanem, K.; Lozovoi, A. Y.; Alavi, A. Unbiasing the initiator approximation in fullconfiguration interaction quantum Monte Carlo. The Journal of chemical physics , ,224108.(48) Sabzevari, I.; Sharma, S. Improved speed and scaling in orbital space variational MonteCarlo. Journal of chemical theory and computation , , 6276–6286.2549) Neuscamman, E. Subtractive manufacturing with geminal powers: making good use of abad wave function. Molecular Physics , , 577–583.(50) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-orderperturbation theory with a CASSCF reference function. Journal of Physical Chemistry , , 5483–5488.(51) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. Introduction ofn-electron valence states for multireference perturbation theory. The Journal of ChemicalPhysics , , 10252–10264.(52) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. N-electron valence state perturbation theory: afast implementation of the strongly contracted variant. Chemical physics letters , ,297–305.(53) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory: Aspinless formulation and an efficient implementation of the strongly contracted and of thepartially contracted variants. The Journal of chemical physics , , 9138–9153.(54) Neuscamman, E.; Yanai, T.; Chan, G. K.-L. A review of canonical transformation theory. International Reviews in Physical Chemistry , , 231–271.(55) Yanai, T.; Kurashige, Y.; Neuscamman, E.; Chan, G. K.-L. Multireference quantumchemistry through a joint density matrix renormalization group and canonicaltransformation theory. The Journal of chemical physics , , 024105.(56) Kurashige, Y.; Yanai, T. Second-order perturbation theory with a density matrixrenormalization group self-consistent field reference function: Theory and application tothe study of chromium dimer. The Journal of chemical physics , , 094104.(57) Sharma, S.; Chan, G. K.-L. Communication: A flexible multi-reference perturbation theoryby minimizing the Hylleraas functional with matrix product states. 2014.2658) Sharma, S.; Alavi, A. Multireference linearized coupled cluster theory for stronglycorrelated systems using matrix product states. The Journal of chemical physics , ,102815.(59) Guo, S.; Watson, M. A.; Hu, W.; Sun, Q.; Chan, G. K.-L. N-electron valence stateperturbation theory based on a density matrix renormalization group reference function,with applications to the chromium dimer and a trimer model of poly (p-phenylenevinylene). Journal of chemical theory and computation , , 1583–1591.(60) Sokolov, A. Y.; Guo, S.; Ronca, E.; Chan, G. K.-L. Time-dependent N-electron valenceperturbation theory with matrix product state reference wavefunctions for large active spacesand basis sets: Applications to the chromium dimer and all-trans polyenes. The Journal ofchemical physics , , 244102.(61) Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochasticheat-bath configuration interaction method: Selected configuration interaction withsemistochastic perturbation theory. Journal of chemical theory and computation , ,1595–1604.(62) Guo, S.; Li, Z.; Chan, G. K.-L. Communication: An efficient stochastic algorithm for theperturbative density matrix renormalization group in large active spaces. The Journal ofchemical physics , , 221104.(63) Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. New implementation of the graphicalunitary group approach for multireference direct configuration interaction calculations. International Journal of Quantum Chemistry , , 91–100.(64) Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfigurationself-consistent field and multireference configuration interaction methods and applications. Chemical reviews , , 108–181. 2765) Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; Müller, T.; Szalay, P. G.;Seth, M.; Kedziora, G. S.; Yabushita, S., et al. High-level multireference methods in thequantum-chemistry program system COLUMBUS: Analytic MR-CISD and MR-AQCCgradients and MR-AQCC-LRT for excited states, GUGA spin–orbit CI and parallel CIdensity. Physical Chemistry Chemical Physics , , 664–673.(66) Saitow, M.; Kurashige, Y.; Yanai, T. Fully internally contracted multireference configurationinteraction theory using density matrix renormalization group: A reduced-scalingimplementation derived by computer-aided tensor factorization. Journal of chemical theoryand computation , , 5120–5131.(67) Paldus, J.; Li, X. Adv. Chem. Phys. , 1–175.(68) Kong, L.; Shamasundar, K.; Demel, O.; Nooijen, M. State specific equation of motioncoupled cluster method in general active space. The Journal of chemical physics , ,114101.(69) Jeziorski, B. Multireference coupled-cluster Ansatz. Molecular Physics , ,3043–3054.(70) ˇCížek, J. On the correlation problem in atomic and molecular systems. Calculationof wavefunction components in Ursell-type expansion using quantum-field theoreticalmethods. The Journal of Chemical Physics , , 4256–4266.(71) ˇCížek, J. On the use of the cluster expansion and the technique of diagrams in calculationsof correlation effects in atoms and molecules. Advances in chemical physics , 35–89.(72) ˇCížek, J.; Paldus, J. Correlation problems in atomic and molecular systems III. Rederivationof the coupled-pair many-electron theory using the traditional quantum chemical methodst.
International Journal of Quantum Chemistry , , 359–379.2873) Paldus, J.; ˇCížek, J.; Shavitt, I. Correlation Problems in Atomic and Molecular Systems. IV.Extended Coupled-Pair Many-Electron Theory and Its Application to the B H 3 Molecule. Physical Review A , , 50.(74) Bartlett, R. J. In Modern Electronic Structure Theory ; Yarkony, D. R., Ed.; World Scientific,1995; Vol. 1; pp 1047–1131.(75) Shavitt, I.; Bartlett, R. J.
Many-body methods in chemistry and physics: MBPT andcoupled-cluster theory ; Cambridge university press, 2009.(76) Kinoshita, T.; Hino, O.; Bartlett, R. J. Coupled-cluster method tailored by configurationinteraction.
The Journal of chemical physics , , 074106.(77) Hino, O.; Kinoshita, T.; Chan, G. K.-L.; Bartlett, R. J. Tailored coupled cluster singles anddoubles method applied to calculations on molecular structure and harmonic vibrationalfrequencies of ozone. The Journal of chemical physics , , 114311.(78) Lyakh, D. I.; Lotrich, V. F.; Bartlett, R. J. The ‘tailored’CCSD (T) description of theautomerization of cyclobutadiene. Chemical Physics Letters , , 166–171.(79) Veis, L.; Antalík, A.; Brabec, J.; Neese, F.; Legeza, O.; Pittner, J. Coupled cluster methodwith single and double excitations tailored by matrix product state wave functions. Thejournal of physical chemistry letters , , 4072–4078.(80) Faulstich, F. M.; Máté, M.; Laestadius, A.; Csirik, M. A.; Veis, L.; Antalik, A.; Brabec, J.;Schneider, R.; Pittner, J.; Kvaal, S., et al. Numerical and theoretical aspects of theDMRG-TCC method exemplified by the nitrogen dimer. Journal of chemical theory andcomputation , , 2206–2220.(81) Vitale, E.; Alavi, A.; Kats, D. FCIQMC-tailored distinguishable cluster approach. Journalof Chemical Theory and Computation , , 5621–5634.2982) Paldus, J.; Planelles, J. Valence bond corrected single reference coupled cluster approach. Theoretica chimica acta , , 13–31.(83) Planelles, J.; Paldus, J.; Li, X. Valence bond corrected single reference coupled clusterapproach. Theoretica chimica acta , , 33.(84) Planelles, J.; Paldus, J.; Li, X. Valence bond corrected single reference coupled clusterapproach. Theoretica chimica acta , , 59.(85) Paldus, J. Externally and internally corrected coupled cluster approaches: an overview. Journal of Mathematical Chemistry , , 477–502.(86) Paldus, J.; Boyle, M. Cluster analysis of the full configuration interaction wave functions ofcyclic polyene models. International Journal of Quantum Chemistry , , 1281–1305.(87) Paldus, J.; Chin, E.; Grey, M. Bond length alternation in cyclic polyenes. II. Unrestrictedhartree–fock method. International journal of quantum chemistry , , 395–409.(88) Paldus, J.; ˇCížek, J.; Takahashi, M. Approximate account of the connected quadruply excitedclusters in the coupled-pair many-electron theory. Physical Review A , , 2193.(89) Piecuch, P.; Tobol, R.; Paldus, J. Approximate account of connected quadruply excitedclusters in single-reference coupled-cluster theory via cluster analysis of the projectedunrestricted Hartree-Fock wave function. Physical Review A , , 1210.(90) Peris, G.; Planelles, J.; Paldus, J. Single-reference CCSD approach employing three-andfour-body CAS SCF corrections: A preliminary study of a simple model. Internationaljournal of quantum chemistry , , 137–151.(91) Li, X.; Peris, G.; Planelles, J.; Rajadall, F.; Paldus, J. Externally corrected singles anddoubles coupled cluster methods for open-shell systems. The Journal of chemical physics , , 90–98. 3092) Peris, G.; Rajadell, F.; Li, X.; Planelles, J.; Paldus, J. Externally corrected singles anddoubles coupled cluster methods for open-shell systems. II. Applications to the low lyingdoublet states of OH, NH2, CH3 and CN radicals. Molecular Physics , , 235–248.(93) Stolarczyk, L. Z. Complete active space coupled-cluster method. Extension ofsingle-reference coupled-cluster method using the CASSCF wavefunction. Chemicalphysics letters , , 1–6.(94) Xu, E.; Li, S. The externally corrected coupled cluster approach with four-and five-bodyclusters from the CASSCF wave function. The Journal of chemical physics , ,094119.(95) Deustua, J. E.; Magoulas, I.; Shen, J.; Piecuch, P. Communication: Approaching exactquantum chemistry by cluster analysis of full configuration interaction quantum MonteCarlo wave functions. The Journal of Chemical Physics , , 151101.(96) Li, X.; Paldus, J. Reduced multireference CCSD method: An effective approach toquasidegenerate states. The Journal of Chemical Physics , , 6257–6269.(97) Li, X.; Paldus, J. Truncated version of the reduced multireference coupled-cluster methodwith perturbation selection of higher than pair clusters. International Journal of QuantumChemistry , , 743–756.(98) Li, X.; Paldus, J. Reduced multireference coupled cluster method with singles and doubles:Perturbative corrections for triples. The Journal of Chemical Physics , , 174101.(99) Li, X.; Paldus, J. Dissociation of N2 triple bond: a reduced multireference CCSD study. Chemical physics letters , , 145–154.(100) Li, X.; Paldus, J. Singlet-triplet splitting in methylene: An accurate description of dynamicand nondynamic correlation by reduced multireference coupled cluster method. Collectionof Czechoslovak chemical communications , , 1381–1393.31101) Li, X.; Paldus, J. Simultaneous handling of dynamical and nondynamical correlation viareduced multireference coupled cluster method: Geometry and harmonic force field ofozone. The Journal of chemical physics , , 2844–2852.(102) Li, X.; Paldus, J. Reduced multireference coupled cluster method: Ro-vibrational spectra ofN 2. The Journal of Chemical Physics , , 9966–9977.(103) Li, X.; Paldus, J. Effect of spin contamination on the prediction of barrier heights bycoupled-cluster theory: F+ H2 → HF+ H reaction.
International Journal of QuantumChemistry , , 281–290.(104) Li, X.; Paldus, J. A truncated version of reduced multireference coupled-cluster method withsingles and doubles and noniterative triples: Application to F 2 and Ni (CO) n (n= 1, 2, and4). The Journal of chemical physics , , 164107.(105) Li, X.; Paldus, J. Singlet–triplet separation in BN and C2: simple yet exceptional systemsfor advanced correlated methods. Chemical physics letters , , 179–184.(106) Li, X.; Paldus, J. Reduced multireference coupled-cluster method: Barrier heights forheavy atom transfer, nucleophilic substitution, association, and unimolecular reactions. TheJournal of Physical Chemistry A , , 11189–11197.(107) Li, X.; Paldus, J. Coupled-cluster approach to spontaneous symmetry breaking in molecules:The linear N3 radical. International Journal of Quantum Chemistry , , 2117–2127.(108) Li, X.; Paldus, J. Full potential energy curve for N2 by the reduced multireferencecoupled-cluster method. The Journal of Chemical Physics , , 054104.(109) Kowalski, K.; Piecuch, P. The method of moments of coupled-cluster equations and therenormalized CCSD[T], CCSD(T), CCSD(TQ), and CCSDT(Q) approaches. J. Chem. Phys. , , 18. 32110) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atomsboron through neon and hydrogen. J. Chem. Phys. , , 1007–1023.(111) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.;McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K. PySCF: thePython-based simulations of chemistry framework. 2017; https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1340 .(112) Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial costalgorithm: A study of the density matrix renormalization group. The Journal of ChemicalPhysics , , 4462–4476.(113) Chan, G. K.-L. An algorithm for large scale density matrix renormalization groupcalculations. The Journal of Chemical Physics , , 3172–3178.(114) Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.-L. Orbital optimization in the densitymatrix renormalization group, with applications to polyenes and β-carotene. The Journal ofChemical Physics , , 144117.(115) Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithmsfor quantum chemistry. The Journal of Chemical Physics , , 124121.(116) Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. SemistochasticHeat-Bath Configuration Interaction Method: Selected Configuration Interaction withSemistochastic Perturbation Theory. J. Chem. Theory Comput. , , 1595–1604.(117) Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-bath Configuration Interaction: Anefficient selected CI algorithm inspired by heat-bath sampling. J. Chem. Theory Comput. , , 3674–3680.(118) Smith, J. E. T.; Mussard, B.; Holmes, A. A.; Sharma, S. Cheap and Near Exact CASSCF33ith Large Active Spaces. Journal of Chemical Theory and Computation , ,5468–5478.(119) Nielsen, I. M. B.; Janssen, C. L. Double-substitution-based diagnostics for coupled-clusterand Møller–Plesset perturbation theory. Chemical physics letters , , 568–576.(120) Magoulas, I.; Gururangan, K.; Piecuch, P.; Deustua, J. E.; Shen, J. Is Externally CorrectedCoupled Cluster Always Better than the Underlying Truncated Configuration Interaction?2021. 34 ppendixA A sweep algorithm converting MPS to CI coefficients with athreshold In DMRG, the electronic wave function is represented by a matrix-product state (MPS) | Ψ (cid:105) = (cid:88) n ,...,n K (cid:88) α ,...,α K − A n α A n α α · · · A n K α K − | n n · · · n K (cid:105) , (A.1) { n } = { vac , ↑ , ↓ , ↑↓} , (A.2)where n i is the occupation of orbital i , | n n · · · n K (cid:105) is the occupation-number representation of adeterminant, and α i is an auxiliary index. Here, (cid:80) α i A n i α i − α i A n i +1 α i α i +1 denotes a matrix product andit is assumed that the dimensions of all the auxiliary indices are the same (the bond dimension M ).The wave function of the ground state can be optimized by the efficient DMRG sweep algorithm.In order to get Q ec3 T and Q ec4 T in Eq. (2), we extract quadruple- and lower-order CI excitationamplitudes from the MPS by a sweep algorithm. To avoid repeated or unnecessary computation,we here describe how to obtain CI coefficients whose values are larger than a threshold in Eq. (3),with a concomitant reduction in computational cost and memory usage from a naive approach.We first start with the MPS in left canonical form. During a sweep to compute the amplitude,at any given point (e.g. at some site p ) one has a set of partial coefficients c α p ( n n . . . n p ) = (cid:80) α ...α p − A n α . . . A n p α p − α p . If (cid:80) α p | c α p ( n n . . . n p ) | < thresh, then this partial coefficient isdropped, as are all determinants involving the occupancy string | n n . . . n p (cid:105) . This is because ifthe MPS is in left canonical form, the above condition on the partial coefficient guarantees that thecoefficient of any determinant generated by the MPS which contains | n n . . . n p (cid:105) as a substringis also less than the threshold in magnitude. In addition, since the orbitals are associated withdefinite hole or particle character, we also drop any coefficient associated with more than 4 holesor 4 particles. Finally, in this process, we can take advantage of the conserved quantum numbers to35nly generate symmetry unique partial coefficients (e.g. if S z = 0 , then the values of c α p +1+1