A polarizable three-layer frozen density embedding/molecular mechanics approach

The correct treatment of ever larger and more complex systems at an amenable computational cost has for long been at the forefront of quantum chemistry research. Alongside the refinement of wave function and density based methods with the purpose of increasing the accuracy of computational simulations, we have witnessed a significant amount of effort in the development of methodologies based on the partitioning of a chemical system into layers that can be treated separately and possibly at different levels of theory.^1–4 

Among the most well-known methods within this class are atomistic approaches based on a Quantum Mechanics (QM)/Molecular Mechanics (MM) paradigm.^5,6 For the latter model, the mutual polarization between the target (QM portion) and the environment (MM portion) may be of particular importance for the modeling of the properties of a molecular system, since in a much simpler electrostatic embedding scheme, the environmental response to a probing electromagnetic field could not be accurately taken into account.^7,8 Many polarizable QM/MM schemes have been presented based on different ways of modeling the polarization arising from the MM portion.^9–19 

Notwithstanding the significant improvements brought by polarizable atomistic models over the commonly used fixed-charge approaches, there are still significant drawbacks that still remain and are difficult to overcome due to the classical nature of these models. Non-electrostatic interactions are not easily integrated within these approaches. In fact, classical force fields model such effects with especially parameterized forms such as the Lennard–Jones potential;²⁰ however, these potentials are not especially useful in the context of QM/MM methods because they have no effect on the electronic density of the QM portion and thus on all properties that depend on it, from dipole moments to higher order response properties and spectroscopies. Ad hoc methods to include non-electrostatic interactions within QM/MM methods have been proposed,^21–23 but another possibility is to resort to QM/QM embedding methodologies to treat close-range interactions between the system and its environment.

Quantum mechanical embedding methods can, indeed, at least in principle, treat all possible types of interactions that may be present in the system.^24–45 Many different types of QM embedding methods have been developed, including projection based methods such as Hartree–Fock (HF)-in-Density Functional Theory (DFT) or DFT-in-DFT,^{32,39,40,46,47} Multilevel HF, DFT, or Coupled Cluster (CC),^{42–45,48–52} and the frozen density embedding (FDE) approach.^41,53–60 The most obvious drawback of quantum embedding schemes is their considerable cost, especially when compared with quantum/classical methods.

To address this problem, we propose a possible solution that consists of the definition of a multiple-layer method in which a small part of the environment and the target molecule are described by means of quantum embedding methods, while the rest of the environment is retained at the classical level only. Such a partitioning yields a three-layer approach, which we here formulate in the context of FDE theory coupled with the polarizable fluctuating charge (FQ) force field to describe the classical portion.⁷ In this way, the FDE layers allow for an accurate ab initio description of short-range interactions, mainly electrostatic and repulsive contributions, with the latter being well defined only at the QM level. Long-range interactions, which are mainly electrostatic in nature, are instead described at the classical level by means of the FQ force field. The purpose of this work is to overcome the limitations of a pure QM/classical embedding scheme when simulating, for instance, the properties of molecular solutes strongly interacting with the closed portion of the environment while at the same time minimizing the computational cost associated with the use of an intermediate QM layer.

Remarkably, the FQ method is a variational embedding approach, which is here chosen because it has been specifically developed to treat the properties of molecular solutes and has found significant success in the case of aqueous solutions, as detailed in recent reviews.^2,7,8 It is, however, worth noting that all variational embedding models are grounded on a common theoretical formulation as recently shown by Nottoli and Lipparini.⁶¹ Therefore, the FQ layer in the QM/FDE/FQ approach can be substituted with any other variational embedding model by keeping the same theoretical formulation that is here presented for the FQ case.

As a test case, QM/FDE/FQ is applied to the particular problem of simulating excitation energies of solutes in aqueous solution, for which an atomistic description of the environment is mandatory to accurately describe directional and specific interactions, such as hydrogen bonding. Although such interactions are dominated by electrostatics, quantum-based forces, such as Pauli repulsion, may play a crucial role in accurately reproducing experimental findings.^23,62,63

This manuscript is organized as follows: In Sec. II, we outline the theory and implementation of the QM/FDE/FQ method and its extension to calculate TDDFT excitation energies within a linear response regime. The methodological section is followed by the testing of the approach on the electronic absorption spectra of acrolein and 4-aminophthalimide in aqueous solution.

