Water is a unique and abundant substance in biological and chemical systems. Considering its importance and ubiquity, numerous water models have been developed to reproduce various properties of bulk water in molecular simulations. Therefore, selecting an appropriate water model suitable for the properties of interest is crucial for computational studies of water systems. The four-point Optimal Point Charge (OPC) and three-point OPC (OPC3) water models were developed in 2014 and 2016, respectively. These models reproduce numerous properties of bulk water with high accuracy, such as density, dielectric constant, heat of vaporization, self-diffusion coefficient, and surface tension. In this study, we evaluated the shear viscosities of the OPC and OPC3 water models at various temperatures ranging from 273 to 373 K using the Green–Kubo formalism to assess their performance. The evaluated viscosities of both models were very close to each other at all the examined temperatures. At temperatures above 310 K, the calculated shear viscosities were in excellent agreement with the experimental results. However, at lower temperatures, the water models systematically underestimated the shear viscosity, with the calculated values at 273 and 298 K being 20% and 10% lower than the experimental values, respectively. Despite this limitation, the OPC and OPC3 water models outperformed other widely used water models.

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.3 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.2 

Currently, the Extended Simple Point Charge (SPC/E)4 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,6 TIP4P/2005, is widely used, which reproduces various experimental properties.7 In particular, for biomolecular simulations, the TIP3P water model has been extensively used within the Amber,8 OPLS-AA,9 and CHARMM10 force fields for an extended period,9 wherein CHARMM uses a slightly modified version of the TIP3P water model.11 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.12 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.3 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.13 Notably, their developmental approach refrained from imposing any geometric constraints on the model, except for symmetry, contrasting with the conventional approach.13 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.12 They suggested that the accuracy limit of the three-point rigid nonpolarizable water model had been reached.12 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.14 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.8 Tempra et al. also reported that the OPC water model accurately reproduces experimental surface tensions and is essential for realistic lipid monolayer simulations.18 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.19 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.

In the Green–Kubo approach, the shear viscosity, η, is related to the autocorrelation function (ACF) of the second-order pressure tensor P,20,
η=VkBT0Pst0+tPst0dt,
(1)
where V is the system volume, kB is the Boltzmann constant, T is the temperature, t is the time, the angle bracket indicates the ensemble average, and Ps is one of the shear components of the pressure tensor, namely, Pxy, Pyz, Pxz, PxxPyy/2, PyyPzz/2, and PxxPzz/2, for an isotropic system.21 The latter three components are obtained by 45° rotations of the pressure tensor about the z-, x-, and y-axes.21 

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 Amber1822 and Amber22,23 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,24 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,25 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 thermostat26 with a collision frequency of 1 ps−1 and the pressure was kept at 1 bar using a Berendsen barostat27 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 method27 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.28 Energies and pressure were saved every 0.01 ps. The Amber MD program computes only the three diagonal components of the pressure tensor, Pxx, Pyy, and Pzz, not the three off-diagonal components during NPT simulations using a Berendsen barostat. To compute the Pxx, Pyy, and Pzz 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 PxxPyy/2, PyyPzz/2, and PxxPzz/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 cpptraj29 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.

FIG. 1.

Shear viscosity, η(t), as a function of the integration time up to 100 (A) and 10 ps (B) with Eq. (1) for the OPC water model at four different temperatures. The values were obtained from 50 independent simulations of 2000 water molecules with the cutoff length of 8 Å for the Lennard-Jones interactions in the NVT ensemble. The error bars resulting from the calculation of η(t) are also shown.

FIG. 1.

Shear viscosity, η(t), as a function of the integration time up to 100 (A) and 10 ps (B) with Eq. (1) for the OPC water model at four different temperatures. The values were obtained from 50 independent simulations of 2000 water molecules with the cutoff length of 8 Å for the Lennard-Jones interactions in the NVT ensemble. The error bars resulting from the calculation of η(t) are also shown.

Close modal

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

TABLE I.

Shear viscosities, η, for the OPC and TIP4P/2005 water models calculated with 2000 and 256 water molecules under the NVT and NVE ensembles at 273 and 283 K.a

