We present a novel quantum mechanical/molecular mechanics (QM/MM) approach in which a quantum subsystem is coupled to a classical subsystem described by the AMOEBA polarizable force field. Our approach permits mutual polarization between the QM and MM subsystems, effected through multipolar electrostatics. Self-consistency is achieved for both the QM and MM subsystems through a total energy minimization scheme. We provide an expression for the Hamiltonian of the coupled QM/MM system, which we minimize using gradient methods. The QM subsystem is described by the onetep linear-scaling DFT approach, which makes use of strictly localized orbitals expressed in a set of periodic sinc basis functions equivalent to plane waves. The MM subsystem is described by the multipolar, polarizable force field AMOEBA, as implemented in tinker. Distributed multipole analysis is used to obtain, on the fly, a classical representation of the QM subsystem in terms of atom-centered multipoles. This auxiliary representation is used for all polarization interactions between QM and MM, allowing us to treat them on the same footing as in AMOEBA. We validate our method in tests of solute-solvent interaction energies, for neutral and charged molecules, demonstrating the simultaneous optimization of the quantum and classical degrees of freedom. Encouragingly, we find that the inclusion of explicit polarization in the MM part of QM/MM improves the agreement with fully QM calculations.

## I. INTRODUCTION

Combining quantum-mechanical (QM) calculations with a classical description is a well-established technique in computational studies of molecular systems.^{1–7} The computational cost of purely QM approaches becomes prohibitive for many systems. One reason is the long-range nature of Coulombic interactions, which makes many properties of interest converge slowly with the size of the system. This makes it necessary to include hundreds of atoms,^{8} if not more,^{9} in the calculation before acceptable accuracy can be reached,^{10} with suitable truncation of larger systems often far from trivial.^{11} Another difficulty arises for systems which cannot be well represented by a single conformer, necessitating a statistical averaging of properties over a large number of configurations.^{12} Linear-scaling approaches^{13–20} can help ameliorate the length scale problem, but they do not address the effects of conformational sampling, which may become more important with the increase of the system size.^{21} While molecular mechanics (MM) approaches can be routinely applied to systems comprising ∼10^{5} to 10^{6} atoms for timescales in the order of 1 *μ*s, the purely classical description is inherently unable to describe intrinsically electronic processes (such as bond breaking and bond reconstruction) and properties (e.g., band gaps or solvent shifts).

In many cases, the properties of interest are localized to a certain part of the system (usually a molecule) that is embedded in a greater environment (such as solvent or solid state matrix). Although typically the presence of the environment cannot be neglected, it is often the case that we are not interested in its detailed properties or behavior. This is the *raison d’être* of focused models such as QM/MM, which seek to describe the subsystem of interest at a desired level of theory, while representing the environment only approximately.^{12}

Hybrid (quantum-classical) calculations fall into two main categories, depending on how they describe this environment. QM/continuum approaches employ an averaged description, with the environment lacking any internal structure. This strategy is best exemplified by widely used implicit solvation techniques, such as PCM,^{22,23} COSMO,^{24,25} or density-dependent models,^{26–28} which help to address the sampling problem by reducing the dimensionality of the system. QM/MM approaches instead adopt a molecular description of the environment, retaining its atomistic detail and thereby preserving directional and specific interactions, such as hydrogen bonds between solvent and solute.^{29} Methodologies combining all three descriptions (e.g., QM/MM/PCM) have recently been reported.^{12,30–32}

The relative ease with which QM/MM methods enable a meaningful compromise between efficiency and accuracy is now well-accepted in biomolecular chemistry,^{33} particularly in applications where electronic excitations are of interest.^{12,30,34} Consequently, a wide gamut of approaches of varying sophistication, and targeting different classes of systems has been proposed^{12,21,29–57} since the pioneering work of Warshel and Levitt.^{1} While even the briefest of reviews is beyond the scope of this paper, an interested reader may benefit from a presentation of model hierarchies,^{58} a review of QM/MM methods for biomolecular systems^{2} and for materials science,^{59} and a recent overview of applications of QM/MM in enzymology.^{60}

In parallel to the progress of QM/MM techniques, the last two decades have witnessed significant developments in the sophistication of force fields,^{61} following the identification of deficiencies in commonly used fixed point charge models. The way in which molecules respond to environmental conditions, such as the presence of a solvent, pH, or ion concentration is, unsurprisingly, difficult to capture without taking polarization into account.^{62} In the absence of explicit polarization terms, fixed point charge models can only describe polarization effects in an averaged fashion, through the reparametrization of the pairwise energy terms that they have at their disposal. This crude description may be adequate under conditions that are sufficiently close to those for which the potential has been parametrized, but, more often than not, it is poorly transferable to different phases or different environments, such as interfaces.^{63}

The non-additive, many-body nature of polarization interactions makes such models more involved and computationally demanding compared to traditional fixed point charge models. In consequence, a wide variety of competing treatments of polarization has emerged (see Refs. 62 and 64 for a review), employing: Drude oscillators,^{65,66} fluctuating charges,^{67,68} and induced point dipoles^{69–75} or higher multipoles.^{76}

Correspondingly, a number of distinct QM/MM models employing polarizable force fields have since been proposed. These models differ in their approach to MM polarization (Drude oscillators,^{50,52} fluctuating charges,^{12,51,55,77} induced dipoles^{21,29–31,33,34,37,48,54,56,78}), and the sophistication of their treatment of the non-polarizable (permanent) part of MM electrostatics (partial charges only,^{12,29–31,37,48,50–52,55,77,78} higher-order multipoles^{21,32,34,56}). Several groups have developed models specifically focused on electronic excitations, using polarizable embedding alongside time-dependent density functional theory (TDDFT),^{12,21,29–33,55} where dynamic mutual polarization^{33} poses an additional challenge. Other models restrict their scope to ground-state density functional theory (DFT). Some of the models further extend their classical description by coupling to a continuum dielectric.^{12,30–32,51,54,78} Other desirable features include analytical gradients,^{31,48,55} the ability to treat covalent bonds spanning the QM/MM interface,^{31,48,50,56} and the availability of post-Hartree-Fock methods for the QM part.^{31,34,53,55,79}

In this work we present a new QM/MM approach, which combines the DFT methodology of onetep ^{80,81} and the polarizable force-field AMOEBA.^{70,72,73} The QM and MM subsystems are coupled electrostatically, and undergo mutual polarization. The electrostatic effect of the MM subsystem is included in the QM Hamiltonian, polarizing the QM subsystem by deforming its electronic charge density. Conversely, the electric field of the QM subsystem is included in the direct field that drives the polarization of the MM subsystem.

A crucial element of our electrostatics model is an auxiliary representation of the QM system in terms of point, atom-centered multipoles, which is used in the calculation of polarization interactions. This representation is obtained through a variant of distributed multipole analysis (DMA),^{82,83} which is a technique for partitioning charge density into single-atom contributions. Invoking this intermediate classical representation for the QM system allows our model to describe QM/MM polarization interactions on the same footing as in AMOEBA, that is using damped, point-multipole electrostatics.

A second distinguishing feature of our approach is the use of linear-scaling DFT^{80,84} to describe the QM subsystem with the aim of, ultimately, undertaking QM/polarizable-MM calculations with QM regions spanning thousands of atoms. To the best of our knowledge, there have been no reports to date of a polarizable QM/MM model with a linear-scaling QM component, although we note in passing that a QM/MM approach with linear-scaling MM (in CPU and memory use) has recently been reported.^{54}

In our approach, van der Waals interactions in the MM subsystem use the unmodified AMOEBA buffered 7-14 formalism,^{73,85} and the same pairwise functional form is used for QM/MM interactions. In the QM subsystem, the repulsive term is naturally accounted for through the exchange-correlation term, while standard empirical corrections^{86} are used to account for dispersion.

We derived a total energy expression for the entire system (QM+MM). This energy is iteratively minimized under suitable constraints using gradient methods^{87} until the simultaneous self-consistency of both subsystems is achieved.

We envision this paper to be the first in a series of continued developments, serving as a proof of concept. For this reason we will primarily focus on a detailed description of theory behind our model, and its validation.

