Molecules in Environments: Towards Systematic Quantum Embedding of Electrons and Drude Oscillators

We develop a quantum embedding method that enables accurate and efficient treatment of interactions between molecules and an environment, while explicitly including many-body correlations. The molecule is composed of classical nuclei and quantum electrons, whereas the environment is modeled via charged quantum harmonic oscillators. We construct a general Hamiltonian and introduce a variational ansatz for the correlated ground state of the fully interacting molecule/environment system. This wavefunction is optimized via variational Monte Carlo and the ground state energy is subsequently estimated through diffusion Monte Carlo. The proposed scheme allows an explicit many-body treatment of electrostatic, polarization, and dispersion interactions between the molecule and the environment. We study solvation energies and excitation energies of benzene derivatives, obtaining excellent agreement with explicit ab initio calculations and experiment.

The challenges of computational electronic-structure methods have evolved from the precise description of short-range interactions involving few electrons to the accurate modeling of longer-range many-electron effects for molecules and materials [1][2][3].This shift has been enabled by the enormous progress in hardware capabilities combined with the development of new effective physical models and numerical techniques.
In particular, embedding approaches, such as the hybrid QM/MM  or QM/QM strategies , based on the partitioning of the system of interest in different subspaces, each treated with a different level of theory and thus of accuracy [1,2], have been devised to study excitation energies, solvation energies or charge transfer between fragments in a complex environment [6,11,18,24,28,38,44,51].
Although harmonic oscillators have been used in classical polarizable force fields [79][80][81], the novelty of our method lies in the quantum description of the full system of atoms and quantum Drude oscillators (QDOs) through one unified Hamiltonian for which we construct a ground state depending on all degrees of freedom.Furthermore, the proposed approach can be used with arbitrary ab initio methods, including CC [62] and QMC [63][64][65].
Within the Born-Oppenheimer approximation, each QDO consists of a pair of charged particles corresponding respectively to the center of the oscillator, with charge q and position R O and to a distinguishable quantum particle called drudon of charge -q, mass µ and position r d .The drudon interacts with its center via a quadratic potential V QDO (r) = 1 2 µω 2 r 2 that depends on the mass µ and on the frequency ω that determine the slope of the quadratic well, and on their relative Euclidean distance r = |r d − R O |.For systems containing more QDOs, all pairs of charged particles (except the drudon and its own center) interact via the Coulomb potential.In such a model the non-covalent interactions, such as dispersion, arise from the quantum fluctuations of the charge densities that are obtained through the coupling of the quantum states of each QDO.Naturally, each QDO depends on a set of well defined parameters {R O i , q i , µ i , ω i } N d i=1 that are obtained by reproducing the leading order polarizabilities and dispersion coefficients of real matter [70] from accurate ab initio or experimental data [67,74,75,82,83].
Once the parameters have been chosen, the analytical expression of the Hamiltonian of a system of N d interacting QDOs, the drudonic Hamiltonian, is written in atomic units as where that include the kinetic energy of the ith drudon and its interaction energies with the oscillators' centers.This drudonic Hamiltonian can be used to model the embedding potential that acts on a molecular system containing N n atomic nuclei and N e electrons, described respectively by the coordinate vectors Rn and re .The total Hamiltonian of electrons and QDOs (El-QDO) takes the form where Ĥe is the standard electronic Hamiltonian describing the interaction between the electrons and the atomic nuclei (that are parameters of the system) and V d−e int is the interaction between the QDOs and the atoms In order to describe molecular systems with static multipole moments of the charge distribution (such as the water molecule) it is necessary to introduce additional point charges in the drudonic Hamiltonian that modify the external potential felt by the drudons and by the electrons.
For example in the QDO model of water introduced by Jones et al. [73] and used in this work, the authors include in the model three additional point charges (centered on the two hydrogen atoms and near the oxygen atom), see Supplemental Material (SM) for details [84].These point charges do not interact with the QDOs of the corresponding molecule but only with the QDOs of the other water molecules, and in our case also with the electronic subsystem.A schematic representation of the interactions that are present in the total Hamiltonian is shown in Fig. 1.As explained below, the QDO model will be additionally corrected by an exchange-repulsion energy term that dominates the short-range interaction region [73,75,85].
The ground state of the many-body Hamiltonian defined in Eq. ( 3) is constructed through an appropriate variational ansatz including the explicit many-body correlations between all quantum particles.In order to integrate and optimize such an ansatz over the set of electronic and drudonic coordinates we choose to work in the framework of QMC methods [63,65,86], as implemented in the QMeCha package [87], although the approach can be generalized to other quantum-chemistry methods.QMC methods are stochastic integration techniques that are able to compute and minimize the energy functional over a chosen trial wave function (see SM [84]).In particular we apply the variational Monte Carlo (VMC) method [63,65,86], optimizing the wave function with the stochastic reconfiguration algorithm [88,89].Subsequently, to improve our energy estimations we use the diffusion Monte Carlo (DMC) projection method [63,65,86] that, within the fixed-node approximation used to solve the electronic sign problem, selects the ground state component of the approximate nodal-surface given by the optimized wave function (see SM [84]).Clearly, in the case of pure QDO systems in which the wave function is always positive and there is no nodal surface, the DMC algorithm converges to the exact solution.
We proceed to build the trial wave function for the system's ground state by expressing it as the product of three contributions which are respectively the electronic ansatz Ψ e (r e ), the drudonic one Ψ d (r d ) and a positive-definite function J d−e (r e , rd ), which contains the correlation effects between the electronic and drudonic systems..In this work, the pure electronic part Ψ e (r e ) is constructed through a Slater determinant of molecular orbitals formed as linear combinations of Gaussian basis sets, multiplied by an electronic Jastrow factor used to explicitly describe the many-body correlation effects between electrons and nuclei (see SM [84]).
The approximate form of the drudonic wave function Ψ d (r d ) is chosen to match the exact solution of a system of QDOs interacting through the dipole potential, which can be used to approximate the Coulomb potential at large distances as in the many-body dispersion MBD method [67,90].For this reason, it is written as the exponential function of a vector-matrix-vector product where rdO = rd − RO is the 3N d dimensional vector containing the distances between each drudon and its corresponding center, and A is a (symmetric) coupling matrix containing 3N d (3N d + 1) /2 independent parameters.The coupling matrix explicitly correlates the drudons subject to the field of their charged centers.
Finally, the electron-drudon factor J e−d (r e , rd ) represents the coupling between the two types of quantum systems.Again we hypothesize that the interactions between the QDO environment and molecular one are dominated by dipole interactions, thus the ansatz will resemble the form in Eq. ( 6) coupling the distance vector rdO of the QDO's with the dipole moment of the molecular part through the rectangular matrix B containing a set of 3 × 3N d free parameters.We remark that the restriction of the QDO wave function to the dipolar approximation is not a significant limitation, since the variational optimization of the wave function parameters is followed by a stochastic DMC calculation of the many-body energy for both electrons and drudons.First, we apply our approach to the Ar [91] and the water dimers [92].While these systems are small, when taken together they are challenging benchmark examples, requiring an accurate description of exchange, dispersion, polarization, and hydrogen bonding.For both Ar and water, a single QDO represents the response of all 8 valence electrons.Accordingly, the QDO model allows a significant reduction in the electronic degrees of freedom.We remark that a bosonic QDO model can be extended to describe exchange [75,93], but this would require a wavefunction with mixed quantum statistics for electrons and QDOs.Instead, in our current embedding scheme we describe exchange via parametrized short-range exponentials (see SM [84]).
In Ar 2 [91], the interaction energy comes from an interplay between correlation in the long range and exchange in the short range.The attractive dispersion energy is correctly described by the QDO-QDO model and essentially coincides with its counterpart computed from electronic symmetry-adapted perturbation theory  (SAPT) [94] (see SM [84]).Noticeably, the electrons-QDO (El-QDO) embedding model reproduces the full interaction curve upon subtracting the exchange contributions from SAPT, highlighting the ability of our method to capture electrostatic and polarization terms.Together with the corresponding repulsive part, both models match the reference CCSD(T) binding curve [91], as shown in Fig. 2.
For the water dimer, where intermolecular interactions are strongly directional, we apply the El-QDO embedding model to study the dissociation energies of its two lowest-energy geometric configurations, C1 and C2, along the O-O axis (Fig. 3).Here we use the extended QDO model of Martyna et al. [73] that contains three point charges and a QDO located at the center-of-mass, substituting the acceptor water molecule (always on the right in Fig. 3).From Fig. 3 one observes that both C1 and C2 binding curves obtained with the QDO-QDO and El-QDO models are in excellent agreement with the CCSD(T) references [92].The models are thus able to describe the different interactions when varying the molecular orientation.Moreover, the comparison of the bare El-QDO interaction curve of the C1 conformer with the SAPT decomposition from Ref. 95 (see SM [84]), shows that our model is again close to the full binding energy curve, upon subtracting the repulsive exchange contributions, thus correctly capturing both dispersion and polarization effects.From the atomic and molecular examples of Ar and water dimers, we conclude that our embedding approach is able to accurately describe the manifold of interactions between an electronic subsystem and an embedding environment, in which the latter is described through a combination of QDOs and point charges.
Among the most intriguing applications of embedding approaches are the effects of solvent on the binding energies between molecules and on their electronic excitations [96,97].Here, we study the binding in Tshaped benzene dimer, prototype for C-H−π interactions [96,97], and the Triplet-Singlet excitation energy (T-S) of ortho-benzyne, which is the non-diradical singlet conformer of benzyne whose excitation energy requires a proper description of dynamical correlation as well as solvent polarization [98].These two systems are studied in cages of water from 4 to 50 molecules (see Fig. 4) with the aim of probing collective polarization and dispersion effects.In order to avoid empirical treatment of exchange in these many-body systems, we impose a minimal distance of ∼ 3 Å between the solute and water.
First we compute the solvation energies, E Solv (X) = E (X Cage )−E (X Vacuum )−E Cage , where X is the T-shaped benzene dimer or its monomers (Table I), in a cage composed of 50 water molecules (structures c) and d) in Fig. 4), or the ortho-benzyne in its singlet (S) and triplet (T) spin states (Table II) in a cage of 4 and 30 water molecules (structures a) and b) in Fig. 4).For all systems the El-QDO DMC results are compared with state-ofthe-art DFT calculations using the PBE0 functional [99] augmented with the Tkatchenko-Scheffler pairwise dispersion energy (PBE0+TS) [74] or the many-body dispersion method (PBE0+MBD) [67,84,90].
From the results reported in Table I and Table II we conclude that El-QDO DMC solvation energies are in excellent agreement with fully electronic DMC calculations (for the 4W cage for which they are feasible) and with PBE0+MBD results.It is known that PBE0+MBD binding energies are often quantitatively accurate for molecular systems, especially when compared to electronic DMC calculations [53].In fact, we note that manybody effects are particularly evident for ortho-benzyne, where the PBE0+TS method overestimates the energies when compared to the El-QDO DMC and PBE0+MBD methods, which both include, in a different manner, the many-body correlation effects of the environment on the solute.Table I.Solvation energies (in kcal/mol) of the T-shaped benzene dimer and monomers in the 50W water cage [100] (Molecular structures c) and d) in Fig. 4).Table II.Solvation energies (in kcal/mol) of orthobenzyne [98] in its singlet (S) and triplet (T) spin states in cages of four (4W) and thirty (30W) water molecules (Molecular structures a) and b) in Fig. 4).Table III.Singlet-Triplet adiabatic excitation energies (S-T) of the ortho-benzyne [98] in vacuum (V), 4-water cage (4W) and 30-water cage (30W) (in kcal/mol), corresponding to the molecular structures a) and b) in Fig.The key role of proper quantum embedding achieved with the El-QDO DMC approach can be observed in the singlet-triplet (S-T) adiabatic excitation energy of ortho-benzyne (Table III).Due to delocalization error, DFT severely underestimates the S-T gap [98], while electronic DMC yields accurate results [102] in vacuum in excellent agreement with the experimental value of 37.5(3) kcal/mol [101].For the smallest cage of four water molecules, the excitation energies computed with electronic DMC and the El-QDO DMC approach are in excellent agreement between each other.Yet, the error bar of El-QDO calculations is much lower compared to electronic DMC computations for an equivalent sampling trajectory.Moreover, as Table III shows, when increasing the size of the water cage, the excitation energy changes by 0.2(1) kcal/mol.This effect, captured by the El-QDO DMC embedding method arises from the higher polarization of the T state with respect to the S state of benzyne.We expect that future extensions of our method will allow more realistic embedding geometries where water molecules can come closer to the solute, yielding more substantial solvation effects.Remarkably, the computational cost of the El-QDO method is essentially independent on the number of QDOs in the embedding region (see Fig. 4) with respect to the cost of the molecule in vacuum.For the 4-water cage, the cost increases by only ∼1%, and by ∼16% for the 30-water cage.This stems from the lack of nodal structure in the QDO wavefunctions and low noise in their DMC sampling with respect to the electronic wavefunctions.
Our study illustrates the need for fully quantum embedding methods when treating binding and excitations energies of solvated systems, and introduces a computationally efficient approach.The El-QDO embedding method is a firm first step in this direction as we demonstrated on the example of prototypical dispersion and hydrogen-bonded dimers and for the solvation of benzene derivates.A future challenge includes the incorporation of a more general representation of the QDO environment and exchange effects, which could be achieved through a synergy with machine-learning methods [103][104][105][106].
MD and MB thank Jorge Alfonso Charry Martinez for participating in the development of the QMeCha code.MB acknowledges financial support from the Luxembourg National Research Fund (IN-TER/DFG/18/12944860).The DFT calculations presented in this paper were carried out using the HPC facilities of the University of Luxembourg [107] (see hpc.uni.lu).An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program.This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.

Figure 1 .
Figure 1.Schematic representation of the embedding model for a single atom and two Drude oscillators.Strings indicate the harmonic interactions between the drudon and their corresponding centers, the dashed lines represent the Coulomb interactions between all other pairs of charged particles (for simplicity, only selected interactions are explicitly shown).

Figure 2 .
Figure 2. Binding energy curve for the Ar dimer obtained using the QDO-QDO model and electrons-QDO (El-QDO) embedding approach (schematically shown in the figure) at the DMC level of theory, with the exponential fit of the short-range repulsion [84].The results are compared to the CCSD(T) [91] reference curve.

Figure 3 .
Figure3.Binding energy curve for the water dimer as a function of the oxygen-oxygen distance in the two energetically lowest geometric configurations (C1 and C2)[92] obtained using the QDO-QDO and El-QDO (with QDO on the right, approximating the acceptor) models at the DMC level of theory[84].The results are compared to the CCSD(T)[92] reference curve.

Figure 4 .
Figure 4. Relative runtimes (see SM [84]) of water-embedded El-QDO DMC calculations with respect to the molecular energies computed with DMC in vacuum.Runtime for electronic DMC is also shown for ortho-benzyne in 4W cage.Inset figures a) to d) correspond to the molecular structures [98, 100].