Water modelEnsembleNumber of watersTemperature (K)η (mPa s)
OPC NVT 2000 273 1.41 ± 0.08 
OPC NVT 256 273 1.44 ± 0.09 
OPC NVE 2000 273b 1.49 ± 0.11 
OPC NVT 2000 283 1.10 ± 0.05 
OPC NVT 256 283 1.11 ± 0.06 
OPC NVE 2000 283c 1.11 ± 0.08 
TIP4P/2005 NVT 2000 273 1.58 ± 0.10 
TIP4P/2005 NVE 2000 273d 1.58 ± 0.13 
Water modelEnsembleNumber of watersTemperature (K)η (mPa s)
OPC NVT 2000 273 1.41 ± 0.08 
OPC NVT 256 273 1.44 ± 0.09 
OPC NVE 2000 273b 1.49 ± 0.11 
OPC NVT 2000 283 1.10 ± 0.05 
OPC NVT 256 283 1.11 ± 0.06 
OPC NVE 2000 283c 1.11 ± 0.08 
TIP4P/2005 NVT 2000 273 1.58 ± 0.10 
TIP4P/2005 NVE 2000 273d 1.58 ± 0.13 
a

All simulations were performed with a cutoff length of 8 Å for the Lennard-Jones interactions. Prior to the NVT and NVE simulations, the simulation systems were equilibrated under NPT conditions at the desired target temperature and pressure of 1 bar.

b

Temperature averaged over the 10 ns NVE run was 274.0 ± 1.9 K.

c

Temperature averaged over the 10 ns NVE run was 284.1 ± 2.3 K.

d

Temperature averaged over the 10 ns NVE run was 274.2 ± 1.9 K.

The effects of temperature control algorithms on the viscosity calculations were investigated by Fanourgakins, Medina, and Prosmiti,33 and Hafner et al.34 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,27 stochastic rescaling,35 and Nosé–Hoover36 thermostats) with a wide range of thermostat coupling strengths.33,34 In contrast, local velocity randomizing thermostats (e.g., the Andersen37 and Langevin26 thermostats) with strong coupling drastically increase the viscosities compared to the NVE results.34 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.33 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.

FIG. 2.

Temperature dependence of the shear viscosities, η, of various water models. The results for the OPC, OPC3, TIP4P/2005, TIP4P-FB, and TIP3P-FB water models were obtained in this work, whereas those for the TIP3P and SPC/E were taken from Refs. 38 and 33, respectively. The experimental values are from Ref. 40. The inset shows an enlargement of the region surrounding 273 K.

FIG. 2.

Temperature dependence of the shear viscosities, η, of various water models. The results for the OPC, OPC3, TIP4P/2005, TIP4P-FB, and TIP3P-FB water models were obtained in this work, whereas those for the TIP3P and SPC/E were taken from Refs. 38 and 33, respectively. The experimental values are from Ref. 40. The inset shows an enlargement of the region surrounding 273 K.

Close modal
TABLE II.

Shear viscosities (mPa s), η, for the OPC, OPC3, TIP4P/2005, TIP4P-FB, and TIP3P-FB water models in the temperature range 273–373 K.a

