Exploring optimization strategies for improving explicit water models: Rigid n-point model and polarizable model based on Drude oscillator

Summary

Rigid n-point water models are widely used in atomistic simulations, but have known accuracy drawbacks. Increasing the number of point charges, as well as adding electronic polarizability, are two common strategies for accuracy improvements. Both strategies come at considerable computational cost, which weighs heavily against modest possible accuracy improvements in practical simulations. In an effort to provide guidance for model development, here we have explored the limiting accuracy of “electrostatically globally optimal” n-point water models in terms of their ability to reproduce properties of water dimer—a mimic of the condensed state of water. For a given n, each model is built upon a set of reference multipole moments (e.g. ab initio) and then optimized to reproduce water dimer total dipole moment. The models are then evaluated with respect to the accuracy of reproducing the geometry of the water dimer. We find that global optimization of the charge distribution alone can deliver high accuracy of the water model: for n = 4 or n = 5, the geometry of the resulting water dimer can be almost within 5⁰ of the ab initio reference, which is half that of the experimental error margin. Thus, global optimization of the charge distribution of classical n-point water models can lead to high accuracy models. We also find that while the accuracy improvement in going from n = 3 to n = 4 is substantial, the additional accuracy increase in going from n = 4 to n = 5 is marginal. Next, we have explored accuracy limitations of the standard practice of adding electronic polarizability (via a Drude particle) to a “rigid base”—pre-optimization rigid n-point water model. The resulting model (n = 3) shows a relatively small improvement in accuracy, suggesting that the strategy of merely adding the polarizability to an inferior accuracy water model used as the base cannot fix the defects of the latter. An alternative strategy in which the parameters of the rigid base model are globally optimized along with the polarizability parameter is much more promising: the resulting 3-point polarizable model out-performs even the 5-point optimal rigid model by a large margin. We suggest that future development efforts consider 3 -⁠ and 4-point polarizable models where global optimization of the “rigid base” is coupled to optimization of the polarizability to deliver globally optimal solutions.

Keywords:

Simulation and modeling – Oxygen – Biochemical simulations – Hydrogen bonding – Monomers – Electrostatics – Dipole moments – Dimers

Introduction

Water molecule has a deceptively simple structure H₂O, yet many anomalies of liquid water are still hard to explain despite the countless studies [1–5]—not surprisingly, water in its liquid state is notoriously difficult to model. Many complex biomolecules (DNA, RNA, proteins, etc.), vital for a living cell’s function, cannot be studied alone without considering their environment—water as the solvent [6]. To study those large biomolecules, atomistic simulations have been widely used, and numerous different water models [7] have been developed to reproduce water properties, including the class of water models most widely used today—simple, rigid, fixed-charge explicit models such as TIP3P [8], TIP4P [9, 10], TIP5P [11] and SPC/E [12]. According to the convention, these classical water models are distinguished by the number of “points”—interaction sites—in them: 3-point models (Fig 1a, with 3 point charges) such as TIP3P are the most common due to what is perceived by many as an acceptable balance between accuracy and computational cost (compared with 4 or 5-point models). For water models with more points, 4-point models (Fig 1b, with 3 point charges and neutral oxygen) such as TIP4P, and especially 5-point models (Fig 1c and 1d, with 4 point charges and neutral oxygen) models such as TIP5P (Fig 1d), the cost considerations become significant, even though these models tend to reproduce water properties better than their 3-point counterparts. However, despite decades of effort by many groups, none of the existing simple water models is perfect. [7, 13–17]

It is known that liquid water properties are determined by a complex network of hydrogen bonds. In these classical water models, hydrogen bonds are mimicked mostly by the electrostatics, complemented by a Lennard-Jones (LJ) potential. The latter—LJ potential—is generally represented by a single site centered on the oxygen, and its corresponding interaction is isotropic and featureless, in contrast to hydrogen bonding, which is directional. Thus, accurately representing electrostatic interactions is paramount for a classical water model to mimic hydrogen bond interactions and reproduce liquid water properties [19]. Therefore, our first question is how accurately can we describe the reference electrostatics within the unavoidable limitations of n-point models, and how far it can get us in terms of water model accuracy, by optimizing the electrostatics alone. This is a non-trivial question since, by construction, these simple n-point rigid models miss many important physics, such as molecular flexibility, electronic polarizability, and charge transfer.

