A comparative computational study of coarse-grained and all-atom water models in shock Hugoniot states

Molecular dynamics (MD) computer simulation technique is widely used to investigate the influence of shock waves on properties of different materials.^1–12 To study how shock waves interact with biological entities like cell membranes or biological fluids, we often start with simulations that employ simplified models of fluids and membranes.^13,14 Although biological membranes are a complex mixture of lipid molecules and proteins, in simulations that study mechanical properties of membranes, they are often considered to be lipid bilayers containing just one component lipid or a mixture of two or three lipid components. Membranes (and other biological entities) are solvated in fluids that mostly contain water; consequently bio-fluids in simulations are often considered to be just water. Therefore, it is important to understand how well water models used in simulations perform when water is exposed to shock waves. There exist a large number of all-atom and also coarse-grained (CG) water models with different levels of sophistication, which describe water for simulation purposes.¹⁵ While all-atom force fields provide a detailed picture of structural changes in water exposed to extreme conditions created by shock wave passage, simulations of biological systems often require use of CG models due to the large size of the systems and need to perform simulations over long time periods.

In order to judge the performance of the force field in the simulations where the shock wave propagates across the material, one often considers the Hugoniot plot that displays how the velocity of shock propagation (u_s) changes as a function of the velocity of a body (piston), u_p, which moves in the material and creates shock waves.^16,17 Such a Hugoniot plot was already obtained for two all-atom force fields describing water. In one case, it was for simulations where the reactive force field (ReaxFF) that can describe bond-breaking and bond-making was used.¹³ It was observed that no bond-breaking or bond-making takes place when u_p has values below 3 km/s, so regular rigid body models of water can be used in these cases. Therefore, in another simulation that studied the influence of cavitation effect on the integrity of lipid bilayers, the rigid SPC model of water was used.¹⁷ The all-atom models of water mentioned above successfully reproduced the Hugoniot plot data obtained from the experiment. For CG models, the Hugoniot plots were constructed only for MARTINI type models. It was observed¹⁴ that Hugoniot plot data obtained for standard MARTINI CG water were in nice agreement with experiment at T₀ = 323 K, when the value of u_p and therefore of u_s are not very large, i.e., when the shock waves are not very intense (u_p has values up to 1.0 km/s). More recently, Adhikari et al.¹⁸ reported that a polarizable MARTINI (PolM) model of water¹⁹ produced a nice fit to an experimental Hugoniot plot (u_s vs. u_p) even for large intensity shock waves. PolM water was developed to accomplish better description of CG simulations using MARTINI for systems when electrical properties of water play an important role. In the standard MARTINI water model (we shall call it non-polarizable MARTINI water, NpolM), four water molecules are replaced by one Lennard-Jones particle. In the PolM water model, the four water molecules from the all-atom model are replaced by three particles. The central particle, W, is neutral and interacts with other W particles through the Lennard-Jones interactions. The other two particles are charged, one positively (WP) and one negatively (WM), so the PolM water is neutral but has a dipole moment. The distances WP–W and WM–W are fixed in PolM, but to allow for dipole fluctuations, a harmonic angular potential for the angle WP–W–WM is also added to the force field. The charged WP and WM particles on different water particles interact just through Coulomb interactions. More details on the PolM model of water and its properties can be found in references.^19,20

In this paper, we focused our attention on the performance of certain water models when influenced by the shock waves. Therefore, we performed MD simulations on shocked water using the standard NpolM CG water model, the PolM CG water model, and also one of the best performing all-atom models, the TIP4P/2005 model.²¹ Different values of u_p were applied to generate the shock waves in Hugoniot shocked states. For each of these states, thermodynamic properties were monitored. Since Hugoniot states are compressed states at high pressures and low specific volumes, our simulations also provide information on the behavior of the water models when they are subjected to these extreme conditions.

We performed MD simulations on three types of systems containing water under shock wave load, and in each type of system, we described water using a different model. Two of the water models we used were coarse-grained, NpolM and PolM, and one of the models was all-atom TIP4P/2005 model. Initially, water molecules (in the case of TIP4P/2005 potential) or particles (in the case of CG potentials) were placed in a rectangular simulation box of size ∼10 × 10 × 100 nm³. These initial systems were equilibrated by performing MD runs for 1 ns using the isothermal-isobaric (NPT) ensemble at 1 bar and 300 K.

