Issue with TIP4P Water Simulation

I am simulating water molecules using the TIP4P model with the 4-point geometry (approach 2, as explained Here. However, when using a larger box (30 Å), the water molecules seem to behave like a solid. They mainly jiggle in place without significant diffusion, even when running the simulation for up to 3 ns.

I tried both NVT and NPT ensembles (with a short NVT equilibration), but the results remain the same, or even worse under NPT, where the mobility is even more restricted. The molecules don’t seem to diffuse like you’d expect in bulk water.

Interestingly, I ran the same simulation in GROMACS using TIP4P water with the OPLS-AA force field under the same conditions (NPT, same timestep), and the water molecules exhibited significantly more movement, behaving more like a liquid with normal diffusion behavior.

I’m wondering if I’m missing something in my LAMMPS setup. I’ve attached both the NVE and NPT input scripts alongside MDP file for GROMCS. I have also attached a comparison trajectory from LAMMPS (NVT) vs GROMACS NPT for 75 ps.

Any advice on whatm could be causing this issue would be appreciated!

NVT:

units real
atom_style charge
atom_modify map array
region box block -15 15 -15 15 -15 15
create_box 3 box

mass 1 15.9994
mass 2 1.008
mass 3 1.0e-100

pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.1550 3.1536
pair_coeff 2 2 0.0 1.0
pair_coeff 3 3 0.0 1.0

fix mol all property/atom mol ghost yes
molecule water tip4p.mol
create_atoms 0 random 909 34564 NULL mol water 25367 overlap 1.33
neigh_modify exclude molecule/intra all

timestep 0.5
fix integrate all rigid/small molecule langevin 300.0 300.0 100.0 2345634

dump 1 all atom 1000 traj.lammpstrj

thermo_style custom step temp press etotal density pe ke
thermo 1000

run 100000


NPT:

units real
atom_style charge
atom_modify map array
region box block -15 15 -15 15 -15 15
create_box 3 box

Set masses based on TIP4P/2005

mass 1 15.9994 # O
mass 2 1.008 # H
mass 3 1.0e-100 # M (tiny non-zero mass for M site)

Set pair style and pair coefficients for TIP4P/2005

pair_style lj/cut/coul/cut 8.0
pair_coeff 1 1 0.1852 3.1589 # Lennard-Jones params for O atoms (TIP4P/2005)
pair_coeff 2 2 0.0 1.0 # No LJ interactions for H atoms
pair_coeff 3 3 0.0 1.0 # No LJ interactions for M atoms

Specify molecular properties for rigid bodies

fix mol all property/atom mol ghost yes
molecule water tip4p.mol # Use the molecule file for TIP4P/2005 geometry
create_atoms 0 random 900 34564 NULL mol water 25367 overlap 1.33
neigh_modify exclude molecule/intra all

Step 1: Use NVT ensemble to stabilize temperature

timestep 0.5
fix nvt all rigid/nvt/small molecule temp 300.0 300.0 100.0

Thermodynamic output for NVT stage

thermo_style custom step temp etotal ke pe
thermo 1000
dump 1 all atom 5000 traj.lammpstrj

Run the NVT ensemble for 50,000 timesteps to stabilize temperature

run 50000

Step 2: Unfix the NVT and switch to NPT ensemble to adjust density

unfix nvt
fix npt all rigid/npt/small molecule temp 300.0 300.0 100.0 iso 1.0 1.0 500.0

Thermodynamic output for NPT stage

thermo_style custom step temp press etotal density pe ke
thermo 1000

run 100000


Gromacs md parameters:

integrator  = md
dt          = 0.0005
tinit       = 0
nsteps      = 150000     ; 100 ps 
nstcomm     = 100
; Output parameters
nstvout             = 5000  ; every 10 ps 
nstfout             = 5000
nstxout-compressed  = 1000
nstenergy           = 5000
; Bond parameters
constraints = h-bonds
constraint-algorithm = LINCS
lincs-order = 4  ; typical for standard MD
lincs-iter = 1   ; usually sufficient for MD
shake-tol = 0.0001
continuation            = no 
; Single-range cutoff scheme
cutoff-scheme   = Verlet
nstlist         = 20
ns_type         = grid 
rlist           = 1.4
rcoulomb        = 1.4
rvdw            = 1.4
; PME electrostatics parameters
coulombtype     = PME
fourierspacing  = 0.12
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
pme_order       = 4
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on in two groups
Tcoupl      = nose-hoover 
tc_grps     = system 
tau_t       = 0.5    
ref_t       = 300  
; Pressure coupling is on
Pcoupl          = C-rescale 
pcoupltype      = isotropic
tau_p           = 1.0       
compressibility = 4.5e-5
ref_p           = 1.0
refcoord_scaling = com
; Generate velocities is on 
gen_vel     = yes 
gen_temp    = 300
; Periodic boundary conditions are on in all directions
pbc     = xyz
; Long-range dispersion correction
DispCorr    = EnerPres

GROMACS trajectory

LAMMPS trajectory

You not only have different short-range cutoffs – you have a completely different electrostatic model between your GROMACS (long-range via PME) and LAMMPS (no long-range electrostatics) inputs.

Dear Shern,

Thank you for your reply. I understand that comparing LAMMPS and GROMACS isn’t directly meaningful since different force fields and electrostatic models are being used. However, my primary concern is not about the differences between the two software packages, but rather about why the water molecules in LAMMPS are not moving as expected and appear to behave like a solid.

I’m following the LAMMPS documentation for the TIP4P water model quite closely, and the only change I made was increasing the size of the simulation box. Despite that, the water molecules are exhibiting minimal diffusion, which is unusual for liquid water.

Could anyone help clarify why this behavior might be occurring?

Yes – I can see now that the LAMMPS documentation at 8.4.4. TIP4P water model — LAMMPS documentation is currently in error. You should simulate TIP4P in LAMMPS with a long-range electrostatics solver. For your setup it would be

pair_style lj/cut/coul/long 8.0
kspace_style pppm 1e-4

I’ll fix up the documentation.

Having said that, it is still very important to know what exactly you are simulating in LAMMPS – and GROMACS for that matter (which would give you the same strange behaviour if you had turned off the long-range electrostatics).

The long-range electrostatic interaction enables long-range interactions between water molecules and thus additional thermalization in the associated vibrational modes. As you have discovered, without those modes, the TIP4P water model’s liquid state evidently loses a lot of vibrational entropy, and thus the solid to liquid transition temperature must correspondingly have been raised past 300K.

Even in Jorgensen’s original paper in 1983(!!), the molecular dynamics of the model were simulated with long-range electrostatics via Ewald summation. Monte Carlo of the static, structural properties was done without long-range electrostatics – that’s reasonable (especially for those times) since those interactions would average to zero structurally.

2 Likes

Dear Shern,
Tthank you so much! Switching to long-range electrostatics with pair_style lj/cut/coul/long and kspace_style pppm fixed the issue—now the water behaves as expected. Your explanation about the importance of long-range electrostatics for proper thermalization makes perfect sense. Thanks also for addressing the documentation!