Multilevel CC2 and CCSD in Reduced Orbital Spaces: Electronic Excitations in Large Molecular Systems

We present efficient implementations of the multilevel CC2 (MLCC2) and multilevel CCSD (MLCCSD) models. As the system size increases, MLCC2 and MLCCSD exhibit the scaling of the lower-level coupled cluster model. To treat large systems, we combine MLCC2 and MLCCSD with a reduced-space approach in which the multilevel coupled cluster calculation is performed in a significantly truncated molecular orbital basis. The truncation scheme is based on the selection of an active region of the molecular system and the subsequent construction of localized Hartree–Fock orbitals. These orbitals are used in the multilevel coupled cluster calculation. The electron repulsion integrals are Cholesky decomposed using a screening protocol that guarantees accuracy in the truncated molecular orbital basis and reduces computational cost. The Cholesky factors are constructed directly in the truncated basis, ensuring low storage requirements. Systems for which Hartree–Fock is too expensive can be treated by using a multilevel Hartree–Fock reference. With the reduced-space approach, we can handle systems with more than a thousand atoms. This is demonstrated for paranitroaniline in aqueous solution.


■ INTRODUCTION
The scaling properties of the coupled cluster hierarchy of methods severely limits the systems for which it is applicable. 1 The methods have polynomial computational scaling, N ( ) n , where N is a measure of system size and n increases with accuracy of the method. Memory and disk space requirements also increase significantly as one moves up through the hierarchy.
The development of reduced cost and reduced scaling coupled cluster methods has been an active topic for decades. Arguably, the most popular approach has emerged from the work of Pulay and Saebø. 2,3 They demonstrated that dynamical electronic correlation could be compactly described using localized orbitals rather than canonical orbitals; specifically, they used localized occupied molecular orbitals (MOs), such as Boys 4 or Pipek-Mezey 5 orbitals, and projected atomic orbitals 2,3 (PAOs) to span the virtual space. Their local correlation approach was later applied to coupled cluster theory by Hampel, Werner, and Schuẗz. 6,7 Other local coupled cluster methods include the local pair natural orbital 8,9 and the orbital-specific-virtual 10 coupled cluster methods. Whereas the success of these local coupled cluster methods in the description of the ground state correlation energy is indisputable, their extension to excited states has turned out to be more complicated. 11−16 A different approach originates from the multireference coupled cluster method of Oliphant and Adamowicz. 17−19 While introduced to describe multireference character, the method is formulated in the framework of single reference coupled cluster theory. An active orbital space is used, and higher order excitation operators (e.g., triple or quadruple excitations) are included with some indices restricted to the active space. Koḧn and Olsen 20 recognized that the method could be used to reduce the cost for single reference systems, and this was further demonstrated by Kaĺlay and Rolik. 21 The multilevel coupled cluster (MLCC) approach, introduced by Myhre et al., 22−24 is closely related to this active space approach.
In MLCC, the goal is to accurately describe excitation energies and other intensive properties, rather than extensive properties such as correlation energies. This is done by restricting the higher order excitation operators to excite within an active orbital space. For example, in the multilevel CCSD (MLCCSD) 22,25 method, the double excitation operator is restricted to excite out of active occupied orbitals and into active virtual orbitals. In this work, we demonstrate the available computational savings of the multilevel CC2 (MLCC2) and MLCCSD models introduced in ref 25; for sufficiently large inactive spaces, we show that the cost is dominated by the lower-level method. This has previously been demonstrated for multilevel CC3 (MLCC3) by Myhre et al. 24 The scaling of the lower-level model cannot, however, be avoided. Therefore, in order to use these methods for large systems, they must be combined with other multilevel or multiscale approaches. For instance, MLCC could be used within a QM/MM 26,27 framework or with the polarizable continnum model. 28,29 Here, we have chosen to perform MLCC calculations in a significantly truncated MO basis. The truncation of the MO basis in coupled cluster calculations is used routinely. For example, the frozen core approximation falls into this category, and there are several examples of truncation of natural orbitals, both of the virtual and occupied spaces. 30−35 The LoFEx 36,37 and CorNFLEx 38 approaches are also notable reduced space coupled cluster approaches that target accuracy in the excited states. In these approaches, a mixed orbital basis consisting of natural transition orbitals (NTOs) 39−41 and localized orbitals is used. The active space is expanded until the excitation energies have converged to within a predefined threshold. One drawback of LoFEx and CorNFLEx is that they are state specific methods, that is, several subsequent calculations with different truncated MO bases must be performed to obtain a set of excitation energies. As a consequence, the calculation of transition moments between excited states is complicated by the fact that the states are nonorthogonal and interacting. An orbital selection procedure similar to that of LoFEx and CorNFLEx has also been used for reduced scaling second-order algebraic diagrammatic construction (ADC(2)) calculations by Mester et al. 35,42 Here, we use a truncation scheme for the MOs where semilocalized Hartree−Fock orbitals (virtual and occupied) are constructed and used to calculate localized intensive properties in large molecular systems. When the region of interest is sufficiently small compared to the full system, the number of MOs (n MO ) in the coupled cluster calculation is much smaller than the number of atomic orbitals (AOs) denoted N AO . This reduced space approach has previously been used with standard coupled cluster models, 43,44 and a very similar approach has been used together with local coupled cluster models. 45 For sufficiently large systems, the cost of Hartree− Fock can become a limiting factor. When this is the case, we handle it by combining the reduced space MLCC approach with a multilevel Hartree−Fock 46,47 (MLHF) reference wave function.
The MLCC2 and MLCCSD implementations are based on Cholesky decomposed electron repulsion integrals. 48,49 We use the two-step Cholesky decomposition algorithm introduced in ref 50. In this algorithm, the Cholesky basis and the Cholesky vectors are determined in two separate steps. We have implemented a direct construction of the Cholesky vectors in the truncated MO basis. This reduces the memory requirement of the vectors from N ( ) AO 3 to N n ( ) AO MO 2 , making it possible to efficiently perform reduced space calculations on systems with several thousands of basis functions. We use n as a measure of the size of the active space, which does not scale with the system. It should be noted that storage of the Cholesky vectors in the AO basis, albeit temporary, can only be avoided in a decomposition algorithm that determines the Cholesky basis and the Cholesky vectors in separate steps. In the Cholesky decomposition, we also use the MO screening procedure that was introduced in ref 50. This MO screening leads to fewer Cholesky vectors, further reducing the memory requirement of the Cholesky vectors to n ( ) 3 .

