Shear viscosity of OPC and OPC3 water models Open Access

Molecular dynamics (MD) simulation is a powerful tool for understanding physical, chemical, and biological phenomena at the atomic level. To ensure the success of simulation studies, the quality of the model should be properly assessed by comparing the simulation results with experimental data. In the case of water models, their ubiquity and importance have led to the development of over 50 explicit water models in the past decade to reproduce the unique properties of water.^1–3 Even though water has a simple chemical composition consisting of merely two hydrogen atoms and one oxygen atom, reproducing all the experimental properties of bulk water simultaneously remains a challenging task.³ Furthermore, since most intracellular reactions involving biological molecules, such as proteins, DNA, and polysaccharides, take place in an aqueous environment, developing more accurate water models and conducting their quality assessments are also important for successful biomolecular simulations.² 

Currently, the Extended Simple Point Charge (SPC/E)⁴ and Transferable Intermolecular Potential three-point (TIP3P)^5,6 models are commonly employed as three-point rigid water models in MD simulations. For four-point rigid water models, the re-parameterized version of the original TIP4P water model,⁶ TIP4P/2005, is widely used, which reproduces various experimental properties.⁷ In particular, for biomolecular simulations, the TIP3P water model has been extensively used within the Amber,⁸ OPLS-AA,⁹ and CHARMM¹⁰ force fields for an extended period,⁹ wherein CHARMM uses a slightly modified version of the TIP3P water model.¹¹ The TIP3P water model accurately reproduces several key features of bulk water at 25 °C and 1 atm: the density, heat of vaporization, and the first peak of the oxygen–oxygen radial distribution functions.^6,9,12 However, for other properties of liquid water, such as the static dielectric constant, self-diffusion coefficient, the temperature of maximum density, and the second peak of oxygen–oxygen radial distribution functions, the performance of the TIP3P water model is relatively poor when compared to that of other recent water models.¹² Since the Amber and CHARMM force fields have been extensively optimized for the TIP3P water models to reproduce experimental data on the thermodynamics and dynamics of biomolecules, some of the deficiencies in the water model should be ameliorated.³ In 2014, Izadi, Anandakrishnan, and Onufriev introduced a nonpolarizable four-point rigid Optimal Point Charge (OPC) water model that exhibited significantly improved accuracy in reproducing a wide range of bulk water properties when compared to commonly used rigid water models across various temperatures.¹³ Notably, their developmental approach refrained from imposing any geometric constraints on the model, except for symmetry, contrasting with the conventional approach.¹³ In 2016, Izadi and Onufriev further introduced a three-point rigid OPC water model (OPC3), which also demonstrated competent performance in reproducing various properties of bulk water, although it was slightly less accurate than the four-point OPC.¹² They suggested that the accuracy limit of the three-point rigid nonpolarizable water model had been reached.¹² The three-point TIP3P-FB and four-point TIP4P-FB rigid water models developed by Wang et al. in 2014 using their automated parameter optimization procedure “ForceBalance” also accurately reproduce a number of experimentally measured physical properties of water.¹⁴ Ion parameters have been optimized for the OPC, OPC3, TIP4P-FB, and TIP3P-FB water models to reproduce their solvation thermodynamics.^15–17 Recently, the Amber development team has recommended the use of the OPC water model in combination with their latest biomolecular force field, Amber ff19SB, to enhance the predictive power for modeling sequence-specific behavior of proteins, protein mutations, and rational protein design.⁸ Tempra et al. also reported that the OPC water model accurately reproduces experimental surface tensions and is essential for realistic lipid monolayer simulations.¹⁸ Therefore, the use of the OPC and OPC3 water models is becoming increasingly common, particularly in the biomolecular modeling and simulation community.

To the best of our knowledge, the viscosities of OPC and OPC3 have not been reported. Thus, to add important information to the performance evaluation of the OPC family of water models, we calculated the shear viscosities of the OPC and OPC3 water models at various temperatures. The shear viscosity is also important for evaluating the self-diffusion coefficients of molecules in an infinite system from MD simulations under periodic boundary conditions since the value of the solvent is required to correct system-size effects in the calculation of diffusion coefficients.¹⁹ In addition, we evaluated other properties, such as the density, self-diffusion coefficient, surface tension, and static dielectric constant, for thoroughness and comparison with other related studies. To verify our evaluation protocol, these properties of the TIP4P/2005, TIP4P-FB, and TIP3P-FB water models were also calculated for comparison with the results in previous reports by others. A summary of these results is provided in the supplementary material.