We note that a good water model intended for biomolecular simulations must reproduce experimental properties of water with high accuracy. There are at least two reasons why existing water models are far from perfect in this respect. One is the limitations of the optimization strategies used in constructing water models. The other reason is on a more fundamental level—the limited physics built into these simple models. As an example of the first kind of limitation, current widely used rigid 3-point models (TIP3P, SPC, SPC/E, etc.) are based on an assumption that the experimental water molecule geometry is somehow optimal, or near optimal, for a classical water models—consequently, these models place the point charges on, or near, the centers of hydrogen and oxygen atoms. While sophistication of the optimization techniques employed to find the optimum has grown tremendously [20], from essentially “guess-and-test” to the complex, state-of-the-art force balance optimizations [21], one crucial aspect of the over-all procedure has not changed until recently: the search for best fit model is performed in the vicinity of the “canonical” water geometry, thus returning a local optimum in the parameter space. It was recently demonstrated [19] that abandoning the restrictions (except fundamental C_2v symmetry) on water model geometry, and performing an exhaustive search for a truly global optimum in the parameter space can result in significant accuracy improvements. The resulting 4-point rigid model—OPC [19] was built without the geometry restriction on point charge placements. It was optimized globally to reproduce bulk water properties as best as possible, without any increase of the computational cost of employing the model in simulations, relative to common 4-point models such as TIP4P-Ew.

Apart from using global optimization, a natural strategy to improve the accuracy of a rigid n-point water model is to consider larger “n”. The key question with this approach is whether the inevitable and substantial increase in computational cost of employing “larger n” water models is justified by significantly better accuracy? If the accuracy gain is large enough, the cost increase may still be well worth it. Answering the question is important both for model developers and practitioners alike. Unfortunately, it is nearly impossible to compare existing n-point water models on the same footing to find out exactly how much gain an increase in n brings. For once, these models are optimized using different criteria, and against non-identical set of water properties. For example, TIP5P-Ew model [22] surpasses many of its 3 -⁠ and 4 -⁠ point predecessors in accuracy of describing water structure, but its accuracy of reproducing the self-diffusion is lower than that of TIP4P-Ew. In addition, optimization protocols vary widely, and, except for OPC-family models, there is no guarantee that the model corresponds to the global optimum in the optimization landscape. Thus, it is entirely possible that a global optimum for a smaller n may yield a more accurate model than a local one for a larger n.

Even if a hypothetically perfect fixed-charge rigid model reproduced a large subset properties of bulk liquid water exactly, the model would still be inherently incapable of responding to the change of polarity of its micro-environment, relevant to biomolecular simulations. Specifically, water is highly polarizable [23]: the experimentally observed change in the dipole moment of real water molecule upon transfer from gas to liquid phase is as large as 1 Debye, while for any rigid fixed-charge model that change is zero by construction. The polarity of micro-environment near a macromolecule can be quite different from that of bulk water; non-polarizable models cannot properly respond to different micro-environments during the course of a simulation [24], e.g. in cross-membrane transport. However, the vast majority of current water models lack polarizability, which must adversely affect the accuracy of biomolecular simulations. This example illustrates the second type of limitation on water model accuracy—the missing physics. To address the lack of electronic polarization, a number of polarizable water models have been developed, for example POL1 [25, 26], FF12POL [27], SWM model family [28] and AMOEBA model family [29–31], see recent reviews for a comprehensive account of the field [7, 32, 33]. There is little doubt that, in principle, availability of highly accurate and efficient polarizable models for routine practical simulations should improve the accuracy of biomolecular simulations significantly, by accounting for effects completely missed by rigid fixed-charge models. However, because of the unavoidable higher complexity of polarizable models relative to the fixed-charge ones, the balance towards wider adoption can only be tipped by distinctly higher accuracy, which does not seem to be the case yet [27, 34]. Since current polarizable water models typically use existing fixed-charge models as their “base”, they likely inherit at least some of the existing flaws of the fixed-charge models, unrelated to their lack of polarizability. That is they are likely represented by local optima in the complex optimization landscape. Indeed, despite undeniably better physical foundation, accuracy of a sophisticated polarizable water model can be lower than that of a globally optimal non-polarizable rigid model [19]. This observation begs the question: how much accuracy gain can the inclusion of polarizability bring, if the above optimization-related limitations were removed, and a truly globally optimal polarizable model was constructed? In other words, is global optimization a potentially useful strategy for polarizable models? For reasons discussed above in the context of comparing non-polarizable models, comparing existing n-point non-polarizable with polarizable models (even with respect to liquid water properties alone) on the same footing is difficult. Moreover, to the best of our knowledge, a truly globally optimal polarizable model has not be constructed yet.

