Piling of supequadric particles

I'm trying to run a simulation like that in Podlozhnyuk et al. (2017) Fig. 7, where superquadric particles are piled up in a container. However, in a small test with 138 I find that the particles fall towards the bottom surface then start to 'dance around' and they never come to rest (see attached animation). I've tried altering the simulation parameters, including adding coefficient of rolling friction, but the problem still occurs. Could anyone suggest what might be incorrect or missing in the script (below).

Many thanks !

Liggghts version 3.8.0

variable box_size equal 48.0e-3
variable cutoff equal 4.0e-3

variable density0 equal 1000

variable N_in equal 138
variable Np_in equal 200

variable vel1_in equal 0.0
variable vel2_in equal 0.0
variable vel3_in equal 0.0

variable shape1 equal 1.6e-3
variable shape2 equal 1.2e-3
variable shape3 equal 1.2e-3

variable blockiness1 equal 2.0
variable blockiness2 equal 2.0

# coefficientRestitution
variable cf_ww equal 0.5
variable cf_pw equal 0.5
variable cf_pp equal 0.5

# coefficientFriction
variable cof_ww equal 0.5
variable cof_pw equal 0.5
variable cof_pp equal 0.5

# coefficientRollingFriction
#variable corf_ww equal 0.0
#variable corf_pw equal 0.0
#variable corf_pp equal 0.0

# youngsModulus
variable Yw equal 2.5e+7
variable Yp equal 2.5e+7

####

variable xmin equal -${box_size}/2
variable xmax equal ${box_size}/2
variable ymin equal -${box_size}/2
variable ymax equal ${box_size}/2
variable zmin equal -${box_size}/2
variable zmax equal ${box_size}/2

####

echo both
log log.liggghts

atom_style superquadric
atom_modify map array
boundary p p m
newton off

communicate single vel yes

units si
#processors 2 2 1

region reg block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} units box
create_box 2 reg

neighbor ${cutoff} bin
neigh_modify delay 0

# Material properties required for granular pair styles

fix m1 all property/global youngsModulus peratomtype ${Yw} ${Yp}
fix m2 all property/global poissonsRatio peratomtype 0.3 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 2 ${cf_ww} ${cf_pw} ${cf_pw} ${cf_pp}
fix m4 all property/global coefficientFriction peratomtypepair 2 ${cof_ww} ${cof_pw} ${cof_pw} ${cof_pp}
#fix m5 all property/global coefficientRollingFriction peratomtypepair 2 ${corf_ww} ${corf_pw} ${corf_pw} ${corf_pp}
fix m6 all property/global characteristicVelocity scalar 2.0
#fix m7 all property/global coefficientRollingViscousDamping peratomtypepair 2 0.1 0.1 0.1 0.1

# pair style
pair_style gran model hertz tangential history #Hertzian without cohesion
pair_coeff * *

timestep 0.000001

# walls
fix zwalls1 all wall/gran model hertz tangential history surface superquadric primitive type 1 zplane ${zmin}
#fix zwalls1 all wall/gran model hertz tangential history rolling_friction epsd2 surface superquadric primitive type 1 zplane ${zmin}

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

###### particle distributions and insertion
region bc block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} units box

fix pts1 all particletemplate/superquadric 15485863 atom_type 2 density constant ${density0} &
shape constant ${shape1} ${shape2} ${shape3} &
blockiness constant ${blockiness1} ${blockiness2}

fix pdd1 all particledistribution/discrete 15485867 1 pts1 1.0

fix ins all insert/pack seed 32452843 distributiontemplate pdd1 &
orientation random &
vel constant ${vel1_in} ${vel2_in} ${vel3_in} &
insert_every once overlapcheck yes all_in yes particles_in_region ${N_in} region bc

# apply nve integration to all particles that are inserted as single particles
fix integr all nve/superquadric integration_scheme 1

# output settings, include total thermal energy
compute rke all erotate/superquadric
thermo_style custom step atoms c_rke cpu time
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

#insert the first particles
run 1
dump dmp all custom/vtk 1000 post/stream_*.vtk id type x y z ix iy iz vx vy vz fx fy fz &
omegax omegay omegaz shapex shapey shapez blockiness1 blockiness2 quat1 quat2 quat3 quat4 &
density mass
dump_modify dmp binary yes

run 1000000