Mutually polarizable QM/MM model with in situ optimized localized basis functions

Molecular dynamics (MD) is a well-established technique for simulating the structure and properties of systems at the atomic scale, with over four decades of applications in biochemistry and materials science, among other fields. The aim of MD is to predict macroscopic behavior from microscopic interactions,¹ and the validity of results strongly depends on how accurately these interactions are described by the molecular mechanics (MM) potential.

The continual increase in available computational power not only extends the scope of MD to larger systems and longer time scales but also enables the refinement of MM models describing inter- and intra-molecular interactions. The last two decades have witnessed the emergence of force fields that directly capture many-body polarization effects, setting out to circumvent well-known deficiencies of pairwise-additive, fixed point charge models.^2–5 Unable to directly account for polarization, fixed point charge force fields struggle to describe, e.g., the interactions of ions with π-electron systems or polar solutes in low-dielectric media,⁶ and they are typically poorly transferable to environments or phases different from those that they were parameterized for, such as interfaces.⁷ 

The non-additive, many-body nature of polarization interactions makes polarizable models more involved and computationally demanding. Consequently, a variety of competing treatments of polarization exists (see Refs. 7–9 for a review): Drude oscillators,^10,11 fluctuating charges,^12,13 induced point dipoles,^14–21 or even induced multipoles of higher order.²² The AMOEBA force field,^15–18 which is of particular significance to this work, describes polarization interactions using damped, induced, point dipoles, while for permanent electrostatics, it employs fixed multipoles up to a quadrupole in lieu of point charges.

Purely classical models, however sophisticated, cannot describe electronic properties, such as bandgaps or solvent shifts, or processes that intrinsically depend on the electronic degrees of freedom, such as bond-breaking. To properly model electronic phenomena, it becomes necessary to employ quantum-mechanical (QM) methods. In practical applications, density functional theory (DFT) is arguably the most commonly used approach,²³ owing to its relatively low computational cost. Even so, length scales (∼10–100 nm) and timescales (∼1 μs) typically used in classical MD simulations remain beyond the scope of DFT today.

QM/MM combines the quantum and classical descriptions, exploiting the fact that the properties of interest are often localized to a part of the system that can be described quantum mechanically, such as a molecule, embedded in an environment that can be described more approximately with MM, e.g., a solvent. Since the seminal work of Warshel and Levitt,²⁴ a profusion of QM/MM approaches has been proposed, targeting different types of systems and varying in the level of sophistication (see, e.g., Refs. 25–36). Even a brief review of QM/MM methods is beyond the scope of this paper; however, we refer interested readers to reviews of QM/MM methods and applications in enzymology,³⁷ biochemistry,³⁸ and materials science.³⁹ 

Recent improvements in force fields promptly become integrated into QM/MM methodologies, and several approaches combining QM with polarizable forcefields (dubbed QM/MMpol) have already been proposed. Many of those approaches employ induced dipoles^40–55 to model polarization, and others use fluctuating charges^56–59 or Drude oscillators.^60,61 Typically, the full QM density is used for all electrostatic QM/MM interactions,^5,52,55 but using auxiliary multipolar representations of QM (for efficiency or convenience) has also been proposed.^51,62 Several groups have developed models specifically focused on electronic excitations, using polarizable embedding alongside time-dependent density functional theory (TDDFT),^{42,43,45,46,50,54,57,59,63} where dynamic mutual polarization poses an additional challenge.⁴⁶ 

We recently presented⁵¹ a novel QM/MM approach (tinktep), which combines the DFT methodology of onetep^64,65 and the polarizable force field AMOEBA,^15,17,18,66 as implemented in tinker.¹⁶ In the tinktep approach, 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 total energy minimization scheme is employed for the Hamiltonian of the coupled QM/MM system. A distinguishing feature of our approach is the use of linear-scaling DFT^65,67 to describe the QM subsystem with the aim of, ultimately, undertaking QM/MMpol calculations with QM regions spanning thousands of atoms.

The main limitation of our first tinktep model, as presented in Ref. 51, was its use of fixed localized orbitals, which represented a trade-off between simplicity and energy accuracy. In this work, we describe an extension of tinktep to the case where the localized orbitals are optimized in situ. The rationale for using optimized orbitals is the near-complete-basis-set accuracy that they offer, even when a minimal basis is used. The resultant accuracy is comparable or superior to even very large bases with fixed orbitals.⁶⁸ 

Incorporating in situ orbital optimization requires computing gradients of all energy terms with respect to the expansion coefficients of the localized orbitals. We present the relevant derivation and describe how the calculation can be implemented to run in linear-scaling time. We subsequently focus on the difficulties that arise as a consequence of using optimized orbitals – QM overpolarization and unphysical charge transfer from QM to MM – and discuss workable solutions to these two problems, using simple QM/MM systems to illustrate our points. We finish by demonstrating the stability, robustness, and accuracy of our model on a number of test cases. We arrive at a robust, mutually polarizable QM/MM model with linear-scaling QM cost, which we show to be more accurate than a non-polarizable QM/GAFF (Generalized amber Force Field) approach, not only in terms of reducing the electrostatic disruption to the QM subsystem but also in terms of improved energetics across the QM/MM interface.

This paper is organized as follows. In Sec. II, we recount the original (fixed-orbital) tinktep approach. In Sec. III, we outline the generalization to in situ optimized orbitals and describe the additional steps that we found to be necessary for obtaining a well-behaved method. The additional steps are best justified using case studies, which, in the interest of clarity of discussion, we relegated to  Appendixes A and  B. Section IV is devoted to validation and demonstration of the utility of the proposed approach on a number of carefully selected systems. Conclusions and closing remarks are found in Sec. V.

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 $μ$⁠, etc.) or Cartesian tensors of rank 2 (e.g., $𝒯 d-dLM$⁠). Matrices with dimensions other than 3 × 3 are typeset with blackboard capitals (e.g., $K$⁠). Indices A, B, and C 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.”

We begin the exposition of the method by briefly recounting the general idea behind tinktep – our first QM/MM approach proposed in Ref. 51. The system is separated into a QM subsystem and an MM subsystem, with the assumption that the separation does not cut through covalent bonds. The total energy of the coupled system is given by

The QM region is described by the density matrix formulation of DFT in the pseudopotential approximation,

with the above terms describing, respectively, the kinetic energy of valence electrons, the Coulombic energy of valence electrons in the (pseudo)potential $vextr$ of the ionic cores, the Hartree energy, the exchange-correlation energy, the mutual Coulombic interaction of N_QM cores having charges $ZI$ and positions $RI$⁠, and empirical dispersion-correction. Open boundary conditions have been assumed. $ρr,r′$ is the density matrix given by

where $φαr$ are non-orthogonal generalized Wannier functions (NGWFs),⁶⁹ which are strictly localized within atom-centered spherical regions. $K=Kαβ$⁠, termed the density kernel, is the matrix representation of the density matrix in the duals of the NGWFs. The electronic (pseudo)density $nr$ is given by

where we assumed a closed-shell system in the interest of brevity. The last term, $EdispQM$⁠, is an empirical dispersion-correction term, which accounts for the well-known deficiency of generalized gradient approximation (GGA) DFT in describing dispersion interactions.⁷⁰ 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⁷¹ formulation, with parameters determined by Hill and Skylaris.⁷⁰ 

The MM subsystem is described by the AMOEBA¹⁸ polarizable force-field, as implemented in the tinker¹⁶ code, with the following general energy expression:

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