For the shock wave generation, we used the momentum mirror method²² and its main idea is illustrated in Fig. 1. According to this method, rigid walls that reflected incoming particles, i.e., momentum mirrors, were placed at both ends of the simulation box at the z-axis. Vacuum layers of 2 and 200 nm thickness were added in the front of the left wall (located at z = 0 nm) and at the right wall, respectively. Particle velocity (u_p) directed toward the left wall at z = 0 nm was added to all particles creating a flow toward this wall. For every system we studied, u_p was considered to have values of 0.5, 1.0, 1.5, and 2.0 km/s. The procedure we employed is equivalent to having an infinitely massive piston pushing water particles toward the right wall and thus creating a shock wave in the system traveling in the +z direction.

The simulations with propagating shock wave were performed using the microcanonical (NVE) ensemble, and periodic boundary condition (PBC) in the xy plane was applied. Electrostatic interactions were calculated by the shifted cutoff and particle mesh Ewald (PME) method²³ for simulations with PolM and TIP4P/2005 water, respectively. The LINCS algorithm²⁴ was applied to satisfy the bond length constraint of the water models, and the time steps of 5 and 1 fs were used for the two CG MARTINI waters and TIP4P/2005 water, respectively.

Since the goal of CG simulations is to speed-up the calculations, we show in Table I that using the coarse-graining potentials NpolM and PolM for water allows us to achieve a speed-up factor of about 148 and 45, respectively, compared with the all-atom TIP4P/2005 potential.

To calculate thermodynamic properties of our three water models, additional MD simulations were performed on bulk water systems corresponding to the presence of 20 000 water molecules in a cubic box. The thermal expansion coefficient, α, isothermal bulk modulus, B, and heat capacity, C_V, were measured by a finite difference approximation

The α and B values were calculated by performing 10-ns NPT simulations, and the C_V values were calculated by performing 10-ns canonical (NVT) simulations at T and P (T and V) values corresponding to Hugoniot states. The values of ΔT and ΔP were set to be 5% of the corresponding Hugoniot temperature and pressure.

All simulations were performed using GROMACS 4.5.7 and GROMACS 2016.3 packages.²⁵ The trajectories and local pressure from the shock wave simulations were analyzed using the MDAnalysis library²⁶ and the customized GROMACS 4.5.5-LS package,^27–29 respectively.

It was observed previously that application of shock waves can ionize water when the pressure reaches values of 15–20 GPa.^30,31 Since we examined here two CG and one all-atom models that are not suitable to describe reactive events, we studied the behavior of our water models in the u_p range of 0.5–2.0 km/s, for which the produced pressure is below the ionization pressure.

Figures 2(a)–2(d) display water density maps along the z-axis after the shock wave propagated across our systems containing water described by our three models. It can be seen from these figures that passage of shock waves compresses water creating a distinct discontinuity in the density and that each compressed region has a homogeneous density. Such discontinuous jumps in density are characteristic of the shock wave. Obviously, the velocity of density discontinuity (u_d) depends on the value of u_p and on the water model. Figure 3 shows Hugoniot plots that characterize the passage of shock waves. Hugoniot plots can be presented either as a plot depicting the dependence of shock velocity u_s on the “piston” velocity u_p, like it is done in Fig. 3(a) or as a plot depicting pressure of Hugoniot states versus density of these states [see Fig. 3(b)].¹⁰ 

The shock velocity (u_s) in Fig. 3(a) was calculated from the sum of the u_d, the velocity of density discontinuity (measured in the time interval from t = 5 to t = 15 ps) and the value of u_p, using the following equation:¹⁴ 

The shock velocity can be also obtained from the Rankine–Hugoniot equation as follows:

where v and v₀ are the specific volumes of the compressed and uncompressed regions, respectively. As we can see from Fig. 3(a), the values of u_s obtained either by Eq. (4) or by Eq. (5) are almost the same, meaning that the momentum mirror method we use for shock wave creation reproduces well Hugoniot predictions made by using continuum arguments. In the case of weak shock (u_p = 0.5 km/s), the u_s versus u_p Hugoniot data for both NpolM and PolM water are in good agreement with that of the all-atom TIP4P/2005 water, and values for all three models display good agreement with the value for the Hugoniot state obtained from the equations recently reported by Myint et al.¹⁶ that describe shock compressed water. However, at u_p ≥ 1.0 km/s, the values of u_s for the NpolM water deviate from the reference values so that at u_p = 2.0 km/s the value of u_s is ∼12% larger than the reference one. It is interesting that for Hugoniot plot depicting u_s versus u_p, the agreement between the results obtained from simulations with PolM water is even slightly better than results obtained by using the all-atom model.

The u_s–u_p Hugoniot plot shows that one can establish a linear relation between u_s and u_p which has a form

where C₀ is the speed of sound in the bulk medium and s is the slope parameter. We can see from Table II that the C₀ values for all three water models are close to the speed of sound in water (1.50 km/s) at 300 K. Because Eq. (5) implies that the initial and Hugoniot specific volumes affect the value of the u_s, the fit to Eq. (6) produces different slope s values for our three water models that have different initial specific volume.

