domain: chim phys

text:
Cristina E González-Espinoza

Christopher A Rumble
email: crumble@psu.edu

Daniel Borgis

Tomasz A Wesolowski
email: tomasz.wesolowski@unige.ch

Quantifying Fluctuations of Average Solvent Environments for Embedding Calculations

In the context of employing embedding methods to study spectroscopic properties, the viability and effectiveness of replacing an ensemble of calculations by a single calculation using an average description of the system of study are evaluated. This work aims to provide a baseline of the expected fluctuations in the average description of the system obtained in the two cases: from calculations of an ensemble of geometries, and from an average environment constructed with the same ensemble. To this end, the classical molecular dynamics simulation of a very simple system was used: a rigid molecule of acetone in a solution of rigid water. We perform a careful numerical analysis of the fluctuations of the electrostatic potential felt by the solute, as well as the fluctuations in the effect on its electronic density, measure through the dipole moment and the atomic charges derived from the corresponding potential. At the same time, we inspect the accuracy of the methods used to construct average environments. Finally, the proposed approach to generate the embedding potential from an average environment density is applied to estimate the solvatochromic shift of the first excitation of 1 acetone. In order to account for quantum-confinement effects that may be important in certain cases, the fluctuations on the shift due to the interaction with the solvent are evaluated using Frozen-Density Embedding Theory. Our results demonstrate that, for normally-behaved environments, the constructed average environment is a reasonably good representation of a discrete solvent environment.

Introduction

Dealing with the electronic properties of molecules requires the use of quantum mechanics. Given the ubiquity of the condensed phase in laboratory chemistry, the calculation of electronic properties in condensed matter is indispensable for connecting theoretical calculations with experimental observables. The large size of such condensed phase systems makes the electronic-structure calculations computationally challenging and approximations are essential for maximizing computational efficiency. A commonly encountered approach for reducing the cost of condensed phase quantum calculations are multi-level methods where the full system is partitioned into two pieces treated with varying levels of complexity. The first partition is a small subsystem of interest, here termed fragment A, which is treated with high-level quantum chemical methods. This subsystem is then embedded within a surrounding environment, here called fragment B, to which various approximations and lower-level methods are applied to minimize computational cost.

An important application of multi-level methods is to the computational study of liquidphase photochemistry, where a quantum mechanical chromophore is embedded in an approximation of the solvent environment. In such studies, one is typically interested in changes in solute electronic excitation energy brought about by the solvent environment, which is characterized by

δν ex = ν emb -ν iso , (1)
where ν iso is the excitation energy of the isolated solute in vacuum, ν emb is the excitation energy of the solute embedded in the solvent environment, and δν ex the 'solvent shift' of the excitation energy. The quality of the excitation energy calculation is usually evaluated by comparing ν emb to the maxima of measured electronic absorption spectra, whereas the shift is compared to the shift in the chromophore's absorption maximum between the solvent of interest and the maximum in a non-polar solvent such as n-hexane, where the non-polar solvent measurement acts as a stand-in for the more demanding gas-phase measurements. Single-reference quantum chemical methods, such as time-dependent Density Functional Theory (TDDFT), [START_REF] Hirata | Time-dependent density functional theory within the Tamm-Dancoff approximation[END_REF][START_REF] Dreuw | Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules[END_REF] configuration interaction (CI), [START_REF] Del Bene | Self-Consistent Molecular Orbital Methods. X. Molecular Orbital Studies of Excited States with Minimal and Extended Basis Sets[END_REF][START_REF] Foresman | Toward a systematic molecular orbital theory for excited states[END_REF] equation of motion coupled-cluster (EOM-CC), [START_REF] Stanton | The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties[END_REF] among others, are used to model the solute electronic transitions, and a number of different approximations can be employed to describe the solvent environment. One method for treating the solvent is to determine a single average solvent environment in which to embed the chromophore. This approach has the obvious benefit of requiring only two calculations for determining δν ex : one for ν iso and one for ν emb . In the polarizable contiuum model (PCM), [START_REF] Tomasi | Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent[END_REF][START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] the environment is described as a polarizable dielectric continuum. In the simplest implementation of PCM, the solvent is modeled as a linear isotropic continuum characterized by the static dielectric constant, , of the bulk solvent. [START_REF] Tomasi | Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent[END_REF] In more sophisticated formulations, the solute cavity is given a molecular shape and the reaction field is described in terms of polarization charges or reaction field factors included in the solute Hamiltonian. This approach allows for self-consistency to be reached between the solute wavefunction and the solvent polarization by means of iterative procedures. In particular, the conductor-like variant of PCM (C-PCM), [START_REF] Barone | Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model[END_REF] based on the COSMO methods, [START_REF] Klamt | COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient[END_REF] has been become very popular due to its ease-of-use, low computational cost, reasonable agreement with experiment, and widespread implementation in quantum chemistry software packages.

