Source code for MDAnalysis.analysis.polymer
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- http://www.MDAnalysis.org
# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein
# and contributors (see AUTHORS for the full list)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""
Polymer analysis --- :mod:`MDAnalysis.analysis.polymer`
=======================================================
:Author: Richard J. Gowers
:Year: 2015
:Copyright: GNU Public License v3
This module contains various commonly used tools in analysing polymers.
"""
from six.moves import range
import numpy as np
import logging
from .. import NoDataError
from ..lib.util import blocks_of
from ..lib.distances import calc_bonds
from .base import AnalysisBase
logger = logging.getLogger(__name__)
[docs]class PersistenceLength(AnalysisBase):
r"""Calculate the persistence length for polymer chains
The persistence length is the length at which two points on the polymer
chain become decorrelated.
Notes
-----
This analysis requires that the trajectory supports indexing
.. versionadded:: 0.13.0
"""
def __init__(self, atomgroups,
start=None, stop=None, step=None):
"""Calculate the persistence length for polymer chains
Parameters
----------
atomgroups : list
List of atomgroups. Each atomgroup should represent a single
polymer chain, ordered in the correct order.
start : int, optional
First frame of trajectory to analyse, Default: 0
stop : int, optional
Last frame of trajectory to analyse, Default: -1
step : int, optional
Step between frames to analyse, Default: 1
"""
self._atomgroups = atomgroups
# Check that all chains are the same length
lens = [len(ag) for ag in atomgroups]
chainlength = len(atomgroups[0])
if not all( l == chainlength for l in lens):
raise ValueError("Not all AtomGroups were the same size")
self._setup_frames(atomgroups[0].universe.trajectory,
start, stop, step)
self._results = np.zeros(chainlength - 1, dtype=np.float32)
def _single_frame(self):
# could optimise this by writing a "self dot array"
# we're only using the upper triangle of np.inner
# function would accept a bunch of coordinates and spit out the
# decorrel for that
n = len(self._atomgroups[0])
for chain in self._atomgroups:
# Vector from each atom to next
vecs = chain.positions[1:] - chain.positions[:-1]
# Normalised to unit vectors
vecs /= np.sqrt((vecs * vecs).sum(axis=1))[:, None]
inner_pr = np.inner(vecs, vecs)
for i in range(n-1):
self._results[:(n-1)-i] += inner_pr[i, i:]
def _conclude(self):
n = len(self._atomgroups[0])
norm = np.linspace(n - 1, 1, n - 1)
norm *= len(self._atomgroups) * self.nframes
self.results = self._results / norm
self._calc_bond_length()
def _calc_bond_length(self):
"""calculate average bond length"""
bs = []
for ag in self._atomgroups:
pos = ag.positions
b = calc_bonds(pos[:-1], pos[1:]).mean()
bs.append(b)
self.lb = np.mean(bs)
[docs] def plot(self):
"""Oooh fancy"""
import matplotlib.pyplot as plt
plt.ylabel('C(x)')
plt.xlabel('x')
plt.xlim([0.0, 40 * self.lb])
plt.plot(self.x, self.results, 'ro')
plt.plot(self.x, self.fit)
plt.show()
[docs]def fit_exponential_decay(x, y):
r"""Fit a function to an exponential decay
.. math:: y = \exp(-x/a)
Parameters
----------
x, y : array_like
The two arrays of data
Returns
-------
a : float
The coefficient *a* for this decay
Notes
-----
This function assumes that data starts at 1.0 and decays to 0.0
Requires scipy
"""
from scipy.optimize import curve_fit
def expfunc(x, a):
return np.exp(-x/a)
a = curve_fit(expfunc, x, y)[0][0]
return a