$EpermMM$ is a sum of purely Coulombic multipolar interactions between atoms in the MM subsystem, with scaling factors¹⁷ used to attenuate or eliminate interactions between first-, second-, third-, and fourth-nearest neighbors (as determined by bond connectivity). The full expression is given in Ref. 18, Eqs. (1) and (10). The $EpermMM$ term is fully local to the MM subsystem that is to say it is insensitive to the presence of the coupling with a QM subsystem.

$EpolMM$ is the polarization energy of the MM subsystem given by [cf. Ref. 51, Eq. (A1)]

where E_L is the direct electric field at site L, $μL$ is the dipole induced at site L in response to the total electric field, and T_LM is a 3 × 3 coupling tensor between sites L and M,

Here, $𝒯LM d-d$ is the Thole-damped, Cartesian dipole-dipole interaction tensor between induced dipoles at sites L and M [cf. Ref. 51, Eq. (28)]. Thole damping^72–74 is a modification to Coulombic electrostatics that helps prevent mutual positive feedback loops involving induced point dipoles, known as a polarization catastrophe.

In a purely MM calculation, the direct electric field E_L is simply the electric field due to the permanent multipoles of MM sites. In a mutually polarizable QM/MM calculation, the direct electric field includes an additional contribution arising from the multipoles representing the QM subsystem. We give the full expression in Ref. 51, Eqs. (26)–(31), where we also explain in detail how a classical multipolar representation of a distributed QM density is obtained (Sec. II D therein). We stress that the non-additive nature of polarization means that the polarization of MM cannot be separated into additive terms due to QM and MM, and so the entire polarization of the MM subsystem is included in $EpolMM$⁠, which explains the absence of an $EpolQM/MM$ term.

$EvalMM$ denotes all short-range valence interactions local to the MM subsystem. The detailed expressions for these terms can be found in Ref. 18, Eqs. (2)–(6), and we shall refrain from recounting them here.

$EvdWMM$ accounts for van der Waals (dispersion-repulsion) interactions local to the MM subsystem. AMOEBA uses the Halgren formulation⁷⁵ of the buffered 14-7 potential,

where ρ_ij = R_ij/R⁰, δ = 0.07, and γ = 0.12. The parameters of the potential are R⁰ and ε. Mixing rules for obtaining pairwise values of the parameters and a description of nuances surrounding hydrogen atoms (“reduction factors”) can be found in Ref. 18, Eqs. (7) and (8).

The final term in (1), $EQM/MM$⁠, accounts for all interactions between the two subsystems, except for mutual polarization. As pointed out earlier, QM contributions to MM polarization have already been included in $EpolMM$ because they are not separable from intra-MM polarization. The effect of MM polarizing QM is automatically included in $EQMnr$ by the deformation of the electronic density $nr$ in response to the electric field of the MM subsystem.

The coupling is described by

where the first term accounts for the electrostatic coupling between QM and permanent MM multipoles and the second term accounts for dispersion-repulsion interactions between QM and MM. In our model, the electrostatic coupling involves the full QM charge density interacting with the Coulombic (not damped) potential of the permanent MM charges, dipoles and quadrupoles,

where

and the expression for $vpMMr$ can be found elsewhere [Ref. 51, Eqs. (37) and (38)].

For the QM/MM van der Waals interaction, $EvdWQM/MM$⁠, our original model uses the same classical, pairwise model that is used for MM/MM [cf. (8)], except the repulsive wall is softened slightly by using δ = 0.21 (cf. Ref. 51, Sec. III.B.3). This has the advantage of being straightforward but has two disadvantages. First, it requires choosing vdW parameter values for atoms in the QM subsystem. Second, and more importantly, this classical form is insensitive to the electronic density of the QM subsystem and as such its contributions to electronic density gradients vanish. This means it fails to provide the Pauli repulsion that would otherwise prevent electrons from unphysically collapsing onto MM atoms. We address this issue in the revised model presented in this paper.

The main limitation of the original model described briefly above, and in detail in Ref. 51, was that the localized orbitals $φα$ were kept fixed and only the density kernel K^αβ was optimized. Allowing $φα$ to be optimized in situ constitutes the main improvement in our revised model.

Optimizing the NGWFs necessitates deriving and implementing functional derivatives of energy terms with respect to the NGWFs. Out of all the energy terms that make up the total energy (1), the following terms are specific to our QM/MM model and have no counterparts in standard onetep: $EvalMM$⁠, $EvdWMM$⁠, $EpermMM$⁠, $EvdWQM/MM$⁠, $EpermQM/MM$⁠, and $EpolMM$⁠. The first three of these are local to the MM subsystem and do not depend on the electronic degrees of freedom, their derivatives with respect to the NGWFs thus vanish. $EvdWQM/MM$ is, so far, described by an electron-independent classical pairwise sum, so its derivative similarly vanishes. We postpone the generalization of this term to an electron-dependent form until later in the text. The remaining terms are $EpermQM/MM$ and $EpolMM$⁠, which we will consider now in sequence.

$EpermQM/MM$ is given by [cf. (10) and (11)]

The potential of permanent MM multipoles $vpMM$ does not depend on the electronic degrees of freedom, so the second term above does not contribute to the derivative. The first term represents an interaction of an electronic density with a fixed potential and so has the same form as the second term in (2). Thus it can be accounted for using usual onetep algorithms by simply adding $vpMMr$ to $vextr$⁠.

The term due to $EpolMM$ is more complicated since it involves differentiating the transformation from the QM density to the set of multipoles representing the QM subsystem that takes part in QM/MM polarization interactions. We refer the reader to Ref. 51, Sec. I.D.2 for a detailed exposition of this transformation, recounting here only the basics needed in the derivation.

All pairs (products) of overlapping NGWFs φ_α and φ_β are expanded in terms of truncated spherical waves $fsr$ centered on both NGWFs, with the coefficients of the expansion given by

where s and t index the spherical waves originating on both atoms (of which there are N_f in total), $VABts$ is an element of the inverse electrostatic overlap matrix between spherical waves originating on atoms A and B, and the notation α ∈ A used in the summations is taken to mean “all NGWFs α belonging to atom A.”

An alternative way to index the spherical waves, and in turn the coefficients $CABs$ and $cαβs$⁠, is via their angular, magnetic, and radial numbers: l, m, and q and a selector for the site on which the spherical wave originates (1 for the first atom from the subscript or 2 for the second atom). This indexing scheme is useful when using the expansion to calculate the spherical multipole moments $Mlm(A)$ that constitute the classical representation of the density,

where $CABlmq,1$ correspond to $CABs$ originating only on atom A and J_lq is a radial, analytical integral given in Ref. 51, Eq. (21). The notation $BSAB≠0$ is taken to mean “atoms B whose NGWFs overlap with those of atom A.”

The interaction energy between all multipoles representing the QM density and the induced MM dipoles is given by

where $P$ is the polarization matrix, with matrix elements α ∈ A, β ∈ B given by

where $wAlm$ captures the electrostatic effect of the entire system of MM induced dipoles on the QM site at R_A: $wA00$ is the Thole-damped electrostatic potential of MM induced dipoles interacting with the charge at R_A, ${wA1,−1,wA1,0,wA1,1}$ is their Thole-damped electric field interacting with the dipole at R_A, and ${wA1,−2,…,wA1,2}$ is their Thole-damped electric field derivative (in spherical representation) interacting with the quadrupole at R_A.

Taking the functional derivative of (15), we obtain [cf. Ref. 76, Eq. (7.25)]

where the Einstein convention has been used for repeated Greek indices.

