# Numerical difference calculation
# of error in forces and virial stress

# adjustable parameters

variable    nsteps index 500    # length of run
variable    nthermo index 10    # thermo output interval 
variable    ndump index 500     # dump output interval
variable    nlat index 3        # size of box
variable    fdelta index 1.0e-4 # displacement size
variable    vdelta index 1.0e-6 # strain size 
variable    temp index 10.0     # temperature

units  	    metal
atom_style  atomic

atom_modify	map yes
lattice 	fcc 5.358000
region 		box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
create_box  	1 box
create_atoms 	1 box
mass 		1 39.903

velocity     all create ${temp} 2357 mom yes dist gaussian

pair_style   lj/cubic
pair_coeff   * * 0.0102701 3.42

neighbor     0.0 bin
neigh_modify every 1 delay 0 check no

timestep     0.001
fix	     nve all nve

# define numerical force calculation

fix	     numforce all numdiff ${nthermo} ${fdelta}
variable     ferrx atom f_numforce[1]-fx
variable     ferry atom f_numforce[2]-fy
variable     ferrz atom f_numforce[3]-fz
variable     ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2
compute	     faverrsq all reduce ave v_ferrsq
variable     fsq atom fx^2+fy^2+fz^2
compute      favsq all reduce ave v_fsq
variable     frelerr equal sqrt(c_faverrsq/c_favsq)
dump errors  all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz

# define numerical virial stress tensor calculation

compute 	myvirial all pressure NULL virial
fix 		numvirial all numdiff/virial ${nthermo} ${vdelta}
variable 	errxx equal f_numvirial[1]-c_myvirial[1]
variable 	erryy equal f_numvirial[2]-c_myvirial[2]
variable 	errzz equal f_numvirial[3]-c_myvirial[3]
variable 	erryz equal f_numvirial[4]-c_myvirial[6]
variable 	errxz equal f_numvirial[5]-c_myvirial[5]
variable 	errxy equal f_numvirial[6]-c_myvirial[4]
variable 	verrsq equal "v_errxx^2 + &
                              v_erryy^2 + &
                              v_errzz^2 + &
                              v_erryz^2 + &
                              v_errxz^2 + &
                              v_errxy^2"
variable 	vsq equal "c_myvirial[1]^2 + &
                           c_myvirial[3]^2 + &
                           c_myvirial[3]^2 + &
		           c_myvirial[4]^2 + &
                           c_myvirial[5]^2 + &
                           c_myvirial[6]^2"
variable     	vrelerr equal sqrt(v_verrsq/v_vsq)

thermo_style 	custom step temp pe etotal press v_frelerr v_vrelerr
thermo 		${nthermo}
run 		${nsteps}