fix shake and angle harmonic

Dear professor,
i want to obtain some transform properties of CO2 use the classical EPM2 model,

a model with rigid bond lengths but with a flexible bond angle potential(harmoonic)

i always cannot obtain the right density and i have read some advice about fix shake(which is not fit to linear molecule but i only shake the bond of CO2 not angle).could you give me some suggestion,thank you !

here are my script.

variable t equal 423

#setup
units real
dimension 3
atom_style full

neighbor 1.7 bin
neigh_modify delay 0 every 1 check yes

read_data data.1000co2

mass 1 12.0107 # C
mass 2 15.9994 # O

group Cgroup type 1

group co2 type 1 2

pair_style lj/cut/coul/long 12.0 8.5

pair_modify tail yes
kspace_style pppm 1.0e-5

pair_coeff 1 1 0.05587 2.757
pair_coeff 1 2 0.0945 2.8917
pair_coeff 2 2 0.1599 3.033

bond_style harmonic

bond_coeff 1 0 1.149

angle_style harmonic

angle_coeff 1 14.77 180

velocity all create $t 22222222 dist gaussian

fix 1 all shake 0.000001 100 0 b 1

fix 3 all nve
fix 4 all temp/berendsen $t t (100*dt)

thermo_style custom step temp ke pe etotal press vol

timestep 1

thermo 10000
run 1000000

unfix 3
unfix 4

reset_timestep 0

fix 5 all npt temp $t t (100dt) iso 246.73 246.73 $(1000dt)
thermo_style custom step temp ke pe etotal press vol
timestep 1

thermo 10000
run 1000000

unfix 5
reset_timestep 0

fix 6 all npt temp $t t (100dt) iso 246.73 246.73 $(1000dt)

thermo_style custom step temp ke pe etotal press density vol

dump 3 Cgroup custom 50 CARBON.lammpstrj id type vx vy vz
dump_modify 3 sort id

timestep 1

thermo 20000
run 3000000

What is the “right” density and how was that determined, and what density are you getting instead?

It is unlikely that any discrepancy may be caused by fix shake. I would first suspect the non-bonded force field parameters. Do you have any information of the error margin to expect. Please also keep in mind, that an MD simulation like LAMMPS is performing is usually employing settings and parameters for condensed systems. So looking reproduce gas phase properties is more challenging and often not even worth doing. In many cases it is better to just set the volume/density to the desired value, run at constant volume and be done with it.

Axel.