Polarizable Density Embedding: A Solution to the Electron Spill-Out

Letter Cite This: J. Phys. Chem. Lett. 2017, 8, 5949−5958

pubs.acs.org/JPCL

Polarizable Density Embedding: A Solution to the Electron Spill-Out Problem in Multiscale Modeling Peter Reinholdt,* Jacob Kongsted,* and Jógvan Magnus Haugaard Olsen* Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark S Supporting Information *

ABSTRACT: We analyze the performance of the polarizable density embedding (PDE) modela new multiscale computational approach designed for prediction and rationalization of general molecular properties of large and complex systems. We showcase how the PDE model very effectively handles the use of large and diffuse basis sets that are otherwise questionabledue to electron spill-out effectsin standard embedding models. Based on our analysis, we find the PDE model to be robust and much more systematic than less sophisticated focused embedding models, and thus outline the PDE model as a very efficient and accurate approach to describe the electronic structure of ground and excited states as well as molecular properties of complex, heterogeneous systems.

O

electron density spill-out from the quantum region into the classical part. This issue is especially problematic when using large or diffuse basis sets in the quantum region, and a pragmatic way of avoiding this problem is to use smaller basis sets in this type of embedding calculations. This is an unsatisfactory solution since this prohibits the usual strategy of improving the predicted results by going toward the basis-set limit, which can be especially important for treatment of excited states and molecular properties. To gain a more realistic reproduction of environmental effects, the embedding potential can be improved in several ways compared to the traditional point-charge electrostatic embedding approach. First of all, the quality of the parameters can be improved by using geometry-specific parameters instead of atom-type parameters, which can suffer from poor transferability. This implies a semiclassical fragment-based approach where the parameters are derived from ab initio calculations on the fragments that make up the environment. Continuing in this direction, the electrostatic interactions can be refined through the use of a distributed multipole expansion of the charge distribution of each fragment.5−16 The next logical improvement is to include polarization in the environment such that mutual polarization between the quantum subsystem and its environment can be included.1,9,10,13−31 Finally, it is possible to go beyond electrostatic interactions and include nonelectrostatic effects, such as exchange repulsion, dispersion, and charge transfer.9,10,16,32,33 The importance of each of the effects depends on the molecular system in question as well as the calculated property. An alternative approach to

ne of the great challenges in theoretical chemistry is an accurate description of large molecular systems aimed at addressing specific details of the underlying electronic structure. For large systems, such as chromophores in solution, embedded in proteins, or more generally in inhomogeneous environments, full quantum-mechanical treatments quickly become intractable due to the steep scaling with respect to system size. In addition, such large systems are often very dynamical in nature meaning that inclusion of finite-temperature effects becomes crucial, thus making full electronic-structure treatments out of reach even with the rapidly advancing linear-scaling methods. As such, development of accurate and efficient alternatives to a full quantum-mechanical treatment is of utmost importance. For systems where a natural division into an active and inactive region is possible, such as, e.g., enzymes and fluorescent proteins, quantum mechanics/molecular mechanics (QM/ MM) methods1−4 can be used as good alternatives to full quantum-mechanical calculations. Within the QM/MM approach, the system under study is partitioned into a smaller region treated using quantum mechanics while the rest of the systemthe environmentis described with classical mechanics. The electrostatic potential due to the environmenthere called the embedding potentialis usually derived from atomcentered point charges. The point-charge electrostatic-embedding QM/MM approach is in many cases good as a first approximation of environmental effects. Indeed, it has proven very useful in many aspects, even if the representation of the environment in terms of charges is rather crude. However, this representation is often inadequate, in particular for the calculation of molecular properties, especially of electronically excited states. A fundamental problem relating to QM/MM approaches is the unphysical violation of the Pauli principle for the quantum−classical interface, and hence the possibility of © XXXX American Chemical Society

Received: October 20, 2017 Accepted: November 27, 2017 Published: November 27, 2017 5949

DOI: 10.1021/acs.jpclett.7b02788 J. Phys. Chem. Lett. 2017, 8, 5949−5958

Letter

The Journal of Physical Chemistry Letters

Figure 1. Difference in the electron density of p-nitrophenolate (pNP) embedded in a water environment including a single sodium ion. Densities were calculated using the CAM-B3LYP exchange-correlation functional and the 6-31G* or 6-31+G* basis sets. Isodensity surfaces were drawn at isovalues of Δρ = ± 0.001 electrons·bohr−3, where the positive isovalue is shown in red and the negative isovalue is shown in blue.

