Hydrogen bond analysis exits with AttributeError: 'list' object has no attribute 'select_atoms'

[vivecalindahl]

Expected behaviour

Get an hbond analysis (even though distance calculations are not taking care of periodicity).

Actual behaviour

Traceback (most recent call last):
File "./reproducehbondissue.py", line 45, in
distance=hbonddist, angle=hbondangle)
File "/usr/local/lib/python2.7/dist-packages/MDAnalysis/analysis/hbonds/hbond_analysis.py", line 612, in init
self._update_selection_2()
File "/usr/local/lib/python2.7/dist-packages/MDAnalysis/analysis/hbonds/hbond_analysis.py", line 798, in _update_selection_2
self._s2_acceptors = self._s2.select_atoms(
AttributeError: 'list' object has no attribute 'select_atoms'

Code to reproduce the behaviour

reproduce.zip

script in reproduce.zip reads gro frame.

Currently version of MDAnalysis:

(run python -c "import MDAnalysis as mda; print(mda.__version__)")
0.15.0

[orbeckst]

I put up your code as https://gist.github.com/orbeckst/eaad2d7cf063901d9b48d5e357ac3e61 – that is easier to browse for people than looking into the zip file. (I didn't want to attach the GRO file t39000.gro in the zip because of its size.)

[orbeckst]

@vivecalindahl I played around with your code a little bit (using 0.16.0-dev). The main problem is that your selection

u.select_atoms(' and '.join(['segid ' + segid, 'resid 197', 'name ND2']))

yields an empty AtomGroup

 <AtomGroup with 0 atoms>

When I looked at the universe itself

>>> u.segments.segids
array(['monomer4'], dtype=object)

but you are trying to select

>>> print(segid)
'monomer1

for the donor group, which is not a segment name in the universe.

HydrogenBondAnalysis should fail with a better error message when one group is empty.

In order for your code to pass through 0.16.0-dev I had to change one thing in segment_monomers(), namely monomer.segids = segid to monomer.segments.segids = segid.

def segment_monomers(u):
    monomers = select_monomers(u)
    segids = ["monomer" + str(i+1) for i in range(len(monomers))]

    for monomer, segid in zip(monomers, segids):
        monomer.segments.segids = segid

    return segids

The problem really lies here, though. Even though the monomer atomgroups consist of different atoms, they still all belong to the same segment (namely, "SYSTEM"). What happens in the loop is that you first name "SYSTEM" as "monomer1". Then you rename "monomer1" to "monomer2"... and finally, it is called "monomer4". But all 4 monomers still belong to the same segment.

I admit that I actually don't know how to split AtomGroups into multiple segments with the new topology system – maybe one of the other @MDAnalysis/coredevs can chime in?

[jbarnoud]

On 04/04/2017 09:48 AM, Oliver Beckstein wrote: @vivecalindahl <https://github.com/vivecalindahl> I played around with your code a little bit (using 0.16.0-dev). The main problem is that your selection u.select_atoms(' and '.join(['segid ' + segid,'resid 197','name ND2'])) yields an empty AtomGroup |<AtomGroup with 0 atoms> | When I looked at the universe itself >>> u.segments.segids array(['monomer4'],dtype=object) but you are trying to select >>> print(segid) 'monomer1 for the donor group, which is not a segment name in the universe. |HydrogenBondAnalysis| *should fail with a better error message when one group is empty*. In order for your code to pass through 0.16.0-dev I had to change one thing in |segment_monomers()|, namely |monomer.segids = segid| to |monomer.segments.segids = segid|. |def segment_monomers(u): monomers = select_monomers(u) segids = ["monomer" + str(i+1) for i in range(len(monomers))] for monomer, segid in zip(monomers, segids): monomer.segments.segids = segid return segids | The problem really lies here, though. Even though the |monomer| atomgroups consist of different atoms, they still all belong to the /same segment/ (namely, "SYSTEM"). What happens in the loop is that you first name "SYSTEM" as "monomer1". Then you rename "monomer1" to "monomer2"... and finally, it is called "monomer4". But all 4 monomers /still belong to the same segment/. I admit that I actually don't know how to split AtomGroups into multiple segments with the new topology system – maybe one of the other @MDAnalysis/coredevs <https://github.com/orgs/MDAnalysis/teams/coredevs> can chime in? — You are receiving this because you are on a team that was mentioned. Reply to this email directly, view it on GitHub <#1268 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABUWuvjmjQmS4vSOyjQGXETSCm9XxWw_ks5rsfXZgaJpZM4MrxWJ>.

I am not aware of an easy way to re-segment a topology. I did it once by creating a new Topology object from the one in the Universe. Fun to do, but convoluted. This is an issue as gomacs files do not give any meaningful segments, especially when working with multimeric proteins.

[richardjgowers]

@orbeckst in 0.16 this is the way to make segments:

import MDAnalysis as mda
from MDAnalysisTests.datafiles import GRO

u = mda.Universe(GRO)

new_seg = u.add_Segment(segid='new')

u.residues[:10].segments = new_seg

This won't work in 0.15 though. I think this is a nice example of where selection strings aren't a good method for defining AtomGroups for analysis modules.

[orbeckst]

@vivecalindahl the segmenting issue is separate from using HydrogneBondAnalysis. If you were to use the selections that you created for the segmentation then it should work. The following code worked for me:

#! /usr/bin/env python2.7

import MDAnalysis as mda
import MDAnalysis.analysis.hbonds

def select_monomers_str(u, nmer=4):
    protein = u.select_atoms("protein")
    return ["protein and bynum {0}:{1}".format(i1, i2)
            for i1 in range(1, protein.n_atoms, protein.n_atoms/nmer)
            for i2 in [i1+protein.n_atoms/nmer - 1]]


if __name__ == "__main__":
    # input
    grofile = "t39000.gro"

    # Create the md universe
    u = mda.Universe(grofile)

    # Analyze monomers separately
    # Need to use selection string for each monomer:
    monomer_selections = select_monomers_str(u)

    hbonddist = 3.0 # distance in Angstrom between donor H and acceptor atom.
    hbondangle = 180-35 # deg (different angle def in the mda hbond tool than for gmx)

    # just look at one monomer for right now;
    # To loop: for monomer_selection in monomer_selections:
    imonomer = 0

    sel1 = "{0} and resid 197 and name ND2".format(monomer_selections[imonomer])
    hbondana = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, \
                    selection1=sel1,
                    selection2='resname SOL and name OW',
                    selection1_type='donor',
                    distance=hbonddist, angle=hbondangle,
                    filter_first=False)
    print 'running hbond ana for monomer {}'.format(imonomer)
    hbondana.run()