Instead of embedding the solute in a dielectric continuum, one can determine an effective solvent potential which acts on the solute. If this potential is constructed using only solvent electrostatics it is called the 'averaged solvent electrostatic potential' (ASEP). [START_REF] Sanchez | Study of solvent effects by means of averaged solvent electrostatic potentials obtained from molecular dynamics data[END_REF] Such a potential was constructed in Ref. 10 by placing point-charges at grid points to approximate the average solvent configuration obtained from molecular dynamics (MD) simulations. The values of the charges are fitted to reproduce the ASEP at the position of the solute by means of the ChELPG procedure, [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] while the position and the number of point-charges is optimized. A second method for constructing the ASEP utilizes scaled charges in all the positions of solvent atoms generated by an MD simulation, as was done by Coutinho et al. [START_REF] Coutinho | An efficient statistically converged average configuration for solvent effects[END_REF] These workers constructed an 'average single configuration' using 60 000 scaled charges from 100 superposed configurations of 200 water molecules centered around a single acetone solute. An important aspect of this procedure, termed as ASEC, is the selection of the configurations, which ensures convergent results by using statistically uncorrelated configurations that extend to outer layers of the solvent. The ASEC method has been used to model the effect of the solvent on spectroscopic properties, [START_REF] De Almeida | Solvent effects on the spectroscopic properties of Damascone derivatives: A sequential Monte Carlo/Quantum Mechanics study[END_REF] such as NMR shielding parameters, absorption spectra, [START_REF] Bistafa | Combining ab initio multiconfigurational and Free Energy Gradient methods to study the π → π * excited state structure and properties of uracil in water[END_REF]15 electronic hyperpolarizabilities, [START_REF] Torres | First hyperpolarizability of isomers of pyridinium N-phenoxide betaine dye in solution using the ASEC-FEG method[END_REF][START_REF] Brandão | Assessing the structure and first hyperpolarizability of Li@B10H14 in solution: a sequential QM/MM study using the ASEC-FEG method[END_REF] and more sophisticated protocols to optimize the free energy of the full system. [START_REF] Orozco-Gonzalez | An Average Solvent Electrostatic Configuration Protocol for QM/MM Free Energy Optimization: Implementation and Application to Rhodopsin Systems[END_REF] Another category of continuum models of the solvent is based on Frozen-Density Embedding Theory (FDET). [START_REF] Wesołowski | Embedding a multideterminantal wave function in an orbital-free environment[END_REF] In these methods, the embedding potential includes a non-linear term in the solute electron density, which is responsible for describing the electronic exchange-correlation and Pauli repulsion between the fragments. Previous studies done with FDETbased methods have used the average solvent density to describe solvatochromic shifts for aminocoumarin, coumarin 153 and acetone in several solvents. [START_REF] Kaminski | Modeling Solvatochromic Shifts Using the Orbital-Free Embedding Potential at Statistically Mechanically Averaged Solvent Density[END_REF][START_REF] Zhou | Multi-scale modelling of solvatochromic shifts from frozen-density embedding theory with non-uniform continuum model of the solvent: the coumarin 153 case[END_REF][START_REF] Zhou | Firstprinciples simulation of the absorption bands of fluorenone in zeolite L[END_REF][START_REF] Shedge | Nonuniform Continuum Model for Solvatochromism Based on Frozen-Density Embedding Theory[END_REF][START_REF] Laktionov | Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations[END_REF] Comparisons between these average environment calculations and experimental spectra are, on the surface, simple and straight-forward. Two single-point calculations are performed to determine δν (calc) ex , two absorption spectra measured to determine δν (exp) ex , and the resulting maxima and shifts compared. In reality, condensed phase spectra are not the result of a single chromophore in a static average environment, but of a fluctuating distribution of flexible chromophore conformations in an ensemble of solvent configurations, each contributing a distinct excitation energy to the observed spectrum. The resulting distribution of excitations, in conjunction with other factors such as vibronic transitions and finite excited state lifetimes, introduces substantial broadening (on the order of 0.34-0.40 eV for coumarin 153, a common solvatochromic probe, in dipolar liquids) [START_REF] Horng | Subpicosecond Measurements of Polar Solvation Dynamics: Coumarin 153 Revisited[END_REF] and asymmetry to the measured spectra. The broadening and asymmetry complicates determination of the experimental 0→0 solute transitions for comparison with theory, hence the use of the spectral maximum as a proxy. The single-point nature of average environment methods does not allow for straightforward determination of spectral broadening and makes it difficult to evaluate differences between δν (calc) ex and δν (exp) ex . One way for accessing these fluctuations is to use molecular dynamics (MD) simulations (or Monte Carlo sampling [START_REF] Acevedo | Quantum and Molecular Mechanical (QM/MM) Monte Carlo Techniques for Modeling Condensed-Phase Reactions[END_REF][START_REF] Kawashima | Monte Carlo microsolvation simulations for excited states using a mixed-Hamiltonian model with polarizable and vibrating waters: Applications to the blueshift of the H2CO excitation[END_REF] ) to generate ensembles of discrete solute/solvent configurations on which to perform an ensemble of single-point quantum mechanical (QM) calculations. [START_REF] Senn | MM Methods for Biomolecular Systems[END_REF][START_REF] Nochebuena | Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials[END_REF] Due to the large number of solvent nuclei needed for simulating realistic dynamics, the solvent cannot be treated explicitly, and approximate descriptions of the discrete solvent environment are needed to make the QM computations tractable. The most simplest method for describing a discrete solvent environment is point-charge embedding (PCE), where the solute is embedded the electric field generated by the surrounding solvent atoms. In order to capture the effects due to mutual polarization or other non-electrostatic interactions, the assignment of charges values for the environment and the embedding process can be made self-consistent. Examples of such polarizable embedding methods [START_REF] Olsen | Excited States in Solution through Polarizable Embedding[END_REF] are the ASEP-MD methods [START_REF] Sanchez Mendoza | A mean field approach that combines quantum mechanics and molecular dynamics simulation: The water molecule in liquid water[END_REF][START_REF] Martín | A multiconfiguration self-consistent field/molecular dynamics study of the (n→π*)1 transition of carbonyl compounds in liquid water[END_REF][START_REF] Sánchez | Solvent Effects by Means of Averaged Solvent Electrostatic Potentials: Coupled Method[END_REF] , the Effective Fragment Potential (EFP) [START_REF] Gordon | Fragmentation Methods: A Route to Accurate Calculations on Large Systems[END_REF] and X-Pol. [START_REF] Gao | Explicit Polarization: A Quantum Mechanical Framework for Developing Next Generation Force Fields[END_REF] An alternative to electrostatic-only embedding is FDET, where the two independent variables, namely: an auxiliary N A -electrons wavefunction Ψ A and the electron density of the environment ρ B (r), are both generated from the individual solute/solvent configurations. [START_REF] Neugebauer | An Explicit Quantum Chemical Method for Modeling Large Solvation Shells Applied to Aminocoumarin C151[END_REF][START_REF] Sapana | Nonuniform Continuum Model for Solvatochromism Based on Frozen-Density Embedding Theory[END_REF] These ensemble methods provide access to the inhomogeneous broadening of spectra, i.e broadening due to the distribution of solute/solvent configurations, and give some insight into the range of δν ex expected in the system that single-point calculations are unable to provide. These approaches are still not perfect, as most applications do not realistically account for homogenous broadening effects such as the aforementioned vibronic transitions and finite excited state lifetimes. Even though they provide valuable information on fluctuations, ensemble calculations can become prohibitively expensive when the embedded system is large or many solute/solvent configurations must be considered, necessitating a fall-back to single-point approaches.

The goal of the present study is assess the validity of substituting a single averaged solvent environment for an ensemble of environments. We wish to determine the fluctuations and uncertainties in properties derived from such calculations, and provide guidance on how to realistically construct an average environment. In essence, we wish to test whether or not an averaged environment can capture the features of an ensemble of calculations, an approximation we summarize as:

Ô Ψ A (R A )[B] ≈ Ô Ψ A ( RA )[ B] . (2)
In this expression, we define Ô[Ψ] as the quantum mechanical expectation value, Ψ| Ô|Ψ of some system property described by Ô. On the left-hand side, we take

Ψ A (R A )[B]
to be the wavefunction of fragment A, with instantaneous nuclear coordinates R A , embedded in instantaneous environment B, which could be described using any of the previously discussed methods. The indicates averaging over an ensemble of realizations of R A and B, and we will refer to this quantity as the ensemble average of Ô. On the right-hand side we have the same operator Ô, but instead of operating on every instantaneous realization of R A and B, Ô is evaluated once for a single wavefunction, Ψ A ( RA ) [ B]. Here, the instantaneous coordinates R A and environment B are replaced by an averaged set of subsystem coordinates RA embedded in the averaged environment B. This quantity will be referred to as the environmental average of Ô. The quantities RA and B may be ensemble averages, but are not required to be, as would be the case in a method such as PCM. If Eq. 2 is a reasonable approximation, we expect the environmental averages to fall well within the fluctuations of the ensemble calculations and be reasonably close to the actual ensemble average. To evaluate Eq. 2, we have performed classical MD simulations of the acetone (ACE) and water solute/solvent system in order to generate a realistic set of solute/solvent configurations for determining ensemble and environmental averages of various system properties as well as the first electronic transition of ACE. We chose to use the ACE/water system because the first excitation of ACE in water has been used extensively by our group and others as a benchmark for solvation models. [START_REF] Kaminski | Modeling Solvatochromic Shifts Using the Orbital-Free Embedding Potential at Statistically Mechanically Averaged Solvent Density[END_REF][START_REF] Laktionov | Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations[END_REF][START_REF] Bernasconi | Time dependent density functional theory study of charge-transfer and intramolecular electronic excitations in acetone-water systems[END_REF][START_REF] Aidas | Coupled Cluster Calculation of the n → π * Electronic Transition of Acetone in Aqueous Solution[END_REF][START_REF] Neugebauer | The merits of the frozen-density embedding scheme to model solvatochromic shifts[END_REF] The ensemble averages will be calculated from N independent MD frames, which can be written using the notation of Eq. 2 as:

Ô Ψ A (R A )[B] = 1 N N i=1 Ô Ψ A (R A )[B i ] = O N , (3)
where B i represents the environment of MD frame i out of a total of N frames, and O N a shorthand for the ensemble average which we will use for the remainder of this work. We consider the ensemble average the 'true' average of Ô and will characterize its fluctuations using the standard deviation, σ:

σ = 1 N N i=1 O N -Ô Ψ A (R A )[B i ] 2 . (4)
In order to keep our analysis focused on Eq. 2 as closely as possible, all molecules will be simulated with fixed molecular geometries. This allows us to replace RA in Eq. 2 with a single R A for the static ACE geometry and all the fluctuations we observe can be attributed to reorientation of the solvent environment. Complicating factors such as chromophore and solvent flexibility will be left for future studies. We then express the environmental averages as

