atom_style	granular
boundary	p p p
newton		off
echo		both
communicate	single vel yes
units		si
processors      1 1 1

region		domain block 0.0 1.0 0.0 1.0 0.0 1.0 units box
create_box	1 domain

neighbor	0.001 bin
neigh_modify	every 1 delay 0

#Materials
hard_particles  yes
fix 		m1 all property/global youngsModulus peratomtype 8.7e9
fix 		m2 all property/global poissonsRatio peratomtype 0.3
fix 		m3 all property/global coefficientRestitution peratomtypepair 1 0.95
fix 		m4 all property/global coefficientFriction peratomtypepair 1 0.0 
fix 		m5 all property/global coefficientRollingFriction peratomtypepair 1 0.0	

#pair_style
pair_style	gran model hertz tangential history
pair_coeff	* * 
timestep	1.0e-7

#distributions for insertion
fix 		pts1 all particletemplate/sphere 123457 atom_type 1 density constant 2500 radius constant 0.15
fix		pdd1 all particledistribution/discrete 15485867 1 pts1 1.0

#region and insertion
group 		nve_group region domain
region		main block 0.3 0.7 0.0 1.0 0.3 0.7 units box

#particle insertion
fix 		ins1 nve_group insert/pack seed 32452867 distributiontemplate pdd1 vel gaussian 0. 0. 0. 5. 5. 5. omega constant 0. 0. 0. &
		insert_every once overlapcheck yes region main ntry_mc 100000000 particles_in_region 5

#integrator
fix		integr nve_group nve/sphere

#Run one step to get particles inserted to avoid initial division by zero in data
thermo_style 	custom step
run    		1

#calculations
compute         ppc all contact/atom
compute         allcontactsinsystem all reduce sum c_ppc
compute 	rot all erotate
compute         peratom nve_group stress/atom
compute         p nve_group reduce/region main sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4]
variable        pyy equal -(c_p[2])
variable        pxy equal (c_p[4])
variable        trans equal 2*ke/(3*mass(nve_group))-(1/3)*(vcm(nve_group,x)^2+vcm(nve_group,y)^2+vcm(nve_group,z)^2)
variable        rotat equal c_rot/(mass(nve_group))-(1/3)*(inertia(nve_group,xx)*omega(nve_group,x)^2+inertia(nve_group,yy)*omega(nve_group,y)^2+inertia(nve_group,zz)*omega(nve_group,z)^2)/mass(nve_group)
variable        num_cont equal c_allcontactsinsystem/2

#print info
thermo_style 	custom step time atoms ke v_trans v_rotat v_pyy v_pxy v_num_cont
thermo		1000
thermo_modify	lost ignore norm no

run             25000 upto