essential to include exchange repulsion when including highaccuracy electrostatics.64−66 In the following, we will illustrate the superior performance of the PDE model and especially that this model leads to a systematic and robust embedding scheme. For this, we will compare its performance to the polarizable embedding (PE) model.14,15 Common for PE and PDE is that environment polarization is very efficiently obtained via distributed anisotropic dipole−dipole polarizabilities, but PE does not model exchange repulsion and uses a distributed multipole expansion to model the electrostatics. Finally, we will also compare the performance of PDE to the simpler pointcharge electrostatic-embedding QM/MM. A further aspect related to the PDE model (and polarizable models in general) is the ability to account for local-field effects associated with an externally applied electric field. When a system is subjected to an external electric field, the local field experienced by the quantum region is modified by the environment polarization. Depending on the local structure of the environment, this polarization might lead to either

the (semi)classical embedding outlined above is to preserve the quantum nature of the environment, i.e., use quantum embedding.34−59 An overview of embedding models can be found in refs 60 and 61. Polarizable density embedding (PDE) is a fragment-based multiscale embedding model that includes accurate electrostatics modeled in terms of electronic densities, polarization based on the induced-dipole approach, and exchange repulsion through a projection operator.62 All quantities are derived from fragment-based ab initio calculations. The combination of density-based electrostatics with classical polarization makes PDE both accurate and efficient. The operator modeling the exchange repulsion is based on the Huzinaga−Cantu approach.63 It scales as the square of the interfragment overlap matrix, which is derived from interfragment overlap in basisfunction space, effectively creating a penalty function that leads to a confinement of the quantum-region electron density. This limits the unphysical violation of the Pauli principle and problems associated with electron spill-out effects. In fact, it is 5950


Letter


Figure 2. Convergence characteristics of the PE and PDE models with respect to basis-set size. The CAM-B3LYP exchange-correlation functional was used. Panel (a) shows the lowest excitation energy, while (b) shows the oscillator strength and (c) shows the 2PA cross section for the lowest transition.

enhancement or reduction of the external field. This direct polarization of the environment by the external field is accounted for in the effective external field (EEF) approach,67 which we combine with the PDE model in this work. The direct environment polarization only affects the strength of the electric field at the quantum region and leaves the frequency unaltered; excitation energies are unaffected by local-field effects and only response and transition properties such as oscillator strengths and two-photon absorption (2PA) strengths are influenced. To illustrate the superior performance of the PDE model with respect to the electron spill-out effect, we compare the electronic ground-state density of p-nitrophenolate (pNP) where an aqueous environment is modeled by the PDE, PE or TIP3P embedding potential (see Figure 1). Red surfaces in Figure 1 correspond to regions where the PE and TIP3P models predict an excess electron density (0.001 electrons· bohr−3) compared to the PDE model, and conversely for blue surfaces. From Figure 1a,c it is seen that rather small density differences appear when using small and compact basis sets, such as the 6-31G* basis, thus indicating similar performance regarding their abilities to predict the density of the solute. However, when including diffuse functions in the basis set, i.e., the 6-31+G* basis set, much larger density differences appear as seen in Figure 1b,d. In particular, a relatively large portion of the electron density is found close to the sodium cation. A similar picture would emerge for other QM/MM-type embedding models as exemplified here using TIP3P. Such electron-density leaks can cause an artificial stabilization of the solute charge density as we will demonstrate in the following analyses. The degree of electron spill-out depends on the properties of the quantum-region molecule as well as its environment. Generally, it is expected that electron spill-out is worst for anionic molecules, least for cationic molecules, and somewhere in between for neutral molecules. To investigate this, we computed ground-state densities of a series of anionic, neutral, and cationic solutes, namely permanganate, acrolein, zwitterionic glycine, pyridinium, and p-nitroanilinium. The

density-difference plots are shown in the Supporting Information. The expected behavior is seen for permanganate, acrolein, and pyridinium, but as exemplified by p-nitroanilinium and glycine, even in neutral and cationic systems, rather large density differences can occur. This is likely related to the charge-separated character of these molecules. The difference in the PE and PDE models ability to predict ground-state densities is mirrored in the quality of their prediction of spectroscopic properties. This can be seen in Figure 2, which features the excitation energy, oscillator strength, and 2PA cross section of the lowest π → π* transition in pNP, calculated from a single structure randomly chosen from a molecular dynamics (MD) trajectory. This excitation is dominated by an electronic transition from the highest-occupied molecular orbital (HOMO) to the lowestunoccupied molecular orbital (LUMO). It is often characterized as an intramolecular charge-transfer (ICT) transition because an electron is transferred from the negatively charged phenolate to the nitro group. The results are based on calculations where the pNP molecule is in the quantum region, and the environment consists of 390 water molecules and one sodium ion. Increasing the number of basis functions leads to a lowering of the excitation energy and enhanced oscillator strength and 2PA cross section (e.g., when going from ccpVDZ to cc-pVTZ). As expected the two models give virtually identical results when a basis set without diffuse functions is used (i.e., 6-31G*, cc-pVDZ, and cc-pVTZ). However, when adding diffuse functions (i.e., 6-31+G*, aug-cc-pVDZ, and augcc-pVTZ), the models yield different results. This is seen most clearly in the excitation energy and 2PA cross section where there is a more substantial decrease and increase, respectively, when the PE model is used. The difference between the excitation energies predicted by the PE model compared to the PDE model can be explained by the differential stabilization of the ground and excited state. First of all, when going from a nondiffuse to a diffuse basis set, the electron density can expand and thus interact more strongly with the environment. Furthermore, because the excited-state density is more diffuse 5951


