problem in MDAnalysis.analysis.hbonds

108 views
Skip to first unread message

Natalia Kanaan Izquierdo

unread,
Sep 19, 2013, 5:46:11 AM9/19/13
to mdnalysis-...@googlegroups.com
Hi all,
I'm using the module MDAnalysis.analysis.hbonds and it works  but, depending on the selection, it doesn't.
For instance, if my selection is:
h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname HSP and resid 151', 'resname TIP3', distance=3.0, angle=120.0)
it works, but if it is:
h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname THY and resid 11', 'resname TIP3', distance=3.0, angle=120.0)
I checked and it understands the selection,  but there is a problem when it is searching for the neighbours.
Does anyone have a clue on what the problem can be?
Thanks a lot in advance,
Natal

Oliver Beckstein

unread,
Sep 19, 2013, 1:43:50 PM9/19/13
to mdnalysis-...@googlegroups.com
Hi,

On 19 Sep, 2013, at 02:46, Natalia Kanaan Izquierdo wrote:

> Hi all,
> I'm using the module MDAnalysis.analysis.hbonds and it works but, depending on the selection, it doesn't.
> For instance, if my selection is:
> h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname HSP and resid 151', 'resname TIP3', distance=3.0, angle=120.0)
> it works, but if it is:
> h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname THY and resid 11', 'resname TIP3', distance=3.0, angle=120.0)

What kind of residue is "THY"?

What are the atoms called? If this is non-standard then you might have to explicitly say which atoms are considered donors or acceptors. See http://pythonhosted.org/MDAnalysis/documentation_pages/analysis/hbonds.html#MDAnalysis.analysis.hbonds.HydrogenBondAnalysis

If you think that your atoms names will be of general use to the community then we can add them to the code.

Oliver

> I checked and it understands the selection, but there is a problem when it is searching for the neighbours.
> Does anyone have a clue on what the problem can be?
> Thanks a lot in advance,
> Natal
>
> --
> You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
> To post to this group, send email to mdnalysis-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/mdnalysis-discussion.
> For more options, visit https://groups.google.com/groups/opt_out.

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst * orbe...@gmail.com

Samuel Jayakanthan

unread,
Mar 20, 2016, 2:27:55 AM3/20/16
to MDnalysis discussion

Dear MDanalysis users,
 I have been using Hydrogen bond analysis module to extract Hbond information for my trajectories. While, I am able to extract several details (like donor acceptor idx, angles, distances.) I do see the digit 3 from the donor or acceptor name TIP3 (for waters) get prefixed onto the folowing column (please see below). The resnames in the PDB and PSF files are "TIP3" as opposed to :TIP 3". I also scrolled down the PSF/PDB files all the way to the bottom to make sure that the column which contains "TIP3" and the adjacent column does have a space. Not sure why MDAnalysis would write that in the output. This becomes a problem because, when I am trying to extract the identities of the waters associated with my protein I get a 3 added to the actual number defining the TIP3.

any pointers to fix this might be very helpful.

thanks
Sam
#TIME donor_idx acceptor_idx donor_resnm donor_resid donor_atom acceptor_resnm acceptor_resid acceptor_atom distance angle
0.04888821 2 69595 GLN 375 HT1 TIP 321598 OH2 2.53951334953 120.064338684
0.04888821 3 70147 GLN 375 HT2 TIP 321782 OH2 1.75176846981 161.643386841
0.04888821 4 70102 GLN 375 HT3 TIP 321767 OH2 1.73476588726 164.393417358
0.04888821 16 69625 GLN 375 HE21 TIP 321608 OH2 2.39797735214 155.587142944
0.04888821 21 70102 ALA 376 HN TIP 321767 OH2 1.9973744154 161.092926025
0.04888821 31 70081 PHE 377 HN TIP 321760 OH2 2.76333880424 163.33543396
0.04888821 79 52942 ARG 379 HE TIP 316047 OH2 2.03155708313 139.219665527
0.04888821 79 54964 ARG 379 HE TIP 316721 OH2 2.73248505592 133.345199585
0.04888821 82 72793 ARG 379 HH11 TIP 322664 OH2 2.69918251038 135.309020996
0.04888821 83 54964 ARG 379 HH12 TIP 316721 OH2 1.89551150799 171.454376221
0.04888821 85 4414 ARG 379 HH21 TIP 3263 OH2 1.85218548775 151.634918213