The main motivation of this work is to address the questions outlines above in a tightly controlled computational experiment. To this end, here we build and evaluate several test water models of two types—rigid n-point globally optimal models (n = 3, 4, 5) and also polarizable models based on these rigid models—that illustrate the two limitations discussed above. We also construct a truly globally optimal polarizable test model. All models are optimized using the same type of protocols, and to the same accuracy level in reproducing the training set. The models are evaluated by examining their ability to reproduce water dimer, which can be considered as the simplest possible mimic of water in condensed state.

Methods

As mentioned above in the Introduction, liquid water properties are determined by a complex network of hydrogen bonds, which is the most challenging part a water model needs to mimic in order to simulate real water. Water dimer, involving only two water molecules and their interactions, is a good minimal starting point to reveal the secrets of hydrogen bonding, and has been studied over half a century [35, 36]. Due to its significance to hydrogen bonding studies, and being the simplest possible mimic of the liquid phase water, the water dimer structure is utilized here to evaluate the quality of the water models.

Overall approach

The overall approach is as follows. Individual water model’s geometry and values of point charges are optimized to match as closely as possible a set of reference multipole moments (Table 1, below). The matching is done in the following precise sense. Through optimal point charge approximation (OPCA, a method to approximate electrostatic charge distributions with a small number of point charges to optimally represent the original charge distribution) [37], each water model has its quadrupole and octupole moments fitted to the reference, and the monomer dipole moment is optimized so that when the two water models form a dimer, the dimer total dipole has a smaller than 0.1% error relative to ab initio calculation—2.68D [38]. Then the water model’s accuracy is evaluated by comparing its dimer’s geometry (Fig 2) with the one obtained from an ab initio calculation [39].

As discussed in the introduction, a water model’s accuracy of simulating liquid water is largely determined by its electrostatic properties, small variations in these properties can lead to large differences in simulation outcomes [40]. Also, the “electrostatic” parameters are where different types of water models differ the most from each other, while the Van der Waals potential is typically simulated in the same manner: by a Lennard-Jones (LJ) site on oxygen. All of model types studied here follow that approach. In reality, optimization of the LJ parameters can be largely decoupled from optimization of the “electrostatic” parameters such as partial charges and their distribution [19]. An accepted approach to construct realistic water models involves optimizing the LJ parameters against reference O-O RDF (radial distribution function), to obtain a close match [19]; we follow essentially the same strategy here, except that we use ab-initio reference oxygen-oxygen distance (r(OO) = 2.91Å [38, 39]) as the target reference for optimization of LJ parameters. With the r(OO) distance used in LJ optimization, the remaining parameters of the dimer geometry—the two angles (θ and φ in Fig 2)—are left as convenient metrics to evaluate the model accuracy, which is mainly determined by the electrostatics.

From the perspective of electrostatics, a model with larger n (number of interaction points) is expected to better reproduce given multipole moments, and thus improve the accuracy of the model [41–44].

For a set of N point charges, in this case—the water model, the Coulomb potential can be written as:

In a Cartesian system, this equation becomes:

Because of the C₂ symmetry of a water molecule, the quadrupole moments and the z-plane elements of octupole moments can be written as:

With the equations all set up, we can then rearrange them to the form where water model parameters (coordinates and charge) are explicitly expressed with multipole moments [19].

