polymer persistence length example fails

[orbeckst]

Expected behavior

When running the example in the docs for calculating the persistence length it should complete without errors.

Actual behavior

I copy and paste (0.20.1 and develop) and get the following error in

sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-2-07a190ddcd6d> in <module>
      7 backbones = [chain.select_atoms('not name O* H') for chain in chains]
      8 # sort the chains, removing any non-backbone atoms
----> 9 sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]
     10 lp = polymer.PersistenceLength(sorted_backbones)
     11 # Run the analysis, this will average over all polymer chains

<ipython-input-2-07a190ddcd6d> in <listcomp>(.0)
      7 backbones = [chain.select_atoms('not name O* H') for chain in chains]
      8 # sort the chains, removing any non-backbone atoms
----> 9 sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]
     10 lp = polymer.PersistenceLength(sorted_backbones)
     11 # Run the analysis, this will average over all polymer chains

~/Projects/Methods/MDAnalysis/mdanalysis/package/MDAnalysis/core/groups.py in check_args(*args, **kwargs)
   4091                                 funcname=func.__name__,
   4092                                 attrs=', '.join(missing)))
-> 4093             return func(*args, **kwargs)
   4094         return check_args
   4095     return require_dec

~/Projects/Methods/MDAnalysis/mdanalysis/package/MDAnalysis/analysis/polymer.py in sort_backbone(backbone)
    112             "Backbone is not linear.  "
    113             "The following atoms have more than two bonds in backbone: {}."
--> 114             "".format(','.join(str(a) for a in branches)))
    115 
    116     caps = [atom for atom in backbone

ValueError: Backbone is not linear.  The following atoms have more than two bonds in backbone: <Atom 7: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 11: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 24: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 28: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 41: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 45: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 58: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 62: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 75: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 79: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 92: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 96: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 109: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 113: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 126: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 130: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 143: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 147: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 160: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 164: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 177: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 181: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 194: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 198: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 211: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 215: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 228: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 232: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 245: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 249: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 262: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 266: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 279: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 283: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 296: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 300: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 313: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 317: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 330: N of type 5 of resname PA, resid 1 and segid PA>,<Atom 334: N of type 5 of resname PA, resid 1 and segid PA>.

Code to reproduce the behavior

from MDAnalysis.tests.datafiles import TRZ_psf, TRZ
import MDAnalysis as mda

from MDAnalysis.analysis import polymer
u = mda.Universe(TRZ_psf, TRZ)
# this system is a pure polymer melt of polyamide,
# so we can select the chains by using the .fragments attribute
chains = u.atoms.fragments
# select only the backbone atoms for each chain
backbones = [chain.select_atoms('not name O* H') for chain in chains]
# sort the chains, removing any non-backbone atoms
sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]
lp = polymer.PersistenceLength(sorted_backbones)
# Run the analysis, this will average over all polymer chains
# and all timesteps in trajectory
lp = lp.run()
print('The persistence length is: {}'.format(lp.pl))
# always check the visualisation of this:
lp.plot()

Specifically, the step

sorted_backbones = [polymer.sort_backbone(bb) for bb in backbones]

fails.

Currently version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)"): 0.20.1 and latest develop
• Which version of Python (python -V)? 3.7
• Which operating system? macOS

[lilyminium]

The is missing after the H in `'not name O H'`. This version continues to work.

Edit: oh, and lp.pl should be lp.lp. I guess that is the cause of the mailing list issue

[orbeckst]

Thanks, that would be an easy GSOC PR…

On Mar 4, 2020, at 6:16 PM, Lily Wang notifications@github.com wrote:

The is missing after the H in 'not name O H'. This version continues to work. https://www.mdanalysis.org/UserGuide/examples/analysis/polymers_and_membranes/polymer.html

[zemanj]

Edit: oh, and lp.pl should be lp.lp. I guess that is the cause of the mailing list issue

I think this is a neat example for a rather suboptimal choice of variable names. IMHO we should replace lp and lb by persistence_length and bond_length and deprecate the old names. Also, I don't understand why results is not a dict containing all results. (If it is just the autocorrelation, it should be named accordingly.)

[orbeckst]

I removed @richardjgowers so that GSoC students don't feel that the issue is taken.

[ss62171]

Can i work on this?

[orbeckst]

By all means, put up a PR. Don't forget to reference this issue.

[siddharthjain1611]

I'd love to commit to this issue.

[ss62171]

Hey @siddharthjain1611 currently I'm working on this and almost done.

[siddharthjain1611]

Okay @ss62171 I'd wait for your patch before I move further with this.

[ss62171]

@orbeckst can you please look at this PR and suggest if some changes required.

[orbeckst]

@ss62171 I reviewed PR #2589 .

@siddharthjain1611 thank you for being polite and letting @ss62171 finish their work. As a general rule for how we have been handling things: we look at the first PR that references the issue. If the PR stalls then we close it and look at the next one.

[siddharthjain1611]

Still wondering, is this issue open?

[orbeckst]

No, see PR #2589