Towards the Efficient Local Tailored Coupled Cluster Approximation and the Peculiar Case of Oxo-Mn(Salen)

TL;DR

This paper presents a local pair natural orbital coupled cluster method using matrix product states, enabling accurate and efficient calculations on larger systems like oxo-Mn(Salen) with high correlation energy recovery.

Contribution

The authors develop a local LPNO-TCCSD implementation with matrix product states, improving scalability and accuracy for complex transition metal systems.

Findings

Recovered over 99.7% of correlation energy for TME

Achieved over 99.8% correlation energy for TME's singlet state

Successfully performed large basis set calculations on oxo-Mn(Salen)

Abstract

We introduce a new implementation of the coupled cluster method tailored by matrix product states wave functions (DMRG-TCCSD), which employs the local pair natural orbital approach (LPNO). By exploiting locality in the coupled cluster stage of the calculation, we were able to remove some of the limitations that hindered the application of the canonical version of the method to larger systems and/or with larger basis sets. We assessed the accuracy of the approximation using two systems: tetramethyleneethane (TME) and oxo-Mn(Salen). Using the default cut-off parameters, we were able to recover over 99.7% and 99.8% of canonical correlation energy for the triplet and singlet state of TME respectively. In case of oxo-Mn(Salen), we found out that the amount of retrieved canonical correlation energy depends on the size of the active space (CAS) - we retrieved over 99.6% for the larger 27 orbital CAS and over 99.8% for the smaller 22 orbital CAS. The use of LPNO-TCCSD allowed us to perform these calculations up to quadruple-$\zeta$ basis set amounting to 1178 basis functions. Moreover, we examined dependance of the ground state of oxo-Mn(Salen) on CAS composition. We found out that the inclusion of 4d$_{xy}$ orbital plays an important role in stabilizing the singlet state at the DMRG-CASSCF level via double-shell effect. However, by including dynamic correlation the ground state was found to be triplet regardless of the size of the basis set or composition of CAS, which is in agreement with previous findings by canonical DMRG-TCCSD in smaller basis.

Figures

Tables

Table 1. Table 1: Energy differences in kcal/mol between LPNO-TCCSD 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
${}_{}{}^1\text{A}$	$2.53$	$4.50$	$4.91$	$6.47$
${}_{}{}^3\text{A}$	$2.47$	$4.55$	$5.51$	$6.50$
$\Delta E_{\text{T-S}}$	$-0.07$	$0.04$	$0.60$	$0.03$

Table 2. Table 2: The singlet and triplet energies of oxo-Mn(Salen) E + 2251 𝐸 2251 E+2251 in atomic units and the difference Δ E T-S = E ( 3 A ) − E ( 1 A ) \Delta E_{\text{T-S}}=E(^{3}\text{A})-E(^{1}\text{A}) in kcal/mol. Results for different active spaces and in various basis sets.

	cc-pVDZ			cc-pVTZ			cc-pVQZ
	${}_{}{}^1\text{A}$	${}_{}{}^3\text{A}$	$\Delta E_{\text{T-S}}$	${}_{}{}^1\text{A}$	${}_{}{}^3\text{A}$	$\Delta E_{\text{T-S}}$	${}_{}{}^1\text{A}$	${}_{}{}^3\text{A}$	$\Delta E_{\text{T-S}}$