0.04888821 4221 401 TIP 3133 H2 PHE 399 O 1.74142551422 176.970474243
0.04888821 4224 1065 TIP 3135 H2 LEU 444 O 1.88931870461 164.345016479
0.04888821 4226 646 TIP 3137 H1 GLN 416 OE1 1.89540481567 165.627960205
0.04888821 4229 2564 TIP 3139 H1 ASP 539 OD1 1.8857613802 178.408721924
0.04888821 4233 2403 TIP 3141 H2 ALA 528 O 2.23189997673 136.599227905
0.04888821 4239 2564 TIP 3145 H2 ASP 539 OD1 2.88705515862 137.274780273
0.04888821 4239 2565 TIP 3145 H2 ASP 539 OD2 1.63372635841 165.801696777
0.04888821 4241 3808 TIP 3147 H1 ASP 622 OD1 1.62873470783 163.808425903
0.04888821 4244 295 TIP 3149 H1 THR 393 OG1 1.9190107584 159.594497681

                                                                              
 
On Thursday, September 19, 2013 at 1:43:50 PM UTC-4, Oliver Beckstein wrote:
Hi,

On 19 Sep, 2013, at 02:46, Natalia Kanaan Izquierdo wrote:

> Hi all,
> I'm using the module MDAnalysis.analysis.hbonds and it works  but, depending on the selection, it doesn't.
> For instance, if my selection is:
> h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname HSP and resid 151', 'resname TIP3', distance=3.0, angle=120.0)
> it works, but if it is:
> h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname THY and resid 11', 'resname TIP3', distance=3.0, angle=120.0)

What kind of residue is "THY"?

What are the atoms called? If this is non-standard then you might have to explicitly say which atoms are considered donors or acceptors. See http://pythonhosted.org/MDAnalysis/documentation_pages/analysis/hbonds.html#MDAnalysis.analysis.hbonds.HydrogenBondAnalysis

If you think that your atoms names will be of general use to the community then we can add them to the code.

Oliver

> I checked and it understands the selection,  but there is a problem when it is searching for the neighbours.
> Does anyone have a clue on what the problem can be?
> Thanks a lot in advance,
> Natal
>
> --
> You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discussion+unsub...@googlegroups.com.

Oliver Beckstein

unread,
Mar 21, 2016, 4:08:51 AM3/21/16
to mdnalysis-...@googlegroups.com
Hi Sam,

I think before blaming h-bond analysis we need to ascertain that the residue names are properly read in the first place.

Which version of MDAnalysis do you run?

print(MDAnalysis.__version__)

How do you load your input, e.g.

u = mda.Universe(PSF, PDB)

?

What does the following show?

water = u.select_atoms("resname TIP*")
print(water.n_residues)
print(water.resnames[:10])
print(water.resids[:10])

Copy and paste the first 10 TIP3P lines from your PSF (something like `cat PSF | grep 'TIP' | head -n 10` might do it)

Oliver
> > To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
> > To post to this group, send email to mdnalysis-...@googlegroups.com.
> > Visit this group at http://groups.google.com/group/mdnalysis-discussion.
> > For more options, visit https://groups.google.com/groups/opt_out.
>
> --
> Oliver Beckstein * orbe...@gmx.net
> skype: orbeckst * orbe...@gmail.com
>
>
> --
> You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
> To post to this group, send email to mdnalysis-...@googlegroups.com.
> Visit this group at https://groups.google.com/group/mdnalysis-discussion.
> For more options, visit https://groups.google.com/d/optout.

Samuel Jayakanthan

unread,
Mar 22, 2016, 9:10:36 AM3/22/16
to MDnalysis discussion


Hi Oliver,
Thanks for getting back. Here are the details you were asking. Should I upgrade to a newer version?
thanks so much
Sam
       


u=Universe(PSF,DCD)

>>> print(MDAnalysis.__version__)
0.8.1



>>> print(water.residues[:10])
<ResidueGroup [<Residue 'TIP3', 1>, <Residue 'TIP3', 3>, <Residue 'TIP3', 5>, <Residue 'TIP3', 7>, <Residue 'TIP3', 9>, <Residue 'TIP3', 11>, <Residue 'TIP3', 13>, <Residue 'TIP3', 15>, <Residue 'TIP3', 17>, <Residue 'TIP3', 19>]>

