2
$\begingroup$

I want to create a polymer chain (in 2D) of a given length such that:
1. First monomer is at {0,0}
2. All other monomers are in the positive x half-plane
3. The distance between two bonded monomers is r0
4. No two non-bonded monomers come closer to each other than rc (rc>r0)

Edit: I found a way to implement these rules:
r0 FoldList[AngleVector,{0, 0},RandomReal[{-1, 1} ArcCos[rc/(2 r0)], n - 1]]
However, this solution does not sample all possible configurations, but only a small subset.

$\endgroup$
  • 1
    $\begingroup$ So... a self-avoiding random walk? Did you try searching the site for random walk implementations? $\endgroup$ – J. M. will be back soon Apr 11 '18 at 13:23
  • $\begingroup$ @J.M.needshelp. It is not exactly a self-avoiding random walk. But I did search for random walk. This is the closest I got to what I want. $\endgroup$ – pukkandan Apr 11 '18 at 13:38
  • $\begingroup$ Your fourth condition sounds like "self-avoiding" to me... $\endgroup$ – J. M. will be back soon Apr 11 '18 at 13:40
  • $\begingroup$ Yes, it has to be self avoiding. But rather than just "not crossing" each other, it has to stay within a certain distance from each other. For example, this is self-avoiding, but does not satisfy rule 4 $\endgroup$ – pukkandan Apr 11 '18 at 13:59
8
$\begingroup$

This could get you started. As always, we use Nearest to determine collision, but put the cheaper collision detection with the left half plane in front. This features also a buffer reservoir for future random steps since creating many random numbers at once is usually much more performant.

r0 = 0.1;
rc = 0.15;
dim = 2;
maxtrials = 200;
reservoircounter = 1;
chaincounter = 2;
maxchainlength = 10000;
reservoirLength = 1000;
getReservoir[n_, r0_] := RandomPoint[Sphere[ConstantArray[0., dim], r0], n];
chain = ConstantArray[0., {maxchainlength, dim}];
reservoir = getReservoir[reservoirLength, r0];
chain[[2]] = x = chain[[1]] + (# Sign[#[[1]]]) &[RandomPoint[Sphere[ConstantArray[0., dim], r0]]];

chaincounter = 2;
While[chaincounter < maxchainlength,
  nf = Nearest[chain[[1 ;; chaincounter - 1]] -> Automatic];
  ncollisions = 1;
  iter = 0;
  While[ncollisions > 0 && iter < maxtrials,
   reservoircounter++;
   iter++;
   If[reservoircounter > reservoirLength,
    reservoircounter = 1;
    reservoir = getReservoir[reservoirLength, r0];
    ];
   xnew = x + reservoir[[reservoircounter]];
   ncollisions = If[xnew[[1]] >= 0., Length[nf[xnew, {∞, rc}]], 1];
   ];
  If[iter >= maxtrials,
   Break[];
   ,
   chain[[chaincounter + 1]] = x = xnew;
   chaincounter++;
   ]
  ];

And some visualization:

Graphics[{
  Line[chain[[1 ;; chaincounter]]], Blue, 
  Point[chain[[1 ;; chaincounter]]],
  Red, Opacity[0.15], Disk[#, rc/2] & /@ chain[[1 ;; chaincounter]],
  Darker@Green, PointSize[0.02], Opacity[1], Point[chain[[1]]], 
  Darker@Red, Point[chain[[chaincounter]]]}
 ]

enter image description here

The algorithm can easily be made working for dimension 3 by setting dim = 3. This way, one can obtain something like this:

Graphics3D[{
  Orange, Specularity[White, 30], 
  Sphere[chain[[1 ;; chaincounter]], rc/2],
  Red, Darker@Green, Sphere[chain[[1]], rc], Darker@Red, 
  Sphere[chain[[chaincounter]], rc],
  }, Lighting -> "Neutral"]

enter image description here

$\endgroup$
  • $\begingroup$ Awesome! Thanks I will keep the question unanswered for now in case anyone else has a more elegant solution $\endgroup$ – pukkandan Apr 12 '18 at 7:37
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Apr 12 '18 at 7:58

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.