Frozen-density embedding theory with average solvent charge densities from explicit atomistic simulations

Abstract

Besides molecular electron densities obtained within the Born–Oppenheimer approximation (ρ_B(r)) to represent the environment, the ensemble averaged density (〈ρ_B〉(r)) is also admissible in frozen-density embedding theory (FDET) [Wesolowski, Phys. Rev. A, 2008, 77, 11444]. This makes it possible to introduce an approximation in the evaluation of the solvent effect on quantum mechanical observables consisting of replacing the ensemble averaged observable by the observable evaluated at ensemble averaged ρ_B(r). This approximation is shown to affect negligibly the solvatochromic shift in the absorption of hydrated acetone. The proposed model provides a continuum type of representation of the solvent, which reflects nevertheless its local structure, and it is to be applied as a post-simulation analysis tool in atomistic level simulations.

---

1 Introduction

Frozen-density embedding theory (FDET) was originally formulated by Wesolowski and Warshel¹ to derive self-consistent expressions for the energy and embedding potential as a universal functional of charge densities: the electrons in the environment (ρ_B(r)), the electron density of the embedded species (ρ_A(r)), and the nuclei in the environment generating the potential (v_B(r)). The FDET approach provides thus an alternative formal framework to the ones used in commonly used QM/MM methods (see classical reviews given in ref. 2–8 or the reviews in volume 115 of Chemical Reviews published in 2015 dedicated entirely to such methods). In QM/MM methods, the embedded wavefunction is obtained using system-specific parameters, whereas the energy and the embedding potential are usually not self-consistent because special terms are a posteriori added. In the original formulation of FDET, a non-interacting reference system was used as a quantum descriptor of the embedded species. Generalizations of FDET to interacting systems have been formulated: for either a multi-determinant reference state⁹ or a one-particle spin-less reduced density matrix used as a quantum mechanical descriptor of the embedded species.¹⁰ In the present work, a non-interacting reference system case is considered. The focus is, however, on the treatment of ρ_B(r) in FDET if the structure of the environment fluctuates. The key feature of the FDET embedding potential is the fact that it is a functional determined uniquely by ρ_B(r), v_B(r), and ρ_A(r):for the definition of the functional E^nad_xct[ρ_A,ρ_B] see the original publications on FDET,^1,9 whereas the approximations used in the present work are given in the corresponding sections of the present work. For the sake of discussion in the present work we will refer to the first two terms in the right-hand-side of eqn (1) as electrostatic components of the embedding potential (v^electr_B(r)).

Usually, the two input quantities ρ_B(r) and v_B(r) needed to evaluate the embedding potential v_emb[ρ_A,ρ_B,v_B](r) are obtained from some approximate method to solve the many-electron problem within the Born–Oppenheimer approximation. Such use of the FDET framework could be called “conventional” as it has been used by many researchers (see our recent review¹¹ or relevant parts of reviews by others,^12–15 for instance). Other authors developed computational methods sharing some elements with FDET.^16–18 One of possible applications of the FDET type of multi-scale method is modeling solvatochromism. In this case, the computational advantages of “conventional” applications originate from the fact that, as in any QM/MM type of modeling, the number of electrons treated explicitly is reduced compared to the full quantum mechanical treatment of the whole system. It is especially important because the embedded wavefunction must be evaluated for many configurations of the environment representing the statistical ensemble. Recently, we proposed an alternative strategy to deal with the solvent effect which concerns the choice made for the input quantities in FDET which are used as the descriptors of the environment: ρ_B(r) and v_B(r).¹⁹ Instead of averaging excitation energies obtained in an ensemble of configurations, the excitation energy is evaluated only once but using the FDET embedding potential corresponding to ensemble averaged charge densities: 〈ρ_B〉(r) and 〈v_B〉(r). Such a strategy hinges, however, on the method to evaluate the average charge densities. In the original work¹⁹ and in the subsequent publications,^20,21 the site-probabilities, i.e., local concentrations of each type of atom of the solvent, dressed up with electrons were used. The site-probabilities were obtained by solving the classical-statistical model of the liquid: 3D-RISM.²² In the present work, a refined model is proposed in which the average charge densities are also generated from site probabilities but the latter are obtained from explicit molecular dynamics simulations at the atomistic level. We aim at the robust computational protocol to evaluate the solvent effect on a specific observable which could be used as a post-simulation analysis of the sample of instantaneous configurations representing the solvent and the solvated species at statistical equilibrium. Such an approach corresponds to the following approximation:

〈ΔÔ[ρ_B]〉_ensemble ≈ ΔÔ[〈ρ_B〉_ensemble], (2)

where Ô represents the operator in question. The significance of the approximation given in eqn (2) is analyzed in the present work using the solvatchromism of acetone in water as a study case. The numerical set-up of the present work is such that it addresses only the approximation given in eqn (2) and not other approximations affecting the calculated excitation energies.

2 FDET based non-uniform continuum model of the solvent

2.1 Frozen-density embedding theory

Frozen-density embedding theory concerns the following constrained search problem in which the number of electrons N_AB, the external potential v_AB(r), and some non-negative function ρ_B(r) are given quantities: If the trial densities in the above search ρ(r) are represented as:where Ψ_A has the form of a N_A-electron function, the condition ∀ rρ(r) ≥ ρ_B(r) is automatically satisfied.

