Reproducing TIP4P/2005 water properties

User discussions
1

GROMACS version: 2024.2 / 2025.2 (2024.2 for production, 2025.2 for analysis)
GROMACS modification: No

Hi,

I am trying to reproduce the TIP4P/2005 water model (ref: DOI 10.1063/1.2121687) properties at ambient conditions (temperature 298 K, 1 bar pressure), to assure that my setup is correct. I got the topology file and initial mdp-file from the Carlos Vega website (catalan.quim.ucm.es) changed mdp-file according to Gromacs recommendations (nstcomm changed to default, nstlist from 1 to 10, changing thermostat to v-rescale with tau-t = 1 ps, barostat to C-rescale with tau-p = 5 ps). I also changed rcoloumb and rvdw to 0.9 (instead of 0.85), as I have seen later publications from the same group use 0.9 (DOI: 10.1080/00268970902784926). Attached topol file and mdp file.

topol_C_Vega.top (1.4 KB)

prodparams_vega_v3.mdp (3.9 KB)

After energy minimization, NVT-equilibration (300 ps) and NPT-equilibration (3 ns), I run three NPT replicates for 300, 100, and 100 ns respectively (each having different starting velocities) with a 5 nm cubic box (4139 molecules), I get the following results (select properties):

Simulation Density(g/cm^3) k_T(10^5 MPa^-1) Cp (cal mol^-1 K^-1) Epsilon Diffusion(10^9 m^2/s)

Rep1 (300 ns) 0.9971 51.4 21.3 57.3 2.16

Rep2 (100 ns) 0.9971 51.9 21.3 57.5 2.19

Rep3 (100 ns) 0.9971 51.6 21.6 57.4 2.14

Reference 0.9979 46.5 21.1 60 2.08

The k_T (isothermal compressibility) seems to me to be too far off. I compute it through the volume fluctuation formula (both through a script of my own, and using Gromacs fluct_props option in gmx energy - they yield similar results). I am aware that fluctuation properties are slow to converge, but surely, with multiple replicates and 100 ns simulation, convergence should not be the issue?
As well, while the original publication used Nosé-Hoover thermostat, I read (DOI: 10.1016/j.molliq.2022.120116) that Nosé-Hoover and V-rescale should yield similar results.

Thanks in advance for any suggested improvements or any spotted errors.

2

I see that this setup uses normal constraints with LINCS. Why is that? SETTLE is more accurate and more efficient. I don’t think this should make a difference, but it would be good to check.

3

Thanks for your reply.

The reason for choosing LINCS was because the reference paper was using LINCS. Will try SETTLE though - if nothing else, I might gain some efficiency.

I have now also tried changing the barostat to Parrinello-Rahman, with tau_p = 2 ps. After a 300 ns simulation (with equilibration steps as before), I get most properties similar to using C-rescale, except that now I get isothermal compressibility of 91.4 (10^5 MPa^-1) - quite different from earlier.

4

Are those compressibility values converged? I would expect both C-rescale and PR to give correct ensembles, if tau_p and nstpcouple are not set to unreasonable values.

5

Concerning convergence:

For C-rescale:
I would refer to the data in the first post - three different replicates, each min. 100 ns long, all yield results within 51-52 (10^5 MPa^-1), which I believe is converged. I can also report that splitting each of those in half, each half differ from each other by no more than 1.2 (10^5 MPa^-1).

For PR:
For a quick-and-dirty convergence check, splitting the trajectory into 5 parts (i.e. each 60 ns), gives the following results:
0-60 ns: 91.6
60-120 ns: 91.2
120-180 ns: 91.7
180-240 ns: 90.6
240-300 ns: 91.7

I think the simulations have a level of convergence good enough to indicate clear differences when using different barostats, even if this analysis is crude.

6

That’s pretty worrying then. Both barostats should give the same results. Are all other mdp options identical?

7

Yes. Attached files for the production runs.

prodparams_vega_v3_PR.mdp (3.8 KB)

prodparams_vega_v3.mdp (3.9 KB)

8

I think that I found the issue. With tau_p=5 ps, the Parrinello-Rahman barostat needs to be updated at least every 10 steps, not the default 100. This is controlled by the nstpcouple mdp option. The error can be very large for small systems and gets smaller as the system size increases.

Currently I don’t know if this is due to inaccuracy of the integration or due to a bug in the code. C-rescale seems to behave correctly.

9

Using settle and nstpcoule = 10, with C-rescale, I get the following results for a 120 ns simulation:

Density = 0.9971, target = 0.9979 (g/cm^2)
Isothermal compressibility = 46.6, target = 46.5 (10^5 MPa^-1)
Heat capacity = 21.3, target = 21.1 (cal mol^-1 K^-1)
Epsilon = 56.86, target = 60
Diffusion = 2.20, target = 2.08 (10^9 m^2/s)

Value for isothermal compressibility is now much better, setting nstpcouple = 10 would seem to be the solution.
Minor deviation for diffusion is probably due to lack of Ye-Hummer correction, which are expected given differences in box sizes.

I consider the issue settled. Thanks a lot for the help!

10

Note that I still get an increase of 2-3% of the RMSF of the volume for a system of 400 water molecules and tau_p=5 and nstpcouple=10. You need to use even lower nstpcouple if want results more accurate than that for small systems. As the system size increases the error goes down.

11

I opened an issue, with high priority: Too large nstpcouple value for Parrinello-Rahman barostat increases volume fluctuations (#5424) · Issues · GROMACS / GROMACS · GitLab

I did tests for the worst case: a small water box. I uploading a fix for the release 2025 branch. This gives nstpcouple=10 with tau_p=5.