We present a novel multilayer polarizable embedding approach in which the system is divided into three portions, two of which are treated using density functional theory and their interaction is based on frozen density embedding (FDE) theory, and both also mutually interact with a polarizable classical layer described using an atomistic model based on fluctuating charges (FQ). The efficacy of the model is demonstrated by extending the formalism to linear response properties and applying it to the simulation of the excitation energies of organic molecules in aqueous solution, where the solute and the first solvation shell are treated using FDE, while the rest of the solvent is modeled using FQ charges.

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;20 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.7 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.61 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:

E=EQM+EFQ+EQM/FQint,
(1)

where EQM and EFQ 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, EFQ can be written as follows:

EFQ=12qλMqλ+qλCQ=12qλJ1λ1λ0qλ+qλχQ,
(2)

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 CQ 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,

EQM/FQint=iqiV[ρ](ri),
(3)

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

V[ρ](ri)=ViN+Vie=α=1NQMZα|riRα|ρ(r)|rir|dr.
(4)

In Eq. (4), ViN and Vie are the nuclear and electronic terms, respectively. Zα are the nuclear charges of NQM 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:

E[ρ,q,λ]=EQM[ρ]+12qλMqλ+qλCQ+qλV(ρ),
(5)

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 λ.13 Such a minimization procedure is equivalent to solving a set of linear equations, which reads as follows:

Mqλ=CQV(ρ).
(6)

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

ρtot=ρI+ρII.
(7)

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:

E=EQM/FDE+EFQ+EQM/FDE/FQint,
(8)

where EQM/FDE, EFQ, 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:

EQM/FDE[ρI,ρII]=Ts[ρI]+Ts[ρII]+Exc[ρI]+Exc[ρII]+(ρI(r)+ρII(r))(vInuc(r)+vIInuc(r))dr+12(ρI(r)+ρII(r))(ρI(r)+ρII(r))|rr|drdr+Tsnadd[ρI,ρII]+Excnadd[ρI,ρII],
(9)

where Ts is the non-interacting kinetic energy, Exc is the exchange–correlation energy, and vnuc 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. Excnadd and Tsnadd are defined as follows:

Excnadd[ρI,ρII]=Exc[ρI+ρII]Exc[ρI]Exc[ρII],Tsnadd[ρI,ρII]=Ts[ρI+ρII]Ts[ρI]Ts[ρII].
(10)

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 NI is conserved. The density obtained can then also be determined through the minimization of the energy functional of this system as follows:

Es[ρI]=Ts[ρI]+ρI(r)veffKSCED[ρI,ρII](r)dr,
(11)

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:

veffKSCED[ρI,ρII](r)=veffKS[ρI](r)+veffemb[ρI,ρII](r),
(12)

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.64 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.67 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:

EQM/FDE/FQint=EI/FQ+EII/FQ,EI/FQ=iqiV[ρI](ri),EII/FQ=iqiV[ρII](ri),
(13)

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

E[ρI,ρII,q,λ]=EQM/FDE[ρI,ρII]+12qλMqλ+qλCQ+qλV(ρI)+qλV(ρII),
(14)

where EQM/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

Mqλ=CQV(ρI)V(ρII)=CQV(ρtot),
(15)

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

veffKSCED[ρI,ρII,q](r)=veffKS[ρI](r)+veffemb[ρI,ρII,q](r)=veffKS[ρI](r)+veffemb, FDE[ρI,ρII](r)iqi|rri|.
(16)

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

ÃB̃Ã*B̃*XY=ω1001XY,
(17)

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:

Ãai,bj=δijδab(εaεi)+aj1|rr|+δ2Exc[ρ]δρ(r)δρ(r)ib+Cai,bjQM/FQ,B̃ai,bj=ab1|rr|+δ2Exc[ρ]δρ(r)δρ(r)ij+Cai,bjQM/FQ,
(18)

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

Cai,bjFQ=pR3φa(r)1rrpφi(r)drqpT(φb,φi),
(19)

where qT are the perturbed FQ charges (placed at positions rp), which arise as a response to the transition density PKT=XK+YK.13 φ 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.68 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 model68–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/FDEFQ, 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 COSMO73,74 approach.

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

FIG. 1.

Molecular structures of acrolein (left) and 4AP (right).

FIG. 1.

Molecular structures of acrolein (left) and 4AP (right).

Close modal

