Near-linear Scaling in DMRG-based Tailored Coupled Clusters: An Implementation of DLPNO-TCCSD and DLPNO-TCCSD(T)

TL;DR

This paper introduces a new, more accurate and scalable implementation of the DMRG-based tailored coupled clusters method (TCCSD) using the DLPNO approach, achieving near-linear scaling and high accuracy in complex systems.

Contribution

The paper presents a novel DLPNO-TCCSD implementation that improves accuracy, scaling, and consistency over previous LPNO versions, including the addition of perturbative triples correction.

Findings

Achieved 99.8-99.9% correlation energy retrieval in tested systems.

Obtained energy differences within 1 kcal/mol of chemical accuracy.

Demonstrated near-linear scaling in the new implementation.

Abstract

We present a new implementation of DMRG-based tailored coupled clusters method (TCCSD), which employs the domain-based local pair natural orbital approach (DLPNO-TCCSD). Compared to the previous LPNO version of the method, the new implementation is more accurate, offers more favorable scaling and provides more consistent behavior across the variety of systems. On top of the singles and doubles, we include the perturbative triples correction (T), which is able to retrieve even more dynamic correlation. The methods were tested on three systems: tetramethyleneethane, oxo-Mn(Salen) and Iron(II)-porphyrin model. The first two were revisited to assess the performance with respect to LPNO-TCCSD. For oxo-Mn(Salen), we retrieved between 99.8-99.9% of the total canonical correlation energy which is the improvement of 0.2% over the LPNO version in less than 63% of the total LPNO runtime. Similar results were obtained for Iron(II)-porphyrin. When the perturbative triples correction was employed, irrespective of the active space size or system, the obtained energy differences between two spin states were within the chemical accuracy of 1 kcal/mol using the default DLPNO settings.

Tables

Table 1. Table 1 : Energy differences in kcal/mol between DLPNO-TCCSD and DLPNO-CCSD(T) calculations with different settings of cut-off parameters and equivalent canonical calculations on oxo-Mn(Salen) in cc-pVDZ.

		CAS(28,22)		CAS(28,27)
		default	TightPNO	default	TightPNO
SD	S	$2.52$	$3.45$	$3.60$	$4.30$
	T	$1.92$	$3.01$	$3.06$	$3.97$
	$\Delta E_{\text{S–T}}$	$0.60$	$0.44$	$0.54$	$0.33$
SD(T)	S	$6.38$	$5.30$	$7.33$	$6.19$
	T	$5.54$	$4.79$	$6.67$	$5.80$
	$\Delta E_{\text{S–T}}$	$0.84$	$0.51$	$0.66$	$0.39$

Table 2. Table 2 : The singlet and triplet energies of oxo-Mn(salen) in cc-pVXZ basis sets. the difference E ( 1 A ) − E ( 3 A ) E(^{1}\text{A})-E(^{3}\text{A}) in kcal/mol. Results for different active spaces and in various basis sets.

	–DZ	–TZ	–QZ
LPNO-TCCSD(28,22)Antalík et al. (2019)	$6.2$	$6.3$	$6.3$
DLPNO-TCCSD(28,22)	$6.9$	$7.3$	$6.9$
DLPNO-TCCSD(T)(28,22)	$5.9$	$6.3$	$5.8$
LPNO-TCCSD(28,27)Antalík et al. (2019)	$3.7$	$3.1$	$2.9$
DLPNO-TCCSD(28,27)	$4.1$	$4.0$	$3.1$
DLPNO-TCCSD(T)(28,27)	$4.0$	$5.0$	$2.9$
NEVPT2(28,22)Sharma et al. (2017)	$-7.4$	$1.6$	$2.4$

Table 3. Table 3 : Energy differences in kcal/mol between DLPNO-TCCSD and DLPNO-CCSD(T) calculations with different settings of cut-off parameters and equivalent canonical calculations on Fe(II)-porphyrin model in def2-SVP.

		CAS(8,12)		CAS(12,16)		CAS(32,34)
		default	TightPNO	default	TightPNO	default	TightPNO
SD	T	$0.85$	$3.01$	$1.10$	$2.27$	$1.42$	$3.28$
	Q	$0.69$	$2.41$	$0.73$	$2.04$	$2.16$	$3.40$
	$\Delta E_{\text{T–Q}}$	$0.16$	$0.60$	$0.37$	$0.23$	$-0.74$	$-0.12$
SD(T)	T	$4.87$	$5.03$	$4.77$	$4.20$	$5.49$	$5.43$
	Q	$4.58$	$4.38$	$4.34$	$3.90$	$6.19$	$5.52$
	$\Delta E_{\text{T–Q}}$	$0.29$	$0.65$	$0.43$	$0.30$	$-0.70$	$-0.09$

Equations

$$|{\Psi }_{\text{TCC}}⟩=e^T|{\Phi }_0⟩=e^{T_{\text{ext}}+T_{\text{CAS}}}|{\Phi }_0⟩=e^{T_{\text{ext}}}{e}^{T_{\text{CAS}}}|{\Phi }_0⟩,$$

$$|{\Psi }_{\text{MPS}}⟩=\sum _{\{\alpha \}}{\mathbf{A}}^{{\alpha }_1}{\mathbf{A}}^{{\alpha }_2}\cdots {\mathbf{A}}^{{\alpha }_k}|{\alpha }_1{\alpha }_2\cdots {\alpha }_k⟩,$$

$$T_{\text{CAS}}^{(1)}$$