This section outlines the fundamentals of the QM/FDE/FQ approach by especially focusing on the coupling between QM/FDE and a layer described by means of the polarizable FQ force field. To this end, we first briefly recall and contextualize the polarizable QM/FQ approach, and then, we discuss the integrated three-layer embedding strategy. In fact, the key in the development of a multi-layer approach as the one proposed in this work is to carefully take into account and accurately describe all possible mutual interactions between the different layers.

In the two-layer QM/FQ approach, the system is divided into two regions, i.e., the target, treated at the QM level, and the environment, which is constituted by a set of molecules (e.g., a solvent) described by means of the FQ force field. The energy of the total system $E$ is written as follows:

where E_QM and E_FQ are the energies of the QM and FQ layers and $EQM/FQint$ is the QM/FQ interaction energy. Each atom belonging to the FQ portion of the system is endowed with a charge q whose value is not fixed (as in standard, non-polarizable force fields) but can vary according to the electronegativity equalization principle, which states that at equilibrium, each atom has the same electronegativity. By exploiting such a principle, E_FQ can be written as follows:

where M is a matrix constituted by the interaction kernel between the FQ charges J and a set of Lagrangian blocks (1_λ) assuring that the charge of each FQ molecule is fixed to Q, i.e., in this basic formulation, charge transfer between different molecules in the FQ region is not allowed. The self-interaction between the FQ charges, i.e., the diagonal elements of the interaction kernel in J, is expressed in terms of atomic chemical hardnesses. In Eq. (2), q_λ is a vector containing FQ charges and Lagrangian multipliers λ, whereas C_Q is a vector containing atomic electronegativities χ and the charge Q, which is maintained fixed on each molecule in the FQ layer.

The mutual polarization between the QM and FQ layers is taken into account by the interaction energy term $EQM/FQint$⁠, which is written in terms of the electrostatic interaction between the FQ charges q and the electrostatic potential generated by the QM density,

where V[ρ](r_i) is the electric potential generated by the QM density ρ on the i-th FQ charge (q_i), i.e.,

In Eq. (4), $ViN$ and $Vie$ are the nuclear and electronic terms, respectively. Z_α are the nuclear charges of N_QM atoms, which are placed at R_α. Notice that the interaction energy in Eq. (3) is written in terms of the electrostatic contribution only, which means that QM–MM Pauli repulsion and dispersion interactions are neglected. Such terms might be included by following an ad hoc approach that has recently been proposed in the literature.^23,62

By assuming the QM region to be described at the DFT level, Eq. (1) can be recast in the following functional:

where the dependence of the energy functional has been explicated. The values of FQ charges are obtained by minimizing Eq. (5) with respect to FQ variables, i.e., q and λ.¹³ Such a minimization procedure is equivalent to solving a set of linear equations, which reads as follows:

The last term in Eq. (3) accounts for the mutual polarization between the QM and FQ portions; therefore, a fully polarizable QM/FQ system is obtained.

FDE describes a QM complex system based on the partitioning of the electronic density into two subsystems I and II.⁶⁰ Each of the two subsystems (which can be two molecules, a solute and the surrounding solvent, or any other partition of a chemical system) is described by a separate electronic density. The total density ρ_tot of the whole system is then expressed as the sum of two densities ρ_I and ρ_II as follows:

The two densities can overlap and are constrained to integrate to an integer number of electrons. In addition to the electronic density, also the corresponding nuclei are apportioned accordingly. The FDE method is based on the assumption that one of the two densities may be kept frozen, while the total energy is variationally optimized with respect to the other density.^{38,41,53–59}

In order to build up a three-layer coupled QM/FDE/FQ system, Eq. (1) can be re-written as follows:

where E_QM/FDE, E_FQ, and $EQM/FDE/FQint$ are the FDE, FQ, and their interaction energies, respectively. In particular, by using Eq. (7), the QM/FDE energy reads as follows:

where T_s is the non-interacting kinetic energy, E_xc is the exchange–correlation energy, and v^nuc is the nuclear potential. Finally, the last two terms account for the non-additive terms that arise from the non-linearity of the DFT energy functionals. E_xc^nadd and T_s^nadd are defined as follows:

At this point, for a given ρ_II value, it is possible to determine the electronic density of subsystem I by minimizing the energy functional with respect to ρ_I itself, with the constraint that the number of electrons N_I is conserved. The density obtained can then also be determined through the minimization of the energy functional of this system as follows:

where $veffKSCED$ is the effective potential that allows the density of the non-interacting system to be the same as the interacting one, which reads as follows:

where $veffKS$ is the Kohn–Sham (KS) effective potential for ρ_I, whereas $veffembd$ is the embedding effective potential of ρ_I embedded in ρ_II.

Generally, FDE is used to model a small subsystem in a quite larger environment. For this reason, in order to reduce the impact of the calculation of the frozen density, the environment is often treated through approximated calculations.⁶⁴ Otherwise, it is possible to use the formalism described thus far to determine the electronic densities of both subsystems, in a manner known as “subsystem DFT.”^65,66

A valid method to achieve convergence uses a number of freeze-and-thaw cycles, in which the subsystems are, in turn, the frozen and non-frozen ones, until convergence is reached.⁶⁷ Alternatively, the equations can be simultaneously solved updating the densities after every cycle.

In order to determine the coupling between QM/FDE and the FQ layers, the interaction term $EQM/FDE/FQint$ in Eq. (8) can be defined, by recalling the partitioning in Eq. (7), as follows:

where the interaction term has been divided into two different contributions arising from both ρ_I and ρ_II. Therefore, Eq. (8) finally reads as follows:

where E_QM/FDE[ρ_I, ρ_II] is reported in Eq. (9). Similar to the basic QM/FQ approach exposed in Sec. II A, FQ charges are obtained by minimizing the functional in Eq. (14) with respect to the FQ variables, i.e., charges q and Lagrangian multipliers λ. The resulting linear system reads

where the linearity of the electrostatic potential with respect to the density has been exploited. Finally, a modified v_eff^KSCED potential can be defined by taking into account the interaction with FQ charges, which enters the definition of the embedding potential,

In addition to the calculation of energy, QM/FDE/FQ can be extended to compute molecular properties and spectroscopies in a similar way to other multilayer approaches.^2,8 As an example, here we focus on the excitation energies of embedded molecules by resorting to Time-Dependent Density Functional Theory (TDDFT). Extension to other properties/spectra is feasible and will, indeed, be the topic of further work by following what we have recently reviewed.⁷ The extension of QM/FDE/FQ to excited states needs to consider the response of the environment to the electronic excitation taking place in the target molecule, i.e., to the molecular density evolving in time as a result of the electronic excitation. The treatment of such a contribution has already been formulated for the QM/FQ model by some of the present authors by explicitly taking into account the FQ response to the molecular transition density.¹³ In particular, FQ charges are adjusted to the transition density associated with the specific electronic excitation under examination.

This results in the definition of modified Casida’s equations, which take into account FQ terms, i.e.,

where i and j correspond to the occupied Kohn–Sham orbitals, a and b correspond to the virtual ones, and ε represents the orbital energies; ω are the excitation energies.

The electrostatic forces between the solute and the FQ environment affect the description of the system by means of the $Ã$ and $B̃$ matrices, which are defined as follows:

where $Cai,bjQM/FQ$ are the additional FQ contributions to the $Ã$ and $B̃$ matrices and are defined as follows:

where q^T are the perturbed FQ charges (placed at positions r_p), which arise as a response to the transition density $PKT=XK+YK$⁠.¹³ φ are the MO orbitals.

In the case of the FDE model, instead, the definition of the charge density response to the time dependent external potential must take into account the partitioning of the system. To this end, the response of the embedding density is not considered, while it has been proven that the response of the target can be expressed in terms of that of the whole system modified by means of a functional for the non-additive kinetic energy term.⁶⁸ In this formalism, the spectroscopic response of the target does not explicitly depend on the environment, but it is indirectly still affected by it by means of its influence on the ground state density. Therefore, in the QM/FDE/FQ approach, the frozen density ρ_II only affects the MOs of subsystem I in the ground state, whereas the FQ response is adjusted to the transition density of the active portion only, as expressed in Eq. (19).