In the case of flexible solutes surrounded by several water molecules, a conformational sampling is needed to properly describe their configuration space.7 Classical molecular dynamics (MD) simulations of acrolein and 4AP in water were carried out as reported in Refs. 63 and 75, respectively, using Gromacs.76 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).

FIG. 2.

Pictorial view of the various levels of theory exploited to study acrolein in aqueous solution.

FIG. 2.

Pictorial view of the various levels of theory exploited to study acrolein in aqueous solution.

Close modal

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 B3LYP77 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 PBE79 and PW91k80 functionals. For the freeze–thaw steps of the embedding density ρII, the PBE79 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 methods63,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.67 

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 

FIG. 3.

Comparison between the QM/FDE/FQ absorption spectra of acrolein in aqueous solution obtained by averaging an increasing number of snapshots. Spectra arising from both nπ* and ππ* transitions are shown.

FIG. 3.

Comparison between the QM/FDE/FQ absorption spectra of acrolein in aqueous solution obtained by averaging an increasing number of snapshots. Spectra arising from both nπ* and ππ* transitions are shown.

Close modal

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

FIG. 4.

Comparison between the QM/FDE/FQ absorption spectra of acrolein in aqueous solution, averaged over the 100 snapshots, obtained with either a frozen or relaxed density ρII. Both nπ* and ππ* transitions are shown.

FIG. 4.

Comparison between the QM/FDE/FQ absorption spectra of acrolein in aqueous solution, averaged over the 100 snapshots, obtained with either a frozen or relaxed density ρII. Both nπ* and ππ* transitions are shown.

Close modal

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/FDEFQ. If an external FQ layer is also included, we denote it as QM/FDEFQ/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.

TABLE I.

Maximum absorption energies and vacuo-to-water solvatochromic shifts (eV) for the two main electronic transitions of acrolein. Experimental values are taken from Refs. 84 and 91.

nπ*ππ*
ModelAbs. ene.ShiftAbs. ene.Shift
In vacuo 3.65 ⋯ 6.40 ⋯ 
Expt. (gas phase) 3.69 ⋯ 6.42 ⋯ 
QM/COSMO 3.95 −0.30 6.30 0.10 
QM/FQ 4.05 −0.40 5.90 0.50 
QM/FDE (3 Å) 3.80 −0.15 6.10 0.30 
QM/FDE/FQ (3 Å) 3.95 −0.30 6.00 0.40 
QM/FDE/FQ (5 Å) 3.90 −0.25 6.05 0.35 
QM/FDEFQ (3 Å) 3.82 −0.17 6.05 0.35 
QM/FDEFQ/FQ (3 Å) 3.94 −0.29 5.94 0.46 
QM/FDEFQ/FQ (5 Å) 3.90 −0.25 5.96 0.44 
Expt. (aqueous solution) 3.94 −0.25 5.90 0.52 
nπ*ππ*
ModelAbs. ene.ShiftAbs. ene.Shift
In vacuo 3.65 ⋯ 6.40 ⋯ 
Expt. (gas phase) 3.69 ⋯ 6.42 ⋯ 
QM/COSMO 3.95 −0.30 6.30 0.10 
QM/FQ 4.05 −0.40 5.90 0.50 
QM/FDE (3 Å) 3.80 −0.15 6.10 0.30 
QM/FDE/FQ (3 Å) 3.95 −0.30 6.00 0.40 
QM/FDE/FQ (5 Å) 3.90 −0.25 6.05 0.35 
QM/FDEFQ (3 Å) 3.82 −0.17 6.05 0.35 
QM/FDEFQ/FQ (3 Å) 3.94 −0.29 5.94 0.46 
QM/FDEFQ/FQ (5 Å) 3.90 −0.25 5.96 0.44 
Expt. (aqueous solution) 3.94 −0.25 5.90 0.52 

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.63 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/FDEFQ 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.19 

FIG. 5.

Comparison between the absorption spectra of acrolein obtained with different embedding approaches.

FIG. 5.

Comparison between the absorption spectra of acrolein obtained with different embedding approaches.

Close modal
FIG. 6.

Comparison between the QM/FDE and QM/FDE/FQ spectra of acrolein in water, with and without the inclusion of the solvent response in the FDE layer (QM/FDEFQ).

FIG. 6.

Comparison between the QM/FDE and QM/FDE/FQ spectra of acrolein in water, with and without the inclusion of the solvent response in the FDE layer (QM/FDEFQ).

Close modal

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.63 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/FDEFQ/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.63 

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