This paper is organized as follows. In Sec. II we describe our approach—starting from a description of the QM and MM components, we follow with an explanation of how they are coupled. Section III is devoted to validation and demonstration of the utility of the proposed approach. Conclusions and closing remarks will be found in Section IV.

## II. METHOD

### A. Conventions and notation

We follow the sign convention where electrons are positively charged. Atomic units are used throughout the text, unless indicated otherwise. Quantities typeset in bold denote Cartesian column vectors (positions **r**, electric fields **E**, dipoles $ \mu $, etc.), Cartesian tensors of rank 2 (e.g., $ T L M d-d $) and of rank 3 ($ T L M d-q $). Matrices with dimensions other than 3 × 3 are typeset with blackboard capitals (e.g., 𝕂). Indices *I* and *J* always refer to atoms in the QM subsystem, and indices *L* and *M* refer to atoms in the MM subsystem. Localized orbitals are indexed with Greek symbols. By van der Waals interactions we will mean the sum of the repulsive and dispersive terms, referring to the attractive term simply as “dispersion.”

### B. QM component

In this section we describe how the QM subsystem is treated in our model. We employ the linear-scaling onetep ^{20,80,88–91} approach, which is a reformulation of Kohn-Sham DFT in terms of the one-electron density matrix, $\rho r , r \u2032 $. The density matrix is expressed in a separable form

in terms of non-orthogonal, localized orbitals (support functions) $ \phi \alpha r $, and a density kernel $K= K \alpha \beta $. The support functions, termed non-orthogonal generalized Wannier functions (NGWFs),^{20} are strictly localized within atom-centered spherical regions. The density kernel is the matrix representation of the density matrix (1) in the duals of the NGWFs.

Linear-scaling is achieved by making the NGWF localization regions finite, and by truncating the density matrix beyond a chosen cutoff distance, exploiting the “nearsightedness principle” of electronic matter.^{92} The assumption of finite localization regions makes the overlap matrix $S= S \alpha \beta $, whose elements are given by

sparse, while the truncation of the density matrix introduces sparsity into the density kernel 𝕂.

Unlike in most linear-scaling approaches, in onetep the NGWFs themselves are expanded in an underlying basis of periodic sinc^{93} (psinc) functions, which are equivalent to plane waves. Gradient methods^{87} are used to minimize the total energy by optimizing not only the elements of the density kernel, but also the expansion coefficients of the NGWFs in terms of psincs. As the NGWFs are optimized *in situ*, using a minimal NGWF basis is sufficient for obtaining high accuracy and systematic convergence of total energies to those of a plane-wave approach with Kohn-Sham molecular orbitals. Alternatively, the NGWFs can be kept fixed following suitable initialization—e.g., to pseudoatomic orbitals (PAOs), or to orbitals that have been pre-optimized in advance. This simpler approach obviates the requirement of deriving and computing energy gradients with respect to the NGWFs. In the interest of simplicity in this communication we limit ourselves to this latter approach.

onetep permits calculations both on periodic systems (in the Γ-point approximation) and on isolated (non-periodic) systems through the use of open boundary conditions (OBCs) and a selection of techniques for eliminating the effects of undesired periodicity.^{94} The model described here uses the latter methodology, and, unless indicated otherwise, we shall assume in the text that follows that OBCs are in effect.

For an isolated QM system, the minimized quantity is the usual Born-Oppenheimer DFT energy of electrons in the potential of clamped nuclei, which in the density-matrix formulation is given by

where $\rho r , r \u2032 $ is the density matrix (1), $ v ext r $ is the external potential of the QM ions, $ E XC QM $ is the exchange-correlation energy functional, and $n r $ is the electronic density, given by

where we assume a closed-shell system for simplicity, with the factor of 2 that is a consequence of spin degeneracy accounted for in the scaling of 𝕂. We note that in practice the pseudopotential approach is used, and so $ v ext r $ would correspond to the pseudopotential of the QM *cores*, and $n r $—to the pseudodensity of *valence* electrons.

The remaining energy terms associated with the QM subsystem are independent of the electronic degrees of freedom. The first of these is the Coulombic repulsion energy of the cores, which, under open boundary conditions is simply

where *N*_{QM} is the number of atoms in the QM subsystem, $ R I I = 1 N QM $ are the positions of the cores, and $ Z I I = 1 N QM $ are their charges. The last term, $ E disp QM $, is an empirical dispersion-correction term, which accounts for the well-known deficiency of generalized gradient approximation (GGA) DFT in describing dispersion interactions.^{86} The exact expression depends on the model used, but the general form is that of a double sum of pairwise terms. This work uses the Elstner^{95} formulation, with parameters determined by Hill and Skylaris.^{86} The total energy associated with the QM subsystem is the sum of the three above-mentioned terms,

### C. MM component

In our model the MM subsystem is described with an unmodified AMOEBA^{73} polarizable force-field, as implemented in the tinker ^{71} code. AMOEBA is a successful polarizable force field offering a consistent treatment of electrostatic interactions through permanent multipoles up to a quadrupole, with polarization realized via damped, induced, point dipoles. The accuracy of the AMOEBA description and parametrization has been demonstrated^{72} for a variety of properties (i.a. heats of vaporization, dimer binding energies, vibrational frequencies, solvation free energies). Its canonical implementation, tinker, provides the means for extracting individual MM energy terms and induced dipoles, which greatly simplifies interfacing with QM. In this section we give a brief account of those aspects of AMOEBA that are relevant to our QM/MM approach and we introduce necessary terminology.

The AMOEBA energy expression can be written in the general form

with the four energy components respectively accounting for permanent electrostatic interactions, polarization, short-range valence interactions, and van der Waals interactions.

#### 1. Permanent electrostatics

In its treatment of electrostatic interactions, AMOEBA uses a point multipole representation. Each atomic site *L* is associated with permanent, point multipoles up to a quadrupole: $ { q L p , \mu L p , Q L p } $. The permanent multipoles are parametrized from *ab initio* calculations,^{96–98} and use a suitable local coordinate frame^{70} to maintain transferability between different molecular conformations. Interactions between permanent multipoles on different sites are purely Coulombic, i.e., they are not damped. However, to ensure a smooth transition between an electrostatic description of interactions at medium range and bonded interactions at short range, permanent electrostatic interactions between nearest and second-nearest neighbors (as determined by bond connectivity) are zeroed, and corresponding interactions between third- and fourth-nearest neighbors are attenuated,^{72} in what is known as *scaling* or *masking*. We shall denote the permanent electrostatic energy of the MM subsystem by $ E perm MM $, readers interested in the full expression can consult Ref. 73, eqs. 1 and 10.

#### 2. Polarization

In addition to permanent multipoles, AMOEBA associates an induced dipole $ \mu L $ with each atomic site *L*. These dipoles are induced primarily by the electric field of the permanent multipoles, termed the direct field, in what is known as *direct induction*. AMOEBA uses an interactive induction scheme, whereby each induced dipole $ \mu L $ will further polarize induced dipoles at other sites *M*≠*L*. Such *mutual induction* continues until the induced dipoles at each site reach convergence.^{73} Induced dipole self-consistency is usually achieved through the use of iterative solver techniques, such as SOR,^{71,99} although a variety of alternative schemes is also in use.^{100–102} For a more detailed discussion the reader is referred to Ref. 103. At convergence the dipole induced at site *L* through linear response is simply

where *α _{L}* is the atomic polarizability of site

*L*, and

**E**

_{L}and $ E L m $ are, respectively, the direct and the mutual electric fields at

**R**

_{L}.

AMOEBA uses a non-additive polarization model, where mutual polarization takes place between all polarizable sites, even those belonging to the same molecule.^{71} This means that the polarizability of a molecule is not a sum of atomic polarizabilities of its constituent atoms; and that intramolecular polarization needs to be accounted for (removed) during the parametrization of permanent multipoles. An important advantage of the non-additive model is that it makes the atomic polarizability tensors isotropic, which is why the polarizability *α _{L}* featuring in (8) is a scalar value.