Letter


Figure 3. Convergence characteristics of the PE and PDE models with respect to basis-set size. The CAM-B3LYP exchange-correlation functional was used. The panels show the excitation energy of the lowest intense ligand-to-metal charge-transfer (LMCT) transition in permanganate, the lowest n → π* and π → π* transitions in acrolein, the lowest n → π* transition in zwitterionic glycine, the lowest π → π* transition in pyridinium, and the lowest π → π* transition in p-nitroanilinium. We note that the 6-31+G* basis set has not been developed for manganese, which is the reason for the missing values for permanganate.

transitions in acrolein, the lowest n → π* transition in zwitterionic glycine, the lowest π → π* transition in pyridinium, and the lowest π → π* transition in pnitroanilinium. These calculations were performed on single structures in aqueous environments taken from MD trajectories. Here we see that the excitation energy of the LMCT transition in permanganate is hardly affected when adding diffuse functions. This is not because there is no electron spill-out but due to the fact that spill-out is already present without diffuse functions in the basis set, which can also be seen in the density-difference plots (see Figure S1 in the SI). From the trends of the excitation energy of the π → π* transition in acrolein and p-nitroanilinium, we see that electron spill-out can cause problems for neutral molecules and even for cationic molecules. The severity of such problems is, perhaps unsurprisingly, very system dependent, and the same type of transition in pyridinium is hardly affected by neither basis set nor embedding model. The excitation energy of the n → π* transition in acrolein and glycine displays the opposite trend compared to π → π* transitions with respect to the addition of diffuse functions in the basis set and with respect to the embedding treatment. The reason for this behavior is that there

for this type of transition, it will interact more strongly with the environment than the ground-state electron density. This explains the general decrease of the excitation energy predicted by both models when diffuse functions are added to the basis set. The different behavior of the two embedding models can finally be explained by the fact that the electron densities predicted by the PE model can expand further into the environment because of the lack of exchange repulsion. This can therefore lead to artificially stabilized densities resulting in a much larger energy decrease compared to the PDE model where exchange repulsion is modeled. It is thus clear that the PDE model can predict more realistic electron densities, which affects both response and transition properties as demonstrated in Figure 2. In order to clearly demonstrate the spill-out effect and the capabilities of the PDE model to prevent it, the analysis presented above was performed using the pNP solute, i.e., an anion where electron spill-out is more likely to occur. Figure 3 shows the basis-set dependence of the excitation energy of selected transitions in the before-mentioned series of solutes, i.e., the lowest ligand-to-metal charge-transfer (LMCT) transition in permanganate, the lowest n → π* and π → π* 5952


Letter


Figure 4. Convergence of the lowest excitation energy (a), its associated oscillator strength (b), and 2PA cross section (c) with respect to the size of the quantum region which is increased by gradually including the nearest molecules (or ion) according to their distance from pNP. The differences to the full cluster results are plotted for 10 different structures which consist of pNP and 40 water molecules or 39 water molecules and one sodium ion. The average error in a property (bold line), as well as minimum and maximum errors (thin lines), are shown.

poses the most challenging case in terms of electron spill-out problems. In Figure 4 we examine the convergence of the spectroscopic properties of pNP with respect to the size of the quantum region. In this and the following analyses we used the 6-31+G* basis set. The convergence analysis is performed using ten randomly chosen structures from a MD trajectory. The structures consist of the pNP ion embedded inside a cluster with either 39 water molecules and one sodium ion or 40 water molecules. The results are plotted relative to the values obtained from full quantum-mechanical calculations on the cluster. The average excitation energy obtained from the full cluster calculations is 3.54 eV, while the average oscillator strength is 0.49 and the average 2PA cross section is 13.6 GM. The PDE model consistently overestimates the excitation energy (Figure 4a) by about 0.06 ± 0.01 eV with only pNP in the quantum region. This is expected and can be explained in terms of intermolecular interactions. PDE models electrostatic, induction, and exchange-repulsion interactions between the quantum fragment and the fragments in its environment. The electrostatic interactions can be overall attractive or repulsive, which, among other things, depends on the orientation of the