The non-vanishing derivative $δKβαδφγ*r$ appearing in the second term results from the use of a so-called purifying transformation in onetep. It has already been derived [cf. Ref. 76, Eq. (7.29) or Ref. 77, Eq. (4.4.6)] and implemented in onetep and as such can be omitted from further discussion. Instead, we focus on the more interesting first term, involving the quantity $δPαβδφγ*r$⁠. From (16), we obtain

where we have used the fact that J_lq is independent of the NGWFs. The quantities $wAlm$ and $wBlm$ involve induced MM dipoles and, as these induce in response to a combined QM + MM electric field, they depend on the NGWFs. However, a zero residual condition at induced dipole self-consistency

obviates the need to calculate $δμLδφγ*r$ [compare Ref. 51, Eq. (33)], allowing us to only consider the dependence of $cαβlmq,1$ and $cαβlmq,2$ on the NGWFs, and so

The remaining functional derivative $δcαβlmq,1δφγ*r$ can be calculated as follows:

where $vtr$ is the potential of a spherical wave (for which an analytical expression is available). $δcαβlmq,2δφγ*r$ is calculated analogously.

By combining (17), (20), and (21) and expressing the sums in the first term explicitly, we obtain

where atom C is the host to the NGWF with respect we differentiate, i.e., γ ∈ C. The presence of both $wClm$ and $wAlm$ in (22) indicates that the gradient with respect to a particular NGWF γ ∈ C depends not only on the electrostatic effect of MM’s induced dipoles at R_C but also at all centres of overlapping NGWFs R_A. This is a consequence of the two-center spherical wave expansion scheme used in onetep.

To maintain linear scaling, an implementation must be able to evaluate (22) in $O1$ time since this calculation must be repeated for all NGWFs γ, and the number of NGWFs is proportional to N_QM. Our implementation in onetep does this by re-ordering (22) as

where

only needs to be evaluated in the intersection of the localization spheres of A and C [cf. (22)]. The cost of evaluating U_AC for a single pair of atoms A-C is system-size independent (⁠$O(1)$⁠) and only depends on the quality of the SW basis set. For any particular atom C, the number of atoms A whose localization spheres overlap with it plateaus at a constant that depends on the density of the system, even if the system size N_QM is increased to arbitrarily large values. This is made explicit by the A-C overlap condition in the first summation of (22). This means that each evaluation of (22) has $O(1)$ cost and, with $O(N)$ such operations, the approach is linear-scaling.

Polarizable force-field models that rely on the induced point dipole approximation have to contend with what is known as a polarization catastrophe. This well-known⁷⁸ artifact consists in an unbounded mutual polarization of two nearby sites through positive feedback and reflects the breakdown of the point-dipole model at a short range. The polarization catastrophe is typically mitigated by replacing Coulombic interactions involving induced dipoles with interactions that are suitably attenuated at a short range using schemes such as Thole damping.^72–74 This is the case in AMOEBA¹⁵ and in our model.⁵¹ The intensity of the damping depends on the polarizabilities of the two atoms – i.e., interactions involving atoms that polarize more readily are more aggressively damped. Beyond several Å, the difference between the Thole-damped and Coulombic quantities (potential, electric field, electric field derivative) becomes negligible and the correct long-range behavior is recovered.

The rationale for using optimized NGWFs in our model, and in onetep in general, is the near-complete-basis-set accuracy that they offer, comparable or superior to even very large bases with fixed orbitals.⁶⁸ However, this additional flexibility results in the basis becoming more diffuse, which is problematic in the context of distributed multipole analysis (DMA^79,80) that we employ to obtain the multipole representation QM^* [cf. (13)–(16)]. More diffuse bases are known^81,82 to engender instabilities in the DMA procedure and lead to increased charge penetration errors (CPE) due to discrepancies between the potential of the original density and that of the multipolar expansion although improved approaches have recently been proposed.^82,83 Indeed, our initial tests revealed that once the NGWFs are no longer fixed and are allowed to change shape during the SCF process, our QM/MM model becomes prone to a QM/MM analog of polarization catastrophe, whereby the QM subsystem becomes excessively polarized by a nearby MM site and vice versa. The problem is particularly severe for MM sites carrying a charge (ions) as they provide a larger initial polarization of QM. We devote a section in the  Appendix A to an elucidation of this mode of failure using a H₂O:Cl⁻ system as an example. In the same section, we show that a simple increase in the damping of QM/MM polarization interactions is sufficient to prevent the QM/MM polarization catastrophe.

Like most QM/MM models, our initial model used a classical, atom-pairwise description of QM/MM dispersion-repulsion (vdW) interactions [cf. (8)]. This strictly classical description has the disadvantage of being insensitive to the electronic degrees of freedom in the QM subsystem that is to say the QM/MM vdW energy only depends on the positions and species of the atoms. The most striking manifestation of this deficiency is that electrons in the QM subsystem do not experience any Pauli repulsion from MM sites. This can be especially problematic when the MM site is a cation, whose electrostatic potential attracts the QM electrons. With no Pauli repulsion to balance this attraction, unphysical charge transfer from QM to MM takes place.

A number of approaches have been proposed to circumvent the problem (see, e.g., Refs. 84–87), but not in the context of linear-scaling QM methods, where it becomes particularly problematic. This is because the spilled electrons accumulate near the peripheries of the localization regions, disrupting SCF convergence, which assumes localized orbitals to be well-decayed at the truncation point.

We refer the reader to the  Appendix B for a case study of this undesired effect on a H₂O:K⁺ system, where we also demonstrate the feasibility and accuracy of an improved model which eliminates this issue. The improvement consists in replacing the repulsive term of Halgren’s vdW potential with a density overlap model⁸⁶ that is sensitive to the QM electronic degrees of freedom, accounting for QM/MM Pauli repulsion. We retain Halgren’s classical description for QM/MM dispersion interactions.

In an overlap model, the Pauli repulsion energy is assumed to be proportional to the overlap between densities, i.e.,

where $nr$ is the QM electronic density, $nLMMr$ is a model density centered on MM atom L, and κ_L is a proportionality constant with a suitable unit.

A reasonable model density is that of a 1s Slater-type function

Instead of working with model densities, we can think of MM atoms as equipped with a model repulsive electrostatic potential, leading to an equivalent energy expression

together with an equivalent MM repulsive potential

characterized by two parameters – a magnitude A and an inverse-width ζ, both of which depend on the chemical species of MM atom L.

The form of (27) is that of a static external potential acting on the electronic density, which means $vMMrepr=∑LNMMvLMMrepr−RL$ can simply be added to $vextr$ in (2), and no new energy gradient expressions need to be derived for this term. Unless stated otherwise, all results presented in this paper have been obtained with the model that includes the MM repulsive potential (and excludes the repulsive contribution from the Halgren QM/MM vdW expression).

Regarding computational efficiency, we point out that the integral in (27) only needs to be computed over the union of localization spheres of the QM subsystem (since $nr$ vanishes elsewhere). Furthermore, and more importantly, since $vLMMrepr$ decays exponentially, only those regions of the QM subsystem that are within a short cutoff radius (say, 5 Å) from any MM atom need to be considered. Generating $vLMMrepr$ in a sphere around R_L, with the sphere radius system-size independent, is an $O(1)$ operation for a single MM atom L. The number of MM atoms within a cutoff radius from the QM subsystem will be proportional to the surface area of the QM subsystem and so to $NQM2/3$⁠. The total cost of evaluating (27) thus scales $O(NQM2/3)$⁠.