In order to avoid a well-known deficiency of point polarizability models known as the “polarization catastrophe”^{104} (an unbounded mutual polarization of nearby dipoles), AMOEBA damps all interactions involving induced dipoles. The damping model devised by Thole^{105} and subsequently revised^{106} modifies the dipole-dipole interaction tensor, ensuring interaction energies approach a finite value rather than becoming infinite as the distance between the two dipoles approaches zero. Analogous modifications^{107} to dipole-charge and dipole-quadrupole interaction tensors make it possible to similarly damp the interactions of induced dipoles with permanent multipoles used in AMOEBA. The relevant modified Cartesian interaction tensors will be given later in the text, cf. Eqs. (27)–(29).

Other modifications to electrostatic interactions involving induced dipoles employed in AMOEBA include the scaling of permanent-induced interactions (which are zeroed for nearest and second-nearest neighbors, as determined by bond connectivity), and the use of polarization groups (permanent multipoles do not induce dipoles within their own polarization groups). A detailed discussion of the AMOEBA electrostatics model is beyond the scope of this paper, and the interested reader is referred to Ref. 70.

The polarization energy $ E pol MM $ is given by the expression (A1) in the Appendix (Appendix 1), where we discuss polarization in more detail.

#### 3. Short-range valence interactions

By $ E val MM $ in (7) we denote all short-range valence interactions local to the MM subsystem. In the AMOEBA force field $ E val MM $ comprises the following terms: bond-stretch, angle-bend, stretch-bend coupling, out-of-plane bend, and torsional rotation. As the detailed expressions for these terms can be found in Ref. 73, Eqs. 2-6, we shall refrain from recounting them here.

#### 4. Van der Waals interactions

The term $ E vdW MM $ in (7) accounts for van der Waals (dispersion-repulsion) interactions local to the MM subsystem. AMOEBA uses the Halgren formulation^{85} of the buffered 14-7 potential,

where $ \rho i j = R i j / R i j 0 $, *δ* = 0.07, *γ* = 0.12. Detailed expressions, along with a description of mixing rules and hydrogen “reduction factors” can be found in Ref. 73, eqs. 7 and 8.

### D. Auxiliary point-multipole representation QM^{∗}

#### 1. Rationale

Like in most plane-wave and plane-wave-like DFT codes, the electronic density in onetep is represented on a uniform, real-space, Cartesian grid, with a typical spacing of ∼0.25 *a*_{0}, which corresponds to commonly used kinetic energy cutoffs in the range of 800 eV. As the localized orbitals in onetep typically extend for ∼8 *a*_{0} from the atomic cores, some MM sites will invariably be found in regions of non-negligible QM electronic density, arbitrarily close to (or even exactly at) grid points on which this density is represented. In the absence of any mitigating measures, such as damping, attempts to directly use the grid representation in calculations of polarization involving point dipoles run into issues of ill-conditioning due to the singularities of the point dipole potential and the sensitivity of energies to the positioning of the grid relative to the point dipole. These issues manifest as quantum-classical counterparts of the *polarization catastrophe*, such as unbounded polarization of the MM point dipole in the direct field of the (discretized) density, or the electronic density being “sucked out” by the unbounded potential of the MM point dipole.

In order to address these difficulties our model uses Thole damping in the calculation of QM/MM polarization, in line with what is done in AMOEBA. Since the Thole damping scheme is pairwise, it requires specifying the polarizabilities of two point dipoles, or in its more general version, of a point dipole and point multipole (cf. (30) and (31)). This requirement makes it impossible to directly use Thole damping for distributed electronic densities. The use of an auxiliary representation of the QM charge density in terms of point multipoles alleviates this problem.

#### 2. Implementation

Here we describe how to obtain an auxiliary representation of the QM subsystem charge density in terms of atom-centered, fixed, point multipoles up to a quadrupole. We shall refer to this representation as QM^{∗}.

The technique whereby a continuous charge distribution is represented in terms of a set of point multipoles is known as distributed multipole analysis (DMA). Atomic centers are usually, although not universally, used as the centers for the multipoles. DMA, first proposed by Rein,^{108} has been pioneered and popularised by Stone^{82} and Alderton.^{83} Distributed multipole analysis is typically performed in a Gaussian basis set.^{109,110} Below we briefly outline how electronic densities represented in a localised (NGWF) basis can be similarly expanded. A more detailed description can be found in Refs. 111 and 112. Our approach belongs to the class of the density-fitting techniques pioneered independently by Baerends *et al.*^{113} and Whitten.^{114}

The first step involves decomposing the electronic density (4) into two-center contributions

The shortcut notation *α* ∈ *I* used in (11) is understood as “NGWFs *α* belonging to atom *I*.” Equation (12) explicitly separates the density into two-center contributions from atomic centers *I* and *J* that have nonzero overlap *S _{IJ}*. The case of

*I*=

*J*, although technically a one-center contribution, can be treated on the same footing for simplicity of notation. By construction, NGWFs which do not overlap do not contribute to density.

We approximate each contribution as a linear combination of one-center contributions, represented in an auxiliary basis set, i.e.,

where $ C I J s $ are the sought coefficients in the expansion, and $ { f s ( r ) } s = 1 N f $ are the functions making up our auxiliary basis set, *N*_{f}/2 of which originate on center *I*, and the remaining *N*_{f}/2 on center *J*.

We subsequently define an electrostatic metric $V= V s t $,

where, for the sake of brevity, we introduced the notation

The use of the electrostatic (Coulomb) metric has long been recognized to be more advantageous compared to the simple overlap metric.^{115,116} Robust schemes that make it possible to maintain variationality in the presence of density fitting have been proposed.^{115–118}

The expansion coefficients are obtained as

by requiring that the electrostatic self-energy of $ n I J r \u2212 n \u0303 I J r $ be minimum. Here *t* indexes auxiliary basis functions originating on centers *I* and *J* in the same way *s* does. *V ^{ts}* are elements of the inverse electrostatic metric matrix.

When calculating gradients with respect to the density kernel, it is useful to separate the atom-pair coefficients $ C I J s $ into combinations of kernel matrix elements and NGWF-pair coefficients $ c \alpha \beta s $. By comparing (11) and (12) we obtain from (16)

We have now decomposed the electronic density of the system into a sum of atom-centered contributions, with each contribution being a linear combination of auxiliary basis functions. For the sake of concreteness, we will now explicitly assume the auxiliary basis functions to be truncated spherical waves (SWs),^{81,119} indexed by their angular momentum number *l*, magnetic quantum number *m*, and number of zeros *q* in the radial part., i.e.,

where *r*_{SW} is the localization radius of the SW, $ j l \u22c5 $ is a spherical Bessel function, and $ Z l m r \u02c6 $ is a real spherical harmonic. The values of the length scale *b _{lq}* are chosen so that $ j l r / b l q $ has exactly

*q*zeros on the interval (0,

*r*

_{SW}], with the last one at

*r*=

*r*

_{SW}for the truncation not to introduce a discontinuity. The index

*q*runs from 1 to a chosen maximum value

*q*

_{max}, which is in the order of 10–20 and in practice is limited by the fineness of the grid. Higher values make the density fitting more accurate at the expense of increased computational cost, which grows as $O N f 2 \u223cO q max 2 l 4 $. The localization radii

*r*

_{SW}of the SWs and of the NGWFs coincide.

We recall that the index *s* runs over auxiliary basis functions originating on *both* centers, while *l*, *m*, *q* describe a truncated spherical wave at **r** = **0**, and so

The sought spherical multipoles associated with an atomic center *I* can be calculated as^{111,120}

where $ C I J l m q $ correspond to $ C I J s $ originating *only* on *I*, and

is a radial integral that can be computed analytically.

While the approach described above is general, in this work we truncate the expansion at quadrupoles (*l* ≤ 2). Charges $ q I QM \u2217 $, Cartesian dipoles $ \mu I QM \u2217 $, and traceless Cartesian quadrupoles $ Q I QM \u2217 $ are obtained from the spherical representation (20) using well-known relations^{121} (eqs. 2.85-2.87). To obtain a representation of the total density, the charges *Z _{I}* of the atomic cores need to be added to $ q I QM \u2217 $.

The procedure described above minimizes the self-interaction energy of the error in the approximate density, with no constraint on the total charge of the entire system $ \u2211 I q I QM \u2217 $. Numerical tests indicate that the total electronic charge obtained from the expansion is within 0.05% of the expected number of valence electrons for reasonable qualities of the auxiliary basis set. Of course this valence charge is to a large degree compensated by a similar, but negative, contribution from the cores, which makes the relative error in the total charge substantially larger (and, infinite, by definition, for neutral systems). In practice we found it sufficient to compensate for this by uniformly rescaling the electronic monopole terms, matching the total with the number of valence electrons. We discuss the effect of this scaling on gradients in Subsection 3 of the Appendix.

