Toward a Unified Modeling of Environment and Bridge-Mediated Contributions to Electronic Energy Transfer: A Fully Polarizable QM/MM/PCM Approach.

Recent studies have unveiled the similar nature of solvent (screening) effects and bridge-mediated contributions to electronic energy transfer, both related to the bridge/solvent polarizability properties. Here, we exploit the similarity of such contributions to develop a fully polarizable mixed QM/discrete/continuum model aimed at studying electronic energy transfer processes in supramolecular systems. In the model, the definition of the three regions is completely flexible and allows us to explore the possibility to describe bridge-mediated contributions by using a polarizable MM description of the linker. In addition, we show that the classical MMPol description of the bridge can be complemented either with an analogous atomistic or with a continuum description of the solvent. Advantages and drawbacks of the model are finally presented and discussed with respect to the system under study.


Introduction
Electronic energy transfer (EET) consists in the radiationless migration of the excitation energy of a sensitized donor (D) to a proximate acceptor (A). In photosynthesis, ultrafast EET among antenna pigment-protein complexes allows to funnel the captured sunlight to reaction centers with very high quantum efficiencies. [1][2][3][4] The same principle is used in organic photovoltaic cells, where the absorbed excitons have to be efficiently transferred to charge-separation interfaces, 5,6 whereas in organic light-emitting diodes (OLED) the processes is reversed, and EET is used to direct the excitons generated upon charge injection. 7,8 The basic design principles of EET processes can be understood based on Förster theory,9 which relates the transfer rate between weakly coupled D/A molecules with the square of the electronic coupling matrix element between the initial and final states of the EET reaction, V , and a spectral overlap factor between normalized donor emission and acceptor absorption lineshapes, In intermolecular EET, the effect of the molecular environment on the electronic coupling V is typically described in terms of a screening effect that significantly slows down the EET process. In Förster theory, this effect is described by a rather crude screening factor s = 1/n 2 , where n is the refractive index of the medium, so that the total coupling is cast in terms of a direct D/A interaction V s and a screening factor s, V = s ·V s . In intramolecular EET processes, however, the D/A units are covalently linked by some molecular spacer that can significantly modulate the coupling between the units, [10][11][12] in so-called through-bond contributions either mediated by superexchange, as in electron transfer, 13,14 or by modifying the D/A Coulombic interaction. 12 In the last decade there have been important efforts in order to accurately model both screening effects [15][16][17][18][19][20][21][22][23] and bridge-mediated effects to EET. [24][25][26][27][28][29][30][31] Screening effects, which arise from the coupling of solvent electronic states with the D/A states, 20 can be modeled either using continuum dielectric models, mixed quantum mechanics/molecular mechanics (QM/MM) methods, or by including them in the QM calculation as in subsystem TDDFT. 21,32 In continuum models, the solvent contribution is dictated by its macroscopic optical dielectric constant, i.e. approximately the refractive index related to the solvent polarizability, 33 whereas in mixed QM/MM models solvent molecules are explicitly considered and described through a classical polarizable force field. 19 Bridge-mediated effects, on the other hand, arising from coupling of the bridge states to the D/A states, are typically modeled by including the linker in the QM calculation. [24][25][26][27][28][29][30][31] Recent studies, however, have pointed out that the leading contributions to bridge-mediated singlet-singlet EET often arise from Coulomb contributions rather than superexchange, 27,28,30,31 thus relating the enhancement or screening of the D/A coupling with the bridge polarizability. 28,30 In the present study, we develop a fully coupled QM/discrete/continuum method for energy transfer that combines a continuum dielectric description of the environment with a flexible definition of QM and polarizable MM regions. This three-level method allows us to explore the possibility to describe bridge-mediated contributions by using a polarizable MM description of the linker combined with a continuum model for solvent effects. In addition, we show that the classical MMPol description of the bridge can be complemented either with an atomistic or a continuum description of the solvent, and discuss some of the advantages and drawbacks of these two options depending on the chemical system under study. Applications to different D-bridge-A systems are reported to determine potentials and limitations of the method.

