Problem with Baker-Hubbard H-bond analysis

[karel10]

Hi, I am a Msc. student of bioinformatics, and I am trying to reproduce the example found at the MDtraj page:
http://mdtraj.org/1.7.2/examples/hbonds.html
I succeeded with the proposed example structure (which seems to be an ensemble of NMR structures), but am not able to successfully complete it with my own trajectories, generated by OpenMM. I have come to the conclusion that there may be some kind of bug which does not allow for the correct determination of H-bonds (with the exception -for some reason-) of .pdb files of NMR-obtained structures.

To try to get it to work in other trajectories/structures I have tried:
a) Using a .dcd trajectory plus a .pdb topology file.
b) Using a .pdb trajectory file (or a single structure file).
c) Using a .hdf5 trajectory file.
d) Looking for other structures at the PDB database and see if they give some result. Apparently the only ones that work are NMR structures. As for example: 2ysh, 1trz, 2n2g, 2n39.

Baker-Hubbard-Test.tar.gz.zip

Two kinds of errors occur. Find next the first kind of error, found when using a .pdb trajectory:

---

IndexError Traceback (most recent call last)
in ()
1 color = itertools.cycle(['r', 'b', 'gold'])
2 for i in [2, 3, 4]:
----> 3 plt.hist(da_distances[:, i], color=next(color), label=label(hbonds[i]), alpha=0.5)
4 plt.legend()
5 plt.ylabel('Freq');
IndexError: index 2 is out of bounds for axis 1 with size 0

Find next the second kind of error, if using a .dcd plus a .pdb topology; or if using a hdf5 file:

---

ValueError Traceback (most recent call last)
in ()
----> 1 hbonds = md.baker_hubbard(t, periodic=False)
2 label = lambda hbond : '%s -- %s' % (t.topology.atom(hbond[0]), t.topology.atom(hbond[2]))
3 for hbond in hbonds:
4 print(label(hbond))
/home/kdepourcq/anaconda3/lib/python3.5/site-packages/mdtraj/geometry/hbond.py in baker_hubbard(traj, freq, exclude_water, periodic)
275 nh_donors = get_donors('N', 'H')
276 oh_donors = get_donors('O', 'H')
--> 277 xh_donors = np.concatenate((nh_donors, oh_donors))
278
279 if len(xh_donors) == 0:
ValueError: all the input arrays must have same number of dimensions

I am using:
Ubuntu 14.04 LTS
mdtraj version 1.7.2
anaconda: 4.0.0-np110py35_0
jupyter: 1.0.0-py35_2
notebook: 4.1.0-py35_1
OpenMM 7.0

While struggling with this, I have taken a .pdb file for which the algorythm works and another for which it doesn't. I have trimmed them down to see if I can find a clear difference indicating the problem. Alas, I have been uncapable of it. I enclose both files, and a reduced, minimal version as well. In those, I have even tried to add the same number of spaces to the "bad" file as are found in the "good" one. Alas, with no success, and finding no obvious difference explaining this behaviour.
I enclose a "fake" .zip file with those .pdb files. It is actually a .taz.gz to which I have just typed an additional .zip extension (the application would not accept the file otherwise)
Baker-Hubbard-Test.tar.gz.zip

EDIT:
I have posted the same issue as well here:
http://discourse.mdtraj.org/t/calculating-h-bonds-problematic-trajectory-files/154
There seems to be at least somebody else with the same issue.

[connorbrinton]

I can replicate this issue, and confirm that it's resolved with pull request #1139. The root problem in the bad trajectory is that there are nitrogen-hydrogen bonds present, but no oxygen-hydrogen bonds. When mdtraj asks numpy to concatenate the two lists together, it can't figure out what shape the empty oxygen-hydrogen bond list should be, so it can't concatenate them.

[karel10]

Thanks a lot, Connor.

I will check it out and see how it works.
Kind and thankful greetings from Barcelona, Spain,

Karel De Pourcq

2016-07-11 19:24 GMT+02:00 Connor Brinton notifications@github.com:

I can replicate this issue, and confirm that it's resolved with pull
request #1139 #1139. The root
problem in the bad trajectory is that there are nitrogen-hydrogen bonds
present, but no oxygen-hydrogen bonds. When mdtraj asks numpy to
concatenate the two lists together, it can't figure out what shape the
empty oxygen-hydrogen bond list should be, so it can't concatenate them.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#1127 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/AO9jbydYJbupd9begATVSdipc6JtkoFrks5qUnw0gaJpZM4IzOVT
.

[karel10]

P.S. It seems that for a week or more the MDtraj official forum is down
(502 bad gateway error).

2016-07-12 10:55 GMT+02:00 Karel karel10@gmail.com:

Thanks a lot, Connor.

I will check it out and see how it works.
Kind and thankful greetings from Barcelona, Spain,

Karel De Pourcq

2016-07-11 19:24 GMT+02:00 Connor Brinton notifications@github.com:

I can replicate this issue, and confirm that it's resolved with pull
request #1139 #1139. The root
problem in the bad trajectory is that there are nitrogen-hydrogen bonds
present, but no oxygen-hydrogen bonds. When mdtraj asks numpy to
concatenate the two lists together, it can't figure out what shape the
empty oxygen-hydrogen bond list should be, so it can't concatenate them.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
#1127 (comment),
or mute the thread
https://github.com/notifications/unsubscribe/AO9jbydYJbupd9begATVSdipc6JtkoFrks5qUnw0gaJpZM4IzOVT
.

[cxhernandez]

It seems that for a week or more the MDtraj official forum is down
(502 bad gateway error).

cc @rmcgibbo

[r08b46009]

May I ask about this issue?
I have the same problem that some of the PDB file like 1JJ2.pdb or 1alk.pdb can't be obtained their hydrogen bond after the 'baker_hubbard' function. How did you solve it?