Expressing the total energy as a functional E^EWF_AB[Ψ_A,ρ_B], which depends explicitly on Ψ_A and ρ_B(r), makes it possible to use the Euler–Lagrange equation in order to find stationary many-electron wave functions:

where λ^I is the Lagrange multiplier associated with the normalization. The index I indicates that the present work concerns not only in the lowest-energy solution but also other stationary states.

The total energy functional denoted as E^EWF_AB[Ψ_A,ρ_B] reads:

where the functional E^nad_xct[ρ_A,ρ_B] comprises three components:

E^nad_xct[ρ_A,ρ_B] = T^nad_s[ρ_A,ρ_B] + E^nad_xc[ρ_A,ρ_B] + ΔF^SC[ρ_A]. (7)

The non-additive bi-functionals T^nad_s[ρ_A,ρ_B] and E^nad_xc[ρ_A,ρ_B], and the functional ΔF^SC[ρ_A] are defined by means of the constrained search procedure (see ref. 9).

The necessary condition for Ψ^I_A to satisfy eqn (5) takes the form resembling the conventional Schrödinger equation for many-electrons:

(Ĥ_A + v̂_emb)Ψ^I_A = ε^IΨ^I_A, (8)

where v̂_emb denotes the “embedding operator” given in the form of the potential which in turn is uniquely determined by the charge densities ρ_A(r) and ρ_B(r) and the fixed potential v_B(r). The functional for this potential is obtained as the functional derivative of the corresponding terms in the total energy functional (eqn (6)) and it is given in eqn (1).

The FDET embedding potential comprises the classical electrostatic components (v^electr_B(r)) and the terms which are defined as the functional derivatives of the density functionals: T^nad_s[ρ_A,ρ_B], E^nad_xc[ρ_A,ρ_B], and ΔF^SC[ρ_A]. The latter are, therefore, also the functionals of ρ_A and ρ_B, whereas v^electr_B(r) (the sum of the first two terms in the RHS of eqn (1)) does not depend on ρ_A(r) at all.

Concerning the functional ΔF^SC(WFT)[ρ], its definition depends on what is used as Ĥ_A and on the search domain. It is neglected in practical applications of FDET based simulations. Moreover, in the case of using a non-interacting reference system for embedded species, as in the present work, it does not contribute to the embedding potential at all.

In multi-scale simulation methods based on FDET, the exact density functionals in eqn (7), T^nad_s[ρ_A,ρ_B], E^nad_xc[ρ_A,ρ_B], and ΔF^SC(WFT)[ρ_A], are replaced by their approximate counterparts: T̃^nad_s[ρ_A,ρ_B], Ẽ^nad_xc[ρ_A,ρ_B], and ΔF̃^SC(WFT)[ρ_A].

2.1.1 Local excitations in embedded species from LR-TDDFT. FDET can be used to obtain excitation energies either as other-then-ground state stationary solutions of the Euler–Lagrange equation for the embedded wavefunction (see eqn (8)) or through the Linear-Response Time-Dependent DFT (LR-TDDFT) approach²³ if the embedded species is described by means of a non-interacting reference system.²⁴ In the present work, the second strategy is used. The fact that the embedding potential depends on ρ_A(r) results in an additional contribution (f^emb[ρ_A,ρ_B](r,r′)) to the response kernel for the isolated species (f^A[ρ_A](r,r′)). This kernel can be conveniently represented as:

f^A(emb)[ρ_A,ρ_B](r,r′) = f^A[ρ_A](r,r′) + f^emb[ρ_A,ρ_B](r,r′), (9)

where the contribution to the kernel due to the embedding potential reads:²⁴ Note that the electrostatic components of the FDET embedding potential do not contribute to f^emb[ρ_A,ρ_B](r,r′) because they do not depend on ρ_A(r).

The above combination of FDET embedding potential with the LR-TDDFT approach to the excitation energy represents a physical approximation consisting of neglecting frequency-dependent response of the environment to the oscillating electric field. For that reason we refer to it as Neglect of Dynamic Response of the Environment (NDRE). Extension going beyond NDRE was formulated by Casida and Wesolowski.²⁵ Neugebauer developed a formalism based on this extension treating efficiently the case of chromophores absorbing at similar wavelengths which are localized in different subsystems.¹⁷

2.1.2 Average 〈ρ_B〉(r) from statistical ensembles for structurally flexible environments. So far ρ_B(r) was assumed to have a quantum-mechanical origin and was associated with some molecular system at a given geometry. The FDET formal framework admits also wider choices for ρ_B(r) such as the ensemble averaged ρ_B(r). One can rewrite the FDET equations replacing ρ_B(r) by the statistical ensemble averaged electron density, which is denoted with 〈ρ_B〉(r) throughout this work and with 〈v_B〉(r) being the ensemble averaged potential generated by the nuclei in the environment. The advantage of such an approach is that while a straightforward calculation of the statistically averaged property 〈O[ρ_B]〉 (left-hand-side of eqn (2)) requires M calculations of O_i[ρ_B] for every configuration i = 1…M in the ensemble, the use of the averaged density enables finding 〈O[ρ_B]〉 in one step (the right-hand-side of eqn (2)).