Ô Ψ A ( RA )[ B] = Ô Ψ A (R A )[ B] = ŌN , (5)
with ŌN being our abbreviation for the environmental averages. The various strategies used for determining B will be discussed as they appear in the text.

In our statistical analysis of Eq. 2, we will consider the ensemble and environmental averages of four different system properties that are influenced by their solvent environment. The most fundamental of these is the electrostatic potential (ESP) generated by the solvent, which we study using the ESP evaluated at the nuclei of the ACE solute. We will then move on to three properties derived from the ESP: the atomic point-charges and dipole moments of ACE embedded in the solvent ESP determined using the ChELPG method, [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] and the first electronic transition of ACE. In each case we will compare the ensemble averages and their fluctuations with the single value generated by an averaged environment calculation. As mentioned previously, the fluctuations we observe will be a consequence of only the rotation and translation of solvent molecules due to our use of rigid solute and solvent molecules in our simulations.

Our results suggest that the ensemble averages of these four properties are well-re-produced by our environmental averaging methods. Although the average properties of the system are captured, properties of a single solvent configuration can become strongly decorrelated when different ESP approximations are applied to the system. We therefore advise caution when comparing the quality of calculations performed on a single solvent configuration under differing sets of approximations, as these differences between different approximations on the same configuration can become much larger than the differences between the ensemble and environmental averages. We conclude that that Eq. 2 is a valid approximation when the fluctuations of the properties in question are normally distributed, as we observed for the rigid solute/solvent simulations performed here.

2 Methods

The Embedding Energy

The embedded quantum system A (the solute) is characterized by the position of its M nuclei, R A , and their momentum, p A . Any instantaneous configuration of the mobile environment B (the solvent) is characterized by the position of all its atoms, R B , and their momentum, p B . The energy of the total A + B system can be written within the Born-Oppenheimer approximation as

E(R A , R B , p A , p B ) =K(p A , p B ) + V B (R B ) + E emb [v emb ](R A , R B ) (6)
with K(p A , p B ) being the classical kinetic energy of the nuclei, V B (R B ) is the self-energy of fragment B, and

E emb [v emb ](R A , R B )
is the energy of embedding A in B. Note that from now on, the dependence on the nuclear coordinates R A and R B will be kept implicit. To be consistent, the integration of the solvent equations of motion should be done with the full Hamiltonian, including both QM and molecular mechanics (MM) parts. [START_REF] Cao | On the Difference Between Additive and Subtractive QM/MM Caltulations[END_REF] Additionally, one should thus include the true AB quantum mechanical coupling using the full electronic density ρ A . This could be done in future studies, but for the time being we apply a simplified strategy. Instead of following the standard QM/MM approaches where both levels are evaluated sequentially or simultaneously with the proper partitioning of the MM and QM energies, [START_REF] Cao | On the Difference Between Additive and Subtractive QM/MM Caltulations[END_REF] we will perform classical MD simulations first, followed by QM calculations on the structures generated by classical simulations. As a consequence, the kinetic energy of the nuclei will be 0, and the total energy will depend only on V B and E emb [v emb ]. A proper way to make the statistics generated this way compatible with the appropriate statistics in the canonical ensemble will be given in Section 3.1.4.

Here, we suppose that the environmental self-energy term, V B , can be reduced to a molecular mechanics force field, i.e, V B = V M M B . Furthermore, to compute the solute-solvent electrostatic interaction, the electron density carried by each solvent molecule is gathered with their nuclei point charges to yield a simplified, effective point charge representation σB (r), that given by the MM force field, i.e.,

σB (r) = N solv j=1 q j δ(r -R B,j ), (7)
where the hat symbol designates discrete distributions. Then, the embedding potential is defined as the electrostatic interaction with these discrete charges

v emb (r) = v elst (r) = dr σB (r ) |r -r | , (8)
which yields the only-electrostatics embedding energy

E emb =H A + V elst = Ψ A | ĤA |Ψ A + drdr (ρ N A (r) + ρ A (r))σ B (r ) |r -r | (9)
with the solute nuclei charge density given by ρN

A (r) = M k=1 Z k δ(r -R A,k ).

Frozen-Density Embedding Theory

The ultimate goal is to go beyond electrostatic-only embedding to describe the interaction between fragment A and B, especially when studying electronic properties, as quantumconfinement effects may play an important role. For that, we turn into methods based on Frozen-Density Embedding Theory (FDET). [START_REF] Wesołowski | Embedding a multideterminantal wave function in an orbital-free environment[END_REF] Contrary to the electrostatics-only case, the embedding potential,

v emb [ρ A , ρ B ; v elst ](r) = v elst (r) + v nad xcT [ρ A , ρ B ](r), (10)
due to the non-additive exchange-correlation and kinetic term, v nad xcT [ρ A , ρ B ](r), is not homogeneous with the density ρ A , i.e., the energy is not equal to the expectation value of the potential. The embedding energy must then be evaluated with Ψ FDET A , the variational solution of the equation

ĤA + v emb Ψ FDET = λΨ FDET . (11)
The FDET embedding potential given in Eq. 10 and the total energy (E FDET v AB [Ψ FDET , ρ B ]) satisfy by construction the condition of self-consistency with the Hohenberg-Kohn energy functional (E HK v AB [ρ]), [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF][START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] meaning that

E FDET v AB [Ψ FDET , ρ B ] = E HK v AB [ρ FDET A + ρ B ] (12)
where ρ FDET A is the density corresponding to the wavefunction obtained in Eq. 11, for any admissible ρ B (r).

FDET interpretation of Eq. 2

In FDET, the expectation value of any quantum mechanical operator evaluated at the geometry

{R A }, {R B } i is a functional of ρ B i : Ψ FDET(R A ,R i B ) A [ρ B i ] Ô Ψ FDET(R A ,R i B ) A [ρ B i ] ≡ ÔR A ,R i B [ρ B i ] (13)
where

Ψ FDET(R A ,R i B ) A [ρ B i ]
denote the optimal embedded FDET wavefunction for a given ρ B i (r) and nuclear positions {R A }, {R B } i . The index i here indicates the i-th configuration of fragment B, from an ensemble of N geometries. The square brackets in Ψ

FDET(R A ,R i B ) A [ρ B i ] indi- cate that,
up to the unitary transformation, the FDET embedded wavefunction is uniquely determined by ρ B i for a given {R A }, {R B } i .

Therefore, B from Eq. 2 is defined in FDET by the set of {R B } i and ρ B i , which determine the local fields (functions of r in 3D), ρ B (r) and ρ nuc B (r). The fields represent the probability to find the negative (electrons) or positive (nuclei) charge in a given position r, respectively. The field ρ nuc B (r) can be easily transformed to the field of the corresponding electrostatic potential which is denoted with v nuc B (r) . Hence, equation 2 takes the following form:

ÔR A ,R i B [ρ B i ] N ≈ ÔR A , v nuc B [ ρ B ] (14)
The approximation in Eq. 14 has been used previously in combination with 3D-RISM and MD simulations to study solvatochromism in absorption [START_REF] Kaminski | Modeling Solvatochromic Shifts Using the Orbital-Free Embedding Potential at Statistically Mechanically Averaged Solvent Density[END_REF][START_REF] Zhou | Multi-scale modelling of solvatochromic shifts from frozen-density embedding theory with non-uniform continuum model of the solvent: the coumarin 153 case[END_REF][START_REF] Laktionov | Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations[END_REF] and emission. [START_REF] Shedge | Nonuniform Continuum Model for Solvatochromism Based on Frozen-Density Embedding Theory[END_REF] For the sake of this publication, we unify the notation and rewrite this key approximation in the same format used in the introduction. In this particular case, as mentioned before, we use a fixed geometry of the environment, which means that Ψ A (R A ) = Ψ A ( RA ), so equation 14 can be expressed as:

Ô Ψ A (R A )[ρ B ] ≈ Ô Ψ A ( RA , vnuc B )[ρ B ] (15)
2.3 Computational Technical Details

