# Randomly generate points on a van der Waals surface

I'm looking to randomly place an atom $\mathbf{A}$ onto a van der Waals surface formed by atoms $\mathbf{B}$, such that the distance between the randomly-generated atom $\mathbf{A}$ and the nearest surface atom $\mathbf{B}_i$ is equal to the sum of the van der Waals radii $r$. Essentially, the problem boils down to a sphere rolling around on an uneven surface.

I can generate a van der Waals surface for graphene with the following code:

latticeGenerate[basisvec_List, numofcell_List, base_List] :=
Module[{x, y, basex, basey}, {x, y} =
Transpose[Tuples[Range /@ numofcell].basisvec];
{basex, basey} = Transpose[base];
Transpose[{Join @@ ((x + #) & /@ basex),
Join @@ ((y + #) & /@ basey)}]]


Defining the C-C bond length and van der Waals radius,

a = 1.4145;
rc = 1.675;

carbonXYZ = Flatten[{#, 0}] & /@
latticeGenerate[{{3 a, 0}, {0, Sqrt[3] a}}, {5,7},
{{0, 0}, {a, 0}, {-a/2, Sqrt[3] a/2}, {3 a/2, Sqrt[3] a/2}}];

vdWsurf = Sphere[carbonXYZ, rc];
Graphics3D[{Gray, vdWsurf}, PlotRangePadding -> 2, Lighting -> "Neutral"]


I'd then like to place an atom with radius 2.5 at random on the top of this surface, e.g.

I did think about a naive approach using FindInstance[], but this just runs for a long time without a result:

ra = 2.5;
minDist[pt_] := Min[Sqrt[(#.#)] & /@ (Transpose[Transpose[carbonXYZ] - pt])]
FindInstance[minDist[{x, y, z}] == (ra + rc) &&
5 < x < 15 && 5 < y < 15 && 0 < z < 5,
{x, y, z}, Reals]


Could the Region[] functionality of Mathematica be exploited here? I'm looking to generate lots of these systems, so speed is of the essence! Perhaps RandomPoint[] can be utilised in this case?

• Aren't you already aware of RandomPoint? – Rahul Dec 20 '15 at 15:54
• For example, set s = DiscretizeRegion[ImplicitRegion[Min[EuclideanDistance[{x, y, z}, #] & /@ carbonXYZ] == rc + ra, {x, y, z}], {{-1, 28}, {-2, 23}, {-5, 5}}, MaxCellMeasure -> 0.1] and then RandomPoint[s] is extremely fast. – Rahul Dec 20 '15 at 17:32
• @Rahul your comment is surely deserving of being an answer! – Oleksandr R. Dec 20 '15 at 17:57
• @Rahul I was aware of RandomPoint, but not of how to apply it! – dr.blochwave Dec 20 '15 at 19:12
• @Rahul as Oleksandr says, if you make your comment an answer, it's exactly what I'm looking for. – dr.blochwave Dec 20 '15 at 20:18

One possibility is to set it up as a minimization problem. For example:

NMinimize[(minDist[{x, y, z}] - (ra + rc))^2, 5<x<15 && 5<y<15 && 0<z<5, {x, y, z}]

{2.43076*10^-21, {x -> 7.32182, y -> 13.4219, z -> 4.14947}}


You can find minima in different locations by changing the constraints:

NMinimize[(minDist[{x, y, z}] - (ra + rc))^2, 5<x<6 && 5<y<15 && 0<z< 5, {x, y, z}]

{7.66876*10^-19, {x -> 5.50514, y -> 13.7025, z -> 4.0791}}


Alternatively, FindMinimum seems faster, and allows you to directly give an initial point:

FindMinimum[{(minDist[{x, y, z}] - (ra + rc))^2, 5<x<15 && 5<y<15 && 0<z<5},
{{x, 10}, {y, 10}, {z, 2.5}}]

{1.21243*10^-12, {x -> 10.2024, y -> 10.4351, z -> 4.11543}}


In line with Oleksandr's comments, you can also try

FindMinimum[(minDist[{x, y, z}] - (ra + rc))^2,
{{x, 8.}, {y, 10}, {z, 3.5}}, Method -> "PrincipalAxis"]


Try out the different methods, some of which may give you warnings (and require adjusting various parameters).

• I think I would try FindMinimum without the constraints if speed is important. If constraints are given it uses the nonlinear interior point algorithm, which is (the only one?) written in top-level code, and slow. Since all of its methods are local optimizers, you can rely on it not going far from the starting point. – Oleksandr R. Dec 20 '15 at 17:54

s = DiscretizeRegion[