# used for offical lammps, developed by deyun
#========================================================================================
variable skindistance equal 0.0001 #say 10% of max particles
variable interval_fix equal 10000
variable interval_dump equal 100000
variable pstart equal 1000000
variable pend equal 1000000
variable pdamp equal 10
variable rate equal -6e-4
#========================================================================================
variable clumpratio equal 1.0 # 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	sphere
dimension	3
newton		on 
echo            both
comm_modify     vel yes 
processors      * * * 
#========================================================================================
#read_data       Particle_coarse.lj
read_restart   restart_filetry1k.52000000
#========================================================================================
neighbor	${skindistance} bin 
group layer1 type 1
group layer2 type 2 
neigh_modify  delay 0 every 1 one 2000 page 20000
#========================================================================================
pair_style granular 
pair_coeff * * hertz/piecewise 3.5e8 7e10 0.25 0.25 1.0 0.005 tangential mindlin_rescale/force 3.2e8 6.4e10 1.0 0.001 0.001 0.005
#pair_coeff * * hertz/material 7e10 0.1 0.2 tangential mindlin_rescale/force NULL 1.0 0.25
#pair_coeff 2 2 hertz/material 7e10 0.1 0.2 tangential mindlin_rescale/force NULL 1.0 0.25

timestep	1e-7
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 layer1 property/atom id radius
variable    layer1vol atom "4.0*PI* c_13[2] * c_13[2] * c_13[2] /3.0"
compute     lay1vol layer1 reduce sum v_layer1vol
compute     14 layer2 property/atom id radius
variable    layer2vol atom "4.0*PI* c_14[2] * c_14[2] * c_14[2] /3.0"
compute     lay2vol layer2 reduce sum v_layer2vol

compute     5 all contact/atom 
compute     connnew all reduce sum c_5  
compute     6 all stress/atom NULL pair
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_lay1vol + c_lay2vol))/(c_lay1vol + c_lay2vol)"
variable    strain equal " (0.0403569 - lz) / 0.0403569 "
#========================================================================================
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/sphere
#========================================================================================
fix 		2 all print ${interval_fix} "${step} ${xxstress} ${yystress} ${zzstress} ${xystress} ${xzstress} ${yzstress} ${xlength} ${ylength} ${zlength}" file try1k_mean_stresses_iso.txt screen no
fix 		6 all print ${interval_fix} "${step} ${voidratio} ${coordNum_new} ${strain}" file try1k_voidratio.txt screen no
#========================================================================================
dump	    1 all local  ${interval_dump} try1k_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} try1k_atom/dump*.atom id type c_5 c_5 xu yu zu 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]
#========================================================================================
restart     1000000 restart_filetry1k
#=========================================================================================
run		    999999999
#run         100000
#=========================================================================================