■ THEORY
In coupled cluster theory, the wave function is defined as where |HF⟩ is the Hartree−Fock reference, X are the cluster operator, x μ are cluster amplitudes, and τ μ are the excitation operators. The standard models within the coupled cluster hierarchy are defined by restricting X to include the excitation operators up to a certain order. In the CCn models, such as CC2 51 and CC3, 52 the nth order excitations are treated perturbatively.
In the following, the indices α, β, γ, ... and p, q, r, ... refer to spatial atomic and molecular orbitals, respectively, and the indices i, j, k, ... and a, b, c, ... refer to occupied and virtual orbitals. The total number of occupied and virtual orbitals are denoted by N o and N v , respectively, and the number of active occupied and active virtual orbitals are denoted by n o a and n v a .

■ MULTILEVEL CC2 AND CCSD
The MLCC2 cluster operator is given by where the single excitation operator, X 1 , is unrestricted, that is, defined for all orbitals, whereas the double excitation operator, S 2 , is restricted to excite within an active orbital space. As in standard CC2, S 2 is treated perturbatively. The MLCC2 ground state equations are given by where Ĥis the X 1 -transformed Hamiltonian and F is the Fock operator. The doubles projection space, {⟨μ 2 S |}, is associated with S 2 . Except for the restriction of S 2 and the projection space, these equations are equivalent to the standard CC2 ground state equations. The MLCC2 equations are solved in a basis where the active-active blocks of the occupied-occupied and virtual-virtual Fock matrices are diagonal. In this semicanonical basis, eq 4 can be solved analytically for the S 2 amplitudes in each iteration. The double amplitudes are inserted into eq 3, which is solved with a DIIS-accelerated 53 quasi-Newton solver 54 to obtain X 1 . The MLCC2 equations are formulated in terms of the Cholesky vectors in the X 1 -basis. See Appendix A for detailed expressions.
If we consider a fixed active space, the overall scaling of the MLCC2 ground state equations is N ( ) 4 : the X 1 -transformation of the Cholesky vectors scales as N ( ) 4 , as does the computation of the Fock matrix in the X 1 -basis and the correlation energy. The construction of Ω scales as N ( ) 2 . The MLCC2 excitation energies are determined as the eigenvalues of the Jacobian matrix, where T 2 τ ν is a double excitation included in T 2 . In addition to the terms of the transformation by A MLCC2 that enter the transformation by A MLCCSD , there are terms which scale as (1) , N ( ) , and N ( ) 2 . Integral construction for the different terms scales, depending on the number of restricted indices, as N ( ) , N ( ) 2 , or N ( ) 3 . See Appendix A for detailed expressions.
Partitioning the Orbital Space. Selecting the active orbital space for a multilevel coupled cluster calculation is not trivial. Generally, the canonical Hartree−Fock orbitals must be transformedthrough occupied-occupied and virtual-virtual rotationsto an orbital basis that can be intuitively partitioned. To determine the type of orbitals to use, both the targeted property and the system must be considered. There are two main approaches to select the active spaces. If the property of interest is adequately described at a lower level of theory, then the information from that lower level can be exploited to partition the orbitals. An example is the use of correlated NTOs (CNTOs). 25,41 If the property of interest is spatially localized, then localized or semilocalized orbitals can be applied. For instance, Cholesky orbitals 44,56 have been used in multilevel coupled cluster calculations by Myhre et al. 23,24,57 The CNTOs are constructed using excitation vectors, R, from a lower-level calculation. The matrices are diagonalized; the matrices that diagonalize M and N are the transformation matrices of the occupied and virtual orbitals, respectively. From eqs 11 and 12 it may seem that the lower level method must include double excitation amplitudes in its parametrization. However, CNTOs can be generated from CCS excitation vectors by constructing approximate double excitation vectors: Here, ω CCS is the CCS excitation energy, and ϵ ij ab = ϵ a + ϵ b − ϵ i −ϵ j , where the ϵ q are orbital energies. The integrals g ̅ aibj are defined as where g pqrs = ∑ J L pq J L rs J is the electronic repulsion integrals in the MO basis and I I I ij ab aibj aibj bjai = + (I aibj are elements of a rank-4 tensor). Equations 13 and 14 were suggested by Baudin and Kristensen 38 and are based on CIS(D). 58 In our previous work, we have found that the CNTOs obtained from a CCS calculation (using eqs 13 and 14) perform well, considering accuracy and cost, compared to CNTOs from a CC2 calculation. 25 It should be noted, however, that these orbitals are not expected to perform well for states dominated by double excitations with respect to the reference.
The active space is selected by considering the eigenvalues of M and N: active orbitals result from the eigenvectors corresponding to the largest eigenvalues. In this work, we either explicitly select the number of active occupied and active virtual orbitals (n o a and n v a ) or we select n o a and let the number of active virtual orbitals be determined from the total fraction of virtual to occupied orbitals; that is, Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Alternatively, one can use the selection criterion given in ref 41. The latter approach is more suitable for production calculations; on the other hand, eq 15 is convenient for testing the models. Several excited states can be considered simultaneously by diagonalizing sums of M and N matrices generated from the individual excitation vectors (eqs 11 and 12). 25 Cholesky orbitals 44,56 are obtained by a restricted Cholesky decomposition of the Hartree−Fock densities (occupied and virtual); the pivots of the decomposition procedure are restricted to correspond to AOs centered on active atoms.
As an alternative to Cholesky orbitals for the virtual space, one can use projected atomic orbitals 3 (PAOs). To construct the PAOs, the occupied orbitals are projected out of the AOs, {χ α }, centered on the active atoms: (16) Here, C is the orbital coefficient matrix, D is the idempotent Hartree−Fock density, and S is the AO overlap matrix. The orbital coefficient matrix for the active PAOs is therefore where S′ is rectangular and contains the columns of S which correspond to AOs on active atomic centers. The PAOs are nonorthogonal and linearly dependent. To remove linear dependence and orthonormalize the active virtual orbitals, we use the Loẅdin canonical orthonormalization procedure. 59 The inactive virtual orbitals are obtained in a similar way: the occupied orbitals, as well as the active virtual orbitals, are projected out of the AOs and the resulting orbitals are finally orthonormalized. After the orbitals have been partitionedregardless of which orbitals are usedwe transform to the semicanonical MO basis that is used in MLCC2 and MLCCSD calculations. This transformation involves block-diagonalizing the virtualvirtual and occupied-occupied Fock matrices such that the active-active and inactive-inactive blocks become diagonal.
Reduced Space Multilevel Coupled Cluster. Recall that MLCC methods exhibit the scaling of the lower-level coupled cluster model. To overcome this limitation, we apply a reduced space approach in which only a subregion of the molecule is described at the coupled cluster level. The orbitals in this subregion are divided into active and inactive sets for the MLCC calculation. The rationale behind this approach is that localized intensive properties can be described by using accurate and expensive methods only for the region of interest. In particular, it is assumed that the effect of the more distant environment is sufficiently well captured through contributions to the Fock matrix. A few numerical results 43,44 indicate that excitation energies can be described accurately with this f rozen Hartree−Fock approach. However, a comprehensive study has not yet been published.
To perform reduced space MLCC calculations, we must first choose the region of the molecular system to be treated with MLCC. After the Hartree−Fock calculation, localized occupied and virtual orbitals are constructed for the active region. Any localization procedure can be employed; however, we use Cholesky orbitals for the occupied space and PAOs for the virtual space. This set of orbitals enters the MLCC calculation. The remaining occupied orbitals enter the equations through their contributions to the Fock matrix, Here, N o f is the number of frozen occupied orbitals and the index I denotes a frozen occupied orbital. The multilevel coupled cluster calculation now has n MO ≪ N AO , but the procedure is otherwise unchanged: the reduced set of MOs is partitioned into active and inactive sets and the MLCC equations are solved. We write n MO (with lower case n) to indicate that the number of MOs does not scale with the system in such calculations.
A multilevel Hartree−Fock 46,47 (MLHF) reference can also be used. As in MLCC, one first determines the active orbitals: a set of active atoms is selected, and the active occupied orbitals are obtained through a partial limited Cholesky decomposition of the initial idempotent density; PAOs can be used to determine the active virtual orbitals. Only the active orbitals are optimized in the Roothan−Hall procedure, which is performed in the MO basis. 47 The inactive orbitals enter the optimization through an ef fective Fock matrix that assumes the same form as in eq 17. The inactive two-electron contribution (F f ) is only computed once at the beginning of the calculation and is subsequently transformed to the updated MO basis in every iteration (for details, see ref 46).
The reduced space MLCC approach relies on the definition of levels of active regions of the system, see Figure 1. We must first select which atoms are active in the Hartree−Fock (HF) calculation. If all atoms are active, we have a standard HF reference. Second, we must determine which atoms enter the MLCC calculation. Lastly, if we use Cholesky/PAOs to partition the orbitals in the MLCC calculation, we must determine which atoms should be treated with the higher level coupled cluster method. This is not necessary when CNTOs are used. Note that the active atom sets for higher level methods are contained within the active atom sets of lower level methods (see Figure 1).
Since these methods rely on selecting active regions, they are especially well suited for solute/solvent systems. They may also be used for other large systems where the region of interest is known.
Integral Handling for Reduced Space Calculations. When n MO ≪ N AO and N AO is large, as is often the case in reduced space calculations, the electron repulsion integrals must be handled carefully to avoid prohibitive scaling with total system size. In the AO basis, the Cholesky vectors, L J , have a storage requirement of N ( ) AO 3 ; as demonstrated by Røeggen and Wisløff-Nilssen, 60 the number of Cholesky vectors, N J , is approximately MN AO when a decomposition threshold of 10 −M is used. For example, with a loose decomposition threshold of 10 −2 , about 28 TB of memory is needed to store the Cholesky vectors of a molecular system with 12000 AOsassuming double precision and no screening.
We have previously suggested a two-step Cholesky decomposition algorithm 50 in which the Cholesky basis (i.e., the set of pivots), , is determined in the first step. The Cholesky vectors are constructed in the second step through an RI-like expression, 18) where the matrix Q is the Cholesky factor of the matrix S KL = (K|L) for K L , ∈ . This two-step algorithm makes it possible to directly construct the Cholesky vectors in the MO basis: We emphasize that it is not possible to avoid storing the AO Cholesky vectors with a one-step Cholesky decomposition of the AO electron repulsion integral matrix. Alternatively, the MO electron repulsion integrals can be constructed from the AO integrals. To reduce the scaling, one can combine screening on the AO integrals and the MO-coefficients.
Below we outline an algorithm to construct and store the vectors directly in the MO basis (see Algorithm 1). This is done after the elements of the basis K ∈ have been determined, S has been constructed and decomposed, and Q has been inverted. When the MO Cholesky factor, L, is too large to store in memory, L pq J is constructed for a maximum number of p indices (resulting in several batches, P 1 , P 2 , ..., P n ). The direct construction of the Cholesky vectors in the MO basis reduces the storage requirement to N n ( ) AO MO 2 . Note that this is linear, rather than cubic, in N AO .
Algorithm 1 is designed to avoid the IO operations involved in temporary storage and reordering of the intermediate X.
Alternatively, X can be constructed and stored on disk before L is constructed in batches over p or q. With the latter approach, the integrals (αβ|K) are never recalculated. It should be noted, however, that when n MO ≪ N AO , batching over p is typically not necessary.
The number of Cholesky vectors, N J , canthrough a method-specific screeningbe made to scale with n MO rather than N AO . Consequently, the storage requirements become n ( ) MO 3 . Method-specific decompositions were first considered by Boman et al. 61 We use the active space screening given in ref 50. In a given iteration of the Cholesky decomposition procedure, the next element of the basis is determined by considering the updated diagonal of the integral matrix Here, the sum is over the current elements of the basis. In the standard decomposition algorithm, the next element of the basis is selected as the K = αβ corresponding to the largest element of D. The decomposition procedure is terminated when where τ is the decomposition threshold. In the spirit of method specific Cholesky decomposition, 61 one can consider the Cholesky decomposition of the matrix with elements The positive semidefiniteness of M follows directly from the positive semidefiniteness of g. The diagonal of M,  (25) and where C a is the MO coefficient matrix of the reduced space MLCC calculation. We can modify the procedure to determine the Cholesky basis. The selection and termination criteria are changed by considering the screened diagonal  (27) where the errors Δ αβpq,γδrs are less than τ.
Finally, let us briefly consider the computational scaling of the decomposition procedure. Except for the initial integral cutoff screening, which scales as N ( ) AO 2 in our implementation, the MO-screened decomposition algorithm scales as n ( ) MO 3 . The prescreening step can be implemented with a lower scaling; however, this step is not time-limiting in any of the reported calculations.

