# Generate a set of 3D coordinates subject to constraints

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

positions = RandomReal[{0, 10}, {10, 3}];
Lighting -> "Neutral"]


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).

• May 23, 2014 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. May 23, 2014 at 19:41
• Do you want to generate arbitrary-length molecules or do you have some upper bound? May 23, 2014 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). May 23, 2014 at 19:46
• Looking forward to your solution (although Andys solution does use While) :D May 23, 2014 at 19:49

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];]
EuclideanDistance, {pts, Nearest[pts][#, 2][[-1]] & /@ pts}]] > minD
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;