In the work of Kaminski et al.¹⁹ 〈ρ_B〉(r) was obtained by means of solving the 3D-RISM equations to obtain the 3D solute-solvent site distribution function g_γ(r) for each solvent site γ. In the next step, each site probability g_γ(r) was dressed up with spherically averaged atomic densities ρ_γ(r − r′). Specifically, the contribution of each atom is obtained subsequently as the convolution of the site probability with a simple (spherically symmetric) model of electron density of the atom ρ_γ(r − r′), where q_γ is the number of electrons associated with each solvent site γ.

In the method proposed in ref. 19, 〈ρ_B〉(r) obtained from the above procedure is used to evaluate the FDET embedding potential and subsequently the embedded wavefunction and its response needed for LR-TDDFT evaluation of the excitation energy. The solvatochromic shifts were obtained from eqn (2) which takes the following form:

〈Δε[ρ_B]〉 ≈ Δε[〈ρ_B〉]. (12)

Despite several successful applications, the method proposed by Kaminski et al. is not free from drawbacks. In particular, the 3D-RISM method is known to suffer from the numerical instability problem, arising from the non-unique partitioning of the charge density among atoms. Additionally, the difference between the shift obtained in the LHS and the RHS of eqn (12) cannot be attributed only to moving the averaging from shifts to densities but also to the closure used in 3D-RISM which represents an additional approximation.

The purpose of our study was to develop a new method to generate 〈ρ_B〉(r) based not on solving 3D-RISM equations but using explicit numerical simulations at the atomistic level such as classical Molecular Dynamics (MD) simulations.

The procedure to generate 〈ρ_B〉_MD(r) corresponding to the ensemble of instantaneous geometries considered in the previous section consists of the following steps: (a) generation of the solvent–solute distribution function i.e. the map in 3D of concentrations of each type of the atom of the solvent and (b) dressing up with electrons the concentrations obtained in the previous step. Compared to ref. 19 the procedure used in the previous work differs in the first step. The solvatchromic shift is evaluated also using the approximation given in eqn (12).

The nuclear part of the FDET embedding potential depends on neither ρ_B(r) nor ρ_A(r) and averaging is straightforward. The steps needed for the evaluation of the above averaged quantities are described in detail in the subsequent sections.

2.1.3 Generation of a 3D solute–solvent distribution function from the MD trajectory. In the first step, a 3D solute–solvent distribution function g_k(r), which gives the number of particles of the given atom type k of solvent in the given volume element d³r, is computed using the ensemble of coordinates of solvent molecules {R^B_j} produced by the MD simulation. For the MD trajectory with M configurations the first geometry of the solute R^A₁ is considered to be a reference whose geometrical center is located at the point where P is the number of atoms of solute, r_i is a coordinate vector of the solute atom i. Following the algorithm of Kabsch et al.,^26,27 for every geometry R^A_j, j = 2…M we compute a translation vector T_j = O_j − O₁ and a rotation matrix U_j and apply the obtained transformation to the geometry R^A_j

R̃^A_j = U_j(R^A_j − T_j), j = 2…M (13)

to get a new set of transformed geometries {R̃^A_j} which minimizes the Root-Mean-Square-Deviation (RMSD) of atomic positions computed for the geometries R^A₁ and R^A_j as follows: As a result the set of matrices U_j and translation vectors T_j contains all information about the translational, vibrational and rotational degrees of freedom of the solute molecule surrounded with solvent.

In the subsequent step, the set of coordinates of the nuclei of solvent {R^B_j} is averaged over the ensemble of M configurations. To do that, the rotation U_j and translation T_j transformations are applied to the ensemble of solvent geometries {R^B_j} to yield a new set of geometries {R̃^B_j}

R̃^B_j = U_j(R^B_j − T_j), j = 2…M (16)

With the two sets of coordinates {R̃^A_j} and {R̃^B_j} sharing the same origin located at the point O₁, one can finally determine the 3D solvent distribution function g_k(r). For that purpose the whole volume V occupied by the system is divided into small boxes and the number n_α,k of occurrences of the given atom type k into each elementary volume v_α is computed. A uniform set of cubic elementary volumes v_α of (0.3)³ size each was used to construct the 3D solvent distribution function g_k(r). Consequently, the 3D solvent distribution function g_k(α) for the given atom type k in the box α reads:where k is an atom type of solvent, α is the index that runs over every elementary volume v_α and N is the total number of nuclei.

2.1.4 Dressing up with electrons the 3D solvent distribution functions to obtain 〈ρ_B〉(r). To begin with, we first assume that each solvent site γ located at r′ is surrounded by an effective spherically averaged electronic density cloud ρ_γ(r − r′) moving together with that nucleus, the general formula for the statistical ensemble averaged electron density 〈ρ_B〉(r) takes on the following form: Another approximation is that the whole electronic density is concentrated in the same point r′ as that of the nucleus with the total integrated charge being equal to the number of electrons q_γ associated with each site γ.

ρ_γ(r − r′) = q_γδ(r − r′) (19)