As was mentioned above, the u_s–u_p Hugoniot plot shown in Fig. 3(a) can be transformed into a plot depicting shocked compressed states at pressure, P as a function of a specific volume, v. We calculated the pressure in the compressed region using two different methods: one uses a calculation of the local stress tensor by using the customized GROMACS-LS 4.5.5 package and the other method uses the Rayleigh line equation obtained by using continuum macroscopic conservation laws,

where P and P₀ are the pressure of the compressed and uncompressed regions, respectively. As shown in Fig. 3(b), the pressures calculated from local stress tensor are also in very good agreement with those obtained by using the continuum approach. Figure 3(b) again shows that results from simulations with the NpolM model do not match well the Hugoniot curve when u_p is larger than 0.5 km/s. It is noteworthy to mention here that the PolM model and especially TIP4P/2005 model reproduce the reference Hugoniot curve reasonably well for values of u_p that are less than or equal to 2.0 km/s.

The PolM model does a better job than the NpolM model in describing the Hugoniot states of water. Is the structure of compressed PolM water very different from that of the NpolM model? Figure 4 shows the W–W radial distribution functions (RDFs) of the two CG models in the compressed region during the shock simulation at t = 10 ps. The RDF plots reveal the structural similarity between the two MARTINI waters at each Hugoniot state, nevertheless the peaks of the RDF from PolM water are slightly taller and the first and second hydration shells of PolM water are shifted to shorter distances. These factors point to an increased water ordering in PolM water and larger coordination numbers when compared to NpolM water (data not shown), thus producing more compressed states for the PolM model, as also could be seen in Fig. 3(b).

According to the Hugoniot equations, the internal energy gain per unit mass, as a result of shock compression is

This equation is obtained by using the laws of conservation of mass, momentum, and energy, and it reduces to simple relationship ΔU ≈ 1/2 u_p² since the Hugoniot pressure P$>>$P₀. Using the MD program, we calculated the total energy gain per unit mass for the NpolM, PolM, and TIP4P/2005 water models at different values of u_p. Indeed, Fig. 5(a) indicates that total energy gain per unit mass, irrespective of the model, follows the value of 1/2 u_p², meaning that the kinetic energy added by the assignment of u_p at the initial stage is converted into the total energy gain during the NVE MD runs. Nevertheless, as Fig. 5(b) shows, the kinetic energy gain varies among the water models, due to the difference in the potential energy functions.

Using the values of the kinetic energy, we calculated the temperature in the compressed region at different values of u_p [see Fig. 5(c)]. The temperature curve we obtained for the TIP4P/2005 water was in line with that predicted by the theory of Gurtman et al.,^32,33 while the curves for CG MARTINI waters displayed much higher temperatures compared to the one from the all-atom model. In particular, the temperature of the NpolM water was very high and it exceeded 5000 K at u_p = 2.0 km/s. This extremely high temperature has also been reported in other MD simulations with CG force fields,^7,34 and it is due to the reduction in the number of degrees of freedom. Considering that there was a small difference in the kinetic energy gain between the NpolM and PolM waters at each compressed state [see Fig. 5(b)], the relationship between temperatures of two MARTINI waters can be explained by using the following expression: 3ΔT_NpolM ≈ 7ΔT_PolM. This expression resulted from the observation that the kinetic part of heat-capacity of NpolM water is just 3/2 R, while it is 7/2 R for PolM water. The higher temperature of PolM water, when compared to temperature in all-atom simulations, is also due to a smaller specific heat-capacity of the CG model (we replaced four water molecules by three particles when coarse-graining in this model). It is possible that a low potential energy of PolM water also contributed to the difference in temperature.

Figure 5(d) represents the velocity map of the three water models at u_p = 2 km/s and at t = 10 ps obtained by averaging the velocity vectors in the yz plane using a grid size of 1 nm. The clear discontinuities in the velocity vectors are observed in the map and these are correlated with discontinuities in density. From Eqs. (4), (5), and (7), we get that

which means that u_d at the same u_p value is proportional to the product of two factors: the specific volume v and change in pressure P − P₀. As one can see from Fig. 3(b), there is a substantial decrease in the specific volume, determining the change in Eq. (9), as one goes from NpolM to TIP4P/2005 and finally to the PolM model. Correspondingly, the u_d and therefore u_s values follow the same trend. We can also observe from Fig. 5(d) that the magnitudes of the velocities in the compressed region reveal their Hugoniot temperature which has a highest value in the case of the NpolM water.