### E. Consistent coupling between QM and MM

Our model uses the following total energy expression:

which we variationally minimize using gradient methods. The minimization is a single SCF process following the approach of Aida *et al.*,^{37} where the linear response equations for induced dipoles are solved at each SCF step. This iterative solution requires only the evaluation of classical multipole interactions, and its computational cost is minor compared to the costs associated with QM terms in the energy. To maintain variationality, this scheme requires the recomputation of the QM^{∗} representation at every SCF step. This can be done efficiently for two reasons. First, the bottleneck here is the evaluation of the electrostatic metric matrix (14), which does not depend on either the electronic or MM degrees of freedom, and can be precomputed. Second, by construction, the expansion coefficients $ c \alpha \beta s $ in (17) only need to be recomputed when the NGWFs change, and not at every SCF step.

In this communication we restrict ourselves to computing gradients with respect to the density kernel 𝕂, which we will derive for every term in turn. In future communications we plan to outline the calculation of gradients with respect to the NGWFs (which will permit optimizing them *in situ*), and with respect to ionic positions (which will pave the way towards geometry optimization, transition state searches, and molecular dynamics).

The superscript of each term identifies it as a property of a particular subsystem (QM, MM) or as a cross energy term (QM/MM). We shall discuss each term in turn. The single-system energy terms $ E QM $ and $ E MM $ are defined by (6) and (7), respectively. Regardless of whether the QM and MM subsystems are isolated or coupled, the functional form of these two terms is the same. That is to say, our approach does not modify the descriptions of the individual subsystems. When coupling is introduced, however, the interpretation of these terms changes: $ E QM $ becomes the energy of the QM subsystem polarized by the electric field of the MM subsystem, and $ E MM $ becomes the energy of the MM subsystem polarized by the electric field of the QM subsystem. To elucidate how these changes arise, we insert (6) and (7) into (22), to obtain the total energy expression for our model

The first five energy terms are insensitive to whether QM and MM are isolated or coupled, and, as they have been defined already, can be omitted from further discussion. The introduction of coupling does, however, change the values of $ E DFT QM $ and $ E pol MM $, which reflects the fact that the two subsystems mutually polarize.

#### 1. Polarization

The change in $ E DFT QM = E DFT QM n r $ is brought about by the deformation of the electronic density $n r $ in response to the electric field of the MM subsystem, and thus accounts for the polarization of QM due to MM. This deformation of $n r $ is driven by the gradient contributions (36) and (41).

The change in $ E pol MM $ is a consequence of the inclusion of the electric field of the QM subsystem in the direct field experienced by the MM induced dipoles, cf. (24). Following the introduction of coupling, $ E pol MM $ accounts for not only the internal polarization of the MM subsystem (cf. Fig. 1, interactions ⑤, ⑥), but also for the polarization of MM due to QM (cf. Fig. 1, interaction ④). The two contributions are non-additive. In our model polarization contributions from QM/MM interactions are damped, consistent with MM/MM polarization contributions owing to the use of the auxiliary QM^{∗} representation (Sec. II D 2). This representation is used in the calculation of $ E pol MM $ and its gradient $d E pol MM /d K \eta \theta $ (cf. (33)). This gradient contribution, apart from being essential for maintaining the variational behavior of our method, enables the electronic degrees of freedom to respond to the induced dipoles of the MM subsystem, thereby capturing the polarization effect of the environment.

Before we proceed with a more detailed account of polarization in our model, we refer the reader to Subsection 1 of the Appendix for a brief review of how polarization energy is calculated for an isolated MM system—this serves as a starting point for the discussion that follows.

In our QM/MM model the direct field experienced by an MM site *L* (which we will denote with $ E L \u2032 $) contains two contributions—the direct electric field of permanent MM multipoles (A5), and the direct electric field of QM^{∗} multipoles

The corresponding polarization energy is

which has the same form as (A6), except for the fact that the dipoles $ \mu L $ are now induced in response to the total direct field $ E L \u2032 $. Here, *N*_{MM} is the number of atoms in the MM subsystem.

The multipole expansion QM^{∗} is truncated at quadrupoles, leading to the following expression (cf. Ref. 121, eq. 2.63) for the electric field it generates at **R**_{L}:

The sum runs over all atomic sites *I* of the QM^{∗} representation. Each site contains a point charge $ q I QM \u2217 $, dipole $ \mu I QM \u2217 $, and quadrupole $ Q I QM \u2217 $.

The Thole-damped, dipole-charge interaction tensor $ T L I d-c $ is given by

the corresponding dipole-dipole interaction tensor is given by

and the dipole-quadrupole interaction tensor $ T L I d-q $ is given by

where *R* is the norm of the vector **R**_{I} − **R**_{L}, and $ R x , R y , R z $ are its Cartesian components.

Here, *λ*_{3}, *λ*_{5}, and *λ*_{7} are the Thole damping factors, given by

*u*is the reduced distance

*a* is a dimensionless width parameter,^{73} and *α _{L}*,

*α*are the (scalar) polarizabilities of the sites

_{I}*L*,

*I*, respectively. Since our model adopts a damped polarization treatment consistent with AMOEBA, it requires specifying classical scalar polarizabilities of QM sites,

*α*.

_{I}During the AMOEBA self-consistency procedure in which $ \mu L $ are determined, the degrees of freedom in QM^{∗} are temporarily clamped, meaning $ E L QM \u2217 $, and in turn $ E L \u2032 $, remain constant. A zero residual condition analogous to (A7) is satisfied at induced dipole self-consistency

As we follow a gradient-based approach, we seek the total derivative of the polarization energy with respect to a density kernel element *K*^{ηθ}:

Once the induced dipoles reach self-consistency, the second term vanishes owing to (32), leading to

which, with $ E L $ being independent of the density kernel, simplifies to

#### 2. Permanent electrostatics

The next term, $ E perm QM / MM $, accounts for the electrostatic interaction between the permanent multipoles in the MM subsystem and the QM charge density (cf. Fig. 1, interactions ①, ②). In our model this interaction is not damped, by analogy to how AMOEBA does not damp interactions between permanent multipoles. The Coulombic potential at **r** due to the permanent MM multipoles (cf. Ref. 121, eq. 2.62) is simply

where $ q L p $, $ \mu L p $, and $ Q L p $ have been defined in Sec. II C, and $ T L r c-c $, $ T L r c-d $, $ T L r c-q $ are, respectively, the Coulombic (undamped) charge-charge, charge-dipole, and charge-quadrupole interaction tensors acting from **R**_{L} to **r**,

*R*= |

**R**

_{L}−

**r**|.

The charge density of the QM subsystem is

and its energy of interaction with the potential of the permanent MM multipoles is given by

where Ω can be restricted to the union of the localization regions of $ \phi \alpha $, as $ n QM r $ is zero elsewhere.

In practice, $n r $ is defined on a uniform Cartesian grid, and the integral (40) is computed as a three-dimensional sum over grid points. Even though $n r $ is a slowly varying (pseudo)density, this approach suffers from similar ill-conditioning issues as those described in Sec. II D 1—since a grid point can be located arbitrarily close to a point MM multipole, the value of *R* in (38) can become arbitrarily small (or indeed zero), making the corresponding tensors, and the resultant potential (37) unbounded, and $ E perm QM / MM $ ill-conditioned.

One solution to this problem would be to resort again to the QM^{∗} representation and to calculate this interaction energy as a pairwise sum of point-multipole–point-multipole interactions. However, abandoning the distributed description of charge density in favor of a multipole expansion would lead to the introduction of charge penetration error. This well-known deficiency of point-multipole models^{122} is a consequence of their poor description of the interaction between extended atomic charge densities at short distances. In this regime the overlap of the two atomic densities becomes significant, leading to a decrease in the shielding of the nuclear charge by its density,^{123} an effect that is not captured by point-multipole models, unless specifically corrected for.

