Different output stress values when using different TIP3P water models

Hi everyone,

This link - https://lammps.sandia.gov/doc/Howto_tip3p.html - shows that different water models (with some different parameters for O and H) can be implemented for TIP3P water (the original 1983 TIP3P model, TIP3P from Mackerell et al. in 1998, and TIP3P from Price et al. in 2004).

I have tried to run the same simulation but with those 3 different water models (the simulated system is a phospholipid bilayer under deformation - with 72 lipids and 9072 water molecules in the system). And it turns out that the output stress values is different when different water models are used (other parameters and settings are the same). Besides, I also tried to run the same simulation with different water models in GROMACS, but there is no significant difference in the stress values, and stress values in GROMACS for all cases are much lower compared to LAMMPS.

So my questions are:

  • Why and how some small differences in the water parameters would significantly affect the output stress values in LAMMP simulation?

  • Is there any explanation for the difference in stress values between LAMMPS and GROMACS (as all other settings are the same)? In other words, how would the water parameters affect the output stress value, leading to higher stress values in LAMMPS compared to GROMACS?

Please accept my apology if these are obvious questions. Thank you very much for your help.

Best regards,

Anh Vo

Hi everyone,

[…]

So my questions are:

  • Why and how some small differences in the water parameters would significantly affect the output stress values in LAMMP simulation?

dense water is a rather incompressible material, thus what may look as small changes in parameters can have significant changes in stress. you can reverse this and, for example, look at the relation of density to pressure. also, the question is what you consider a significant difference. typically people look at systems around atmospheric pressure (or “zero pressure”) where relative changes are particularly large.

  • Is there any explanation for the difference in stress values between LAMMPS and GROMACS (as all other settings are the same)? In other words, how would the water parameters affect the output stress value, leading to higher stress values in LAMMPS compared to GROMACS?

if there would be exactly the same physics model implemented, there must be the same results (the math does not change). however, what you may consider “the same” is likely not the exactly the same.

Axel.

Dear Dr. Axel,

Thank you for your quick response.

As I’m simulating a lipid bilayer under deformation, the von Mises stress-strain curve is plot for analysis. So the difference is the different stress values in the stress-strain curve when different water models are used. Specifically, In LAMMPS, the output stress values is the highest in simulation with Price water model, then decreases in simulation with the original 1983 TIP3P model, and the stress is lowest in simulation with Mackerell water model. I also run the same simulation in GROMACS, but the von Mises stress in stress-strain curve remain similar when different water models are used. (Please see the attachment for the plots).

In both LAMMPS and GROMACS, the system contains 72 POPC and 9072 water molecules. The two equilibrated bilayer structures in LAMMPS and GROMACS have similar number of atoms, as well as periodicity and membrane dimensions. The temperature is 310K and the pressure is 1 atm during deformation simulation, and other parameters in the input file are set to be the same.

So I don’t understand why von Mises stress is different in LAMMPS (but remain similar in GROMACS) when different water models are used?

I’m sorry for not being clear on what the output stress and what the difference is in my previous email. Thank you very much for your help.

Best regards,

Anh Vo

StressStrainCurve_differentWaterModels_LAMMPS_vs_GROMACS.tif (111 KB)

Dear Dr. Axel,

Thank you for your quick response.

[…]

So I don’t understand why von Mises stress is different in LAMMPS (but remain similar in GROMACS) when different water models are used?

while some of the settings you enter are the same, there are differences in a number of details between Gromacs and LAMMPS.

  • there are a variety of lennard-jones variants with shifted or truncated or smoothly truncated functional forms
  • you didn’t mention what kind of thermostat and time integrator you are using and if or not you are using a barostat
  • you didn’t say if you use tail correction for the LJ contribution to pressure
  • you didn’t say whether you are using a single or double precision version of Gromacs, similarly whether you are using anything similar in LAMMPS (e.g. from USER-INTEL or GPU).
  • Gromacs and LAMMPS have different ways of building their neighbor lists and decide differently what neighbors to include
  • You didn’t state which convergence criteria you are using for long-range electrostatics. Please also note, that those numbers cannot be directly compared between LAMMPS and Gromacs and they do use different methods (PPPM vs. SPME).
  • LAMMPS in general has many more options and settings and thus also more opportunities to overlook something
  • you didn’t say, if you are using the default settings for tabulation with bitmask lookup for long-range coulomb in LAMMPS

So to make certain that you have absolutely comparable settings and parameters, you have to make sure that everything that would affect the stress calculation has to be well converged (e.g. a PPPM convergence of 1.0e-4 is not well converged. it is merely sufficient to have forces that are sufficiently converged for constant volume simulations. good convergence for pressure needs better settings than that, and you are probably also better off to not use tabulation.

So as a starting point, you should run tests with just 2 atoms of opposite charge or no charge at different distances and then determine what settings in both LAMMPS and Gromacs are necessary to have well converged forces and energies and stresses for coulomb and LJ interactions separately and then together. And carefully choose the pair style in LAMMPS, too.

In summary, there is much more that you have to keep in mind, if you want to have an apples to apples comparison.

axel.