Generate a set of 3D coordinates subject to constraints

Question

Generating a set of random 3D coordinates is ok, e.g.

positions = RandomReal[{0, 10}, {10, 3}];
radius = 1;
Graphics3D[{Red, Sphere[positions, radius]}, PlotRangePadding -> 2, 
 Lighting -> "Neutral"]

Some random 3D points

Now I want to add a constraint, so that the distance between any point and its nearest neighbour is fixed.

A good way of thinking about it would be like generating a virtual molecule, where the atoms can't be too close to each other, but nor can they be too far apart otherwise they wouldn't be bonded.

I had a prototype made in C++ a while back that just iterated over and over, adding a new point each time the randomly generated coordinates happened to fit the constraint. For example, first I placed a single point with random coordinates and called it i. Then I created a new point, j, with random coordinates. I then defined a distance between point i and point j:

euclideandistance = Sqrt[
    (positions[[i,1]]-positions[[j,1]])^2 + 
    (positions[[i,2]]-positions[[j,2]])^2 + 
    (positions[[i,3]]-positions[[j,3]])^2
]

Then if the constraint, based on bondlength,

0.9*bondlength < euclideandistance < 1.1*bondlength

was satisfied, I'd add point j and start over with a new point k, comparing the distance to its nearest neighbour (which of course could be i or j now) and repeating...

This seems somewhat inefficient...it was a quick-and-dirty prototype after all! Coming from a C++ background to Mathematica, is there a good way to do this in Mathematica that doesn't involve a While loop?

Update

Courtesy of Yves Klett, this question (in 2D and with only a lower bound) Efficient way to generate random points with a predefined lower bound on their pairwise Euclidean distance is what I had in mind, I didn't spot it in my search before posting this.

I'll go away and have a play myself when I can tomorrow to apply it to my scenario, and post what I come up with when I've done so (though any answers in the meantime are fine! It all adds to the learning experience for me).

Related (?): Efficient way to generate random points with a predefined lower bound on their pairwise Euclidean distance — Yves Klett, May 23 '14 at 19:39
Yep I think that is related, in 2D but that's no problem. I hadn't spotted that question when I was searching, but I guess the title isn't what I'd thought of. — dr.blochwave, May 23 '14 at 19:41
Do you want to generate arbitrary-length molecules or do you have some upper bound? — Yves Klett, May 23 '14 at 19:44
It doesn't matter too much to me, as I'm just looking to generate some random shapes with the constraint I mentioned (that's the key thing). — dr.blochwave, May 23 '14 at 19:46
Looking forward to your solution (although Andys solution does use While) :D — Yves Klett, May 23 '14 at 19:49
Yes, I realised it as I was looking! But it's nice and quick, which is the important thing. Thanks for your help! — dr.blochwave, May 23 '14 at 20:05

score 3 · Accepted Answer · edited Apr 13 '17 at 12:55

So here's my edit to the answer given in Efficient way to generate random points with a predefined lower bound on their pairwise Euclidean distance

First define the modified function, which now has the upper bound included, as well as a third coordinate:

findPoints = 
  Compile[{{n, _Integer}, {low, _Real}, {high, _Real}, {minD, _Real},{maxD, _Real}}, 
   Block[{data = RandomReal[{low, high}, {1, 3}], k = 1, rv, temp}, 
    While[k < n, rv = RandomReal[{low, high}, 3];
     temp = Transpose[Transpose[data] - rv];
     If[Min[Sqrt[(#.#)] & /@ temp] > minD && 
       Min[Sqrt[(#.#)] & /@ temp] < maxD, data = Join[data, {rv}];
      k++;];];
    data]];

Next job is to specify the parameters to try it out:

npts = 50;
minD = 3;
maxD = 4;
low = 0;
high = 100;

Then benchmarking, checking the constraints etc. shows it's pretty quick (as the original answer was), and meets the min & max constraints I've set.

AbsoluteTiming[pts = findPoints[npts, low, high, minD, maxD];]
Min[MapThread[
   EuclideanDistance, {pts, Nearest[pts][#, 2][[-1]] & /@ pts}]] > minD
Min[MapThread[
   EuclideanDistance, {pts, Nearest[pts][#, 2][[-1]] & /@ pts}]] < maxD
Length[pts] == npts

(* Out={0.008019, Null} *)
(* Out=True *)
(* Out=True *)
(* Out=True *)

Then plotting it with spheres, as I was in my question:

radius = 1;
Graphics3D[{Red, Sphere[pts, radius]}, PlotRangePadding -> 2, 
 Lighting -> "Neutral"]

Gives the following result:

enter image description here

And if I increase the upper bound, e.g. to 50, I get a much more dispersed solution (as expected!):

enter image description here

Generate a set of 3D coordinates subject to constraints

1 Answers1

Linked