Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to calculate the persistence length (Lp) of a polymer chain using Mathematica. The first and simplest method is to compare the EuclideanDistance between two points on the polymer and arc length separating the points. This works for me. I'll describe how I implemented this now and then give the other method based on tangent-tangent angle correlations which should also work but I'm having trouble with. Hopefully someone might know about this stuff and be able to help. This article describes what I'm trying to do. I've put some more references and a helpful figure at the end. Two of the linked articles at the end are free/open access and I can provide copies of the other linked articles if people are interested and cannot retrieve them.

You can get a sample dataset here, units are in microns but it doesn't really matter:

data = Uncompress@Import["http://pastebin.com/raw.php?i=60aazeuc"]

In the literature traced curves from pixel-based images are often smoothed using a weighted moving average

smoothdata = MovingAverage[data, {1, 2, 4, 2, 1}];

To calculate the arc length lets interpolate along the curve with a cubic spline to a parametrized curve (starting at 0).

f = Interpolation[Transpose[{Range[0, Length[smoothdata] - 1], smoothdata}], 
                  Method -> "Spline", InterpolationOrder -> 3];

To find the length along the curve I based it on @Szabolcs answer.

v = FunctionInterpolation[Norm[f'[t]], {t, 0, Length[smoothdata] - 1}]
dist = Derivative[-1][v]
totallength = dist[Length[smoothdata] - 1];
invdist = FunctionInterpolation[InverseFunction[dist][d], 
                               {d, 0, totallength}]

The invdist function returns the parameter along the curve for a given length.

Now I made a function that takes a window of arc length x (in microns) and translates it along the curve (in 5 nm steps) and then returns the mean of the square of the EuclideanDistance between the end points of all the subsets of the curve. This is repeated for a range of arc lengths and then the result

distance[x_] := 
  Mean[
    Table[
        EuclideanDistance[f[invdist[t]], f[invdist[t + x]]]^2, 
        {t, 0, totallength - x, 0.005}
    ]
  ];

Now we can run the function:

tofit = ParallelTable[{x, distance[x]}, {x, 0, totallength, 0.1}]

and fit the result (dropping the first (0,0) point) to the model from the literature:

fit = NonlinearModelFit[
        Rest@tofit, {4 Lp s (1 - (2 Lp /s) (1 - Exp[-s/(2 Lp)]))}, {Lp}, s,
        MaxIterations -> 1000, Method -> "NMinimize"
      ];

Ta da! It works. Lets plot the data and the fit.

Show[
  ListPlot[tofit, Joined -> True, PlotRange -> All], 
  Plot[fit[x], {x, 0, totallength}, PlotStyle -> {Red, Dashed}]
]

Mathematica graphic

In this case Lp is 0.79 microns, sounds reasonable.

The other method to measure the persistence length is based on the angle between tangents on the curve separated by a certain arc length. See: this graphic stolen from the linked paper. Using the cubic spline fit its easy to calculate the tangent vectors at a point. This next function similarly translates a window of arc length x along the curve except this time it finds the dot product (cosine of the angle) of the unit tangent vectors at the endpoints of the window and returns the mean of all the dot products (cosines of the angles).

cosinecorrelation[x_] := 
  Mean[
    Table[
          Dot[Normalize[f'[invdist[t]]], Normalize[f'[invdist[t + x]]]], 
         {t, 0, totallength - x, 0.005}
    ]
  ];

This should result in an exponential decay <cos(theta)> = Exp[-x/Lp], However, I get basically a linear function, also the values drop below zero. As far as I can tell that shouldn't happen... Does something look wrong with what I did?

ListPlot[
   ParallelTable[{x, cosinecorrelation[x]}, 
                 {x, 0, totallength, 0.1}
   ], 
   Joined -> True, PlotRange -> All
]

Mathematica graphic

The next method is to look at the mean square angle (instead of the cosine of the angle), <theta^2> = -x/Lp which should be a linear function of arc length, x. However, (warning, potentially stupid question follows), when I use ArcCos I always get a positive angle (same with VectorAngle), even when the dot product is negative. This must have to do with how Mathematica defines ArcCos, how, can I make it so negative dot products give the equivalent negative angle. Here is my broken function, based on the MMA ArcCos

anglesquaredcorrelation[x_] := 
  Mean[
    Table[
      (ArcCos[Dot[Normalize[f'[invdist[t]]], Normalize[f'[invdist[t + x]]]]])^2, 
      {t, 0, totallength - x, 0.005}
    ]
  ];

An alternative implementation, based on VectorAngle, which confusingly doesn't seem to give the same values for larger x.

vectorangle[x_] := 
 Mean[
   Table[
      VectorAngle[f'[invdist[t]], f'[invdist[t + x]]]^2, 
      {t, 0, totallength - x, 0.005}
 ];

Any help would be greatly appreciated! For interested parties much of the theory was presented in this 1979 paper based on earlier work by Landau and Lifshitz (Landau, L. D. & Lifshitz, E. M. (1958) Statistical Physics, Course of Theoretical Physics, Vol. 5, Pergamon, London.)

Here is a helpful set of plots and equations from another paper:

Dropbox graphic

Open access article this time: PDF, and another open access paper PDF.

share|improve this question
1  
Sqrt[#.#] &[...] could be written as Norm[...] –  Sjoerd C. de Vries Aug 1 '12 at 19:26
    
The ArcCos doc page says: For real z between -1 and 1, the results are always in the range 0 to pi. I believe this is the standard convention. (see also Wikipedia) –  Sjoerd C. de Vries Aug 1 '12 at 19:42
    
@SjoerdC.deVries, I added an image borrowed from a paper that illustrates the bending angle and how they have negative angles. –  s0rce Aug 1 '12 at 20:04
2  
Your cosinecorrelation looks perfectly reasonable, judging from the data. To see an exponential decay you would either need a much wigglier chain, or to average over many data sets. Your vectorangle and anglesquaredcorrelation gave exactly the same results when I tried them. Try re-running your code with a fresh kernel. If you want a signed angle between tangents, I suggest instead of ArcCos[a.b] you use ArcTan@@a - ArcTan@@b –  Simon Woods Aug 2 '12 at 12:10
    
@s0rce is your data 2D $(X,Y)$, or 3D $(\phi,\theta)$? In my opinion MovingAverage introduces an artificial correlation on your data, which is actually what you need to measure. –  rhermans Feb 20 at 18:58

1 Answer 1

up vote 1 down vote accepted

You say your data is a list of $(X,Y)$ so it looks like this:

ListLinePlot[data, Frame -> True, FrameLabel -> {"X", "Y"}]

Mathematica graphics

Far from been a random coil where the definition of persistent length makes sense. Let see, we want to calculate the average $\cos(\theta) = \hat{v_1} \cdot \hat{v_1}$, as a function of the contour distance, so I use Dot[v1, v2]/(Norm[v1] Norm[v2]) versus Total[Norm /@ Differences[data[[i ;; i + j]] for all possible combinations of $\vec{v_1}$, $\vec{v_2}$

cosThetaList = Flatten[
   Table[
    Table[
     With[
      {
       v1 = data[[i]] - data[[i + 1]]
       , v2 = data[[i + j]] - data[[i + j + 1]] 
       }, 
      N[{Total[Norm /@ Differences[data[[i ;; i + j]]]], Dot[v1, v2]/(
        Norm[v1] Norm[v2])}]]
     , {j, 1, Length[data] - i - 1}]
    , {i, 1, Length[data] - 2}
    ], 1];

Now for the average value

avgCosTheta = MovingAverage[Sort[cosThetaList, #1[[1]] > #2[[1]] &], 100];

ListLinePlot[
 avgCosTheta
 , PlotRange -> All
 , Frame -> True
 , FrameLabel -> {"Contour length", "<Cos[\[Theta]]>"}
 ]

Mathematica graphics

Which resembles your result. In my opinion, is that the randomness implicit in the model that defines the persistent length is not present in your data, probably because the length of your polymer is too short or the physics is different from a random walk with partial memory.

Probably a way to define $P$ regardless the non-exponential shape of the curve, would be to find the value of the contour distance at which the cosine crosses $\exp(-1)$, in this case approximately $P\approx1.07$ in whatever the units your data was on.

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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