Temperature (K)OPCOPC3TIP4P/2005TIP4P-FBTIP3P-FB
8 Å cutoff12 Å cutoff8 Å cutoff12 Å cutoff8 Å cutoff8 Å cutoff8 Å cutoffExperimentb
273 1.41 ± 0.07 1.42 ± 0.07 1.45 ± 0.07 1.43 ± 0.07 1.58 ± 0.10 1.71 ± 0.09 1.63 ± 0.09 1.79 
283 1.11 ± 0.05 1.09 ± 0.06 1.11 ± 0.05 1.10 ± 0.06 1.18 ± 0.05 1.28 ± 0.07 1.25 ± 0.06 1.31 
293 0.88 ± 0.04 0.88 ± 0.04 0.88 ± 0.04 0.89 ± 0.05 0.93 ± 0.04 0.99 ± 0.05 0.96 ± 0.06 1.01 
298 0.79 ± 0.04 0.80 ± 0.04 0.81 ± 0.04 0.79 ± 0.04 0.82 ± 0.04 0.89 ± 0.04 0.88 ± 0.04 0.90 
303 0.73 ± 0.03 0.72 ± 0.03 0.72 ± 0.04 0.73 ± 0.04 0.75 ± 0.03 0.80 ± 0.04 0.79 ± 0.04 0.80 
313 0.62 ± 0.03 0.62 ± 0.03 0.61 ± 0.03 0.61 ± 0.03 0.62 ± 0.03 0.66 ± 0.04 0.66 ± 0.03 0.66 
333 0.46 ± 0.03 0.46 ± 0.02 0.45 ± 0.03 0.46 ± 0.03 0.45 ± 0.03 0.48 ± 0.03 0.49 ± 0.02 0.47 
353 0.37 ± 0.02 0.37 ± 0.02 0.36 ± 0.02 0.36 ± 0.01 0.35 ± 0.02 0.37 ± 0.02 0.38 ± 0.02 0.36 
373 0.30 ± 0.02 0.30 ± 0.01 0.29 ± 0.01 0.29 ± 0.02 0.28 ± 0.01 0.30 ± 0.01 0.30 ± 0.01 0.28 
Temperature (K)OPCOPC3TIP4P/2005TIP4P-FBTIP3P-FB
8 Å cutoff12 Å cutoff8 Å cutoff12 Å cutoff8 Å cutoff8 Å cutoff8 Å cutoffExperimentb
273 1.41 ± 0.07 1.42 ± 0.07 1.45 ± 0.07 1.43 ± 0.07 1.58 ± 0.10 1.71 ± 0.09 1.63 ± 0.09 1.79 
283 1.11 ± 0.05 1.09 ± 0.06 1.11 ± 0.05 1.10 ± 0.06 1.18 ± 0.05 1.28 ± 0.07 1.25 ± 0.06 1.31 
293 0.88 ± 0.04 0.88 ± 0.04 0.88 ± 0.04 0.89 ± 0.05 0.93 ± 0.04 0.99 ± 0.05 0.96 ± 0.06 1.01 
298 0.79 ± 0.04 0.80 ± 0.04 0.81 ± 0.04 0.79 ± 0.04 0.82 ± 0.04 0.89 ± 0.04 0.88 ± 0.04 0.90 
303 0.73 ± 0.03 0.72 ± 0.03 0.72 ± 0.04 0.73 ± 0.04 0.75 ± 0.03 0.80 ± 0.04 0.79 ± 0.04 0.80 
313 0.62 ± 0.03 0.62 ± 0.03 0.61 ± 0.03 0.61 ± 0.03 0.62 ± 0.03 0.66 ± 0.04 0.66 ± 0.03 0.66 
333 0.46 ± 0.03 0.46 ± 0.02 0.45 ± 0.03 0.46 ± 0.03 0.45 ± 0.03 0.48 ± 0.03 0.49 ± 0.02 0.47 
353 0.37 ± 0.02 0.37 ± 0.02 0.36 ± 0.02 0.36 ± 0.01 0.35 ± 0.02 0.37 ± 0.02 0.38 ± 0.02 0.36 
373 0.30 ± 0.02 0.30 ± 0.01 0.29 ± 0.01 0.29 ± 0.02 0.28 ± 0.01 0.30 ± 0.01 0.30 ± 0.01 0.28 
a

The shear viscosities were estimated from 50 independent simulations under NVT conditions. Prior to the NVT simulations, the simulation systems were equilibrated under NPT conditions at the desired target temperature and pressure of 1 bar.

b

Experimental values are from Ref. 40.

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).

FIG. 3.

Oxygen–oxygen radial distribution function, g(r), for various water models calculated from the NVT simulations at 273 K under 1 bar. All simulations were performed using 2000 water molecules in cubic boxes with a cutoff length of 8 Å for Lennard-Jones interactions. The bin size is 0.1 Å. The curves are natural cubic spline of the data points. (A) g(r) for oxygen–oxygen distances less than 8 Å, and (B) magnification of the region containing the first peak of the g(r).

FIG. 3.

Oxygen–oxygen radial distribution function, g(r), for various water models calculated from the NVT simulations at 273 K under 1 bar. All simulations were performed using 2000 water molecules in cubic boxes with a cutoff length of 8 Å for Lennard-Jones interactions. The bin size is 0.1 Å. The curves are natural cubic spline of the data points. (A) g(r) for oxygen–oxygen distances less than 8 Å, and (B) magnification of the region containing the first peak of the g(r).

Close modal

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.

