# Define variables
# Temperature:
variable T equal 300
# Length of sides of cubic box
variable L equal 200
# Spring bond (0) or rigid bond (1)
variable rigidbond equal 0

units real
#    mass = grams/mole
#    distance = Angstroms
#    time = femtoseconds
#    Energy = Kcal/mole
#    velocity = Angstroms/femtosecond
#    force = Kcal/mole-Angstrom
#    temperature = Kelvin
#    pressure = atmospheres
#    density = gram/cm^dim

atom_style full
region box block 0 $L 0 $L 0 $L
create_box 1 box bond/types 1 extra/bond/per/atom 1
pair_style lj/class2 10.0
pair_coeff * * 0.05980 3.8008
mass * 14.0
bond_style class2
bond_coeff 1 1.0977 1651.3730 -4069.3178 5984.9629
molecule N2 N2.molecule
create_atoms 0 random 500 12345 box mol N2 1234
#make sure no explosion happens
minimize 1e-4 1e-4 100 100

#Make the N-N bond rigid if rigidbond=1
if "${rigidbond} == 1" then &
"fix shake all shake 0.0001 20 0 b 1"

neigh_modify every 1 delay 0 check yes

velocity all create $T 3928459
timestep 2

# Equilibrate the system at constant T
fix nvt all nvt temp $T $T 1000.0
thermo          100
run             1000

#Start heat addition and run to stabilize
unfix nvt
fix nve all nve
variable eFlux equal 0.001
fix heat all heat 1 ${eFlux} region box
run 1000

#Define new thermodynamic output and trajectory dump and start production run
thermo_style custom step temp epair etotal press 
thermo 100
dump 1 all custom 100 dump.lammpstrj id x y z vx vy vz
run 100000