4.3.1. Hydrogen Bond Analysis — MDAnalysis.analysis.hydrogenbonds.hbond_analysis

Author:

Paul Smith

Year:

2019

Copyright:

GNU Public License v3

New in version 1.0.0.

This module provides methods to find and analyse hydrogen bonds in a Universe.

The HydrogenBondAnalysis class is a new version of the original MDAnalysis.analysis.hbonds.HydrogenBondAnalysis class from the module MDAnalysis.analysis.hbonds.hbond_analysis, which itself was modeled after the VMD HBONDS plugin.

4.3.1.1. Input

Required:

  • universe : an MDAnalysis Universe object

Options:

  • donors_sel [None] : Atom selection for donors. If None, then will be identified via the topology.

  • hydrogens_sel [None] : Atom selection for hydrogens. If None, then will be identified via charge and mass.

  • acceptors_sel [None] : Atom selection for acceptors. If None, then will be identified via charge.

  • d_h_cutoff (Å) [1.2] : Distance cutoff used for finding donor-hydrogen pairs

  • d_a_cutoff (Å) [3.0] : Distance cutoff for hydrogen bonds. This cutoff refers to the D-A distance.

  • d_h_a_angle_cutoff (degrees) [150] : D-H-A angle cutoff for hydrogen bonds.

  • update_selections [True] : If true, will update atom selections at each frame.

4.3.1.2. Output

  • frame : frame at which a hydrogen bond was found

  • donor id : atom id of the hydrogen bond donor atom

  • hydrogen id : atom id of the hydrogen bond hydrogen atom

  • acceptor id : atom id of the hydrogen bond acceptor atom

  • distance (Å): length of the hydrogen bond

  • angle (degrees): angle of the hydrogen bond

Hydrogen bond data are returned in a numpy.ndarray on a “one line, one observation” basis and can be accessed via HydrogenBondAnalysis.results.hbonds:

results = [
    [
        <frame>,
        <donor index (0-based)>,
        <hydrogen index (0-based)>,
        <acceptor index (0-based)>,
        <distance>,
        <angle>
    ],
    ...
]

4.3.1.3. Example use of HydrogenBondAnalysis

The simplest use case is to allow HydrogenBondAnalysis to guess the acceptor and hydrogen atoms, and to identify donor-hydrogen pairs via the bonding information in the topology:

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

u = MDAnalysis.Universe(psf, trajectory)

hbonds = HBA(universe=u)
hbonds.run()

It is also possible to specify which hydrogens and acceptors to use in the analysis. For example, to find all hydrogen bonds in water:

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

u = MDAnalysis.Universe(psf, trajectory)

hbonds = HBA(universe=u, hydrogens_sel='resname TIP3 and name H1 H2', acceptors_sel='resname TIP3 and name OH2')
hbonds.run()

Alternatively, hydrogens_sel and acceptors_sel may be generated via the guess_hydrogens and guess_acceptors. This selection strings may then be modified prior to calling run, or a subset of the universe may be used to guess the atoms. For example, find hydrogens and acceptors belonging to a protein:

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

u = MDAnalysis.Universe(psf, trajectory)

hbonds = HBA(universe=u)
hbonds.hydrogens_sel = hbonds.guess_hydrogens("protein")
hbonds.acceptors_sel = hbonds.guess_acceptors("protein")
hbonds.run()

Slightly more complex selection strings are also possible. For example, to find hydrogen bonds involving a protein and any water molecules within 10 Å of the protein (which may be useful for subsequently finding the lifetime of protein-water hydrogen bonds or finding water-bridging hydrogen bond paths):

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

u = MDAnalysis.Universe(psf, trajectory)

hbonds = HBA(universe=u)

protein_hydrogens_sel = hbonds.guess_hydrogens("protein")
protein_acceptors_sel = hbonds.guess_acceptors("protein")

water_hydrogens_sel = "resname TIP3 and name H1 H2"
water_acceptors_sel = "resname TIP3 and name OH2"