The presence of a frozen density layer that separates the QM system from the FQ layer results in a system where two layers that dynamically respond to the external time-dependent perturbation (the solute and the FQ) are separated by a layer that does not. This may result in an unbalanced picture and affect the quality of the results, particularly for those electronic transitions for which the environment response is especially important. A possible solution would be to calculate the response of the FDE layer as well by resorting to a coupled FDE model^68–70 and coupling the FQ layer to both QM portions and their response functions. We propose a different strategy: solely for the response equations, we embed the FDE layer with fluctuating charges just like the FQ layer and calculate their response using the same method. This amounts to modifying Eq. (19) by extending the set of charges to run over both the FQ and FDE atoms. This way the entire environment responds to the external perturbation. This choice has the advantage of not needing any expansion in Eq. (17) because the response of the frozen density layer is included automatically within the $Ã$ and $B̃$ matrices and not through the explicit excitations of its orbitals. The calculation of terms in Eq. (19) is very cheap; therefore, this does not affect the overall cost of the calculation in the slightest, and the resulting scheme, which we denote QM/FDE_FQ, offers the best of both worlds, i.e., allowing for the full inclusion of electrostatic and quantum repulsion effects of the environment upon the system in the ground state, without giving up the inclusion of electrostatic effects in the response calculations.

The QM/FDE/FQ approach has been implemented in a modified version of the Amsterdam Density Functional (ADF) 2019.3 program.^71,72 The ADF exploits Slater type orbitals, instead of the commonly used Gaussian type orbitals, and, as recalled in the name, it is designed in the DFT framework. For this reason, contrary to all previous implementations of the QM/FQ approach, due to the specificity of the ADF program, all the integrals that are involved in Eqs. (6), (4), (15), (16), and (19), i.e., the evaluation of the electrostatic potentials, are numerically calculated in the DFT grid. Furthermore, we remark that in our ADF implementation, we only deal with the density function on a set of grid points; thus, the density matrix of the QM part is never explicitly constructed. In addition, thanks to the modular structure of the ADF code, the FQ method has been implemented in a separate Fortran 90 module, which guarantees the coupling with existing features in the ADF package. Such a feature has been exploited in this work to allow for the coupling of the FQ layer with the QM/FDE description.

All calculations have been performed using a modified version of the ADF program.^71,72 In order to test the new coupled QM/FDE/FQ approach, the results have been obtained using varying partitions of the system. In addition, for the sake of comparison, some calculations employed a continuum modeling of the aqueous solution by means of the COSMO^73,74 approach.

The systems tested are aqueous solutions of acrolein and 4-Aminophthalimide (4AP) (see Fig. 1).

In the case of flexible solutes surrounded by several water molecules, a conformational sampling is needed to properly describe their configuration space.⁷ Classical molecular dynamics (MD) simulations of acrolein and 4AP in water were carried out as reported in Refs. 63 and 75, respectively, using Gromacs.⁷⁶ In the case of acrolein in water, 100 snapshots were extracted from the MD simulation and elaborated to be employed in QM/FDE/FQ calculations. In particular, the QM solute shell has always been composed of acrolein only, and solvent molecules have been chosen cutting a sphere of a fixed radius around it. For the QM/FQ calculations, the FQ layer was made up by water molecules that had at least one atom less distant than 17 Å from acrolein. For the QM/FDE calculations, subsystem I (i.e., the active QM) is constituted by the acrolein moiety, whereas only the water molecules placed at a distance lower than 3 Å are retained and treated as subsystem II (i.e., the frozen density). In the case of the QM/FDE/FQ approach, water molecules within 3 or 5 Å from acrolein were treated as subsystem II, while all the remaining solvent molecules constituted the FQ layer. With these choices, subsystem II contained, on average, 10.6 water molecules for the 3 Å threshold and 39.0 water molecules for the 5 Å threshold (see also Fig. 2).

For 4-Aminophthalimide (4AP), electronic excitation energies were calculated in the gas phase and in aqueous solution, describing the solvation phenomenon by means of QM/COSMO and QM/FDE/FQ approaches. In addition, in this case, for QM/FDE/FQ, 100 uncorrelated snapshots extracted from a classical MD simulation were considered. Solvation shells were defined by means of different cutting radii: 4AP and the water molecules closer than 3 Å from 4AP were treated as FDE, whereas all the remaining solvent molecules that are within 14 Å of the solute were described at the FQ level. With this partition, there were, on average, 17.2 water molecules in the frozen density layer.

