# Initialize Simulation
units si
atom_style granular
atom_modify map array
boundary f f f
newton off
communicate single vel yes

# Define Region and Create Box
region box_region block 0 0.038 0 0.038 0 0.080 units box
create_box 1 box_region

# Neighbor Settings
neighbor 0.001 bin
neigh_modify delay 0

# Material Properties for Hertz-Mindlin Contact Model
fix m1 all property/global youngsModulus peratomtype 45e6
fix m2 all property/global poissonsRatio peratomtype 0.27
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
fix m4 all property/global coefficientFriction peratomtypepair 1 0.4
#fix m5 all property/global coefficientRollingFriction peratomtypepair 1 0.05

# Define Contact Model (Hertz-Mindlin)
pair_style gran model hertz tangential history
pair_coeff * *

# Set Timestep
timestep 0.000005

fix			grav all gravity 9.81 vector 0.0 0.0 -1.0

# Create Fixed Walls
fix wall_x1 all wall/gran model hertz tangential history primitive type 1 xplane 0.0
fix wall_x2 all wall/gran model hertz tangential history primitive type 1 xplane 0.038
fix wall_y1 all wall/gran model hertz tangential history primitive type 1 yplane 0.0
fix wall_y2 all wall/gran model hertz tangential history primitive type 1 yplane 0.038
fix wall_bottom all wall/gran model hertz tangential history primitive type 1 zplane 0.0
fix wall_top all wall/gran model hertz tangential history primitive type 1 zplane 0.08

# Define Group for Insertion
group nve_group region box_region

# Particle templates with unique prime seeds for 50 particles
fix pt01 all particletemplate/sphere 15485863 atom_type 1 density constant 2700 radius constant 7.165000e-05
fix pt02 all particletemplate/sphere 15485867 atom_type 1 density constant 2700 radius constant 1.213664e-04
fix pt03 all particletemplate/sphere 15485917 atom_type 1 density constant 2700 radius constant 1.710827e-04
fix pt04 all particletemplate/sphere 15485927 atom_type 1 density constant 2700 radius constant 2.207990e-04
fix pt05 all particletemplate/sphere 15485933 atom_type 1 density constant 2700 radius constant 2.705153e-04
fix pt06 all particletemplate/sphere 15485941 atom_type 1 density constant 2700 radius constant 3.202316e-04
fix pt07 all particletemplate/sphere 15485959 atom_type 1 density constant 2700 radius constant 3.699480e-04
fix pt08 all particletemplate/sphere 15485989 atom_type 1 density constant 2700 radius constant 4.196643e-04
fix pt09 all particletemplate/sphere 15485993 atom_type 1 density constant 2700 radius constant 4.693806e-04
fix pt10 all particletemplate/sphere 15486013 atom_type 1 density constant 2700 radius constant 5.190970e-04
fix pt11 all particletemplate/sphere 15486041 atom_type 1 density constant 2700 radius constant 5.688135e-04
fix pt12 all particletemplate/sphere 15486047 atom_type 1 density constant 2700 radius constant 6.185295e-04
fix pt13 all particletemplate/sphere 15486059 atom_type 1 density constant 2700 radius constant 6.682460e-04
fix pt14 all particletemplate/sphere 15486071 atom_type 1 density constant 2700 radius constant 7.179620e-04
fix pt15 all particletemplate/sphere 15486101 atom_type 1 density constant 2700 radius constant 7.676785e-04
fix pt16 all particletemplate/sphere 15486139 atom_type 1 density constant 2700 radius constant 8.173950e-04
fix pt17 all particletemplate/sphere 15486157 atom_type 1 density constant 2700 radius constant 8.671110e-04
fix pt18 all particletemplate/sphere 15486173 atom_type 1 density constant 2700 radius constant 9.168275e-04
fix pt19 all particletemplate/sphere 15486181 atom_type 1 density constant 2700 radius constant 9.665440e-04
fix pt20 all particletemplate/sphere 15486193 atom_type 1 density constant 2700 radius constant 1.016260e-03
fix pt21 all particletemplate/sphere 15486209 atom_type 1 density constant 2700 radius constant 1.065976e-03
fix pt22 all particletemplate/sphere 15486221 atom_type 1 density constant 2700 radius constant 1.115693e-03
fix pt23 all particletemplate/sphere 15486227 atom_type 1 density constant 2700 radius constant 1.165409e-03
fix pt24 all particletemplate/sphere 15486241 atom_type 1 density constant 2700 radius constant 1.215126e-03
fix pt25 all particletemplate/sphere 15486257 atom_type 1 density constant 2700 radius constant 1.264842e-03
fix pt26 all particletemplate/sphere 15486259 atom_type 1 density constant 2700 radius constant 1.314558e-03
fix pt27 all particletemplate/sphere 15486277 atom_type 1 density constant 2700 radius constant 1.364275e-03
fix pt28 all particletemplate/sphere 15486281 atom_type 1 density constant 2700 radius constant 1.413991e-03
fix pt29 all particletemplate/sphere 15486283 atom_type 1 density constant 2700 radius constant 1.463708e-03
fix pt30 all particletemplate/sphere 15486287 atom_type 1 density constant 2700 radius constant 1.513424e-03
fix pt31 all particletemplate/sphere 15486347 atom_type 1 density constant 2700 radius constant 1.563140e-03
fix pt32 all particletemplate/sphere 15486421 atom_type 1 density constant 2700 radius constant 1.612856e-03
fix pt33 all particletemplate/sphere 15486433 atom_type 1 density constant 2700 radius constant 1.662572e-03
fix pt34 all particletemplate/sphere 15486437 atom_type 1 density constant 2700 radius constant 1.712289e-03
fix pt35 all particletemplate/sphere 15486451 atom_type 1 density constant 2700 radius constant 1.762005e-03
fix pt36 all particletemplate/sphere 15486469 atom_type 1 density constant 2700 radius constant 1.811722e-03
fix pt37 all particletemplate/sphere 15486481 atom_type 1 density constant 2700 radius constant 1.861438e-03
fix pt38 all particletemplate/sphere 15486487 atom_type 1 density constant 2700 radius constant 1.911154e-03
fix pt39 all particletemplate/sphere 15486511 atom_type 1 density constant 2700 radius constant 1.960871e-03
fix pt40 all particletemplate/sphere 15486517 atom_type 1 density constant 2700 radius constant 2.010587e-03
fix pt41 all particletemplate/sphere 15486533 atom_type 1 density constant 2700 radius constant 2.060303e-03
fix pt42 all particletemplate/sphere 15486557 atom_type 1 density constant 2700 radius constant 2.110020e-03
fix pt43 all particletemplate/sphere 15486571 atom_type 1 density constant 2700 radius constant 2.159736e-03
fix pt44 all particletemplate/sphere 15486589 atom_type 1 density constant 2700 radius constant 2.209452e-03
fix pt45 all particletemplate/sphere 15486649 atom_type 1 density constant 2700 radius constant 2.259168e-03
fix pt46 all particletemplate/sphere 15486671 atom_type 1 density constant 2700 radius constant 2.308884e-03
fix pt47 all particletemplate/sphere 15486673 atom_type 1 density constant 2700 radius constant 2.358601e-03
fix pt48 all particletemplate/sphere 15486703 atom_type 1 density constant 2700 radius constant 2.408317e-03
fix pt49 all particletemplate/sphere 15486707 atom_type 1 density constant 2700 radius constant 2.458033e-03
fix pt50 all particletemplate/sphere 32452843 atom_type 1 density constant 2700 radius constant 2.507750e-03