hbonds.hydrogens_sel = f"({protein_hydrogens_sel}) or ({water_hydrogens_sel} and around 10 not resname TIP3})"
hbonds.acceptors_sel = f"({protein_acceptors_sel}) or ({water_acceptors_sel} and around 10 not resname TIP3})"
hbonds.run()

To calculate the hydrogen bonds between different groups, for example a protein and water, one can use the between keyword. The following will find protein-water hydrogen bonds but not protein-protein or water-water hydrogen bonds:

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import (
  HydrogenBondAnalysis as HBA)

u = MDAnalysis.Universe(psf, trajectory)

hbonds = HBA(
  universe=u,
  between=['resname TIP3', 'protein']
  )

protein_hydrogens_sel = hbonds.guess_hydrogens("protein")
protein_acceptors_sel = hbonds.guess_acceptors("protein")

water_hydrogens_sel = "resname TIP3 and name H1 H2"
water_acceptors_sel = "resname TIP3 and name OH2"

hbonds.hydrogens_sel = f"({protein_hydrogens_sel}) or ({water_hydrogens_sel})"
hbonds.acceptors_sel = f"({protein_acceptors_sel}) or ({water_acceptors_sel})"

hbonds.run()

It is further possible to compute hydrogen bonds between several groups with with use of between. If in the above example, between=[[‘resname TIP3’, ‘protein’], [‘protein’, ‘protein’]], all protein-water and protein-protein hydrogen bonds will be found, but no water-water hydrogen bonds.

One can also define hydrogen bonds with atom types:

from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
hbonds = HBA(
             universe=u,
             donors_sel='type 2',
             hydrogens_sel='type 1',
             acceptors_sel='type 2',
            )

In order to compute the hydrogen bond lifetime, after finding hydrogen bonds one can use the lifetime function:

...
hbonds.run()
tau_timeseries, timeseries = hbonds.lifetime()

It is highly recommended that a topology with bond information is used to generate the universe, e.g PSF, TPR, or PRMTOP files. This is the only method by which it can be guaranteed that donor-hydrogen pairs are correctly identified. However, if, for example, a PDB file is used instead, a donors_sel may be provided along with a hydrogens_sel and the donor-hydrogen pairs will be identified via a distance cutoff, d_h_cutoff:

import MDAnalysis
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import (
  HydrogenBondAnalysis as HBA)

u = MDAnalysis.Universe(pdb, trajectory)

hbonds = HBA(
  universe=u,
  donors_sel='resname TIP3 and name OH2',
  hydrogens_sel='resname TIP3 and name H1 H2',
  acceptors_sel='resname TIP3 and name OH2',
  d_h_cutoff=1.2
)
hbonds.run()

4.3.1.4. The class and its methods

class MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis(universe, donors_sel=None, hydrogens_sel=None, acceptors_sel=None, between=None, d_h_cutoff=1.2, d_a_cutoff=3.0, d_h_a_angle_cutoff=150, update_selections=True)[source]

Perform an analysis of hydrogen bonds in a Universe.

Set up atom selections and geometric criteria for finding hydrogen bonds in a Universe.

Hydrogen bond selections with donors_sel , hydrogens_sel, and acceptors_sel may be achieved with either a resname, atom name combination, or when those are absent, with atom type selections.

