I want to simulate two gases contained in an adiabatic cylinder and separated by a moving piston (something like this http://www.opensourcephysics.org/items/detail.cfm?ID=10260), and I want the piston made of particles. All the interactions beetwen particles (gas and piston) are of the LJ type and, to keep togheter the particles composing the piston, I should introduce for each particle an external potential (eg harmonic) with the equilibrium position depending on the macroscopic position of the piston (eg the center of mass of the particles composing the piston), thus the equilibrium positions of each potential will depend on time as the piston is free to move along the axis of the cylinder.

Do you know how I can do this in lammps?


As a first step, you can treat the piston as a rigid body

and constrain its motion to the x direction. See the
fix rigid command. Once that works, you can
add non-rigid atoms to each piston atom that are
attached by harmonic bonds, and turn off the rigid

piston atoms with neigh_modify so they are
not seen by the gas atoms.



thanks for the reply. Created the rigid body but now I have to fix new atoms to the rigid ones.
I had a look in the documentation and some examples but it still not clear how to do this.

This is the script I use now for the gas+ rigid piston

# 2d points

units lj
dimension 2
atom_style atomic
boundary f f p

read_data data.atomic

group pistone type 2
group gas type 1

mass 1 1.0
mass 2 1.0

fix xwalls all wall/reflect xlo EDGE xhi EDGE
fix ywalls all wall/reflect ylo EDGE yhi EDGE

velocity gas create 1.44 87287 loop geom
velocity pistone create 0 87287 loop geom

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0
pair_coeff 1 2 1.0 1.0
pair_coeff 2 2 0.5 1.0

neighbor 0.3 bin

fix 1 all nve
fix 2 all enforce2d

timestep 0.003
thermo_style custom step time temp ke pe etotal press vol
thermo 500
run 10000


Quoting Steve Plimpton <[email protected]>:

Like any system with fixed bonds, you need to
create a data file that lists the bonds, as well
as the coords of each pair of bonded atoms.

I.e. the initial geometry and topology of your


So should I use

atom_style bond

instead of atomic and define all the bonds in the data file?

So should I use

atom_style bond

instead of atomic and define all the bonds in the data file?

lammps cannot store (and use) information about bonds unless the atom
style supports it. so you do need atom style bond, angle, molecular,
or full.

and defining bonds in a data file, is usually the simplest way to let
LAMMPS know which bonds you want to go where.



I modified the script and now I can simulate my system.

I have 3 groups of atoms, 1 for the gas, 2 for the rigid piston and 3 for the flexible piston attached to the rigid one, and one kind of bond, between atoms of type 2 and 3 (harmonic bond). The pairwise interaction is lj/cut.

I fix the rigid piston with

fix 3 piston rigid single force * off off off torque * off off off

for this test it cannot move at all. Then I use

neigh_modify exclude type 1 2
neigh_modify exclude type 2 2

and I integratethe gas+flexible piston with

fix 1 gas_pistonf nve

But then I saw in the output that the total energy (etotal) increases systematically. And I cannot explain why.

Quoting Axel Kohlmeyer <[email protected]>:

Could be any # of reasons for increasing energy.
E.g. too big a timestep.