For both the QM in vacuo and the QM/COSMO approaches, calculations were carried out on the optimized geometry of the molecules each at the respective level of theory.

In QM calculations, the hybrid B3LYP⁷⁷ exchange–correlation energy functional and the DZP basis set were used, while the Adiabatic Local Density Approximation (ALDA) kernel was employed for TDDFT calculations. The parameters for the FQ force field were taken from Ref. 78. Non-additive exchange–correlation and kinetic energies were evaluated in terms of PBE⁷⁹ and PW91k⁸⁰ functionals. For the freeze–thaw steps of the embedding density ρ_II, the PBE⁷⁹ exchange–correlation functional was used. In order to plot absorption spectra, the TDDFT excitations were convoluted using Gaussian lineshapes with a σ value of 0.25 eV for acrolein and 0.2 eV for 4AP.

In this section, we report the results of the application of the QM/FDE/FQ approach to describe a couple of organic molecules in aqueous solution. As it was anticipated in the Introduction, we focus on electronic excitation energies, which are especially suited as test properties for any focused approach. In fact, they arise from a specific portion of the system (the solute in our case) and are modified but not determined by the surrounding environment. The three-layer QM/FDE/FQ method is compared to alternative schemes ranging from full FDE computations to two-layer purely continuum QM/COSMO or QM/FQ approaches.

For the first application of the method, we selected acrolein (Fig. 1, left), and we focused on solvatochromic shifts due to water. This is a chemical system that has been studied extensively using quantum-chemical methods^63,81–85 because of the presence of different types of transitions in its electronic spectrum and the fact that it is composed of a few atoms, which limit the cost of QM calculations. In our case, the small size of the system allows us to test the approximations involved in different approaches. In particular, the aqueous solution was described by means of QM/COSMO, QM/FQ, QM/FDE, and QM/FDE/FQ approaches, and we also performed an in vacuo calculation to use as a reference and compare the results. Moreover, in calculations involving the FQ model coupled to the FDE and subsystem DFT methods, the solvent density has been treated as frozen and relaxed by means of freeze-and-thaw cycles, respectively.⁶⁷ 

Before reporting on the results of the different embedding schemes, we note that an important step in all atomistic spectral simulations is to ensure that the final spectra are at convergence with respect to the number of structures extracted from the MD. Therefore, the full acrolein–water configuration space has been sampled; otherwise, results risk being biased by an unrepresentative mean solvation environment. To this end, the spectra of acrolein in aqueous solution have been calculated using the QM/FDE/FQ model on an increasing number of snapshots to make sure that the spectrum does not significantly change upon addition of more snapshots to the ensemble.

As Fig. 3 shows, for the π$→$π^* transition, the energy of the maximum absorption peak does not considerably vary beyond 50 snapshots, while moving from 50 to 125 snapshots, there is a slight decrease in the intensity of the peak. A different behavior is, instead, reported for the n$→$π^* transition, where 50 snapshots are not appropriate to achieve convergence. In fact, at least 100 structures yield a converged spectrum, in terms of both energy (wavelength) and intensity. In both cases, the results obtained considering 125 snapshots do not appreciably differ from those obtained with 100 snapshots, so all the spectra presented in the following have been obtained from the average of 100 structures. Note that such findings in terms of the number of snapshots do not substantially differ from what we have reported for similar properties investigated at the QM/FQ level.^86–88 