For a 3 or 4-point water model, the rearranged equations become:

A 4-point water model is shown in Fig 3b. With the C₂ symmetry, only 4 parameters are needed to construct the 4-point water model—y₁, z₁, z₂ and electric charge q. Both 3 and 4-point models have 3 point charges, so the calculation process are the same (Eqs 9–12) except for 3-point models the condition z₁ = 0 needs to be imposed (Fig 3a).

For a 5-point model, there are 5 parameters to be determined (q, x₁, y₁, z₁, z₂, see Fig 3c).

From Eqs (3)–(6), we have:

With Eqs (13)–(16), a series of expression (x₁, y₁, z₁, z₂) = f(q, μ, Q_t, Q₀, Ω_t) can be derived:

After choosing a set of multipole data (μ, Q_t, Q₀, Ω_t, Ω₀), insert μ, Q_t, Q₀, Ω_t values into Eqs (18)–(21) so that we can then express the coordinates (x, y, z₁, z₂) with one argument—q the charge value. And now, what is left is a very simple optimization problem: finding the optimal q that makes the Ω₀ error of the model as small as possible with exhaustive search. Then with the optimal q and its corresponding coordinates (x, y, z₁, z₂), a 5-point water model is then constructed and ready for evaluation.

Evaluation of the built water models starts with the energy minimization calculation of two identical water molecules that form a water dimer, using AMBER 2019 [45]. Steepest descent, conjugate gradient method (for 3, 4-point rigid models) and Limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) quasi-Newton algorithm (for 5-point rigid models and 3-point Drude models) are used in this minimization calculation. With the energy minimization process, the LJ parameters on the oxygen of each model are determined so that the minimized water dimer has an oxygen-oxygen distance matching the reference (2.91Å). With the oxygen-oxygen distance controlled, the dimer system has a very limited number of variables (rotation angles), therefore finding the global minimum is straightforward. To assert it, we ran multiple minimizations with random starting coordinates, and found that all the minimized states were the same. Thus, after minimization, the dimer system of the constructed water model reaches its global lowest energy state. The total dipole moment of the minimized water dimer is then calculated and compared with the ab initio result(2.68D [38]). If the relative error of the total dipole moment is greater than 0.1%, the monomer dipole moment will be changed, we then redo the model building and the dimer minimization. This optimization process is continued until a monomer dipole value is found when the total dipole moment of the minimized dimer has a relative error within 0.1%.

As a metric to evaluate the geometry of water dimer, two parameters are used—θ and φ. θ and φ are the angles between water molecule’s C₂ axis and the O-O line for the two molecules respectively. θ is of the water molecule coplanar with O-O, φ is of the other one non-coplanar with O-O (Fig 2). We use an ab initio calculation of a water dimer as the reference geometry, where θ₀ = 57.9°, φ₀ = 55.6° and r(OO) = 2.91Å. [38, 39] This ab initio calculation is done by Klopper et al. in 2000 [39], optimized at the level of CCSD(T) theory with IO249 basis. It is seen as a benchmark result of theoretical studies on water dimer. [36]

Reference multipole moment sets

To ascertain robustness of our methodological conclusions, we use three different sets multipole moments as references to which we fit our water models (Table 1). The monomer dipole, quadrupole and octupole moments of gas phase QM calculation (CCSD with aug-cc-pCVTZ basis) are from Anandakrishnan et al [37]. The gas phase experimental data are from Clough et al [46]. When used to construct water models, the octupole moments from the gas QM set are used in the gas phase experimental set (whose octupole moments are not available). The third set of reference multipole moments is of liquid phase water, a QM calculation done by Niu et al [47].

Each set of reference multipole moments is used to construct a set of optimal water models comprising of 3-point, 4-point and 5-point rigid models. These three sets of water models are called “gas exp.”, “gas QM” and “liquid MP2/4MM” respectively, based on the dataset used for model construction. The “gas exp” reference multipole moments are also used to construct polarizable models in this study.

n-point rigid models