Molecular Dynamics

All classical molecular dynamics simulations were performed using DL_POLY_4.10. [START_REF] Todorov | New Dimensions in Molecular Dynamics Simulations via Massive Parallelism[END_REF][START_REF] Bush | A DAFT DL_POLY Distributed Memory Adaptation of the Smoothed Particle Mesh Ewald Method[END_REF] The intermolecular interactions were treated using a simple Lennard-Jones plus Coulomb force-field, and the simulated system consisted one acetone (ACE) molecule dissolved in 1000 SPCE waters. [START_REF] Berendsen | The Missing Term in Effective Pair Potentials[END_REF] The ACE structure and point charges were taken from Aidas et. al [START_REF] Aidas | Coupled Cluster Calculation of the n → π * Electronic Transition of Acetone in Aqueous Solution[END_REF] and the Lennard-Jones parameters from the ketone and methyl group parameters of the OPLS-AA force-field. [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] Lennard-Jones parameters for cross-species interactions were calculated using the Lorentz-Berthelot mixing rules. Electrostatics were handled using an Ewald summation and electrostatics and Lennard-Jones forces were cutoff at 15 Å. Trajectories were integrated using the velocity Verlet algorithm with a 2 fs timestep. All molecules were held rigid using the SHAKE algorithm with the default tolerance settings of DL_POLY_4.10. Initial conformations were generated using the fftool and PACKMOL 48 packages. The systems contained one ACE and 1000 SPCE molecules in an initial 32 Å box with cubic periodic boundary conditions. The following equilibration procedure was then followed. First, forces were minimized using the conjugate gradient method followed by a 25 ps simulation in the NVT ensemble. The temperature was 298.15 K and controlled using a Nosé-Hoover thermostat with a relaxation time of 0.5 ps. Next, a 50 ps simulation in the NPT ensemble was conducted in order to relax the system volume using a Nosé-Hoover barostat with 1 atm of pressure a relaxation time of 1.5 ps. Following the NPT simulation, the volume of the system was then fixed followed by a final 25 ps equilibration simulation in the NVT ensemble. A production simulation of 5 ns was then performed with frames saved every 1 ps, resulting in 5000 frames for later analysis. The production frames were then centered and wrapped around the ACE center-of-mass using the MDAnalysis package. [START_REF] Michaud-Agrawal | A Toolkit for the Analysis of Molecular Dynamics Simulations[END_REF] All images generated from the MD simulations were created using the VMD 1.9.3 package 50 and the Tachyon ray tracing library. [START_REF] Stone | An Efficient Library for Parallel Ray Tracing and Animation[END_REF] Scheme 1: Space-filling models and atom labels of all simulated systems.

Embedding Calculations

For the FDET calculations in this work, the energy and density of the core fragment A is evaluated with PySCF 52 using HF and the aug-cc-pVDZ basis set. The local-density approximation (LDA) was used for the non-additive functional and functional derivatives,

E nad xcT [ρ A , ρ B ], v nad xcT [ρ A , ρ B ](r)
, meaning the Thomas-Fermi kinetic energy functional, [START_REF] Thomas | The calculation of atomic fields[END_REF][START_REF] Fermi | Eine statistische Methode zur Bestimmung einiger Eigenschaften des Atoms und ihre Anwendung auf die Theorie des periodischen Systems der Elemente[END_REF] the Slater exchange functional [START_REF] Slater | The Theory of Complex Spectra[END_REF] and the Vosko-Wilk-Nusair fit of Eq. 4.4 for the correlation energy of the uniform electron gas (VWN5). [START_REF] Vosko | Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis[END_REF] The AO contraction of the embedding potential was done using functions from PySCF and our Python libraries FDET-Average 57 and FDETaco. 58 The ChELPG charges discussed in subsection 3.1.3 were calculated with the new single-fragment module in Q-Chem 5.4 [START_REF] Epifanovsky | Software for the frontiers of quantum chemistry: An overview of developments in the Q-Chem 5 package[END_REF] from densities obtained from doing linearized-FDET [START_REF] Prager | First time combination of frozen density embedding theory with the algebraic diagrammatic construction scheme for the polarization propagator of second order[END_REF] embedding with HF/aug-cc-pVDZ and adding a Møller-Plesset first-order correction to them. Finally, the excitation energies were obtained also with the linearized-FDET approximation, with the reference density, ρ ref A , being the ground-state HF/aug-cc-pVDZ density of the isolated acetone. For frame-wise calculations with FDET, the density of the environment was obtained from the superposition of molecular densities at the HF/augcc-pVDZ level of theory. All excitation energies were computed with the ADCC python module, [START_REF] Herbst | adcc: A versatile toolkit for rapid development of algebraic-diagrammatic construction methods[END_REF] following the FDE-ADC(2) procedure as described in Prager et al . [START_REF] Prager | First time combination of frozen density embedding theory with the algebraic diagrammatic construction scheme for the polarization propagator of second order[END_REF] 3 Results and Discussion

Electrostatic Properties

Solvent Systems and Environmental Averaging Methods

We take this opportunity to not only test Eq. 2, but to also evaluate different sets of approximations used for building the solvent environment from the MD simulation frames which we will call Direct/Cube, Direct/Supercube, Grid/Cube, and Grid/Sphere. In the Direct/Cube method, each MD frame is first centered around the ACE center-ofgeometry1 and the solute coordinates wrapped so that all atoms of the cubic system are within a box centered at ACE and maintaining the system's periodicity. Next, we directly calculate the ESP at the nucleus of each ACE atom using a simple pair-wise sum over all solvent atoms:

v elst (r n ) = N solv j=1 q j |r n -r B,j | , ( 16
)
where r n in this case is the position of a particular solute nucleus, N solv the number of solvent atoms, and q j the charge of solvent atom j at position r B,j . This process is repeated for every frame and the values of v elst (r n ) from each frame are used to build probability distributions and ensemble averages. Using this method, the ensemble and environmental average ESPs are by definition exactly the same. Because there are no approximations made in calculating the ESP in this method, we consider Direct/Cube our reference method.

In order to determine if the 1000 SPCE water system was sufficiently large for the ESP to converge, we created the Direct/Supercube system, in which 26 copies of the solvent box were tiled around a central box. The ESP experienced by the central ACE molecule was then calculated and compared to that of the single 1000 solvent molecule box. The pair-wise calculation in Direct/Cube can be extremely expensive due to poor scaling of Eq. 16 with N solv and the number of grid-points {r n }, particularly when evaluating v elst on a dense integration grid, as it is done in the next subsections to calculate derivative properties. To address this limitation, we add one approximation in what we call the Grid/Cube method. The MD frames are prepared in the same manner as in Direct/Cube, but are then rotated such that each frame is aligned to the ACE nuclei. Next, solvent atoms and their charges are compacted inside 0.025 nm cubic volume elements using a 3D histogram, illustrated in Scheme 2. The ESP at the solute nuclei is then calculated according to

v elst (r n ) = N grid l=1 q l |r n -r l | , ( 17
)
where N grid is the number of grid points, q l the average of all solvent charges compacted into into volume element l, and r l the center of the volume element. This cubic grid ESP is then extrapolation over the molecular Becke-Lebedev grid using the nearest-neighbor method.

The cubic grid is dense enough such that each point typically contains only 1 atom per MD frame, therefore no efficiency is gained using the grid method for determining the ensemble average. On the other hand, the environmental average ESP calculation is much more efficient, as all solvent atoms from the 5000 frames are compacted into grid points which significantly reduces the total number of electrostatic interactions to calculate. There are two main differences in this approximation compared to the that in Ref. 10. First, we build an evenly distributed cubic grid instead of using the minimal number of charges placed on spherical shells around the solute with optimized radius. Secondly, we do not optimize the charges, they are only compacted according to their position.