$$T_{\text{CAS}}^{(2)}$$

$$⟨{\Phi }_i^a|H{e}^{T_{\text{ext}}}{e}^{T_{\text{CAS}}}|{\Phi }_0{⟩}_c$$

$$⟨{\Phi }_{ij}^{ab}|H{e}^{T_{\text{ext}}}{e}^{T_{\text{CAS}}}|{\Phi }_0{⟩}_c$$

$$|{\tilde{\mu}}\rangle=\Big{(}1-\sum_{i}|{i}\rangle\langle{i}|\Big{)}|{\mu}\rangle$$

$${\text{(DOI)}}_{i\mu }=\sqrt{{\int }_{{\mathbb{R}}^3}|{\varphi }_i{|}^2|{\varphi }_{\tilde{\mu }}{|}^2{\mathrm{d}}^{3}r}$$

$$\begin{gathered} \Delta E_{\text{CutPairs}}=\sum _{ij}^{\text{weak}}{\epsilon }_{ij}^{\text{PAO}} \\ +\sum _{ip}^{\text{weak}}{\epsilon }_{ip}^{\text{PAO}}+\sum _{pq}^{\text{weak}}{\epsilon }_{pq}^{\text{PAO}}, \end{gathered}$$

$${\mathbf{S}}^{ij}={\mathbf{I}}_{N_{\text{CAS}}}\oplus {\mathbf{d}}^{ij},$$

$${\mathbf{D}}^{pp}={\mathbf{t}}_{pp}^{(1)†}{\mathbf{t}}_{pp}^{(1)},$$

$${\mathbf{t}}_{pp}^{(1)}=\sum _{{\tilde{\mu }}_{pp}{\tilde{\nu }}_{pp}}\frac{⟨p{\tilde{\mu }}_{pp}||p{\tilde{\nu }}_{pp}⟩}{2f_{pp}-{\epsilon }_{{\tilde{\mu }}_{pp}}-{\epsilon }_{{\tilde{\nu }}_{pp}}},$$

$$\begin{array}{rl}\Delta E_{\text{CutPNO}}=& \sum _{ij}^{\text{strong}}({\epsilon }_{ij}^{\text{PAO}}-{\epsilon }_{ij}^{\text{PNO}})\\ +& \sum _{ip}^{\text{strong}}({\epsilon }_{ip}^{\text{PAO}}-{\epsilon }_{ip}^{\text{PNO}})\\ +& \sum _{pq}^{\text{strong}}({\epsilon }_{pq}^{\text{PAO}}-{\epsilon }_{pq}^{\text{PNO}}),\end{array}$$

$$⟨{\Phi }_i^{\bar{a}}|H{e}^{{\bar{T}}_{\text{ext}}^{(1)}+{\bar{T}}_{\text{ext}}^{(2)}}{e}^{T_{\text{CAS}}}|{\Phi }_0{⟩}_c=0\{i,a\}\not\subset \text{CAS},$$

$$⟨{\Phi }_{ij}^{\bar{a}\bar{b}}|H{e}^{{\bar{T}}_{\text{ext}}^{(1)}+{\bar{T}}_{\text{ext}}^{(2)}}{e}^{T_{\text{CAS}}}|{\Phi }_0{⟩}_c=0\{i,j,a,b\}\not\subset \text{CAS}.$$

$$\mathbf{D}^{ijk}=\frac{1}{3}\Big{(}\mathbf{D}^{ij}+\mathbf{D}^{jk}+\mathbf{D}^{ik}\Big{)},$$

Full text

domain-based local pair natural orbital approximation, tailored coupled clusters, density matrix renormalization group

I Introduction

Since its introduction to quantum chemistry Čížek (1966), the coupled cluster (CC) approach has become one of the most widely used methods for the accurate calculations of dynamic correlation. It offers numerous favorable properties, such as compact description of the wave function, size-extensivity, invariance to rotations within occupied or virtual orbital subspaces and also a systematic hierarchy of approximations converging towards the full configuration interaction (FCI) limit Gauss (1998). For instance, the CCSD(T) method Raghavachari et al. (1989), which includes connected single-, double- and perturbative triple excitations, is notorious referred to as the gold standard of quantum chemistry Gauss (1998).

Although the CC method performs well for single reference molecules, it becomes fairly inaccurate or breaks down completely for systems with strongly correlated electrons. Such systems are multireference in nature since they include quasi-degenerate frontier orbitals. This situation is common during dissociation processes, in diradicals, or compounds containing transition metals. Over the years, numerous efforts to generalize the CC ansatz and thus overcome this drawback gave rise to a broad family of multireference CC methods (MRCC) Bartlett and Musiał (2007); Tew et al. (2010); Lyakh et al. (2012).