that results in the following expression for statistical ensemble averaged densitywhere q_γ is the number of electrons associated with each site γ. The assumption 19 is valid when the ρ_γ(r − r′) is a short-range function confined in the elementary volume v_α that enables us to replace the integral over the space coordinates by the sum over the elementary volumes α so that 20 can be rewritten as followswhere q_k is the number of electrons associated with each type of atoms of solvent k and g_k(α) is the 3D solvent distribution function giving the average number of atoms of the type k in the elementary volume v_α.

2.1.5 Generation of the electrostatic part of the embedding potential. Once the ensemble averaged density 〈ρ_B〉(r) is determined, the ensemble averaged electrostatic potential of the environment is computed aswhere r_α′ denotes a center of the volume v_α, q_k is the number of electrons associated with each type of atom of solvent k, Z_k is a nuclear charge of the atom type k and g_k(α) is a 3D solvent distribution function. It is worthwhile to recall that a similar approximation was introduced in ref. 28 where the averaged electrostatic potential was used. Compared to the method used in the present work, we average not the potential but the charge-densities. Averaging potential is equivalent to averaging the charge density only in the case of embedding potentials which are order-one homogeneous in ρ_B(r) such as the electrostatic potential averaged in ref. 28. The FDET embedding potential (either exact or approximated) is, however, not order-one homogeneous in ρ_B(r) because it contains the non-electrostatic component (the last term in eqn (1)).

2.2 Computational details

In order to validate the approximation given in eqn (2), either the left-hand-side or the right-hand-side of eqn (2) was evaluated using the method introduced in ref. 24 which uses the FDET embedding potential in the general framework of the LR-TDDFT approach to the excited states.²³ The solvatochromic shift for the lowest n → π* transition in an acetone molecule was calculated as the difference of the respective excitation energies in gas and in the solvent Δε = ε_H2O − ε_gas. Subjected to the extensive theoretical^29,30 and experimental study,^31–34 acetone is a very convenient model system to test novel computational approaches for prediction of solvatochromic shifts. To evaluate either sides of eqn (12), the M = 121 instantaneous geometries taken from the QM/MM MD simulation by Adias et al.³⁵ were used. The simulation reported in ref. 35 was performed at 298.15 K for the fixed geometry acetone molecule embedded in 511 explicit water molecules. The acetone molecule was treated at the QM level, while the rest of the system was treated classically. The used instantaneous geometries correspond to a 10 ps long equilibrated trajectory. The excitation energies were evaluated using the ADF code³⁶ featuring the implementation of LR-TDDFT FDET equations by Wesolowski²⁴ subsequently improved by Jacob et al.³⁷ The Slater Type AUG-TZP basis sets were used. The Statistical Averaged Orbital Potential (SAOP) approximation for the exchange–correlation potential in Ĥ^A_o was used in calculations for both isolated and embedded chromophores, the local-density approximation was used for both T^nad_s[ρ_A,ρ_B] and E^nad_xc[ρ_A,ρ_B]. Note that the same approximations (basis sets and functionals) were used in solving LR-TDDFT FDET equations for evaluating both sides of eqn (12). As a result, any differences between the obtained numerical values could be attributed only to the approximation eqn (12) investigated in this work and not to other factors.

For atomistic simulations (the left-hand-side of eqn (2)) the used computational protocol is essentially the same as the one used in the paper of Neugebauer et al.²⁹ applied also for evaluation of absorption band shapes in ref. 38 and 20. Here we report only its key elements and the details specific for the present calculations. The corresponding absorption spectra were computed for M = 121 conformations sampled during the QM/MM MD simulations using the LR-TDDFT FDET method. To make LR-TDDFT FDET calculations affordable we reduce the number of water solvent atoms by introducing a 8 cut off radius whose origin is situated at the center of mass of the acetone. The selected molecular cluster comprises 68 water molecules contributing to ρ_B(r). At each instantaneous geometry, ρ_B(r) is evaluated as a superposition of atomic densities (for the performance of FDET with such construction of ρ_B(r) see ref. 39):

where i refers to the atom in embedding subsystem B, N is the total number of atoms of the solvent and ρⁱ_B(r) is the normalized spherically symmetric electron density of the solvent, Zⁱ_B is the total charge equal to the atomic number of the atom i and nⁱ_B is the net charge of the atom i.

As in ref. 19, the solvent net oxygen and hydrogen charges were taken to be n^O_B = − 0.80e and n^H_B = +0.40e respectively. Those charges correspond to the formal ionic charges. The STO-DZP basis sets were used for the calculation of ρ_B(r). Given the M cluster models, the 8 lowest excitation energies and their corresponding oscillator strengths were computed and averaged over the whole statistical ensemble of M configurations. The number of excitations chosen could be expected to affect the final excitation energies for two reasons. First of all, the n → π* transition is optically forbidden and its small oscillator strength in the solvent is the result of breaking the symmetry in the instantaneous geometry. It is possible therefore that the higher more intense excitations could affect the position of the maximum of the lowest absorption band. In the considered case it is not the case because of the large gap between the two relevant excitation energies. The mean energy of the lowest excitation equals 4.78 eV (the oscillator strength of this transition does not exceed f ≈ 7 × 10⁻⁵). The lowest excitation energy varies between 4.6 eV and 5.0 eV in the considered ensemble. The energy of the second lowest excitation, which has a significantly higher oscillator strength f ≈ 0.01, varies between 6.14 eV and 7.12 eV. The other reason why the number of explicitly considered excitations could affect the calculated solvatochromic shifts is the fact that the quality of higher excitation energies deteriorates due to the wrong asymptotic behavior of the used approximation for the exchange–correlation functional and due to the use of the Davidson diagonalization procedure to solve Casida's equations. The ultimate verification of the importance of the above factors was made by repeating the calculations using different number of excitations (6, 8, or 10). The average excitation energies evaluated for the considered trajectories using any of the chosen number of excitations are practically the same (see Table 1).