Naturally, physically reasonable values for A_L and ζ_L need to be determined for all species of interest appearing in the MM subsystem. As a proof of concept, in Sec. IV A, we show how suitable values can be found for Cl⁻, K⁺, and H₂O.

In this section, we demonstrate the accuracy and viability of the proposed approach on a number of systems. In all QM calculations, we used the Perdew-Burke-Ernzerhof (PBE)⁸⁸ exchange-correlation functional, with an empirical dispersion correction in the Elstner⁷¹ formulation, with parameters determined by Hill and Skylaris.⁷⁰ The NGWF localization radius was set to 3.7 Å.

We begin by examining the interaction energy curves of three simple systems: H₂O:K⁺, H₂O:Cl⁻, and a water dimer. The latter two systems were studied in our earlier studies,^5,52 using a different QM/MM model and using energy decomposition analysis (EDA) to compare AMOEBA against a high-quality DFT functional ωB97X-V.⁸⁹ For each of the systems, we compare the predictions of the QM/MM model that is the focus of this paper and those of AMOEBA, against reference results obtained from fully QM calculations (i.e., PBE-D as described above). All QM calculations used a kinetic energy cutoff of 1290 eV.

By performing a parameter scan in the space of $A,ζ$⁠, we established MM repulsive potential parameters for K⁺ that, for this system, are optimal in the sense of minimizing the mean squared difference between the interaction energy curves from QM/MM and fully QM calculations. The values we obtained are $AK+=230Ha/e$ and $ζK+=1.379a0−1$⁠. The interaction energy curves are compared in Fig. 1. AMOEBA (green curve) is seen to model this interaction faithfully, with the position of the minimum accurate to 0.007 Å and only very slight underbinding (less than 1 kcal/mol). This degree of agreement is expected since the charge density of K⁺ is tightly localized and thus well-approximated by a point multipole model, with very little charge penetration error. The QM H₂O/MM K⁺ description is in even better agreement with full QM – the position of the minimum is accurate to 0.002 Å and the energy is no further than 0.3 kcal/mol from the fully QM result for all interatomic separations.

We now turn our attention to the H₂O:Cl⁻ system, which we expect to be more difficult for a polarizable point dipole model due to the larger electronic delocalization of Cl⁻, which increases the charge penetration error. By following the same protocol as for the H₂O:K⁺ system, we established optimal parameters for the repulsive MM potential for Cl⁻: $ACl−=250Ha/e$ and $ζCl−=1.140a0−1$⁠, which, compared to K⁺, represent a marginally stronger and somewhat less localized potential, consistent with expectations.

The interaction energy curves are compared in Fig. 2. Compared to our QM reference, AMOEBA is seen to underbind at all interatomic separations, particularly at short distances, where the magnitude of the error increases from ≈1 kcal/mol to over 5 kcal/mol. A large fraction of this error can be attributed to the neglect of charge transfer. In our reference, fully QM calculations as much as 0.22e is transferred from the Cl⁻ ion to the water molecule at the shortest studied separation (2.8 Å), corresponding to a stabilizing effect of ≈−5 kcal/mol. As the separation is increased, this charge transfer becomes less pronounced – at 4 Å only 0.05e is transferred and the corresponding change in energy is only about −0.2 kcal/mol.

However, we must acknowledge the fact that the QM reference curve against which AMOEBA is benchmarked is too a result of a physical model. DFT calculations involve a number of approximations, chief among which are the use of an approximate exchange-correlation functional, the pseudopotential approximation, the use of a finite basis, and – in linear-scaling DFT – the use of finite radii for the localized orbitals. Different choices for these parameters will lead to slightly, but noticeably, different interaction energy curves, particularly since the water molecule is well-known to be difficult to describe with GGA DFT (see, e.g., Refs. 90–92). For instance, when instead of PBE, ωB97X-V⁸⁹ is used for the same system (as reported by some of us in Ref. 5, cf. Fig. 3 therein), the reference curve shifts upwards by ≈1 kcal/mol in the long range (practically matching AMOEBA), and by as much as ≈3.2 kcal/mol at 2.8 Å, reducing AMOEBA’s perceived underbinding at the shortest separation studied here to 2 kcal/mol. We thus caution against treating all differences between MM and QM reported here strictly as deficiencies of the MM model.

Naturally, the charge penetration error is also expected to be more significant for Cl⁻ than for K⁺. AMOEBA does not explicitly model charge transfer or account for charge penetration and must resort to approximating these effects through polarization and vdW interactions. More severe underbinding at short distances leads to a shift in the position of the minimum, which, compared to our QM reference, AMOEBA overestimates by 0.08 Å (or 0.05 Å against the QM reference of Ref. 5). Our QM/MM model achieves better agreement with fully QM results, underbinding by less than 2 kcal/mol, with the magnitude of the error being almost independent of the distance between Cl⁻ and the water molecule. Thus, the predicted interaction energy curve is very similar in shape to the reference one, only shifted by a constant, and the position of the minimum is predicted very accurately (to 0.002 Å), showing the QM/MM interface to be well-balanced in this scenario.

We now turn our attention to the H₂O dimer. In earlier work⁵² on the same system, we showed that charge penetration is significant at the equilibrium distance and below it, making this system challenging for AMOEBA, which has to compensate for CPE by artificially softening the repulsive vdW wall, relying on cancellation of errors to model the hydrogen bond. Thus (cf. Fig. 3), the agreement between AMOEBA and a fully QM calculation worsens at short separations, where AMOEBA underbinds by as much as 4 kcal/mol (2.9 kcal/mol against the QM reference of Ref. 5), but is still remarkably good at the equilibrium distance and beyond, where AMOEBA underbinds by only ≈0.5 kcal/mol. The rms error across the entire curve is 1.1 kcal/mol. The position of the minimum is also predicted accurately (to 0.004 Å).

Determining suitable parameters for our model’s repulsive MM potential for O and H atoms is more challenging than in the previous two cases. First, there are four parameters to be simultaneously optimized (A_H, ζ_H, A_O, ζ_O), making the parameter scan more involved. Second, if our model is to be well-transferable, it must accurately describe both the situation where the hydrogen bond donor is described by QM (and the acceptor by MM), and the situation where QM is used to describe the hydrogen bond acceptor (and MM – the donor).

After a thorough parameter scan, we determined the following suitable values for the parameters of the MM repulsive potential: A_H = 35 Ha/e, $ζH=2.40a0−1$⁠, A_O = 550 Ha/e, $ζO=1.58a0−1$⁠. With these values, the interaction energy curves predicted by our QM/MM model (Fig. 3, blue curves) are in very good agreement with fully QM results. When QM is used to model the hydrogen bond donor (solid blue curve), the rms error in energy across the entire curve is 0.7 kcal/mol and the position of the minimum is accurate to 0.002 Å. When QM is used to model the hydrogen bond acceptor (dashed blue curve), the rms error in energy is only 0.4 kcal/mol, but the position of the minimum is predicted less accurately and is underestimated by 0.026 Å. Crucially, in both cases, the interaction energy curve stays within 1 kcal/mol from the reference curve obtained with fully QM calculations, even at separations below the equilibrium distance. This indicates that the QM/MM interface in our model is well-balanced even in the presence of hydrogen bonds.

The above examination of the performance of our QM/MM model for three representative systems (MM cation, MM anion, and MM neutral molecule with a hydrogen bond spanning the QM/MM interface) can be considered a proof of concept. We showed that our QM/MM model is stable for all studied intermolecular separations, even well below the equilibrium distance, and that it gives reasonable predictions for interaction energy profiles, which we find remarkable given that in the studied systems, the crucial interactions crossed the QM/MM boundary.

