#  define units
units       real

#  specify periodic boundary conditions
boundary p p p

#  define atom_style
#  full covers everything
atom_style  full 

#  define simulation volume 
#  If I want N = 512 atoms 
#  and I want a density of rho = 0.5 atoms/lj-sigma^3
# in real units density is rho = 0.012665393 atoms/Angstrom^3
#  Then I can determine the size of a cube by 
#  side = (N/rho)^(1/3)
variable side      equal 29.34340048
region      boxid block 0.0 ${side} 0.0 ${side} 0.0 ${side}
create_box  1 boxid

#  specify initial positions of atoms
#  sc = simple cubic
#  set lattice parameter to match density above
#  target number of atoms N = 512
#  lattice parameter = N^(1/3)*side = 512^(1/3)*29.34340048
#  lattice parameter = 3.6679 Angstroms + epsilon
lattice     sc  3.668

#  place atoms of type 1 in boxid
create_atoms    1 box

#   define mass of atom type 1
mass        1 39.948

#  specify initial velocity of atoms
#  group = all
#  temperature is T = 119.8 K 
#  seed for random number generator
#  distribution is gaussian (e.g. Maxwell-Boltzmann)
velocity    all create 119.8 87287 dist gaussian

#  specify interaction potential
#  pairwise interaction via the Lennard-Jones potential with a cut-off at 2.5 lj-sigma (8.5125 Angstroms)
pair_style  lj/cut 8.5125

#  specify parameters between atoms of type 1 with an atom of type 1
#  epsilon = 0.23806, sigma = 3.405, cutoff = 8.5125
pair_coeff  1 1 0.23806 3.405 8.5125

# add long-range tail correction
pair_modify tail yes

#  specify parameters for neighbor list 
#  rnbr = rcut + 0.3 (1.0215 Angstrom)
neighbor    1.0215 bin

#  specify thermodynamic properties to be output
#  pe = potential energy
#  ke = kinetic energy
#  etotal = pe + ke
#  temp = temperature
#  press = pressure
#  density = number density
#  output every thousand steps
#  norm = normalize by # of atoms (yes or no)
thermo_style custom step pe ke etotal temp press density

# report instantaneous thermo values every 100 steps
thermo 100

# normalize thermo properties by number of atoms (yes or no)
thermo_modify norm no

#   save configurations
#   dumpid = 1
#   all atoms
#   atomic symbol is Ar
#   save positions every 100 steps
#   filename = output.xyz
#
dump    1       all xyz 100 output.xyz
dump_modify 1 element Ar


#
#  FIRST EQUILIBRATION:  MELT Simple Cubic Crystal in NVE
#

#  specify ensemble
#  fixid = 1
#  atoms = all
#  ensemble = nve or nvt
fix     1 all nve

#  define time step (fs)
timestep 2.0

# run 5000 steps in the NVE ensemble
# (this equilibrates positions)
run 5000

#  stop fix with given fixid
#  fixid = 1
unfix 1

#
#  SECOND EQUILIBRATION:  Equilibrate to set temperature in NVT
#

#  specify ensemble
#  fixid = 2
#  atoms = all
#  ensemble = nvt
#  temp = temperature
#  initial temperature = 119.8
#  final temperature = 119.8
#  thermostat controller gain = timestep*100 = 200.0 (units of time, bigger is less tight control)
fix     2 all nvt temp 119.8 119.8 200.0 

# run 5000 steps in the NVT ensemble
# (this equilibrates thermostat)
run     5000

#
#  DATA PRODUCTION:  run at set temperature in NVT
#

# run 50,000 more steps in the NVT ensemble
# (this is data production) 
run     50000