are two environment-induced stabilization effects that come into play for this type of transition. As for π → π* transitions, we expect a more diffuse excited-state density, which thus has a stabilizing effect. However, for n → π* transitions, electron density is moved from an oxygen lone pair, where it had a considerable stabilizing effect in the ground-state, thus, counteracting the stabilizing effect of the more diffuse excited-state density. Overall, this results in an excited-state density that is not stabilized to the same degree as for a π → π* transition and the excitation energy is therefore mainly determined by the stabilization of the ground-state density. Adding diffuse functions to the basis set will therefore stabilize the ground-state density to a higher degree than the excitedstate density, which explains the increase of the excitation energy when the basis set is extended. Furthermore, since the PE model will stabilize the ground-state density to a higher degree than the PDE model this leads a higher excitation energy for the PE model. Clearly, the degree of electron spillout is not only dictated by the charged state of the solute, but also by other factors such as the character of the transition and the specific interactions between the solute and solvent. The following analyses are based solely on the pNP solute because it 5953


Letter


Figure 5. Distributions of the five lowest excitation energies based on 120 structures from a 1.2 ns MD trajectory when including no solvent in the quantum region are shown for (a) PE and (b) PDE. The error distributions relative to a reference in which 10 solvent molecules are treated fully quantum mechanically are shown in panels (c) and (d).

Figure 6. Distributions of the oscillator strengths across the trajectory are shown in (a) and (b), while the error distributions relative to reference calculations are shown in (c) and (d).

attractive, e.g., dispersion and charge transfer, and lead to a lower excitation energy. Therefore, it stands to reason that the inclusion of electrostatic, induction, and exchange repulsion, such as in the PDE model, should either give the correct excitation energy (if the remaining interactions are negligible) or higher, which is exactly what is observed. Furthermore, since dispersion is not modeled by the density functional approximation used here, it seems that charge-transfer interactions play a role and are needed to get a better agreement with the full cluster results. It should be noted, however, that conventional functionals are known to over-

interacting fragments. On average, we would expect that the electrostatic interactions would lead to a higher excitation energy for the ICT transition, because the ground state, which has a larger dipole moment, will have more favorable interactions with the polar water molecules and the sodium ion (when present). The induction interactions are always attractive and will therefore stabilize the ground state and, to a higher degree, the excited state, because of its more diffuse character, thus lowering the excitation energy. The opposite holds for the exchange-repulsion interactions. The other interactions that might be relevant in this case would be 5954


Letter


Figure 7. Distributions of the 2PA cross section are shown in (a) and (b). The error distributions relative to reference calculations are shown in (c) and (d).

estimate charge-transfer interactions in certain systems,68−71 which could affect our reference values and possibly explain the somewhat slow convergence that we observe for the PDE model. In any case, we observe smooth convergence with respect to the size of the quantum region. For example, with 10 molecules added to the quantum region the errors drop to approximately 0.03 ± 0.01 eV, and to around 0.02 ± 0.01 eV with 15 molecules in the quantum region (in addition to pNP). Importantly, the span of the errors does not change significantly for the smaller quantum regions, and all errors are positive, which shows that the PDE model is both stable and consistent. This is further demonstrated in Figure 5 where the distributions of the five lowest excitation energies are shown (Figure 5b) together with the error distributions relative to calculations with ten molecules in the quantum region (Figure 5d). These results are from calculations performed on 120 structures and include all water molecules within 12 Å of the pNP ion and one sodium if it is within the cutoff. Similar distributions of the oscillator strength and 2PA cross section of the lowest transition are presented in Figures 6 and 7, respectively. From the error distributions of the excitation energies (Figure 5d) it is seen that the PDE model yields excitation energies with average errors around 0.05 eV or less compared to the larger quantum region, and is thus equally applicable to the higher-lying excitations. The PE model, which on average reproduces the full cluster calculation quite well, yields a much larger error span for small quantum regions, where the errors range from about −0.03 to 0.05 eV (see Figure 4). From Figure 5c we see that this is even worse for higher-lying excitations. The average errors are small for the first and second transition but increase for the higherlying excitations, where also the error distributions are rather broad. Going back to Figure 4, we see that the convergence behavior of the PE model is rather irregular. On average, no significant change is visible until about five molecules are included in the quantum region, after which the energy tends to increase slightly. As discussed above, the reason for this

behavior is the differential stabilization of the ground- and excited-state electron densities mainly due to the induction interactions and the lack of exchange repulsion, which leads to electron spill-out and thus exacerbates the stabilization effects. In other words, there is no proper confinement of the quantum-region electron density by the surrounding molecules in the calculations with a small quantum region, which is present in the larger quantum regions where pNP is surrounded by quantum molecules. Incidentally, this limitation in the PE model leads to excitation energies that correspond quite well to the full cluster calculations even for small quantum regions. However, if other interactions were to be included, but not exchange repulsion, then this would lead to an underestimation of the excitation energy (cf. discussion above regarding intermolecular interactions). Included in Figure 4 are also calculations based on an embedding potential that consists of TIP3P point charges, with which only electrostatic interactions are modeled, though it may be argued that some polarization is implicitly included via excessive charges. The excitation energies calculated using the TIP3P embedding potential deviate from the full cluster results by 0.10 eV on average for the smallest quantum region. This result is consistent with the explanation in terms of intermolecular interactions given above. The span of the errors covers a range from 0.03 to 0.13 eV, which is most likely due to a varying degree of electron spill-out. Interestingly, the convergence behavior with respect to the size of the quantum region is quite smooth, which further indicates that the main reason for the low average errors produced by the PE model is the use of polarizable embedding potentials without exchange repulsion. The convergence of the oscillator strength with respect to the size of the quantum region is shown in Figure 4b. The PDE model without including the EEF overestimates the oscillator strength, with errors that initially increase as the quantum region is expanded (from about 0.04 to 0.08), before converging rather slowly toward the full cluster results. The 5955


