Return to Answer

Here's a variation that attempts to do more than one bounded random walk, where the starting point of the $n$th random walk is a randomly chosen point from the $(n-1)$-th random walk:

sphrPak2[s_, r_, n_Integer?Positive, walks_Integer?Positive, its_Integer?Positive] := Module[{sh = s/2, j, k, p0, pc, sphList}, While[ p0 = RandomReal[{-sh, sh}, 3]; And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]]; sphList = {p0}; j = 0; While[Length[sphList] < n && j < walks, k = 0; While[Length[sphList] < n && k < its, pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]]; If[(PolyhedronData["Cube", "RegionFunction"] @@ ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc]) && (And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]), k = 0; AppendTo[sphList, pc]; p0 = pc, ++k]; ]; If[Length[sphList] < n, p0 = RandomChoice[Map[EuclideanDistance[#, sphList[[1]]] &, sphList] -> sphList]]; ++j; ]; If[j == walks, Print[StringForm["Only `1` spheres were generated.", Length[sphList]]]]; Graphics3D[Sphere[sphList, r], Axes -> Automatic, PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]] 

Try it out:

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"]; sphrPak2[20, 0.7, 1500, 1, 1000]] Only 662 spheres were generated. 

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"]; sphrPak2[20, 0.7, 1500, 20, 1000]] Only 1335 spheres were generated. 

sphrPak[s_, r_, n_Integer?Positive, lim_Integer?Positive] := Module[{sh = s/2, k, p0, pc, sphList}, (* generate initial sphere within box *) While[ p0 = RandomReal[{-sh, sh}, 3]; And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]]; sphList = {p0}; k = 0; While[Length[sphList] < n && k < lim, (* random direction *) pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]]; If[(* is center within the cube? *) PolyhedronData["Cube", "RegionFunction"] @@ ScalingTransform[ConstantArray[1/s, 3]][pc] && (* is the sphere within the cube? *) (And @@ Thread[Map[Min[sh - #, sh + #] &, pc] >= r]) && (* does the sphere not overlap with other spheres? *) (And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]), k = 0; AppendTo[sphList, pc]; p0 = pc, (* else *) ++k]; ]; If[k == lim, Print[StringForm["Only `1` spheres were generated.", Length[sphList]]]]; Graphics3D[Sphere[sphList, r], Axes -> Automatic, PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]]; sphrPak[s_, r_, n_Integer] := sphrPak[s, r, n, Quotient[2 n, 3]] 

sphrPak[s_, r_, n_Integer?Positive, lim_Integer?Positive] := Module[{sh = s/2, k, p0, pc, sphList}, (* generate initial sphere within box *) While[ p0 = RandomReal[{-sh, sh}, 3]; And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]]; sphList = {p0}; k = 0; While[Length[sphList] < n && k < lim, (* center for new sphere chosen in random direction *) pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]]; If[(* is sphere within the cube? *) PolyhedronData["Cube", "RegionFunction"] @@ ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc] && (* does the sphere not overlap with other spheres? *) (And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]), k = 0; AppendTo[sphList, pc]; p0 = pc, (* else *) ++k]; ]; If[k == lim, Print[StringForm["Only `1` spheres were generated.", Length[sphList]]]]; Graphics3D[Sphere[sphList, r], Axes -> Automatic, PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]]; sphrPak[s_, r_, n_Integer] := sphrPak[s, r, n, Quotient[2 n, 3]] 

More attempts:

BlockRandom[SeedRandom[4092, Method -> "MersenneTwister"]; sphrPak[20, 0.7, 1500]] Only 805 spheres were generated. 

BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"]; sphrPak[20, 0.7, 1500]] Only 932 spheres were generated.