# define unit set and class of atomic model
units metal
atom_style atomic

# periodic boundary conditions along all three dimensions
boundary p p p

# create a supercell with cubic lattice (fcc, bcc, sc, or diamond)
# using 10*10*10 conventional (orthogonal) unit cells
lattice       bcc  4.42748749256
region        simbox  block 0 10 0 10 0 10 units lattice
create_box    1 simbox
create_atoms  1 box

# specify which KIM Model to use, letting LAMMPS compute the virial/pressure
pair_style    kim LAMMPSvirial Pair_Morse_Shifted_GirifalcoWeizer_LowCutoff_Na__MO_707981543254_001
pair_coeff    * * Na 
mass          1 22.98977

# create initial velocities consistent with the chosen temperature
velocity      all create 333.15 17 mom yes rot yes

# set the time step (in picoseconds)
timestep      0.001

# set NPT ensemble for all atoms
fix           ensemble all npt temp 333.15 333.15 0.1 iso 0.0 0.0 1

# compute the time averages of pressure, temperature, and volume, respectively
# ignore the first 5000 timesteps
variable      myTemp  equal temp
variable      myPress equal press 
variable      myVol   equal vol
fix           avgmyTemp  all ave/time 5 20 100 v_myTemp  ave running start 5000
fix           avgmyPress all ave/time 5 20 100 v_myPress ave running start 5000
fix           avgmyVol   all ave/time 5 20 100 v_myVol   ave running start 5000

variable      T equal f_avgmyTemp  # cannot use f_avgmyTemp directly in `if' statement
variable      P equal f_avgmyPress 
variable      V equal f_avgmyVol 

# set error bounds for temperature and pressure
variable      T_low equal "333.15 - 0.2"
variable      T_up  equal "333.15 + 0.2"
variable      P_low equal "0.0 - 0.2"
variable      P_up  equal "0.0 + 0.2"

# print to logfile every 1000 timesteps 
thermo        1000 

# Run a simulation for at most 2000*1000 timesteps. At each 1000th time step, check 
# whether the temperature and volume have converged. If yes, break.  
label top 
variable a loop 2000
run 1000
if "${T}>${T_low} && ${T}<${T_up} && ${P}>${P_low} && ${P}<${P_up}" then "jump SELF break"
print "flag: Temp = ${T}, Press = ${P}"
next a
jump SELF top 
label break

# Write final averaged volume to file if temperature and volume have converged; otherwise wirte a 
# flag to indicate non-convergence.
variable      myStep equal step 
if "${myStep} < 2000000" then "print '${V}'  file output/vol_T333.15.out"  &
else  "print 'not_converged'  file output/vol_T333.15.out"

print "ALL DONE"