Optimized against each of the three multipole moment data sets (“gas exp.”, “gas QM” and “liquid MP2/4MM”), three rigid models were constructed with 3, 4 and 5 interaction points respectively, following the process shown in Fig 4. Parameters of these rigid models are listed in Table 2.

For each n-point rigid water model built, its monomer dipole moment μ value is varied for optimization. When the monomer dipole moment changes, the resulting water model parameters (coordinates and charges) change, as a result, the water dimer of this model will be different. An optimal monomer dipole moment is chosen when the total dipole of the corresponding dimer is within 0.1% relative error from the ab initio calculation [38].

Polarizable water model with Drude oscillator

There are several different approaches available to add electronic polarizability to a “rigid base” water model, adding the Drude oscillator [48, 49] being one of the simplest methods. The Drude oscillator is easy to implement in most MD simulation packages. Critically, it only adds one more degree of freedom for optimization, which will not impair the search for the global optimum. Compared with other approaches, the Drude approach also has the benefit of high computational efficiency [32, 50–54]. Here the polarizable water model is constructed by adding a Drude particle to the oxygen in a 3-point rigid water model (Fig 5, the Drude particle has no interactions with the 2 hydrogens in the same molecule). It is referred to as a “3-point polarizable model”. The polarizability of this model is

Based on “gas exp.” data set (Table 1), we built two polarizable models: the “limited optimal” one with only polarizability optimized for fair comparison with rigid models; and the “globally optimal” one with both polarizability and gas phase monomer dipole moment optimized for practical evaluation.

Optimizing the 3-point polarizable water model

For the limited optimal 3-point polarizable model, we fixed the monomer dipole moment to experimental value–1.86D (thus the rigid base is fixed), and optimized it by varying polarizability until its dimer dipole reproduces experimental result with in 0.1% error, the procedure is shown in Fig 6. The optimization and the resulting model are shown in Fig 7 and Table 3, respectively. The resulting polarizability of the limited optimal polarizable model is 0.94ų, 35% lower than the gas phase experiment value–1.44ų [55].

For the globally optimal 3-point polarizable model, we modify the model by changing two input parameters—monomer dipole moment μ and polarizability α (both in the vicinity of their gas phase experimental value, μ_exp = 1.86D and α_exp = 1.44ų), to make the model’s dimer total dipole moment fall within 0.1% error from reference [39]. Error distributions of the intermediate models during this optimization process are shown in Fig 8, with respect to the two varying parameters.

The resulting globally optimal 3-point polarizable model is constructed with polarizability α = 2.2900ų and monomer dipole moment μ = 1.7258D (Table 3). This model has a total angle error of 3.19° when simulating a water dimer. The monomer dipole moment (1.7258D) of this polarizable model is relatively close to the gas phase experimental data (1.86D), while its polarizability (2.275ų) is larger than the gas phase experimental value (1.44ų).

Results and discussion

For biomolecular simulations, the most important state of water is the liquid state at ambient temperature and pressure. However, the number of experimental parameters used to characterize liquid water is vast, and existing practical procedures used to build water models to match a subset of these parameters are highly diverse, non-trivial and expensive [7]. The process is still somewhat of an art. For a strict, limiting-case apples-to-apples comparison we want to present here, using liquid state properties is thus not ideal, if not completely impossible. We therefore resort to using a mimic of condensed state of water—water dimer. Water dimer is the smallest form of a water cluster. It has a hydrogen bond, which greatly contributes to numerous anomalies in water properties. As a starting point to study water properties, water dimer has been used extensively [35, 36, 56]. It has also been employed as a standard reference to test water models [8, 9, 11, 12, 19, 27–30].

Our over-all strategy is as follows (see Methods for details). For n-point rigid models, we build each model to match, as closely as possible, the reference quadrupole, octupole moments of water monomer (Table 1), and a monomer dipole moment optimized for water dimer. For polarizable models, its rigid “base” is built first to match the reference multipole moments including dipole, and then the polarizable model is built by adding a Drude oscillator with an optimal polarizability value—optimized for water dimer—to the rigid “base”. Parameters of the “base” are not re-optimized. The optimization of the last 3-point polarizable model is slightly different: it involves varying both the monomer dipole and the polarizability, such that both its “base” and Drude oscillator are optimized. Regardless of what parameters are optimized, the standard of the optimization is the same for all these models—the total dipole moment of the water model’s dimer matches the ab initio reference value for water dimer dipole [39] within a given tolerance.

