# used for offical lammps
#========================================================================================
variable skindistance equal 0.0001 #say 10% of max particles
variable interval_fix equal 1000
variable interval_dump equal 1000000
variable EWall equal  7e7
#========================================================================================
units		si # (This means m, kg, s for base units.)
boundary	p p f  # (x and y are p; z is the rigid wall)
atom_style	sphere
dimension	3
newton		on 
echo            both
comm_modify     vel yes 
processors      * * * 
#========================================================================================
read_data       single.data
#========================================================================================
neighbor	${skindistance} bin 
group rubber type 1
group sand   type 2
neigh_modify  delay 0 every 1 one 2000 page 20000
#========================================================================================
pair_style granular 
pair_coeff * * hertz/material 7e7 1.0 0.5 tangential mindlin_rescale/force NULL 1.0 1.0
#pair_coeff * * hooke/piecewise 1e4 1e7 0 0.0001 tangential linear_nohistory 1.0 0.4 damping mass_velocity
#tangential mindlin_rescale/force NULL 1.0 1.0

timestep	1e-9
run_style	verlet
#========================================================================================
compute     4 all property/atom id radius x y z vx vy vz fx fy fz
compute     12 all property/local patom1 patom2 cutoff radius
#=====================================================
fix 3 all wall/gran granular hertz/material ${EWall} 1 0.5 tangential mindlin_rescale/force NULL 1.0 1.0 zplane -0.0100000001 NULL wiggle z 0.1 0.1 contacts
fix 4 all wall/gran granular hertz/material ${EWall} 1 0.5 tangential mindlin_rescale/force NULL 1.0 1.0 zplane NULL 0.01 contacts
fix 5 all nve/sphere
#======================================================
variable  step equal step
variable  bottomforcez equal f_3[4]
variable  topforcez equal f_4[4]
variable  bottomz equal f_3[7]
variable  topz equal f_4[7]
compute   bottomwallforce all reduce sum f_3[4]
compute   topwallforce all reduce sum f_4[4]
compute   bottomtag all reduce sum f_3[1]
compute   toptag    all reduce sum f_4[1]
compute   bottompos all reduce sum f_3[7]
compute   toppos    all reduce sum f_4[7]
variable  bottom equal c_bottomwallforce
variable  top    equal c_topwallforce
variable  bottomwall equal c_bottompos
variable  topwall    equal c_toppos
variable  bottomtagg equal c_bottomtag  
variable  toptagg    equal c_toptag

#variable  bottomwallpos equal "c_bottompos / c_bottomtag"
#variable  topwallpos equal "c_toppos / c_toptag"
#========================================================================================
fix         2 all print ${interval_fix} "${step} ${bottom} ${top} ${bottomwall} ${topwall}  ${bottomtagg}  ${toptagg} " file pd3_force.txt screen no

#dump        2 all custom ${interval_dump} pd3_atom/dump*.atom id x y z diameter c_4[9] c_4[10] c_4[11]
#dump        6 all local  ${interval_dump} pd3_bond/dump*.bond c_12[1]  c_12[2] c_11[1] c_11[2] c_11[3] c_11[4] c_11[5] c_11[6] c_11[7] c_11[8] 
#dump        7 all custom ${interval_dump} pd3_contact/dump*.contact f_3[1] f_3[2] f_3[3] f_3[4] f_3[5] f_3[6] f_3[7] f_3[8] f_4[1] f_4[2] f_4[3] f_4[4] f_4[5] f_4[6] f_4[7] f_4[8]
#=========================================================================================
run		    500000000
#=========================================================================================