Persistence length of a polymer chain

Question

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}]
]

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
]

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:

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

Answer 1

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

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

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]]>"}
 ]

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.