FIG. 7.

Comparison between the absorption spectra of 4AP obtained with different levels of theory. The top panel also shows the stick spectrum, i.e., the unbroadened ensemble of the oscillator strengths obtained from each snapshot at their respective peak positions.

FIG. 7.

Comparison between the absorption spectra of 4AP obtained with different levels of theory. The top panel also shows the stick spectrum, i.e., the unbroadened ensemble of the oscillator strengths obtained from each snapshot at their respective peak positions.

Close modal

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/FDEFQ/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/FDEFQ/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 –NH2 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.92 

FIG. 8.

Comparison between the calculated absorption spectra of 4AP in water with different models and the experimental spectrum.

FIG. 8.

Comparison between the calculated absorption spectra of 4AP in water with different models and the experimental spectrum.

Close modal
FIG. 9.

Frontier orbitals of 4AP in water evaluated at the B3LYP/COSMO level of theory.

FIG. 9.

Frontier orbitals of 4AP in water evaluated at the B3LYP/COSMO level of theory.

Close modal

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 –NH2 group. The other transitions can be likewise assigned, with both visible transitions being characterized as ππ* excitations. Overall, we can conclude that QM/FDEFQ/FQ can effectively account for both long-range electrostatic interactions and short-range non-electrostatic effects, with visible differences in the predicted spectra.

TABLE II.

Excited states of 4AP evaluated at the B3LYP/COSMO level of theory. H: HOMO and L: LUMO.

StateEnergy (eV)f/10−2Orbitals
3.46 7.033 H → L 
4.22 0.012 H − 2 → L 
4.42 0.430 H − 1 → L, H → L + 1 
4.86 0.001 H − 4 → L 
5.14 41.810 H → L + 1, H − 1 → L 
5.52 10.830 H − 3 → L 
StateEnergy (eV)f/10−2Orbitals
3.46 7.033 H → L 
4.22 0.012 H − 2 → L 
4.42 0.430 H − 1 → L, H → L + 1 
4.86 0.001 H − 4 → L 
5.14 41.810 H → L + 1, H − 1 → L 
5.52 10.830 H − 3 → L 

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

