# 4.4.1. Hydrogen Bond Analysis — `MDAnalysis.analysis.hydrogenbonds.hbond_analysis`

Author:

  Paul Smith

Year:

  2019

Copyright:

  Lesser GNU Public License v2.1+

Added 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.4.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.4.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.4.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.4.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)

  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.

  Added 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

  Changed in version 2.8.0: Introduced `get_supported_backends()` allowing for parallel execution on `multiprocessing` and `dask` backends.

  Changed in version 2.10.0: The donors_sel, hydrogens_sel, and acceptors_sel are stored as properties now and update the internal _donor, _hydrogens, and _acceptors when they are modified.

  results.hbonds

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

    Added 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.

  property acceptors_sel

    Selection string for the hydrogen bond acceptor atoms.

    Added in version 2.10.0.

  count_by_ids()

    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()

    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()

    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.

  property donors_sel

    Selection string for the hydrogen bond donor atoms.

    Added in version 2.10.0.

  classmethod get_supported_backends()

    Tuple with backends supported by the core library for a given class. User can pass either one of these values as `backend=...` to `run()` method, or a custom object that has `apply` method (see documentation for `run()`):

    > - ‘serial’: no parallelization

    > - ‘multiprocessing’: parallelization using multiprocessing.Pool

    > - ‘dask’: parallelization using dask.delayed.compute(). Requires installation of mdanalysis[dask]

    If you want to add your own backend to an existing class, pass a `backends.BackendBase` subclass (see its documentation to learn how to implement it properly), and specify `unsupported_backend=True`.

    Returns:

      names of built-in backends that can be used in `run(backend=...)()`

    Return type:

      tuple

    Added in version 2.8.0.

  guess_acceptors(select='all', max_charge=-0.5)

    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)

    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)

    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

  property hydrogens_sel

    Selection string for the hydrogen bond hydrogen atoms.

    Added in version 2.10.0.

  lifetime(tau_max=20, window_step=1, intermittency=0)

    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

  property parallelizable

    Boolean mark showing that a given class can be parallelizable with split-apply-combine procedure. Namely, if we can safely distribute `_single_frame()` to multiple workers and then combine them with a proper `_conclude()` call. If set to `False`, no backends except for `serial` are supported.

    Note

    If you want to check parallelizability of the whole class, without explicitly creating an instance of the class, see `_analysis_algorithm_is_parallelizable`. Note that you setting it to other value will break things if the algorithm behind the analysis is not trivially parallelizable.

    Returns:

      if a given `AnalysisBase` subclass instance is parallelizable with split-apply-combine, or not

    Return type:

      bool

    Added in version 2.8.0.

  run(start: int = None, stop: int = None, step: int = None, frames: Iterable = None, verbose: bool = None, n_workers: int = None, n_parts: int = None, backend: str | BackendBase = None, *, unsupported_backend: bool = False, progressbar_kwargs=None)

    Perform the calculation

    Parameters:

      - start (int, optional) – start frame of analysis

      - stop (int, optional) – stop frame of analysis

      - step (int, optional) – number of frames to skip between each analysed frame

      - frames (array_like, optional) –

        array of integers or booleans to slice trajectory; `frames` can only be used instead of `start`, `stop`, and `step`. Setting both `frames` and at least one of `start`, `stop`, `step` to a non-default value will raise a `ValueError`.

        Added in version 2.2.0.

      - verbose (bool, optional) – Turn on verbosity

      - progressbar_kwargs (dict, optional) – ProgressBar keywords with custom parameters regarding progress bar position, etc; see `MDAnalysis.lib.log.ProgressBar` for full list. Available only for `backend='serial'`

      - backend (Union[str, BackendBase], optional) –

        By default, performs calculations in a serial fashion. Otherwise, user can choose a backend: `str` is matched to a builtin backend (one of `serial`, `multiprocessing` and `dask`), or a `MDAnalysis.analysis.results.BackendBase` subclass.

        Added in version 2.8.0.

      - n_workers (int) –

        positive integer with number of workers (processes, in case of built-in backends) to split the work between

        Added in version 2.8.0.

      - n_parts (int, optional) –

        number of parts to split computations across. Can be more than number of workers.

        Added in version 2.8.0.

      - unsupported_backend (bool, optional) –

        if you want to run your custom backend on a parallelizable class that has not been tested by developers, by default False

        Added in version 2.8.0.

    Changed in version 2.2.0: Added ability to analyze arbitrary frames by passing a list of frame indices in the frames keyword argument.

    Changed in version 2.5.0: Add progressbar_kwargs parameter, allowing to modify description, position etc of tqdm progressbars

    Changed in version 2.8.0: Introduced `backend`, `n_workers`, `n_parts` and `unsupported_backend` keywords, and refactored the method logic to support parallelizable execution.