In order to preserve the advantages that having access to the full density $ n QM r $ offers, we choose to avoid the ill-conditioning by smoothing the potential (37) in the vicinity of every MM multipole. To this effect we replace Coulombic interaction tensors (38) with their Thole-damped counterparts, but instead of using atomic polarizabilities in the denominator of (31), we use a fixed value of 0.2 *a*_{0} for the characteristic length. The influence of such smearing is negligible already for *R* > 0.5 *a*_{0}, but it effectively removes singularities as *R* → 0. We find the results to be practically insensitive to the particular choice of this value.

Since $ v p MM r $ does not depend on the electronic degrees of freedom, the corresponding energy gradient is simply a matrix element of the potential in the NGWF basis

This gradient term enables the electronic degrees of freedom to respond to the permanent dipoles of the MM subsystem, thereby capturing the direct polarization of QM by MM.

#### 3. Van der Waals interactions

In the calculation of $ E vdW QM / MM $ we use the same formalism as is used for MM/MM van der Waals interactions, i.e., the pairwise Halgren formulation^{85} of the buffered 14-7 potential.^{73} The calculation of this term is performed entirely in tinker, with the QM subsystem treated as an embedding *inactive region*, which avoids calculating QM/QM contributions that are already accounted for in $ E QM $.

Adopting such classical approach requires choosing suitable vdW parameters for atoms in the QM subsystem. In contrast to the fine-grained system of atom types that AMOEBA uses for electrostatics, van der Waals interactions are parametrized using a broader notion of *atom classes*, meaning the vdW parameters are, to a large extent, shared by atoms whose chemical species and hybridization are identical. With this in mind, we simply use unmodified AMOEBA parameters for atom classes deemed to be nearest matches for atoms in the QM subsystem, avoiding the need for parametrizing the QM subsystem altogether.

Since $ E vdW QM / MM $ is independent of the electronic degrees of freedom, the corresponding gradient contribution vanishes

### F. The TINKTEP implementation

tinktep is a software package that enables self-consistent, mutually polarizable QM/MM calculations combining onetep and tinker. The theory outlined above has been implemented entirely within onetep, within an infrastructure of general-use software modules: spherical wave resolution of identity (SWRI), spherical wave expansion (SWX), distributed multipole analysis (DMA), and polarizable embedding (PE).

Our approach does not require any intervention into the general mechanisms of tinker in order to actualize the coupling. Only a small, well-contained set of trivial adjustments to tinker’s input/output routines is necessary to increase the numerical precision of certain inputs and outputs, and to adjust the maximum line length. This is realized through a patch that is distributed with onetep.

tinktep itself comprises several core scripts and about a dozen utility scripts. The role of the core scripts is to provide a user-friendly environment in which onetep and tinker are executed, and to oversee and synchronize their execution, providing user feedback, error logging, and graceful abort in case of errors. The utility scripts facilitate the set-up of QM/MM calculations, file format conversions, translations of frame of reference etc.

tinktep is available to all users of onetep v4.5 and later. In addition to polarizable QM/MM calculations, tinktep supports QM/MM calculations with fixed point-charge force fields, such as the General Amber Force Field (GAFF).

## III. DEMONSTRATION OF METHODOLOGY

We demonstrate the feasibility of our QM/MM approach using two test setups. First, we briefly validate the correctness of the total (QM+MM) energy gradient by examining the convergence of total energy for a diphenylhydramine molecule embedded in a water sphere. Second, we closely examine how accurately the binding energy of six solutes to water solvent is reproduced by a number of MM and QM/MM approaches. Having selected larger, small, neutral, and charged solutes, and investigating the binding energy curves for increasing sizes of solvent spheres, we can elucidate the advantages and deficiencies of each approach. We demonstrate that in terms of accuracy our approach outperforms QM/MM schemes that do not take mutual polarization into account.

### A. Correctness of total energy gradient

We begin with basic validation of the correctness of our implementation. We compared total (QM+MM) energy gradients obtained by finite differences (FD) with analytical gradients, ensuring they were in agreement for a variety of systems, including charged systems, where mutual polarization becomes significant. In all cases we were able to converge the embedded system to the same thresholds as the purely QM system, and we obtained expected agreement between analytical and FD gradients. Below (Fig. 2) we restrict ourselves to demonstrating that the convergence of the total energy during density matrix optimization proceeds very similarly regardless of the type of embedding used. We demonstrate this on a representative system, where the QM subsystem comprises a diphenylhydramine molecule with 75 surrounding H_{2}O molecules, and the MM subsystem comprises 256 H_{2}O molecules.

### B. Binding energies

To assess the accuracy and robustness of our QM/MM approach we employed it to calculate the interaction energy of six solutes with progressively larger shells of explicit water molecules. Three of the solutes were chosen from the SAMPL4 blind challenge^{124}—these were (a) (–)-menthol, (b) diphenylhydramine, and (c) 2-chloro-4-hydroxy-3,5-dimethoxybenzaldehyde. These moderately sized molecules (31, 40, and 23 atoms, respectively) encompass a number of chemical features: a cyclohexane ring (a), an ether group (b), an aromatic ring (b), an amine group (b), a halogen atom (c), and an aldehyde group (c). The remaining three molecules were (d) ammonia (NH_{3}), (e) the ammonium ion ($ NH 4 + $), and (f) the cyanide ion (CN^{−})—which we chose to verify if our model correctly describes small and charged solutes.

#### 1. Computational set-up

Each of the six solutes has been solvated in approx. 660 explicit H_{2}O molecules under periodic boundary conditions (PBCs). Classical polarizable MD trajectories in the *NpT* ensemble (*p* = 1 atm, *T* = 298 K) were then obtained using the dynamic program from the tinker suite. The dynamics was run for 50 ps with a timestep of 1 fs and the final configuration was used in subsequent calculations. AMOEBA parametrization for solutes (a)-(c) has been taken from Ref. 96, solutes (d)-(f) used parameters natively available in AMOEBA09^{72} and for water molecules we used the AMOEBA 2003 water model. Long-range electrostatics used the particle mesh Ewald approach,^{125} with a real-space cutoff of 8 Å. van der Waals interactions were cut off at 9 Å, and a long-range correction was applied.

The above procedure yields configurations suitable for calculations with periodic boundary conditions (cf. Fig. 3, panel (a)). At present our QM/MM implementation only supports calculations with open boundary conditions, which necessitates augmenting the simulation protocol with an intermediate step. First, we discard all water molecules which are not entirely contained in a sphere with a radius of half of the simulation box size (cf. Fig. 3, panel (b)). No longer part of the bulk, the water molecules at the surface of the sphere are now misoriented, generating a non-negligible, spurious dipole moment. As has been recently shown by Lever *et al.*,^{11} this can lead to an unphysical lowering of the HOMO-LUMO gap in *ab initio* calculations, making it crucial to address this effect. We mitigate the problem by allowing the outermost water molecules to relax to a local energy minimum. This is illustrated in Fig. 3, panel (c). The relaxation is performed with tinker’s optimize program, using open boundary conditions (OBC) and an infinite interaction cutoff. The solute and water molecules whose geometric centers are within 12 Å from the geometric center of the solute are restrained during the relaxation. As expected, in the course of the relaxation we observe a lowering of the dipole moment of the system and the final solute-solvent binding energy is close to that of the original PBC system. The resultant configurations consist of the solute and 325-400 H_{2}O molecules and are suitable for calculations with OBC.

We examined the binding energies between the solute and the surrounding water molecules for a single snapshot for each of the solutes, investigating how the solute-solvent interaction energy converged as the number of surrounding H_{2}O molecules was increased. Water molecules were added in the order of increasing distance from the center of the solute. Binding energies were calculated as

with atomic positions taken from the solute+solvent configuration.

Our comparison involved four computational approaches:

Fully QM calculations with no embedding (entire system treated at the DFT level of theory), which serve as reference.

QM calculations using a purely electrostatic embedding, where

**the QM subsystem encompassed only the solute**, and H_{2}O molecules were described with fixed partial charges. In this set-up only a fixed, external potential is included in the QM Hamiltonian; we emphasize the neglect of van der Waals interactions between the QM and the embedding.QM/MM calculations with either a fixed point-charge embedding (GAFF v1.5

^{126}) or a polarizable embedding (AMOEBA). Here too**the QM subsystem encompassed only the solute**, and all water molecules were described by a classical force field. What is also different from (b) is that van der Waals interactions between the solvent and solute were included at the MM level of theory (Lennard-Jones potential for GAFF embedding, Halgren’s 7-14 potential (cf. (9)) for AMOEBA embedding).Fully MM calculations, where the entire system was treated with classical MM (GAFF or AMOEBA).

