Analysis of environment response effects on excitation energies within subsystem-based time-dependent density-functional theory

We analyze the ability of subsystem time-dependent density-functional theory (sTDDFT) to describe environmental response effects. To this end, we utilize the recently proposed “ exact ” version of sTDDFT relying on projection-based embedding (PbE), which so far was applied only for the special case of two subsystems. We confirm that PbE-sTDDFT in combination with supersystem bases yields results equivalent to those of supermolecular TDDFT calculations for systems solvated by many solvent molecules, using the previously studiedsystem of methylene-cyclopropene (cid:1)(cid:1)(cid:1) (H 2 O) 17 as anexample. By means of this exact reference embedding framework, we are able to disentangle solvent effects introduced in terms of the embedding potential from those caused by solvent response couplings, both for the PbE variant and for sTDDFT with approximate non-additive kinetic energy functionals. Furthermore, we show that the use of a monomer basis introduces significant errors for the environmental response contribution. Employing avirtual-orbital localization strategy on top of PbE-sTDDFT, we can also directly assess the impact of inter-subsystem charge-transfer excitations on the entire solvent effect, which turn out to play a significant role for the environmental response. Finally, we analyze the response effects introduced by the individual solvent molecules and their interdependence, and show that a simple, pair-wise additive correction for solvent response yields excellent results in the present example.

usually obtained by solving a truncated eigenvalue problem, starting from the results of monomer TDDFT calculations within the related framework of (uncoupled) frozen-density-embedding TDDFT (FDEu-TDDFT). This neglects all couplings with the remaining monomer excitations; in particular, environmental response effects (through "off-resonance excitonic couplings") are missing in this treatment. As recently shown in our lab, all three problems can be solved by switching to a projection-based version of sTDDFT in combination with a supermolecular basis set if all orbital transitions are considered. [12,13] While this approach is rather inefficient if used without further approximations, it can provide a reference point for benchmarking more approximate versions of subsystem TDDFT.
Previous studies have indicated that environmental response effects on excitation energies can, in principle, be captured by approximate sTDDFT versions. [3,14] But this requires explicit inclusion of a very large number of environmental transitions in the coupling procedure. The study in ref. [14] in fact explored certain strategies to reduce the computational effort. However, no full equivalence with a supermolecular approach could be demonstrated with this approximate scheme. Furthermore, the implementation used for sTDDFT calculations in refs. [3,14] was intended for exciton-coupling-like problems and hence always involved a full diagonalization of the (truncated) response matrix, assuming that this step is not the bottleneck of the calculation. But if a large number of environmental transitions are included, this type of solution becomes prohibitively expensive. A different route toward the determination of environmental effects in the context of approximate sTDDFT has been developed by the Pavanello group. In their real-time sTDDFT [15] framework, upon employment of periodic boundary conditions, many-body (coupling-) effects on excitation energies have been investigated both for liquid water [16] and molecules at surfaces. [17] A large variety of other quantum chemical analyses and practical solution strategies have been presented to model environmental response effects with different degrees of sophistication (besides classical models, which will not be considered explicitly here). For instance, Schwabe has formulated an approach starting from long-range perturbation theory, in which nonresonant excitonic coupling effects are related to dynamic polarizabilities of the environment, [18] similar to earlier developments in refs. [19,20]. Megow et al include transition-density couplings in terms of dispersive site-energy shifts for aggregates of dye molecules, [21][22][23] following the work by Renger et al. [24] Within the fragment molecular orbital method, an incremental ansatz using pair corrections to monomer excitation energies has been suggested. [25] In the context of wavefunction-in-DFT (WF-in-DFT) embedding, state-specific embedding potentials have been explored to mimic the environmental response. [26,27] Ricardi et al have shown that at least part of the electronic polarization of the environment in such calculations can implicitly be modeled by allowing the active-system density to delocalize into the environment through the use of supermolecular bases. [28] As an alternative, Goodpaster and co-workers have employed supermolecular TDDFT to approximately include environmental response into an embedded WF calculation [29] in an ONIOM-like fashion. Finally, we would like to mention that combinations of density-based embedding strategies with classical, polarizable embedding have been suggested to cover environmental response. [30,31] Since environmental response effects have, in the past years, seen a revived interest (see, for instance, the examples for chromophore-protein interactions in refs. [32,33]), we will, in the following, analyze the characteristics of environmental response effects and possible simplifications within sTDDFT, starting from an exact embedding strategy. We chose microsolvated methylene-cyclopropene (MCP) with up to 17 water molecules in the solvation shell as a test system, since (a) our previous work on state-specific embedding potentials has shown that large environment effects occur for MCP-water complexes, [27] (b) the system contains sufficiently many water molecules to explore approximations for efficiency increase, and (c) it is still small enough to easily facilitate full supermolecular reference calculations.

