Statistical Physics Simulation
### DEFINE SIMULATION BOX GEOMETRY ###
variable density equal 0.8
lattice sc ${density}
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 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
reset_timestep 0
unfix nvt
fix nve all nve
### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp density atoms vol
variable volume equal vol
variable dens equal density
variable dens2 equal density*density
variable energy equal etotal
variable energy2 equal etotal*etotal
variable atoms2 equal atoms*atoms
variable temp equal temp
variable temp2 equal temp*temp
fix aves all ave/time 100 1000 100000 v_dens v_temp v_dens2 v_temp2 v_energy v_energy2
run 100000
variable cp equal ${atoms2}*(f_aves[6]-f_aves[5]*f_aves[5])/f_aves[4]
variable avedens equal f_aves[1]
variable avetemp equal f_aves[2]
variable errdens equal sqrt(abs(f_aves[4]-f_aves[1]*f_aves[1]))
variable errtemp equal sqrt(abs(f_aves[5]-f_aves[2]*f_aves[2]))
print "Averages"
print "--------"
print "Heat Capacity: ${cp}"
print "Volume: ${volume}"
print "Density: ${avedens}"
print "Stderr: ${errdens}"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