DMRG-CASSCF(28,22)Wouters et al. (2014)	$-0.7509$	$-0.7593$	$-5.3$
DMRG-CASSCF(26,21)Stein and Reiher (2016)	$-0.7963$	$-0.7954$	$0.6$
DMRG-CASSCF(28,22)Sharma et al. (2017)	$-0.7991$	$-0.8002$	$-0.7$	$-0.9957$	$-0.9926$	$1.9$	$-1.0449$	$-1.0418$	$1.9$
DMRG-CASSCF(28,22)	$-0.7502$	$-0.7582$	$-5.0$	$-0.9402$	$-0.9473$	$-4.5$	$-0.9891$	$-0.9961$	$-4.4$
DMRG-CASSCF(28,27)	$-0.8803$	$-0.8791$	$0.8$	$-1.0717$	$-1.0697$	$1.2$	$-1.1210$	$-1.1187$	$1.5$
TCCSD(34,25)Veis et al. (2016) 111These values were obtained with split-localized ROHF orbitals.	$-2.7273$	$-2.7330$	$-3.6$
MRLCC(28,22)Sharma et al. (2017)	$-3.2830$	$-3.2600$	$14.4$	$-4.1303$	$-4.1310$	$-0.5$
NEVPT2(28,22)Sharma et al. (2017)	$-3.0109$	$-2.9990$	$7.4$	$-3.8437$	$-3.8463$	$-1.6$	$-4.1441$	$-4.1481$	$-2.4$
LPNO-TCCSD(28,22)	$-3.1455$	$-3.1554$	$-6.2$	$-3.9698$	$-3.9798$	$-6.3$	$-4.2479$	$-4.2578$	$-6.3$
LPNO-TCCSD(28,27)	$-3.1491$	$-3.1550$	$-3.7$	$-3.9749$	$-3.9798$	$-3.1$	$-4.2531$	$-4.2578$	$-2.9$

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$$

$${}^{\alpha \alpha }{t}_{ij}^{ab}=-\frac{K_{ij}^{ab}-K_{ij}^{ba}}{f_{aa}+f_{bb}-f_{ii}-f_{jj}},$$

$${}^{\alpha \beta }{t}_{ij}^{ab}=-\frac{K_{ij}^{ab}}{f_{aa}+f_{bb}-f_{ii}-f_{jj}},$$

$${\mathbf{D}}_{\alpha \alpha }^{ij}$$

$${}^{(\alpha )}{\mathbf{D}}_{\alpha \beta }^{ij}$$

$${}^{(\beta )}{\mathbf{D}}_{\alpha \beta }^{ij}$$

$${\mathbf{D}}^{ij}{\mathbf{d}}_{\bar{a}}^{ij}=n_{\bar{a}}^{ij}{\mathbf{d}}_{\bar{a}}^{ij},$$

$${\tilde{\mathbf{D}}}^{ij}={\mathbf{D}}_{\text{CAS}}^{ij}\oplus {\mathbf{D}}_{\text{ext}}^{ij},$$

$${\mathbf{D}}_{\text{CAS}}^{ij}$$

$$(D_{\text{ext}}^{ij}{)}_{ab}$$

$$⟨{\Phi }_{ij}^{\bar{a}\bar{b}}|H{e}^{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},$$

Full text

local pair natural orbital approximation, tailored coupled clusters, density matrix renormalization group, oxo-Mn(Salen)

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 notoriously referred to as the gold standard of quantum chemistry Gauss (1998).

Although the CC method performs well for single reference molecules, it becomes fairly inaccuarate or breaks down completely for systems with strongly correlated electrons. Such systems are multireference in nature since they include quasi degenerate frontier orbitals, which are 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 Musial (2007); Tew et al. (2010); Lyakh et al. (2012).

One such approach, aiming to include 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, 2001, 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); Adamovicz et al. (1998); Piecuch (2010); Piecuch et al. (2002, 2002); Kowalski and Piecuch (2002); WŁoch et al. (2006); Lodriguito et al. (2006); Piecuch et al. (2004); Kowalski and Piecuch (2001, 2000). The retrieved information can be then 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, in which the amplitudes corresponding to single and double excitations are split into two parts. The active part is treated by complete active space configuration interaction (CAS-CI) and 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).

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, 2008); Marti and Reiher (2010); Chan and Sharma (2011); Wouters and Van Neck (2014); Szalay et al. (2015); Yanai et al. (2014). As a numerical approximation to full configuration interaction (FCI), it can handle significantly larger active spaces compared to the conventional method. 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 attemps has been made to tackle this limitation for example DMRG-CASPT2Kurashige and Yanai (2011), Cholesky decomposition DMRG-NEVPT2Freitag et al. (2017), DMRG-icMRCISaitow 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. , and also our aforementioned CC tailored by MPS wave functions (DMRG-TCCSD)Veis et al. (2016).