We then test each model’s ability to reproduce the geometry of a water dimer formed by the two water model molecules. The reference dimer geometry is from the benchmark ab initio calculation by Klopper et al. [38, 39] For reasons detailed in “Methods”, we keep the O-O distance in the dimer fixed to its reference value, thereby excluding oxygen LJ parameters from the consideration.

Rigid n-point water models

First we explore accuracy of rigid non-polarizable water models with 3, 4 and 5 point charges. The model parameters are given in Table 2. Three sets of n-point (n = 3, 4, 5) rigid water models are included. These three sets of water models are constructed based on 3 different sets of multipole moments data, see “Methods”.

Accuracies of these models with respect to reproducing our main reference are compared in Fig 9. For reference, TIP3P-, OPC3 -⁠ and OPC-based dimers are also tested and included in Fig 9. The water dimer geometry error values of the 4 -⁠ and 5-point models are considerably smaller than that of 3-point models. The remaining errors are at ∼5° level, well below the ±10° experimental error margin of the reference values [38]. Note that this is the result of optimizing the electrostatics only. Additionally, our models have smaller errors than the existing models (TIP3P, OPC3 and OPC, Fig 9), which is not unexpected because the latter ones were optimized for liquid water instead of the water dimer. This relatively high accuracy of the “electrostatics only” optimization emphasizes the critical importance of getting right the distribution of charge in a water model. The result also suggests that one can improve a rigid water model significantly by optimizing its electrostatics, even though some physics is still missing from the model.

The dimers of 4-point water model show notably better agreement with the reference data than the dimers based on 3-point model—error of the two key angles reduced to ∼5° from ∼20°. This improvement of accuracy from 3-point to 4-point model is robust, seen for all reference sets used to build these models. However, further increasing n from 4 to 5 only has minimal effect on improving the agreement with the reference dimer geometry. The 5-point models with “gas exp.” and “gas QM” data improve from their 4-point counterparts with marginal changes (<1°). The “liquid MP2/4MM” 5-point model even has a slightly inferior accuracy compared with its 4-point model (error is 1° larger). Thus, adding more points beyond n = 4, which increases the accuracy of reproducing multipole moments, n-point rigid models still cannot improve much from what a 4-point model can already do. As seen from Fig 9, all of the above conclusions are independent of the specific reference multipole set used.

Limited optimal 3-point polarizable model—only polarizability optimized

One obvious physical feature missing from rigid n-point water models is electronic polarizability [29, 49, 57–62]. This missing physics may be the reason why there is little, if any, improvement in the accuracy of the rigid model in reproducing water dimer as n is increased beyond n = 4, Fig 9. How much of an accuracy improvement can be achieved by adding the missing polarizability to the model? To answer this question, we compare two different approaches of improving rigid water models.

So far, in the rigid models we built, the geometry and point charge values are optimized (with OPCA method [37], by varying the monomer dipole moment of the water model) specifically for water dimer to compensate the missing gas-to-condensed-phase polarization. Now, we add polarizability to a 3-point model, whose monomer multipole moments (including the dipole) fit the reference for the monomer, then optimize the added polarizability against water dimer to the same standard—error of dimer total dipole moment smaller than 0.1%. This way, an apple-to-apple comparison is established between the two approaches to improve a rigid water model: 1) adding the missing physics of electronic polarizability; 2) compensating for its absence by optimizing the electrostatics. Either approach has only one parameter for optimization so they are competing on a level field.

Parameters of the resulting polarizable model are shown in Table 3. The accuracy of the resulting 3-point limited optimal polarizable model is compared to that of our rigid models in Fig 10.

As illustrated in Fig 10, this limited globally optimal 3-point polarizable model has ∼5° smaller total angle error than the 3-point rigid model. This means that under the same optimization standard, adding and optimizing polarizability (while keeping the monomer dipole moment fixed to its gas phase value) is more beneficial than compensating the lack of polarization by only optimizing the dipole moment of the rigid model.