# now the PSD distribution
fix psd1 all particledistribution/discrete/numberbased 49979693  50 pt01 0.119474 pt02 0.095764 pt03 0.078340 pt04 0.065553 pt05 0.055887 pt06 0.048371 pt07 0.042384 pt08 0.037520 pt09 0.033501 pt10 0.030133 pt11 0.027277 pt12 0.024829 pt13 0.022712 pt14 0.020866 pt15 0.019246 pt16 0.017814 pt17 0.016542 pt18 0.015405 pt19 0.014385 pt20 0.013466 pt21 0.012633 pt22 0.011878 pt23 0.011189 pt24 0.010559 pt25 0.009982 pt26 0.009451 pt27 0.008962 pt28 0.008510 pt29 0.008092 pt30 0.007704 pt31 0.007343 pt32 0.007007 pt33 0.006694 pt34 0.006401 pt35 0.006127 pt36 0.005869 pt37 0.005628 pt38 0.005401 pt39 0.005188 pt40 0.004986 pt41 0.004796 pt42 0.004617 pt43 0.004447 pt44 0.004286 pt45 0.004134 pt46 0.003990 pt47 0.003853 pt48 0.003722 pt49 0.003598 pt50 0.003480


region insert_region block 0.003 0.035 0.003 0.035 0.000 0.070 units box



# Particle insertion using that distribution
fix insert1 nve_group insert/pack seed 67867979 distributiontemplate psd1 maxattempt 10000 insert_every 100 overlapcheck yes all_in yes vel constant 0.0 0.0 -1.0 region insert_region particles_in_region 5000 ntry_mc 10000



# Apply NVE integration to all inserted particles
fix integr all nve/sphere


# Output settings to include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no

# Insert the first particles to avoid an empty dump
run 100000

unfix insert1
dump dmp all custom/vtk 800 post/TX_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius

# Insert particles
run 100000