Another point that is worth analyzing in greater detail comes from the fact that the FDE method in its basic form describes the environment of the system using a frozen density that is converged separately from the density of the solute.⁶⁷ This description may be refined by “thawing” the density through a number of freeze-and-thaw cycles in which the two layers switch role, with one being optimized and the other being kept frozen, until self-consistency is achieved. This approach, also known as subsystem DFT, allows one to take into account the polarization of the environment density caused by the system (e.g., acrolein in this case) at the expense of an increase in computational cost caused by the fact that several SCF cycles are needed to achieve convergence. While this approach would certainly seem more accurate, it has been frequently reported in the literature that this optimization does not improve the calculated properties of non-charged embedded species.^89,90 This is mainly due to the inaccuracies in describing the non-additive kinetic energy potential, which are enhanced by the relaxation of ρ_II. Nevertheless, because our three-layer QM/FDE/FQ method has never been examined form this point of view, it is worth investigating this aspect of the problem as well. As a matter of fact, spectra obtained for acrolein using the two aforementioned approaches do not significantly differ from one another, as can be seen in Fig. 4, for both investigated transitions. For this reason, the following discussion will only focus on spectra obtained by treating the solvent density as frozen. In the supplementary material, we also report the values of the energies and oscillator strengths for the two models (with the frozen and thawed density ρ_II) for each snapshot to show that the differences between the values are small on a per-snapshot basis rather than as a consequence of the averaging.

We finally come to the analysis of the results given by the different multilayer embedding models. Table I reports maximum absorption energies and solvatochromic shifts obtained from the averaged spectra of the 100 uncorrelated snapshots for both n → π^* and π → π^* transitions. We have applied QM/FQ, QM/FDE (including only water molecules within 3 Å of the solute’s atoms), and QM/FDE/FQ (with two different partitions at 3 and 5 Å from each atom for the FDE layer). At this stage, we also analyze the effect of including the solvent response within the FDE layer by endowing it with fluctuating charges as explained in Sec. II C. Such data are denoted as QM/FDE_FQ. If an external FQ layer is also included, we denote it as QM/FDE_FQ/FQ, and the radius used to select the boundary of the FDE region is also indicated in parenthesis. Table I also shows values obtained in vacuo and by using the continuum QM/COSMO approach and experimental data taken from Refs. 84 and 91.

It is worth noticing that water has opposite effects on the two peaks: while for the n → π^* transition, the intermolecular forces cause a small increase in the energy gap between the two states, for the π → π^* transition, water acts to reduce this gap. The inclusion of the solvent response in the FDE layer does not qualitatively alter this picture; however, it is immediately noticeable that it affects the shift of the two types of transitions in very different ways: there is almost no effect for the forbidden n → π^* transition, while a comparatively larger change is observed in the case of the bright π → π^* transition. This is to be expected since a brighter transition is expected to induce a much more potent electrostatic effect upon the solvent.

Analyzing more in detail each peak, Fig. 5 (top panel) shows the comparison of the results obtained with different models on the first absorption peak. All the models show a blueshift of the peak with respect to the vacuo result due to a larger stabilization of the n state with respect to the π^* one. Moreover, the QM/COSMO and QM/FQ models, which both only consider the electrostatic and polarization forces between acrolein and water, yield higher excitation maxima than the other models. In fact, QM/FQ also gives a higher excitation energy than QM/FDE since the latter lacks long-range solvent electrostatic effects. This behavior respects the physical interpretation of the system since the introduction of non-electrostatic forces is expected to stabilize the π^* state and so to result in a lower excitation energy, as can be found with the QM/FDE and QM/FDE/FQ models.⁶³ It is also worth comparing the QM/FDE/FQ model at 3 Å with the QM/FDE model, as the latter lacks the long-range polarization effects brought by FQ. As already noted, when comparing solvatochromic shifts, the inclusion of the solvent response within the FDE layer does not significantly affect this particular transition. The spectra, depicted in the inset of Fig. 6, show that this is true for both excitation energies and oscillator strengths, which are essentially the same for both the QM/FDE and QM/FDE_FQ models, and in fact, in the case where the boundary of the FDE region is 5 Å, the two spectra, with and without the response of the FDE region, are almost identical, while tiny differences are visible in the other cases. Note that the picture could be further improved by overcoming the limitations of a pure-linear response and instead adopt a corrected-linear response formalism for the evaluation of the solvent’s response.¹⁹ 