All QM calculations were performed using the onetep linear-scaling package. Since our QM/MM approach does not currently support optimizing the NGWFs *in situ*, we used a fixed minimal NGWF basis pre-optimized in vacuum, and only optimized the density kernel 𝕂. Open boundary conditions were imposed by using direct Coulombic summation (5) for $ E core - core QM $, and through the use of the cut-off Coulomb^{127} technique in the calculation of electronic interactions. The kinetic energy cutoff was set to 1000 eV. The PBE^{128} exchange-correlation functional was used. Empirical dispersion correction in the formulation of Elstner,^{95} with parameters determined by Hill and Skylaris^{86} was employed to correct the deficiencies of GGA DFT description of dispersion interactions.

All MM calculations were performed with tinker’s analyze program. Open boundary conditions were used, with an infinite cutoff for all interactions. Induced dipoles were converged to a tight rms threshold of 10^{−11} D in order to make the residual (A7) negligible.

QM/MM calculations used the above settings for the QM and MM subsystems, respectively. QM/MM coupling was effected through the use of the tinker-onetep interface, tinktep (cf. Sec. II F), that has been implemented as part of this work. In QM/polarizable-MM calculations the QM^{∗} expansion (cf. Sec. II D 1) was truncated after quadrupoles. In fixed point charge QM/MM calculations and in QM calculations with fixed point-charge embedding we used partial charges of 0.417 *e* for H atoms and −0.834 *e* for O atoms, which are identical to the TIP3P^{129} model used in GAFF. In this case the QM^{∗} expansion was limited to atom-centered charges only—in the absence of polarization, the QM^{∗} representation is only of relevance to the tests of the refined model described in Sec. III B 3, and in Fig. 6, curve (b) in particular.

#### 2. Initial results

We first focus on one of the molecules, (–)-menthol, to illustrate in detail how the accuracy of MM, QM, and QM/MM approaches compares when no adjustments are made to any of the models, and to elaborate on the metrics we used for judging accuracy.

Figure 4 shows the binding energy between the (–)-menthol molecule and the molecules of the solvent. By the time the water shell comprises 130 H_{2}O molecules, the reference QM energy is converged to within 1 kcal/mol of the value obtained for the largest shell. The purely classical AMOEBA results closely track the QM trend, consistently underbinding by as little as 2 kcal/mol for the larger systems. GAFF’s predictions are not as accurate, with the error accumulating quickly in the short range and plateauing at about 10 kcal/mol by the time 50 H_{2}O molecules are reached.

The behaviors of the studied approaches are easier to interpret if we adopt the QM results as a baseline, and examine *errors* in the binding energy, understood as energy differences from the reference. To this effect we replot the same data in Figure 5, immediately re-establishing that AMOEBA’s predictions are very good for all system sizes, and that GAFF’s long-range behavior is correct, but its predictions are plagued by short-range error.

We now examine the predictions of QM/MM models and the of the purely electrostatic embedding approach. Electrostatic (point-charge) embedding severely underbinds the system by accumulating as much as 20.5 kcal/mol of error in the short- and medium-range. To a large extent this is an expected consequence of not taking solute-solvent (QM/MM) van der Waals interactions into account. Since the partial charges used in calculations with electrostatic embedding and QM/MM with GAFF are identical, and the van der Waals contribution in the latter is calculated classically and so does not affect the electronic degrees of freedom, the entire difference between electrostatic embedding and QM/MM with GAFF is due to missing van der Waals interactions. Here these interactions are strongly attractive and their neglect accounts for the majority (12.5 kcal/mol) of electrostatic embedding’s error (cf. Fig. 5(a)). The remaining error is the same as the error in QM/MM with GAFF and amounts to about 8 kcal/mol in the long range (cf. Fig. 5(b)).

In the calculation of binding energies, all intra-MM terms (valence, van der Waals, and electrostatic) cancel out between the solute-solvent “complex” and the solvent-only calculation. The classical intra-QM dispersion correction similarly cancels out between calculations on the complex and on the solute. The only remaining energy terms are QM/MM electrostatic and van der Waals terms, and intra-QM electronic energies (i.e., the polarization response of the solute). Of these three, only point-charge electrostatic interactions are long-ranged, and the behavior of the QM/MM GAFF curve—which becomes almost flat beyond 100 H_{2}O molecules—suggests that in this case the point-charge description is sufficient at long range. We will later demonstrate that this is no longer the case for charged solutes. Here, however, the errors in both MM GAFF and QM/MM with GAFF are already accumulated (to values of ∼10 kcal/mol and ∼8 kcal/mol, respectively) by the time the solvent shell comprises ∼60 H_{2}O molecules—a regime where the mean distance between a solute atom and a solvent atom is 6.8 Å. This indicates the unsurprising breakdown of the fixed point charge description, and/or GAFF’s van der Waals model, at short range, where AMOEBA’s polarizable model is seen to cope very well. The curves obtained from MM with GAFF and QM/MM with GAFF have very similar shapes, differing mostly by an almost constant offset of 2.5 kcal/mol. This implies that the majority of the error (the remaining ∼8 kcal/mol) is due to a short-range deficiency in the shared components of the two models, i.e., in GAFF’s treatment of solute-solvent vdW interactions and the limitations of point-charge description—rather than any serious deficiency of the QM/MM interface. Aida *et al.* reached similar conclusions for solute-solvent interactions with their QM/MM-pol-vib model.^{37}

The opposite is true for the QM/MM with AMOEBA model proposed in this work. A comparison of the QM/MM with AMOEBA curve with the MM AMOEBA curve reveals that their shapes differ significantly only at short range (up to ∼40 H_{2}O molecules), where the QM/MM calculation already underbinds by 8.5 kcal/mol (cf. Fig. 5(c)) compared to MM AMOEBA. From this point on the offset between the two curves remains practically constant, indicating almost identical, correct medium-, and long-range behavior. Given that at short range the purely MM AMOEBA description is in remarkable agreement with the DFT result (cf. Fig. 5(d)), we must conclude that in this case the QM/MM interface itself is responsible for most of the error by which our QM/MM scheme underbinds (–)-menthol.

#### 3. Model refinement

We will now elucidate the sources of this error, and propose a simple measure to remediate the problem. The main differences between the treatment of interactions in QM/MM with AMOEBA and MM with AMOEBA are the following:

Quantum-mechanical treatment of the QM subsystem.

Absence of QM-side charge penetration error in the treatment of QM/permanent-MM electrostatics (cf. (40)).

Inclusion of polarization contributions in $ n QM r $ appearing in (40) in the interaction of QM with permanent MM multipoles. The entire QM contribution is thus not damped, whereas it would be more consistent with AMOEBA methodology to separate this term into a permanent part ($ n vac QM r $) and an induced part ($ n QM r \u2212 n vac QM r $), damping the latter.

Of these three differences, the first two are, of course, desired features of the model, introduced to improve upon the MM description. The last one is an unwelcome simplification resulting from the desire to avoid dealing with induced partial charges and induced quadrupoles in the QM^{∗} representation, as these would not be directly compatible with the AMOEBA model, and having to perform a separate QM calculation in vacuum.

Regardless of their origin, each of these changes affects the behavior of electrostatic interactions in the system, with no corresponding change in the repulsive part of van der Waals interactions that AMOEBA leverages to balance strongly attractive electrostatics at short range. Similarly, since our approach changes the description of permanent interactions, while retaining AMOEBA’s damped multipolar scheme for polarization, we potentially disrupt the balance between permanent and induction interactions. That is to say, even when the description of electrostatics is improved, e.g., by the elimination of charge penetration error on the QM side, the balance between electrostatics and van der Waals interactions can easily become disrupted, necessitating adjustments to parametrization or functional forms for interactions crossing the QM/MM interface. Similar conclusions have been reached independently by Aida, Yamataka, and Dupuis.^{37} Carnimeo *et al.*^{55} have also found it necessary to adjust van der Waals parameters in the QM/MM interface. Ultimately, the classical description of van der Waals interactions, with its inherent neglect of the coupling to an atom’s local electronic structure,^{130} may altogether prove unsatisfactory for describing QM/MM interactions. More refined models, where van der Waals parameters could be made density-dependent through atomic volumes,^{131} or where exchange and dispersion energies are made charge-dependent in a many-body formulation^{130} might alleviate this issue in the long term.