Two systems were employed to simulate liquid water: 2000 and 256 water molecules in a cubic box with a side length of ∼40 and 20 Å, respectively. For the MD simulations of the former system, we used the GPU version of the pmemd program from Amber18²² and Amber22,²³ running on NVIDIA GeForce RTX 3090 graphics cards. For the latter system, the CPU version of the pmemd program was used. Periodic boundary conditions were applied in all directions. Long-range electrostatic interactions were calculated using the particle mesh Ewald method,²⁴ whereas short-range electrostatic and Lennard-Jones (LJ) interactions were evaluated using a cutoff method. All intra-water geometries were constrained using the SHAKE algorithm,²⁵ and a time step of 2 fs was used. The simulation system was equilibrated for 2 ns under isobaric–isothermal (NPT) conditions after energy minimization, where the temperature was controlled using a Langevin thermostat²⁶ with a collision frequency of 1 ps⁻¹ and the pressure was kept at 1 bar using a Berendsen barostat²⁷ with a relaxation time of 1 ps. After equilibration, a production run was performed for 10 ns under isochoric–isothermal (NVT) or microcanonical (NVE) conditions to evaluate the shear viscosity of water. For the NVE simulations, a time step of 1 fs was used to minimize energy drift. To control the temperature during the production run, the Berendsen weak-coupling method²⁷ with a relaxation constant of 1 ps was used to minimize the artificial effects of the thermostat algorithm on the transport properties of the liquid molecules.²⁸ Energies and pressure were saved every 0.01 ps. The Amber MD program computes only the three diagonal components of the pressure tensor, P_xx, P_yy, and P_zz, not the three off-diagonal components during NPT simulations using a Berendsen barostat. To compute the P_xx, P_yy, and P_zz during NVT or NVE simulations in the Amber programs, the functionality of the Berendsen barostat was turned on, but with a compressibility parameter of zero to prevent rescaling the coordinates and the box sides during the simulations. In viscosity calculations using Eq. (1), three independent shear components of the pressure tensors $Pxx−Pyy/2$⁠, $Pyy−Pzz/2$⁠, and $Pxx−Pzz/2$ were included to improve the statistical accuracy.