| THEORETICAL BACKGROUND
Subsystem TDDFT calculations following the formulation of ref. [2] proceed in two steps. In the first step (FDEu-TDDFT, see above), the local (uncoupled) excitations of the individual subsystems are calculated. Usually, only a few of the lowest transitions are computed in this step. Technically, this proceeds in a manner similar to a TDDFT calculation carried out for the individual subsystems, but with two important differences: (a) the orbitals and orbital energies used for the TDDFT step are obtained in subsystem DFT calculations, that is, under the influence of an embedding potential, rather than Kohn-Sham-DFT calculations on the isolated subsystems, and (b) the effective response kernel in that case involves terms arising from the non-additive exchange-correlation and kinetic energy contributions to the sDFT energy functional. [34] In the second step of the sTDDFT procedure (also called coupled FDE-TDDFT or FDEc-TDDFT), the response problem is solved in the truncated space of the local excitations obtained in the first step. This method was originally intended for systems of interacting chromophores, in which the splitting effect due to excited-state interactions of local excitations is often of great interest (see, for example, refs. [35,36]). In a minimal model, such a splitting can be described in terms of just one excited state per chromophore. Usually, a few more states per monomer may be needed to quantitatively converge the splitting for aggregates of chromophores. But converging this splitting for systems in responsive environments may require to explicitly couple the excitonically interacting states with a large number of environmental excitations. [14] Already in ref. [14], an attempt was made to reduce the number of coupled excitations required to describe an environmental screening effect on excitonically coupled chromophores based on their contributions to the polarizability of the system. However, the implementation was limited by the fact that always a full diagonalization of the response matrix was performed. Our new implementation, by contrast, allows us to perform Davidson-type subspace iteration solutions starting from a diagonalization of a truncated response problem on the basis of monomer excitations.
None of the previous tests was able to demonstrate equivalence of subsystem and supersystem results for chromophores in responsive environments. This requires accurate methods to handle (or circumvent) the NAKE contribution as well as the usage of a supersystem basis. In fact, we have recently demonstrated the equivalence of sTDDFT and supersystem TDDFT for small molecular dimers. [12] In this work, we extend that work by producing reference sTDDFT results for systems explicitly solvated by many solvent molecules. Subsequently, these results will be used as a starting point for an analysis of possible approximations.
In addition, we are decomposing the response within our exact sTDDFT calculations into local response (LR) and charge-transfer response (CTR). This is achieved via an additional virtual-orbital localization procedure (see below). It has to be kept in mind that projection-based embedding (PbE)-sDFT calculations with supermolecular basis lead to virtual-orbital spaces for each subsystem that contain the entire supersystem virtual space. The virtual-orbital localization then allows us to go back to a subsystem-based picture by separating virtual orbitals located on different subsystems, even though they formally appear for each system. As illustrated in Figure 1, this allows the separation of local and intersubsystem charge-transfer (CT) excitations in the FDEu calculations. In the coupling procedure, these selected orbital-transition spaces can then be used to analyze the individual environment response contributions.
Several virtual-orbital localization procedures were proposed in the literature (cf. refs. [37][38][39][40][41]). Here, we employ a simple strategy based on an occupied-orbital localization procedure proposed by Claudino et al [42] for the virtual-orbital space. Note that the following procedure is strictly applicable only for PbE in a supersystem basis. In this case, the virtual-orbital space for each subsystem is identical. [13] In a first step, the virtualorbital coefficient matrix of a particular subsystem A is orthogonalized: In the next step, the molecular orbital overlap matrix in the atomic orbital basis of subsystem A is constructed, Via a singular value decomposition of S A the rotated (localized) virtual orbitals are obtained asC Individual orbitals are then selected based on their singular values Σ A . For singular values >0.5, the orbital is assigned to be located on the particular subsystem, while a value <0.5 refers to an orbital located on the environment. Because of the fact that the intra-subsystem Fock matrix is not diagonalized by the rotated orbitals, an additional diagonalization is necessary. Note that this procedure is similar to the orthogonality-constrained basis-set expansion used in the embedding context to keep orbitals of different subsystems orthogonal. [43] Transforming the embedded Fock matrix of subsystem A (F A(Env) ) in the space spanned by the selected virtual orbitals yields a new set of eigenvectors  A and eigenvalues via diagonalization of  A . In a last step, the localized quasi-canonical virtual orbitals of subsystem A are obtained from the transformation ofC virt with  A : Because of this new set of virtual orbitals, additional coupling-matrix contributions related to external orthogonality, that is, the mutual orthogonality of orbitals in different subsystems, need to be taken into account. [44] Therefore, additional terms corresponding to virtual-orbital off-diagonal Lagrange multipliers δ ij F ab must be considered (similar to the external-orthogonality contributions for occupied orbitals derived in ref. [13]).