As for the second absorption peak, instead, it represents a π → π^* transition. In this case, Fig. 5 (bottom panel) shows that the introduction of the solvent causes a redshift of the peak, which is reported by each model. In addition, for this transition, both the electrostatic forces and the non-electrostatic ones seem to play an important role in the energies of the states involved. In fact, the increasing accuracy in the description of solute–solvent interactions, which is reached moving from the QM/COSMO to the QM/FDE/FQ model, results in an accentuated redshift of the peak. QM/FQ also yields better results than QM/COSMO; however, it should be emphasized that this model also benefits from cancellation of errors resulting from the absence of Pauli repulsion in the model, which lowers the agreement once included.⁶³ An estimation of long-range electrostatic effects can be made by comparing QM/FDE/FQ at 3 Å with the data obtained without the inclusion of the FQ layer. In fact, these effects shift the two transitions in opposite directions, as can be seen in Table I. Contrary to what was observed for the forbidden n$→$π^* peak, the inclusion of the solvent response in the FDE layer (Fig. 6) has a very visible effect for this transition notwithstanding the small number of water molecules included in the FDE layer. In addition to the redshift that was already discussed, we also see a significant increase in the oscillator strength, which brings the spectrum closer to the QM/FQ spectrum. Therefore, we can conclude that in general, the response of the solvent should not be neglected even when the FDE layer is made to be as small as possible and includes only those solvent molecules that are closest to the solute.

Moving to the comparison between computed and experimental data (see Table I), it is possible to see that the calculated shifts quite accurately reproduce the experimental trend for both peaks. In particular, for the n → π^* transition, the best agreement of the solvatochromic shift is provided by the QM/FDE/FQ and QM/FDE_FQ/FQ approaches in the case of a 5 Å FDE shell, which give similar results for the reasons discussed above. With these approaches, in fact, both electrostatic and non-electrostatic effects are considered since both are essential in this case.⁶³ 

On the other hand, the shift observed for the π → π^* transition seems to be more accurately reproduced by the QM/FQ model, while the models that include FDE seem to be less accurate. In addition, note that QM/FQ results often benefit from a lucky cancellation of errors since the underlying parameters that define the model (i.e., atomic electronegativities and hardnesses) can be optimized to yield the best results by decreasing the solvent polarization in order to compensate for the lack of non-electrostatic effects. For instance, when quantum repulsion effects are added back into the model, one finds that the electrostatic polarization should, indeed, be increased, as was demonstrated in previous studies.^19,63,86 The QM/FDE_FQ/FQ model includes both long-range electrostatic and short-range quantum repulsion effects and thus does not need to rely on cancellation of errors. A novel parameterization for the electrostatic part of the FQ may be desirable to be specifically applied to the combined three-layer model.

In order to further test the QM/FDE_FQ/FQ model, we applied it to the simulation of the electronic excitations of 4-Aminophthalimide (4AP) (Fig. 1, right), and the results were compared with in vacuo, QM/COSMO, and QM/FQ. FQ spectra were obtained by averaging the results over 100 uncorrelated snapshots extracted from the classical MD.

Finally, since the spectra obtained with the frozen and relaxed approaches for the embedding density do not significantly differ from one another, as demonstrated in Sec. IV A and also reported in the literature,^89,90 for this molecule, calculations were performed by only treating the water density as fully frozen.

In Fig. 7, spectra obtained at different levels of theory are plotted. The top panel shows the ensemble of all QM/FDE/FQ results arising from each of the 100 snapshots superposed with their convolution, showing how this approach naturally produces the shape of each band. Comparing the results in the bottom panel, the importance of the representation of the solvent is clearly evident. In fact, introducing the solvent into the description, even by means of an approximate method such as QM/COSMO, we observe a hypochromic effect and a redshift on both the calculated peaks. Since this approach only considers mean electrostatic solute–solvent forces, the shifts obtained can be read as indicators of the importance of these interactions for the 4AP–water system.

As the structure itself can suggest, with its nitrogen and oxygen atoms that can be involved in hydrogen bonding with the solvent and its diffuse π electron cloud, other interactions can also play a significant role. This observation can be easily confirmed in the spectra obtained with the QM/FQ and QM/FDE_FQ/FQ models, which, if compared to the QM/COSMO and in vacuo calculations, reveal notable shifts in the energies of the bands. These differences are clearly due to the non-electrostatic solute–solvent interactions and hydrogen bonding that a model such as QM/FDE_FQ/FQ can properly take into account and were neglected in the continuum solvation results. Comparing the results obtained with and without FDE, it can be seen that the inclusion of the intermediate quantum embedding layer leads to higher excitation energies closer to the vacuum values, as was observed in the case of acrolein.

