Hydrogen bond calculation

[horlust]

Hi everyone,

I have a simulation box containing 10 molecules of a solute called B00 in aqueous solution. I've tried two methods of computing the hydrogen bonds of the solute SOL with water, and I thought they would give the same results, but it's not the case, and I wanted to understand why.

First method: simple calculation.

u = mda.Universe(f'prod_10_molecules_{T}K.tpr', f'prod_10_molecules_{T}K.xtc')
hbonds = HBA(universe=u, between=['resname SOL', 'resname B00'], d_a_cutoff=distance, d_h_a_angle_cutoff=angle)
hbonds.run(verbose=True)

Second method: defining selections.

u = mda.Universe(f'prod_10_molecules_{T}K.tpr', f'prod_10_molecules_{T}K.xtc')
hbonds = HBA(universe=u, between=['resname SOL', 'resname B00'], d_a_cutoff=distance, d_h_a_angle_cutoff=angle)

hydrogen_donors_sel = f"resname B00"
hydrogens_sel = hbonds.guess_hydrogens("resname B00")
 acceptors_sel = hbonds.guess_acceptors("resname B00")
    
water_hydrogens_sel = "resname SOL and name HW1 HW2"
water_acceptors_sel = "resname SOL and name OW"
water_hydrogen_donors_sel = "resname SOL and name OW"

hbonds.hydrogens_sel = f"({hydrogens_sel}) or ({water_hydrogens_sel})"
hbonds.acceptors_sel = f"({acceptors_sel}) or ({water_acceptors_sel})"
hbonds.donors_sel = f"({hydrogen_donors_sel}) or ({water_hydrogen_donors_sel})"

hbonds.run(verbose=True)

I've noticed that by defining selections the calculations run faster, and there are more hydrogen bonds per frame. Why are these calculation results so different?

Thanks! :)

[orbeckst]

I would compare the auto-guessed donors and acceptors to your manual selections. The auto-guessing is mostly geared towards water and biomolecules such as proteins. If your B00 is a random organic molecule then you likely need to define all of this yourself.

Maybe @p-j-smith has some additional insights?

[p-j-smith]

As, @orbeckst said, it's worth checking which atoms were selected in each case, alhtough as far as I can see the selections should be the same. A couple other points:

• as you're using a tpr file for your topology, you can omit setting hbonds.donors_sel. This way, donor-hydrogen pairs will be found based on the topology. If you set a donor selection, the donor-hydrogen pairs are determined based on distances instead, which is slower and less reliable.
• Also, as your atom selections don't update at each frame, you can pass update_selections=False when initialising HBA. This will avoid finding donor-hydrogen pairs at each frame and so the analysis should run faster.