| COMPUTATIONAL DETAILS
All calculations were performed with a development version of the quantum chemistry program SERENITY, [45] using a def2-SVP basis set [46] and the B3LYP [47,48] exchange-correlation functional. For embedding calculations, either PbE (not to be confused with the exchange-correlation functional PBE [49,50] ) or an approximate NAKE functional was employed. The notation used throughout this study is as follows: XCFUNC[/NAXC]/ NAKE/basis-set(monomer/supermolecular). These are, respectively, the intra-subsystem exchange-correlation functional XCFUNC, the nonadditive exchange-correlation functional (NAXC, only specified if different from XCFUNC), the NAKE functional, the type of basis set, and a label to distinguish between monomer and supersystem basis. In case of PbE, the acronym "PbE" substitutes the NAKE functional and refers to bottom-up PbE (see below), unless explicitly stated otherwise. Furthermore, PbE calculations make use of the Huzinaga projection operator, [51,52] while those with an approximate NAKE functional utilize the LLP91 [53] functional.
The choice of the non-additive exchange-correlation functional raises a problem: For consistency reasons, B3LYP should be used throughout for exchange-correlation contributions. However, this is straightforward only in the context of PbE with externally orthogonal orbitals, which is why we chose the BLYP [54,55] functional as a substitute if external orthogonality is not ensured. Even though a monomer basis is generally preferred for efficiency reasons, PbE with external orthogonality can only be guaranteed in supersystem basis. Hence, we employed BLYP as a non-additive exchange-correlation functional for PbE in monomer basis here as well. An alternative would be to ignore the external-orthogonality violation in the PbE (monomer basis) calculations and use B3LYP as a non-additive exchange-correlation functional simply assuming perfect external orthogonality. Test calculations have indicated that this leads to very similar excitation energies with deviations smaller than 0.011 eV w.r.t. the B3LYP/BLYP/PbE (monomer basis) calculations. Note that in PbE calculations with monomer bases, kernel contributions possibly arising from small non-orthogonality effects were neglected.
In PbE calculations, artificially shifted occupied environment orbitals [56] are removed from the orbital-transition space in sTDDFT calculations in supersystem basis. Throughout this study, full coupling procedures were carried out (involving all transitions in all subsystems). These were started from a set of uncoupled guess vectors and utilized a Davidson-like root-following algorithm. [57] In top-down [56] calculations, the intrinsic bond orbital (IBO) [58] localization procedure was employed to localize the occupied orbital space. In case of bottom-up embedding calculations, four freeze-and-thaw (FaT) cycles were performed in order to relax the subsystem electron densities. The number of cycles was chosen after a careful convergence test of excitation energies within the different embedding schemes such that the lowest uncoupled excitation energies changed by less than 0.0001 eV with further FaT cycles.
The structure of the MCP surrounded with up to 17 water molecules was taken from ref. [27] and is displayed in Figure 2. Substructures MCPÁ Á Á(H 2 O) n are defined as those clusters including MCP and water molecules W01 till W n in this study.

| RESULTS
The presentation of the results is structured as follows: First, we will validate our coupled PbE approach for excited states utilizing MCP surrounded by several water molecules, and address the issue of the ambiguity in the underlying uncoupled excitations. Thereafter, uncoupled excitation energies from PbE and an approximate NAKE functional will be analyzed, followed by comparisons of coupled excitation energies from different approximate embedding schemes. Lastly, we present an incremental strategy for efficient description of environmental response. An overview of the response contributions covered by the methods used in this study can be found in Table 1.