We can then move on to comparing the calculated spectra with experiment, as shown in Fig. 8. The experimental spectrum, taken from Ref. 92, shows three absorption bands. Computational results seem to only reproduce two of them, one around 3.4 eV and the other around 4.9 eV, but not the smaller band at 4.1 eV. In fact, as seen in Fig. 7, the stick spectrum shows a faint grouping of signals precisely where this band should appear. In the gas-phase calculated spectrum, this band is not visible nor does it appear in any of the solvated results. We, therefore, tentatively attribute this discrepancy in the inability of the underlying electronic structure method to accurately reproduce the intensity of this band. Errors in the absolute excitation energies and intensities are to be expected due to the systematic errors introduced by the DFT description of the electronic structure, though environmental effects may still be accurately reproduced by the model. The analysis of excited states (Fig. 9) suggests that the first transition is due to the interaction between the electron rich region of the –NH₂ group and the electron-poor carbonyl group. This is experimentally confirmed comparing the 4-aminophthalimide and phthalimide spectra since in the latter, this peak is not found.⁹² 

This is also consistent with the assignment of the transitions based on the MO analysis (see Table II). The first transition is a HOMO–LUMO excitation, and we see in Fig. 9 that HOMO has a strong component on the –NH₂ group. The other transitions can be likewise assigned, with both visible transitions being characterized as π$→$π^* excitations. Overall, we can conclude that QM/FDE_FQ/FQ can effectively account for both long-range electrostatic interactions and short-range non-electrostatic effects, with visible differences in the predicted spectra.

We have presented a three-layer computational model in which the chemical system is divided into a central moiety treated quantum mechanically, while its close-range environment is modeled through the Frozen Density Embedding theory, and long-range interactions are included via the polarizable QM/MM model based on fluctuating charges (QM/FDE/FQ). This approach aims to overcome the main limitation of the quantum/classical QM/FQ description that is based solely on electrostatic interactions by including a middle layer that is also described quantum mechanically and therefore retains a description of non-electrostatic interactions. The FQ layer allows one to reduce the computational cost by reducing the size of the FDE region that must be included since long-range electrostatic interactions are accurately described by the classical polarizable embedding. One noteworthy feature of the presented model is the ability to include the environment’s response due to the FDE layer by using the classical FQ model. This is done by endowing each atom in the frozen density layer with a fluctuating charge, which only comes into play within the response equations. Thus, the ground state density of the quantum system is evaluated at the QM/FDE/FQ level of theory once the three layers have been defined, and then, response equations that are needed to simulate the electronic excitations of the system are solved by modeling the entire environment (FDE+FQ) using the FQ method. This approach allows us to consider the full electrostatic effects of the environment without the need to include the FDE layer in the response equations, which would significantly increase the computational cost while retaining all the beneficial features of the QM/FDE approach for the ground state, including the ability to treat non-electrostatic effects of quantum origin, which would be neglected if only the FQ model were employed.

The model was implemented and tested by simulating the electronic absorption spectra of organic molecules in water solution. For this type of system, the cybotactic region is treated using FDE, while the rest of the solvent is modeled via FQ. We have shown that the model can indeed reproduce solvatochromic effects, and therefore, it can be of great help whenever non-electrostatic effects, which are generally discarded in QM/classical calculations, play an important role.

While for the first application of the method, we opted to simulate the absorption of a single molecule in aqueous solution, we believe that the model possesses the versatility needed for the modeling of much more complex systems, where the presence of an intermediate layer modeled quantum mechanically is especially crucial. One example could be, for instance, molecules that interact with biochemical systems in solution or that are adsorbed on nanostructured plasmonic materials.^93–98 Other avenues that will be considered for exploration are simulations involving more complex response properties and spectroscopies, such as those involving magnetic perturbations like Nuclear Magnetic Resonance (NMR) or chiral properties such as circular dichroism.

See the supplementary material for comparison of the energies and oscillator strengths of acrolein computed with the QM/FDE/FQ method with and without relaxing the density of the FDE layer for each snapshot.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 818064). T.G. acknowledges funding from the Research Council of Norway through its grant TheoLight (Grant No. 275506).