■ RESULTS AND DISCUSSION
The MLCC2 and MLCCSD methods have been implemented in a development version of the e T program. 43 The following thresholds are applied, unless otherwise stated: the Hartree− Fock equations are solved to within a gradient threshold of 10 −8 ; the Cholesky decomposition threshold is 10 −3 ; the coupled cluster amplitude equations are solved such that |Ω| < 10 −6 ; the excited state equations are solved to within a residual threshold of 10 −4 ; and occupied Cholesky orbitals are constructed using a threshold of 10 −2 on the pivots. The frozen  n o a and n v a are the number of active occupied and virtual orbitals, and ω is the lowest excitation energy. The wall times to solve the ground and excited state equations (t gs and t es ) and to construct the CNTOs (t CNTO ) are also given. The calculations were performed on two Intel Xeon E5-2699 v4 processors using 44 threads. The calculations were performed with 500 GB memory available. Peak memory usage (PMU) is given in GB. n o a and n v a are the number of active occupied and virtual orbitals, and ω i is the ith excitation energy. The wall times to solve the ground and excited state equations (t gs and t es ) and to construct the CNTOs (t CNTO ) are also given. The calculations were performed on two Intel Xeon Gold 6138 processors with 40 threads and 355 GB memory available. Peak memory usage (PMU) is given in GB. core approximation is used throughout. All geometries are available from ref 62.
Performance and Scaling. The MLCC2 and MLCCSD methods can be used to obtain excitation energies of CC2 and CCSD quality, at significantly reduced cost. This is demonstrated for rifampicin and adenosine, see Figure 2. For rifampicin, the lowest excitation energy is calculated at the MLCC2/aug-cc-pVDZ and CC2/aug-cc-pVDZ levels of theory. For adenosine, the three lowest excitation energies are calculated at the MLCCSD/aug-cc-pVDZ and CCSD/augcc-pVDZ levels of theory. We have used CNTOs to partition the orbitals. The results are given in Tables 1 and 2, respectively. These show that the error in the MLCCSD and MLCC2 excitation energies with respect to CC2 and CCSD is smaller than the expected error of CC2 and CCSD. 63,64 Furthermore, the cost is drastically reduced in all cases.
In Tables 1−3, we have given the available memory and peak memory used in these calculations. Note that the calculations may be performed with less memory since the models are implemented with batching for the memory intensive terms to use no more storage than the specified available memory.
The lowest excitation energy of rifampicin was also calculated with MLCCSD/aug-cc-pVDZ, see Table 3. Since the system has 1806 MOs, a full CCSD calculation would be demanding; therefore, we do not present a reference CCSD calculation. However, the variation of the excitation energy is less than 0.05 eV for the different active spaces and can therefore be considered converged. In our experience, MLCCSD excitation energies converge smoothly to the CCSD values. 25 Note that the MLCC2 and MLCCSD timings cannot be compared as the calculations were performed on different processors.
To demonstrate the scaling properties, we consider a system of PNA and water molecules. The size of the active space is fixedwith 36 occupied and 247 virtual orbitalsand the system size is increased by adding water molecules (see Figure  3). We use the aug-cc-pVDZ basis set.
In Figure 4, we show the overall wall times of the Hartree− Fock calculation, the CNTO construction, and the MLCC ground and excited state calculations. The steep N ( ) 5 scaling of the CNTO construction is apparent: for the largest system, it is the most expensive step. The ground and excited state MLCC equations scale as N ( ) 4 ; however, for the larger systems we have considered, the Hartree−Fock calculation is seen to be more expensive. This must be understood in the context of system size and the use of an augmented basis set. For sufficiently large inactive spaces, the N ( ) 4 terms of MLCC2 and MLCCSD will become more expensive than Hartree−Fock.
In Figure 5, we present a timing breakdown of an iteration to solve the MLCC ground state equations. The iteration is dominated by the N ( ) 4 step to construct the X 1 -transformed Cholesky vectors. The calculation of the energy, and the necessary blocks of the Fock matrix in the X 1 -basis, also scale as N ( ) 4 , but the prefactor is lower for these operations. The construction of the Ω-vector scales as N ( ) 2 . In MLCCSD, the Ω-vector contains additional contractions, compared to MLCC2, that scale as (1) or N ( ) (see Appendix A). n o a and n v a are the number of active occupied and virtual orbitals, and ω is the lowest excitation energy. The wall times to solve the ground and excited state equations (t gs and t es ) and to construct the CNTOs (t CNTO ) are also given. The calculations were performed on two Intel Xeon Gold 6138 processors with 40 threads and 355 GB memory available. Peak memory usage (PMU) is given in GB.  In Figure 6, we plot the wall time of the Jacobian matrix transformation together with the time spent on terms that arise at the CCS, CC2, and CCSD level of theory. The CCS terms scale more steeply ( N ( ) 4 ), and for MLCC2, we see that these terms dominate when the inactive space is sufficiently large. For MLCCSD, the CCS terms are significant, but they do not dominate for any of the systems.
Reduced Space Calculations. We now consider a larger PNA-in-water system. The geometry is extracted from a single snapshot of a molecular dynamics simulation taken from ref 65. The PNA-in-water system is restricted to a sphere centered on PNA with a 15 Å radius and includes 499 water molecules, see Figure 7.
To assess the accuracy of the MO screening procedure of eqs 25 and 26, we consider the lowest MLCCSD-in-HF excitation energy of the system, which corresponds to a charge transfer process in PNA. We compare the MO-screened Cholesky decomposition with the standard Cholesky decomposition. Note that we use the partitioned Cholesky decomposition (PCD) algorithm, described in ref 50, with two batches. In these MLCCSD calculations, the atoms within a sphere of 5 Å are included in the MLCC region (r CCS = 5 Å) and the atoms within a sphere of radius 3.5 Å are defined as active at the CCSD level of theory (r CCSD = 3.5 Å). The orbitals are partitioned with the Cholesky/PAO approach. For the CCSD/CCS/HF levels of theory, we use the aug-cc-pVDZ/cc-pVDZ/STO-3G basis sets. The total number of basis functions is 3971, and in the MLCCSD-in-HF calculations, we have n o CCSD = 90, n v CCSD = 287, n o CCS = 57, and n v CCS = 219, that is, n MO = 653. The results are given in Table 4.
The MO screening yields significantly fewer Cholesky vectors without introducing large errors in the excitation energies. As expected, the number of Cholesky vectors, N J , is seen to be on the same order of magnitude as N AO and n MO for the standard and MO-screened decomposition algorithms, respectively. Fewer Cholesky vectors reduces the cost of the coupled cluster calculation, for which the Cholesky vectors are either used to construct the integrals or applied directly in Cholesky vector-based algorithms. Moreover, the decomposition time is reduced when the MO screening is employed; for instance, with a threshold of 10 −4 , the decomposition time was 1161 s without screening and 286 s with screening. In any case, the decomposition time is not a bottleneck in any of these calculations.
The largest error in the approximated AO integral matrix, ϵ, is also given in Table 4. For standard PCD, the errors are comparable to the decomposition threshold. With MO screening, ϵ is large because AO integrals that do not contribute to the MO integrals are not described by the Cholesky vectors. Without MO screening, a Cholesky decomposition threshold of 10 −2 or 10 −3 is typically Figure 5. Timing breakdown of the MLCC2 (top) and MLCCSD (bottom) ground state iteration for PNA and water. N w is the number of water molecules, t Energy is the time to compute the MLCC correlation energy, t F is the time to construct the necessary blocks of the Fock matrix in the X 1 -basis, t Ω is the time to construct the Ωvector, and t X 1 is the time to X 1 -transform the Cholesky vectors. The calculations were performed on two Intel Xeon E5−2699 v4 processors using 44 threads and with 1.4 TB memory available.  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article sufficient. 50 For MLCC2 or MLCCSD in a reduced space calculation, the MO screening can be used and a threshold of 10 −4 seems suitable. In the calculation with MO screening and a threshold of 10 −2 , the MLCCSD calculation did not converge.
We have also performed MLCC calculations on the PNA-inwater system in Figure 7 with larger basis sets. In Table 5, we present timings for MLCC2-in-HF/MLHF and MLCCSD-in-HF/MLHF calculations with r CCS = 6.0 Å and r CC2/CCSD = 3.5 Å. The aug-cc-pVDZ basis is used for all atoms included in the CC active region, and cc-pVDZ is used on the remaining atoms. In total, there are 12669 AOs and 1498 MOs in the coupled cluster calculation. The Cholesky decomposition is performed with MO screening using a threshold of 10 −4 . For the reference calculations, a gradient threshold of 10 −6 is used.
Comparing Tables 4 and 5, we see that the MLCCSD-in-HF excitation energies do not change significantly with a larger basis and an increased r CCS . For the calculations presented in Table 5, the reference calculation is the most expensive step. Since the active region of the MLHF calculation is large (10.0 Å), we do not obtain large savings using an MLHF reference. However, this can be achieved by reducing r HF . Furthermore, MLHF is applicable for systems where standard Hartree−Fock is not computationally feasible. The CC2-in-HF calculation for this system, with a CC2 radius of 6 Å, yields ω = 3.732 eV. Hence, the error of using MLCC2, compared to CC2, is approximately 0.1 eV. The effect of extending the CCS radius to r CCS = 8.0 Å is to increase the excitation energy by 0.003 eV to ω = 3.824 eV.
Solvation effects can be estimated by performing calculations on a series of snapshots from a molecular mechanics simulation, for instance using the QM/MM approach for the individual snapshots, such as in ref 65. The calculations in this paper demonstrate that a fully quantum mechanical approachMLCC-in-HF and MLCC-in-MLHFcan be used to determine such solvation effects. For the former of these approaches, the Hartree−Fock calculation is likely to be the time limiting step.
The MLCC-in-HF and MLCC-in-MLHF approaches are not only applicable to solute−solvent systems. They can also be used for large molecules. As a proof of concept, we present MLCC-in-HF calculations for the lowest excitation in rifampicin in Table 6. In Figure 8, we have plotted the NTOs from a CCS/aug-cc-pVDZ calculation. The excitation is seen to be located in a subregion of the molecule. It can therefore be treated with CC-in-HF or MLCC-in-HF. In Figure 8, we have also plotted the Hartree−Fock density of the active occupied orbitals treated with CC or MLCC. The active atoms in the CC and MLCC calculations were selected by hand by inspecting the NTOs. CNTOs were used to partition the orbitals in the case of MLCC2/MLCCSD. The shift observed by going from X to X-in-HF is about 0.1 eV in all the presented calculations. It should be noted that this system is too small to be suitable for CC-in-HF and MLCC-in-HF, but is chosen because the reference CC2 calculations are available. The MLCC2 or MLCCSD methods are preferable for systems of this size since, as can be seen from Tables 1 and 3, these calculations can be performed with ease. The PCD algorithm is used. The threshold, τ, the number of Cholesky vectors, N J , and the largest error in the approximated matrix in the AO basis, ϵ, are given. There are 3971 basis functions. The atoms within a radius of 6 Å are included in the MLCC calculation, and the atoms within a radius of 3.5 Å are treated with the higher level coupled cluster method (CC2 or CCSD). In the MLHF reference calculation, the atoms within a radius of 10 Å are active. The aug-cc-pVDZ basis is used on the atoms that are included in the MLCC calculation, and cc-pVDZ is used on the remaining atoms. The wall times for the reference calculation (t ref ) and the MLCC calculation (t MLCC ) are also given. The calculations were performed on two Intel Xeon Gold 6152 processors with 44 threads and 1.4 TB memory available. The peak memory usage (PMU) is given in GB. a The aug-cc-pVDZ basis is used, unless otherwise stated. The active atoms in the CC and MLCC calculations were selected by hand (see Figure 8). CNTOs were used to partition the orbitals in the MLCC calculations. b aug-cc-pVTZ on active atoms. . With an additional screening during the Cholesky decomposition algorithm, we further reduce this requirement to n ( ) MO 3 , making the storage requirement independent of the size of the environment. Exploiting the Cholesky factorization in this manner, we can handle systems with several thousand basis functions using existing MLCC implementations. The MLCCin-HF/MLHF framework is therefore suited to accurately model solvation effects on intensive properties on the solute. It can also be used for chromophores in biomolecules.

