4.4.2. Hydrogen bond autocorrelation — MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel

Author:

Richard J. Gowers

Year:

2014

Copyright:

Lesser GNU Public License v2.1+

Added in version 0.9.0.

Changed in version 2.0.0: Module moved from MDAnalysis.analysis.hbonds.hbond_autocorrel to MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel.

4.4.2.1. Description

Calculates the time autocorrelation function, \(C_x(t)\), for the hydrogen bonds in the selections passed to it. The population of hydrogen bonds at a given startpoint, \(t_0\), is evaluated based on geometric criteria and then the lifetime of these bonds is monitored over time. Multiple passes through the trajectory are used to build an average of the behaviour.

\[C_x(t) = \left \langle \frac{h_{ij}(t_0) h_{ij}(t_0 + t)}{h_{ij}(t_0)^2} \right\rangle\]

The subscript \(x\) refers to the definition of lifetime being used, either continuous or intermittent. The continuous definition measures the time that a particular hydrogen bond remains continuously attached, whilst the intermittent definition allows a bond to break and then subsequently reform and be counted again. The relevent lifetime, \(\tau_x\), can then be found via integration of this function

\[\tau_x = \int_0^\infty C_x(t) dt`\]

For this, the observed behaviour is fitted to a multi exponential function, using 2 exponents for the continuous lifetime and 3 for the intermittent lifetime.

\[C_x(t) = A_1 \exp( - t / \tau_1) + A_2 \exp( - t / \tau_2) [+ A_3 \exp( - t / \tau_3)]\]

Where the final pre expoential factor \(A_n\) is subject to the condition:

\[A_n = 1 - \sum\limits_{i=1}^{n-1} A_i\]

For further details see [1].

4.4.2.2. Input

Three AtomGroup selections representing the hydrogens, donors and acceptors that you wish to analyse. Note that the hydrogens and donors selections must be aligned, that is hydrogens[0] and donors[0] must represent a bonded pair. For systems such as water, this will mean that each oxygen appears twice in the donors AtomGroup. The function find_hydrogen_donors() can be used to construct the donor AtomGroup

import MDAnalysis as mda
from MDAnalysis.analysis import hydrogenbonds
from MDAnalysis.tests.datafiles import waterPSF, waterDCD
u = mda.Universe(waterPSF, waterDCD)
hydrogens = u.select_atoms('name H*')
donors = hydrogenbonds.find_hydrogen_donors(hydrogens)

Note that this requires the Universe to have bond information. If this isn’t present in the topology file, the MDAnalysis.core.groups.AtomGroup.guess_bonds() method can be used as so

import MDAnalysis as mda
from MDAnalysis.analysis import hydrogenbonds
from MDAnalysis.tests.datafiles import GRO
# we could load the Universe with guess_bonds=True
# but this would guess **all** bonds
u = mda.Universe(GRO)
water = u.select_atoms('resname SOL and not type DUMMY')
# guess bonds only within our water atoms
# this adds the bond information directly to the Universe
water.guess_bonds()
hydrogens = water.select_atoms('type H')
# this is now possible as we guessed the bonds
donors = hydrogenbonds.find_hydrogen_donors(hydrogens)

The keyword exclusions allows a tuple of array addresses to be provided, (Hidx, Aidx),these pairs of hydrogen-acceptor are then not permitted to be counted as part of the analysis. This could be used to exclude the consideration of hydrogen bonds within the same functional group, or to perform analysis on strictly intermolecular hydrogen bonding.

Hydrogen bonds are defined on the basis of geometric criteria; a Hydrogen-Acceptor distance of less then dist_crit and a Donor-Hydrogen-Acceptor angle of greater than angle_crit.

The length of trajectory to analyse in ps, sample_time, is used to choose what length to analyse.

Multiple passes, controlled by the keyword nruns, through the trajectory are performed and an average calculated. For each pass, nsamples number of points along the run are calculated.

4.4.2.3. Output

All results of the analysis are available through the HydrogenBondAutoCorrel.solution attribute. This is a dictionary with the following keys

• results: The raw results of the time autocorrelation function.

• time: Time axis, in ps, for the results.

• fit: Results of the exponential curve fitting procedure. For the continuous lifetime these are (A1, tau1, tau2), for the intermittent lifetime these are (A1, A2, tau1, tau2, tau3).

• tau: Calculated time constant from the fit.

• estimate: Estimated values generated by the calculated fit.

The results and time values are only filled after the HydrogenBondAutoCorrel.run() method has completed, fit, tau and estimate are filled after the HydrogenBondAutoCorrel.solve() method has been used.

4.4.2.4. Worked Example for Polyamide

This example finds the continuous hydrogen bond lifetime between N-H..O in a polyamide system. This will use the default geometric definition for hydrogen bonds of length 3.0 Å and angle of 130 degrees. It will observe a window of 2.0 ps (sample_time) and try to gather 1000 sample point within this time window (this relies upon the trajectory being sampled frequently enough). This process is repeated for 20 different start points to build a better average.

import MDAnalysis as mda
from MDAnalysis.analysis import hydrogenbonds
from MDAnalysis.tests.datafiles import TRZ_psf, TRZ
import matplotlib.pyplot as plt
# load system
u = mda.Universe(TRZ_psf, TRZ)
# select atoms of interest into AtomGroups
H = u.select_atoms('name Hn')
N = u.select_atoms('name N')
O = u.select_atoms('name O')
# create analysis object
hb_ac = hydrogenbonds.HydrogenBondAutoCorrel(u,
            acceptors=O, hydrogens=H, donors=N,
            bond_type='continuous',
            sample_time=2.0, nsamples=1000, nruns=20)
# call run to gather results
hb_ac.run()
# attempt to fit results to exponential equation
hb_ac.solve()
# grab results from inside object
tau = hb_ac.solution['tau']
time = hb_ac.solution['time']
results = hb_ac.solution['results']
estimate = hb_ac.solution['estimate']
# plot to check!
plt.plot(time, results, 'ro')
plt.plot(time, estimate)
plt.show()

4.4.2.5. Functions and Classes

MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel.find_hydrogen_donors(hydrogens)

Returns the donor atom for each hydrogen

Parameters:

hydrogens (AtomGroup) – the hydrogens that will form hydrogen bonds

Returns:

donors – the donor atom for each hydrogen, found via bond information

Return type:

AtomGroup

Added in version 0.20.0.

class MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel.HydrogenBondAutoCorrel(universe, hydrogens=None, acceptors=None, donors=None, bond_type=None, exclusions=None, angle_crit=130.0, dist_crit=3.0, sample_time=100, time_cut=None, nruns=1, nsamples=50, pbc=True)

Perform a time autocorrelation of the hydrogen bonds in the system.

Results are stored in solution dict.

universeUniverse

MDAnalysis Universe that all selections belong to

hydrogensAtomGroup

AtomGroup of Hydrogens which can form hydrogen bonds

acceptorsAtomGroup

AtomGroup of all Acceptor atoms

donorsAtomGroup

The atoms which are connected to the hydrogens. This group must be identical in length to the hydrogen group and matched, ie hydrogens[0] is bonded to donors[0]. For water, this will mean a donor appears twice in this group, once for each hydrogen.

bond_typestr

Which definition of hydrogen bond lifetime to consider, either ‘continuous’ or ‘intermittent’.

exclusionsndarray, optional

Indices of Hydrogen-Acceptor pairs to be excluded. With nH and nA Hydrogens and Acceptors, a (nH x nA) array of distances is calculated, exclusions is used as a mask on this array to exclude some pairs.

angle_critfloat, optional

The angle (in degrees) which all bonds must be greater than [130.0]

dist_critfloat, optional

The maximum distance (in Angstroms) for a hydrogen bond [3.0]

sample_timefloat, optional

The amount of time, in ps, that you wish to observe hydrogen bonds for [100]

nrunsint, optional

The number of different start points within the trajectory to use [1]

nsamplesint, optional

Within each run, the number of frames to analyse [50]

pbcbool, optional

Whether to consider periodic boundaries in calculations [True]

solution

Results are saved in this dictionary under the following keys:

• results: The raw results of the time autocorrelation function.

• time: Time axis, in ps, for the results.

• fit: Results of the exponential curve fitting procedure. For the continuous lifetime these are (A1, tau1, tau2), for the intermittent lifetime these are (A1, A2, tau1, tau2, tau3).

• tau: Calculated time constant from the fit.

• estimate: Estimated values generated by the calculated fit.

The results and time values are only filled after the run() method has completed, fit, tau and estimate are filled after the solve() method has been used.

Type:

dict

When you use this function please cite [1].

Changed in version 1.0.0: save_results() method was removed. You can instead use np.savez() on HydrogenBondAutoCorrel.solution['time'] and HydrogenBondAutoCorrel.solution['results'] instead.

run(force=False)

Run all the required passes

Intermediate results are stored in the attribute solution (a dict). In order to perform the exponential fits, run the solve() method.

Parameters:

force (bool, optional) – If previous results exist and force=True, then the analysis will run and will overwrite previous results; in the same situation with force=Fales, this method will return immediately without making any changes to solution.

Raises:

ValueError – If all runs start without hydrogen bonds (no hydrogen bonds at the first frame of any run), so that the autocorrelation is undefined.

Changed in version 2.11.0: Runs without hydrogen bonds at their first frame are skipped with a warning instead of giving NaN. A ValueError is raised if all runs are skipped.

solve(p_guess=None)

Fit results to an multi exponential decay and integrate to find characteristic time

Parameters:

p_guess (tuple of floats, optional) – Initial guess for the leastsq fit, must match the shape of the expected coefficients

Continuous defition results are fitted to a double exponential with scipy.optimize.leastsq(), intermittent definition are fit to a triple exponential.

The results of this fitting procedure are saved into the fit, tau and estimate keys in the solution dict.

• fit contains the coefficients, (A1, tau1, tau2) or (A1, A2, tau1, tau2, tau3)

• tau contains the calculated lifetime in ps for the hydrogen bonding

• estimate contains the estimate provided by the fit of the time

  autocorrelation function

In addition, the output of the scipy.optimize.leastsq() function is saved into the solution dict

• infodict

• mesg

• ier

References

[1] (1,2)

Richard J. Gowers and Paola Carbone. A multiscale approach to model hydrogen bonding: the case of polyamide. The journal of Chemical Physics, 142(22):224907, 2015. doi:10.1063/1.4922445.