A possible approach for including static correlation in the CC scheme is to employ a different method like complete active space self-consistent field (CASSCF) or multireference configuration interaction (MRCI) in order to extract the information about the most important excitations Li and Paldus (1997); Li (2001); Li and Paldus (2001, 2001); Paldus and Planelles (1994); Piecuch et al. (1996); Li et al. (1997); Kinoshita et al. (2005); Lyakh et al. (2011); Melnichuk and Bartlett (2012, 2014); Piecuch et al. (1993); Piecuch and Adamowicz (1994); Adamowicz et al. (1998); Piecuch (2010); Piecuch et al. (2002, 2002); Kowalski and Piecuch (2002); łoch et al. (2006); Lodriguito et al. (2006); Piecuch et al. (2004); Kowalski and Piecuch (2001, 2000). The retrieved information can then be introduced to a CC calculation as an external correction. One of such methods is tailored CC with single and double excitations (TCCSD) proposed by Kinoshita et al. Kinoshita et al. (2005), which draws on the split-amplitude ansatz by Piecuch et al.Piecuch et al. (1993), in which the cluster operator corresponding to single and double excitations is split into two parts. The active part is imported from a complete active space configuration interaction (CAS-CI) and kept fixed, while the external amplitudes are iterated using the standard CCSD framework. We recently extended this approach by using the density matrix renormalization group (DMRG) method to obtain the active space amplitudes Veis et al. (2016, 2017). Related externally corrected CC approaches employ fixed $T_{3}$ and $T_{4}$ amplitudes obtained from $MRCI$ Li and Paldus (1997) or from stochastic CI Deustua et al. (2018) and iterate all singles and doubles in their presence. Compared to TCC, this has the advantage that the active space $T_{1}$ and $T_{2}$ can reflect the dynamic correlation outside of the active space, but one pays the price of much larger number of the $T_{3}$ and $T_{4}$ amplitudes involved.

The DMRG method, which originated in solid-state physics White and Noack (1992); White (1992, 1993), is nowadays well established in quantum chemistry for the treatment of strongly correlated systems White and Martin (1999); Chan and Head-Gordon (2002); Legeza et al. (2003, ); Marti and Reiher (2010); Chan and Sharma (2011); Wouters and Neck (2014); Szalay et al. (2015); Yanai et al. (2014). As a numerical approximation to FCI, it can handle significantly larger active spaces compared to the conventional methods. However, even then the prohibitive scaling does not allow to include dynamic correlation and it is therefore necessary to employ some "post-DMRG" procedure. Many different attempts has been made to tackle this limitation for example with the complete active space perturbation theory (CASPT2)Kurashige and Yanai (2011), Cholesky decomposition n-electron valence state perturbation theory (NEVPT2)Freitag et al. (2017), MRCI using cumulant reconstruction with internal contraction of DMRG wave functionSaitow et al. (2013), canonical transformationNeuscamman et al. (2010), matrix product state (MPS) based formulation of multireference perturbation theorySharma and Chan (2014), DMRG pair-density functional theorySharma et al. (2019), and also our aforementioned CC tailored by MPS wave functions (DMRG-TCCSD)Veis et al. (2016); Faulstich et al. (2019, 2019).

Even though the DMRG-TCCSD method offers a reasonably efficient treatment of both static and dynamic correlationVeis et al. (2018), its application to larger systems is hampered by the unfavorable scaling of the CCSD part of the calculation. With such a steep scaling, even massive parallelization is not sufficient to make the method applicable to molecules with hundreds of atoms. To overcome this issue, Pulay proposed to exploit the locality of the electron correlationPulay (1983); Sæbø and Pulay (1985). Due to its short range character for non-metallic systems, it is possible to take advantage of sparsity of the Hamiltonian matrix by employing the basis of localized orbitals.

Although the occupied orbital space can be easily localized using a variety of appropriate methods Foster and Boys (1960); Pipek and Mezey (1989); Knizia (2013), for the virtual space things get slightly more complicated. In their first works on locality, Pulay and Sæbø used projected atomic orbitals (PAOs) Sæbø and Pulay (1987, 1988), which were also used by Werner and Schütz in the local CC method Hampel and Werner (1996); Schütz and Werner (2001); Schütz (2002); Werner and Schütz (2011); Schütz (2002). In this approach, each localized occupied orbital is assigned a domain of PAOs obtained by projecting out the occupied orbital components from atomic orbitals. The pairs of occupied orbitals are subsequently classified according to their real space distance and treated at either coupled cluster level (strong pairs), perturbative level (weak and distant pairs), or neglected altogether (very distant pairs).

Another idea how to make use of locality is based on the concept of dividing a large system into smaller segments and performing the calculations on each of these subsystems separately. Among such approaches belong the divide-expand-consolidate method Kristensen et al. (2011); Høyvik et al. (2012), the divide-and-conquer method Kobayashi and Nakai (2008), the incremental method Stoll (1992), the local natural orbital method Rolik and Kállay (2011), and the fragment molecular orbital method Fedorov and Kitaura (2005). The closely related cluster-in-molecule method Li et al. (2001) is based on energy decomposition into contributions corresponding to individual occupied orbitals.

Possibly the most effective way of virtual space truncation is, however, the use of pair natural orbitals (PNOs), which are known to provide compact parametrization of the virtual space. They were first used in the 1960s by Edmiston and Kraus Edmiston and Krauss (1965) and later by Meyer Meyer (1971, 1974); Werner and Meyer (1976); Botschwina and Meyer (1977); Rosmus and Meyer (1978), Ahlrichs and Kutzelnigg Ahlrichs and Kutzelnigg (1968); Ahlrichs et al. (1975). For many years, the progress in the area of PNO-based methods stalled, until their revival in 2009, when the new local pair natural orbital (LPNO) variants of CEPA and CCSD methods were introduced by Neese et alNeese et al. (2009, 2009); Hansen et al. (2011); Huntington et al. (2012). The cornerstone of this approach, the idea to use PNOs in combination with localized occupied orbitals, was later developed into the more advanced domain based local pair natural orbital (DLPNO) methodsRiplinger and Neese (2013); Pinski et al. (2015); Riplinger et al. (2016); Saitow et al. (2017).