Moreover, because the monomer dipole moment is fixed to the gas phase value in this limited optimization process, the resulting polarizable model reproduces correct gas phase multipole moments when in gas phase, unlike rigid models whose dipole moment is changed for condensed phase optimization. Also, the optimized polarizability is 0.94ų, 35% lower, but still qualitatively similar to the gas phase experimental value—1.44ų.

However, this limited global optimal 3-point polarizable model is still inferior to the 4 -⁠ and 5-point rigid models with quite a margin (∼10° higher total angle error). Therefore, adding polarizability alone cannot fix the defects of suboptimal electrostatics of a 3-point model.

With the computational cost in mind (in real simulations, 3-point polarizable model is even slower than the corresponding 4-point rigid models), this result is not very satisfactory. It suggests that the optimization protocol in which the fixed-charge base model is ignored and only polarizability is optimized has significant limitations.

Globally optimal 3-point polarizable model—Both the polarizability and the “rigid base” are optimized

So far, we have focused on purely conceptual questions about water models’ limitations in reproducing liquid water, mimicked here by water dimer. In practical MD simulations, water models are mostly used to simulate the solvent in liquid phase. Therefore from practical perspective of water model design, we do not have to keep the rigid base of a polarizable model unoptimized. Also, the results of the limited optimal 3-point polarizable model make it clear that only optimizing polarizability is not sufficient to improve a water model’s accuracy beyond the rigid n-point model’s limitations. Therefore, here we build a globally optimal 3-point polarizable model, in which two parts are both optimized for the water dimer scenario: the geometry and point charge values of the rigid base (through OPCA [37], varying the monomer dipole); and the polarizability introduced by the Drude oscillator.

As discussed in “Methods” section, under the same optimization standard (reproducing dimer total dipole) a polarizable model has one more degree of freedom—the added polarizability—than the rigid models. By optimizing both its rigid base and its polarizability, we are utilizing this additional degree of freedom afforded by the polarizable model. In other words, this 3-point polarizable model is optimized with two parameters: polarizability(α) and gas phase monomer dipole moment(μ). The result is shown in Fig 10. In comparison with the 3, 4, 5-point rigid models and the limited optimal polarizable model, this globally optimal polarizable model’s result is shown along side with them as the rightmost column.

With both polarizability and gas phase monomer dipole moment optimized, the dimer geometry error of this fully optimal polarizable model shows a significant improvement from the limited optimal polarizable model, and is more accurate than the corresponding 4 -⁠ and 5-point rigid models—the φ error is very close to zero and the θ error is at the same level or smaller than that of all other models. Not only the accuracy of the globally optimal polarizable model is considerably higher than that of all rigid models and the other limited optimal polarizable model, it is worth noting that this globally optimal polarizable model has a gas phase monomer dipole moment of 1.7258D, not far from the gas phase experimental value of 1.86D. The resulting polarizability is 2.29ų, which is considerably larger than the experimental 1.44ų. This deviation of the optimized polarizability value points to remaining deficiencies of the specific polarizable water model, likely that the model is too simplistic or/and that other physical effects not explicitly considered here (e.g. charge transfer) play a noticeable role. For the 3-point polarizable water model optimized for water dimer, decreasing the monomer dipole moment by 7% (from 1.86D to 1.7258D) resulted in a 144% polarizability increase (from 0.94ų to 2.29ų). This implies that the water dimer properties are less sensitive to the model’s polarizability than to its monomer dipole moment, which agrees with the study of Soetens and Millot [63]. Although the water dimer properties are less sensitive to the water model’s polarizability, the polarizability is still very important in improving the model’s accuracy, since a 3-point polarizable model outperforms a 5-point non-polarizable model by quite a margin.

The combined results show the full potential of a polarizable water model when properly optimized. Under the same optimization standard (reproducing dimer total dipole) and with the additional degree of freedom utilized, the 3-point polarizable model outperforms 3, 4 and even 5-point rigid models by achieving the smallest error in reproducing water dimer geometry, while still having electrostatic properties in gas phase close to experimental values. At the same time, even this globally optimal polarizable model is not perfect.