Theory
We implemented a fully polarizable mixed QM/discrete/continuum model, where a solute, described at QM level, is surrounded by a set of polarizable MM atoms and by a polarizable structureless medium. The latter is treated with the Polarizable Continuum Model (PCM), 34 implemented according to the Integral Equation Formalism (IEF-PCM 35,36 ), while the MM atoms are represented as a set of charges and induced dipoles, using the MMPol model recently developed by some of the present authors. 19 This MM layer may be employed to describe both solvent molecules and non-QM parts of the solute itself. Here below we first report the main elements of the general implementation (which is similar to the one presented by Steindal et al. 37 ), and successively we detail on its extension to describe EET processes within a TD-DFT formalism.

The electrostatic problem
In the PCM model, the polarization of the solvent is represented by placing a set of "apparent" charges, {q PCM i }, on the mesh of the molecular cavity's surface. These charges are obtained from the electrostatic potential of the QM solute, plus eventually the potential generated by the MM distribution, using the equation: where Q Q Q is the vector containing the {q PCM i } set of PCM charges on each of the N tes surface finite elements (generally called "tesserae"), V V V a vector with the value of the potential on the same tesserae, and T T T and R R R are matrices depending on geometrical factors and on the solvent permittivity ε (either the static or the optical one, depending on the case).
On the other hand, the MM molecules are described as a set of fixed charges, {q MM i }, centered on the atoms, and a set of induced dipoles, {µ µ µ i }, whose positions are generally referred to as polarizable sites. The induced dipoles can be obtained, similarly to the PCM case, from the electrostatic field due to the QM solute and the PCM and MM distributions, calculated on the polarizable sites: where M M M is the vector containing the induced dipoles {µ µ µ i } on the N pol polarizable sites, E E E the value of the electric field on the same sites and α α α the MMPol matrix containing the polarizability tensors {α α α i } and whose expression depends on the MMPol model used.
Eq. Eq. (2) and Eq. (3) differ from the standard PCM and MMPol equations in that the potential and electric field also contain "mixed" contributions from both the MMPol and PCM distributions, and therefore the two problems cannot be separated.
In more detail, the elements of the total potential of Eq. (2) can be written as: where indices t, c and p label tesserae, MM charge sites and polarizable sites, respectively.
Similarly, the electric field on each polarizable site, appearing in Eq. (3), can be expressed as: where 3(µ µ µ p ·r r r p,p )r r r p,p − µ µ µ p r r r p − r r r p 3 = a a a p + ∑ p =p P P P p,p µ µ µ p and Substituting the total potential and electric field into Eq. (2) and Eq. (3), and using a more compact matrix form, one can write: where the Eq. (8) must be solved simultaneously. Noting that a a a and b b b are respectively the electric field and potential due to the MM charges, and that S S S = −U U U † , the system of equations can be rearranged as: The induced dipoles and PCM charges are therefore obtained simultaneously by solving, either iteratively or by matrix inversion, Eq. (9). The MM-PCM matrix is dimensioned (N tes + 3N pol ) × (N tes + 3N pol ); matrix U U U is responsible for the MMPol/PCM mixing, and setting it to zero would decouple the problem.

The QM problem
Within the Born-Oppenheimer approximation, the QM solute satisfies the Schrödinger equation, where the Hamiltonian is divided into an unperturbed (gas-phase) termĤ 0 and an perturbation termĤ env accounting for the coupling with the environment.
The perturbation term includes all the energetic terms arising from the presence of a set of PCM where:Ĥ In the last equations,V andÊ E E are the potential and electric field operators, respectively, due where the label 'QM+c' indicates, as in Eq. (9), that the interaction is with the QM system (nuclei and electron density) plus the MM charges.
The Fock matrix contains extra 1-electron and 2-electron terms due to the environment (plus the constant charge-charge interaction), so that it reads: where h h h 0 and G G G 0 are the usual gas-phase 1-and 2-electron terms (plus the eventual contribution from the MM charges), while h h h env and X X X env arise from the interaction of the QM electron density with the PCM charges and MM charges and induced dipoles. In particular, the 1-electron term h h h env accounts for the interaction with both the PCM and the MM charges and the dipoles induced by the fixed charges (nuclei and MM charges), while X X X env accounts for the interaction with those induced by the electron distribution, i.e. those found by solving Eq. (9) with the electric field and potential due to the QM electronic wavefunction. These PCM charges and induced dipoles therefore depend themselves on the wavefunction, and this is stressed by explicitly indicating the term dependence on the density matrix P P P.