I then get as output

>>> hbondana.timeseries
[[[14165,
   80196,
   14164,
   80195,
   'ASN197:HD21',
   'SOL10202:OW',
   1.7456797,
   149.12985]]]

NOTE: I also needed to set filter_first=False.

If I use the default filter_first=True I get the error

AttributeError: 'list' object has no attribute 'select_atoms'
> /Volumes/Data/oliver/Biop/Library/python/mdanalysis/package/MDAnalysis/analysis/hbonds/hbond_analysis.py(819)_update_selection_2()
    817         self._s2_acceptors = {}
    818         if self.selection1_type in ('donor', 'both'):
--> 819             self._s2_acceptors = self._s2.select_atoms(
    820                 'name {0}'.format(' '.join(self.acceptors)))
    821             self.logger_debug("Selection 2 acceptors: {0:d}".format(len(self._s2_acceptors)))

For reasons that I do not understand, the filtering does not seem to be finding that one water molecule that is actually close by (according to the result above, SOL10202:OW is within 1.7 Å of the donor) --- you can see this in the log file (mda.start_logging()):

2017-04-04 02:12:21,108 MDAnalysis.analysis.hbonds DEBUG    Size of selection 1: 1 atoms
2017-04-04 02:12:21,109 MDAnalysis.analysis.hbonds DEBUG    Selection 1 donors: 1
2017-04-04 02:12:21,109 MDAnalysis.analysis.hbonds DEBUG    Selection 1 donor hydrogens: 1
2017-04-04 02:12:21,122 MDAnalysis.analysis.hbonds DEBUG    Size of selection 2 before filtering: 20920 atoms
2017-04-04 02:12:21,144 MDAnalysis.analysis.hbonds DEBUG    Size of selection 2: 0 atoms
2017-04-04 02:12:21,144 MDAnalysis.analysis.hbonds WARNING  Selection 2 "resname SOL and name OW" did not select any atoms.

Maybe something with the NeighborSearch in hbond_analysis.py L809 is not working well???

[vivecalindahl]

@orbeckst What you find in the developer code is different from what I find in 0.15. It could be that I've misunderstood how the segments work... but I get

print u.segments.segids
['monomer1' 'monomer2' 'monomer3' 'monomer4' 'SYSTEM']

and

sel=u.select_atoms(' and '.join(['segid ' + segids[0], 'resid 197', 'name ND2']))
<AtomGroup with 1 atoms>

so the issue is different than what you described. How am I using segments incorrectly?

