I am trying to monitor the average stress on each atom as a graphene sheet gets pulled apart. I planned on doing this by dividing the force on each atom by the volume of the atoms. However the voronoi compute used to get the volume of each atom brings back an error. Is there a work arround for adding computes that might require additional packages?
Below is the code ran:
# ------------------------ INITIALIZATION ----------------------------
units real #real units, makes time femtoseconds
dimension 3
boundary p p f
atom_style charge
# ----------------------- ATOM DEFINITION ----------------------------
region whole block 0 21.4 0 21.4 -5 5 #expanded box due to greater interatomic spacing
create_box 1 whole
read_data SLGSarmcharge.txt add merge shift 13.5 13 0
# ------------------------ FORCE FIELDS ------------------------------
pair_style reax/c NULL checkqeq yes
pair_coeff * * ffield.reax.mattsson C
neighbor 0.3 bin
neigh_modify delay 5
# ------------------------- SETTINGS ---------------------------------
compute peratom all pe/atom
compute sv all stress/atom NULL #maybe add virial
compute vol all voronoi/atom
log log.txt
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
######################################
# EQUILIBRATION
reset_timestep 0
timestep 0.001
velocity all create 300 12465 mom yes rot no
fix 2 all nvt temp 300 300 $(100*dt)
min_style sd
minimize 1.0e-5 1.0e-6 1000 10000
# Set thermo output
thermo 1000
thermo_style custom step lx ly lz press pxx pyy pzz pe temp
unfix 2
# Store final cell length for strain calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"
######################################
# DEFORMATION
reset_timestep 0
fix 2 all nvt temp 300 300 $(110*dt)
variable srate equal 1.7e10
variable srate1 equal "v_srate / 1.0e12"
fix 3 all deform 1 x erate ${srate1} units box remap x
# Output strain and stress info to file
# for units real, pressure is in [atm] = 0.000101325 [GPa]
# sxx is in GPa, vol is in nm^3
variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "v_strain"
variable vol equal "lx*ly*lz"
variable s atom "c_sv[1]/c_vol[1]" #converts units press*vol to press
compute sxx all reduce ave v_s #ave the per atom stresses
variable sxx equal "c_sxx*0.000101325"
fix def1 all print 100 "${p1} ${sxx} ${vol}" file armchair_tensile.def1.txt screen no title "strain stress_x (Gpa) vol (nm^3)"
# Use cfg for Ovito
dump 1 all cfg 500 armchair.tensile_*.cfg mass type xs ys zs c_peratom fx fy fz
#Display thermo
thermo 1000
thermo_style custom step v_strain v_sxx temp ke pe press
run 60000
unfix 2
unfix 3
######################################
# SIMULATION DONE
print "All done"