TD-DFT and Electronic Energy Transfer
The QM/MMPol/PCM method has been implemented within the TD-DFT Linear Response scheme, in the same way it had been done for the PCM and MMPol independently. In a basis of single excitations between KS orbitals, solving Eq. (18) provides the excitation frequencies ω n and the corresponding eigenstates X X X n Y Y Y n † : where A A A and B B B form the Hessian of the electronic energy. The inclusion of both PCM and MMPol effects in this picture is simply a combination of the two effects, 19,33,38

so that matrices A A A and B B B
can be written as: where indices i and j label occupied orbitals and a and b virtual orbitals, and where matrices K K K, C C C and D D D are written in terms of the orbitals as follows: K ai,b j = dr r rdr r r ψ * a (r r r)ψ i (r r r) 1 |r r r − r r r | + g xc (r r r, r r r ) ψ * b (r r r )ψ j (r r r ); The term in Eq.
where the indices D and A refer to the donor and acceptor, and S S S includes the overlap between donor and acceptor orbitals. When the D/A interaction is turned on, the energies are no longer degenerate, and the electronic coupling between donor and acceptor is defined as half the resulting energy splitting. If the interaction is treated as a perturbation, Eq. (24) can be solved approximately and the coupling can be obtained from the first-order perturbed energies as: 15 Note that matrices K K K DA , C C C PCM DA and K K K MMPol DA differ from those of Eq. (21) -Eq. (23) in that the orbitals labelled i and a belong to the donor, while those labelled j and b belong to the acceptor; the same holds for the overlap matrix S S S DA . When these matrices are contracted with W W W † D and W W W A , the resulting contributions to the first-order coupling can be written in terms of the diagonal elements of the one-particle transition density matrices ("transition densities"),ρ D andρ A : V Coul = dr r rdr r r ρ † A (r r r) 1 |r r r − r r r |ρ D (r r r ); V xc = dr r rdr r r ρ † A (r r r)g xc (r r r, r r r )ρ D (r r r ); The definitions reported in Eq. (26) show that the first-order electronic coupling between two chromophores can be calculated from the transition densities of the non-interacting chromophores.
In practice, a coupling calculation is performed in two stages: first the transition densities of the

Computational details
All the QM calculations were run at the TD-DFT level, employing the CAM-B3LYP functional 39 and the 6-31G(d) basis set, using a locally modified version of the Gaussian09 suite of codes. 40 The choice of the functional and the basis set was done mainly to maintain the consistency with previous studies on the same systems.
As for MMPol parameters, the fixed MM charges were obtained from a fit of the electrostatic potential of the molecule or fragment, using the Merz and Kollman method. 41 The polarizable sites coincide with the MM atoms. The values of the isotropic polarizabilities placed on each atom depend on the MM treatment of polarization chosen. We adopted the Thole model, which avoids intramolecular overpolarization problems by using a smeared dipole-dipole interaction tensor. 44 Atomic isotropic polarizability values were taken from the fit of experimental molecular polarizabilities performed by van Duijnen and Swart using the linear version of Thole dipole-dipole tensor. 45 We also tested the polarizabilities derived recently in the context of the Amber force field, 46 always finding little dependence of the results on the model used.
The continuum solvent description was based on the IEF formulation of the PCM model, 35,36 using the discretization procedure into point charges available in Gaussian09 asking for the Gas-sian03 defaults.
In order to obtain representative structures of the distribution of toluene solvent molecules around the PDI-TDI system, we also performed classical molecular dynamics simulations of the dyad, in which the structure of the PDI-TDI molecule was kept frozen. The simulation system was built by solvating the dyad in a toluene box (buffer zone of 20 Å) using the Leap module of the Amber9 suite of programs. 47 Both the toluene solvent and the PDI-TDI dyad were described using the GAFF force field. 48 The system was initially thermalized by running a 100-ps MD simulation at constant volume in order to increase the temperature from 100 K to 298 K. Then, a 100-ps equilibration at constant pressure (1 atm) and temperature (298 K) was performed using standard coupling schemes in order to reach an appropiate density of the toluene solution. Finally, the simulation was extended for 2 ns for production purposes. All runs were performed with Amber9 using an integration time step of 1 fs, periodic boundary conditions, the Particle Mesh Ewald approach to deal with long-range electrostatics, and a nonbonded cutoff equal to 10 Å. For pressure regulation, an isothermal compressibility of 92 TPa −1 was used for toluene.

