Triaxial Test

I would be really happy that any of you share a triaxial code with me as I am beginner in this field.
I have already wrote this code :

#####################################################################Triaxial Test######################################################################

#################################################################### Granular model ####################################################################

atom_style granular

atom_modify map array

boundary m m m

newton off

communicate single vel yes

units si

#################################################################### Simulation box ####################################################################

region reg block -0.01 0.01 -0.01 0.01 0.0 0.05 units box

create_box 2 reg

neighbor 0.002 bin

neigh_modify delay 0

hard_particles yes

################################################### Material properties required for new pair styles ####################################################

fix m1 all property/global youngsModulus peratomtype 2.6e7 6.9e13
fix m2 all property/global poissonsRatio peratomtype 0.1 0.1
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.7 0.6 0.6 0.8
fix m4 all property/global coefficientFriction peratomtypepair 2 0.24 0.24 0.24 0.24
fix m5 all property/global RollingFriction peratomtypepair 2 0.1 0.1 0.1 0.1
fix m6 all property/global SlidingFriction peratomtypepair 2 0.166 0.166 0.166 0.166
#fix m7 all property/global characteristicVelocity scalar 2.0

###################################################################### pair style #######################################################################

pair_style gran model hertz tangential history

pair_coeff * *

timestep 0.00001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

######################################################################## boudary ########################################################################

fix xwalls1 all wall/gran model hertz tangential history primitive type 2 xplane -0.005
fix xwalls2 all wall/gran model hertz tangential history primitive type 2 xplane 0.005
fix ywalls1 all wall/gran model hertz tangential history primitive type 2 yplane -0.005
fix ywalls2 all wall/gran model hertz tangential history primitive type 2 yplane 0.005
fix zwalls1 all wall/gran model hertz tangential history primitive type 2 zplane 0.00

################################################################## region of packing ####################################################################

region bc block -0.00498 0.00498 -0.00498 0.00498 0.0 0.023 units box

################################################################### particle packing ####################################################################

fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 1631 radius constant 0.0005
fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 1631 radius constant 0.0006
fix pdd1 all particledistribution/discrete/numberbased 86028157 2 pts1 0.8 pts2 0.2

fix ins all insert/pack seed 32452843 distributiontemplate pdd1 vel constant 0.0 0.0 -0.35 &
insert_every once overlapcheck yes all_in yes particles_in_region 800 region bc

run 20000

#################################################################### Wall Properties ####################################################################

fix cad3 all mesh/surface/stress/servo file cad3.stl type 1 scale 0.001 stress on com 0.0 -0.005 0.005 ctrlPV force axis 0.0 -1.0 0.0 target_val -1000 vel_max 0.025
fix cad4 all mesh/surface/stress/servo file cad4.stl type 1 scale 0.001 stress on com 0.0 0.005 0.005 ctrlPV force axis 0.0 1.0 0.0 target_val -1000 vel_max 0.025
fix cad5 all mesh/surface/stress/servo file cad5.stl type 1 scale 0.001 stress on com 0.005 0.0 0.005 ctrlPV force axis -1.0 0.0 0.0 target_val -1000 vel_max 0.025
fix cad6 all mesh/surface/stress/servo file cad6.stl type 1 scale 0.001 stress on com -0.005 0.0 0.005 ctrlPV force axis 1.0 0.0 0.0 target_val -1000 vel_max 0.025

fix geometry1 all wall/gran model hertz tangential history mesh n_meshes 4 meshes cad3 cad4 cad5 cad6

####################################################### apply nve integration to all particles ##########################################################

fix integr all nve/sphere

################################################### output settings, include total thermal energy #######################################################

compute peratom all stress/atom

compute rke all erotate/sphere

thermo_style custom step atoms ke c_rke vol

thermo 1000

thermo_modify lost ignore norm no

fix zwalls2 all wall/gran model hertz tangential history primitive type 2 zplane 0.03

run 20000

unfix ins
unfix geometry1

fix cad1 all mesh/surface/stress/servo file cad1.stl type 1 scale 0.001 stress on com 0.0 0.0 0.0 ctrlPV force axis 0.0 0.0 1.0 target_val 1000 vel_max 0.025
fix cad2 all mesh/surface/stress/servo file cad2.stl type 1 scale 0.001 stress on com 0.0 0.0 0.01 ctrlPV force axis 0.0 0.0 -1.0 target_val 1000 vel_max 0.025
fix geometry2 all wall/gran model hertz tangential history mesh n_meshes 6 meshes cad1 cad2 cad3 cad4 cad5 cad6

run 20000

dump dmp all custom/vtk 800 post/newModels_*.vtk c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
dump dmpBOX all mesh/stl 800 post/Disc*.stl

######################################################################### run ############################################################################

run 200000

but I want to make a comparison.

Thank you !