So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] := With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points, d = (2 r)}, c + d*Normalize@(c - p) ]; generateRandomWalking[] := With[{p = RandomChoice[points], \[Theta] = RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)}, p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]], d*Cos[\[Theta]]} ]; tryPoint[s_] := If[AllTrue[s, r < # < m - r &] && AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]]; 

Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7; m = 20; points = {}; tryPoint@RandomReal[{r, m - r}, 3]; With[{now = Now}, Monitor[While[Length[points] < 1825, tryPoint@generateRandomAnchoring[] ], Column@{Now - now, Length@points}] ] 

It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False] 

And this is what we see:

Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

###Edit

Note that one could also add a weight-tracking scheme for the random choice of points (in the second of the methods). I.e., points that have been chosen and for which the addition failed decrease in weight by some percentage. The likelihood that this will be a significant improvement is low, but it could eke out some performance gains.

So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] := With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points, d = (2 r)}, c + d*Normalize@(c - p) ]; generateRandomWalking[] := With[{p = RandomChoice[points], \[Theta] = RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)}, p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]], d*Cos[\[Theta]]} ]; tryPoint[s_] := If[AllTrue[s, r < # < m - r &] && AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]]; 

Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7; m = 20; points = {}; tryPoint@RandomReal[{r, m - r}, 3]; With[{now = Now}, Monitor[While[Length[points] < 1825, tryPoint@generateRandomAnchoring[] ], Column@{Now - now, Length@points}] ] 

It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False] 

And this is what we see:

Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

Edit

Note that one could also add a weight-tracking scheme for the random choice of points (in the second of the methods). I.e., points that have been chosen and for which the addition failed decrease in weight by some percentage. The likelihood that this will be a significant improvement is low, but it could eke out some performance gains.

So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] := With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points, d = (2 r)}, c + d*Normalize@(c - p) ]; generateRandomWalking[] := With[{p = RandomChoice[points], \[Theta] = RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)}, p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]], d*Cos[\[Theta]]} ]; tryPoint[s_] := If[AllTrue[s, r < # < m - r &] && AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]]; 

Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7; m = 20; points = {}; tryPoint@RandomReal[{r, m - r}, 3]; With[{now = Now}, Monitor[While[Length[points] < 1825, tryPoint@generateRandomAnchoring[] ], Column@{Now - now, Length@points}] ] 

It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False] 

And this is what we see:

Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

So using another random walk-like method (well I tried two which I realized are equivalent) I was able to do this pretty quickly and compactly:

generateRandomAnchoring[] := With[{p = RandomReal[{r, m - r}, 3], c = RandomChoice@points, d = (2 r)}, c + d*Normalize@(c - p) ]; generateRandomWalking[] := With[{p = RandomChoice[points], \[Theta] = RandomReal[2 \[Pi]], \[Phi] = RandomReal[2 \[Pi]], d = (2 r)}, p + {d*Cos[\[Phi]] Sin[\[Theta]], d*Sin[\[Theta]] Sin[\[Phi]], d*Cos[\[Theta]]} ]; tryPoint[s_] := If[AllTrue[s, r < # < m - r &] && AllTrue[points, (Norm[# - s] >= 2 r) &], AppendTo[points, s]]; 

Both of these take a random pregenerated point and attach a point to it at a random position, just in one case it generates a spherical polar coordinate and in the other it just attaches a randomly generated point. They're entirely equivalent as best I can reason out.

Then we'll just initialize the state and update the point generating code from wxffles answer so that it also spits out the time elapsed:

r = .7; m = 20; points = {}; tryPoint@RandomReal[{r, m - r}, 3]; With[{now = Now}, Monitor[While[Length[points] < 1825, tryPoint@generateRandomAnchoring[] ], Column@{Now - now, Length@points}] ] 

It runs for about a minute before hitting the 1825 mark. And we'll use the same viewer code as wxffles too:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; Graphics3D[{cube, Sphere[#, 0.7] & /@ points}, Boxed -> False] 

And this is what we see:

Then just adjust the number of spheres and timing to your liking and you have a your packed cube.

###Edit

Note that one could also add a weight-tracking scheme for the random choice of points (in the second of the methods). I.e., points that have been chosen and for which the addition failed decrease in weight by some percentage. The likelihood that this will be a significant improvement is low, but it could eke out some performance gains.