Hydrogen bond distances

[horlust]

Hello,

I slightly modified the hbond_analysis file in order to obtain the atom ids that form an hydrogen bond for each frame:

def list_of_sets(self):
        # Extract the hydrogen bonds IDs only in the format
        # [set(superset(x1,x2), superset(x3,x4)), ..]
        found_hydrogen_bonds = [set() for _ in self.frames]
        for frame_index, frame in enumerate(self.frames):
            for hbond in self.results.hbonds[self.results.hbonds[:, 0] == frame]:
                found_hydrogen_bonds[frame_index].add(frozenset(hbond[2:4]))
                
        return found_hydrogen_bonds

Then, for a system containing water, I performed this analysis and got an array of arrays containing the atom ids:

# define universe with MDAnalysis
u = mda.Universe('md_2.tpr', 'md_2.trr')

hbonds = HBA(universe=u, d_a_cutoff=3.5, d_h_a_angle_cutoff=30)
hbonds.run()

from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import *
list_of_sets = hbonds.list_of_sets()

I wanted to see the evolution of the distance between pairs of atoms throughout the simulation. I took 2 atom ids that were bonded in the first frame and analyzed the evolution of the distance:

from MDAnalysis.analysis.distances import distance_array

id_1 = 1172
id_2 = 1158

time = np.arange(100)
distances = []

for t in time:
    coordinates_1 = u.trajectory[t].positions[id_1]
    coordinates_2 = u.trajectory[t].positions[id_2]
    distance = MDAnalysis.analysis.distances.distance_array(coordinates_1, coordinates_2)
    distances.append(float(distance))

plt.plot(time, distances)
plt.show()

However, from what I can see the distance between these two atoms (which were bonded in the first frame, according to list_of_sets) in the first frame is not < 3.5 A:

I tried other atom ids that were bonded in the first frame and I often get a distance bigger than 3.5 A. What am I doing wrong?

Thanks in advance! :-)

[lunamorrow]

Hi @horlust, I am new to MDA myself but I have been doing some work on the hydrogen_analysis class and I think I can help you. Since you haven't provided a donor selection, the _get_dh_pairs method will guess the donor atoms for the hydrogen bonds based off of all the bonds in Universe, so from what I can see the d_a_cutoff you passed in is not considered until the first frame is run, which would be frame 2. So while it is identifying bonds correctly from your TPR file, they may not be within the cutoff you want initially. Otherwise, you should be able to apply d_a_cutoff straight away if you supply a selection string for both the donors and hydrogens args. The second option is probably best, and is straightforward for a system with just water.

u = mda.Universe('md_2.tpr', 'md_2.trr')

# the selection strings may vary depending on how atoms are defined in your files
hbonds = HBA(universe=u, donors_sel="name OH2", hydrogens_sel="name H1 H2", d_a_cutoff=3.5, d_h_a_angle_cutoff=30)
hbonds.run()

from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import *
list_of_sets = hbonds.list_of_sets()

Can you check the distances between bonded pairs identified at frame 2, and/or pass in the selection strings when you make the HBA object and check the distances at frame 1? If you can post them as a list of numbers that would be great for me to see. I hope this helps! :)