4021 WATA     1        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4022 WATA     1        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4023 WATA     1        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4024 WATA     3        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4025 WATA     3        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4026 WATA     3        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4027 WATA     5        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4028 WATA     5        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4029 WATA     5        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4030 WATA     7        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4031 WATA     7        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02

4804 SOLV     1        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4805 SOLV     1        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4806 SOLV     1        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4807 SOLV     2        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4808 SOLV     2        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4809 SOLV     2        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4810 SOLV     3        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4811 SOLV     3        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4812 SOLV     3        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4813 SOLV     4        TIP3     OH2         3  -0.834000       15.9994           0   0.00000     -0.301140E-02
      4814 SOLV     4        TIP3     H1          1   0.417000       1.00800           0   0.00000     -0.301140E-02
      4815 SOLV     4        TIP3     H2          1   0.417000       1.00800           0   0.00000     -0.301140E-02
> > To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discussion+unsub...@googlegroups.com.
> > To post to this group, send email to mdnalysis-...@googlegroups.com.
> > Visit this group at http://groups.google.com/group/mdnalysis-discussion.
> > For more options, visit https://groups.google.com/groups/opt_out.
>
> --
> Oliver Beckstein * orbe...@gmx.net
> skype: orbeckst  * orbe...@gmail.com
>
>
> --
> You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discussion+unsub...@googlegroups.com.

Tyler Reddy

unread,
Mar 22, 2016, 11:08:25 AM3/22/16
to mdnalysis-...@googlegroups.com
Hi Sam,

That's a pretty old version of MDAnalysis. I'd suggest using 0.14 (see: http://www.mdanalysis.org/pages/installation_quick_start/)

Tyler

> > To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
> > To post to this group, send email to mdnalysis-...@googlegroups.com.
> > Visit this group at http://groups.google.com/group/mdnalysis-discussion.
> > For more options, visit https://groups.google.com/groups/opt_out.
>
> --
> Oliver Beckstein * orbe...@gmx.net
> skype: orbeckst  * orbe...@gmail.com
>
>
> --
> You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
> To post to this group, send email to mdnalysis-...@googlegroups.com.
> Visit this group at https://groups.google.com/group/mdnalysis-discussion.
> For more options, visit https://groups.google.com/d/optout.

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst  * orbe...@gmail.com

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

Oliver Beckstein

unread,
Mar 22, 2016, 12:52:25 PM3/22/16
to mdnalysis-...@googlegroups.com
Yes, please. If it's still a problem then we dig deeper. 

--
Oliver Beckstein

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

Samuel Jayakanthan

unread,
Mar 23, 2016, 10:49:39 AM3/23/16
to MDnalysis discussion
Hi Tyler and Oliver,
I installed the newest version. and re-ran my script.

In [3]: print(MDAnalysis.__version__)
0.14.0

I still do see the same "aliasing/overlapping" effect shown above.

My script writes out the data as prescribed by the hbond tutorial on the MDAnalysis website. please let me know if you require more information..


                f1.write(str(table1[i][0]))
                f1.write(" ")
                f1.write(str(table1[i][1]))
                f1.write(" ")
                f1.write(str(table1[i][2]))
                f1.write(" ")
                f1.write(str(table1[i][3]))
                ............................... and so on.

thanks for taking the time to trouble shoot. really appreciate it.
-Sam

Oliver Beckstein

unread,
Mar 23, 2016, 12:44:25 PM3/23/16
to mdnalysis-...@googlegroups.com
Hi Sam,

At first glance this looks like a bug. Can you please file an issue in the issue tracker http://issues.mdanalysis.org , and also add a link to this discussion. We'll continue there. 

Oliver


--
Oliver Beckstein

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

Samuel Jayakanthan

unread,
Mar 23, 2016, 4:44:39 PM3/23/16
to MDnalysis discussion
Hi Oliver,
 I just posted the details (including the script) on issue tracked. I will log on from time to time...and look at the issue tracker.
thanks so much
Sam
Reply all
Reply to author
Forward
0 new messages