I am trying to recreate the density of TIP4P/2005 water at 298 K and 1 bar per this paper. Simulation details are in Section II. I am observing an average density of 0.98 as opposed to the predicted value of 0.9971. This is a relative error of 1.7%. This isn’t too bad but it is large enough to question my inputs.

The system: I have 360 water molecules in a cubic box with dimensions 22x22x22 angstrom^3 and periodic boundary conditions. All water molecules are held rigid and planar using the SHAKE algorithm. I am running the simulation for 15 million steps with a timestep of 0.25 fs. LJ and Coulombic cutoffs are at 8.5 angstroms. LJ, Coulombic, and geometric parameters were taken from the LAMMPS TIP4P page.

I have posted my data file (tip4p_2005.data), settings file (tip4p_2005.settings), and input file (in.tip4p).

With SHAKE you should be able to run with a timestep of 2fs.

What timestep did the authors in your reference paper use and which software?
When running with NPT the kspace convergence matters. What does the reference paper use?

Oh, sorry for not mentioning the LAMMPS version. I am using LAMMPS Sep 29 2021.

Per the paper, they run 15 million steps at 0.25 fs. They used DLPOLY.

A time step of 0.25 fs ensures energy conservation within a
0.05% in a 15 million steps run. For the calculation of the static dielectric constant and the self-diffusion coefficient we have used the molecular-dynamics package DLPOLY.

This is what they state about their NPT simulations:

Unless otherwise stated, the simulations were carried out using the Monte Carlo method at constant pressure and temperature NpT. Isotropic NpT simulations are adequate for the liquid phase while anisotropic Monte Carlo simulations.Parrinello-Rahman-type are required for the solid phases.

Regarding k-space methods, they mention:

In our simulations, the LJ potential was truncated at 8.5 Å. Standard long-range corrections to the LJ energy were added. The Ewald summation technique has
been employed for the calculation of the long-range electrostatic forces. For the real space cutoff we also employed 8.5 Å. The screening parameter and the number of vectors in the reciprocal space considered was carefully selected for each phase.

I am using pppm/tip4p as opposed to Ewald, but that seems to be the version available to LAMMPS.

That is a rather old version. For a “new” study, I would always start with a recent release.

You are comparing Monte Carlo with Molecular Dynamics. Not sure how well those would agree specifically on NPT. Monte Carlo always favors more disordered geometries, while Molecular Dynamics can have a long “memory” of previously ordered structures.

This is contradicting the statement about MC vs. MD. In MC you do not have a timestep. And a timestep of 0.25fs would suggest using fix rigid instead of fix shake.

This conflicts with your choice.

pair_modify mix arithmetic shift yes

Long-range corrections to LJ are enabled with

pair_modify tail yes

What your input does is shift the energy so that it is zero at the cutoff. That is different.

For some reason, I assumed something about their barostat was Monte Carlo as opposed to the whole simulation being Monte Carlo. Strange.

Thank you for the correction regarding the LJ correction. Glad to know that there was no issue with the implementation of SHAKE. I appreciate the prompt responses!

Also note that as stated in the LAMMPS doc, Nosé-Hoover equation require adjustment of the \tau_T ans \tau_P parameters with regard to the timestep. In your case they end up larger by a factor 4 compared to the recommended values.

With a relative error of 1.7\%, and such a small system, fluctuations come into play. You might want to compare the standard error and standard deviation from your data to check how far you are from you reference value.

Because having correlated data, you will also have to look into statistical efficiency. The standard error may be underestimated. For that you compute the average and standard error for the entire data set, then two halves, four quarters and so on. That will also give you a gauge if there are some slow moving trends. There is a section in the Allen and Tildesley book discussing this.