Even though the DMRG-TCCSD method offers a reasonably efficient treatment of both static and dynamic correlationVeis et al. , its applications to larger systems is hampered by the infavorable scaling of the CCSD part of the calculation. To remove this restriction, we decided to implement the method using a local approach. Out of many possibilities how to exploit locality Flocke and Bartlett (2004); Fedorov and Kitaura (2005); Kobayashi and Nakai (2008); Scuseria and Ayala (1999); Li et al. (2001); Stoll (1992); Kristensen et al. (2011); Li et al. (2009), we opted for pair natural orbitals (PNO) based methods, in particular, the local pair natural orbital (LPNO) approachNeese et al. (2009); Li et al. (2009); Hansen et al. (2011).

The PNOs were introduced by Edmiston and Krause and were shown to provide a compact parametrization of the virtual spaceEdmiston and Krauss (1965). Over the years, several correlation methods that made use of advantages the approach offered were developed Meyer (1971, 1974); Taylor (1981); Fink and Staemmler (1993), but the full potential of PNOs could have been unleashed only due to more recent advances in modern hardware and modern integral transformation technologies, particularly the density fitting or resolution of the identity methodologyVahtras et al. (1993). The approach was further extended to domain LPNO (DLPNO), with even more favourable scaling Riplinger and Neese (2013). Apart from single-reference methods, the LPNO and DLPNO methodologies were also succesfully applied to multireference CC techniques Demel et al. (2015); Lang et al. (2017); Brabec et al. (2018); Lang et al. (2019)

In this paper, we contribute to these efforts by implementing the LPNO version of DMRG-TCCSD. We demonstrate the properties of the method on two systems which were previously studied with its canonical implementationVeis et al. (2016, ). First, we used tetramethyleneethane (TME) as a benchmark system to estimate the amount of correlation energy possible to retrieve by LPNO approach compared to the canonical version of the method. Then, we performed a similar test on oxo-Mn(Salen) in a double-$\zeta$ basis set and subsequently performed calculation using up to a quadruple-$\zeta$ basis set, far beyond the capabilities of the canonical implementation. This way, we were able to investigate the effect of dynamic correlation in the basis sets of size previously unfeasible. Moreover, we explored the impact of active space composition on the spin state ordering of oxo-Mn(Salen) in order to shed light into previously varying claims about its ground state.

In the rest of this paper, we will use the acronym TCCSD($e$,$o$) to denote a DMRG-TCCSD calculation, in which the active space of DMRG consists of $e$ electrons in $o$ orbitals. In the same manner, LPNO-TCCSD($e$,$o$) denotes the calculation with the CC part perfomed employing the LPNO approach.

II Theory and Implementation

II.1 DMRG-based Tailored Coupled Clusters

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

[TABLE]

