Open Access Article

This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Andrey
Laktionov
*,
Emilie
Chemineau-Chalaye
and
Tomasz A.
Wesolowski
*

Université de Genève, Département de Chimie Physique 30, quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland. E-mail: tomasz.wesolowski@unige.ch; Fax: +41-22-379-6518; Tel: +41-22-379-6101

Received
22nd January 2016
, Accepted 3rd March 2016

First published on 3rd March 2016

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) |

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^{11} 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).^{19} 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^{19} 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.^{22} 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) |

(3) |

(4) |

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:

(5) |

The total energy functional denoted as E^{EWF}_{AB}[Ψ_{A},ρ_{B}] reads:

(6) |

E^{nad}_{xct}[ρ_{A},ρ_{B}] = T^{nad}_{s}[ρ_{A},ρ_{B}] + E^{nad}_{xc}[ρ_{A},ρ_{B}] + ΔF^{SC}[ρ_{A}]. | (7) |

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

(Ĥ_{A} + _{emb})Ψ^{I}_{A} = ε^{I}Ψ^{I}_{A}, | (8) |

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: ^{nad}_{s}[ρ_{A},ρ_{B}], Ẽ^{nad}_{xc}[ρ_{A},ρ_{B}], and Δ^{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^{23} if the embedded species is described by means of a non-interacting reference system.^{24} 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:

where the contribution to the kernel due to the embedding potential reads:^{24}

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).

f^{A(emb)}[ρ_{A},ρ_{B}](r,r′) = f^{A}[ρ_{A}](r,r′) + f^{emb}[ρ_{A},ρ_{B}](r,r′), | (9) |

(10) |

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.^{25} Neugebauer developed a formalism based on this extension treating efficiently the case of chromophores absorbing at similar wavelengths which are localized in different subsystems.^{17}

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)).

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.

In the work of Kaminski et al.^{19} 〈ρ_{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 γ.

(11) |

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) |

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^{3}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}_{1} 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_{1} and a rotation matrix U_{j} and apply the obtained transformation to the geometry R^{A}_{j}

to get a new set of transformed geometries {^{A}_{j}} which minimizes the Root-Mean-Square-Deviation (RMSD) of atomic positions computed for the geometries R^{A}_{1} 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.

With the two sets of coordinates {^{A}_{j}} and {^{B}_{j}} sharing the same origin located at the point O_{1}, 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)^{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.

^{A}_{j} = U_{j}(R^{A}_{j} − T_{j}), j = 2…M | (13) |

(14) |

(15) |

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 {^{B}_{j}}

^{B}_{j} = U_{j}(R^{B}_{j} − T_{j}), j = 2…M | (16) |

(17) |

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 γ.

that results in the following expression for statistical ensemble averaged density

where 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 follows

where 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_{α}.

(18) |

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

(20) |

(21) |

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 as

where 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)).

(22) |

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.^{29} 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):

(23) |

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^{−5}). 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).

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)).

(24) |

(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.

The recent reported theoretical value of Δε = 0.20 eV^{30} 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).

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. 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.

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. 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. |

〈ρ^{symm}_{B}〉(r′) = 〈ρ^{symm}_{B}〉(r′′) = 〈ρ^{symm}_{B}〉(r′′′) = 〈ρ^{symm}_{B}〉(r) | (26) |

〈ρ^{symm}_{B}〉(r) = 0.25(〈ρ_{B}〉(r) + 〈ρ_{B}〉(r′) + 〈ρ_{B}〉(r′′) + 〈ρ_{B}〉(r′′′)) | (27) |

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^{−12} 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^{−6}.

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.,^{28} 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^{19} 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^{40} 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,^{9} 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.

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

This journal is © the Owner Societies 2016 |