[lammps-users] Simulation of ice slabs and molecules

Dear All,

I am trying to simulate H2 molecules interacting with a slab of water ice.

So, I first prepare a slab of ice, which works fine. But, then when I introduce
even a single molecule of H2 into the system, the who systems blows up.
In the sense that, the water molecules are no more in a slab configuration, but
fill up the entire box and eventually i get an error.

I have checked my H2+ H20 interaction potential n number of times, but couldn’t
get to solve this problem. I guess, most of you simulate similar systems. So
has anyone encountered such problem? Is there anything in particular I need to
look into.


Could be any number of reasons, most of them involving
a bad model or bad params in your input script. Without
details it's impossible to tell.


Thanks Steve for a quick reply. I have added the input script in this mail.
In one of my test, I even commented out all the code in the new potential
and compiled it. Although there are no forces on the H2, but the water molecules
still expand and go out of the box. So, that means, it was not the interaction
of H2-H20 (as the new code was completely commmented out!) but something else.

Any inputs will be highly appreciated.
Thanks a ton!


tip4p water and H2 MD calculation

units real
atom_style full
boundary p p f
echo both
pair_style hybrid lj/cut/coul/long/tip4p 1 2 1 1 0.125 15.0 lj/cut/coul/long/h2water 1 2 3 1 1 1 0.125 15.0
pair_modify shift yes mix arithmetic
kspace_style pppm/tip4p 0.00001
kspace_modify slab 3.0
bond_style harmonic
angle_style harmonic
read_data data2.slab

pair_coeff 1 1 lj/cut/coul/long/tip4p 0.16275 3.16435
pair_coeff 1 2 none
pair_coeff 2 2 none

pair_coeff 1 3 lj/cut/coul/long/h2water 1.1425 3.10
pair_coeff 2 3 lj/cut/coul/long/h2water 0.5112 2.15
pair_coeff 3 3 none

bond_coeff 1 80.0 0.9572
#bond_coeff 2 80.0 0.7400

angle_coeff 1 80.0 104.52

neighbor 2.0 bin
neigh_modify delay 5
timestep 1.0
thermo_style custom step temp epair press vol lx ly lz
thermo 0

group water type 1 2
group H2 type 3

reset_timestep 0

#velocity H2 create 10.0 238946237 mom yes dist gaussian
#velocity water create 10.0 238946237 mom yes dist gaussian

fix 1 water nvt temp 10 10.0 5.0
fix 2 water shake 0.00001 200 0 b 1 a 1
fix 3 H2 rigid/nve single

thermo_style custom step temp epair etotal press vol lx ly lz
thermo 100
thermo_modify flush yes

dump 2 H2 xyz 100 testh.xyz
dump 1 all xyz 1000 test.xyz
run 5000
write_restart test.restart

I'm not sure what you're asking. I would model the water
by itself, first with NVE, then with NVT, then with SHAKE
and see if you can get a stable system in each case for
a longish run.


Hi Steve,

Thanks for your previous replies.

I have modeled TIP4P water by first with NVE, then with NVT and so on to get a stable
system. I could get the properties of the ice slab which I was looking for. This I have been
able to do earlier too.

But, the moment I add a H2 particle into this system. I use pair_style hybrid, the
pressure of the system blows up. The energy and the temperature of the system are good
but the pressure values are way higher, eventually breaking away the slab.

I even turned off all the interaction terms between H2 and H2O molecules, but this didn’t help either.

I looked into a very old thread of the users group
where a similar issue was reported, although I am not sure how to solve it.

Let me know if I am missing something or you need further details about the problem.


Can you post input scripts (and data files) for small
problems where you see a dramatic pressure difference
with and without a single H2 molecule?