[lammps-users] introducing an interface

Dear LAMMPS users,

I am interested in simulating a liquid/vapor interface, but I am running into a strange problem. My system consists of N=256 rigid molecules, each molecule consisting of three Lennard-Jones atoms. Atoms within a molecule have no interactions with each other. Only different molecules interact with each other.

I have successfully equilibrated my system in a cubic volume, but when I introduce an interface by expanding the zlo, zhi values and run fix rigid, my system breaks down: the total energy is no longer conserved the kinetic energy increases without bound. (Sometimes the total energy decreases before increasing later.) The potential energy, however, does not change much. When I take the same input configuration file but change the zlo and zhi values to their original values, the system is stable.

I've tried the following:

- I've checked for rigid molecules across the z-boundaries, and after fixing that so that the molecule is 'whole' in the configuration file, I find that the system blows up
- I've simulated a system of 256 LJ monomers, and this system is stable after the expanding the zlo,zhi values.
- I've compiled lammps on a different architecture machine (64bit versus 32 bit) and the systems blow up

I've actually looked at several systems, with both charged and uncharged atoms, and of various system sizes (N=256, 500, 864). What is confusing is that my systems do not always blow up; only some blow up but others can be perfectly stable, exhibiting the expected behavior after the introduction of an interface. I basically cannot predict whether a system will blow up or not upon the introduction of an interface, despite the many different experiments I've run.

Based on my studies, I think either there is some implicit assumption somewhere in lammps that I am violating by extending zlo and zhi (or any boundary since I've also tried changing x and y boundaries), or maybe there is a bug somewhere??

Attached are the configuration file and input file for the simplest system that behaves strangely when an interface is introduced.

I would really appreciate any input on this! Thanks in advance,

in.sample (606 Bytes)

sampleConfig.txt (72.6 KB)

HI Joyce. If the product of particle speed and timestep size gets too
big in any MD simulation, unstable integration can occur and energy is
not conserved. Trying to move particles too far in a single timestep
is bad. After playing with your problem quite a bit, it looks to me
like that is the problem you have here. When I reduce the particles'
speed or the timestep size, the simulation conserves energy much
better. For example, try this input script where I've done both:

units lj
atom_style full
pair_style lj/cut 2.5
read_data joyce.data
velocity all create 0.1 4928459
timestep 0.0005
thermo 100
thermo_style custom step temp pe ke etotal press
fix 1 all rigid molecule
run 10000

I didn't change your data file at all. Of course, you don't want to
overdo it and use a tiny timestep because you'll waste CPU time. You
have to find that balance between speed and good energy conservation.

I'd also recommend using shake to make your molecules rigid instead of
using fix rigid. You'll have to create bonds and angles and then apply
a shake fix on the bonds and the angles. This will likely help achieve
more stable and faster integration.