In all examples so far, we used the H₂O molecule for the QM subsystem and so the question of whether our QM/MM model is transferable, particularly concerning the parameterization of the MM repulsive potential, remains open. In the text that follows, we will examine the model’s performance for a number of different molecules, both neutral and charged, demonstrating that it is indeed transferable as its predictions remain accurate.

The need to determine suitable parameters for the repulsive MM potential of all MM species of interest can be seen as a weakness of our model. Intuitively, one would hope that the parameters A_L and ζ_L could be derived from corresponding classical vdW parameters ε_L and $RL0$ – e.g., we expected A_L ∼ ε_L and $ζL∼1/RL0$⁠. However, we found this not to be the case. For instance, for the parameters, we determined $AHAO≈15.7$⁠, whereas $εHεO≈8.1$⁠, that is, our model uses a substantially weaker potential on H. Similarly, we have $ζHζO≈1.52$⁠, whereas $1/RH01/RO0≈1.28$⁠, meaning the potential on H used in our model is also slightly tighter.

One reason is that while $vLMMrepr$ [cf. (28)] is linear in A_L just like $EvdWMMRij$ [cf. (8)] is linear in ε_L, the energy expression (27) for $ErepQM/MM$ is not linear in A_L. This is because $nr$ implicitly depends on A_L, that is to say, the electronic density responds to the MM repulsive potential by deforming accordingly. Thus, not only is $ErepQM/MM$ not linear in A_L, but also other energy terms are indirectly influenced by the repulsive MM potential through the change in $nr$⁠. Another reason is that the vdW parameters adopted in AMOEBA have been fitted to partially compensate for the deficiencies in the classical treatment of electrostatics, some of which are no longer present in our QM/MM formulation. Finally, AMOEBA employs additional “tweaks” in its vdW formulation, for instance, the repulsive sites of H atoms are slightly offset from the actual atomic sites (“reduction factor”).

While we plan to investigate routes for automatically obtaining A_L and ζ_L in future work, in this paper, we will focus on MM water, for which we have obtained good parameters already. In this way, we can apply our model to a large class of systems that is of practical interest – QM solutes embedded in MM water. We defer applications with different MM species to a later time.

We now set out to demonstrate the transferability of our model, turning our attention to a number of QM solutes embedded in spherical shells of MM water. We will use the same systems and the same methodology as in our earlier work⁵¹ – the QM subsystem will only encompass the solute, and we will study the behavior of the QM/MM system as the size of the MM H₂O shell is increased (cf. Fig. 3 in Ref. 51).

Three of the solutes were chosen from the SAMPL4 blind challenge:⁹³ (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₃), (e) the ammonium ion (⁠$NH4+$⁠), and (f) the cyanide ion (CN⁻) – were chosen with the aim of verifying if our model correctly describes small and charged solutes.

We compared 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₂O molecules were described with fixed partial charges. In this setup, only a fixed, external potential is included in the QM Hamiltonian; we emphasize the neglect of vdW interactions between the QM and the embedding;

• QM/MM calculations with either a fixed point-charge embedding (GAFF v1.5⁹⁴) 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. For the fixed point-charge (GAFF) embedding, vdW interactions between the solvent and solute were included at the MM level of theory (Lennard-Jones potential). Thus, the MM repulsive potential introduced in Sec. III C was not used in this case. Similarly, polarization damping (Sec. III B) was not relevant here as the force field was not polarizable. For the polarizable embedding (AMOEBA), we used the final, refined QM/MM model, as described in Secs. III B–III C;

• Fully MM calculations, where the entire system was treated classically (with GAFF or AMOEBA).

In fixed point charge QM/MM calculations and in QM calculations with fixed point-charge embedding, we used partial charges of 0.417e for H atoms and −0.834e for O atoms, which are identical to the TIP3P⁹⁵ model used in GAFF. All QM calculations used a kinetic energy cutoff of 1000 eV. The configurations were prepared by solvating the solutes in approximately 660 explicit H₂O molecules under periodic boundary conditions. Classical polarizable MD trajectories in the NpT ensemble (p = 1 atm, T = 298 K) were then obtained, and the final configurations after 50 ps were used. Final densities were 1.00 ± 0.01 g cm⁻³. For details on how the configurations were processed for calculations under open boundary conditions, details of cutoffs assumed, and the models and parameters used for the solutes and water, see Ref. 51, Sec. 3.B.1.

We begin by comparing the interaction energies between the QM solute and the MM water shell as a function of the size of the shell (number of H₂O molecules). To better elucidate the long-range behavior, in Fig. 4, we only plot the error in the energy with respect to the fully QM calculation that we use as reference. Even though the systems studied here are the same as in our earlier work,⁵¹ we point out that the energy error curves are not directly comparable with those of Ref. 51, since in current work, we used a more refined QM approach with in situ optimized NGWFs both for QM/MM calculations and for the fully QM reference.

Our first observation is that for all six systems, purely MM calculations with AMOEBA (green circles) clearly outperform GAFF (gray diamonds). With the exception of NH₃ and CN⁻, typical errors in the GAFF description are at least three times larger than their AMOEBA counterparts. Even in cases where GAFF fares relatively well (NH₃), or where AMOEBA’s error is rather large (CN⁻), long-range behavior is clearly much better described by AMOEBA. This can be appreciated from the much flatter profiles of the AMOEBA curves, which indicate an almost constant energy shift from the reference QM calculation. The energy changes from adding subsequent water molecules are more erratic for GAFF, and the convergence with the number of H₂O molecules is much worse. As expected, this is particularly pronounced for charged solutes – in the case of CN⁻, for instance, the binding energy is not well converged even at 400 MM H₂O molecules. We attribute this to polarization partly compensating for the boundary effects that result from truncating the water shells.

The behavior of point-charge electrostatic embedding (QM + EE, red diamonds) is best understood by comparing it against QM/MM with GAFF (orange crosses) since these two approaches only differ by the neglect of QM/MM vdW interactions in the former. This neglect leads to a very rapid accumulation of error at a short range, particularly for larger solutes (a)–(c), where this interaction is more significant. At longer range, QM/MM vdW interactions are well-decayed, which is reflected in the almost identical profiles of QM/MM-GAFF and QM + EE curves starting at approximately 150 H₂O molecules for large solutes (a)–(c) and as early as approximately 50 H₂O molecules for small solutes (d)–(f). As expected, the neglect of QM/MM vdW interactions makes the QM + EE approach inadequate for calculating interaction energies between the QM solute and embedding although occasionally (e.g., for NH₃ and CN⁻) the error fortuitously cancels out some of the errors in the electrostatics.

It is also worthwhile to compare the results of QM/MM-GAFF (orange crosses) against purely MM GAFF calculations (gray diamonds) because it reveals the effect of treating the solute at the QM level of theory, all other components of the two models being identical. For all six solutes, QM/MM-GAFF is more accurate, and, as the long-range profiles of the two curves are almost identical, it is clear that this gain in accuracy is due to a much improved description of short-range interactions, i.e., the ability of the QM subsystem to realistically polarize in response to the MM environment.