1.
B.
Guillot
,
J. Mol. Liq.
101
,
219
(
2002
).
2.
A. V.
Onufriev
and
S.
Izadi
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1347
(
2018
).
3.
S. P.
Kadaoluwa Pathirannahalage
,
N.
Meftahi
,
A.
Elbourne
,
A. C. G.
Weiss
,
C. F.
McConville
,
A.
Padua
,
D. A.
Winkler
,
M.
Costa Gomes
,
T. L.
Greaves
,
T. C.
Le
,
Q. A.
Besford
, and
A. J.
Christofferson
,
J. Chem. Inf. Model.
61
,
4521
(
2021
).
4.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
5.
W. L.
Jorgensen
,
J. Am. Chem. Soc.
103
,
335
(
1981
).
6.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
,
J. Chem. Phys.
79
,
926
(
1983
).
7.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
8.
C.
Tian
,
K.
Kasavajhala
,
K. A. A.
Belfon
,
L.
Raguette
,
H.
Huang
,
A. N.
Migues
,
J.
Bickel
,
Y. Z.
Wang
,
J.
Pincay
,
Q.
Wu
, and
C.
Simmerling
,
J. Chem. Theory Comput.
16
,
528
(
2020
).
9.
W. L.
Jorgensen
and
J.
Tirado-Rives
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
6665
(
2005
).
10.
R. B.
Best
,
X.
Zhu
,
J.
Shim
,
P. E. M.
Lopes
,
J.
Mittal
,
M.
Feig
, and
A. D.
MacKerell
,
J. Chem. Theory Comput.
8
,
3257
(
2012
).
11.
W. H.
Reiher
III
, Ph.D. thesis,
Harvard University
,
1985
.
12.
S.
Izadi
and
A. V.
Onufriev
,
J. Chem. Phys.
145
,
074501
(
2016
).
13.
S.
Izadi
,
R.
Anandakrishnan
, and
A. V.
Onufriev
,
J. Phys. Chem. Lett.
5
,
3863
(
2014
).
14.
L. P.
Wang
,
T. J.
Martinez
, and
V. S.
Pande
,
J. Phys. Chem. Lett.
5
,
1885
(
2014
).
15.
Z.
Li
,
L. F.
Song
,
P. F.
Li
, and
K. M.
Merz
,
J. Chem. Theory Comput.
16
,
4429
(
2020
).
16.
Z.
Li
,
L. F.
Song
,
P. F.
Li
, and
K. M.
Merz
,
J. Chem. Theory Comput.
17
,
2342
(
2021
).
17.
A.
Sengupta
,
Z.
Li
,
L. F.
Song
,
P. F.
Li
, and
K. M.
Merz
,
J. Chem. Inf. Model.
61
,
869
(
2021
).
18.
C.
Tempra
,
O. H. S.
Ollila
, and
M.
Javanainen
,
J. Chem. Theory Comput.
18
,
1862
(
2022
).
19.
I. C.
Yeh
and
G.
Hummer
,
J. Phys. Chem. B
108
,
15873
(
2004
).
20.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
2017
).
21.
B. L.
Holian
and
D. J.
Evans
,
J. Chem. Phys.
78
,
5147
(
1983
).
22.
D. A.
Case
,
I. Y.
Ben-Shalom
,
S. R.
Brozell
,
D. S.
Cerutti
,
T. E.
Cheatham
III
,
V. W. D.
Cruzeiro
,
T. A.
Darden
,
R. E.
Duke
,
D.
Ghoreishi
,
G.
Giambasu
,
T.
Giese
,
M. K.
Gilson
,
H.
Gohlke
,
A. W.
Goetz
,
D.
Greene
,
R.
Harris
,
N.
Homeyer
,
Y.
Huang
,
S.
Izadi
,
A.
Kovalenko
,
R.
Krasny
,
T.
Kurtzman
,
T. S.
Lee
,
S.
LeGrand
,
P.
Li
,
C.
Lin
,
J.
Liu
,
T.
Luchko
,
R.
Luo
,
V.
Man
,
D. J.
Mermelstein
,
K. M.
Merz
,
Y.
Miao
,
G.
Monard
,
C.
Nguyen
,
H.
Nguyen
,
A.
Onufriev
,
F.
Pan
,
R.
Qi
,
D. R.
Roe
,
A.
Roitberg
,
C.
Sagui
,
S.
Schott-Verdugo
,
J.
Shen
,
C. L.
Simmerling
,
J.
Smith
,
J.
Swails
,
R. C.
Walker
,
J.
Wang
,
H.
Wei
,
L.
Wilson
,
R. M.
Wolf
,
X.
Wu
,
L.
Xiao
,
Y.
Xiong
,
D. M.
York
, and
P. A.
Kollman
, “
AMBER 2019
” (
University of California
,
San Francisco
,
2019
).
23.
D. A.
Case
,
H. M.
Aktulga
,
K.
Belfon
,
J. T.
Ben-Shalom
,
S. R.
Berryman
,
S. R.
Brozell
,
D. S.
Cerutti
,
T. E.
Cheatham
III
,
G. A.
Cisneros
,
V. W. D.
Cruzeiro
,
T. A.
Darden
,
R. E.
Duke
,
G.
Giambasu
,
M. K.
Gilson
,
H.
Gohlke
,
A. W.
Goetz
,
R.
Harris
,
S.
Izadi
,
S. A.
Izmailov
,
K.
Kasavajhala
,
M. C.
Kaymak
,
E.
King
,
A.
Kovalenko
,
T.
Kurtzman
,
T. S.
Lee
,
S.
LeGrand
,
P.
Li
,
C.
Lin
,
J.
Liu
,
T.
Luchko
,
R.
Luo
,
M.
Machado
,
V.
Man
,
M.
Manathunga
,
K. M.
Merz
,
Y.
Miao
,
O.
Mikhailovskii
,
G.
Monard
,
H.
Nguyen
,
K. A.
O’Hearn
,
A.
Onufriev
,
F.
Pan
,
S.
Pantano
,
R.
Qi
,
A.
Rahnamoun
,
D. R.
Roe
,
A.
Roitberg
,
C.
Sagui
,
S.
Schott-Verdugo
,
A.
Shajan
,
J.
Shen
,
C. L.
Simmerling
,
N. R.
Skrynnikov
,
J.
Smith
,
J.
Swails
,
R. C.
Walker
,
J.
Wang
,
J.
Wang
,
H.
Wei
,
R. M.
Wolf
,
X.
Wu
,
Y.
Xiong
,
Y.
Xue
,
D. M.
York
,
S.
Zhao
, and
P. A.
Kollman
, “
AMBER 2022
” (
University of California
,
San Francisco
,
2022
).
24.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
,
10089
(
1993
).
25.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J. C.
Berendsen
,
J. Comput. Phys.
23
,
327
(
1977
).
26.
R. J.
Loncharich
,
B. R.
Brooks
, and
R. W.
Pastor
,
Biopolymers
32
,
523
(
1992
).
27.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
Vangunsteren
,
A.
Dinola
, and
J. R.
Haak
,
J. Chem. Phys.
81
,
3684
(
1984
).
28.
J. E.
Basconi
and
M. R.
Shirts
,
J. Chem. Theory Comput.
9
,
2887
(
2013
).
29.
D. R.
Roe
and
T. E.
Cheatham
,
J. Chem. Theory Comput.
9
,
3084
(
2013
).
30.
P. J.
Daivis
and
D. J.
Evans
,
J. Chem. Phys.
103
,
4261
(
1995
).
31.
O. A.
Moultos
,
Y.
Zhang
,
I. N.
Tsimpanogiannis
,
I. G.
Economou
, and
E. J.
Maginn
,
J. Chem. Phys.
145
,
074109
(
2016
).
32.
K. S.
Kim
,
M. H.
Han
,
C.
Kim
,
Z.
Li
,
G. E.
Karniadakis
, and
E. K.
Lee
,
J. Chem. Phys.
149
,
044510
(
2018
).
33.
G. S.
Fanourgakis
,
J. S.
Medina
, and
R.
Prosmiti
,
J. Phys. Chem. A
116
,
2564
(
2012
).
34.
R.
Hafner
,
G.
Guevara-Carrion
,
J.
Vrabec
, and
P.
Klein
,
J. Phys. Chem. B
126
,
10172
(
2022
).
35.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
36.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
37.
H. C.
Andersen
,
J. Chem. Phys.
72
,
2384
(
1980
).
38.
M. A.
González
and
J. L.
Abascal
,
J. Chem. Phys.
132
,
096101
(
2010
).
39.
Y. D.
Qiu
,
P. S.
Nerenberg
,
T.
Head-Gordon
, and
L. P.
Wang
,
J. Phys. Chem. B
123
,
7061
(
2019
).
40.
E. W.
Lemmon
,
I. H.
Bell
,
M. L.
Huber
, and
M. O.
McLinden
,
NIST Chemistry WebBook
,
NIST Standard Reference Database Number 69
, edited by
P. J.
Linstrom
and
W. G.
Mallard
(
National Institute of Standards and Technology
,
Gaithersburg, MD
,
2023
).
Published open access through an agreement with Tokyo Rika Daigaku

Supplementary Material