Parameters:

  • universe (Universe) – MDAnalysis Universe object

  • donors_sel (str) – Selection string for the hydrogen bond donor atoms. If the universe topology contains bonding information, leave donors_sel as None so that donor-hydrogen pairs can be correctly identified.

  • hydrogens_sel (str) – Selection string for the hydrogen bond hydrogen atoms. Leave as None to guess which hydrogens to use in the analysis using guess_hydrogens. If hydrogens_sel is left as None, also leave donors_sel as None so that donor-hydrogen pairs can be correctly identified.

  • acceptors_sel (str) – Selection string for the hydrogen bond acceptor atoms. Leave as None to guess which atoms to use in the analysis using guess_acceptors

  • between (List (optional),) – Specify two selection strings for non-updating atom groups between which hydrogen bonds will be calculated. For example, if the donor and acceptor selections include both protein and water, it is possible to find only protein-water hydrogen bonds - and not protein-protein or water-water - by specifying between=[“protein”, “SOL”]`. If a two-dimensional list is passed, hydrogen bonds between each pair will be found. For example, between=[[“protein”, “SOL”], [“protein”, “protein”]]` will calculate all protein-water and protein-protein hydrogen bonds but not water-water hydrogen bonds. If None, hydrogen bonds between all donors and acceptors will be calculated.

  • d_h_cutoff (float (optional)) – Distance cutoff used for finding donor-hydrogen pairs. Only used to find donor-hydrogen pairs if the universe topology does not contain bonding information

  • d_a_cutoff (float (optional)) – Distance cutoff for hydrogen bonds. This cutoff refers to the D-A distance.

  • d_h_a_angle_cutoff (float (optional)) – D-H-A angle cutoff for hydrogen bonds, in degrees.

  • update_selections (bool (optional)) – Whether or not to update the acceptor, donor and hydrogen lists at each frame.

Note

It is highly recommended that a universe topology with bond information is used, as this is the only way that guarantees the correct identification of donor-hydrogen pairs.

New in version 2.0.0: Added between keyword

Changed in version 2.4.0: Added use of atom types in selection strings for hydrogen atoms, bond donors, or bond acceptors

results.hbonds

A numpy.ndarray which contains a list of all observed hydrogen bond interactions. See Output for more information.

New in version 2.0.0.

hbonds

Alias to the results.hbonds attribute.

Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use results.hbonds instead.

count_by_ids()[source]

Counts the total number hydrogen bonds formed by unique combinations of donor, hydrogen and acceptor atoms.

Returns:

counts – Each row of the array contains the donor atom id, hydrogen atom id, acceptor atom id and the total number of times the hydrogen bond was observed. The array is sorted by frequency of occurrence.

Return type:

numpy.ndarray

Note

Unique hydrogen bonds are determined through a consideration of the hydrogen atom id and acceptor atom id in a hydrogen bond.

count_by_time()[source]

Counts the number of hydrogen bonds per timestep.

Returns:

counts – Contains the total number of hydrogen bonds found at each timestep. Can be used along with HydrogenBondAnalysis.times to plot the number of hydrogen bonds over time.

Return type:

numpy.ndarray

count_by_type()[source]

Counts the total number of each unique type of hydrogen bond.

Returns:

counts – Each row of the array contains the donor resname, donor atom type, acceptor resname, acceptor atom type and the total number of times the hydrogen bond was found.

Return type:

numpy.ndarray

Note

Unique hydrogen bonds are determined through a consideration of the resname and atom type of the donor and acceptor atoms in a hydrogen bond.

guess_acceptors(select='all', max_charge=-0.5)[source]

Guesses which atoms could be considered acceptors in the analysis.

Acceptor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types.

Parameters:

  • select (str (optional)) – Selection string for atom group from which acceptors will be identified. (e.g., (resname X and name O1) or type 2)

  • max_charge (float (optional)) – The charge of an acceptor atom must be less than this value.

Returns:

potential_acceptors – String containing the resname and name of all atoms that potentially capable of forming hydrogen bonds.

Return type:

str

Notes

Acceptor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types.

This function makes use of and atomic charges to identify which atoms could be considered acceptor atoms in the hydrogen bond analysis. If an atom has an atomic charge less than max_charge then it is considered capable of participating in hydrogen bonds.

If acceptors_sel is None, this function is called to guess the selection.

Alternatively, this function may be used to quickly generate a str of potential acceptor atoms involved in hydrogen bonding. This str may then be modified before being used to set the attribute acceptors_sel.

Changed in version 2.4.0: Added ability to use atom types

guess_donors(select='all', max_charge=-0.5)[source]