Of greatest interest to this paper is, of course, the comparison between QM/MM-AMOEBA (blue squares) and QM/MM-GAFF (orange crosses). In terms of absolute errors in energy, our model outperforms QM/MM-GAFF in all cases except for NH₃. Moreover, the long-range behavior of QM/MM-AMOEBA is much better (flatter curves), particularly for charged solutes, where all fixed-point charge approaches (MM-GAFF, QM/MM-GAFF, QM + EE) clearly suffer from neglecting polarization. Out of all five models, the QM/MM-AMOEBA model has the lowest maximum error in the long range (4 kcal/mol for CN⁻, compared with 9 kcal/mol of AMOEBA, and maximum errors in excess of 10 kcal/mol for the other approaches). We summarize these results in Table I, from which it also becomes clear that, when averaged over all six systems, our approach has the lowest signed and unsigned errors of all the considered approaches.

Satisfied that the energetics of our QM/MM interface is accurate, we now focus on how the QM solute is affected by the presence of the QM/MM interface. Naturally, we would like the electronic structure of the QM solute in the presence of MM embedding to resemble the electronic structure of the full QM system as much as possible, i.e., for the MM embedding to faithfully mimic QM. Since we cannot compare electronic energy levels between QM/MM and full QM and comparing electronic densities would require density partitioning, we will use the total dipole moment of the solute as a proxy.

In Fig. 5, we plot the magnitudes of the solute dipole moment for the same six solutes as a function of the size of the water shell. In QM/MM calculations, the solute (QM) dipoles are immediately available. In fully QM calculations, the solute dipoles were obtained from DMA analysis. In fully MM calculations, the solute dipoles are either obtained by a suitable vector summation of permanent dipoles with induced dipoles (AMOEBA) or, in the absence of polarization, are simply constant (GAFF). For neutral systems, the dipole moment is invariant to the choice of the reference point. For charged systems, and in fully QM calculations where charge transfer between the solute and solvent can make the total solute charge non-zero, we chose the centroid of the molecule as the reference point.

Our first observation is that, particularly for larger solutes, the solute dipole moment is rather sensitive to the environment and can change abruptly depending on where subsequent H₂O molecules are added. The qualitative behavior of this sensitivity is captured to a similar degree by all models, except of course MM-GAFF, which does not permit solute polarization. The accuracy of the constant dipole moment of the GAFF solute is hit-and-miss – e.g., GAFF’s prediction is excellent for (−)-menthol, severely underestimated for diphenylhydramine and $NH4+$⁠, and severely overestimated for CN⁻.

The predictions of QM + EE (red diamonds) and QM/MM-GAFF (orange crosses) are expected to be identical since the two approaches only differ by the absence/presence of a classical QM/MM vdW term that does not affect the electronic degrees of freedom. In practice, we observe very small differences (<0.1 D) that are the consequence of different smearing of the singularities of the Coulombic permanent fixed-point charges of the embedding on the Cartesian grid on which electronic density is evaluated in onetep.

For all six solutes, the predictions of our QM/MM model are more accurate than those of AMOEBA, indicating the expected superiority of a QM description of the solute (this is most striking for $NH4+$⁠, which is underpolarized with AMOEBA). For three of the six solutes ((−)-menthol, NH₃, CN⁻), the predictions of our model are closest to the fully QM results in absolute terms. For the remaining three molecules, our model is slightly less accurate than QM/MM-GAFF, but not much so. Furthermore, this only happens when AMOEBA itself fares worse (diphenylhydramine, 2-Cl-4-OH-3,5-dimethoxy-BALD, and $NH4+$⁠), possibly implicating the polarizable water model, rather than the QM/MM interface, as the culprit. When the errors are averaged over all the systems, our QM/MM model yields the lowest root mean square error (RMSE). Details are summarized in Table II.

We now turn our attention to the interaction between three solutes: H₂O, Cl⁻, Na⁺, and their first solvation shells. The three solutes are meant to be representative of neutral, anionic, and cationic species, respectively. To benchmark our QM/MM approach, we model only the solute at the QM level of theory, while the solvent (water) is being described by AMOEBA. We calculate solute-solvent interaction energies, comparing the performance of our approach against a non-polarizable model (QM/MM-GAFF) and purely classical models (where the entire system is described with GAFF or AMOEBA). Our aim is to verify whether our QM/MM interface correctly reproduces the purely QM results obtained at the same level of theory. Thus, our reference is pseudopotential DFT + D with PBE, the limitations of which we acknowledged earlier. We caution the reader against interpreting the discrepancies between MM results and our reference as “failures” of MM in absolute terms.

To obtain meaningful statistics, we performed calculations for 100 configurations (for each solute) obtained from MD runs, where each solute was solvated in 215 H₂O molecules. In each MD snapshot, all but N_1st H₂O molecules closest to the solute were then stripped, leaving only the first solvation shell. The values of N_1st were 4, 8, and 6 for H₂O, Cl⁻, and Na⁺, respectively. The configurations studied are the same as used in our earlier work.⁵² For a more detailed description of how the configurations were obtained and the rationale for choosing N_1st, we refer the reader to Ref. 52, Sec. III.C.

The calculated solute-solvent interaction energies are plotted in Fig. 6, while Table III reports crucial statistics. For the H₂O solute, our QM/AMOEBA model performs well, with a general trend of underbinding by about 1 kcal/mol, and comparing favorably against QM/GAFF in all five metrics (rms error, mean signed error, maximum error, correlation coefficient, and slope of linear fit to the model’s interaction energy vs. reference). We find this rather remarkable since QM/MM models are typically not very good at reproducing interaction energies spanning the interface. The largest error for our model is 3.7 kcal/mol, which is better than pure AMOEBA (4.6 kcal/mol), and much better than pure GAFF (6.2 kcal/mol) or QM/GAFF (7.3 kcal/mol). The correlation between the model and reference is also very good (r = 0.95) although pure AMOEBA does marginally better (r = 0.96).

The Cl⁻ solute exposes the weaknesses of purely MM treatments. GAFF is particularly inaccurate here, with rms and mean signed errors in excess of 10 kcal/mol and a maximum error of over 40 kcal/mol, which is not surprising, given the likely importance of polarization effects in this system, for both the solute and solvent. These results were obtained using the vdW parameters proposed by Fox and Kollman:⁹⁶ R^* = 1.948 Å, ε = 0.265 kcal/mol, which we used for consistency with previous work. We note in passing that GAFF performs slightly better here when vdW parameters tailored specifically for Cl⁻ in TIP3P water are used (Joung and Cheatham:⁹⁷ R^* = 2.513 Å, ε = 0.0356 kcal/mol). The improved values are shown in parentheses in Table III.

AMOEBA’s predictions are better (an rms error of 6.6 kcal/mol), but it does not avoid occasional embarrassments (a max error of 20.7 kcal/mol). Both QM/MM models perform significantly better, which highlights the importance of treating the Cl⁻ ion at the QM level of theory, in order to be consistent with the latter. Our QM/AMOEBA model correlates better with pure QM (r = 0.98 against r = 0.92 for QM/GAFF, a linear slope of 0.97 against 1.23 for QM/GAFF), but it is seen to underbind slightly across the board (rms error of 4.1 kcal/mol, compared to 3.9 kcal/mol for QM/GAFF). Its maximum error is 7.1 kcal/mol, which is rather large, but still better than QM/GAFF (9.4 kcal/mol) and much better than the double-digit errors of MM models.

