#
#  define units
#
units       lj

#
#  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
#  Then I can determine the size of a cube by 
#  size = (N/rho)^(1/3)
#

variable l      equal 10.0793684
region      boxid block 0.0       ${l}      0.0       ${l}      0.0       ${l}
create_box  2 boxid

#
#  Split the simulation into two regions
#
variable     half_l equal ${l}/2.0
region      boxa  block 0.0       ${half_l} 0.0 ${l} 0.0 ${l} 
region      boxb  block ${half_l} ${l}      0.0 ${l} 0.0 ${l} 

#
#  specify initial positions of atoms
#  sc = simple cubic
#  0.5 = density in lj units
#
lattice     sc 0.50
#  place atoms of type 1 in boxa
#  place atoms of type 2 in boxb
create_atoms    1 region boxa
create_atoms    2 region boxb

#
#   define mass
#
#   mass of atom type 1
mass        1 1.0
mass        2 1.0

#
#  specify initial velocity of atoms
#  group = all
#  reduced temperature is T = 1.0 = lj-eps/kb 
#  seed for random number generator
#  distribution is gaussian (e.g. Maxwell-Boltzmann)
#
velocity    all create 1.0 87287 dist gaussian

#
#  specify interaction potential
#  pairwise interaction via the Lennard-Jones potential with a cut-off at 2.5 lj-sigma
#
pair_style  lj/cut 2.5
#  specify parameters between atoms of type 1 with an atom of type 1
#  epsilon = 1.0, sigma = 1.0, cutoff = 2.5
pair_coeff  1 1 1.0 1.0 2.5
pair_coeff  2 2 1.2 0.9 2.5

#
# add long-range tail correctio

#
pair_modify tail yes

#
#  specify parameters for neighbor list 
#  rnbr = rcut + 0.3
#
neighbor    0.3 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
thermo 100
thermo_modify norm no

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

timestep 0.001

run 1000
#
#  stop fix with given fixid
#  fixid = 1
#
unfix 1
#
#  specify ensemble
#  fixid = 2
#  atoms = all
#  ensemble = nvt
#  temp = temperature
#  initial temperature = 1.0
#  final temperature = 1.0
#  thermostat controller gain = 0.1 (units of time, bigger is less tight control)
#
fix     2 all nvt temp 1.0 1.0 0.1 

run     1000

#
#   save configurations
#   dumpid = 1
#   filename = output.xyz
#
dump    1       all xyz 100 output.xyz
dump_modify 1 element Ar Ne


run     1000