1.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
, “
Quantum mechanical continuum solvation models
,”
Chem. Rev.
105
,
2999
3094
(
2005
).
2.
C.
Cappelli
, “
Integrated QM/polarizable MM/continuum approaches to model chiroptical properties of strongly interacting solute-solvent systems
,”
Int. J. Quantum Chem.
116
,
1532
1542
(
2016
).
3.
A.
Warshel
and
M.
Levitt
, “
Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme
,”
J. Mol. Biol.
103
,
227
249
(
1976
).
4.
M. J.
Field
,
P. A.
Bash
, and
M.
Karplus
, “
A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations
,”
J. Comput. Chem.
11
,
700
733
(
1990
).
5.
H. M.
Senn
and
W.
Thiel
, “
QM/MM methods for biomolecular systems
,”
Angew. Chem., Int. Ed.
48
,
1198
1229
(
2009
).
6.
H.
Lin
and
D. G.
Truhlar
, “
QM/MM: What have we learned, where are we, and where do we go from here?
,”
Theor. Chem. Acc.
117
,
185
199
(
2007
).
7.
T.
Giovannini
,
F.
Egidi
, and
C.
Cappelli
, “
Molecular spectroscopy of aqueous solutions: A theoretical perspective
,”
Chem. Soc. Rev.
49
,
5664
5677
(
2020
).
8.
T.
Giovannini
,
F.
Egidi
, and
C.
Cappelli
, “
Theory and algorithms for chiroptical properties and spectroscopies of aqueous systems
,”
Phys. Chem. Chem. Phys.
22
,
22864
22879
(
2020
).
9.
C.
Curutchet
,
A.
Muñoz-Losa
,
S.
Monti
,
J.
Kongsted
,
G. D.
Scholes
, and
B.
Mennucci
, “
Electronic energy transfer in condensed phase studied by a polarizable QM/MM model
,”
J. Chem. Theory Comput.
5
,
1838
1848
(
2009
).
10.
J. M. H.
Olsen
and
J.
Kongsted
, “
Molecular properties through polarizable embedding
,”
Adv. Quantum Chem.
61
,
107
143
(
2011
).
11.
D.
Loco
,
É.
Polack
,
S.
Caprasecca
,
L.
Lagardère
,
F.
Lipparini
,
J.-P.
Piquemal
, and
B.
Mennucci
, “
A QM/MM approach using the AMOEBA polarizable embedding: From ground state energies to electronic excitations
,”
J. Chem. Theory Comput.
12
,
3654
3661
(
2016
).
12.
T.
Giovannini
,
A.
Puglisi
,
M.
Ambrosetti
, and
C.
Cappelli
, “
Polarizable QM/MM approach with fluctuating charges and fluctuating dipoles: The QM/FQFμ model
,”
J. Chem. Theory Comput.
15
,
2233
2245
(
2019
).
13.
F.
Lipparini
,
C.
Cappelli
, and
V.
Barone
, “
Linear response theory and electronic transition energies for a fully polarizable QM/classical Hamiltonian
,”
J. Chem. Theory Comput.
8
,
4153
4165
(
2012
).
14.
F.
Lipparini
,
C.
Cappelli
,
G.
Scalmani
,
N.
De Mitri
, and
V.
Barone
, “
Analytical first and second derivatives for a fully polarizable QM/classical Hamiltonian
,”
J. Chem. Theory Comput.
8
,
4270
4278
(
2012
).
15.
A. H.
Steindal
,
K.
Ruud
,
L.
Frediani
,
K.
Aidas
, and
J.
Kongsted
, “
Excitation energies in solution: The fully polarizable QM/MM/PCM method
,”
J. Phys. Chem. B
115
,
3027
3037
(
2011
).
16.
T.
Schwabe
,
J. M. H.
Olsen
,
K.
Sneskov
,
J.
Kongsted
, and
O.
Christiansen
, “
Solvation effects on electronic transitions: Exploring the performance of advanced solvent potentials in polarizable embedding calculations
,”
J. Chem. Theory Comput.
7
,
2209
2217
(
2011
).
17.
D.
Loco
,
S.
Jurinovich
,
L.
Cupellini
,
M. F. S. J.
Menger
, and
B.
Mennucci
, “
The modeling of the absorption lineshape for embedded molecules through a polarizable QM/MM approach
,”
Photochem. Photobiol. Sci.
17
,
552
560
(
2018
).
18.
T.
Giovannini
,
L.
Grazioli
,
M.
Ambrosetti
, and
C.
Cappelli
, “
Calculation of IR spectra with a fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles
,”
J. Chem. Theory Comput.
15
,
5495
5507
(
2019
).
19.
T.
Giovannini
,
R. R.
Riso
,
M.
Ambrosetti
,
A.
Puglisi
, and
C.
Cappelli
, “
Electronic transitions for a fully polarizable QM/MM approach based on fluctuating charges and fluctuating dipoles: Linear and corrected linear response regimes
,”
J. Chem. Phys.
151
,
174104
(
2019
).
20.
J. E.
Lennard-Jones
, “
Cohesion
,”
Proc. Phys. Soc.
43
,
461
(
1931
).
21.
H.
Gökcan
,
E. G.
Kratz
,
T. A.
Darden
,
J.-P.
Piquemal
, and
G. A.
Cisneros
, “
QM/MM simulations with the Gaussian electrostatic model, a density-based polarizable potential
,”
J. Phys. Chem. Lett.
9
,
3062
3067
(
2018
).
22.
M. S.
Gordon
,
M. A.
Freitag
,
P.
Bandyopadhyay
,
J. H.
Jensen
,
V.
Kairys
, and
W. J.
Stevens
, “
The effective fragment potential method: A QM-based MM approach to modeling environmental effects in chemistry
,”
J. Phys. Chem. A
105
,
293
307
(
2001
).
23.
T.
Giovannini
,
P.
Lafiosca
, and
C.
Cappelli
, “
A general route to include Pauli repulsion and quantum dispersion effects in QM/MM approaches
,”
J. Chem. Theory Comput.
13
,
4854
4870
(
2017
).
24.
M. S.
Gordon
,
Q. A.
Smith
,
P.
Xu
, and
L. V.
Slipchenko
, “
Accurate first principles model potentials for intermolecular interactions
,”
Annu. Rev. Phys. Chem.
64
,
553
578
(
2013
).
25.
M. S.
Gordon
,
L.
Slipchenko
,
H.
Li
, and
J. H.
Jensen
, “
The effective fragment potential: A general method for predicting intermolecular interactions
,”
Annu. Rep. Comput. Chem.
3
,
177
193
(
2007
).
26.
Q.
Sun
and
G. K.-L.
Chan
, “
Quantum embedding theories
,”
Acc. Chem. Res.
49
,
2705
2712
(
2016
).
27.
G.
Knizia
and
G. K.-L.
Chan
, “
Density matrix embedding: A strong-coupling quantum embedding theory
,”
J. Chem. Theory Comput.
9
,
1428
1432
(
2013
).
28.
D. V.
Chulhai
and
J. D.
Goodpaster
, “
Projection-based correlated wave function in density functional theory embedding for periodic systems
,”
J. Chem. Theory Comput.
14
,
1928
1942
(
2018
).
29.
D. V.
Chulhai
and
J. D.
Goodpaster
, “
Improved accuracy and efficiency in quantum embedding through absolute localization
,”
J. Chem. Theory Comput.
13
,
1503
1508
(
2017
).
30.
X.
Wen
,
D. S.
Graham
,
D. V.
Chulhai
, and
J. D.
Goodpaster
, “
Absolutely localized projection-based embedding for excited states
,”
J. Chem. Theory Comput.
16
,
385
398
(
2020
).
31.
F.
Ding
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Embedded mean-field theory with block-orthogonalized partitioning
,”
J. Chem. Theory Comput.
13
,
1605
1615
(
2017
).
32.
J. D.
Goodpaster
,
T. A.
Barnes
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Density functional theory embedding for correlated wavefunctions: Improved methods for open-shell systems and transition metal complexes
,”
J. Chem. Phys.
137
,
224113
(
2012
).
33.
J. D.
Goodpaster
,
T. A.
Barnes
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Accurate and systematically improvable density functional theory embedding for correlated wavefunctions
,”
J. Chem. Phys.
140
,
18A507
(
2014
).
34.
F. R.
Manby
,
M.
Stella
,
J. D.
Goodpaster
, and
T. F.
Miller
 III