To get further insight into the performance of our water models under high pressure and density, we calculated the Grüneisen parameter, which is a useful quantity for description of the thermomechanical properties of liquids,³⁶ as well as of solid materials.^37,38 This parameter is defined as follows:

where α is the thermal expansion coefficient, B is the isothermal bulk modulus, and C_V is the heat capacity. We observed freezing of PolM water at u_p = 0.5, 1.0, and 1.5 km/s in the additional MD simulations although experimental studies do not confirm this.^16,39 Therefore, we did not calculate the Grüneisen parameters in these cases. We also did not study the thermodynamic properties of NpolM water at u_p = 0.5 km/s because we observed that in this case, the sample was very close to its freezing point and thermodynamic parameters given by Eqs. (1)–(3) could not be obtained reliably.

The values of the Grüneisen parameter as a function of the specific volume obtained from calculations of α, B, and specific C_V [Figs. 6(a)–6(c)] are shown in Fig. 6(d) for the three models of water. Two reference curves reported by Gurtman et al. and Myint et al. are also presented in Fig. 6(d) together with another reference point taken from the IAPWS-95 formulations³⁵ at 1 GPa and 340 K since the IAPWS-95 is valid only up to 1 GPa. These pressure and temperature are similar to the Hugoniot state of TIP4P/2005 at u_p = 0.5 km/s (1.25 GPa and 340 K). Indeed, the values of α, B, and C_V of TIP4P/2005 at u_p = 0.5 km/s are in good agreement with this reference point. (We need to emphasize here that all reference data are not experimental but are obtained from theory.) Figures 6(a)–6(c) show that the values of C_V for the CG models were much smaller than for the all-atom model, as expected. Also, the values of α were smaller for CG models, again because of a large change in temperature along Hugoniot curves. The values of B are not very sensitive to the model although CG models predict slightly larger compressibility. We can see from Fig. 6 that the NpolM water model strongly mis-estimates the thermomechanical properties of water at Hugoniot states, mostly due to its incorrect temperature behavior. It is interesting that despite that the total number of degrees of freedom was reduced for the PolM water model compared to the all-atom model, it is still successfully predicting thermomechanical properties of a compressed state when u_p = 2.0 km/s. Unfortunately, this model of water has a freezing problem in bulk simulations corresponding to compressed states we obtained when u_p = 0.5, 1.0, and 1.5 km/s.

Finally, Fig. 6(d) shows that the TIP4P/2005 model can reproduce Grüneisen parameters for liquid water at high pressures and high temperatures, as we see from comparison with the reference Grüneisen curves.

Understanding the behavior of water models in simulations performed under strong compression is of special interest when we study the effect of shock waves on aqueous systems including biological systems. Therefore, we carried out MD simulations on three different water models to study how these models perform when shock waves pass through water, creating extreme conditions of high pressure and density. Two of the three models were CG models used in the MARTINI force field: one regular non-polarizable water model (NpolM) and the other is a polarizable version of water in MARTINI (PolM). The third model was an all-atom model, TIP4P/2005, which is considered to be among the best all-atom models for water.

In our simulations, we observed that both PolM and TIP4P/2005 water models well reproduced the pressure and density of the reference shocked Hugoniot states, if the passage of shock waves occurs over short time periods. The NpolM model is fairly accurate only when the shock wave intensity is weak (less than 1.5 GPa). On the other hand, the two CG MARTINI waters have intrinsic limitations in reproducing the Hugoniot temperature and thermodynamic properties, resulting from the reduction in the number of degrees of freedom. The CG water models display extremely high temperatures, low thermal expansion coefficient, slightly high compressibility, and low heat capacity at their Hugoniot states, when compared with the all-atom water and the referenced data. We also observed the solidification of PolM water at elevated pressures, but less than 7 GPa, even with overpredicted temperatures, when simulations were run for a longer period of time. It indicates that equilibrium phase of PolM water at those Hugoniot states is incorrect. For higher pressure corresponding to u_p = 2.0 km/s, the temperature of the PolM water sample is high and this prevents its solidification. No solidification is observed in simulations that used TIP4P/2005 all-atoms model, and the thermodynamic data obtained from simulations with this model show good agreement with the reference data.

In summary, we conclude that MARTINI type CG water models should be used with caution in simulations that study the effects of shock waves although they can be used when the passage of shock waves occurs over very short time periods and the temperature of the compressed water region is not of interest. Judging by the behavior of the Grüneisen parameter and the Hugoniot curves, the TIP4P/2005 model seems to perform very well under shocked Hugoniot states (up to ∼11 GPa) when compared to reference data.

The support by a Grant No. N00014-17-1-2164 from the Office of Naval Research is gratefully acknowledged.