where 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, with $|\Phi_{0}\rangle$ being the reference wave function.

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 matrix product state (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 relations between CI and CC coefficients

[TABLE]

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

[TABLE]

analogous 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.

II.2 The LPNO Approach for DMRG-TCCSD

The first step in LPNO methods is the localization of internal orbitals and subsequent transformation of virtual space to a PNO basis. This process can be divided into three distinct steps. First, based on the MP2 calculation, pair energies and MP2 amplitudes are calculated. From these, pair density matrices are constructed for pairs of occupied orbitals and by their subsequent diagonalization PNOs are acquired.

Since approximations are made during this procedure, the accuracy of the method is controlled mainly by two cut-off parameters. The first, $T_{\text{CutPairs}}$, limits the number of occupied orbital pairs chosen for CCSD correlation treatment based on their respective pair energies. The remaining weak pairs are then treated solely at the MP2 level. The second, $T_{\text{CutPNOs}}$, determines the truncation of the PNO expansion for a given pair, based on the PNO occupation numbers.

In order to make the LPNO-TCCSD method function properly, we considered crucial to maintain the properties of the canonical TCCSD method as well as the behavior of the original LPNO-CCSD method. Just like the latter, it should maintain smooth dependance of retreived canonical correlation energy with respect to a change in the cut-off parameters . Moreover, it is necessary to perform the PNO transformation in such a way that the active orbitals in the new basis exactly align with the orbitals in the original MO basis. Failing to do so would result in a mismatch between the orbitals and imported amplitudes. For the same reason, it is required for active orbitals to pass both $T_{\text{CutPairs}}$ and $T_{\text{CutPNOs}}$ screenings.

Because we implemented LPNO-TCCSD in the spin-unrestricted version of the LPNO-CCSD code, the following steps will be described in the respective formalism Hansen et al. (2011). Also, for more compact notation we will write ${}^{\sigma\sigma^{\prime}}\!t_{ij}^{ab}=t_{i_{\sigma}j_{\sigma^{\prime}}}^{a_{\sigma}b_{\sigma^{\prime}}}$.

First, for each pair of occupied orbitals $ij$, we construct the MP2 amplitudes

[TABLE]

from exchange integrals $K_{ij}^{ab}$ and orbital energies $f_{pp}$. The $\beta\beta$ case is constructed analogously as the $\alpha\alpha$ case. Subsequently, we replace the active amplitudes, for which $\{i,j,a,b\}\subset\text{CAS}$, with the amplitudes imported from a DMRG calculation. At this point, the first cut-off parameter comes into play. For each pair, we calculate an MP2 pair energy and if it is larger than $T_{\text{CutPairs}}$, the pair is kept for further correlation treatment. To ensure that none of the active amplitudes is discarded, we circumvent this screening for all active pairs and keep them automatically.

In the next step, the PNOs are generated. First, for given pair $ij$ a pair density matrix is built from the matrix $\mathbf{T}^{ij}_{\sigma\sigma^{\prime}}$ containing the amplitudes ${}^{\sigma\sigma^{\prime}}\!t_{ij}^{ab}$

[TABLE]

Since two sets of PNOs are needed for $\alpha\beta$ pairs, we use the superscripts $(\alpha)$ and $(\beta)$ to distinguish between them. Once again, the $\beta\beta$ case is analogous to the $\alpha\alpha$ case.

For an inactive pair $\{i,j\}\not\subset\text{CAS}$, we proceed directly to diagonalization of the pair matrix to solve

[TABLE]

and obtain set of PNOs $\mathbf{d}^{ij}_{\bar{a}}$ and their respective occupation numbers $n^{ij}_{\bar{a}}$, where barred index refers to the PNO basis. This PNO expansion is then truncated based on the second cut-off parameter. Only PNOs with occupation numbers larger than $T_{\text{CutPNOs}}$ are kept and the remaining orbitals are discarded.

However, the process gets slightly more complicated for active pairs $\{i,j\}\subset\text{CAS}$. To maintain the alignment between the original active orbitals in MO basis and the active orbitals in PNO basis, it is necessary to keep their coefficients untouched during the PNO transformation. We achieve this by setting the active-external elements of a pair density matrix to zero and replacing the active-active part with an identity matrix. In order to also preserve the correct order of the active orbitals, we add "infinitesimaly" small positive numbers $\varepsilon_{a}$ to the active diagonal, for which holds that $\varepsilon_{a}>\varepsilon_{a+1}$. Thus, the resulting matrix has the block form (see Figure 1)

[TABLE]

where

[TABLE]

This way, we make sure that after solving (12) all active orbitals have the largest eigenvalues, which are in given order at the beginning, and therefore pass the $T_{\text{CutPNOs}}$ screening.

The resulting equation for singly excited amplitudes remains formally the same as the equation for the canonical method (5). On the other hand, the equation for doubly excited amplitudes (6) now becomes

[TABLE]

with the active amplitudes formally in PNO basis.

III Computational Details

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

In case of TME, we used CASPT2(6,6)/cc-pVTZ geometries for seven values of the dihedral angle from our previous work Veis et al. . The orbitals were prepared by CASSCF(6,6) calculation with the active space containing six 2pz orbitals on carbon atoms.

In case of 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 T.H. 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 the Results section. The orbitals were then split-localized using the Pipek-Mezey algorithm Pipek and Mezey (1989) in the following orbital subspaces: internal, active doubly occupied, active singly occupied and active virtual.

The orbitals for DMRG were ordered using the Fiedler method Barcza et al. (2011); Fertitta et al. (2014) combined with some manual adjustments. All DMRG runs were initialized by CI-DEAS procedure Legeza and Sólyom (2003); Szalay et al. (2015). We employed the dynamical block state selection (DBSS) procedure Legeza et al. (2003); Legeza and Sólyom (2004) to control the accuracy of the larger oxo-Mn(Salen) calculations with the truncation error criterion set to $10^{-6}$. This resulted in block dimension varying between 1000 up to 2500 block states for CAS(28,22) and up to 8200 in case of CAS(28,27). 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. Auxiliary basis sets cc-pVQZ/C and cc-pV6Z/C were used for the resolution of the identity approximation for oxo-Mn(Salen) and TME respectively Weigend et al. ; Bross et al. (2013). The default LPNO cut-off parameters were set to $T_{\text{CutPNO}}=3.33\cdot 10^{-7}$, $T_{\text{CutPairs}}=10^{-4}$ and $T_{\text{CutMKN}}=10^{-3}$ and these were used unless otherwise stated. The production runs of oxo-Mn(Salen) were performed with ORCA’s TightPNO settings i.e. the cut-off parameters set to $T_{\text{CutPNO}}=10^{-7}$, $T_{\text{CutPairs}}=10^{-5}$ and $T_{\text{CutMKN}}=10^{-4}$. For calculations which purpose was to estimate the dependance of LPNO-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 LPNO approach with reference to a DMRG-TCCSD energy calculated with the traditional TCCSD implementation.

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. . Moreover, it was already a subject of our previous study with the canonical DMRG-TCCSD method Veis et al. , so it only seems natural to use this system to test the performance of the LPNO approach to TCCSD. For this purpose, we investigate the behavior of the approximation with respect to different geometries corresponding to the rotation about its central C–C bond (see Figure 2) and different values of the cut-off parameters.

We only present results for a small active space corresponding to six electrons in six 2p orbitals. This decision followed an effort to perform the performance evaluation on three active spaces of different sizes. However, because of the small localization subspaces stemming from a small number of occupied orbitals, many orbitals remained rather delocalized. This resulted in large numbers of PNOs necessary to maintain the accuracy, even for looser cut-off parameters, which ultimately rendered the LPNO approximation useless due to its low efficiency. Therefore, we compare the performance of LPNO-TCCSD with different sized active spaces on oxo-Mn(Salen), which is better suited for this purpose.

At this point, we investigate the amount of correlation energy retrieved by the LPNO-TCCSD method compared to the canonical TCCSD method.

Figure 3 shows the dependence of the retrieved correlation energy with respect to the cut-off parameter $T_{\mathrm{CutPairs}}$, which controls the number of pairs treated by CCSD based on their estimated pair correlation energies. Both singlet and triplet state show the desired behaviour and converge towards a certain value. For settings looser than the default value of $10^{-4}$, the increasing number of pairs treated by MP2 results in overestimation of the correlation energy. As can be seen, for the triplet state the method recovers consistently about 0.1% of correlation energy less than for the singlet state. This means that the energy difference between the two states is therefore off by about 0.7 kcal/mol compared to a canonical calculation. Regarding the consistency of the approximation across different geometries, the accuracy is mostly consistent, but slight discrepancies are visible around the default value of the cut-off parameter. However, the largest difference in energies at this value is less than 0.03% canonical correlation energies, which corresponds to an energy difference smaller than 0.2 kcal/mol. Note that these differences are relevant for singlet or triplet calculation alone.

Figure 4 shows the dependence on the second studied cut-off parameter $T_{\mathrm{CutPNOs}}$. From $10^{-4}$, the recovered energy grows gradually when tightening the parameter resulting in 99.83% and 99.72% of correlation energy recovered at the default value for the singlet and triplet states respectively. Again the difference between the geometries is less than 0.03% for each state alone. For the smallest values the singlet calculation retrieves more than 100% of correlation energy, which is caused by using the default setting for $T_{\mathrm{CutPairs}}$. This means that more pairs are treated only by perturbative correction which overshoots the contributions to correlation energy.

The observed difference in accuracy for these two cases then might arise from the neglect of some of the terms in the current LPNO implementation. However, this behavior should be eradicated in the future DLPNO implementation of the method.

IV.2 oxo-Mn(Salen)

The oxo-Mn(Salen) molecule has been a subject of numerous computational studies motivated mainly by its role in catalysis of the enantioselective epoxidation of unfunctional olefins Zhang et al. (1990); Irie et al. (1990). Moreover, its closely lying singlet and triplet states pose a considerable challenge for multireference methods. Over the years, several multireference studies has 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 recently the first DMRG results with dynamic correlation treatment were presented Veis et al. (2016); Sharma et al. (2017). Our aim was to contribute to these efforts by exploring the effect of the active space and basis set dependence. on the character of the ground state. With our LPNO implementation we were able to study the effect of dynamic correlation up to the quadruple-$\zeta$ basis set. This would not be possible without the LPNO approximation, since the cc-pVQZ basis for this systems amounts to 1178 basis functions.

In order to assess the accuracy of the LPNO-TCCSD method with respect to active spaces of different size and investigate the different ground states reported at the CASSCF level, we selected two active spaces. In accordance with the study by Wouters et al. Wouters et al. (2014), the smaller CAS(28,22) consists of ten $\pi$ orbitals on equatorial conjugated rings (C, N and O atoms), five 3d orbitals on Mn, three 2p orbitals on the axial O atom and four 2p orbitals on equatorial N and O atoms forming $\sigma$ bonds with the Mn atom. On top of these, we added extra five orbitals on Mn resulting in CAS(28,27), namely 4dxy, 4dyz, and 4px, 4py and 4s, which form $\sigma^{*}$ bonds with Mn. The effect of inclusion of these particular orbitals is discussed further in the text. On top of that, we also tried to add 3p orbitals on Cl to the active space, since these were included in some of the studies Olivares-Amaya et al. (2015); Veis et al. (2016) but based on the results of entanglement analysis (one-orbital entropies) we concluded that their effect was negligible.

First, we investigate the behavior of the LPNO approximation in the smallest cc-pVDZ basis. The dependence of recovered correlation energy for both spin states and active spaces with respect to LPNO cut-off parameters is shown in graphs in Figure 6. When varying $T_{\mathrm{CutPairs}}$, the curves for all four calculations converge smoothly. The singlet and triplet curves exhibit an excellent behavior, with the errors stemming from the approximation canceling out perfectly for both states. In case of the smaller CAS(28,22), this is valid even for less conservative values.

The amount of correlation energy with respect to $T_{\mathrm{CutPNOs}}$ parameter changes smoothly towards 100% with smaller values. With the default settings, we were able to recover over 99.85% correlation energy for the smaller and over 99.78% for the larger active space. The method overshoots for $T_{\mathrm{CutPNOs}}$ over $3.33\cdot 10^{-8}$ for the same reason as with TME i.e. the fixed $T_{\mathrm{CutPairs}}$ parameter. Also the small consistent gap between singlet and triplet curves disappears with smaller cut-off for estimated pair correlation energies. As can be seen from Table 1, although the absolute energies acquired with tighter cut-offs seem to be slightly worse, these settings noticeably improve the reproduction of the canonical singlet-triplet gap.

For this system, the treatment of different spin states is balanced due to more effective cancellation of errors compared to TME. However, a difference in accuracy arises between the larger and smaller CAS, with slightly better retrieval of correlation energy in case of the latter. Similarly to TME, this can be attributed to a neglect of certain terms in the LPNO implementation of the CC method.

After the accuracy assessment, we approached the actual system. The first step was to obtain the energies of Wouters et al. Wouters et al. (2014) with our smaller CAS(28,22) considering that this study reports on the composition of the active space in detail. We successfully reproduced their results establishing the triplet ground state, but since several later studies reported significantly lower energies of both states Stein and Reiher (2016); Sharma et al. (2017) with the CAS of the same size or smaller, we wanted to investigate this further.

Therefore we included, on top of three orbitals forming the $\sigma^{*}$ bonds, additional 4dxy, 4dyz to examine how much the double-shell effect influences the stability of these states. With this active space, we obtained significantly lower energies (partially due to the larger number of variational parameters) and more importantly, the two states switched resulting in the singlet ground state. The energies can be found in Table 2.

Based on the diagrams of mutual information plotted in Figure 7, we assume that the orbitals responsible for considerable decrease in energy are mainly the antibonding $\sigma^{*}$ and partially the 4dyz orbitals. As expected, these are strongly correlated to the respective bonding orbitals. Since their effect is the same for both states of interest, the remaining 4dxy has to cause the change in the ground state. This change can be then easily explained when we consider that in singlet state the valence 3dxy orbital is fully occupied and therefore the 4dxy orbital stabilizes the singlet via the double-shell effect. This claim is supported by the fact, that mutual information between the two orbitals is fairly large in case of the singlet state. Even though its value in the diagram is of order of magnitude $10^{-2}$, out of the remaining values of the same order it is in fact one of the largest, namely 0.086, while in case of the triplet state its value is merely a fraction i.e. 0.026. Moreover, this orbital was also present in the study of Stein and Reiher Stein and Reiher (2016), who obtained the singlet ground state as well, with very similar energy difference between the two states.

When dynamic correlation is added to the calculation by the LPNO-TCCSD method, a significant shift in energies occurs. Looking at the energies of the triplet state we can observe that the extra orbitals in the larger CAS seems to contribute solely to the dynamic correlation at the CASSCF level, since the LPNO-TCCSD energies are virtually the same for both active spaces. Singlet energies, on the other hand, differ depending on the active space, which can be attributed to the effect of aforementioned static correlation from the inclusion of 4dxy orbital. Nonetheless, irrespective of the active space, we found out the ground state to be a triplet with slightly larger singlet-triplet gap for the smaller CAS. This is in agreement with our previously published results Veis et al. (2016), in which we used canonical TCCSD, but only with orbitals coming from an ROHF calculation.

Furthermore, these results remain consistent with larger basis sets. Going from DZ to QZ basis at CASSCF level, the ground state of the particular CAS slightly stabilizes. Although no significant change occurs for the singlet-triplet gap by including dynamic correlation in the smaller CAS, the size of the gap decreases for larger CAS. This change roughly corresponds to what was observed at CASSCF level.

Comparing these results to those published by Sharma et al. Sharma et al. (2017), we conclude that employing larger basis set does not qualitatively affect the ground state. Since we have shown that the system is quite sensitive to composition of the active space and apart from DZ, their energies agree very reasonably with ours, we suppose there might have been some inconsistency in the active space for DMRG-CASSCF in the smallest studied basis set. Especially our best result LPNO-TCCSD(28,27) is in excellent agreement with the best available NEVPT2(28,22) in cc-pVQZ basis.

Regarding the future endeavors with this system, we would suggest the geometry reoptimization. Although the used geometry serves as a useful benchmark for the reason that it is frequently used by several different groups, since the introduction of dynamic correlation, the optimization at CASSCF(10,10)/6-31G* appears to be fairly inadequate.

V Conclusions

We introduced a new version of DMRG-TCCSD method, which employs the local pair natural orbital approach. The method has been implemented in ORCA presently at the singles and doubles level.

We performed accuracy assessment of the method employing two systems, which were previously studied by the canonical TCCSD method. Regarding tetramethyleneethane, we were able to retrieve over 99.7% for triplet and over 99.8% for singlet state, while using the default settings of cut-off parameters. For oxo-Mn(Salen), the amount of retrieved correlation was dependent on the size of the active space used, ranging from 99.6% for the larger CAS(28,27) to 99.8% for smaller CAS(28,22). Despite this dependence, an excellent agreement was achieved concerning the accuracy between the spin states. Using the default settings resulted in singlet-triplet gap being off by 0.6 kcal/mol and with tighter cut-offs only 0.04 kcal/mol compared to canonical calculation.

Furthermore, we investigated previously unexplored problem of varying reports of different ground states of oxo-Mn(Salen) at CASSCF level. We discovered that this inconsistency most likely originates from the composition of CAS. In particular, we found out that the orbital responsible for stabilizing the singlet state is the 4dxy orbital via the double-shell effect. However, by employing the dynamic correlation treatment with LPNO-TCCSD, regardless of the basis set the ground state was unambiguously found to be a triplet.

Regarding the future of the method, we would like to implement the DLPNO version of the TCCSD, which we hope to further enhance capabilities of the method and also include the perturbative triples correction for even more accurate results.

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. 16-12052S and 18-18940Y), Czech Ministry of Education, Youth and Sports (Project No. LTAUSA17033), Charles University (Project GA UK No. 376217) and by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project „IT4Innovations National Supercomputing Center – LM2015070“.