[lammps-users] Simulation of CO2 with a piston

Hello everybody,

Since few weeks I try to fix my simulation and I’m stacked.


Bulk of rigid CO2 molecules with temperature of 400K under pressure 250 bars with size around 118 x 118 x 118 A^3 (supercritical conditions). The system as whole is big in X direction (up to 10000A in both directions)(attached figure). The bulk is constrained by a LJ wall in positive X direction which prevent molecules to fly away behind certain region (purple one). On the other side of CO2 bulk I made a dense wall of carbon atoms. The piston is suppose to keep constant pressure when I remove the LJ wall from CO2 to let it defuse. So the piston should move slowly keeping the constant pressure of 250 bars no matter what I’m doing with CO2.


units metal
boundary p p p
atom_style full

read_data system_plus_wall_38_atoms.data
timestep 0.001

#pair_style lj/cut/coul/long 8.775 10.0
#pair_coeff 2 2 0.00681 3.05 9.15
#pair_coeff 1 1 0.00233 2.80 8.40
#pair_coeff 1 2 0.00398 2.925 8.775

pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 2 2 0.00681 3.05
pair_coeff 1 1 0.00233 2.80
pair_coeff 1 2 0.00398 2.925

pair_coeff 2 3 0.00398 2.925 3.283
pair_coeff 1 3 0.00233 2.80 3.143
pair_coeff 3 3 0.00233 2.80 3.143

kspace_style pppm 0.001

group O type 2
group C type 1

group CO2 type 1 2
group wall type 3

neighbor 2.0 bin
neigh_modify every 5 delay 10 check yes

compute Tr CO2 temp/com
compute Tp wall temp/com

fix right_cont_CO2 all wall/lj126 xhi 309.0 0.00398 2.925 3.2 units box pbc yes

fix NVT CO2 rigid/nvt/small molecule temp 400.0 400.0 0.4
fix PISTON wall rigid group 1 wall

velocity CO2 create 400.0 4928459 #zero angular

fix moving wall aveforce 2.2 0.0 0.0 #moving the system at certain direction
fix pNVE wall nve

The rest is dump and run commands.


I calculated with high attention the force I applied on piston which is 2.2 ev/A according to metal units. The piston is freezed due to the fix command. When I run the simulation (with optimized CO2 or not) the piston is pushing the gas really fast like there is no resistance and no pressure from CO2 molecules with LJ wall at the end of CO2 box or without. I simulated CO2 box separately and everything was fine about the pressure and equilibration went well. Could you tell me what can be wrong here? Or maybe I do not use a command as I should. I tried many option and non of them works well. I also was following all the pitson problems discussion and examples you gave but it doesn’t work for me.

Thank you in advance.

All the best,
Katarzyna Walczak.


the force you add with fix aveforce is a per-atom force, so if 2.2 eV/A would be the intended force on the piston, you would have to divide it by the number of atoms to get the value per atom. something like this:

fix moving wall aveforce $(2.2/count(wall)) 0.0 0.0


Thank you Axel for your reply. It seems to work for now. I misunderstood the documentation of LAMMPS. I thought that this force is applied on a defined group of atoms and it is already divided among them and I was looking for a problem somewhere else.

All the best,