We now lay out a very simple two-step remediation measure that we find to be sufficient to obtain reasonable accuracy for both neutral and charged solutes.

In the first step we resign ourselves to approximately reintroducing charge penetration error (CPE) into QM/permanent-MM interactions in the hope of restoring some of the balance between electrostatics and the repulsive van der Waals term. We achieve that by “correcting” the converged total energy—removing the CPE-mitigated interaction (40) and replacing it with the CPE-afflicted interaction energy of QM^{∗} with the permanent MM multipoles, as calculated by tinker. This is not equivalent to reintroducing the entirety of CPE—we emphasize that the QM density is still optimized under the original CPE-mitigated Hamiltonian, and the “correction” is done *a posteriori*. Furthermore, the CPE-mitigated interaction only accounted for charge penetration on the QM side, with the MM subsystem still represented by point multipoles. Nevertheless, this improves the agreement of our model with the QM reference (cf. Fig. 6, (a)), not only for (–)-menthol, but for all six studied molecules. Interestingly, analogous reintroduction of CPE into QM/MM with GAFF leads to marked worsening of obtained binding energies, again for all six molecules (shown for (–)-menthol in Fig. 6, (b)). It remains to be determined whether the difference in behaviors is due to the multipolar nature of AMOEBA, the fact that it is a polarizable model (and thus has to balance permanent electrostatics and polarization), or the fact that its balance between electrostatics and van der Waals interactions is perhaps somewhat more fragile. Intuitively we would expect switching back to the multipolar description to bring the QM/MM description closer to the MM model, and this is indeed what happens—binding becomes stronger in QM/MM with AMOEBA (compare Fig. 6, (a) with Fig. 5, blue and green curves), and weaker in the case of QM/MM with GAFF (compare Fig. 6, (b) with Fig. 5, orange and grey curves). A detailed study of the electrostatics of our model, supported by energy decomposition analysis, will be the subject of a future communication. In the following we will of course *not* include this detrimental adjustment in the results obtained with QM/MM with GAFF.

The second step of our refinement consists of a simple adjustment of the steepness of the repulsive wall in the buffered 7-14 van der Waals potential (9) used in AMOEBA. We stress that this is only done for interactions crossing the QM/MM interface. Moreover, this adjustment does not involve any changes to the position or value of the minimum of the potential, and is independent of chemical species. That is, we keep using the original parametrization of AMOEBA, only applying a uniform *a posteriori* adjustment to the shape of the repulsive wall. Our change amounts in replacing the value of *δ* = 0.07 in (9) with *δ* = 0.21, which achieves the effect illustrated in Fig. 7. This adjustment is done with the aim of attenuating a small number of severely repulsive van der Waals interaction pairs corresponding to a scenario where a QM atom and an MM atom are drawn close together by very favorable electrostatic interactions, as is the case for hydrogen bonds crossing the QM/MM interface. In AMOEBA the energetics of such pairs relies on a delicate balance between strongly attractive electrostatic interactions and strongly repulsive van der Waals interactions. This balance, disrupted by the differences in the treatment of electrostatics in our QM/MM model, can be restored, to a large degree, by attenuating only the most excessive van der Waals interactions. The significant improvement to the short-range behavior of our QM/MM model resulting from this simple adjustment can be seen in Fig. 6, (c). Each of the remaining solutes, with the sole exception of $ NH 4 + $, benefits from this final adjustment to our model. We point out that a similar adjustment is neither possible, nor necessary for QM/MM with GAFF. This is because the Lennard-Jones potential does not offer the flexibility of the buffered 7-14 potential, and, respectively, because in GAFF the balance in question can be more robust owing to the absence of polarization. The choice of *δ* = 0.21 minimises the extremely short range error (best examined, for the (–)-menthol test case, by comparing curves (a) and (c) in Fig. 6) across the six tested molecules, while ensuring long-range behavior is not adversely affected—as excessive changes to *δ* will gradually influence the attractive regime.

#### 4. Final results

The performance of each of the approaches is compared in Table I, where we report errors of each approach (with DFT results as baseline) obtained for the systems with the largest numbers of solvent molecules. Error values were averaged over systems with 200 or more H_{2}O molecules to give an idea of the feasibility of each approach in practical scenarios, where the computational effort associated with the MM embedding is small compared with the effort of the QM calculation, permitting the use of large embedding regions.

. | MM . | MM . | QM/MM . | QM/MM . | QM/MM . | QM/MM . |
---|---|---|---|---|---|---|

Molecule . | GAFF . | AMOEBA . | point-charge . | GAFF . | AMOEBA (initial) . | AMOEBA (refined) . |

(–)-menthol | 10.4 | 1.9 | 20.5 | 8.0 | 10.5 | −1.3 |

diphenylhydramine | 13.7 | −0.1 | 41.3 | 14.9 | 15.3 | 2.8 |

2-Cl-4-OH-3,5-dimethoxy-BALD | 4.7 | −0.6 | 24.4 | 6.8 | 6.9 | 0.3 |

NH_{3} | 2.7 | 1.5 | 3.7 | 4.1 | 7.6 | −1.8 |

$ NH 4 + $ | 5.1 | 1.1 | −5.0 | 4.5 | 1.2 | −4.1 |

CN^{−} | −7.7 | 7.2 | 0.9 | −1.8 | 18.2 | 2.9 |

RMSE | 8.3 | 3.1 | 21.5 | 7.9 | 11.4 | 2.5 |

MSE | 4.8 | 1.8 | 14.3 | 6.1 | 10.0 | −0.2 |

. | MM . | MM . | QM/MM . | QM/MM . | QM/MM . | QM/MM . |
---|---|---|---|---|---|---|

Molecule . | GAFF . | AMOEBA . | point-charge . | GAFF . | AMOEBA (initial) . | AMOEBA (refined) . |

(–)-menthol | 10.4 | 1.9 | 20.5 | 8.0 | 10.5 | −1.3 |

diphenylhydramine | 13.7 | −0.1 | 41.3 | 14.9 | 15.3 | 2.8 |

2-Cl-4-OH-3,5-dimethoxy-BALD | 4.7 | −0.6 | 24.4 | 6.8 | 6.9 | 0.3 |

NH_{3} | 2.7 | 1.5 | 3.7 | 4.1 | 7.6 | −1.8 |

$ NH 4 + $ | 5.1 | 1.1 | −5.0 | 4.5 | 1.2 | −4.1 |

CN^{−} | −7.7 | 7.2 | 0.9 | −1.8 | 18.2 | 2.9 |

RMSE | 8.3 | 3.1 | 21.5 | 7.9 | 11.4 | 2.5 |

MSE | 4.8 | 1.8 | 14.3 | 6.1 | 10.0 | −0.2 |

For each of the molecules, purely classical calculations with AMOEBA are in significantly better agreement with DFT than GAFF is, with CN^{−} being the only solute that is not captured to within 2 kcal/mol, presumably due to the complex electronic structure of this ion (a *σ*-donor, *π*-acceptor ligand). Calculations with point-charge embedding perform very poorly due to their neglect of solute-solvent van der Waals interactions. With the exception of $ NH 4 + $ and NH_{3}, these interactions are attractive, which means their neglect leads to consistent underbinding, which is particularly severe for larger molecules—this is well evidenced in the large, positive mean average error (14.3 kcal/mol) plaguing this approach. The inclusion of van der Waals interactions at the GAFF level of theory (QM/MM with GAFF embedding) improves accuracy, but only to a level that is on par with purely MM GAFF calculations. As pointed out earlier, we believe the errors here are mostly attributable to the fixed point charge description, rather than to an unbalanced QM/MM interface.