In these, PNOs were expressed as a linear combination of PAOs in a pair-domain, which ultimately removed the bottleneck of previous PNO methods and achieved genuine linear scaling. For instance, the resulting DLPNO-CCSD method is applicable to systems with hundreds of atoms and thousands of basis functions, which renders the prior SCF calculation possibly computationally more demanding than the actual correlation treatment.

Moreover, the PNO-based approaches possess many desirable properties, which allow them to be used in black box fashion. They provide a very compact description of the virtual space, which makes it computationally feasible to use sufficiently large domains of PAOs; something that would be too costly for a purely PAO-based approach. They use a limited number of cut-off parameters, which do not involve any real space distance and the calculated correlation energy smoothly depends on the values of these parameters.

Currently, PNO-based methods are developed in a number of groups including WernerWerner et al. (2015); Menezes et al. (2016); Schwilk et al. (2017); Tew et al. (2011) and Hättig Helmich and Hättig (2013); Schmitz et al. (2013); Guo et al. (2016) and are widely employed in the context of various systems of chemical interestAntony et al. (2011); Anoop et al. (2010); Liakos and Neese (2011); Zade et al. (2011); Kubas et al. (2012); Ashtari and Cann (2011); Zhang and Dolg (2014); Minenkov et al. (2015); Sparta and Neese (2014); Liakos et al. (2015). Apart from single-reference methods, the methodology was also successfully applied to multireference CC techniques Demel et al. (2015); Lang et al. (2017); Brabec et al. (2018); Lang et al. (2019).

In this article, we build upon the previous investigation of our LPNO-TCCSD methodAntalík et al. (2019) and introduce the newly implemented DLPNO-TCCSD and DLPNO-TCCSD(T) methods. Based on our experience, we revisit the molecule of tetramethyleneethane (TME) and address the drawbacks of the former method. The performance of the new approximation is then assessed using two benchmark systems, namely, oxo-Mn(Salen) and a model of Fe(II)-porphyrins (FeP).

II Theory and Implementation

II.1 DMRG-based Tailored Coupled Clusters

The tailored coupled cluster methodKinoshita et al. (2005), which belongs to the class of externally corrected methods, employs the split-amplitude wave function ansatz proposed by Piecuch et al. Piecuch et al. (1993)

[TABLE]

where $|\Phi_{0}\rangle$ is the reference wave function and the cluster operator $T$ is split into two parts: $T_{\text{CAS}}$ which contains the active amplitudes obtained from an external calculation and $T_{\text{ext}}$ which contains the external amplitudes i.e. the amplitudes with at least one index outside the CAS space. Another way to justify this ansatz is the formulation of CC equations based on excitation subalgebras recently introduced by Kowalski Kowalski (2018); Bauman et al. (2019).

In our implementation, we employed the DMRG method to obtain the active amplitudes. Using the DMRG algorithm we first optimize the wave function, which is provided in the MPS form

[TABLE]

where $\alpha\in\{|-\rangle,|\downarrow\rangle,|\uparrow\rangle,|\downarrow\uparrow\rangle\}$ and $\mathbf{A}^{\alpha_{i}}$ are MPS matrices. These are then contracted to obtain CI coefficients for single and double excitations $C$ Moritz and Reiher (2007); Boguslawski et al. (2011). Using the intermediate normalization and the relations between CI and CC coefficients

[TABLE]

we are able to acquire the respective CC amplitudes, which are subsequently introduced into the CC calculation. At this point, these active amplitudes are kept frozen, while the remaining amplitudes in $T_{\text{ext}}$ are optimized by solving the equations

[TABLE]

analogously to the standard CCSD equations. This way, the active amplitudes account for static correlation and by optimizing the external amplitudes, we are able to recover the remaining dynamic correlation. The effect of dynamic correlation on $T_{\text{CAS}}$ is neglected, which is an approximation inherent in the TCC method.

Due to the two-body Hamiltonian, TCC recovers the DMRG energy for $T_{\text{ext}}=0$. In the limit of CAS including all orbitals, FCI energy is recovered, although the TCC energy does not behave monotonously when extending the CAS spaceFaulstich et al. (2019). Nevertheless, in practice the optimal CAS size related to the energy minimum can be determined with low cost DMRG calculationsFaulstich et al. (2019). In addition, a quadratic error bound valid for DMRG-TCC methods is also derivedFaulstich et al. (2019).

On top of the TCCSD routine, the perturbative triples correction can be appliedLyakh et al. (2011). However, in order to prevent the double counting of static correlation, it is necessary to omit the terms that include the active amplitudes. This can be straightforwardly achieved by setting all active single and double amplitudes to zero during the calculation of the (T) correction.

II.2 Domain Local Pair Natural Orbital Approximation for TCCSD

The presented method is based on the open-shell DLPNO-CCSD code as implemented in ORCANeese (2011). In this section, we briefly outline DLPNO-CCSD and describe the modifications that were made to accommodate the tailored version of the method.

As in all PNO-based local methods, the whole process starts with the localization of the occupied orbitals. Based on our previous experienceAntalík et al. (2019), we opt for the split-localization scheme, where the orbitals are separately localized within four distinct orbital subspaces: doubly occupied external, doubly occupied active, singly occupied active and active virtual. Such choice has been shown to provide a set of orbitals which yields a reasonable convergence behavior of the DMRG procedureOlivares-Amaya et al. (2015), without influencing the DMRG energy.