Scheme 2: The Grid method for averaging. For each volume element, the compacted charge q k used to evaluate the averaged environment density in Eq. 25, is the combination of all the charges that appear inside the volume element v α throughout the trajectory. The size of the cube allows only one atom per volume element in each frame.

Unfortunately, the rotation required to align each frame in Grid/Cube results in a nonuniform solvent density for distances more than 1/2 a box length away from the ACE nuclei while also destroying the periodicity of the system. Therefore, physically meaningful environmental averages cannot be extracted from the Grid/Cube system, and consequently, they are not calculated. To alleviate this problem, we employ the Grid/Sphere method, were we only consider a sphere of water molecules whose oxygen atom is ≤ 1 /2 a box length from the ACE center-of-geometry. Our hope is that these spherical clusters exhibit the same properties as the Grid/Cube reference system. Illustrations of Direct/Cube, Di-rect/Supercube, and Grid/Sphere are shown in Scheme 3. Note that all environmental averages require alignment of the individual frames, and for the Direct approaches, the electrostatic potential is evaluated as in Ref. 12, using the charges of all water molecules over the 5000 frames, scaled by a factor of 1/5000.

Electrostatic Potentials

Distributions of the ESP at the ACE nuclei, r n = R A,n , and the corresponding environmental averages from the Direct/Supercube, Direct/Cube, Grid/Cube, and Grid/Sphere methods are shown in Figure 1. The ensemble averages, standard deviations, and environmental averages are plotted in Figure 2 and tabulated in Table ??. These data demonstrate that the ESP distributions of all the ACE nuclei are roughly Gaussian and all chemically equivalent atoms have the same distributions and ensemble average. We also observe negligible differences between the Direct/Supercube and Direct/Cube methods, allowing us to conclude that the 1000 water molecule box was sufficiently large to converge the system. The standard deviations of the Grid/Cube distributions are roughly 30% larger than the Direct/Cube reference system, whereas the Grid/Sphere distributions are only 20% larger. We attribute this discrepancy to the zero-th order extrapolation scheme that is used to project the discrete solvent charges on the 3D-grid; higher order schemes like B-splines would be more appropriate. The environmentally averaged ESPs are in strong agreement with their ensemble counterparts and fall well within the standard deviation of the ensemble averages. We also found that the ESP distributions, ensemble averages, and environmental averages of the Grid/Sphere method were all shifted by about -0.29 AU from those of the Direct/Cube and Grid/Cube methods. This offset can be attributed to the differences in boundary conditions between the periodic cubic and non-periodic spherical systems, which results in a constant shift in the ESP. As we shall see, the ESP offset is constant, i.e. it is the same at every position, explaining why this difference does not affect any of the derivative properties. See the Appendix A for a complete explanation of the origin of the shift. These results demonstrate that there is no significant difference between the ensemble and environmentally averaged ESPs derived from Direct/Cube reference method and those of the Grid/Sphere method results, which allows us to conclude that both the grid approximation and the spherical cluster are acceptable approximations for determining the ESP experienced by a solute. These approximations effect only a minor increase in the width of the ESP probability distributions while giving the same ensemble averages. The environmental averaging techniques we employed resulted in average ESPs that are equivalent to the full ensemble calculations. Scheme 3: Cubic (top) and spherical (bottom) systems.

ChELPG Charges and Molecular Dipole Moments

The first derivative properties we will consider are the ChELPG charges and dipole moments of the ACE solute determined in the presence of the solvent ESP. The solute electron densities used to evaluate the electrostatic potential in the ChELPG procedure are obtained from a Hartree-Fock optimization of ACE in the presence of the solvent ESP, v elst (r), as in Eq. 9. We then performed the ChELPG fitting procedure for each individual frame of the MD simulations in order to determine the probability distributions and ensemble averages,

Q ChELPG A
, of the point-charge of each ACE atom. This procedure was carried out for ESPs derived using the Direct/Cube Grid/Cube and Grid/Sphere methods described in Section 3.1.2. The environmentally averaged charges, QChELPG A , were determined from a single ChELPG calculation with ACE embedded in the environmentally averaged ESP determined using either the Direct/Cube or Grid/Sphere methods. ??. As with the ESPs, all ChELPG charge distributions are Gaussian and chemically equivalent atoms have the same distributions of point charges. The ensemble and environmental point charges from all systems are practically identical, just as was seen for the ESPs, but the standard deviation of the Grid/Cube distributions are, for most atoms, a factor of 2 larger than Direct/Cube. This is in contrast to the Grid/Sphere standard deviations which are nearly the same as Direct/Cube. The offset in the ESP observed in the Grid/Sphere method did not have an effect on the distribution of ChELPG charges. We now can be confident that ESPs derived from the grid and spherical solvent cluster approximations do not siginficantly change a solute's ensemble and environmental average point-charges compared to a direct ESP evaluation. As another test of the quality of the ChELPG method and our averaging procedures, we calculated the dipole moments generated from the quantum mechanical (QM) charge density resulting from the Hartree-Fock optimization of ACE and the classical (CL) dipole moments calculated from the fitted ChELPG point charges. The distribution of dipole moments from the ensemble of calculations and the corresponding environmental averages using the Direct/Cube Grid/Cube and Grid/Sphere methods are shown in Figure 5. The ensemble averages, standard deviations, and environmental averages can be found in Figure 6 and Table 1. The ensemble and environmental average dipole moments are nearly identical in all systems, whereas the standard deviation of the Grid/Sphere distributions are 15% narrower than their Direct/Cube and Grid/Cube counterparts. We can conclude from these data that the ChELPG method does an excellent job of reproducing the average quantum mechanical solute dipoles as well as their distributions. Additionally, the grid and sphere approximations do not significantly affect the dipole moment of the solute. Although the QM and CL dipole moment averages and distributions are the same for all three systems, that does not guarantee that the dipole moments of individual frames are also the same. Confirming this correlation is critical for evaluating differences between calculations performed with different approximations on the same MD frame. The top panel of Figure 7 shows the correlation between the QM and CL dipole moments for Direct/Cube, Grid/Cube, and Grid/Sphere. These data confirm that the ChELPG procedure faithfully reproduces the dipole moments for every frame/ESP method combination considered here. On the other hand, correlations between frames treated with different ESP approximations are not as strong (bottom panel Figure 7). These data show that, for a single frame, the Direct/Cube Grid/Cube and Grid/Sphere methods result in significantly different dipole moments and, by inference, ChELPG charges and ESPs. We attribute this decorrelation to two factors: 1) the slight movement of solvent atoms from their 'true' positions when the gridding procedure is applied and 2) the removal of solvent molecules when constructing the spheres from the cubes. When either procedure 1) or 2) are applied to the original solvent cube they produce a 'new' set of solvent coordinates that are not exactly the same as the original cube of solvent. Due to the sensitive nature of the ESP, frame-to-frame correlation is lost, which manifests later in decorrelated ChELPG charges and dipole moments. Even though the within-frame correlation is lost, no effect is seen in the average properties or their distributions, suggesting that the ensemble properties of the system are conserved even though individual frame correlation is not. This effect must always be kept in mind when applying different approximations to individual frames. Although these approximations do not affect the ensemble and environmental averages, evaluating approximations by comparing results from a single MD frame can be misleading. 5 and bars in Figure 6). The columns labeled 'σ / D' are the standard deviations of the ensemble averages (error bars in Figure 6). Rows labeled 'QM' are the dipole moments calculated from the quantum mechanical charge density and 'CL' from the CHELPG charges.

Direct/Cube

Grid/Cube Grid/Sphere

|M | / D σ / D | M | / D |M | / D σ / D | M | / D |M | / D σ / D | M | / D QM 4

The Separated QM/MM Time-Averaging Scheme

The electrostatic interaction energy defined in Eq. 9 can be written in two equivalent ways, as a function of the embedding electrostatic potential created by B on A, v els (r) in our previous notations, or reversely the external electrostatic potential created by A on B, v A (r; [ρ A ]), i.e.,

