# modified by LIU for reference 
#  define input parameters 
variable skindistance equal 0.00005 # !!!not sure, say 10% of max particles

variable shearmodulus equal 29.17e9
variable poissonratio equal 0.2
variable frictioncoeff equal 0.1


variable maxrate equal 1.0
variable targetstress equal 100000 #pa
variable localDamp equal 0.0

variable     interval_fix equal 10000
variable     interval_dump equal 100000

# define the contacts used for rose diagram



# Initialisation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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   * * * # new thing (!!!not identified in XIN version)

# Create particles ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#read_data     Particle_ToyouraSandCoordinates_e1.1.lj
read_restart    restart_filerandom.384200000


# Pre-fix settings ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
neighbor	${skindistance} bin 
neigh_modify	delay 0 every 1 one 200000 page 2000000

pair_style granular 
pair_coeff 1 1 hertz 40000000.0 150.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
pair_coeff 2 2 hertz 20000000.0 100.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall
pair_coeff 1 2 hertz 30000000.0 125.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall


timestep	0.00000001 # only related to the mass of particle, very consertative
run_style	verlet


compute 	2 all pair/local dist dx dy dz force fx fy fz 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 

variable        volume atom "4.0*PI* c_4[2] * c_4[2] * c_4[2] /3.0"
compute         vol all reduce sum v_volume


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 


fix         9 all press/berendsen x 100000 100000 100 y 100000 100000 100 z 100000 100000 100

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_vol)/(c_vol)"

fix 		2 all print ${interval_fix} "${step} ${xxstress} ${yystress} ${zzstress} ${xystress} ${xzstress} ${yzstress} ${xlength} ${ylength} ${zlength}" file  random_mean_stresses_iso.txt screen no


fix 		4 all print ${interval_fix} "${step} ${coordNum_new}" file random_coord.txt screen no

fix 		6 all print ${interval_fix} "${step} ${voidratio}" file random_voidratio.txt screen no

fix		3 all nve/sphere

# fix          5 all momentum/gran  


# !!!!!!!!WARNING!!!£¬ CHANGE the steps and folder name!!!!
# dump files

 dump		1 all local ${interval_dump} random_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]  

# remove the "id" in order to be consistent with TOM

dump         2 all custom ${interval_dump} random_connectivity/dump*.connectivity_1 id c_5 x y z diameter

dump         3 all atom ${interval_dump}   random_atom/dump*.atom
dump_modify  3 scale no

dump         4 all custom ${interval_dump}   random_stress/dump*.stress_1 id diameter c_6[1] c_6[2] c_6[3] c_6[4] c_6[5] c_6[6]

dump         5 all custom ${interval_dump} random_velocity/dump*.all_velocities c_4[1] c_4[2] c_4[3] c_4[4] c_4[5] c_4[6] c_4[7] c_4[8] c_4[9] c_4[10] c_4[11]


#========================================================================================

# !!!!!!!!WARNING!!!£¬ restart name!!!!
restart         10000 restart_filerandom

#=========================================================================================
run		999999999