# Define variables
# Starting temperature:
variable T equal 1
# Number of molecules:
variable N equal 500
# Density:
variable rho equal 0.01
# Length of sides of cubic box
variable L equal (${N}/${rho})^(1/3)
# How much heat to add per timestep
variable eFlux equal 1
# Spring bond (0) or rigid bond (1)
variable rigidbond equal 0

# Use Lennard-Jones units, periodic boundaries and atomic model
units lj
boundary p p p
atom_style full

# Define the space and put atoms in
dimension 3
region simbox block 0 $L 0 $L 0 $L
create_box 1 simbox bond/types 1 extra/bond/per/atom 1

# Define properties of the molecules
bond_style class2
bond_coeff 1 0.7 35.0 -40.0 59.0
molecule LJdiat diat.molecule
create_atoms 0 random $N 12345 simbox mol LJdiat 1234
mass 1 1.0
velocity all create $T 87287
pair_style lj/class2 2.5
pair_coeff 1 1 1.0 1.0 2.5

#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"
#timestep 2

# Speed up simulation with neighbour lists
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

# Define thermodynamic ensemble for thermalization
fix nvt all nvt temp $T $T 1

#Thermalization run
run             100000

# Define thermodynamic ensemble for long run
unfix nvt
fix nve all nve

# Start heat addition and run to stabilize
fix heat all heat 1 ${eFlux} region simbox
run 1000

#Define frequency of 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