Letter


in order to systematically improve the quality of results. Getting reliable predictions of transition strengths is significantly improved by the inclusion of local field effects. Indeed, we have demonstrated that accurate results for oscillator strengths and 2PA cross sections can be obtained with the PDE model in combination with EEF. Thus, the use of effective and focused embedding models in quantum chemistry is a powerful and cost-effective approach to describing the electronic structure of complex, heterogeneous systems.

PE model has essentially the same errors and convergence behavior except for the smallest quantum regions, where there is a small insignificant difference. Including the EEF has a dramatic effect on the oscillator strength, yielding errors that are significantly smaller for all quantum regions. In fact, even when using the smallest quantum region, the oscillator strengths calculated using the PDE model with the EEF only deviate between −0.02 and 0.02 from the full cluster results with an average error just below 0.02. The PE model including EEF produces slightly larger errors, but already with about five additional molecules in the quantum region, both models have essentially the same performance. Overall, the average errors do not change significantly as the quantum region is expanded, which is also confirmed by the error distributions shown in Figure 6. There is, however, a slight tendency for the errors to increase for the smallest quantum regions (as seen in Figure 4b), which could explain why the error distributions, shown in Figure 6c,d, are centered slightly to the negative side, i.e., the oscillator strengths are slightly larger when 10 molecules are added in the quantum region. The 2PA cross section is a challenging and interesting property to investigate in the context of multiscale embedding models because it is highly sensitive to the environment; it relies on an accurate description of not only the directly involved molecular orbitals but, in principle, all molecular orbitals. The molecular orbitals that contribute most to the 2PA cross section of the lowest transition in the solvated pNP molecule are most likely those that are localized on the pNP itself. It is evident from the analysis of the oscillator strength that the EEF is crucial to get good agreement with the full cluster result. Indeed, without the EEF the trends for the 2PA cross section are essentially the same as for the oscillator strength. On that account, we will focus on the results that include the EEF. In both cases, we again see a substantial decrease in the errors, showing that the EEF is indeed also very important for the 2PA cross section. However, where both models have an essentially equivalent performance for the oscillator strength, the situation is different for the 2PA cross section, as can be seen in Figure 4c. The average error of PDE(EEF) is around zero GM for all quantum regions. PE(EEF) initially has an average error of around one GM, which then gradually decreases and the models become nearly equivalent with around ten additional molecules in the quantum region. Both models have an error range of approximately ±1 GM for the smaller quantum region. The improved performance of the PDE model is evident in Figure 7 where the effect of the EEF is also apparent. Here we see that the error distribution (Figure 7d) is centered close to zero and that most errors are below one GM. This result is indeed a significant improvement compared to the PE model (Figure 7c), which has a rather broad error distribution centered at approximately one GM. This indicates that the molecular orbitals that are important for the 2PA cross section are better reproduced by the PDE model. The presented analyses demonstrate the importance of including effects beyond electrostatics and polarization when describing transition properties of solute−solvent systems using focused multiscale embedding models. In particular, we have shown that the PDE model, which includes exchange repulsion, gives a robust and systematic prediction of excitation energies without relying on cancellation of errors. This is imperative if one is to apply the usual strategy of increasing the basis-set size

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b02788. Electron density difference plots of a series of solutes in aqueous solvent, together with computational details (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected]. ORCID

Peter Reinholdt: 0000-0003-2406-700X Jacob Kongsted: 0000-0002-7725-2164 Jógvan Magnus Haugaard Olsen: 0000-0001-7487-944X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Erik Rosendahl Kjellgren for providing the structures of pNP in water solution. Computation/simulation for the work described in this paper was supported by the DeIC National HPC Centre, SDU. J.K. thanks the Danish Council for Independent Research (Grant ID DFF-0602-02122B) for financial support. J.M.H.O. thanks the Carlsberg Foundation (Grant ID CF15-0823) for financial support.

■

REFERENCES