Table 1 The ensemble averaged energy of the lowest excitation for the solvated acetone evaluated using different numbers of excitations in solving Casida's equations by means of Davidson's diagonalization procedure

Number of states	6	8	10
Energy, eV	4.7823043	4.7823043	4.7823044

---

All the solvatochromic shifts in the n → π* transition reported in the present work use the reference value of 4.606 eV for this excitation. It was evaluated using the same approximations in LR-TDDFT equations as the ones for the embedded species and the same geometry as the one used in the original MD simulations. After that the electronic spectra of acetone in solution were convoluted with Gaussian functions to mimic the broadening of the electronic transitions due to the vibronic degrees of freedom and the thermal fluctuations of the solvent. For a given transition γ of the frequency ω_γ the molar extinction coefficient , which shows the ability of a chemical substance to absorb light at a given wavelength ω, reads (eqn (24)).

In eqn (24), M is the number of the solute molecular configurations, ω_iγ and f_iγ are the excitation energy and the corresponding oscillator strength of the transition γ computed for the molecular configuration i, D = 3.4819 × 10⁻⁵ mol cm eV⁻¹ is a constant and b is the full width at half maximum (FWHM) chosen to be 0.08 eV. The molar extinction coefficient computed for the first 8 excited state transitions is given by eqn (25)

For each volume element, ρ_B(r) was evaluated as a sum of delta-functions eqn (19) and was computed using eqn (21) that reads as 〈ρ_B〉_MD(r) = (−0.6e)g_H(r) + (−8.8e)g_O(r), where g_H(r) and g_HO(r) are the 3D distribution functions for the oxygen and hydrogen atoms, q_H = −0.6e and q_O = −8.8e being the respective effective net charges.

3 Results

3.1 The absorption band shape from atomistic simulations

Fig. 6 shows the absorption band shape obtained using atomistic simulations. The solvatochromic shifts for the lowest n → π* transition in acetone in comparison with available theoretical and experimental data are collected in Table 2. While the atomistic simulations lead to Δε = 0.18 eV, the experimental shift of the maxima of the absorption band is reported to lie between 0.19–0.22 eV.^31–34

Table 2 Solvatochromic shifts from the water to gas phase of the lowest n → π* transition (in eV)

LR-TDDFT/atomistic	CCSD/MD^30	Exp.^31–34
0.18	0.20	0.19–0.22

---

The recent reported theoretical value of Δε = 0.20 eV³⁰ obtained by means of the combined linear response coupled cluster/molecular mechanics agrees also very well with the result of our atomistic simulations which justifies to use the obtained value and the used approximations (for functionals, basis sets, and the construction of ρ_B(r)) as a reference for analysis of the approximation given in eqn (2).

3.2 Solvatochromic shift from the nonuniform-continuum model: Δε[〈ρ_B〉_MD]

Fig. 1 shows the map of 〈ρ_B〉(r) featuring a low density region around the solute indicating that the solvent density 〈ρ_B〉(r) does not penetrate the interior of the acetone molecule. The radius of the cavity surrounding the acetone nuclei is about 1.5–2. A first solvation shell is present and has a pronounced diffuse character.

Fig. 1 Statistical ensemble averaged embedding density 〈ρ_B〉(r) given in a.u.

Fig. 2 shows the map of the net charge. Being the sum of positively charged nuclei and negatively charged electron density distributions the net density takes on both positive and negative values. Similar to the map of 〈ρ_B〉(r) shown in Fig. 1 and Fig. 2 also features a cavity representing the solute. Two positively charged clouds located in vicinity to the acetone oxygen atom are associated with the protons of molecules of water and indicate the existence of two hydrogen bonds. The statistical ensemble averaged electrostatic potential of the environment v^electr_B[〈ρ_B〉] is shown in Fig. 3. The blue region near the negatively charged ion of the oxygen atom of acetone signifies the field generated by positive charges. On the other hand, negative values of the electrostatic potential prevail in the bottom of the picture that signifies the polarization of the solvent. The acetone's oxygen is situated in a region where the magnitude of the electrostatic potential is close to zero.

Fig. 2 Statistical ensemble averaged net density 〈ρ_Net〉(r) given in a.u.

Fig. 3 Statistical ensemble averaged embedding electrostatic potential of the environment v^electr_B[〈ρ_Net〉_MD](r) given in a.u. The transparent contour represents the Kohn–Sham electronic density of the acetone with isovalue ρ_A = 0.001 a.u.