The next step is the construction of orbital domains. Using the idea of Werner et al.Adler and Werner (2011, 2009), we first construct PAOs

[TABLE]

by projecting out the localized occupied and active orbitals $|i\rangle$ from the original set of atomic orbitals $|\mu\rangle$. The acquired orbitals are subsequently normalized and the orbital domains are constructed based on the differential overlap integrals

[TABLE]

between PAOs and the set of occupied and active orbitals. Here, the first prescreening parameter comes into play – for a given occupied orbital $i$, only PAOs for which $\text{(DOI)}_{i\mu}>T_{\text{CutDO}}$ are included within its domain. Also, if an atom contains at least one PAO, all PAOs belonging to this atom are considered for further domain construction.

At this point, the main difference of the TCC approach compared to the conventional DLPNO-CCSD is that all active indices including the active virtuals are formally treated as singly occupied. This means that during the creation of domains they share the same domain. Moreover, in the following dipole prescreening, it is ensured that every active occupied pair (i.e. a pair with both active indices) automatically survives the dipole prescreening.

After the prescreenings, unrestricted MP2 pair energies $\varepsilon_{ij}^{\text{PAO}}$ are calculated and used to categorize occupied orbital pairs $ij$ into three classes, based on the preset cut-off parameters. Specifically, the active pairs and the pairs with energies larger than $T_{\text{CutPairs}}$ are classified as strong, while the rest is either weak or neglected according to a related parameter. The whole process is executed in two consecutive steps. First, a crude screening is performed in smaller domains using loose thresholds, followed by the second screening in which the strong and weak pairs are again distributed between the three categories, but this time with finer thresholds. The strong pairs are then passed to the next stage, while the remaining pair energies are stored as a correction to the final energy

[TABLE]

where $i$, $j$ are indices of doubly and $p$, $q$ of singly occupied orbitals. If the perturbative triples correction is invoked, the final weak pairs with energies higher than $0.01\cdot T_{\text{CutPairs}}$ are also saved for later use.

Afterwards, using the non-redundant PAOs, the NEVPT2 pair densities $\mathbf{D}^{ij}$ are constructed for the surviving pairs which do not contain any explicit information about the tailored CAS space. These are diagonalized to obtain PNO expansions $\mathbf{d}^{ij}$, which are then truncated based on the cut-off parameter $T_{\text{CutPNO}}$. Only PNOs with occupation numbers larger than its value are kept and the remaining orbitals are discarded. The final PAO/PNO transformation matrix for a given pair is then obtained by enlarging the former transformation matrix by a unit matrix

[TABLE]

where $N_{\text{CAS}}$ is the number of active SOMOs and virtual orbitals, see Figure 1. While this artificially increases the final number of PNOs, it barely increases the computational time. Singles PNOs are obtained in the same manner from the densities which are computed as

[TABLE]

where the amplitudes are

[TABLE]

where $f_{pp}$ are diagonal Fock matrix elements. Furthermore, the pair energies are recalculated in the new truncated PNO basis and the differences between them and the original MP2 estimates in the PAO basis

[TABLE]

are stored as a correction to the final energy.

The resulting equations for singly excited amplitudes (5) now become

[TABLE]

where the barred index $\bar{a}$ indicates PNO basis. Similarly, the equations for doubly excited amplitudes (6) become

[TABLE]

Once these equations are solved, the stored corrections (9) and (13) are added to the acquired energy. Except the fact that the active amplitudes are ‘frozen’ during the CCSD iterations, these equations are identical to single-reference DLPNO-CCSD as implemented in ORCASaitow et al. (2017).

II.3 Perturbative triples correction to DLPNO-TCCSD

To calculate the triples correction using the DLPNO approachRiplinger et al. (2013); Guo et al. (2018), it is first necessary to identify relevant compact parameterization of virtual space for every triple $ijk$. Domains for these triples are created as a union of $i$, $j$ and $k$ domains and the $ij$, $jk$ and $ik$ pair densities are created within these triples domains.

Next, the triples densities are constructed by averaging over the respective pair densities

[TABLE]

where the triple $ijk$ is composed either of three strong pairs or of two strong and a weak pair, as it has been previously shown that using only the strong pairs is insufficientRiplinger et al. (2013). Since the amplitudes for weak pairs are not known prior to the triples calculation, the approximate MP2 amplitudes are used. The process then continues analogously to the construction of PNOs in CCSD.

Once the densities are constructed, they are diagonalized in order to obtain the triple natural orbitals (TNO) and their corresponding natural occupation numbers. The TNO expansion is then truncated based on the occupation numbers and the cut-off parameter $T_{\text{CutTNO}}$ and from the orbitals that passed the screening a transformation matrix is formed. This matrix is enlarged by a unit matrix of dimension $N_{\text{CAS}}$ and its final form is used to transform the integrals, amplitudes and subsequently to calculate the energy correction. As in the canonical version of the method, all active single and double amplitudes are set to zero during the calculation of the correction to prevent double-counting.

III Computational Details

The DMRG calculations were performed by the Budapest QC-DMRG code Legeza et al. . The DLPNO-TCCSD and DLPNO-TCCSD(T) methods were implemented in the ORCA program package Neese (2011), which was also used to prepare the orbitals.