(1) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227−249. (2) Lin, H.; Truhlar, D. G. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185−199. (3) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (4) Brunk, E.; Rothlisberger, U. Mixed Quantum Mechanical/ Molecular Mechanical Molecular Dynamics Simulations of Biological Systems in Ground and Electronically Excited States. Chem. Rev. 2015, 115, 6217−6263. (5) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. II. Overlap Populations, Bond Orders, and Covalent Bond Energies. J. Chem. Phys. 1955, 23, 1841−1846. (6) Karlström, G. In Proceedings of the 5th Seminar on Computational Methods in Quantum Chemistry; van Duijnen, P. T., Nieuwpoort, W. C., Eds.; Max-Planck-Institut für Physik und Astrophysik: Groningen, The Netherlands, 1981; p 353. (7) Stone, A. J. Distributed multipole analysis, or how to describe a molecular charge distribution. Chem. Phys. Lett. 1981, 83, 233−239. (8) Stone, A. J.; Alderton, M. Distributed multipole analysis. Mol. Phys. 1985, 56, 1047−1064.

5956


Letter

The Journal of Physical Chemistry Letters (9) Day, P. N.; Jensen, J. H.; Gordon, M. S.; Webb, S. P.; Stevens, W. J.; Krauss, M.; Garmer, D.; Basch, H.; Cohen, D. An effective fragment method for modeling solvent effects in quantum mechanical calculations. J. Chem. Phys. 1996, 105, 1968−1986. (10) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. The Effective Fragment Potential Method: A QM-Based MM Approach to Modeling Environmental Effects in Chemistry. J. Phys. Chem. A 2001, 105, 293−307. (11) Gagliardi, L.; Lindh, R.; Karlström, G. Local properties of quantum chemical systems: The LoProp approach. J. Chem. Phys. 2004, 121, 4494−4500. (12) Stone, A. J. Distributed Multipole Analysis: Stability for Large Basis Sets. J. Chem. Theory Comput. 2005, 1, 1128−1132. (13) Söderhjelm, P.; Husberg, C.; Strambi, A.; Olivucci, M.; Ryde, U. Protein Influence on Electronic Spectra Modeled by Multipoles and Polarizabilities. J. Chem. Theory Comput. 2009, 5, 649−658. (14) Olsen, J. M.; Aidas, K.; Kongsted, J. Excited States in Solution through Polarizable Embedding. J. Chem. Theory Comput. 2010, 6, 3721−3734. (15) Olsen, J. M. H.; Kongsted, J. Advances in Quantum Chemistry; Elsevier: Amsterdam, 2011; pp 107−143. (16) Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Accurate First Principles Model Potentials for Intermolecular Interactions. Annu. Rev. Phys. Chem. 2013, 64, 553−578. (17) Singh, U. C.; Kollman, P. A. A combined ab initio quantum mechanical and molecular mechanical method for carrying out simulations on complex molecular systems: Applications to the CH3Cl + Cl− exchange reaction and gas phase protonation of polyethers. J. Comput. Chem. 1986, 7, 718−730. (18) Thompson, M. A.; Schenter, G. K. Excited States of the Bacteriochlorophyll b Dimer of Rhodopseudomonas viridis: A QM/ MM Study of the Photosynthetic Reaction Center That Includes MM Polarization. J. Phys. Chem. 1995, 99, 6374−6386. (19) Thompson, M. A. QM/MMpol: A Consistent Model for Solute/Solvent Polarization. Application to the Aqueous Solvation and Spectroscopy of Formaldehyde, Acetaldehyde, and Acetone. J. Phys. Chem. 1996, 100, 14492−14507. (20) Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580−10594. (21) Bryce, R. A.; Buesnel, R.; Hillier, I. H.; Burton, N. A. A solvation model using a hybrid quantum mechanical/molecular mechanical potential with fluctuating solvent charges. Chem. Phys. Lett. 1997, 279, 367−371. (22) Jensen, L.; van Duijnen, P. T.; Snijders, J. G. A discrete solvent reaction field model within density functional theory. J. Chem. Phys. 2003, 118, 514−521. (23) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Coupled Cluster/Molecular Mechanics Method: Implementation and Application to Liquid Water. J. Phys. Chem. A 2003, 107, 2578−2588. (24) Illingworth, C. J. R.; Gooding, S. R.; Winn, P. J.; Jones, G. A.; Ferenczy, G. G.; Reynolds, C. A. Classical Polarization in Hybrid QM/ MM Methods. J. Phys. Chem. A 2006, 110, 6487−6497. (25) Nielsen, C. B.; Christiansen, O.; Mikkelsen, K. V.; Kongsted, J. Density functional self-consistent quantum mechanics/molecular mechanics theory for linear and nonlinear molecular properties: Applications to solvated water and formaldehyde. J. Chem. Phys. 2007, 126, 154112. (26) Geerke, D. P.; Thiel, S.; Thiel, W.; van Gunsteren, W. F. Combined QM/MM Molecular Dynamics Study on a CondensedPhase SN2 Reaction at Nitrogen: The Effect of Explicitly Including Solvent Polarization. J. Chem. Theory Comput. 2007, 3, 1499−1509. (27) Curutchet, C.; Muñoz-Losa, A.; Monti, S.; Kongsted, J.; Scholes, G. D.; Mennucci, B. Electronic Energy Transfer in Condensed Phase Studied by a Polarizable QM/MM Model. J. Chem. Theory Comput. 2009, 5, 1838−1848. (28) Lipparini, F.; Barone, V. Polarizable Force Fields and Polarizable Continuum Model: A Fluctuating Charges/PCM Approach. 1. Theory and Implementation. J. Chem. Theory Comput. 2011, 7, 3711−3724.