I tested adding filter_first=False which does remove the original error. hbondana.timeseries now returns an empty list (for the first monomer/segment). This could be the correct answer but it is not the same as the result you posted. Actually, I get the results you posted for the fourth monomer.

[kain88-de]

We have a almost complete rewrite of the core of MDAnalysis in the development version. You can read more here. This is why segments work slightly different now.

[orbeckst]

Am Apr 4, 2017 um 08:58 schrieb vivecalindahl ***@***.***>: How am I using segments incorrectly?

You're not, you're doing the right thing with 0.15.0. Sorry for the confusion (but good to keep in mind for upcoming 0.16.0).

I tested adding filter_first=False which does remove the original error. hbondana.timeseries now returns an empty list (for the first monomer/segment). This could be the correct answer but it is not the same as the result you posted.

Because mine effectively looked at monomer4, as you found.

Actually, I get the results you posted for the fourth monomer.

The issue that you report is a bug in the code. It should handle the case of no atoms within the 3*distance gracefully. This needs to be fixed. Thanks for reporting. In the meantime please try without the filter_first optimization --- sorry.

[vivecalindahl]

@orbeckst Ok, thanks for looking into it.

[kain88-de]

If I use the default filter_first=True I get the error

This is already fixed in develop. I don't get an error any longer. I still do get it with 0.16.0 though.

@orbeckst can you confirm and close this if the issue is solved or provide more information about the error.

[orbeckst]

@vivecalindahl 's case is now running for me; I am using reproducehbondissue-0.16.1.py . I am detecting 0, 7, 7, 7 (monomer 1-4) waters in the 2nd selection group when using filter_first=True. This agrees with the visualization in VMD (I did a selection like same residue as (name OW and within 9 of (resid 197 and name ND2)); the residue 197 in each of the four monomers are shown and the ND2 atom is highlighted as blue VDW sphere):

There are three h-bonds detected

2017-05-21 00:07:30,008 MDAnalysis.analysis.hbonds DEBUG    S1-D: 6780 <-> S2-A: 109349 1.981566 A, 153.637421 DEG
2017-05-21 00:07:30,256 MDAnalysis.analysis.hbonds DEBUG    S1-D: 10472 <-> S2-A: 105461 2.199296 A, 163.300568 DEG
2017-05-21 00:07:30,501 MDAnalysis.analysis.hbonds DEBUG    S1-D: 14164 <-> S2-A: 80195 1.745680 A, 170.232483 DEG

I am almost certain that this got fixed in d5bc8f5 (PR #1325 for #1327 ).

I am closing this for now. We can reopen if @vivecalindahl says that there is still a problem with the latest develop version.

[orbeckst]

@kain88-de many thanks for checking... I had totally forgotten that @xiki-tempula had encountered a similar problem and fixed #1327 !

[omkantnirala]

Hi everyone,
i am beginner in MDAnalysis , can anyone help regarding this error.

TypeError Traceback (most recent call last)
in ()
8 H = u.select_atoms('byres(name OW HW1 HW2 and around 6.3 name SOL)')
9 hb_ac = hbonds.HydrogenBondAutoCorrel(u,acceptors=A, hydrogens=H, donors=D ,number_time =2, nruns =10,
---> 10 nsample =10, dist_crit=2.4,angle_crit=130.0,pbc=True)
11 hb_ac.run()
12 hb_ac.solv()

TypeError: init() got an unexpected keyword argument 'nsample'

Thanks.

[orbeckst]

Hi, Please don’t use the issue tracker for questions: subscribe to the mailing list https://groups.google.com/forum/#!forum/mdnalysis-discussion and ask there. This is where we answer questions. If you look at the docs for HydrogenBondAutoCorrel <https://www.mdanalysis.org/docs/documentation_pages/analysis/hbond_autocorrel.html#MDAnalysis.analysis.hbonds.hbond_autocorrel.HydrogenBondAutoCorrel> then you see that the name of the argument is nsamples and not “nsample”. Oliver

…

On Apr 16, 2019, at 2:38 AM, omkantnirala ***@***.***> wrote: Hi everyone, i am beginner in MDAnalysis , can anyone help regarding this error. TypeError Traceback (most recent call last) in () 8 H = u.select_atoms('byres(name OW HW1 HW2 and around 6.3 name SOL)') 9 hb_ac = hbonds.HydrogenBondAutoCorrel(u,acceptors=A, hydrogens=H, donors=D ,number_time =2, nruns =10, ---> 10 nsample =10, dist_crit=2.4,angle_crit=130.0,pbc=True) 11 hb_ac.run() 12 hb_ac.solv() TypeError: init() got an unexpected keyword argument 'nsample' Thanks.