2
$\begingroup$

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

$\endgroup$
  • 2
  • $\begingroup$ 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. $\endgroup$ – dr.blochwave May 23 '14 at 19:41
  • $\begingroup$ Do you want to generate arbitrary-length molecules or do you have some upper bound? $\endgroup$ – Yves Klett May 23 '14 at 19:44
  • $\begingroup$ 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). $\endgroup$ – dr.blochwave May 23 '14 at 19:46
  • 1
    $\begingroup$ Looking forward to your solution (although Andys solution does use While) :D $\endgroup$ – Yves Klett May 23 '14 at 19:49
3
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.