Results
There are various instances when a coupled fully-polarizable QM/MM/PCM model is desirable.
We present and discuss here the study of two different cases. Fückel et al. 28 and Chen et al. 30 ), and that the bridge effect is in general a complex phenomenon dependent on the bridge orientation, and possibly characterized by anisotropic polarizabilities. 12 For this reason, if one wants to carry out a mixed study on DBA systems, where the chromophores are described at a QM level and the bridge at a lower level, a model retaining structural information of the bridge and accounting for its anisotropic polarizability should be employed, rather than averaged models. At the same time, full-QM models, where both the chromophores and the bridge are treated quantum-mechanically, may be cumbersome and computationally expensive. The polarizable MMPol model 19 seems to be an ideal compromise, as it greatly reduces the computational requirements of a full QM description, but at the same time correctly describes the bridge charge distribution and polarizability. The first aim of this section is to investigate whether such a classical polarizable description of the bridge correctly reproduces the enhancement of the D/A coupling which has been observed in some DBA systems, and which has been ascribed to mediated Coulomb interactions. At the same time, we expect our model to fail whenever the bridge effect is quantum-mechanical in character (e.g. in case of a superexchange effect).
Two families of DBA systems, characterized by a significant effect of the bridge on the electronic coupling, have been studied. The systems have been selected so to allow a direct comparison with the "exact" benchmark, e.g. a full QM description. The general approach followed here consists in comparing the coupling results obtained when the bridge is treated quantum-mechanically to those obtained when it is treated classically, with the MMPol model. The solvent is then considered as a PCM charge distribution, placed on the molecular cavity.

PDI-TDI and PMI-TDI systems
Previous studies on electronic energy transfer in a perylene diimide -terrylene diimide (PDI-TDI) dyad, separated by a substituted terphenyl spacer (see Figure 1) The couplings obtained are reported in Table 1 (exchange-correlation and overlap terms are not reported as they are negligible). The total couplings are shown in Figure 2. Red bars refer to calculations in vacuo, blue ones to those in toluene.
Comparing the QM models, M0 through Mc, we notice the strong coupling enhancement due to the presence of the bridge between the chromophores; such effect is estimated in Table 1  Assuming that the contributions to the coupling arise from both through-space polarization effects and through-bond QM effects, and considering that the latter are completely disregarded by the classical MMPol model, we can interpret the discrepancies observed between the two models and estimate the through-space polarization effect of the bridge to be approximately ranging from 68% to 60% of the total coupling (as obtained from results in vacuo and in toluene, respectively).    When the solvent is included in the picture, the resulting total coupling values (plotted as blue bars in Figure 2) show an evident decrease due to the overall screening effect of the solvent, which reduces the total coupling to up to 50% of the value in vacuo. In particular, see also the explicit PCM term, V PCM , in Table 1 As Figure 5 shows, the coupling enhancement effect induced by the pPh spacer is very large, but it is almost completely lost in MMPol calculations, as expected. The results do not qualitatively change when we include the solvent (toluene) but what we see is just a reduction of the coupling due to the explicit PCM term.