V elst = dr v elst (r) (ρ N A (r) + ρ A (r)) (18) = dr v A (r; [ρ A ]) σB (r) dr v A (r; Q ChELP G k [ρ A ]) σB (r) (19)
N j=1 q j M k=1 Q ChELP G k |R B,j -R A,k | . ( 20
)
In the near equality, the electrostatic potential due to the embedded system charge density has been approximated by the one created by a set of point charges located on the nuclei, that are best fitted according to an appropriate scheme. As mentioned above, we have chosen the electrostatic-potential based ChELPG scheme. Note that if the QM method used to obtain Ψ A is variational, it can be shown easily (Appendix B) that the appropriate molecular dynamics scheme to handle the Hamiltonian defined by Eq. 9, with energy conservation only limited by the near equality of Eq. 19, is typically that for a polarizable system, obtained with the following procedure: 1) for a given configuration R B , compute Ψ A and ρ A by minimization of Eq. 9 using Eq. 18, monitor the 'gas phase + polarization' energy Ψ A | ĤA |Ψ A , and extract the set of charges

{Q ChELP G k [ρ A ]}, 2
) compute the forces on the environment atoms using those point charges and Eq. 20, and 3) update the environment positions R B and iterate. Such algorithm is indeed standard in the literature, but its implementation in our analysis will be left for future investigations.

In the case of a completely separated MM and QM calculation scheme, as in this work, the MM simulations are performed with a reference non-polarizable system to generate afterwards the correct statistics for the polarizable system itself. To this end, one can define a reference, fixed set of charges {Q 0 k }, and a reference Hamiltonian H 0 where the electrostatic interaction seen by the environment in Eq. 20 is replaced by that with the fixed charges {Q 0 k }. The polarizable system Hamiltonian can be written as

H = H 0 + ∆V elst (21) = H 0 + N j=1 q j M k=1 Q ChELP G k -Q 0 k |R B,j -R A,k | . ( 22
)
Then, in the canonical thermodynamic ensemble, the canonical average value of any observ-able of the polarizable system can be computed as

Ô(R B ) = Ô(R B ) exp(-β∆V elst (R B )) 0 exp(-β∆V elst (R B )) 0 , (23)
where the subscript 0 indicates a dynamics performed with the reference Hamiltonian H 0 , i.e., with the same self-consistent scheme as above but with the forces computed with the fixed charges {Q 0 k } in step 2. The biased ensemble averages (calculated while neglecting the exponentials in Eq. 23) and unbiased ensemble average (calculating using Eq. 23) for the ChELPG charges of all ACE atoms from the 5000 frame Direct/Cube and Grid/Sphere systems provided in Table 2, in addition to the initial ChELPG charges used in the MD simulations. The variations we see here between the biased and unbiased ChELPG charges are smaller than the fluctuations due to the environment. Therefore, we are confident in only considering the biased averages in future calculations. This conclusion by be different, though, if the initial ChELPG charges were drastically different from the final charges.

Conclusions

We find that the grid method for determining the ESP and spherical solvent cluster approximation have no significant effect on the ensemble and environmental average ESP, ChELPG charges, or solute dipole moments. On the other hand, these approximations are less effective when applied within the same frame, which was attributed to differences in the solvent cluster size and shifting of solvent nuclei after applying the grid procedure. This means that the quality of the approximations can only be evaluated by looking at ensembles of calculations. Single-frame agreement should be improved by using a higher-order charge extrapolation scheme and/or finer grids. This would be at the expense of an increased computation time for the 3D electrostatic embedding potential acting on the solute. Additionally, unbiasing the trajectories did not effect the average system properties, but this may not hold for systems in which the simulated solute charges are drastically different than those of the embedded solute. Therefore, we conclude that Eq. 2 holds for the ESP, ChELPG charges, and dipole moments and that Grid/Sphere set of approximations reasonably reproduces the electrostatic properties of the Direct/Cube reference system. We can now confidently use ESPs determined using the Grid/Sphere method for computing electronic structure properties.

Acetone Electronic Transitions

Next we will test the applicability of Eq. 2 to solute electronic transitions by calculating the first electronic transition of ACE embedded in different descriptions of the solvent environment. For all transition energy calculations, all obtained with the ADC(2) method, [START_REF] Dreuw | The algebraic diagrammatic construction scheme for the polarization propagator for the calculation of excited states[END_REF] the reference is the first excitation energy of an isolated ACE molecule in vacuum.

One set of calculations will be performed using point charge embedding (PCE), where the ACE solute is embedded in the ESP generated by the solvent point charges taken directly from the MD forcefield (that of SPCE water). For the ensemble averages, the embedding ESP is determined from a sphere of solvent molecules with a radius of 1 /2 the box length. The ESP is then determined using the pair-wise method of Eq. 16. This method combines the direct pair-wise evaluation of the ESP in the Direct/Cube method of Section 3.1.1 with the spherical solvent cluster approximation. Excitation energies for each frame are then calculated. For the environmental average, the embedding ESP was calculated using the Grid/Sphere method and a single-point excitation energy calculation. These procedures for determining the ensemble and environmental averages were chosen in order to compare our computationally cheap Grid/Sphere method of ESP determination with a very expensive, yet very exact, implementation of the PCE procedure.

We also calculated second set of ACE excitation energies using FDET embedding. This method combines the electrostatic contribution from the solvent with the non-electrostatic (i.e. Pauli exclusion) portion of the energy. The ensemble averages are determined by first taking a sphere of solvent molecules with a radius of 1 /2 a box length centered at the ACE solute. A quantum mechanical charge density is then determined at the HF/aug-cc-pVDZ level of theory from an inner-sphere of solvent molecules 1 /4 of a box length from the ACE. This superposed density is then used in the evaluation of the FDET embedding potential. All solvent molecules farther than 1 /4 a box length distance from the ACE only contribute to the electrostatic portion of the embedding energy. This procedure is performed for each of the 5000 MD frames in order to construct the ensemble average FDET excitation energy. For the FDET environmental averages, the embedding charge density is evaluated using only the positions of the solvent nuclei and their forcefield charges; no quantum mechanical embedding density is used in order to maximize computational efficiency. We begin by gridding the solvent nuclei as in the Grid/Sphere method of ESP determination. The electronic density is obtained by dividing the effective electronic charge q e k , that is the effective number of electrons associated to the atom type k, by the volume element formed by the grid. As we need to account for all atoms over the ensemble of solvent geometries, instead of tracking individual atoms we group them into atom types. Here, it is important to stress that q j from the definition of σB (Eq. 7), is the total effective charge of the j-th atom, taken from the force field. So, if the j-th atom is of the type k, then q e j(k) = -q j + Z j(k) (24) where Z j(k) is the nuclear charge of the atom j, of type k. To simplify the notation, from now on we only use k.

The average solvent density at the position of the point with index α (r α = {x α , y α , z α }), as in our previous work, [START_REF] Laktionov | Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations[END_REF] is evaluated as

ρ B (α) = k q e k g k (α), (25)
here, the division over the volume elements is contained inside the solvent distribution function g k (α). For a given atom type k in the volume element corresponding to the point

α: g k (α) = n α,k N ν α ; k,α g k (α)ν α = N nuc , ( 26
)
where N is the number of configurations (MD-frames), N nuc being the total number of nuclei in the box, n α,k the number of occurrences of an atom of type k in the volume element ν α , and the total volume of the box V = α ν α . Therefore, the average electrostatic potential is computed as

v elst [ ρ B ](r) = k,α (-q e k + Z k )g k (α)ν α |r -r α | (27)
We emphasize that this ESP is equivalent to the Grid/Sphere method described in Section 3.1.1. Whereas for the PCE the embedding potential is just v elst (r), for the FDET calculations the embedding potential is:

v emb [ρ A , ρ B ; v elst ](r) =v elst [ ρ B ](r) + v nad(lin) xcT [ρ A , ρ B ](r). (28)
The environmental average excitation from FDET is then obtained from a single-point calculation in the presence of the potentials of Eq. 28.

