I am trying to calculate surface tension for Argon liquid but I still have fluctuation in the surface tension results it is unstable even after 2 million timestep this is the script :

# Argon Surface tension

units real

atom_style hybrid atomic sphere

lattice fcc 5.26

boundary p p p

region box block 0 1 0 1 0 1

create_box 1 box

create_atoms 1 box

# forcefield

mass 1 39.948 # Ar

pair_style lj/cut 8.7675

pair_coeff 1 1 0.233 3.507 8.7675

replicate 10 10 10

variable xLo equal -0.5*lx
variable xHi equal 0.5*lx

variable yLo equal -0.5

*ly*

variable yHi equal 0.5ly

variable yHi equal 0.5

variable zLo equal -0.5

*lz*

variable zHi equal 0.5lz

variable zHi equal 0.5

# Recenter system coords about the origin (0,0,0):

displace_atoms all move {xLo} {yLo} ${zLo} units box

# Adjust box dimensions:

change_box all x final {xLo} {xHi} units box

change_box all y final {yLo} {yHi} units box

change_box all z final {zLo} {zHi} units box

# Increase z edge to yield “Vacuum|Argon|Vacuum” system:

variable zLoNew equal -2.5*lz
variable zHiNew equal 2.5*lz

change_box all z final {zLoNew} {zHiNew} units box

variable Nblock equal 10000

variable Nrun equal 50*{Nblock}
variable Nr equal {Nblock}/100

variable A_in_m equal 1e-10 # Angstrom in meter

variable atm_in_Pa equal 101325 # note: 1 Pa = 1 N/m^2

variable N_in_mN equal 1e3 # Newton in milliNewton

variable convFac equal {A_in_m}*{atm_in_Pa}*${N_in_mN}

variable Text equal 100

group Argon type 1

velocity all create ${Text} 1234

neighbor 2.0 bin

neigh_modify every 1 delay 0 check yes

fix 1 all nve/sphere

fix thermostat all langevin {Text} {Text} 100 1234 zero yes

fix removeMomentum all momentum 1 linear 1 1 1

compute T all temp

fix TAve all ave/time 10 100 1000 c_T

variable P equal press

fix PAve all ave/time 10 100 1000 v_P

variable xPress equal c_thermo_press[1]

variable yPress equal c_thermo_press[2]

variable zPress equal c_thermo_press[3]

# Evaluate and average surface tension in mN/m:

variable st equal 0.5*lz*(v_zPress-0.5*(v_xPress+v_yPress))*${convFac}

fix st all ave/time 10 100 1000 v_st

variable xyArea equal lx*ly

thermo_style custom step temp f_TAve press f_PAve etotal v_xyArea lz f_st

thermo_modify norm yes flush yes

thermo 1000

dump trj all custom 1000 Argon.trj id type x y z

dump_modify trj sort id pad 6

timestep 1

run ${Nrun}