, “
A simple, exact density-functional-theory embedding scheme
,”
J. Chem. Theory Comput.
8
,
2564
2568
(
2012
).
35.
J. D.
Goodpaster
,
N.
Ananth
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Exact nonadditive kinetic potentials for embedded density functional theory
,”
J. Chem. Phys.
133
,
084103
(
2010
).
36.
K.
Zhang
,
S.
Ren
, and
M.
Caricato
, “
Multi-state QM/QM extrapolation of UV/Vis absorption spectra with point charge embedding
,”
J. Chem. Theory Comput.
16
,
4361
4372
(
2020
).
37.
P.
Ramos
,
M.
Papadakis
, and
M.
Pavanello
, “
Performance of frozen density embedding for modeling hole transfer reactions
,”
J. Phys. Chem. B
119
,
7541
7557
(
2015
).
38.
M.
Pavanello
and
J.
Neugebauer
, “
Modelling charge transfer reactions with the frozen density embedding formalism
,”
J. Chem. Phys.
135
,
234103
(
2011
).
39.
P. K.
Tamukong
,
Y. G.
Khait
, and
M. R.
Hoffmann
, “
Density differences in embedding theory with external orbital orthogonality
,”
J. Phys. Chem. A
118
,
9182
9200
(
2014
).
40.
P. K.
Tamukong
,
Y. G.
Khait
, and
M. R.
Hoffmann
, “
Accurate dissociation of chemical bonds using DFT-in-DFT embedding theory with external orbital orthogonality
,”
J. Phys. Chem. A
121
,
256
264
(
2017
).
41.
T. A.
Wesolowski
, “
Hydrogen-bonding-induced shifts of the excitation energies in nucleic acid bases: An interplay between electrostatic and electron density overlap effects
,”
J. Am. Chem. Soc.
126
,
11444
11445
(
2004
).
42.
S.
Sæther
,
T.
Kjærgaard
,
H.
Koch
, and
I.-M.
Høyvik
, “
Density-based multilevel Hartree–Fock model
,”
J. Chem. Theory Comput.
13
,
5282
5290
(
2017
).
43.
I.-M.
Høyvik
, “
Convergence acceleration for the multilevel Hartree–Fock model
,”
Mol. Phys.
118
,
1626929
(
2020
).
44.
G.
Marrazzini
,
T.
Giovannini
,
M.
Scavino
,
F.
Egidi
,
C.
Cappelli
, and
H.
Koch
, “
Multilevel density functional theory
,”
J. Chem. Theory Comput.
17
,
791
(
2021
).
45.
T.
Giovannini
and
H.
Koch
, “
Energy-based molecular orbital localization in a specific spatial region
,”
J. Chem. Theory Comput.
17
,
139
150
(
2021
).
46.
S. J.
Bennie
,
B. F. E.
Curchod
,
F. R.
Manby
, and
D. R.
Glowacki
, “
Pushing the limits of EOM-CCSD with projector-based embedding for excitation energies
,”
J. Phys. Chem. Lett.
8
,
5559
5565
(
2017
).
47.
S. J.
Lee
,
M.
Welborn
,
F. R.
Manby
, and
T. F.
Miller
 III