Initial results obtained with QM/MM with AMOEBA indicate that the QM/MM interface becomes unbalanced by the mitigation of charge penetration error and other changes that our model applies to the treatment of electrostatics. As a result, the initial model consistently underbinds due to strongly repulsive van der Waals terms that are no longer sufficiently balanced by strong electrostatic attraction. This effect manifests most significantly for hydrogen bonds crossing the QM/MM interface. Once the model is adjusted by re-introducing some of the charge penetration error and softening the repulsive wall of the buffered 7-14 potential, we observe a dramatic improvement. While our refined model cannot yet boast chemical accuracy, its rms error is the lowest of all the studied approaches, and the very low mean average error (−0.2 kcal/mol) is a good indication of a well-balanced approach. Given that in this work only the solute molecules are treated at the DFT level of theory, we consider achieving an rms error 2.5 kcal/mol against full DFT calculations on 1000+ atoms a success.

Full binding energy error curves illustrating how the studied approaches compare at all system sizes are shown in Fig. 8. For QM/MM with AMOEBA only the results obtained with the refined model are shown for clarity. The conclusions we have reached using (–)-menthol as a test case are seen to be generally applicable to the remaining molecules. Below we briefly comment on a number of differences observed for particular molecules.

Earlier we used the similarity in the trends between QM/MM with GAFF and MM GAFF and to argue that the majority of the error in both approaches is likely due to a shared component of the models (GAFF’s description of solute-solvent interactions). If so, we would expect the short-range part of the two curves to differ more for solutes that are more difficult to describe with a pure MM description (where the difference between the QM and MM treatments of the solute would be highlighted). This is indeed the case for 2-chloro-4-hydroxy-3,5-dimethoxybenzaldehyde.

For (–)-menthol we reasoned about good long-range behavior of MM GAFF from the flat shape of the energy error curve at long range. We point out that this is no longer the case for charged solutes, for which the fixed-point charge model is clearly insufficient. This is evidenced by large oscillations in the MM GAFF, QM/MM with GAFF and QM with point charge embedding curves for $ NH 4 + $ and CN^{−}. The same curves are much flatter for MM AMOEBA and QM/MM with AMOEBA, indicating that, as expected, polarization of the solvent needs to be taken into account for charged solutes.

Finally, we point out that the largest errors in our QM/MM model also appear for charged solutes, and have opposite signs for a cation and an anion. This suggests directions for further refinement of the model through improving the description of electrostatics, and polarization in particular.

## IV. CONCLUSIONS AND CLOSING REMARKS

We have presented an implementation of a QM/MM approach in which the quantum subsystem described by DFT is coupled to a classical subsystem described by the AMOEBA polarizable force field. The two components mutually polarize one another within a total energy minimization scheme which achieves self-consistency for both the MM and QM subsystems. We have derived an expression for the Hamiltonian of the coupled QM/MM system, which we minimize using gradient methods.

We describe the QM subsystem with the onetep linear-scaling DFT program, which makes use of localized orbitals expressed in a set of periodic sinc basis functions equivalent to plane waves. We have interfaced onetep with the tinker code, which is a prototypical implementation of the AMOEBA force field, used in our model to describe the MM subsystem. We have put great emphasis on treating polarization interactions consistently between the MM and QM subsystems, particularly with regard to damping, which is a crucial mechanism in polarizable point dipole approaches.

We have carried out tests to validate our method, demonstrating the simultaneous optimization of the quantum and classical degrees of freedom. We identified and remediated the sources of inaccuracy in the QM/MM interface that stem from a disruption of the balance between (improved) electrostatics and van der Waals interactions.

This is a proof-of-principle implementation, as we have not yet implemented the *in situ* optimization of the local orbitals of onetep. Despite this limitation, which we plan to address in a future communication, our results indicate that our approach offers superior convergence and accuracy compared to conventional QM/MM methods. Future work also will be devoted to investigating suitable reparametrization of interactions crossing the QM/MM interface, and refining the treatment of electrostatics in our model.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the coordinates of solvated test systems referred to in Figs. 2, 4–6, and 8.

## Acknowledgments

J.D. and C.K.S. acknowledge the support of the Engineering and Physical Sciences Research Council (EPSRC Grant Nos. EP/J015059/1 and EP/K039156/1), and of the UKCP consortium (EPSRC Grant No. EP/K013556/1) for access to the ARCHER national supercomputer and the IRIDIS High Performance Computing Facility of the University of Southampton. T.H.G., M.H.G., and Y.M. acknowledge the support by Grant No. CHE-1363320 from the U.S. National Science Foundation. JWP would like to acknowledge NSF Award No. CHE1265712. Some of the calculations in this work were carried out at the TASK Computer Centre in Gdansk, Poland (supported in part by the PL-Grid Infrastructure, Grant No. POIG.02.03.00-00-096/10).

The authors are grateful to Alex Albaugh, Richard Bradshaw, Bernard Brooks, David Case, Jon Essex, David O’Regan, Jean-Philip Piquemal, and Joshua Rackers for fruitful discussions. Parametrizations and initial structures of solutes used in this work were kindly provided by Richard Bradshaw.

### APPENDIX: POLARIZATION ENERGY AND ITS DERIVATIVE WITH RESPECT TO THE DENSITY KERNEL

#### 1. Polarization in an isolated MM system

We start from the expression for the polarization energy of a set of point dipoles as formulated by Simmonett *et al.*^{100} (eq. 9), which we rewrite using explicit summations over *N*_{MM} polarizable dipole sites

where **E**_{L} is the *direct* electric field at site *L*, $ \mu L $ is the dipole induced at site *L* in response to the total (direct and mutual) electric field, and **T**_{LM} is a 3 × 3 coupling tensor between sites *L* and *M*,

$ T L M d-d $ is the Thole-damped, Cartesian dipole-dipole interaction tensor (dipole field tensor) between induced dipoles at sites *L* and *M* (cf. (28)). Authors preferring the block matrix notation refer to the 3*N* × 3*N* block matrix $T= T L M $ describing the entire system as the coupling tensor^{100} or the relay matrix.^{106}

We identify the three terms in (A4) as: the work cost of assembling the set of dipoles, the mutual interaction energy of induced dipoles, and the interaction of induced dipoles with the direct field.

By expressing the direct field as a difference between the total field and the mutual field

we can rewrite (A3) simply as

which is the same as the result given in Ref. 72, eq. 5.

In the standard AMOEBA formulation the induced dipoles $ \mu L $ are determined through an iterative procedure, with the zero residual condition satisfied at selfconsistency^{100}

#### 2. Gradients of the auxiliary representationQM^{∗} with respect to the density kernel

Here we finalize the derivation presented in Sec. II E, where we sought to calculate the gradient of the polarization energy with respect to the density kernel. We continue from Sec. II E, where the components that remain to be derived are the derivatives of Cartesian traceless multipole moments with respect to the density kernel, $ \u2202 q I QM \u2217 \u2202 K \eta \theta $, $ \u2202 Q I QM \u2217 \u2202 K \eta \theta $, and $ \u2202 \mu I QM \u2217 \u2202 K \eta \theta $.

where

The shortcut notation *α* ∈ *I* used in (A8) is understood as “NGWFs *α* belonging to atom *I*.” Since $ d \alpha \beta l m $ are independent of the density kernel, it follows from (A8) that

This means that a multipole at a QM atomic site *I* only has non-vanishing derivatives with respect to density kernel elements coupling NGWFs on atom *I* to NGWFs on its neighbors *J*. By neighbors we mean atoms *J* (including *J* = *I*), whose overlap with *I* is non-zero, *S _{IJ}*≠0.

With this in place, the only operation left is the conversion between spherical multipoles and Cartesian traceless multipoles

#### 3. Charge scaling

In Sec. II D 2 we briefly mentioned that the total charge obtained from the DMA procedure is not constrained to be an integer, and typically needs to be scaled to match the number of valence electrons. This is achieved by replacing all electronic monopoles $ q I QM \u2217 $ with $ q \u0303 I QM \u2217 =\lambda q I QM \u2217 $, where *λ* is the ratio between the expected valence charge and the total charge obtained from DMA, i.e.,

The scaling is very modest, with typical values of $\lambda \u2208 0 . 9995 , 1 . 0005 $. Nevertheless, this additional dependence on 𝕂 needs to be taken into account in the gradients. To this effect we correspondingly replace (A11) with

where we used