4
$\begingroup$

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.

enter image description here

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?

$\endgroup$
  • 1
    $\begingroup$ Aren't you already aware of RandomPoint? $\endgroup$ – Rahul Dec 20 '15 at 15:54
  • 2
    $\begingroup$ 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. $\endgroup$ – Rahul Dec 20 '15 at 17:32
  • 1
    $\begingroup$ @Rahul your comment is surely deserving of being an answer! $\endgroup$ – Oleksandr R. Dec 20 '15 at 17:57
  • $\begingroup$ @Rahul I was aware of RandomPoint, but not of how to apply it! $\endgroup$ – dr.blochwave Dec 20 '15 at 19:12
  • $\begingroup$ @Rahul as Oleksandr says, if you make your comment an answer, it's exactly what I'm looking for. $\endgroup$ – dr.blochwave Dec 20 '15 at 20:18
4
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Dec 20 '15 at 17:54
2
$\begingroup$

Answer courtesy of Rahul:

s = DiscretizeRegion[
        ImplicitRegion[
          Min[EuclideanDistance[{x, y, z}, #] & /@ carbonXYZ] == rc + ra, 
        {x, y, z}], 
       {{-1, 28}, {-2, 23}, {-5, 5}}, 
       MaxCellMeasure -> 0.1] 

RandomPoint[s]
$\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.