ZnFbB(CH 3 ) 4 system
Another system studied is a porphyrin dimer. Several studies showed that dimers of Zinc (Zn) and metal-free (Fb) porphyrins, linked by semirigid bridges, have interesting properties; in particular, some are characterized by a rapid and highly efficient energy transfer from the Zn-porphyrin (donor) to the metal-free one (acceptor). [50][51][52] We present here the results on one of these dimers, having a diaryl-ethyne linker, referred to as ZnFbB(CH 3 ) 4 in the paper by Strachan et al.; 52 its structure is shown in Figure 6.
In order to assess the influence of the bridge on the coupling, and whether a classical description is appropriate for it, we have built four models following what did before for the other two systems.
As before, model M0 completely neglects the bridge whereas Mc includes the entire bridge at QM level. We note that, in Mc, the partition into donor and acceptor can be done either associating the  Table 2 reports the excitation energies obtained in vacuo for the four models.  For what concerns the coupling results, we present in Table 3 the couplings between all possi-ble pairs of the bright states reported in Table 2 (the third and fourth excited states). The oscillator strengths of these B-states are between 1.2 and 1.6, while those of the lower Q-states are approximately between 0.001 and 0.01. Table 3: Electronic coupling terms in cm −1 , in vacuo, between the strong B-states of the ZnFb porphyrin dimer studied. V i, j indicates the electronic coupling between the donor's S i state and the acceptor's S j state.
59. For models Mx, the only important contribution to the coupling comes from the Coulomb interaction, since exchange-correlation and overlap contributions (not reported) are negligible. In order to compare the different models one has to take into account that shifts between the (quasi) degenerate states are possible moving from one model to the other. In particular, from the analysis of the data reported in Table 3 it is evident that introducing the full bridge induces a shift between the states 3 and 4 in the donor. As a result V 3,i of Mc corresponds to V 4,i of M0 and M1 and viceversa: this shift is exactly reproduced by the MMPol description.
To have a more direct analysis, in Figure 7 we report a graph of the total couplings obtained as a sum of all the couplings listed in Table 3. Both results in vacuo and in toluene are shown. We

Description of the solvent
Discrete and continuum models follow almost opposite approaches in describing the effect of the solvent on QM solutes, and they may be thought as complementary for what concerns their advantages and disadvantages. Discrete QM/MM models have the advantage of retaining information on the solvent structure, and therefore of being able to describe short-range, specific solute-solvent interactions with good accuracy. At the same time, in order to get a correct picture of the dynamic solute-solvent interaction, it is necessary to sample a large number of configurations from a molecular dynamics simulation. Conversely, continuum solvation models describe the solvent in an averaged way, and therefore do not need configuration sampling. However, the solvent structural information is lost, and specific solute-solvent interactions are neglected.
The idea beyond a mixed discrete/continuum modeling of the solvent is then to try and combine the two models in order to exploit their complementary strengths: a discrete representation of the solvent molecules at short range may correctly describe the specific interactions, even using a relatively small number of configurations. For what concerns the average, long-range interactions, a continuum model may be employed beyond the discrete shell of solvation.
With this in mind, we have carried out a study to verify the advantages of such mixed approach in the calculation of the chromophore properties and electronic coupling of the PDI-TDI dyad in toluene. First, we ran a classical MD simulation of the system and extracted 11 structures from the trajectory. For each structure, several cutoff radii were applied, so that toluene molecules located further than the cutoff from the closest QM atom were discarded. Cutoff values of 0 (with no explicit solvent molecules), 3, 5, 8, 10, 15 and 30 Å were used. The average number of solvent molecules, n solv , included in each case is reported in Table 4.
For each configuration, at each cutoff value, two coupling calculations are run, both describing the PDI-TDI dyad at QM level (using the previous Mc model) and the explicit solvent molecules at MMPol level. In one set, however, an extra PCM layer is considered to account for the long-range solvent effect beyond the cutoff. We will refer to this set as QM/MMPol/PCM, while the set that does not include the continuum solvent description on top of the discrete one will be referred to as It is evident that each single QM/MMPol/PCM calculation is generally much more computationally expensive than each QM/MMPol one. However, the data reported in (Table 4) does not include any possible optimization of the PCM description when combined to the MMPol. As a matter of fact, it is reasonable to assume that when the first solvation shells are treated explicitly, the additional PCM description can be treated with lower numerical accuracy than the case where a pure PCM description is used. In particular, the quality of the mesh used can be significantly reduced with a corresponding important reduction of N tes . This and other numerical aspects of the implementation have not been investigated in this paper but surely need to be optimized to obtain the required efficiency: efforts in this direction are in progress. Finally, Figure 10(c) reports the solvent screening factor, s, which can be thought, in a simple dipole-dipole picture, as the effective screening of the Coulomb interaction due to the solvent. It is defined as: s = V TOT V Coul , so that V TOT = sV Coul . In Förster theory, s is calculated as the inverse of the solvent refractive index squared, but other approximations exist, based on different assumptions.
The Förster screening factor is indicated in Figure 10(c) as a horizontal black line, and is in good agreement with our pure-PCM result (horizontal red line), around 0.45. The limit of the QM/MM and QM/MM/PCM calculations is also close, at around 0.5. This confirms our previous conclusion, that the solvent effect is mostly an average, long-range effect, which is well reproduced by a pure PCM model. 22