For Na⁺ and its first solvation shell, AMOEBA performs very well, while all the remaining models are rather inaccurate. Since Na⁺ is a compact, barely polarizable ion, it is well-described by MM methods. This explains why GAFF results are almost identical to QM/GAFF results. What is significantly more important in this system is the description of the water solvent. GAFF’s water model cannot capture the polarization of the solvent, which is highly relevant here, due to the charge on the Na⁺ solute. Thus GAFF and QM/GAFF both yield a poor description of the whole system, with rms errors above 10 kcal/mol and maximum errors of almost 20 kcal/mol. This is despite using vdW parameters for Na⁺ that were specifically tailored for a sodium ion in TIP3P water (Ref. 97, R^* = 1.369 Å, ε = 0.087439 kcal/mol). AMOEBA, by contrast, performs very well, with an rms error of only 1 kcal/mol and good correlation with purely QM results (r = 0.96), highlighting the importance of a polarizable description of the water solvent. Since our QM/AMOEBA model shares its description of the solvent with AMOEBA, one would expect it to yield a similarly good description. However, this is not the case. While the correlation with purely QM results is good (r = 0.95, slope of 0.99), there is significant overbinding for all snapshots, leading to large errors in energy, dominated by a mean signed error of −13.3 kcal/mol. This almost constant shift points to a deficiency of our QM/MM interface in handling cationic solutes, presumably due to the repulsive MM potential having been parameterized only using H₂O–H₂O interactions. We attribute the observed overbinding to an insufficient repulsion between the compact Na⁺ and nearby MM oxygen atoms.

We will now briefly investigate the effect of using in situ optimized orbitals on the quality of QM/AMOEBA. We calculated the interaction energies of the three systems from Fig. 6 using fixed pseudoatomic orbitals [single-zeta (SZ), double-zeta and polarization (DZP), triple-zeta and polarization (TZP)] and compared them with results obtained using an in situ optimized minimal basis. We report the results in Fig. 7 and Table IV. As expected, for a compact Na⁺ cation, a minimal (SZ) basis is sufficient, and increasing the flexibility of the basis makes very little difference, with the slope and correlation coefficients practically unchanged, and an essentially rigid shift of the interaction energies by 0.5–0.7 kcal/mol compared to an in situ optimized basis. Since our model systematically overbinds this system, the resultant shift actually makes the fixed-basis results marginally better, owing to cancellation of a small fraction of the error. For a diffuse Cl⁻ anion, the effect of the basis size is more pronounced, with a clear systematic improvement of about 1.5 kcal/mol in the mean signed error (MSE) for each time the basis set quality is increased although the optimized basis “overshoots” by a small amount, underbinding by 3.9 kcal/mol. Nevertheless, QM/MM with an in situ optimised basis yields the most accurate results under all metrics, except for the slope (where the differences between basis set qualities are marginal). For H₂O, the effect of using an optimized basis is dramatic, which is expected, since we anticipate the orbitals in a molecule to be poorly described with small basis sets. Consequently, a minimal fixed basis (SZ) yields entirely wrong results, consistently predicting large and positive interaction energies, with an MSE as large as 31.2 kcal/mol. The addition of polarization functions improves results dramatically, but convergence with the size of the basis set is slow – DZP yields an MSE of 4.2 kcal/mol, and TZP yields 3.6 kcal/mol, with correlation coefficients of only ∼0.8. Only when in situ optimized orbitals are used do the results improve markedly – the MSE falls below 1 kcal/mol, the correlation coefficient exceeds 0.95, and maximum error diminishes by a factor of 2.7 compared to TZP. The optimized orbital formulation of QM/AMOEBA is a clear winner in this case, underscoring the advantages of using an in situ optimized basis set for the QM subsystem.

Having examined the energetics of our model, in the last step of our analysis, we will verify how the QM subsystem (solute) is affected by the presence of the QM/MM interface. Following the rationale of Sec. IV C, we will compare solute dipole moments as a proxy for the shape of the electronic density. As reference, we use fully QM calculations, where we partitioned the dipole moments into atomic contributions using DMA.⁹⁸ We compare the magnitudes of the total solute dipole between the two QM/MM approaches (where the solute dipole is simply the total QM dipole) and the two purely MM approaches [where the solute dipole is directly obtained from atomic dipoles (including induced dipoles in the case of AMOEBA) and charges].

We begin our discussion with the water pentamer system. While an isolated water molecule is charge-neutral, and its dipole moment is position-independent, in fully QM calculations of a solvated water molecule, we observe moderate charge transfer (below ±0.1e) between the central H₂O molecule and the solvation shell, which makes the solute dipole position-dependent, and necessitates choosing a reference point. For consistency with the rest of this work, we chose the centroid of the solute as the point, where we evaluate the dipole moment. The total dipole moments, as well as individual atom-centered dipole moments (for O and H atoms separately) are plotted in Fig. 8 and summarized in Table V.

Qualitatively, the total dipole moment of the central H₂O molecule (solute) is described similarly by both MM and QM/MM approaches – with substantial scatter and overpolarization (of about 0.6 D) relative to the purely QM results. We attribute this to the differences in the description of the four surrounding H₂O molecules – here only the reference uses a DFT description, while both MM and QM/MM approaches use a classical description. Given that standard DFT GGA models are known to struggle to correctly describe the structure and properties of water (cf., e.g., Refs. 90–92), and that we expect the dipole moment to be larger than the gas phase value (1.85 D) and lower than the value for bulk water (≈2.7 D),⁹⁹ we believe that it is in fact the reference calculation that underpolarizes the solvated H₂O molecule. Neither of the MM or QM/MM models seems to have a particular advantage in this case although QM/AMOEBA correlates slightly better with the reference.

It is more interesting to examine the individual atom-centered dipoles. Here we expect MM and QM/MM results to differ substantially because both MM approaches are disadvantaged by the constraint of fixed charge on the atoms, while in QM/MM, the charge density is free to transfer between atoms. GAFF, being non-polarizable, additionally yields zero atom-centered dipoles by construction. Thus, it is only meaningful here to benchmark the two QM/MM approaches against one another. For the O atom, we find QM/AMOEBA to be superior to QM/GAFF in terms of correlation with the reference result and a lower rms error (0.059 D vs 0.099 D). For the H atom, QM/GAFF does not reproduce the change in the dipole moment between the snapshots at all, while QM/AMOEBA shows the right trend although it mostly overpolarizes.

We finish with an examination of the two ionic solutes (Fig. 9, Table V). We evaluate the dipole at the position of the ion. For Cl⁻, we find all approaches, except for the non-polarizable GAFF, to generally overpolarize the ion, which is again expected, given the differences between classical and DFT treatments of the water solvent. AMOEBA is the least accurate possibly because of the large magnitude of charge penetration error for a diffuse Cl⁻ ion. QM/AMOEBA performs better than QM/GAFF (an rms error of 0.330 D vs. 0.584 D). For Na⁺, which barely polarizes at all (μ < 0.1 D), QM/AMOEBA turns out be the most accurate (an rms error of 0.007 D vs. 0.013 D for QM/GAFF) although all models are qualitatively correct.

We conclude that, at least for small hydrated solutes, our model is superior to QM/GAFF not only in how the QM subsystem is affected by the QM/MM interface but also in better energetics across the interface.

We presented and benchmarked a new mutually polarizable QM/MM model, where the QM subsystem is described using DFT with in situ optimized, localized orbitals (non-orthogonal generalized Wannier functions, NGWFs), and the MM subsystem is described using the AMOEBA force field. By implementing our model in the onetep linear-scaling DFT framework, we pave the way for affordable large-scale QM/MM calculations, with QM subsystems spanning thousands of atoms. However, in this work, we only studied small QM subsystems (up to 40 atoms), which are outside the linear-scaling regime.

