Need help in LAMMPS: MD simulation on Water

Dear Altruist,
I am trying to work on finding out the pressure vs density relationship for water molecules. Here, I used a data file for water molecules (containing 1500 molecules and density of 1000Kg/m3). The potential forcefield was chosen to lj interaction with coulombic range. But, when I run the simulation with nvt, then I found that pressure was changing, but density was remained constant. Also, the temperature was not found constant at 300K. Would you like to help how could I get good pressure vs density curve from this input script? Where is the problem in this input script? I really need help in this case as I am new to learning simulation.
Sincerely,
Akash Talapatra
MS student, Virginia Tech
#Input script

units real
atom_style full

boundary p p p
read_data water(1000).data

replicate 1 1 1

variable dt equal dt
variable ext_temp equal 300
mass 1 15.9994
mass 2 1.008
variable Tdamp equal 100.0*${dt}

pair_style lj/cut/coul/long 9.0 9.0
pair_modify tail yes
kspace_style pppm 1.0e-4
pair_coeff 1 1 0.15539421659476232 3.16555789
pair_coeff * 2 0.0000 0.0000

bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none

bond_coeff 1 5000.00 1.000
angle_coeff 1 500.0 109.47
#Create initial velocity distribution
velocity all create ${ext_temp} 432567 dist uniform

#Reset the time step counter
reset_timestep 0

#SHAKE Parameters to preserve bondlengths
fix 1 all shake 0.0001 20 0 b 1 a 1
#Ensemble set-up: NVT ensemble, thermostatted by Nose-Hoover a chained thermostat, set at temperature ext_temp fix 2 all nvt temp {ext_temp} {ext_temp} {Tdamp}
#timestep 0.5
#Output Style
variable Nevery equal 20
variable Nrepeat equal 50
variable Nfreq equal {Nevery}*{Nrepeat}
variable PotentialEnergy equal epair
variable Pressure equal press
variable Temperature equal temp
variable Density equal density

fix 3 all ave/time {Nevery} {Nrepeat} ${Nfreq} v_Temperature v_PotentialEnergy v_Pressure file ave.dens_1000kgm3.out format %.8g
dump st all custom 100 dump.lammpstrj id type x y z
thermo 1000
thermo_style custom step temp density epair press
#Run it!
run 60000

#Outputs
Step Temp Density E_pair Press
0 450.05002 1.0000059 -9823.1202 15649.381
1000 297.78229 1.0000059 -16549.086 -681.74344
2000 296.56455 1.0000059 -16723.059 386.22762
3000 297.68996 1.0000059 -16781.364 462.93578
4000 297.96643 1.0000059 -16689.808 531.10791
5000 300.97047 1.0000059 -16693.015 -419.22291
6000 302.47401 1.0000059 -16838.386 131.47509
7000 298.78537 1.0000059 -16798.136 124.73786
8000 296.31094 1.0000059 -16759.38 80.834287
9000 296.9679 1.0000059 -16783.12 -302.43309
10000 305.12358 1.0000059 -16772.815 152.49701
11000 301.0554 1.0000059 -16749.773 81.422811
12000 295.94107 1.0000059 -16727.506 -111.04608
13000 293.41613 1.0000059 -16711.442 -543.91158
14000 303.78818 1.0000059 -16728.467 -337.74409
15000 302.84376 1.0000059 -16702.389 -70.632047
16000 301.97565 1.0000059 -16729.029 494.58834
17000 304.80025 1.0000059 -16698.838 599.31846
18000 298.86411 1.0000059 -16781.149 405.64783
19000 303.71063 1.0000059 -16716.012 89.099738
20000 303.65877 1.0000059 -16814.3 -195.30294
21000 306.9728 1.0000059 -16768.676 146.53848
22000 305.16509 1.0000059 -16767.722 -330.66466
23000 311.26889 1.0000059 -16642.254 135.18624
24000 303.30852 1.0000059 -16807.744 67.820395
25000 310.80295 1.0000059 -16790.206 868.86566
26000 300.21857 1.0000059 -16601.377 20.002035
27000 297.02305 1.0000059 -16534.666 371.72524
28000 293.81809 1.0000059 -16635.014 -136.15918
29000 295.67957 1.0000059 -16731.441 493.57746
30000 293.33474 1.0000059 -16693.485 -248.27554
31000 302.17025 1.0000059 -16748.443 331.17401
32000 298.32418 1.0000059 -16629.347 161.55044
33000 302.9666 1.0000059 -16681.582 -63.021334
34000 300.54374 1.0000059 -16817.865 -67.415523
35000 305.49678 1.0000059 -16828.5 62.712364
36000 296.90462 1.0000059 -16726.777 235.69704
37000 299.64712 1.0000059 -16855.124 314.17953
38000 289.29864 1.0000059 -16753.823 67.643625
39000 298.25035 1.0000059 -16696.994 376.87236
40000 289.85982 1.0000059 -16814.177 327.49061
41000 293.26169 1.0000059 -16689.921 379.2272
42000 296.04682 1.0000059 -16816.177 -602.71629
43000 307.06156 1.0000059 -16735.833 187.46977
44000 303.15135 1.0000059 -16829.363 417.48567
45000 304.24218 1.0000059 -16712.185 -435.83222
46000 296.84282 1.0000059 -16734.585 138.69402
47000 295.51977 1.0000059 -16744.186 231.43523
48000 296.63454 1.0000059 -16698.147 -750.76937
49000 303.13461 1.0000059 -16703.513 -369.22278
50000 292.71464 1.0000059 -16686.245 433.57969
51000 301.67637 1.0000059 -16617.251 332.72543
52000 303.97802 1.0000059 -16722.127 681.26606
53000 297.76704 1.0000059 -16750.351 -337.23038
54000 303.08515 1.0000059 -16722.017 -236.79643
55000 298.90843 1.0000059 -16777.384 193.38446
56000 293.29691 1.0000059 -16776.642 -55.164427
57000 302.65835 1.0000059 -16727.362 67.004336
58000 296.20336 1.0000059 -16742.291 -63.663928
59000 301.9545 1.0000059 -16698.477 -116.76893
60000 304.42752 1.0000059 -16597.112 -537.5654

The problem is more a problem of understanding MD and statistical mechanics than a LAMMPS problem.
Besides, pressure and temperature fluctuations have been discussed many times before, thus please search through the old posts for such discussions and enlightenment.

Dear Sir,
Thank you for your response. Pressure difference was expected, but why density was not changing? Is there any problem with my ensemble selection?

You should discuss with your adviser. Not a LAMMPS issue. The simulation does what you asked for.

You don’t select an ensemble, you select an integrator. And again, this has been discussed many times before. Please have the courtesy to check the (extensive) collection of old topics reaching back to 2005.

Thank you for your suggestion.