Prior to DMRG calculations, we split-localized the orbitals by the Pipek-Mezey methodPipek and Mezey (1989) in the following orbital subspaces: internal, active doubly occupied, active singly occupied and active virtual. Orbital ordering was subsequently optimized using the Fiedler method Barcza et al. (2011); Fertitta et al. (2014) combined with minor manual adjustments. All DMRG runs were initialized by the CI-DEAS procedure Legeza and Sólyom (2003); Szalay et al. (2015) and employed the dynamical block state selection (DBSS) procedure Legeza et al. (2003); Legeza and Sólyom (2004) to control the accuracy of the calculation, with the truncation error criterion set to $10^{-6}$. The convergence threshold was set to energy difference between two subsequent sweeps smaller than $10^{-6}$ a.u.

The core electrons were kept frozen throughout all coupled cluster calculations. The default set of DLPNO cut-off parameters $T_{\text{CutPNO}}=3.33\cdot 10^{-7}$, $T_{\text{CutPairs}}=10^{-4}$, $T_{\text{CutMKN}}=10^{-3}$ and $T_{\text{CutDO}}=10^{-2}$ was employed unless otherwise stated. The other settings, referred to as TightPNO, were used for production runs: $T_{\text{CutPNO}}=10^{-7}$, $T_{\text{CutPairs}}=10^{-5}$, $T_{\text{CutMKN}}=10^{-3}$ and $T_{\text{CutDO}}=5.0\cdot 10^{-3}$. For perturbative triples the default value of relevant cut-off parameter is $T_{\text{CutTNO}}=10^{-9}$ for both sets of settings. To estimate the dependence of DLPNO-TCCSD energies on these parameters, one parameter was varied with remaining parameters fixed to the default value. We assess the amount of retrieved correlation energy by DLPNO approach with reference to the DMRG-TCCSD energy calculated with the canonical TCCSD implementation.

For TME, we used CASPT2(6,6)/cc-pVTZ geometries for seven values of the dihedral angle from our previous work Veis et al. (2018). The orbitals were prepared by CASSCF(6,6) calculation with the active space containing six 2pz orbitals on carbon atoms. The cc-pV6Z/C auxiliary basis set was used for the resolution of the identity (RI) approximation Bross et al. (2013).

For oxo-Mn(Salen), we used the singlet CASSCF(10,10)/6-31G* optimized geometry by Ivanic et al. Ivanic et al. (2004). The orbitals were optimized using the DMRG-CASSCF method Ghosh et al. (2008); Zgid and Nooijen (2008); Yanai et al. (2009) in Dunning’s cc-pVXZ $\text{X}\in${D,T,Q} basis sets Dunning (1989); Woon and Dunning (1993); Balabanov and Peterson (2005). The optimization was carried out with fixed bond dimension $M=1024$ for the smaller CAS(28,22) and $M=2048$ for CAS(28,27). The composition of these active spaces is discussed further in detail in the paperAntalík et al. (2019). The cc-pVQZ/C auxiliary basis set was used for RI Weigend et al. (2002).

Finally, the calculations on FeP were performed at the triplet geometry from the previous study of Li Manni and Alavi Manni and Alavi (2018), with CASSCF/def2-SVP orbitals, for three distinct active spaces, which were selected based on the entanglement analysis . For the largest CAS the DMRG-CASSCF was applied, with the fixed bond dimension $M=1024$. The def2-SVP and def2-TZVP basis setsWeigend and Ahlrichs (2005) were used with the def2-TZVPP/C auxiliary basis setHellweg et al. (2007) for RI.

IV Results and Discussion

IV.1 Tetramethyleneethane

Although small, the tetramethyleneethane molecule is a challenging system due to its complex electronic structure. To correctly describe the character of its singlet state, one needs to employ a theory with a balanced description of both static and dynamic correlation combined with a reasonably large basis set. This is the reason, why it often serves as a benchmark system for multireference methods Pittner et al. (2001); Bhaskaran-Nair et al. (2011); Chattopadhyay et al. (2011); Pozun et al. (2013); Demel et al. (2015); Veis et al. (2018). We previously studied the system with the canonical DMRG-TCCSD method Veis et al. (2018), as well as with its LPNO versionAntalík et al. (2019). In the latter, we encountered an issue with different accuracy for singlet and triplet states, which we attributed to the neglect of some terms in the LPNO-CCSD implementation in ORCA. For this reason, we revisit the system to investigate the improvement by the new DLPNO method, which should offer more robust approximation.

We performed calculations in seven different geometries corresponding to the rotation about the central C–C bond of the molecule (see Figure 2a) and different values of the cut-off parameters. These benchmarks were performed only for $T_{\text{CutPairs}}$ and $T_{\text{CutPNO}}$ parameters, since $T_{\text{CutDO}}$ does not affect the results, unless an extremely small value is chosen. This behavior corresponds with the observations from previous studiesLang et al. (2019); Brabec et al. (2018); Saitow et al. (2017).

The left plot in Figure 3 shows the dependence of retrieved canonical correlation energy with respect to estimated pair correlation energies cut-off $T_{\text{CutPairs}}$ averaged over the geometries. Both methods converge in a similar fashion, but DLPNO-TCCSD recovers 0.1–0.3% more correlation energy than LPNO-TCCSD. The right plot in the same figure shows the dependence on the second cut-off parameter, that is PNO occupation number $T_{\text{CutPNO}}$, which is again averaged over the geometries. Here, the DLPNO-TCCSD shows noticeably faster convergence to the canonical value compared to LPNO-TCCSD, with more accurate correlation energies even for higher values of the parameter. For instance, when the default values are used, DLPNO-TCCSD extracts over 99.91% of the canonical correlation energy.