In this study, two different cutoff distances, 8 and 12 Å, as well as nine different temperatures, 273, 283, 293, 298, 303, 313, 333, 353, and 373 K, were examined. To enhance the statistical accuracy, 50 independent simulations were performed with different random seeds in the equilibration phase for each temperature. The analysis was conducted using cpptraj²⁹ and our in-house programs. Some of our analysis tools and input files for the MD simulations used in this study are available on GitHub (https://github.com/Ando-Group/Water_Viscosity_Calculation). The average and standard deviation were obtained from 50 independent simulations, and the latter is reported as a statistical uncertainty.

Figure 1 shows the cumulative integrals of the pressure tensor ACF, η(t), in Eq. (1) up to t = 100 ps obtained from the NVT simulations of 2000 OPC water molecules with a cutoff of 8 Å at four different temperatures. At all temperatures, η(t) reached a plateau within 10 ps and remained almost constant up to 100 ps. The shear viscosities calculated from the other simulation models had the same tendencies. Because a longer integration time increased the statistical errors, an upper limit of 10 ps was imposed on the ACF integral in this study.

Table I lists the shear viscosities for the OPC and TIP4P/2005 water models calculated from the simulations with 2000 and 256 water molecules under NVT and NVE ensembles at 273 and 283 K. The η values of the OPC water with 2000 and 256 molecules were very close within the statistical error at 273 and 283 K, respectively. Therefore, finite size effects are negligible in the viscosity calculations, which have already been reported by several groups.^19,30–32

The effects of temperature control algorithms on the viscosity calculations were investigated by Fanourgakins, Medina, and Prosmiti,³³ and Hafner et al.³⁴ Their studies showed that the shear and bulk viscosities estimated for the NVE ensemble are statistically indistinguishable from those obtained in the NVT ensemble using global velocity scaling thermostats (e.g., the Berendsen,²⁷ stochastic rescaling,³⁵ and Nosé–Hoover³⁶ thermostats) with a wide range of thermostat coupling strengths.^33,34 In contrast, local velocity randomizing thermostats (e.g., the Andersen³⁷ and Langevin²⁶ thermostats) with strong coupling drastically increase the viscosities compared to the NVE results.³⁴ Similar to their conclusion, our estimates of η in the NVT ensemble simulations using a Berendsen thermostat were also statistically indistinguishable from those in the NVE ensemble simulations for the OPC and TIP4P/2005 water models as listed in Table I. In addition, our estimation of the shear viscosity for the TIP4P/2005 at 273 K in NVT and NVE ensembles with 2000 water systems was η = 1.58 ± 0.10 and 1.58 ± 0.13 mPa s, respectively, which is in good agreement with the value of η = 1.65 ± 0.05 mPa s in the NVE simulation of 256 TIP4P/2005 water molecules at 273 K reported by Fanourgakins and Prosmiti.³³ Hereafter, we present the results obtained from the NVT simulations of 2000 water molecules unless otherwise stated.

Figure 2 shows the shear viscosities of the OPC, OPC3, TIP4P/2005, TIP4P-FB, and TIP3P-FB water models calculated at various temperatures. The figure also includes the calculated shear viscosities reported in the literature for other commonly used water models, namely, TIP3P and SPC/E, for comparison. All calculated shear viscosities for the OPC and OPC3 models, and the TIP4P/2005, TIP4P-FB, and TIP3P-FB models, including the statistical errors, are listed in Table II. The results for the TIP4P/2005, TIP4P-FB, and TIP3P-FB models were in good agreement with those reported in Refs. 14, 38, and 39. The OPC and OPC3 models yielded similar values over the examined temperature range. The difference in the cutoff length had negligible effects on the shear viscosity for these models, although the surface tension was strongly dependent on the cutoff length for the LJ interactions (Fig. S3 in the supplementary material and Ref. 18). For the OPC and OPC3 models, the calculated shear viscosities at higher temperatures corresponded well with the experimental results. However, a deviation from the experimental values was noted at temperatures below 313 K. At 298 K, the OPC and OPC3 models predict a η of ∼0.8 mPa s, 10% lower than the experimental value of 0.89 mPa s. At 273 K, the predicted η is 20% lower than the experimental value. Compared to other water models, the performance of the OPC and OPC3 models in terms of predicting the shear viscosity was lower than that of TIP4P/2005, TIP4P-FB, and TIP3P-FB, but only for temperatures below 293 K. However, it was notably better than the performance of the SPC/E and TIP3P models for the entire temperature range examined in this study.

The reason for the underestimation of the shear viscosities of the OPC and OPC3 water models at low temperatures was investigated by calculating the water oxygen–oxygen radial distribution function, g(r), at 273 K for the water models examined in this study, and these results are shown in Fig. 3. The higher g(r) value indicates a higher probability of finding a water molecule at the distance r. The height of the first peak in the g(r) decreased in the order TIP4P-FB > TIP4P/2005> TIP3P-FB > OPC3 > OPC. This order almost corresponds to the decreasing order of the measured shear viscosity at 273 K, TIP4P-FB > TIP3P-FB > TIP4P/2005 > OPC3 ≈ OPC (see Fig. 2), indicating that the shear viscosity is correlated with the height of the first peak of g(r). The first peak for the TIP3P-FB appeared at a slightly shorter distance than that of the TIP4P/2005. For the OPC water model, the first peak was positioned at a slightly longer distance than that of the OPC3. This result suggests that the shear viscosity would decrease with the first peak position in the g(r).

In this work, the shear viscosities of the OPC and OPC3 water models at various temperatures ranging from 273 to 373 K were calculated to assess their performance. Based on the results presented here, along with those from previous reports by other researchers, we can conclude that the OPC and OPC3 are among the best nonpolarizable water models, accounting for the various static and dynamic properties of water.

The densities, self-diffusion coefficients, surface tensions, and static dielectric constants of the OPC, OPC3, TIP4P/2005, TIP4P-FB, and TIP3P-FB water models for temperatures ranging from 273 to 373 K are provided as the supplementary material.

The author has no conflicts to disclose.

Tadashi Ando: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.