### DEFINE SIMULATION BOX GEOMETRY ### variable density equal 0.2 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 energy equal etotal variable energy2 equal etotal*etotal variable atoms2 equal atoms*atoms variable dens2 equal density*density 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 heatcap 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(f_aves[4]-f_aves[1]*f_aves[1]) variable errtemp equal sqrt(f_aves[5]-f_aves[2]*f_aves[2]) print "Averages" print "--------" print "Heat Capacity: ${heatcap}" print "Volume: ${volume}" print "Density: ${avedens}" print "Stderr: ${errdens}" print "Temperature: ${avetemp}" print "Stderr: ${errtemp}"