■ APPENDIX A
We use the following notation for the MLCC2 and MLCCSD equations: indices a, b, c, ... denote active virtual orbitals; A, B, C, ..., unrestricted virtual orbitals; i, j, k, ..., active occupied orbitals; I, J, ..., unrestricted occupied orbitals; p, q, r, ..., general active orbitals; and P, Q, R, ..., general unrestricted orbitals. The index K is used to denote Cholesky vectors. We also define n o and n v as the number of active occupied and active virtual orbitals, respectively, and N o and N v as the total number of occupied and virtual orbitals. We also use N MO and N J for the number of MOs and Cholesky vectors, respectively. We have adopted the Einstein notation with implicit summation over repeated indices. In the following screening considerations, we assume a fixed active space and expanding inactive space.
The electron repulsion integrals are Cholesky decomposed, The Cholesky vectors are stored in both the MO and the X 1transformed basis. In general, X 1 -transformed quantities are denoted with tilde; for example, ( 29) When no indices are restricted to the active space, the construction of g or g̃from L or L, respectively, is an N ( )

operation.
The Ground State Equations and the Correlation Energy. Solving the projected coupled cluster equations, entails the iterative construction of the Ω-vector, the iterative construction of the Fock matrix in the X 1 -transformed basis (F), and the calculation of the correlation energy. The correlation energy is computed in every iteration; however, this is not necessary as convergence can be determined purely from the norm of Ω. The Fock matrix in the X 1 -transformed basis is given by ( 31) and its construction, in terms of the Cholesky vectors, scales as N MO 2 N O N J ( N ( ) 4 ). However, depending on the coupled cluster model, only certain subblocks of F̃are needed to solve eq 30.
In MLCCSD and MLCC2, the correlation energy is given by where we have introduced L PQRS = 2g PQRS − g PSRQ . The last term in eq 32 is restricted to the active space. The first term is calculated according to which scales as N V N O 2 N J ( N ( ) 4 ), avoiding the N ( ) 5 integral constructions.

The MLCC2 Ω-vector
In MLCC2, the cluster operator is given by X X S 1 2 = + The Ω-vector becomes  where F PQ are elements of the Fock matrix. The Ω-vector is coded as . All orbital indices are restricted to the active space and only the integral construction scales with the system (linear scaling, N ( ) ).

Jacobian Transformation
The linear transformation by the Jacobian matrix, Ac σ = (43) must be calculated in order to obtain excitation energies in coupled cluster theory.

MLCC2 Jacobian Transformation
The MLCC2 Jacobian matrix is given by δ ϵ μ μ ν in the semicanonical basis, where ϵ aibj = F aa + F bb − F ii − F jj .
The terms of the singles part of the transformed vector are The contractions in eq 51 scale as either (0) , N ( ) , or N ( ) 2 . Additionally, the integral constructions scale as N ( ) , N ( ) 2 , or N ( ) 3 . The intermediates are calculated once before the iterative loop, and are defined as