For comparison with the single-point environmental averages, we calculate the ensemble average transition energies using two approaches. First, we take the weighted average of the transition energies:

ν ex = N i=1 w i ν (i) ex , (29)
where ν

(i)
ex is the excitation energy of frame i, and the weight, w i , is determined from the oscillator strength, f i , according to

w i = f i N i=1 f i . ( 30
)
The weighted standard deviation of ν ex is calculated as:

σ wt = N i=1 w i ν i -ν 2 . (31)
Note that Eq. 3.2 is only valid when i w i = 1. We also calculate the unweighted average of ν ex ,

ν ex d = 1 N N i=1 ν (i) ex . ( 32
)
which is the version of an ensemble average closest to the environmental average, where a single calculations is performed. Note that this expression assumes that the excitation energies can be decoupled from their oscillator strength (hence the subscript d). If this assumption holds we expect to observe ν ex = ν ex d .

Figure 9: Electronic transitions of ACE in SPCE water evaluated using PCE and FDET embedding methods. Due to the symmetry of the transition, the oscillator strength of the isolated ACE and environmentally averaged calculations is effectively 0. The ACE transitions from these three sets of calculations are plotted in Figure 9, and all transition energy averages and single-point calculation results are given in Table 3. The environmental average transition energies are in excellent agreement with the ensemble aver-ages, to within 0.2%. One can also observe that the difference between the PCE and FDET methods is rather small. This is consistent with previous observations which found that the magnitude of the shift can be primarily attributed to the solute-solvent electrostatic interaction [START_REF] Laktionov | Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations[END_REF][START_REF] Aidas | Coupled Cluster Calculation of the n → π * Electronic Transition of Acetone in Aqueous Solution[END_REF] . Additionally, the small difference between ν ex d and ν ex shows that decoupling the oscillator strengths and excitation energies is a good approximation, applicable to model the average behavior of the transitions. However, in terms of frame-to-frame correlation, similar to the dipole moments, the single-frame oscillator strengths and transition energies are not well correlated with each other when different sets of approximations are applied (see Figure 10).

Sample Size Dependence

At the outset of this study, we arbitrarily chose to use 5000 MD frames in the hope that this was a sufficiently large sample for properly converging the ensemble and environmental averages. In future applications involving high level theories, 5000 frames may be computationally intractable and a smaller number of calculations must be considered. We now look to determine the number of calculations needed in order to be 'reasonably' confident in the average excitation shift generated by an ensemble of calculations. For this section, we will consider only the solvatochromic shift, δν ex = ν emb -ν iso , i.e. the shift between the isolated and embedded transition energies. We use this quantity, instead of the excitation energies, as we are interested in the accuracy of the description of the effect of the environment, and not the accuracy of the QM method itself.

First we will consider the number of frames required to converge an ensemble average shift, δν ex . We begin by splitting the 5000 frame FDET/Sphere trajectory into n samp samples consisting of n frame sequential frames, where n samp * n frame = 5000. The sample's mean excitation shifts were calculated as in Section 3.2, together with the weighted standard deviation. We will call these individual sample calculations the 'sample means' and 'sample standard deviations.' After processing each sample, the mean of the sample means and its 95% confidence interval were calculated, which we will refer to as the 'aggregate mean' and 'aggregate standard deviation.' The aggregate standard deviation gives insight into the uncertainty of δν ex when it is calculated from n frame individual calculations. Although we have confirmed that the ensemble calculations behave normally, this does not necessarily mean that the same conclusion can be made for the environmental averages.To characterize the uncertainty in νex , we broke down the 5000 frame MD trajectory into samples of 100, 500, and 1000 sequential frames. An average environment using the Grid/Sphere method from each sample and an FDET calculation was performed in the same manner as in Section 3.2. Figure 13 shows the resulting environmentally averaged transition energy shift from each of these samples, and the aggregate means and standard deviations of δν ex are given in Table 4. We find that in order to ensure an uncertainty in δν ex of ≤ 5 meV, a minimum of 500 MD frames must be used to construct the average solvent environment. Note that this is a lower bound on the uncertainty, as we have not included any solute or solvent flexibility or polarizability in our simulations. Our result is roughly in agreement with the coupled-cluster with single-and double-excitations reported in Ref. 39, where it was shown that from 400 MD configurations, the standard deviation of ensemble averages of the first excitation energy of acetone was smaller then 100 cm -1 (approx. 12 meV).

Conclusions

The results of our study demonstrate that Eq. 2 is a valid approximation for a rigid solute dissolved in a rigid solute. The the three different methods used for treating the solvent environment, Direct/Cube Grid/Cube and Grid/Sphere, were found to give the same ensemble averages of the ESP, ChELPG charges, and solute dipole moments. Strong correspondence was found between these ensemble averages and their corresponding environmental averages. When reducing the system from a cube of solvent to a sphere, we observed a uniform shift of -0.29 eV in the ESP due to the change in the boundary conditions. Because this shift was constant for all points in space, it did not affect any properties derivative of the ESP. We found that the most efficient method to get environmental averages, Grid/Sphere, is valid for determining averages, but frame-to-frame correlation is lost due to the shifting of solvent atoms inherent in the grid approximation. This correspondence can only be recovered by using a very fine grid, so that the potential converges to the exact solution. Therefore, the Grid/Sphere method is not recommended when looking at properties that depend on the instantaneous fluctuations of the solvent or the solute. Because we chose to perform uncoupled QM/MM calculations, the effect of the choice in the charges used for the classical simulation was explored. We found that there was no significant difference in the ensemble averages before and after an unbiasing procedure was applied.

When studying the first vertical excitation of acetone, as in this particular case, the shift in the excitation energy due to the solvent is mainly the result of the electrostatic interactions, it was not a surprise finding that the results from PCE and FDET were very similar. The ensemble and environmental averages were found to be very close to each other in both cases: 0.233 eV and 0.223 eV with PCE, and 0.223 eV and 0.227 eV with FDET, proving once more the validity of Eq. 2. Finally, we looked at the dependence of the shift in the excitation energy on the number of independent frames used for constructing the ensemble and environmental averages. We found that at least 775 frames are needed for the ensemble averages to have an uncertainty of < 5 meV, whereas the environmental averages required at least 500 frames.

The immediate following steps in this project are to: 1) explore other approaches to construct environmental averages, for instance molecular density functional theory, and 2) investigate the applicability of environmental averages in cases where the solute is flexible and where the construction of an average geometry is required. [START_REF] Ricardi | Embedding non-rigid solutes in an averaged environment: a case study on rhodopsins[END_REF] changing here E in H, and knowing indeed that there is a parametric dependence of ρ A (r) on both R A and R B , i.e. ρ A (r) ≡ ρ A (r; R A , R B ). V M M (R A , R B ) is the MM force field that describes the mutual interaction of the solvent molecules and their non-electrostatic interaction with the solute, modelled by site-site Lennard-Jones contributions. The force on any atom j of the solvent can be written as

∂E ∂R B,j = -dr δH A δρ A (r) + v elst (r) ∂ρ A (r) ∂R B,j (37)
-dr (ρ A (r) + ρN A (r))

∂v elst (r) ∂R B,j - ∂V M M ∂R B,j
The bracket in the first term is just the Euler-Lagrange equation for the quantum system and it vanishes for all r's. So does the whole integral. The second term is the electrostatic contribution to the force, F elst B,j . Accounting from the equalities 17-19 in the main text it can be rewritten as

F elest B,j dr M k=1 Q ChELP G k δ(r -R A,k ) N j=1 q j (R B,j -r) |r -R B,j | 3 = N j=1 q j M k=1 Q ChELP G k (R B,j -R A,k ) |R A,k -R B,j | 3 = - ∂ ∂R B,j N j=1 q j M k=1 Q ChELP G k |R A,k -R B,j | = - ∂ ∂R B,j N j=1 q j v A (R B,j ; Q ChELP G k [ρ A ]) (38)
This implies that once the quantum Hamiltonian has been minimized, yielding the ground state electronic density ρ min A , the total energy of the QM/MM system is given by