Guesses which atoms could be considered donors in the analysis. Only use if the universe topology does not contain bonding information, otherwise donor-hydrogen pairs may be incorrectly assigned.

Parameters:

  • select (str (optional)) – Selection string for atom group from which donors will be identified. (e.g., (resname X and name O1) or type 2)

  • max_charge (float (optional)) – The charge of a donor atom must be less than this value.

Returns:

potential_donors – String containing the resname and name of all atoms that are potentially capable of forming hydrogen bonds.

Return type:

str

Notes

Donor selections may be achieved with either a resname, atom name combination, or when those are absent, atom types.

This function makes use of and atomic charges to identify which atoms could be considered donor atoms in the hydrogen bond analysis. If an atom has an atomic charge less than max_charge, and it is within d_h_cutoff of a hydrogen atom, then it is considered capable of participating in hydrogen bonds.

If donors_sel is None, and the universe topology does not have bonding information, this function is called to guess the selection.

Alternatively, this function may be used to quickly generate a str of potential donor atoms involved in hydrogen bonding. This str may then be modified before being used to set the attribute donors_sel.

Changed in version 2.4.0: Added ability to use atom types

guess_hydrogens(select='all', max_mass=1.1, min_charge=0.3, min_mass=0.9)[source]

Guesses which hydrogen atoms should be used in the analysis.

Parameters:

  • select (str (optional)) – Selection string for atom group from which hydrogens will be identified. (e.g., (resname X and name H1) or type 2)

  • max_mass (float (optional)) – The mass of a hydrogen atom must be less than this value.

  • min_mass (float (optional)) – The mass of a hydrogen atom must be greater than this value.

  • min_charge (float (optional)) – The charge of a hydrogen atom must be greater than this value.

Returns:

potential_hydrogens – String containing the resname and name of all hydrogen atoms potentially capable of forming hydrogen bonds.

Return type:

str

Notes

Hydrogen selections may be achieved with either a resname, atom name combination, or when those are absent, atom types.

This function makes use of atomic masses and atomic charges to identify which atoms are hydrogen atoms that are capable of participating in hydrogen bonding. If an atom has a mass less than max_mass and an atomic charge greater than min_charge then it is considered capable of participating in hydrogen bonds.

If hydrogens_sel is None, this function is called to guess the selection.

Alternatively, this function may be used to quickly generate a str of potential hydrogen atoms involved in hydrogen bonding. This str may then be modified before being used to set the attribute hydrogens_sel.

Changed in version 2.4.0: Added ability to use atom types

lifetime(tau_max=20, window_step=1, intermittency=0)[source]

Computes and returns the time-autocorrelation (HydrogenBondLifetimes) of hydrogen bonds.

Before calling this method, the hydrogen bonds must first be computed with the run() function. The same start, stop and step parameters used in finding hydrogen bonds will be used here for calculating hydrogen bond lifetimes. That is, the same frames will be used in the analysis.

Unique hydrogen bonds are identified using hydrogen-acceptor pairs. This means an acceptor switching to a different hydrogen atom - with the same donor - from one frame to the next is considered a different hydrogen bond.

Please see MDAnalysis.lib.correlations.autocorrelation() and MDAnalysis.lib.correlations.intermittency() functions for more details.

Parameters:

  • window_step (int, optional) – The number of frames between each t(0).

  • tau_max (int, optional) – Hydrogen bond lifetime is calculated for frames in the range 1 <= tau <= tau_max

  • intermittency (int, optional) – The maximum number of consecutive frames for which a bond can disappear but be counted as present if it returns at the next frame. An intermittency of 0 is equivalent to a continuous autocorrelation, which does not allow for hydrogen bond disappearance. For example, for intermittency=2, any given hydrogen bond may disappear for up to two consecutive frames yet be treated as being present at all frames. The default is continuous (intermittency=0).

Returns:

  • tau_timeseries (np.array) – tau from 1 to tau_max

  • timeseries (np.array) – autcorrelation value for each value of tau