Random periodic packing with different volume fractions

I am trying to generate random periodic packing of spheres in a cube with different volume fraction. I am looking for random packing of spheres with volume fraction as one parameter. I have used yade to generate random dense packing using the script below. Is there any easy way to generate packing for different volume fractions in yade?

from yade import pack
pack.randomPeriPack(10,(100,100,100), rRelFuzz = 0.1).toSimulation()
O.step() # to initialize bounding boxes
xImages = [] # list of "edge" particles

for b in O.bodies:
 xmin,xmax = O.cell.wrap(b.bound.min)[0], O.cell.wrap(b.bound.max)[0] # wrap is important as physically the particles can be anywhere

 if xmin > xmax: # this means that bounding box crosses periodic cell. You can define various different conditions..
  xImages.append(b)

pr = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in xImages] # list or [pos,radius] of "edge" particles. Now pos is always inside the cell

for i,(pos,r) in enumerate(pr):
 shift = Vector3(O.cell.size[0],0,0) * (1 if pos[0]<0.5*O.cell.size[0] else -1) # determine shift to create a periodic image +x... direction if particles is near -x face and vice versa
 pr[i][0] += shift

# saves images into a file
f = open("/home/nikita/1_1_x_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr)
f.close()

yImages = []

for b in O.bodies:
 ymin,ymax = O.cell.wrap(b.bound.min)[1], O.cell.wrap(b.bound.max)[1]
 if ymin > ymax:
  yImages.append(b)

pr1 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in yImages]

for i,(pos,r) in enumerate(pr1):
 shift = Vector3(0,O.cell.size[1],0) * (1 if pos[1]<0.5*O.cell.size[1] else -1)
 pr1[i][0] += shift

# saves images into a file
f = open("/home/nikita/1_1_y_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr1)
f.close()

zImages = []

for b in O.bodies:
 zmin,zmax = O.cell.wrap(b.bound.min)[2], O.cell.wrap(b.bound.max)[2]
 if zmin > zmax:
  zImages.append(b)

pr2 = [[O.cell.wrap(b.state.pos),b.shape.radius] for b in zImages]

for i,(pos,r) in enumerate(pr2):
 shift = Vector3(0,0,O.cell.size[2]) * (1 if pos[2]<0.5*O.cell.size[2] else -1)
 pr2[i][0] += shift

# saves images into a file
f = open("/home/nikita/1_1_z_1.dat","w")
f.write("x y z r\n")
f.writelines("%g %g %g %g\n"%(pos[0],pos[1],pos[2],r) for pos,r in pr2)
f.close()

# saving normal spheres
f = open("/home/nikita/normal_1.dat", "w")
f.write('x y z r\n')
for b in O.bodies:
 p,r = O.cell.wrap(b.state.pos),b.shape.radius
 f.write('%g %g %g %g\n' %(p[0],p[1],p[2],r))
f.close()