### DEFINE SIMULATION BOX GEOMETRY ###
variable d equal 0.2
lattice sc ${d}
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0
variable p equal 2.75
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
unfix nvt
fix nve all nve
reset_timestep 0

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp atoms vol
variable volume equal vol
variable temp equal temp
variable N2 equal atoms*atoms
variable temp2 equal temp*temp
variable energy equal etotal
variable energy2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_energy v_temp v_energy2 v_temp2
run 100000

variable avenergy equal f_aves[1]
variable avtemp equal f_aves[2]
variable avenergy2 equal f_aves[3]
variable avtemp2 equal f_aves[4]
variable cv equal ${N2}*(${avenergy2}-${avenergy}*${avenergy})/(${avtemp2})


print "Averages"
print "--------"
print "Heat_Capacity: ${cv}"
print "Heat_Capacity/Vol: ${cv}/${volume}"
print "Volume: ${volume}"
print "Energy: ${avenergy}"
print "Energy^2: ${avenergy2}"
print "Temperature: ${avtemp}"
print "Temperature^2: ${avtemp2}"