failure to fit persistence length across the whole PDB trajectory

[Adolfo Poma]

Dear all,

I have a PDB trajectory (502 Frances and one liner chain with 198 beads connected) and each single monomer is represented by one bead (name CA). The script below seems to read and calculate the persistence length for the first 29 frames and then suddenly it prints out the following message:

Taken from attached message)lp.txt:

The persistence length is 11.7583 A

Traceback (most recent call last):

  File "lp.py", line 22, in <module>

    p.run()

  File "/Users/apoma/miniconda2/envs/mdaenv/lib/python3.6/site-packages/MDAnalysis/analysis/base.py", line 200, in run

    self._conclude()

  File "/Users/apoma/miniconda2/envs/mdaenv/lib/python3.6/site-packages/MDAnalysis/analysis/polymer.py", line 231, in _conclude

    self._perform_fit()

  File "/Users/apoma/miniconda2/envs/mdaenv/lib/python3.6/site-packages/MDAnalysis/analysis/polymer.py", line 254, in _perform_fit

    self.lp = fit_exponential_decay(self.x, self.results)

  File "/Users/apoma/miniconda2/envs/mdaenv/lib/python3.6/site-packages/MDAnalysis/analysis/polymer.py", line 307, in fit_exponential_decay

    a = scipy.optimize.curve_fit(expfunc, x, y)[0][0]

  File "/Users/apoma/miniconda2/envs/mdaenv/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 756, in curve_fit

    raise RuntimeError("Optimal parameters not found: " + errmsg)

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 400.

could somebody help me to fix this issue.

best regards,

Adolfo

------ the persistence length script -----------------------------

from __future__ import print_function

import MDAnalysis as mda

from MDAnalysis.analysis.polymer import PersistenceLength

import matplotlib.pyplot as plt

#--- select the data to read from ---#

u = mda.Universe('chainA.pdb')

print('We have a universe: {}'.format(u))

print('We have {} chains'.format(len(u.residues)))

print('Our atom types are: {}'.format(set(u.atoms.types)))

print('Trajectory frames are: {}'.format(len(u.trajectory)))

#--- select desire atoms/index ---#

ags=[u.atoms.select_atoms('name CA')]

print(ags)

#--- calculate Lp ---#

with open('output.dat', 'w') as out:

    for i in range(len(u.trajectory)):

        p = PersistenceLength(ags,start=i, stop=i+1)

        p.run()

        p.perform_fit()

        print("The persistence length is {:.4f} A".format(p.lp))

        out.write('frame {} lp {}'.format(i, p.lp))

------ the persistence length script -----------------------------

[Oliver Beckstein]

Hi Adolfo,

The calculation depends on having sufficient data to perform a fit. If there aren't enough data or the data are weird in some way it might fail. I think you are calculating the persistence length for single frames of the trajectory. That would explain why you might not have enough statistics. You can either increase your window over which you average or simply ignore the frames for which you don't have enough data:

p = PersistenceLength(ags,start=i, stop=i+1)

try:

  p.run()

except RuntimeError:

  continue  # start a new iteration of the loop

# do something with p

or something similar. However, I am not sure how much I would trust the final time series of persistence lengths if there aren't enough data points.

Oliver