(29) Lipparini, F.; Cappelli, C.; Barone, V. Linear Response Theory and Electronic Transition Energies for a Fully Polarizable QM/ Classical Hamiltonian. J. Chem. Theory Comput. 2012, 8, 4153−4165. (30) Boulanger, E.; Thiel, W. Solvent Boundary Potentials for Hybrid QM/MM Computations Using Classical Drude Oscillators: A Fully Polarizable Model. J. Chem. Theory Comput. 2012, 8, 4527−4538. (31) Mennucci, B. Modeling environment effects on spectroscopies through QM/classical models. Phys. Chem. Chem. Phys. 2013, 15, 6583. (32) Ö hrn, A.; Karlström, G. A theoretical study of the solvent shift to the transition in formaldehyde with an effective discrete quantum chemical solvent model including non-electrostatic perturbation. Mol. Phys. 2006, 104, 3087−3099. (33) Giovannini, T.; Lafiosca, P.; Cappelli, C. A General Route to Include Pauli Repulsion and Quantum Dispersion Effects in QM/MM Approaches. J. Chem. Theory Comput. 2017, 13, 4854−4870. (34) Inglesfield, J. E. A method of embedding. J. Phys. C: Solid State Phys. 1981, 14, 3795−3806. (35) Barandiarán, Z.; Seijo, L. The ab initio model potential representation of the crystalline environment. Theoretical study of the local distortion on NaCl:Cu+. J. Chem. Phys. 1988, 89, 5739−5746. (36) Cortona, P. Self-consistently determined properties of solids without band-structure calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 44, 8454−8458. (37) Wesołowski, T. A.; Warshel, A. Frozen density functional approach for ab initio calculations of solvated molecules. J. Phys. Chem. 1993, 97, 8050−8053. (38) Wesołowski, T. A.; Weber, J. Kohn-Sham equations with constrained electron density: an iterative evaluation of the groundstate electron density of interacting molecules. Chem. Phys. Lett. 1996, 248, 71−76. (39) Kotliar, G.; Savrasov, S. Y.; Haule, K.; Oudovenko, V. S.; Parcollet, O.; Marianetti, C. A. Electronic structure calculations with dynamical mean-field theory. Rev. Mod. Phys. 2006, 78, 865−951. (40) Swerts, B.; Chibotaru, L. F.; Lindh, R.; Seijo, L.; Barandiaran, Z.; Clima, S.; Pierloot, K.; Hendrickx, M. F. A. Embedding Fragment ab Initio Model Potentials in CASSCF/CASPT2 Calculations of Doped Solids: Implementation and Applications. J. Chem. Theory Comput. 2008, 4, 586−594. (41) Zgid, D.; Chan, G. K.-L. Dynamical mean-field theory from a quantum chemical perspective. J. Chem. Phys. 2011, 134, 094115. (42) Höfener, S.; Gomes, A. S. P.; Visscher, L. Molecular properties via a subsystem density functional theory formulation: A common framework for electronic embedding. J. Chem. Phys. 2012, 136, 044104. (43) Höfener, S.; Visscher, L. Calculation of electronic excitations using wave-function in wave-function frozen-density embedding. J. Chem. Phys. 2012, 137, 204120. (44) Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller, T. F. A Simple, Exact Density-Functional-Theory Embedding Scheme. J. Chem. Theory Comput. 2012, 8, 2564−2568. (45) Knizia, G.; Chan, G. K.-L. Density Matrix Embedding: A Simple Alternative to Dynamical Mean-Field Theory. Phys. Rev. Lett. 2012, 109, 186404. (46) Höfener, S.; Gomes, A. S. P.; Visscher, L. Solvatochromic shifts from coupled-cluster theory embedded in density functional theory. J. Chem. Phys. 2013, 139, 104106. (47) Knizia, G.; Chan, G. K.-L. Density Matrix Embedding: A StrongCoupling Quantum Embedding Theory. J. Chem. Theory Comput. 2013, 9, 1428−1432. (48) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 325−362. (49) Libisch, F.; Huang, C.; Carter, E. A. Embedded Correlated Wavefunction Schemes: Theory and Applications. Acc. Chem. Res. 2014, 47, 2768−2775. (50) Myhre, R. H.; Sánchez de Merás, A. M. J.; Koch, H. Multi-level coupled cluster theory. J. Chem. Phys. 2014, 141, 224105. (51) Inglesfield, J. E. The Embedding Method for Electronic Structure; IOP Publishing: Bristol, U.K., 2015; pp 2053−2563. 5957