One can notice the aforementioned discrepancy in accuracy for LPNO-TCCSD method, which can be explained by the neglected terms in the UHF-LPNO formalism. DLPNO-TCCSD does not exhibit this behavior and describes both spin states equally well. This means that the difference of the recovered correlation energy is less than 0.01% compared to LPNO-TCCSD, for which this percentage is an order of magnitude worse.

Figure 4 shows the dependence of non-parallelity error (NPE) on $T_{\text{CutPNO}}$ for both LPNO and DLPNO version of TCCSD. It can be seen that the NPE is no larger than 0.17 kcal/mol for either spin state, with very similar value for TightPNO settings, as well as for the LPNO version.

From the chemical point of view, the most interesting property is the behavior of the singlet-triplet gap with respect to the cut-off parameters. It can be observed from the presented data that this property is well described at the default and even looser settings, which means that it can be calculated by DLPNO-TCCSD with virtually no loss in accuracy.

IV.2 oxo-Mn(Salen)

The oxo-Mn(Salen) molecule (Figure 2b) has been a subject of numerous computational studies motivated mainly by its role in catalysis of the enantioselective epoxidation of unfunctional olefines Zhang et al. (1990); Irie et al. (1990). Due to its closely lying singlet and triplet states it is also considered to be a challenging system even for multireference methods. Over the years, many multireference studies have been published Ivanic et al. (2004); Sears and Sherrill (2006); Ma et al. (2011), some of which employed the DMRG method Wouters et al. (2014); Olivares-Amaya et al. (2015); Stein and Reiher (2016) and DMRG results with dynamic correlation treatment were presented as wellVeis et al. (2016); Sharma et al. (2017); Antalík et al. (2019). This includes our previous studies, in which we used the system as a benchmark for the predecessors of the current method, namely DMRG-TCCSD and its LPNO version.

Again, we first assessed the percentage of the recovered canonical TCCSD correlation energy with respect to the cut-off parameter $T_{\text{CutPairs}}$. The curve shown in the left plot of Figure 5 follows the same trend as for TME and quickly converges. This means that the retrieved correlation energy is already converged by the value $T_{\text{CutPairs}}=10^{-5}$. Even though the correlation energy is not yet stable for the default value $10^{-4}$, the difference in accuracy between two spin states is minimal. Since we are interested in the width of the singlet-triplet gap, it appears that even this setting provides quite reasonable accuracy.

The right plot of Figure 5 shows the dependence on the cut-off parameter $T_{\text{CutPNO}}$. These curves smoothly near towards 100% as the parameter tightens, the behavior as observed for conventional DLPNO-CCSD. However, the energies overshoot for the very conservative values of the parameter, which is caused by overcompensation of the neglected pairs with the MP2 pair energy. For this reason, we calculated the same dependence also with the tighter setting of $T_{\text{CutPairs}}=10^{-5}$ (the thinner dotted lines in the plot). Since for this value the correlation energies are basically converged with respect to $T_{\text{CutPairs}}$, we can observe that the curves now smoothly converge towards 100%. The singlet-triplet gap seems reasonably accurate for the default value of $T_{\text{CutPNO}}$ and for more conservative values the difference in accuracy completely disappears.

We also examined the effect of these cut-offs on the perturbative triples correction recovered by DLPNO-TCCSD(T), which are shown in Figure 6. At the default value of $T_{\text{CutPairs}}$ the energies seem to be almost converged and are fully converged at the TightPNO setting, the actually relevant relative accuracy does not change. In case of $T_{\text{CutPNO}}$, the curves seem to slowly converge towards 100%, which is even more apparent when tighter cut-off on pair energies is in place. Note, that these values represent only the percentage of the canonical triples correction not the total correlation energy and for this particular system, the triples correction amounts to less than 3% of the total correlation energy. On top of that, we investigated the sensitivity of overall correlation energy on the $T_{\text{CutTNO}}$ with the conclusion that the default value is more than adequate, see Figure 7.

Table 1 contains the differences between DLPNO and canonical energies for the singlet-triplet gap. Although the errors in energies for singlet or triplet state alone are substantial (especially for the larger CAS), the errors for the gap are well within the chemical accuracy of 1 kcal/mol. When the perturbative triples correction is invoked, tighter cut-offs are preferable. The same can be observed in Figure 8, where the chemical accuracy is achieved even with looser than default settings. Also the behavior of the DLPNO approximation seems quite stable with respect to the size of the active space.

When we compare these results to the previous LPNO-CCSD calculations Antalík et al. (2019), the DLPNO-TCCSD has smaller error in absolute energies, but larger in gaps. This can be attributed to the rather fortunate cancellation of errors in the LPNO case, but otherwise the DLPNO implementation is more reliable due to smaller errors in absolute energies and smooth convergence towards the canonical correlation energies.

Moreover, we performed calculations with up to cc-pVQZ basis set. We found the results to be consistent with the previous LPNO-TCCSD calculations and NEVPT2(28,22), with perturbative triples correction being responsible only for a minor shift in energy, see Table 2. Even though the calculation in cc-pVQZ, which amounts to 1178 basis functions, is perfectly within the possibilities of the LPNO approach, it took 60% longer to finish under the same conditions (8 cores, 30GB of memory per core) compared to the DLPNO calculation.

IV.3 Iron(II) Porphyrin Model