, “
Projection-based wavefunction-in-DFT embedding
,”
Acc. Chem. Res.
52
,
1359
1368
(
2019
).
48.
R. H.
Myhre
,
A. M. J.
Sánchez de Merás
, and
H.
Koch
, “
Multi-level coupled cluster theory
,”
J. Chem. Phys.
141
,
224105
(
2014
).
49.
R. H.
Myhre
and
H.
Koch
, “
The multilevel CC3 coupled cluster model
,”
J. Chem. Phys.
145
,
044111
(
2016
).
50.
S. D.
Folkestad
and
H.
Koch
, “
Multilevel CC2 and CCSD methods with correlated natural transition orbitals
,”
J. Chem. Theory Comput.
16
,
179
189
(
2019
).
51.
S. D.
Folkestad
and
H.
Koch
, “
Equation-of-motion MLCCSD and CCSD-in-HF oscillator strengths and their application to core excitations
,”
J. Chem. Theory Comput.
16
,
6869
6879
(
2020
).
52.
S. D.
Folkestad
,
E. F.
Kjønstad
,
L.
Goletto
, and
H.
Koch
, “
Multilevel CC2 and CCSD in reduced orbital spaces: Electronic excitations in large molecular systems
,”
J. Chem. Theory Comput.
17
,
714
(
2021
).
53.
J.
Neugebauer
,
M. J.
Louwerse
,
E. J.
Baerends
, and
T. A.
Wesolowski
, “
The merits of the frozen-density embedding scheme to model solvatochromic shifts
,”
J. Chem. Phys.
122
,
094115
(
2005
).
54.
T. A.
Wesolowski
,
S.
Shedge
, and
X.
Zhou
, “
Frozen-density embedding strategy for multilevel simulations of electronic structure
,”
Chem. Rev.
115
,
5891
5928
(
2015
).
55.
S.
Fux
,
C. R.
Jacob
,
J.
Neugebauer
,
L.
Visscher
, and
M.
Reiher
, “
Accurate frozen-density embedding potentials as a first step towards a subsystem description of covalent bonds
,”
J. Chem. Phys.
132
,
164101
(
2010
).
56.
C. R.
Jacob
,
J.
Neugebauer
, and
L.
Visscher
, “
A flexible implementation of frozen-density embedding for use in multilevel simulations
,”
J. Comput. Chem.
29
,
1011
1018
(
2008
).
57.
C. R.
Jacob
and
L.
Visscher
, “
Calculation of nuclear magnetic resonance shieldings using frozen-density embedding
,”
J. Chem. Phys.
125
,
194104
(
2006
).
58.
C. R.
Jacob
,
J.
Neugebauer
,
L.
Jensen
, and
L.
Visscher
, “
Comparison of frozen-density embedding and discrete reaction field solvent models for molecular properties
,”
Phys. Chem. Chem. Phys.
8
,
2349
2359
(
2006
).
59.
T. A.
Wesołowski
, “
Embedding a multideterminantal wave function in an orbital-free environment
,”
Phys. Rev. A
77
,
012504
(
2008
).
60.
T. A.
Wesolowski
and
A.
Warshel
, “
Frozen density functional approach for ab initio calculations of solvated molecules
,”
J. Phys. Chem.
97
,
8050
8053
(
1993
).
61.
M.
Nottoli
and
F.
Lipparini
, “
General formulation of polarizable embedding models and of their coupling
,”
J. Chem. Phys.
153
,
224108
(
2020
).
62.
P.
Arora
,
L. V.
Slipchenko
,
S. P.
Webb
,
A.
DeFusco
, and
M. S.
Gordon
, “
Solvent-induced frequency shifts: Configuration interaction singles combined with the effective fragment potential method
,”
J. Phys. Chem. A
114
,
6742
6750
(
2010
).
63.
T.
Giovannini
,
M.
Ambrosetti
, and
C.
Cappelli
, “
Quantum confinement effects on solvatochromic shifts of molecular solutes
,”
J. Phys. Chem. Lett.
10
,
5823
5829
(
2019
).
64.
T. A.
Wesolowski
and
A.
Warshel
, “
Ab initio free energy perturbation calculations of solvation free energy using the frozen density functional approach
,”
J. Phys. Chem.
98
,
5183
5187
(
1994
).
65.
P.
Cortona
, “
Self-consistently determined properties of solids without band-structure calculations
,”
Phys. Rev. B
44
,
8454
8458
(
1991
).
66.
C. R.
Jacob
and
J.
Neugebauer
, “
Subsystem density-functional theory
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
325
362
(
2014
).
67.
T. A.
Wesolowski
and
J.
Weber
, “
Kohn-Sham equations with constrained electron density: An iterative evaluation of the ground-state electron density of interacting molecules
,”
Chem. Phys. Lett.
248
,
71
76
(
1996
).
68.
M. E.
Casida
and
T. A.
Wesołowski
, “
Generalization of the Kohn–Sham equations with constrained electron density formalism and its time-dependent response theory formulation
,”
Int. J. Quantum Chem.
96
,
577
588
(
2004
).
69.
J.
Neugebauer
, “
Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory
,”
J. Chem. Phys.
126
,
134116
(
2007
).
70.
J.
Neugebauer
, “
Photophysical properties of natural light-harvesting complexes studied by subsystem density functional theory
,”
J. Phys. Chem. B
112
,
2207
2217
(
2008
).
71.
G.
te Velde
,
F. M.
Bickelhaupt
,
E. J.
Baerends
,
C.
Fonseca Guerra
,
S. J. A.
van Gisbergen
,
J. G.
Snijders
, and
T.
Ziegler
, “
Chemistry with ADF
,”
J. Comput. Chem.
22
,
931
967
(
2001
).
72.
ADF 2020, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, http://www.scm.com.
73.
A.
Klamt
, “
The COSMO and COSMO-RS solvation models
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
699
709
(
2011
).
74.
C. C.
Pye
and
T.
Ziegler
, “
An implementation of the conductor-like screening model of solvation within the amsterdam density functional package
,”
Theor. Chem. Acc.
101
,
396
(
1999
).
75.
T.
Giovannini
,
M.
Macchiagodena
,
M.
Ambrosetti
,
A.
Puglisi
,
P.
Lafiosca
,
G.
Lo Gerfo
,
F.
Egidi
, and
C.
Cappelli
, “
Simulating vertical excitation energies of solvated dyes: From continuum to polarizable discrete modeling
,”
Int. J. Quantum Chem.
119
,
e25684
(
2019
).
76.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
, “
GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers
,”
SoftwareX
1-2
,
19
25
(
2015
).
77.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
, “
Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density
,”
Phys. Rev. B
37
,
785
789
(
1988
).
78.
S. W.
Rick
,
S. J.
Stuart
, and
B. J.
Berne
, “
Dynamical fluctuating charge force fields: Application to liquid water
,”
J. Chem. Phys.
101
,
6141
6156
(
1994
).
79.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
80.
A.
Lembarki
and
H.
Chermette
, “
Obtaining a gradient-corrected kinetic-energy functional from the Perdew-Wang exchange functional
,”
Phys. Rev. A
50
,
5328
(
1994
).
81.
A. V.
Marenich
,
C. J.
Cramer
,
D. G.
Truhlar
,
C. A.
Guido
,
B.
Mennucci
,
G.
Scalmani
, and
M. J.
Frisch
, “
Practical computation of electronic excitation in solution: Vertical excitation model
,”
Chem. Sci.
2
,
2143
2161
(
2011
).
82.
I.
Duchemin
,
C. A.
Guido
,
D.
Jacquemin
, and
X.
Blase
, “
The Bethe–Salpeter formalism with polarisable continuum embedding: Reconciling linear-response and state-specific features
,”
Chem. Sci.
9
,
4430
4443
(
2018
).
83.
C.
Bistafa
,
L.
Modesto-Costa
, and
S.
Canuto
, “
A complete basis set study of the lowest n–π* and π–π* electronic transitions of acrolein in explicit water environment
,”
Theor. Chem. Acc.
135
,
129
(
2016
).
84.
K.
Aidas
,
A.
Møgelhøj
,
E. J. K.
Nilsson
,
M. S.
Johnson
,
K. V.
Mikkelsen
,
O.
Christiansen
,
P.
Söderhjelm
, and
J.
Kongsted
, “
On the performance of quantum chemical methods to predict solvatochromic effects: The case of acrolein in aqueous solution
,”
J. Chem. Phys.
128
,
194503
(
2008
).
85.
L.
Goletto
,
T.
Giovannini
,
S. D.
Folkestad
, and
H.
Koch
, “
Combining multilevel Hartree Fock and multilevel coupled cluster with molecular mechanics: A study of electronic excitations in solutions
,”
Phys. Chem. Chem. Phys.
23
,
4413
4425
(
2021
).
86.
G.
Marrazzini
,
T.
Giovannini
,
F.
Egidi
, and
C.
Cappelli
, “
Calculation of linear and non-linear electric response properties of systems in aqueous solution: A polarizable quantum/classical approach with quantum repulsion effects
,”
J. Chem. Theory Comput.
16
,
6993
7004
(
2020
).
87.
T.
Giovannini
,
M.
Ambrosetti
, and
C.
Cappelli
, “
A polarizable embedding approach to second harmonic generation (SHG) of molecular systems in aqueous solutions
,”
Theor. Chem. Acc.
137
,
74
(
2018
).
88.
T.
Giovannini
,
M.
Olszówka
,
F.
Egidi
,
J. R.
Cheeseman
,
G.
Scalmani
, and
C.
Cappelli
, “
Polarizable embedding approach for the analytical calculation of Raman and Raman optical activity spectra of solvated systems
,”
J. Chem. Theory Comput.
13
,
4421
4435
(
2017
).
89.
M.
Humbert-Droz
,
X.
Zhou
,
S. V.
Shedge
, and
T. A.
Wesolowski
, “
How to choose the frozen density in frozen-density embedding theory-based numerical simulations of local excitations?
,”
Theor. Chem. Acc.
133
,
1405
(
2013
).
90.
G.
Fradelos
,
J. J.
Lutz
,
T. A.
Wesołowski
,
P.
Piecuch
, and
M.
Włoch
, “
Embedding vs supermolecular strategies in evaluating the hydrogen-bonding-induced shifts of excitation energies
,”
J. Chem. Theory Comput.
7
,
1647
1666
(
2011
).
91.
A.
Moskvin
,
O.
Yablonskii
, and
L.
Bondar
, “
An experimental investigation of the effect of alkyl substituents on the position of the K and R absorption bands in acrolein derivatives
,”
Theor. Exp. Chem.
2
,
469
472
(
1966
).
92.
L.
Bucsiová
and
P.
Hrdlovič
, “
Medium effect of polymer matrices on spectral properties of 4-aminophthalimide and 4-dimethylaminophthalimide
,”
J. Macromol. Sci., Part A
44
,
1047
1053
(
2007
).
93.
T.
Giovannini
,
M.
Rosa
,
S.
Corni
, and
C.
Cappelli
, “
A classical picture of subnanometer junctions: An atomistic Drude approach to nanoplasmonics
,”
Nanoscale
11
,
6004
6015
(
2019
).
94.
T.
Giovannini
,
L.
Bonatti
,
M.
Polini
, and
C.
Cappelli
, “
Graphene plasmonics: Fully atomistic approach for realistic structures
,”
J. Phys. Chem. Lett.
11
,
7595
7602
(
2020
).
95.
S. M.
Morton
and
L.
Jensen
, “
A discrete interaction model/quantum mechanical method to describe the interaction of metal nanoparticles and molecular absorption
,”
J. Chem. Phys.
135
,
134103
(
2011
).
96.
S. M.
Morton
,
D. W.
Silverstein
, and
L.
Jensen
, “
Theoretical studies of plasmonics using electronic structure methods
,”
Chem. Rev.
111
,
3962
3994
(
2011
).
97.
S. M.
Morton
and
L.
Jensen
, “
A discrete interaction model/quantum mechanical method for describing response properties of molecules adsorbed on metal nanoparticles
,”
J. Chem. Phys.
133
,
074103
(
2010
).
98.
L.
Jensen
,
C. M.
Aikens
, and
G. C.
Schatz
, “
Electronic structure methods for studying surface-enhanced Raman scattering
,”
Chem. Soc. Rev.
37
,
1061
1073
(
2008
).

Supplementary Material