| Validation
In our previous work, we have shown that PbE-sTDDFT is able to reproduce supermolecular results if a supersystem basis is employed and all excitations of all subsystems are coupled. [12,13]  which is of π ! π* type, are shown in Figure 3. HOMO and LUMO will be abbreviated as H and L, respectively, in the following. In the top-down calculations, the water environment is treated as one subsystem, while in the bottom-up calculations each water molecule is considered as a separate subsystem. Both types of calculations, though technically very different, reproduce the supermolecular results to very high accuracy. This confirms that the equivalence of subsystem and conventional TDDFT can be established in general cases with more than two subsystems and that sTDDFT is able to treat environmental response effects in precisely the same way as supermolecular TDDFT. However, especially the bottom-up calculations in supersystem basis become costly quite quickly with increasing system size. Hence the question remains as to how more efficient yet also more approximate versions of sTDDFT perform in terms of accuracy compared to this reference approach.
Before we address this question in more detail, we would like to highlight another important aspect of sTDDFT calculations. The abovementioned equivalence to supermolecular results holds only for the FDEc-TDDFT excitations. The intermediately obtained uncoupled (FDEu) excitations are, by contrast, not uniquely defined, [60] which is related to the nonunique density partitioning in sDFT. [61] This is illustrated for 17 cluster used throughout this study, taken from ref. [27]. Graphics: VMD (ref. [59]) Note: Charge-transfer response is denoted analogously as CTR a , CTR e , and CTR e, e , respectively. Entries in parentheses indicate contributions that are formally included, but suffer from the incorrect asymptotics of the embedding potential (Section 4.3). A detailed description of the FDEcc method is presented in Section 4.4.
MCPÁ Á ÁH 2 O by the data in Tables 2 and 3, which contain coupled and uncoupled excitation energies, respectively, from different PbE-TDDFT calculations. In top-down calculations, the density partitioning relies on an explicit orbital localization step, starting from a supermolecular calculation.
The localized orbitals are then assigned to one of the subsystems. The bottom-up scheme typically starts from isolated system orbitals, which are intrinsically subsystem-localized, and iteratively modifies them under the influence of the embedding potential. It can thus be regarded as yet another way of obtaining subsystem-localized orbitals that correctly reproduce the electron density of the supermolecule. The fully coupled results in Table 2 show that top-down PbE with different localization strategies [IBO and Foster-Boys [62,63] ] as well as bottom-up PbE yield virtually the same excitation energies for the excitations studied here, with deviations of 0.0001 eV or less. Also, the deviations between the different types of PbE-TDDFT and the supermolecular results are very small (0.0007 eV or less); the remaining discrepancies are mainly due to numerical aspects like the chosen grid accuracy or the predefined number of FaT cycles in the bottom-up procedure.
In case of uncoupled top-down or bottom-up results, the situation changes: The bottom-up PbE data (Table 3)  can even reach more than 0.1 eV for transitions that show a total solvent effect of the same order of magnitude.
As mentioned above, the reason for these discrepancies lies in the so-called nonuniqueness of the density partitioning in sDFT in general [61] (see, however, the work by Wasserman and co-workers in the related field of partition (TD)DFT [64,65]  different combinations of subsystem densities can sum up to the same (correct) supersystem densities. In the present case, the different PbE procedures lead to rather pronounced differences in the HOMO of MCP (Figure 4): In the case of top-down PbE with the IBO and Foster-Boys localization procedure, we observe a rather sharp polarization feature in the π system of MCP, directed toward the closer of the two hydrogen atoms of water. Contrary to that, the bottom-up calculation does not show such a feature, but rather has a small p-type contribution at the water oxygen, which resembles the shape of the corresponding feature in the HOMO of the supermolecular calculation. The HOMO from bottom-up PbE is thus much closer to the canonical HOMO of the supermolecular reference calculation than to those from the top-down calculations. This also results in different orbital energies, which then enter the TDDFT calculation in terms of orbital energy differences. In fact, the differences in the HOMO energies (Table 4) explain the major part of the differences in the excitation energies. Again, bottom-up PbE shows the smallest deviations from the reference for the HOMO energy, which in turn explains why bottom-up PbE yields uncoupled excitation energies in good agreement with the reference.
From these results we conclude that the bottom-up approach leads to orbitals that are in better agreement with the supermolecular orbitals, while the top-down procedure with its initial localization leads to larger deviations that translate into more drastic changes in the corresponding orbital energies and uncoupled excitation energies.
In the following, we will systematically introduce approximations in terms of basis, NAKE functionals, and consideration of couplings, as a function of the number of water molecules for substructures of the MCPÁ Á Á(H 2 O) 17 system.

| Uncoupled excitation energies
From here on, we will compare uncoupled excitation energies for the H MCP ! L MCP transition from different FDEu-TDDFT embedding procedures to the supermolecular reference. The corresponding data are shown in Figure 5. There, the uncoupled excitation energies from both sTDDFT with approximate NAKE functionals and from PbE-sTDDFT in supersystem and monomer basis are shown.
The reference excitation energies initially increase linearly until the fourth water molecule is added. Thereafter, the excitation energies show a less regular, partially oscillatory behavior upon addition of water molecules. Although all the four methods tested here are able to capture the general trend and change in excitation energy for an increasing water shell of up to 10 water molecules, the relative deviations differ substantially.
For PbE in monomer basis, uncoupled excitation energies are consistently too low and oscillate around 4.5-4.6 eV after the addition of W12. A similar trend can be observed for the uncoupled excitation energies obtained with an approximate NAKE functional in monomer basis. However, the excitation energies are consistently shifted to higher energies and level off at around 4.7 eV.
In supersystem basis, the picture changes. Here, the trend in the excitation energies for approximate NAKE functionals and PbE is similar.
Both systematically overestimate the reference excitation energy. In direct comparison, the two methods differ by a maximum value of 0.039 eV

| Coupled excitation energies
In the next step, the uncoupled excitations evaluated with the different embedding methods are fully coupled to the excitations of the water molecules in the environment. In general, the coupling of excitations located on different subsystems is intended to remedy deviations to the supermolecular reference due to yet neglected environmental response. In Figure 6, the coupled excitation energies for PbE and approximate NAKE functionals in both monomer and supersystem basis are shown. As already demonstrated (Figure 3), reference excitation energies are restored if we employ PbE in combination with a supersystem basis. For PbE in monomer basis, excitation energies systematically underestimate the reference, and the overall trend is almost identical to the excitation energies obtained in uncoupled calculations. Similar trends can also be observed for the coupled calculations obtained with approximate NAKE functionals, even though these, again, overestimate the reference.
The small influence of the environmental water molecules on the coupled excitation energies in case of embedding calculations in monomer basis indicates that important response contributions are missing. What is still missing in these calculations are the incoming and outgoing CT excitations, caused by the restriction of the basis set. Therefore, besides the intrinsic error introduced by the method itself, CT-related response effects seem to play an important role in order to restore reference excitation energies. In case of approximate NAKE functionals, CT response can, in principle, be captured. However, as already noted in ref. [66], the wrong distance-dependent behavior of approximate NAKE functionals can lead to artificially low-lying virtual orbitals on the environment. Therefore, a wrong description of CT response is rather likely in that case. We would like to note that even in PbE calculations, the CT-response effect might be overrated due to the fact that the exchange-correlation functional employed (B3LYP) only includes 20% of exact exchange and is therefore not able to correctly describe CT excitations. [67] More specifically, CT excitations are usually energetically underestimated with this functional and thus are closer in energy to the MCP valence excitations. This, in turn, could mean that their coupling effect is overestimated. Nevertheless, test calculations employing CAM-B3LYP [68] show a similar influence of the CT response on the excitation energy of the MCP, so that the effect seems to be present also with functionals better suited for CT excitations.
In order to separate the influence of local and CT response, we employ the orbital-localization procedure introduced in Section 2 to decom-

| Incremental strategies for environmental response
The full coupling procedure for a large number of subsystems, especially those of large size, is still very expensive. For this reason, we apply a pairwise coupling procedure in which our system of interest-MCP in the present case-is coupled only to a single water molecule at a time, while all other water molecules are treated as uncoupled environment. Subsequently, the effects of all pairwise couplings are added up to give the total coupling contribution, which is then added to the uncoupled excitation energy. This reduces the number of coupled subsystems per calculation step to exactly 2, and also the number of coupled states decreases significantly. Moreover, this approach is trivially parallelizable by distributing calculations for the different coupled water molecules. In this way, we can separate the influence of the individual water molecules on the supersystem excitation energy. We will refer to this approximation as cumulatively coupled frozen density embedding (FDEcc).
Individual contributions to the total FDEcc shifts together with the total FDEcc response contribution as well as the supermolecular reference for the PbE approach are shown in Figure 10.  9 ], the coupling effect of others is almost negligible. However, these water molecules can still have a strong indirect influence on water molecules closer to the MCP due to ground-state polarization.
We would like to address the role of the closest water molecule (W01) in more detail. It can be seen from Figure 10 that W01 contributes about −0.03 eV to the total environmental response, and this value is fairly independent of the total number n of water molecules considered (between −0.038 eV for n = 3 and −0.026 eV for n = 12). The only exception is the structure with n = 4, for which a contribution of −0.078 eV is observed. This can be rationalized as follows: In the full structure with 17 water molecules, W01 exhibits two hydrogen bonds to other water molecules, namely W04 and W05. Neither W04 nor W05 have a significant direct coupling effect. But while the addition of W04 changes the electronic structure of W01 in such a way that the coupling effect of the latter increases, the effect of W05 brings the coupling back to the expected range. Among the uncoupled transitions of W01, the H W01 ! L W01 excitation is the transition that couples most strongly to the H MCP ! L MCP excitation considered here. Figure 11 shows the H W01 ! L W01 excitation energy for all solvated substructures. The H W01 ! L W01 excitation In a next step, the FDEcc approximation was applied for coupling calculations with approximate NAKE functionals in supersystem basis.
The corresponding graph is shown in Figure 12. First of all, the magnitude of the contributions of each water molecule to the cumulative coupling is much smaller than in the case of PbE ( Figure 10). Furthermore, an almost constant shift per molecule is observed for the individual water subsystems. Especially, the prominent change in the coupling magnitude for W01 with the addition of more water molecules in the solvation shell cannot be observed. But also with approximate NAKE functionals we observe a close agreement of cumulatively coupled and fully coupled FDE excitation energies. Figure 13 summarizes these findings for approximate NAKE functionals and PbE and shows the resulting excitation energies. In the PbE case, the differences between monomer-basis and supersystem-basis calculations are even larger, and are significant already for the smallest clusters. This is not very surprising, as a truncation of the basis set here necessarily affects the projection operator, leading to increasingly strong violations of external orthogonality. Instead, for n = 17 and employment of an approximate NAKE functional in monomer basis and F I G U R E 1 2 Environmental response contributions ΔE exc to the excitation energy (in eV) of the lowest methylene-cyclopropene (MCP) excitation in MCPÁ Á Á(H 2 O) n substructures. Shown are contributions of the individual water molecules (in color code; see Figure 7) to the total FDEcc response contribution [dashed line; B3LYP/ BLYP/LLP91/def2-SVP (supersystem basis)]. For comparison, also both the fully coupled FDEc (gray linepoints) and the supermolecular reference (solid black line) are shown Supermolecular F I G U R E 1 3 Fully coupled and cumulatively coupled excitation energies E exc (in eV) for the LLP91 NAKE functional (blue circles and orange squares, respectively) and for bottom-up PbE (red triangles and green diamonds, respectively). Calculations were performed in supersystem basis supersystem basis, the errors are only 0.10 eV and 0.05 eV, respectively. We therefore want to emphasize the usefulness of calculations with approximate NAKE functionals using monomer bases for weakly interacting solvent systems if the goal of the calculation is a good semiquantitative description of environmental effects with high efficiency rather than an in-depth analysis of different contributions to the solvent shifts. This method reproduces the total (micro-)solvation effect in this system (0.88 eV) to within 11% at low computational cost. But, of course, it is not suited at all to study environment response effects.
Subsystem TDDFT also easily allows to selectively switch the response couplings for specific molecules on or off, which enables a detailed analysis of the contributions of individual solvent molecules to the overall solvent response shift for a given structure. This is particularly simple in the context of the cumulative evaluation of pairwise couplings between MCP and the different water molecules. In this framework, the total response of the environment is approximated as the sum of contributions of individual solvent molecules, neglecting the inter-solvent response couplings. We have shown that this approximation leads to very good agreement for the total environmental response compared to the reference calculations.
One conclusion we want to emphasize is the importance of a proper description of CT response in order to capture environmental effects onto the active system, which we inferred from our virtual-orbital localization procedure. At least for the example of MCP in water studied here, these effects are crucial for a quantitative modeling of environment response effects. However, most of the fragment-based approaches, which try to include the environment contributions through point charges, induced dipoles, or similar strategies, do not include CT response, which can lead to severe discrepancies between the actual and the modeled environmental contributions. This was also already recognized in the context of inter-chromophore CT response by Mennucci and co-workers. [69][70][71] A further analysis of this aspect is part of our ongoing work. In particular, it needs to be studied to which extent this observation is related to the typical error of approximate TDDFT for CT excitations.