How to obtain interstitial packing

Thank you Chris, I was thinking about the packing due to gravity drop. I used the pour command and got the particles in to a cylindrical region. But I have another problem. Iam not able to delete atoms in a new small cylindrical region inside the first cylindrical region. The input command is like this

### Initialization
# Preliminaries
units si
atom_style sphere
atom_modify map array
boundary f f f
newton off
communicate single vel yes

#processors 2 2 3

# Declare domain

region domain block -0.017 0.033 -0.017 0.033 -0.001 0.011 units box
create_box 1 domain

# Neighbor listing

neighbor 0.003 bin
neigh_modify delay 0

### Setup

# Material and interaction properties required

fix m1 all property/global youngsModulus peratomtype 5e6 5e6
fix m2 all property/global poissonsRatio peratomtype 0.31 0.31
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.1 0.1 0.1 0.1
fix m4 all property/global coefficientFriction peratomtypepair 2 0.2 0.2 0.2 0.2

# Gravity

fix grav all gravity 9.81 vector 0.0 0.0 -1.0

region void cylinder z 0. 0. 0.015 0. 0.009 units box
region delete_region cylinder z 0. 0. 0.006 0. 0.009 units box
delete_atoms region delete_region

fix ins all pour/legacy 10000 1 12345 &
region void vol 0.5 100 rate 1. diam 10e-06 200e-06

# Define the physics

pair_style gran/hertz/history
pair_coeff * *

### Detailed settings

# Integrator

fix integrate all nve/sphere

# Time step

timestep 0.00000001

# Thermodynamic output settings

thermo_style custom step atoms ke cpu
thermo 10000
thermo_modify norm no lost ignore

# Check time step and initialize dump file

fix ctg all check/timestep/gran 1 0.01 0.01
run 1
unfix ctg

# Create imaging information

dump dmp all custom 10000 Desktop/new6/dump.1 id type type x y z ix iy iz vx vy vz &
fx fy fz omegax omegay omegaz radius

run 50000

What could be my mistake ? Please help.

Regards
Hari