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 boxmass 1 15.9994
mass 2 1.008
mass 3 1.0e-100pair_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.0fix 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 alltimestep 0.5
fix integrate all rigid/small molecule langevin 300.0 300.0 100.0 2345634dump 1 all atom 1000 traj.lammpstrj
thermo_style custom step temp press etotal density pe ke
thermo 1000run 100000
NPT:
units real
atom_style charge
atom_modify map array
region box block -15 15 -15 15 -15 15
create_box 3 boxSet 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 atomsSpecify 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 allStep 1: Use NVT ensemble to stabilize temperature
timestep 0.5
fix nvt all rigid/nvt/small molecule temp 300.0 300.0 100.0Thermodynamic output for NVT stage
thermo_style custom step temp etotal ke pe
thermo 1000
dump 1 all atom 5000 traj.lammpstrjRun 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.0Thermodynamic output for NPT stage
thermo_style custom step temp press etotal density pe ke
thermo 1000run 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