# used for offical lammps, developed by deyun
#========================================================================================
variable skindistance equal 0.0001 #say 10% of max particles
variable interval_fix equal 100
variable interval_dump equal 100
variable pstart equal 0 # use 0 
variable pend equal 30000 
variable pdamp equal 1000
variable rate equal -1.04e-2 # real velocity unit should be m
#========================================================================================
variable clumpratio equal 0.199327526 # be careful, this is estimated from matlab, should be case by case
#========================================================================================
units		si # (This means m, kg, s for base units.)
boundary	p p p 
atom_style	bpm/sphere
dimension	3
newton		on off
echo            both
comm_modify     vel yes 
processors      * * * 
#========================================================================================
read_data       rubber_sand15_12.data
#read_restart   restart_fileverify.100
#========================================================================================
neighbor	${skindistance} bin 
group rubber type 1
group sand   type 2
#neigh_modify exclude molecule/intra rubber delay 0 every 1 one 200 page 2000
neigh_modify delay 0 every 1 one 200 page 2000
bond_style      bpm/rotational
bond_coeff      1 1e4 1e4 5e5 5e5 1.0e12 1.0e12 1.0e12 1.0e12 1 1 0.0004 0.0004
special_bonds   lj 0 1 1 coul 1 1 1
#========================================================================================
pair_style granular 
pair_coeff 2 2 hertz/material 7.0e10 0.1 0.2 tangential mindlin_rescale/force NULL 1.0 0.001 
pair_coeff 1 1 hertz/material 7e6 0.5 0.5 tangential mindlin_rescale/force NULL 1.0 1.0    



timestep	1e-8
run_style	verlet
#========================================================================================
compute     2 all pair/local dist dx dy dz force fx fy fz p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 cutoff radius
compute     3 all property/local patom1 patom2 cutoff radius
compute     4 all property/atom id radius x y z vx vy vz fx fy fz 
compute     13 sand property/atom id radius 
compute     14 rubber property/atom id radius
variable    sandvol atom "4.0*PI* c_13[2] * c_13[2] * c_13[2] /3.0"
variable    rubbervol atom "4.0*PI* c_14[2] * c_14[2] * c_14[2] /3.0"
compute     svol sand reduce sum v_sandvol
compute     rvol rubber reduce sum v_rubbervol 
compute     5 all contact/atom 
compute     connnew all reduce sum c_5  
compute     6 all stress/atom NULL
compute     7 all pressure NULL pair 
compute 	aveconn all reduce ave c_5 
compute     10 all nbond/atom
compute     11 all bond/local dist force fx fy fz dx dy dz
compute     12 all property/local batom1 batom2 cutoff radius
#========================================================================================
variable 	step equal step
variable	xlength equal lx
variable	ylength equal ly
variable	zlength equal lz
variable    xxstress equal c_7[1]
variable    yystress equal c_7[2]
variable    zzstress equal c_7[3]
variable    xystress equal c_7[4]
variable    xzstress equal c_7[5]
variable    yzstress equal c_7[6]
variable    coordNum_new  equal  c_aveconn
variable    voidratio equal "(lz * lx * ly - (c_svol + c_rvol * v_clumpratio))/(c_svol + c_rvol * v_clumpratio)"
#========================================================================================
fix         9 all press/berendsen x ${pstart} ${pend} ${pdamp} y ${pstart} ${pend} ${pdamp} z ${pstart} ${pend} ${pdamp}
#fix        10 all deform 10 z vel ${rate} units box 
fix	        3 all nve/bpm/sphere
#========================================================================================
fix 		2 all print ${interval_fix} "${step} ${xxstress} ${yystress} ${zzstress} ${xystress} ${xzstress} ${yzstress} ${xlength} ${ylength} ${zlength}" file verify_mean_stresses_iso.txt screen no
fix 		6 all print ${interval_fix} "${step} ${voidratio} ${coordNum_new}" file verify_voidratio.txt screen no
#========================================================================================
dump	     1 all local  ${interval_dump} verify_contact/dump*.contact c_3[1] c_3[2] c_2[1] c_2[2] c_2[3] c_2[4] c_2[5] c_2[6] c_2[7] c_2[8] c_2[9] c_2[10] c_2[11] c_2[12] c_2[13] c_2[14] c_2[15] c_2[16] c_2[17] c_2[18] c_2[19] c_2[20]     
dump         2 all custom ${interval_dump} verify_atom/dump*.atom id type c_5 c_10 x y z diameter c_6[1] c_6[2] c_6[3] c_6[4] c_6[5] c_6[6] c_4[6] c_4[7] c_4[8] c_4[9] c_4[10] c_4[11]
dump         6 all local  ${interval_dump} verify_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] 
#========================================================================================
restart     100 restart_fileverify
#=========================================================================================
run		    999999999
#=========================================================================================