Fe(II) porphyrin derivatives play important roles in reactions related to material science and biological processes due to their closely lying spin states. The chosen model has previously been a subject of several large scale CASSCF studiesManni and Alavi (2018); Manni et al. (2019) and we therefore consider it to be an interesting system to test the efficiency of the method for several active spaces of different size.

The left plot in Figure 9 shows the dependence on the $T_{\text{CutPairs}}$ cut-off parameter. As for oxo-Mn(Salen), the curves converge smoothly and at $10^{-5}$ they are practically converged. The accuracy is consistent for triplet and quintet spin states and it differs maximally by 0.05% of the retrieved canonical correlation energy. This is true regardless of the size of the active space, although the overall accuracy is slightly lower for the largest CAS.

The right plot in Figure 9 shows the dependence on the $T_{\text{CutPNO}}$ cut-off parameter. For very loose values, the discrepancy in accuracy for different spin states is apparent, especially for the larger CAS. However, it disappears with default and TightPNO settings of the cut-off. The method once again overestimates the correlation energy for the default value of $T_{\text{CutPairs}}$, but for tighter value $T_{\text{CutPairs}}=10^{-5}$ it converges smoothly towards 100%.

Regarding the triples correction, the convergence behavior is similar as for oxo-Mn(Salen). The percentages of retrieved triples energy with respect to both $T_{\text{CutPairs}}$ and $T_{\text{CutPNO}}$ are shown in Figure 10. The DLPNO approximation for this correction is again the most sensitive to changes in $T_{\text{CutPNO}}$ parameter, but shows clear convergence towards 100%. Dependence of total DLPNO correlation energy on $T_{\text{CutTNO}}$ is shown in Figure 11.

Finally, we assess the accuracy of the triplet-quintet gap, which is the actual value of interest. Table 3 lists the energy differences between TCCSD, TCCSD(T) and their DLPNO versions. At the CCSD level, the errors for absolute energies grow with growing CAS. These errors are larger for the TightPNO settings, which is understandable given the method overshoots for looser $T_{\text{CutPairs}}$ thresholds. However, the tighter cut-offs significantly improve the accuracy in energy of the triplet-quintet gap, with the exception of the smallest active space. Although the decrease in accuracy of the absolute energies is observed, once the perturbative triples correction is included, the errors in gaps are basically the same as in the CCSD case. Moreover, the dependence of the DLPNO error on $T_{\text{CutPNO}}$ is presented in Figure 12. Even for the default cut-offs, the error is well under 1 kcal/mol and is even smaller for the TightPNO settings, with the previously discussed exception of CAS(8,12). Interestingly, the accuracy of the quintet state is consistently lower than for the triplet state for the largest CAS, while it is the opposite for the smaller active spaces.

V Conclusions

We introduced a new version of DMRG-TCCSD method, which employs the domain-based local pair natural orbital approach. The method was implemented in ORCA and should be available in the next release, while the Budapest DMRG code is used for the DMRG part.

We performed several accuracy assessments of the method employing three systems, two of which have been previously studied by the canonical and LPNO versions of the TCCSD method.

For tetramethyleneethane, we were able to retrieve over 99.9% of the canonical correlation energy, while using the default settings of cut-off parameters, with negligible non-parallelity error with respect to its dihedral rotation. Also, the drawbacks present in the LPNO version were eliminated.

For oxo-Mn(Salen), the amount of retrieved correlation was dependent on the size of the used active space, ranging from 99.8% for the larger CAS(28,27) to 99.9% for smaller CAS(28,22), which is an improvement about 0.2% with respect to LPNO-TCCSD. Using the default settings resulted in singlet-triplet gap being off by 0.5-0.6 kcal/mol and with tighter cut-offs only 0.3-0.4 kcal/mol compared to canonical calculation. These results are slightly worse than those of LPNO-TCCSD, which we believe is a consequence of a fortunate cancellation of errors in LPNO-TCCSD, since DLPNO-TCCSD are consistently better for the remaining two systems. Moreover, new implementation offers significantly better timing than the LPNO one.

For Fe(II)-porphyrin, the dependence on cut-off parameters was very similar as for oxo-Mn(Salen). With the default settings, we were able to retrieve more than 99.9% canonical correlation energy. Irrespective of the active space size, the method determined the triplet-quintet gap with error safely within the chemical accuracy of 1kcal/mol even with the default settings.

To summarize, we believe that DLPNO-DMRG-TCCSD(T) is a possible approach to treat systems with a moderate multireference character. In the near future, we would like to implement the DLPNO version of multireference TCCSD, which we hope to further enhance capabilities of this method.

Acknowledgment

We would like to thank Prof. Frank Neese for providing us with access to ORCA source code, as well as for helpful discussions, Dr. Frank Wennmohs for technical assistance with the ORCA code and Dr. Ondřej Demel for helpful discussions.

This work has been supported by the Czech Science Foundation (Grants No. 18-24563S, 16-12052S and 18-18940Y), the Hungarian National Research, Development and Innovation Office (NKFIH) through (Grant No. K120569), the Hungarian Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001), by the Pacific Northwest National Laboratory (Contract No. 462268), and by the Czech Ministry of Education, Youth and Sports (Project No. LTAUSA17033 and the Large Infrastructures for Research, Experimental Development and Innovations project „IT4Innovations National Supercomputing Center – LM2015070“. Ö. Legeza acknowledges financial support from the Alexander von Humboldt foundation. Mutual visits with members of the Hungarian group have been partly supported by the Hungarian-Czech Joint Research Project MTA/19/04.