Hello,
I’m working on a NPT simulation where a few hydrogen molecules (represented as beads with a mass of 2) are solvated in TIP4P water (rigid, 4-point geometry). Since my system involves a hybrid model with both rigid bodies (water) and non-rigid particles (hydrogen beads), I followed the LAMMPS documentation for the fix rigid command (link) and applied the following commands for temperature and pressure control:
# Apply fix rigid/npt/small to water molecules, dilating all atoms when box changes
fix rigid_water water rigid/npt/small molecule temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 dilate all
# Apply fix nvt to hydrogen beads
fix nvt_h2 hydrogen_beads nvt temp 300.0 300.0 100.0
The issue is that after a few hundred timesteps, the simulation becomes unstable. The error I get varies depending on the number of CPU cores or ensemble I use, but they typically include:
ERROR on proc 0: Non-numeric box dimensions - simulation unstable (src/KSPACE/pppm.cpp:1803)
or
ERROR: Out of range atoms - cannot compute PPPM
or
ERROR on proc 1: Non-numeric box dimensions - simulation unstable
I’ve checked the velocities of the particles (as Axel recommended Here, and it seems that the hydrogen beads have significantly higher velocities compared to the water molecules.
So far, I have tried:
- Reducing the timestep to very small values.
- Using an NVT ensemble instead of NPT (as per the documentation), but the issue persists.
I’m wondering if anyone can point out what might be causing the instability.
Any suggestions or insights would be greatly appreciated!
I’ve attached the input file for reference.
Thanks in advance!
Ehsan
# LAMMPS Input Script for Water and Hydrogen Beads Simulation with Explicit 4-Point TIP4P
units real
atom_style charge
atom_modify map array
boundary p p p
# Define the simulation box
region box block 0 15 0 15 0 15
create_box 4 box
# Define masses
mass 1 15.9994 # Oxygen in water
mass 2 1.008 # Hydrogen in water
mass 3 1.0e-6 # M site (tiny non-zero mass)
mass 4 2.01588 # Hydrogen bead (H2 molecule)
# Define pair style and kspace style
pair_style lj/cut/coul/long 10.0
kspace_style pppm 1.0e-4
# Define pair coefficients
pair_coeff 1 1 0.1550 3.1536
pair_coeff 2 2 0.0 1.0
pair_coeff 3 3 0.0 1.0
pair_coeff 4 4 0.06796 2.96
# Cross interactions between hydrogen beads and water molecules
pair_modify mix arithmetic
# Read the water molecule file (ensure the path is correct)
molecule water tip4p.mol
# Use fix property/atom mol to store molecule IDs
fix prop all property/atom mol
# Create hydrogen beads first to avoid overlaps with water molecules
create_atoms 4 random 10 67890 box overlap 2.5 maxtry 10000
# Now create water molecules using the molecule template
create_atoms 0 random 100 12345 box mol water 25367 overlap 2.5 maxtry 10000
# Define groups
group water type 1 2 3
group hydrogen_beads type 4
# Exclude intra-molecular interactions within water molecules
neigh_modify exclude molecule/intra water
# Initialize velocities (adjust temperature as needed)
velocity all create 300.0 5463576
# Set timestep
timestep 0.5
# Apply fix rigid/npt/small to water molecules, dilating all atoms when box changes
fix rigid_water water rigid/npt/small molecule temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 dilate all
# Apply fix nvt to hydrogen beads
fix nvt_h2 hydrogen_beads nvt temp 300.0 300.0 100.0
# Output settings
thermo_style custom step temp press etotal pe density
thermo 100
# Run the simulation
run 200000