Fig. 4 shows the embedding potential v_emb(r) generated using the nonuniform-continuum model. To compare the contribution of the non-additive terms E^nad_xct[ρ_A,ρ_B] in eqn (7) and that of the electrostatic term v_B[〈ρ_Net〉_MD] to the embedding potential, the map of the embedding potential was plotted on the same scale as the one used for electrostatic potential (Fig. 5). A zero-line passes between the atoms of carbon and oxygen and slightly changes its location upon the addition of the non-additive terms E^nad_xct[ρ_A,ρ_B]. A comparison of Fig. 5 with Fig. 3 shows that a highly attractive embedding potential, which overlaps the electronic density of the solute, appears in the proximity of the oxygen atom of acetone. The magnitude of embedding potential within the cavity comprising the acetone molecule in Fig. 3 remains unchanged as compared to the Fig. 5 due to the absence of the embedding density near the solute.

Fig. 4 The FDET embedding potential (in a.u.) evaluated at average charge densities of the environment for acetone in water.

Fig. 5 The FDET embedding potential (in a.u.) evaluated at average charge densities of the environment for acetone in water together with the isovalue of ρ_A = 0.001 of ρ_A(r).

The solvatochromic shifts obtained with various methodologies are collected in Table 3.

Table 3 Solvatochromic shifts (in eV) from the water to gas phase for the lowest n → π* transition

Methodology	Solvatochromic shift
a Present work. b Ref. 19.
FDET/non-uniform continuum^a	0.20
FDET/atomistic^a	0.18
FDET/non-uniform continuum^b	0.21
QM/MM^30	0.20
Exp.^31–34	0.19–0.22

---

The solvatochromic shift obtained using the non-uniform continuum model and the average charge densities obtained from explicit atomistic simulations (0.20 eV) are in excellent agreement with both of the theoretical (0.20 eV) and the experimental (0.19–0.22 eV) values. It is also close to that obtained using atomistic simulations described in the previous section (0.18 eV) which justifies the applicability of the approximation given in eqn (2). A reduction of the computational expense (one calculation of the excitation energy instead of more than one hundred) is accompanied by the relative change of the calculated shift by only about 10%. Compared to the shifts obtained using an alternative FDET based non-uniform model using 3D-RISM site probabilities to obtain average charge densities, the model presented in this work also leads to a very similar value. The difference between 0.21 eV (3D-RISM) and 0.18 eV (this work) can be attributed to different force field parameters which are used in both methods and to the closure relation needed to solve 3D-RISM equations.

The n → π* excitation involves the HOMO and LUMO orbitals in the transition as shown in Fig. 7 and 8.

Fig. 6 Simulated absorption spectra of acetone in water, from 4.4 eV to 5.4 eV.

Fig. 7 The highest occupied molecular orbital for the solvated acetone at one of the conformations of the solvent used for averaging.

Fig. 8 The lowest unoccupied orbital for the solvated acetone at one of the conformations of the solvent used for averaging.

3.3 Approaching the statistical limit for 〈ρ_B〉(r)

The quality of the calculated solvatchromic shift in the approach proposed in the present work and discussed in the previous sections depends critically on the representability of the sample of the used configurations. The averaged density considered in the previous section was obtained using a relatively small sample of instantaneous configurations. A straightforward way to approach the statistical limit would consist of using a longer equilibrated MD trajectory. This is not always possible as we propose the current approach as a post-simulation tool to be used independently of the simulations providing the instantaneous configurations. In the present section, we address the issue of representability of the used sample by exploiting two features of the investigated system and FDET: the C_2v symmetry of the chromophore, which results in the same symmetry of 〈ρ_B〉(r) in the statistical limit, and the fact that 〈ρ_B〉(r) (input quantity in FDET equations) is a real function in 3D which can be a subject of manipulations. We notice that the averaged densities and potentials shown in Fig. 1–4 are not symmetric and look very “grainy”. That indicates that the used sample does not represent the statistical limit. In the following, we estimate the effect of this distortion of the symmetry on the calculated shifts. To this end, the asymmetric density 〈ρ_B〉(r) can be symmetrized in a straightforward manner:

〈ρ^symm_B〉(r′) = 〈ρ^symm_B〉(r′′) = 〈ρ^symm_B〉(r′′′) = 〈ρ^symm_B〉(r) (26)

where

〈ρ^symm_B〉(r) = 0.25(〈ρ_B〉(r) + 〈ρ_B〉(r′) + 〈ρ_B〉(r′′) + 〈ρ_B〉(r′′′)) (27)

and where r, r′, r′′, and r′′′, are related by symmetry operations. Compared to their symmetric counterparts the symmetrized charge densities and potentials shown in Fig. 9 and 10, which feature the right symmetry, are smoother because of the smearing out the deviations from the symmetry.

Fig. 9 Symmetrized statistical ensemble averaged embedding density 〈ρ^symm_B〉(r) given in a.u.

Fig. 10 Symmetrized statistical ensemble averaged embedding electrostatic potential of the environment v^electr_B[〈ρ^symm_Net(r)〉_MD](r) given in a.u.