Conclusion

In order to study possible novel avenues for optimization strategies aimed at improving accuracy of explicit water models for biomolecular simulations, different “toy” water models were constructed and examined with respect to their ability to reproduce properties of water dimer—a mimic of the condensed state of water. Specifically, we constructed rigid models with 3, 4 and 5 points, and two different 3-point polarizable models with a single Drude particle to mimic electronic polarizability The models were built to match reference multipole moments, including ab-initio, and were then optimized to reproduce, as closely as possible, the total dipole moment of water dimer. The ability of the models to reproduce the water dimer geometry was used as the metric of the model accuracy.

First, we conclude that optimizing the “electrostatics” (charge distribution) alone can deliver high accuracy of the water model: the geometry of the resulting water dimer is essentially within 5° of the ab initio reference, this remaining error is almost half of the reported experimental error margin (±10°) on the dimer geometry. Not unexpectedly, the resulting water models show smaller errors than the existing models (TIP3P, OPC3, OPC), which were optimized for liquid water instead of the water dimer. This result reinforces the notion that optimizing electrostatics is key to water model quality, and an “electrostatically globally optimal” water model can be quite accurate.

Second, we conclude that for rigid n-point water models, increasing the number of interaction points from n = 3 to n = 4 can easily bring better accuracy, while by further increasing n from 4 to 5 only a marginal improvement can be achieved. Considering the steep increase with n of the computational cost associated with performing simulations based on n-point water model, a very convincing justification is needed for renewed efforts to build rigid water models with 5 or even more point charges. Justification for such models may need to be based on expected or demonstrated improvements in specific areas, e.g. where the possibility that non-planarity of a water models may be critical, such as phase transition between ice and liquid water that TIP5P model has a melting point much closer to experiment than most 4-point models [64]. The above conclusions are robust to the specific type—out of three different ones—of the reference multipole set used to build the models.

Next, we have investigated the effect of adding electronic polarizability to the rigid base model. Our first exercise was aimed at quantifying how much extra accuracy one can gain by having a dedicated mechanism (Drude model) to account for the polarizability. To make an apples-to-apples comparison, we compared two optimization options that employed just one adjustable parameter. The first option is the method we already applied to build our rigid models where point charge placements are adjusted to compensate for the polarizability in an average sense; the second option is to add and optimize polarizability of the Drude oscillator without changing the rigid “base” model. The second option bears similarity to optimization strategies used in practice to construct Drude polarizable models. The same reference multipole moments and dimer geometry accuracy metric were used in both options. The results show that for the same 3-point “base”, having the dedicated polarizability component does achieve a better accuracy than the purely rigid model, which is an argument in favor of polarizable models, despite their extra complexity. At the same time, the accuracy of this limited optimal 3-point polarizable model falls significantly short of the optimal 4-point rigid non-polarizable model. Apparently, the correct physics added to the rigid model via electronic polarizability is mostly “wasted” on correcting the deficiencies of the 3-point “base” model.

To reveal full potential of polarizable models, we have explored a 3-point polarizable model in which both the base charge distribution and the polarizability are globally optimized simultaneously—to the best of our knowledge, this strategy has not yet been used to build existing polarizable models used in practical simulations. The resulting truly globally optimal 3-point polarizable model easily outperforms 4 -⁠ and 5-point globally optimal rigid models. Not only that the globally optimal polarizable model is more accurate in reproducing the water dimer, its unpolarized monomer dipole moment (1.7258D) is quite close to the correct gas phase value (1.86D), indicating that the model is a decent mimic for both the condensed phase and gas phase of water.

We stress that all of the water models developed in this work are not intended for use in actual atomistic simulations. Instead, these “toys” models serve only one purpose—to compare various model optimization strategies on the same footing and suggest avenues to consider for future practical optimization efforts to improve realistic water models. In the future, we plan to build and explore a truly electrostatically globally optimal polarizable model, following the over-all optimization strategy outlined in this work.

Supporting information

S1 Protocol [zip]
AMBER files for each water model we built in this study.