Letter

The Journal of Physical Chemistry Letters (52) Fornace, M. E.; Lee, J.; Miyamoto, K.; Manby, F. R.; Miller, T. F. Embedded Mean-Field Theory. J. Chem. Theory Comput. 2015, 11, 568−580. (53) Chulhai, D. V.; Jensen, L. Frozen Density Embedding with External Orthogonality in Delocalized Covalent Systems. J. Chem. Theory Comput. 2015, 11, 3080−3088. (54) Wesołowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891−5928. (55) Chibani, W.; Ren, X.; Scheffler, M.; Rinke, P. Self-consistent Green’s function embedding for advanced electronic structure methods based on a dynamical mean-field concept. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93, 165106. (56) Wouters, S.; Jiménez-Hoyos, C. A.; Sun, Q.; Chan, G. K.-L. A Practical Guide to Density Matrix Embedding Theory in Quantum Chemistry. J. Chem. Theory Comput. 2016, 12, 2706−2719. (57) Pernal, K. Reduced density matrix embedding. General formalism and inter-domain correlation functional. Phys. Chem. Chem. Phys. 2016, 18, 21111−21121. (58) Ding, F.; Tsuchiya, T.; Manby, F. R.; Miller, T. F. LinearResponse Time-Dependent Embedded Mean-Field Theory. J. Chem. Theory Comput. 2017, 13, 4216−4227. (59) Sæther, S.; Kjærgaard, T.; Koch, H.; Høyvik, I.-M. DensityBased Multilevel Hartree-Fock Model. J. Chem. Theory Comput. 2017, 13, 5282. (60) Gomes, A. S. P.; Jacob, C. R. Quantum-chemical embedding methods for treating local electronic excitations in complex chemical systems. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2012, 108, 222. (61) Sun, Q.; Chan, G. K.-L. Quantum Embedding Theories. Acc. Chem. Res. 2016, 49, 2705−2712. (62) Olsen, J. M. H.; Steinmann, C.; Ruud, K.; Kongsted, J. Polarizable Density Embedding: A New QM/QM/MM-Based Computational Strategy. J. Phys. Chem. A 2015, 119, 5344−5355. (63) Huzinaga, S.; Cantu, A. A. Theory of Separability of ManyElectron Systems. J. Chem. Phys. 1971, 55, 5543−5549. (64) Fradelos, G.; Wesołowski, T. A. Importance of the Intermolecular Pauli Repulsion in Embedding Calculations for Molecular Properties: The Case of Excitation Energies for a Chromophore in Hydrogen-Bonded Environments. J. Phys. Chem. A 2011, 115, 10018−10026. (65) Fradelos, G.; Wesołowski, T. A. The Importance of Going beyond Coulombic Potential in Embedding Calculations for Molecular Properties: The Case of Iso-G for Biliverdin in Protein-Like Environment. J. Chem. Theory Comput. 2011, 7, 213−222. (66) Nåbo, L. J.; Olsen, J. M. H.; List, N. H.; Solanko, L. M.; Wüstner, D.; Kongsted, J. Embedding beyond electrostatics−The role of wave function confinement. J. Chem. Phys. 2016, 145, 104102. (67) List, N. H.; Jensen, H. J. Aa.; Kongsted, J. Local electric fields and molecular properties in heterogeneous environments through polarizable embedding. Phys. Chem. Chem. Phys. 2016, 18, 10070− 10080. (68) Neugebauer, J.; Louwerse, M. J.; Baerends, E. J.; Wesołowski, T. A. The merits of the frozen-density embedding scheme to model solvatochromic shifts. J. Chem. Phys. 2005, 122, 094115. (69) Lange, A.; Herbert, J. M. Simple Methods To Reduce ChargeTransfer Contamination in Time-Dependent Density-Functional Calculations of Clusters and Liquids. J. Chem. Theory Comput. 2007, 3, 1680−1690. (70) Jakobsen, S.; Kristensen, K.; Jensen, F. Electrostatic Potential of Insulin: Exploring the Limitations of Density Functional Theory and Force Field Methods. J. Chem. Theory Comput. 2013, 9, 3978−3985. (71) Isborn, C. M.; Mar, B. D.; Curchod, B. F. E.; Tavernelli, I.; Martínez, T. J. The Charge Transfer Problem in Density Functional Theory Calculations of Aqueously Solvated Molecules. J. Phys. Chem. B 2013, 117, 12189−12201.

5958


Polarizable Density Embedding: A Solution to the Electron Spill-Out

Recommend Documents