The symmetrization affects negligibly the calculated excitation energy (ε[ρ^symm_B(r)] = 4.8053 eV vs. ε[〈ρ_B(r)〉] = 4.8048 eV). As expected, the oscillator strength of the considered excitation approaches zero (it is actually f[〈ρ^symm_B(r)〉] = 0.2836 × 10⁻¹² due to the used numerical procedures) due to the fact that the lowest n → π* state in acetone (1A2 state in the C_2v symmetry) is optically forbidden. Interestingly, the oscillatory strength in the asymmetric case although small is non-zero f[〈ρ_B(r)〉] = 0.3023 × 10⁻⁶.

4 Discussion and conclusions

Replacing the averaging the quantum mechanical observable evaluated in an ensemble of frozen densities (ρ_B(r)_i) by evaluating the observable at the ensemble average of the charge density of the environment (〈ρ_B〉(r)) is the principal approximation analyzed in this work. This approximation has a compact mathematical form expressed in eqn (2). Physically, it represents replacing the explicit atomistic level of description for the statistical ensemble representing the solvent by two continuous fields: 〈ρ_B〉(r) and the average potential generated by the nuclei 〈v_B〉(r).

In practice, approximation given in eqn (2) leads to huge potential computational savings. With explicit description of the solvent, the evaluation of the observable in question (vertical excitation energy in the present study) must be performed for each configuration representing the statistical ensemble. The use of the continuum fields involves, however, just one evaluation of the quantum mechanical observable. The exact savings depend on the case: costs of generating the statistical ensemble used to generate 〈ρ_B〉(r), the costs of the evaluation of the FDET embedding potential as a functional of charge densities, for instance. In the present study, the CPU timings needed for right-hand-side and left-hand-side of the eqn (2) differ by almost two orders of magnitude.

The idea of using an average embedding potential is not new, to our knowledge, it was introduced by Coutinho et al.,²⁸ where only the electrostatic component of the embedding potential was considered. The present work concerns the complete local embedding potential, which includes also the non-electrostatic term, which as shown in ref. 1 and 9 has the form of a universal functional of charge-densities. Similar to the method introduced in one of our earlier works¹⁹ and opposite to the approach of Coutinho et al., it is not the potential but the charge density of the environment which is averaged. The distinction is not necessary in the case of electrostatic-only embedding because the embedding potential depends linearly on the charge density of the environment. The non-electrostatic term (the functional derivative of E^nad_xct[ρ_A,ρ_B]) can be associated with the quantum confinement effects. In this work, as in other applications of FDET, E^nad_xct[ρ_A,ρ_B] is approximated by an analytic expression Ẽ^nad_xct[ρ_A,ρ_B]. The fact that E^nad_xct[ρ_A,ρ_B] is order-one homogeneous⁴⁰ in neither ρ_A nor ρ_B lies at the origin of the fact that eqn (2) is an approximation and not an exact relation. The same concerns Ẽ^nad_xct[ρ_A,ρ_B]. The present results show that the inhomogeneity in ρ_B of the used approximated density functional for E^nad_xct[ρ_A,ρ_B] results in a negligible error in the solvatochromic shift (in the 0.02 eV range). Concerning the inhomogeneity in ρ_A, it leads to coupling between the electronic structure of the embedded species and the structure of the solvent if the observable is evaluated form the left-hand-side of eqn (2). For each geometry of the solvent, the corresponding ρ_A(r) is different, because the embedding potential depends on ρ_A(r). In the right-hand-side of eqn (2) such couplings are completely neglected.

The actual results presented in this work concern a simple and rigid chromophore solvated in liquid water. In such a case, different configurations representing the statistical ensemble differ but the differences concern mainly the solvent and not the solute. The FDET based averaging approach proposed in the present work can be expected to be adequate used for studies of rigid chromophores in structurally disordered and fluctuating environments. It can be used without further modifications also for studies of a flexible solute provided its structure oscillates around a common global geometry. In the case of a solute occurring in several conformations (cis/trans forms, different tautomeric forms, etc.), the approach proposed in the present work can also be adopted in a straightforward manner. To this end eqn (2) shall be applied for each conformer separately. Such applications are beyond the scope of this article and are subject of our current research.

Using an explicit ρ_A,ρ_B-dependent analytic formula to approximate the functional E^nad_xct[ρ_A,ρ_B] leads also to the corresponding analytic approximated expressions in the embedding potential. Replacing ρ_B by 〈ρ_B〉 can always be made in the analytic formula. In the case of the exact functional, the situation is more involved. Whereas E^nad_xct[ρ_A,〈ρ_B〉] is well defined for any pair of N-representable densities,⁹ the existence of the corresponding potential hinges on a more stringent condition – v-representability of either ρ_A(r) or ρ_A(r) + 〈ρ_B〉(r). Similar to any other density functional theory formulations, this condition cannot be assured in practice. We notice, however, that 〈ρ_B〉(r) is a significantly smoother function that any instantaneous ρ_B(r) taken from the statistical ensemble. It becomes even uniform for non-structured liquids.

The embedding potential (v^FDET_emb[ρ_A,〈ρ_B〉,〈v_B〉](r)) evaluated as a functional of the non-uniform continuum descriptors of the solvent: ρ_A(r), 〈ρ_B〉(r), and 〈v_B〉(r) derived from the statistical ensemble of conformations of the solvent molecules obtained in atomistic simulations was used in the present work to evaluate the solvatochromic shifts in absorption within the general framework of LR-TDDFT. The embedded wavefunction (both interacting or non-interacting) can be, however, used to evaluate any quantum mechanical observable associated with an operator Ô according to the formula given in eqn (2). Although the adequacy of the approximation given in eqn (2) shall be subject to numerical validation similar to that presented in the present work, the presented numerical results indicate that the prospects of using this approximation for other observables are promising.