H(R A , R B , p B ) =K(p B ) + H A [ρ min A (r)] + N j=1 q j v A (R B,j ; Q ChELP G k [ρ min A ]) + V M M (R A , R B ), (39)
and the force are given consistently by Eq. 38 above. This leads to the self-consistent molecular dynamics algorithm described in the text. As mentioned there too, one can define at this point a reference classical Hamiltonian H 0 for which the energy is defined exactly as above but with a fixed set of solute charges {Q 0 k }. One can then decompose the total energy into H = H 0 + ∆V elst as in Eqs. 20-21. The following Eq. 22 is a standard relationship, exact in the canonical thermodynamic ensemble, which enables to compute the property of the system under study from a dynamics performed with the reference Hamiltonian. Even if exact in theory, the formula is only applicable in practice if the reference system remains close to the exact one (here if the reference charges are a good guess with respect to the ChELPG charges measured on average), or if the statistics in ∆V elst is Gaussian (a seemingly good approximation in our study). In that case, the following formula can be used

Ô(R B ) = Ô(R B ) 0 -2β δ∆V elst (R B )δ Ô(R B ) 0 , (40)
where δX = X = X indicate the fluctuation.

Figure 1 :

 1 Figure 1: Probability distributions and and environmental averages (sticks) of the ESP experienced by the ACE nuclei for the Direct/Supercube, Direct/Cube, Grid/Cube, and Grid/Sphere systems.

Figure 2 :

 2 Figure 2: Ensemble and environmental averages of the ESP at all ACE nuclei. Error bars are the standard deviation.

Figure 3 :

 3 Figure 3: Probability distributions and and environmental averages (vertical lines) of the ChELPG charges of all ACE atoms for the Direct/Cube, Grid/Cube, and Grid/Sphere systems.

Figure 4 :

 4 Figure 4: Ensemble and environmental averages of the ChELPG charges for all ACE nuclei. Error bars are the standard deviation.

Figure 5 :

 5 Figure 5: Distribution of the ACE quantum mechanical (QM, solid lines) dipole moments and the classical dipole moment derived from the ChELPG point-charges (CL, dashed lines) of the ACE nuclei from the Direct/Cube, Grid/Cube, and Grid/Sphere systems. The distributions have been offset vertically for clarity. Vertical lines are the environmental averages.

Figure 6 :

 6 Figure 6: Ensemble and environment averages of the quantum mechanical (QM) and classical (CL) dipole moments. Error bars are the standard deviation.

Figure 7 :

 7 Figure 7: Top: Correlation between the QM and CL dipole moments in the Direct/Cube, Grid/Cube, and Grid/Sphere systems. Bottom: Correlation between the Grid/Sphere CL dipole moments and the Direct/Cube and Grid/Cube systems.

Figure 8 :

 8 Figure 8: System used in the frame-wise FDET-embedding calculations. The internal sphere of water molecules (ball-and-stick models) are treated as densities and used to evaluate the FDET-embedding potential. Water molecules beyond this layer are treated as point-charges.

Figure 10 :

 10 Figure 10: The correlation between the PCE and FDET oscillator strengths (top) and transition energies (bottom) for the first ACE electronic excitation. The black line represents a perfect 1:1 correlation.

Figure 11 :

 11 Figure 11: Break down of the full 5000 FDET/Sphere transition energy shifts, ν = N i=1 w i ν i , into smaller sequential samples of varying size.

Figure 12 :

 12 Figure 12: Top: Aggregate means of the ensemble average transition energy shift δν ex vs. n samp . Error bars are the aggregate standard deviation. Bottom: Aggregate standand deviations vs. n samp . The red curve is a fit to a/ √ n samp where a = 139 meV.

Figure 13 :

 13 Figure 13: Environmentally averaged transition energy shifts, δν ex = νemb -ν iso , constructed from a varying n frame number of sequential simulation frames.

Table 1 :

 1 Ensemble, |M | , and environmental, |M |, average ACE dipole moments (vertical lines in Figure

Table 2 :

 2 Initial, biased, and unbiased ensemble average ChELPG charges of ACE from the Direct/Cube and Grid/Sphere systems.

				Direct/Cube			Grid/Sphere
	atom Initial / e Biased / e Unbiased / e % Diff. Biased / e Unbiased / e % Diff.
	CC	0.7800	0.8759	0.8883	1.4%	0.8743	0.8823	0.9%
	CT1	-0.4926	-0.4885	-0.4706	-3.7%	-0.4890	-0.4743	-3.0%
	CT2	-0.4926	-0.4880	-0.5210	6.8%	-0.4878	-0.5148	5.5%
	HI1	0.1338	0.1437	0.1398	-2.7%	0.1441	0.1385	-3.9%
	HI2	0.1338	0.1440	0.1503	4.4%	0.1437	0.1506	4.8%
	HO1	0.1423	0.1412	0.1344	-4.8%	0.1412	0.1350	-4.4%
	HO2	0.1423	0.1410	0.1550	9.9%	0.1408	0.1442	2.5%
	HO3	0.1423	0.1407	0.1281	-8.9%	0.1404	0.1404	0.0%
	HO4	0.1423	0.1409	0.1439	2.1%	0.1404	0.1558	11.0%
	OA	-0.6317	-0.7508	-0.7481	-0.4%	-0.7482	-0.7577	1.3%

Table 3 :

 3 Ensemble and environmental first excitation energies of ACE. The ensemble averages (Eq. 3) are weighted by the oscillator strength. Shifts are calculated with respect to the isolated ACE. All values are in eV.

		ν ex d δν ex d	ν ex	δν ex	σ wt	νex	δν ex
	Isolated	-	-	-	-	-	4.270	-
	PCE	4.480 0.210 4.503 0.233 0.105 4.493 0.223
	FDET 4.478 0.208 4.493 0.223 0.093 4.497 0.227

Table 4 :

 4 Aggregate mean and standard deviation of the environmentally averaged transition energy shifts.n samp n frame δν ex Agg. Mean / meV δν ex Agg.

	Stdev / meV

Although our solute is rigid, it is permitted to freely translate and rotate during the simulation

Acknowledgement D.G. thanks the Department of Physical Chemistry at University of Geneva for invited Professor position. All calculations were performed on the Baobab and Yggdrasil clusters at the University of Geneva.

Appendix A: Potential shift due to the change from P-type to M-type boundary conditions Solvent molecules located at far distances from the solute molecule perspective, where the solute electric field is small, appear to be rotating freely. Therefore, on average, solvent molecules at the borders can be seen as averaged charges or charged spheres. For instance, water molecules are reduced to an isotropic quadrupole (γ = 2q H d 2 ) and higher isotropic multipoles.

In the language of Kastenholz, 64 the periodic boundary conditions (PBC) used in MD and Molecular Density Functional Theory (MDFT) are labeled as P-boundary type. In this case, the average charge density at every point in the box is zero, and there are no discontinuities at the borders. On the contrary, M-type boundary conditions, where molecules are taken as a whole at the boundaries, give a different picture -due to the split of molecules there is a surface polarization P M at the borders. According to Ref. 64, the polarization at the surface

where ρ is the liquid density, creates inside the box a uniform shift ∆φ of the potential,φ M = φ out -∆φ. In the infinite fluid case,

For SPCE water model, d = 1 Å, ρ = 0.03328 mol /Å 2 , q H = 0.4238 e, and γ = 0.8476, therefore ∆φ = 2π 3 ργ ≈ 0.059 e/Å ≈ 0.031 e/Bohr.

Appendix B: The Separated QM/MM Time-Averaging Scheme Supposing for simplicity that the solute A is kept immobile and the embedding potential reduced to its electrostatic component, and adopting here a DFT framework, the total energy of the A + B system in Eq. 6 can be expressed as