Enforce2d not work for compressing rigid

Dear LAMMPS users,

When I use “rigid/npt” to compress the box containing rigid particles (an ellisoid + two spheres) in 2D with “fix enforce2d”, some particles rotate to z-axis direction like the figures show. The interaction I use is gayberen potential.

The version of LAMMPS I use is 2 Aug 2023. And below is my input script. Is there anything I miss?

units           lj
atom_style      hybrid molecular ellipsoid
dimension       2


read_data       init.data

set             type 1 mass 1.0
set             type 1 shape 3.0 1.0 1.0
set             type 2 mass 1.0
set             type 2 shape 0.6 0.6 0.6
velocity        all create 1.0 87287 loop geom

pair_style      gayberne 1.0 1.0 1.0 5
pair_coeff      1 1 1.0 1.0 1 1 1 1 1 1
pair_coeff      1 2 1.0 1.0 1 1 1 1 1 1 
pair_coeff      2 2 1.0 0.6 1 1 1 1 1 1 

neighbor        4.0 bin
neigh_modify    delay 0 every 1 check yes


fix             1 all rigid/npt/small molecule temp 1.0 1.0 1.0 iso 1.0 5.0 1.0 
fix             2 all enforce2d


neigh_modify    exclude molecule/intra all

compute         mytemp all temp/asphere
compute         mypressure all pressure mytemp
fix_modify      1 temp mytemp press mypressure
compute         0 all property/atom quatw quati quatj quatk shapex shapey shapez

dump           1 all custom 10 dump.ellipsoid.equil id type x y z c_0[*]
dump_modify    1 colname c_0[1] quatw colname c_0[2] quati colname c_0[3] quatj colname c_0[4] quatk colname c_0[5] shapex colname c_0[6] shapey colname c_0[7] shapez

thermo_style    custom step temp epair etotal press
thermo_modify   temp mytemp press mypressure
thermo          100
timestep        0.005
run             10000


init.data (562.7 KB)

Dear all,

I add the keyword “torque” in the “fix rigid”. It works.

But if I use “fix rigid/npt/small”, I cannot use the keyword “torque”. So I change to “fix rigid/npt”.

I hope it can be helpful.