Conclusions
We have presented the formulation, the computational implementation and some applications of a fully polarizable QM/MM/PCM approach to describe EET processes in (supra)molecular systems in condensed phase. The model allowed us to assess the possibility to describe solvent (screening) and bridge-mediated contributions to EET in solvated bridged bichromophoric systems based on a classical polarizable description of the solvent and bridge regions, as suggested by recent theoretical studies. In order to retain the structural information of the bridge and its anisotropic polarizability, the molecular linker is described through a classical polarizable MM force field.
The solvent, conversely, can be described either using a classical MM description or a continuum solvation model. Moreover, a mixed solvent description can also be adopted, where only the first solvation shell is described explicitly, in order to retain specific solute-solvent interactions.
We have shown that the model is able to describe the bridge-mediated enhancement of the EET coupling in different dyads (such as perylene linked by phenyl-based spacers and porphyrins with a diaryl-ethyne linkers). We have also illustrated the limitations of the model which clearly cannot be applied to describe the bridge effect in systems where the bridge is strongly coupled to the donor moiety via charge transfer interactions.
Moreover, we have shown that both MM and PCM models provide a good description of solvent screening effects, and this choice should be considered depending on the system under study.
While it is advantageous to model solvation through PCM in homogeneous media such as a standard solvent, an atomistic MM force field can improve the description of the environment in heterogeneous systems, such as organic crystals, polymers, or biological macromolecules. Indeed, in organic materials the separation between "bridge" and "environment" regions can be less straightforward, and thus a unified modeling of the whole system at the MM level allows to tackle both effects on an equal footing. Even in a homogeneous solvent, it can be interesting to include the first solvation shells at the MM level if strong specific solute-solvent interactions are expected to be important. An atomistic MM description can also be used to provide insights on the role of different structural regions of the system on the enhancement/screening of EET interactions, because the explicit environment-mediated contribution to the coupling can be dissected into atomic or group contributions, as we showed recently for a photosynthetic antenna complex. 22 On the other hand, a dielectric continuum description of the environment (solvent) can be advantageous whenever detailed structural information of the environment is missing, and prevents the need to perform molecular dynamics simulations of the system and subsequently a number of EET calculations in order to account for configurational sampling. Moreover, the easiness of continuum models to be extended to very different environments such as interfaces, membranes, as well as composite systems including metal nanoparticles, 53 makes this integration with a polarizable MM approach a very promising strategy to describe EET processes in systems of increasing complexity.