Acknowledgements

This research was supported by grants from Swiss National Science Foundation (200021_152779). The authors are grateful to Prof. Kurt Mikkelsen for providing the trajectories used for averaging the quantities discussed in this work.

References

1. T. A. Wesolowski and A. Warshel, J. Phys. Chem., 1993, 97, 8050–8053.
2. J. Åquist and A. Warshel, Chem. Rev., 1993, 93, 2523–2544.
3. D. Bakowies and W. Thiel, J. Phys. Chem., 1996, 100, 10580–10594.
4. C. J. Cramer and D. G. Truhlar, Chem. Rev., 1999, 99, 2161–2200.
5. J. Sauer, P. Ugliengo, E. Garrone and V. R. Sounders, Chem. Rev., 1994, 94, 2095.
6. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093.
7. H. M. Senn and W. Thiel, Angew. Chem., Int. Ed., 2009, 48, 1198–1229.
8. J. Gao, Reviews in Computational Chemistry, John Wiley & Sons, Inc., 2007, pp. 119–185.
9. T. A. Wesolowski, Phys. Rev. A: At., Mol., Opt. Phys., 2008, 77, 012504.
10. K. Pernal and T. A. Wesolowski, Int. J. Quantum Chem., 2009, 109, 2520–2525.
11. T. A. Wesolowski, S. Shedge and X. Zhou, Chem. Rev., 2015, 115, 5891–5928.
12. J. Neugebauer, Phys. Rep., 2010, 489, 1–87.
13. C. R. Jacob and J. Neugebauer, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 325–362.
14. A. S. P. Gomes and C. R. Jacob, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem., 2012, 108, 222–277.
15. F. Libisch, C. Huang and E. A. Carter, Acc. Chem. Res., 2014, 47, 2768–2775.
16. N. Govind, Y. A. Wang, A. J. R. da Silva and E. A. Carter, Chem. Phys. Lett., 1998, 295, 129–134.
17. J. Neugebauer, J. Chem. Phys., 2007, 97, 134116.
18. M. Pavanello, J. Chem. Phys., 2013, 138, 204118.
19. J. W. Kaminski, S. Gusarov, T. A. Wesolowski and A. Kovalenko, J. Phys. Chem. A, 2010, 114, 6082–6096.
20. X. Zhou, J. W. Kaminski and T. A. Wesolowski, Phys. Chem. Chem. Phys., 2011, 13, 10565–10576.
21. S. V. Shedge and T. A. Wesolowski, ChemPhysChem, 2014, 15, 3291–3300.
22. A. Kovalenko, Molecular Theory of Solvation, Springer Netherlands, 2003, pp. 169–275.
23. M. E. Casida, Recent Advances in Computational Chemistry, World Scientific, Singapore, 1995, vol. 1, part 1, p. 155.
24. T. A. Wesolowski, J. Am. Chem. Soc., 2004, 126, 11444–11445.
25. M. E. Casida and T. A. Wesolowski, Int. J. Quantum Chem., 2004, 96, 577–588.
26. W. Kabsch, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1976, 32, 922–923.
27. W. Kabsch, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1978, 34, 827–828.
28. K. Coutinho, H. C. Georg, T. L. Fonseca, V. Ludwig and S. Canuto, Chem. Phys. Lett., 2007, 437, 148–152.
29. J. Neugebauer, M. J. Louwerse, E. J. Baerends and T. A. Wesolowski, J. Chem. Phys., 2005, 122, 094115.
30. K. Aidas, J. Kongsted, A. Osted, K. V. Mikkelsen and O. Christiansen, J. Phys. Chem. A, 2005, 109, 8001–8010.
31. I. Renge, J. Phys. Chem. A, 2009, 113, 10678–10686.
32. W. P. Hayes and C. J. Timmons, Spectrochim. Acta, 1965, 21, 529–541.
33. N. S. Bayliss and E. G. McRae, J. Phys. Chem., 1954, 58, 1002–1006.
34. C. W. Porter and C. Iddings, J. Am. Chem. Soc., 1926, 48, 40–44.
35. K. Aidas, K. V. Mikkelsen, B. Mennucci and J. Kongsted, Int. J. Quantum Chem., 2011, 111, 1511–1520.
36. ADF 2014.06, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, http://www.scm.com.
37. C. R. Jacob, J. Neugebauer and L. Visscher, J. Comput. Chem., 2008, 29, 1011–1018.
38. J. Neugebauer, C. R. Jacob, T. A. Wesolowski and E. J. Baerends, J. Phys. Chem. A, 2005, 109, 7805–7814.
39. M. Humbert-Droz, X. Zhou, S. V. Shedge and T. A. Wesolowski, Theor. Chem. Acc., 2014, 133, 1405.
40. A. Zech, F. Aquilante and T. A. Wesolowski, Mol. Phys., 2015 DOI:10.1080/00268976.2015.1125027.

---