The rationale for optimizing NGWFs is the near-complete-basis-set accuracy that they enable even with a minimal basis. This high accuracy is comparable or superior to even very large bases with fixed orbitals.⁶⁸ As part of this work, we derived and implemented the necessary gradients of the total QM + MM energy with respect to the NGWFs, enabling their in situ optimization also in the context of QM/MM calculations. We demonstrated how the additional flexibility of an in situ optimized basis exacerbates known problems of polarizable QM/MM methods – catastrophic overpolarization of the QM region (particularly in the presence of MM ions) and unphysical charge transfer (“charge spilling”) from QM to MM sites (particularly in the presence of MM cations).

We developed, presented, and validated conceptually straightforward solutions to both of these issues. We demonstrated that QM overpolarization can be mitigated by modifying the value of the Thole damping parameter only for polarization interactions spanning the QM/MM interface. This suitably attenuates QM/MM polarization at a very short range while having negligible medium-range and long-range effect. We addressed the charge spilling by introducing a more refined model of QM/MM Pauli repulsion interactions. This refinement replaces classical Halgren vdW repulsion with an electrostatic repulsive potential originating on MM atoms, which is parameterize to mimic Pauli repulsion. This approach is functionally equivalent to a density-overlap-based Pauli repulsion energy model with fixed, species-dependent densities placed on MM atoms, and, crucially, actual QM densities. This formulation is sensitive to the electronic degrees of freedom, which prevents the electronic density from excessive spilling.

Our modified approach requires parameterizing the repulsive potential with two values per MM species. We demonstrated how suitable values could be determined for K⁺ and Cl⁻ ions and for the atomic components of water, giving us confidence that this can be performed in principle. We did not identify a simple relation that would enable us to easily derive the sought parameters from classical vdW parameters, but we plan to investigate this further in future work.

We performed extensive tests to evaluate the transferability and reliability of our model with focus on the MM treatment of water. Using a variety of molecules, from small ions and neutral systems to larger molecules (up to 40 atoms), we showed that our model is, in general, although not universally, superior to nonpolarizable QM/MM and to purely MM approaches. This was reflected in lower disruptive effect of the QM/MM interface on the QM subsystem (which we assessed by comparing the dipole moments against a fully QM reference) but also by better energetics, when calculating interaction energies between the QM and MM subsystems. We find the latter particularly promising, as QM/MM models are typically very poor at describing interaction energies across the QM/MM boundary.

The authors are grateful to Alex Albaugh and Richard Bradshaw for fruitful discussions. AMOEBA parameterizations and initial structures of solutes used in this work were kindly provided by Richard Bradshaw. Structures of water dimers and solvated water molecules and ions were kindly provided by Yuezhi Mao. 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. and M.H.-G. acknowledge the support by Grant No. CHE-1665315 from the U.S. National Science Foundation. The majority of the calculations in this work were carried out at the TASK Academic Computer Centre in Gdansk, Poland.

In this section, we elaborate on the QM/MM polarization catastrophe (cf. Sec. III B) that can arise when NGWFs in the QM subsystem are optimized in situ. We use an H₂O:Cl⁻ system as an example. Figure 10 shows the magnitudes of atom-centered dipoles computed from the DMA procedure for QM atoms and the induced dipole of the MM Cl⁻ ion (whose permanent dipole is of course zero) in the course of SCF optimization. A mutual positive feedback can be observed to intensify at about step 80, quickly leading to absurdly large dipole moments (in excess of 1000 D). The sharp, step-like changes to the dipole values correspond to NGWF optimization steps and the relatively flatter parts of the graph – to the density kernel optimization steps. The classical degrees of freedom (MM induced dipoles) are fully optimized for each energy evaluation. The expected dipole moments on all atoms (as calculated from a fully QM reference calculation), shown with dashed lines, are well below 1 D. It is clear that the MM site, in particular, is overpolarized from the very first step of the optimization.

The main underlying reason for this is the inaccuracy of the multipolar approximation of QM density at short range, i.e., the charge penetration error (in this system, the Cl⁻ ion is 1.8 Å away from the H atom, while localized NGWFs extend for 3.7 Å). This can be appreciated in Fig. 11, where we plot the relevant component of the electric field due to the QM subsystem along the line joining the MM Cl⁻ ion and the leftmost H atom in the QM subsystem.

At the first glance, the field from the full electronic density (solid red line) seems to agree rather well with the field from the multipole approximation (dashed red line), up to ≈1.5 Å, where the multipole expansion starts to diverge. However, this field is to a large degree cancelled out by the field of the QM core, making the relative error in the total (blue) much more pronounced. Additionally, the shoulder to the left of the Cl⁻ ion (a result of a small fraction of the electronic density being attracted there by the electrostatic potential dip from Cl⁻’s induced dipole) cannot be well-represented by the multipole expansion and contributes to the charge penetration error. The issue is compounded by the fact that NGWF optimisation is driven by a gradient expression [cf. (22)] that, for consistency with a multipolar energy expression, has itself been derived via an intermediate step of a multipolar expansion. In consequence, the NGWF gradient is only sensitive to the MM potential, field and field derivative at the QM atom centres, rather than in the entire localization sphere [cf. comment directly below (22)].

One way to avoid the QM/MM polarization catastrophe is to simply reduce the NGWF localization region (e.g., we found 3 Å to be sufficient) – this makes the orbitals less diffuse and the point-multipole approximation more accurate. However, this would sacrifice some accuracy in the QM calculation – onetep calculations typically use localization radii of 3.5 − 5 Å.

What we propose instead is to slightly attenuate QM/MM polarization interactions, leaving permanent QM/MM interactions and MM/MM polarization unchanged. We retain the Thole functional form of the damping, but we reduce Thole’s a parameter by a factor of 2.45, which has the effect of attenuating the interactions at a very short range, having negligible effect elsewhere. This modification can be affected without any changes to onetep or tinker simply by rescaling the apparent polarizabilities of QM atoms as seen by tinker by a square of the above factor. The value of 2.45 has been found by numerical experiments on several small QM/MM systems, and we do not claim it to be optimal (indeed for the system studied here it leads to slight underpolarization, as seen in Fig. 12).

The practicability of the proposed solution is demonstrated in Fig. 12, which shows that all dipoles now converge to reasonable, finite values, and in Fig. 13 which shows how the total QM (electronic + core) electric field is now much better approximated by point multipoles, up to well below 1 Å. The unexpected charge transfer from QM to the left of the Cl⁻ ion seen in Fig. 11 disappears, owing to the dipole on Cl⁻ now being well-behaved.

In this section, we elaborate on the unphysical charge transfer from QM to MM that can manifest when QM/MM Pauli repulsion is not adequately taken into account. We will demonstrate that our improved model (cf. Sec. III C) addresses this issue satisfactorily. We use an H₂O:K⁺ system as an example.

Figure 14, panel (a) shows an isosurface of the electronic density at one point in the SCF optimization, where the unphysical charge transfer is apparent. The density shown is not the converged density, as the calculation fails to converge. This is because the spilled electrons accumulate near the peripheries of the orbital localization regions, and a localized-orbital formulation of DFT that assumes the orbitals to be well-decayed by the time they can be truncated cannot cope well with this situation. This is illustrated in Fig. 15, where a radial cross section of one of the NGWFs on the QM oxygen atom (red curve) is seen to differ markedly from its counterpart in a fully QM calculation (black curve).

Once a more physically sound model is used for QM/MM Pauli repulsion interactions (what we propose in Sec. III C), the charge transfer is prevented [Fig. 14, panel (b)] and the orbitals reacquire the correct shape (Fig. 15, green curve). The parameters used in this demonstration were $AK+=230Ha/e$ and $ζK+=1.379a0−1$⁠.