4.3.3. Water Bridge analysis — MDAnalysis.analysis.hydrogenbonds.WaterBridgeAnalysis
¶
Author: | Zhiyi Wu |
---|---|
Year: | 2017-2018 |
Copyright: | GNU Public License v3 |
Maintainer: | Zhiyi Wu <zhiyi.wu@gtc.ox.ac.uk>, @xiki-tempula on GitHub |
Given a Universe
(simulation
trajectory with 1 or more frames) measure all water bridges for each
frame between selections 1 and 2.
Water bridge is defined as a bridging water which simultaneously forms
two hydrogen bonds with atoms from both selection 1 and selection 2.
A water bridge can form between two hydrogen bond acceptors.
e.g. -CO2-:···H−O−H···:-O2C-
A water bridge can also form between two hydrogen bond donors.
e.g. -NH···:O:···HN- (where O is the oxygen of a bridging water)
A hydrogen bond acceptor and another hydrogen bond donor can be bridged by a water.
e.g. -CO2-:···H−O:···HN- (where H−O is part of H−O−H)
A higher order water bridge is defined as more than one water bridging hydrogen bond acceptor and donor. An example of a second order water bridge:
e.g. -CO2-:···H−O:···H−O:···HN- (where H−O is part of H−O−H)
The following keyword arguments are important to control the behaviour of the water bridge analysis:
- water_selection (
resname SOL
): the selection string for the bridging water- order the maximum number of water bridging both ends
- donor-acceptor distance (Å): 3.0
- Angle cutoff (degrees): 120.0
- forcefield to switch between default values for different force fields
- donors and acceptors atom types (to add additional atom names)
4.3.3.1. Theory¶
This module attempts to find multi-order water bridge by an approach similar to breadth-first search, where the first solvation shell of selection 1 is selected, followed by the selection of the second solvation shell as well as any hydrogen bonding partner from selection 1. After that, the third solvation shell, as well as any binding partners from selection 2, are detected. This process is repeated until the maximum order of water bridges is reached.
4.3.3.2. Output as Network¶
Since the waters connecting the two ends of the selections are by nature a
network. We provide a network representation of the water network. Water bridge
data are returned per frame, which is stored in
WaterBridgeAnalysis.results.network
. Each frame is represented as a
dictionary, where the keys are the hydrogen bonds originating from selection
1 and the values are new dictionaries representing the hydrogen bonds coming
out of the corresponding molecules making hydrogen bonds with selection 1.
As for the hydrogen bonds which reach the selection 2, the values of the corresponding keys are None. One example where selection 1 and selection 2 are joined by one water molecule (A) which also hydrogen bond to another water (B) which also hydrogen bond to selection 2 would be represented as
# (selection 1)-O:···H-O(water 1):···H-(selection 2)
# | :
# H·············O-H(water2)
# H
{(sele1_acceptor, None, water1_donor, water1_donor_heavy, distance, angle):
{(water1_acceptor, None, sele2_donor, sele2_donor_heavy,
distance, angle): None},
{(water1_donor, water1_donor_heavy, water2_acceptor, None,
distance, angle):
{(water2_acceptor, None, sele2_donor, sele2_donor_heavy,
distance, angle): None}
},
}
The atoms are represented by atom index and if the atom is hydrogen bond donor,
it is followed by the index of the corresponding heavy atom
(donor_proton, donor_heavy_atom)
.
If the atom is a hydrogen bond acceptor, it is followed by none.
4.3.3.3. Output as Timeseries¶
For lower order water bridges, it might be desirable to represent the
connections as WaterBridgeAnalysis.results.timeseries
. The results
are returned per frame and are a list of hydrogen bonds between the selection
1 or selection 2 and the bridging waters. Due to the complexity of the higher
order water bridge and the fact that one hydrogen bond between two waters can
appear in both third and fourth order water bridge, the hydrogen bonds in the
WaterBridgeAnalysis.results.timeseries
attribute are generated in a
depth-first search manner to avoid duplication. Example code of how
WaterBridgeAnalysis.results.timeseries
is generated:
def network2timeseries(network, timeseries):
'''Traverse the network in a depth-first fashion.
expand_timeseries will expand the compact representation to the
familiar timeseries representation.'''
if network is None:
return
else:
for node in network:
timeseries.append(expand_timeseries(node))
network2timeseries(network[node], timeseries)
timeseries = []
network2timeseries(network, timeseries)
An example would be.
results = [
[ # frame 1
[ <donor index>, <acceptor index>,
(<donor residue name>, <donor residue number>, <donor atom name>),
(<acceptor residue name>, <acceptor residue number>,
<acceptor atom name>),
<distance>, <angle>],
....
],
[ # frame 2
[ ... ], [ ... ], ...
],
...
]
Using the WaterBridgeAnalysis.generate_table()
method one can reformat
the results as a flat “normalised” table that is easier to import into a
database or dataframe for further processing.
4.3.3.4. Detection of water bridges¶
Water bridges are recorded if a bridging water simultaneously forms hydrogen bonds with selection 1 and selection 2.
Hydrogen bonds are detected based on a geometric criterion:
- The distance between acceptor and hydrogen is less than or equal to distance (default is 3 Å).
- The angle between donor-hydrogen-acceptor is greater than or equal to angle (default is 120º).
The cut-off values angle and distance can be set as keywords to
WaterBridgeAnalysis
.
Donor and acceptor heavy atoms are detected from atom names. The current defaults are appropriate for the CHARMM27 and GLYCAM06 force fields as defined in Table Default atom names for water bridge analysis.
Hydrogen atoms bonded to a donor are searched based on its distance to the donor. The algorithm searches for all hydrogens (name “H*” or name “[123]H” or type “H”) in the same residue as the donor atom within a cut-off distance of 1.2 Å.
group | donor | acceptor | comments |
---|---|---|---|
main chain | N | O, OC1, OC2 | OC1, OC2 from amber99sb-ildn (Gromacs) |
water | OH2, OW | OH2, OW | SPC, TIP3P, TIP4P (CHARMM27,Gromacs) |
ARG | NE, NH1, NH2 | ||
ASN | ND2 | OD1 | |
ASP | OD1, OD2 | ||
CYS | SG | ||
CYH | SG | possible false positives for CYS | |
GLN | NE2 | OE1 | |
GLU | OE1, OE2 | ||
HIS | ND1, NE2 | ND1, NE2 | presence of H determines if donor |
HSD | ND1 | NE2 | |
HSE | NE2 | ND1 | |
HSP | ND1, NE2 | ||
LYS | NZ | ||
MET | SD | see e.g. [Gregoret1991] | |
SER | OG | OG | |
THR | OG1 | OG1 | |
TRP | NE1 | ||
TYR | OH | OH |
element | donor | acceptor |
---|---|---|
N | N,NT,N3 | N,NT |
O | OH,OW | O,O2,OH,OS,OW,OY |
S | SM |
Donor and acceptor names for the CHARMM27 force field will also work for e.g. OPLS/AA or amber (tested in Gromacs). Residue names in the table are for information only and are not taken into account when determining acceptors and donors. This can potentially lead to some ambiguity in the assignment of donors/acceptors for residues such as histidine or cytosine.
For more information about the naming convention in GLYCAM06 have a look at the Carbohydrate Naming Convention in Glycam.
The lists of donor and acceptor names can be extended by providing lists of
atom names in the donors and acceptors keywords to
WaterBridgeAnalysis
. If the lists are entirely inappropriate
(e.g. when analysing simulations done with a force field that uses very
different atom names) then one should either use the value “other” for
forcefield to set no default values or derive a new class and set the
default list oneself:
class WaterBridgeAnalysis_OtherFF(WaterBridgeAnalysis):
DEFAULT_DONORS = {"OtherFF": tuple(set([...]))}
DEFAULT_ACCEPTORS = {"OtherFF": tuple(set([...]))}
Then simply use the new class instead of the parent class and call it with
`forcefield` = "OtherFF"
. Please also consider contributing the list of
heavy atom names to MDAnalysis.
References
[Gregoret1991] | L.M. Gregoret, S.D. Rader, R.J. Fletterick, and F.E. Cohen. Hydrogen bonds involving sulfur atoms in proteins. Proteins, 9(2):99–107, 1991. 10.1002/prot.340090204. |
4.3.3.5. How to perform WaterBridgeAnalysis
¶
All water bridges between arginine and aspartic acid can be analysed with
import MDAnalysis
import MDAnalysis.analysis.hbonds
u = MDAnalysis.Universe('topology', 'trajectory')
w = MDAnalysis.analysis.hbonds.WaterBridgeAnalysis(u, 'resname ARG',
'resname ASP')
w.run()
The maximum number of bridging waters detected can be changed using the order keyword.
w = MDAnalysis.analysis.hbonds.WaterBridgeAnalysis(u, 'resname ARG',
'resname ASP', order=3)
Thus, a maximum of three bridging waters will be detected.
An example of using the WaterBridgeAnalysis
would be
detecting the percentage of time a certain water bridge exits.
Trajectory u
has two frames, where the first frame contains a water
bridge from the oxygen of the first arginine to one of the oxygens in the
carboxylic group of aspartate (ASP3:OD1). In the second frame, the same water
bridge forms but is between the oxygen of the arginine and the other oxygen in
the carboxylic group (ASP3:OD2).
# index residue id residue name atom name
# 0 1 ARG O
# 1 2 SOL OW
# 2 2 SOL HW1
# 3 2 SOL HW2
# 4 3 ASP OD1
# 5 3 ASP OD2
print(w.results.timeseries)
prints out
[ # frame 1
# A water bridge SOL2 links O from ARG1 to the carboxylic group OD1 of ASP3
[[0,2,('ARG',1, 'O'),('SOL',2,'HW1'), 3.0,180],
[3,4,('SOL',2,'HW2'),('ASP',3,'OD1'), 3.0,180],
],
# frame 2
# Another water bridge SOL2 links O from ARG1 to the other oxygen of the
# carboxylic group OD2 of ASP3
[[0,2,('ARG',1, 'O'),('SOL',2,'HW1'), 3.0,180],
[3,5,('SOL',2,'HW2'),('ASP',3,'OD2'), 3.0,180],
],
]
4.3.3.6. Use count_by_type
¶
We can use the count_by_type()
to
generate the frequence of all water bridges in the simulation.
w.count_by_type()
Returns
[(0, 3, 'ARG', 1, 'O', 'ASP', 3, 'OD1', 0.5),
(0, 4, 'ARG', 1, 'O', 'ASP', 3, 'OD2', 0.5),]
You might think that the OD1 and OD2 are the same oxygen and the aspartate has
just flipped and thus, they should be counted as the same type of water bridge.
The type of the water bridge can be customised by supplying an analysis
function to count_by_type()
.
The analysis function has two parameters. The current and the output. The
current is a list of hydrogen bonds from selection 1 to selection 2, formatted
in the same fashion as WaterBridgeAnalysis.network
, and an example will
be
[
# sele1 acceptor idx, , water donor index, donor heavy atom idx, dist, ang.
[ 0, None, 2, 1, 3.0,180],
# water donor idx, donor heavy atom idx, sele2 acceptor idx, distance, angle.
[ 3, 1, 4, None, 3.0,180],]
Where current[0]
is the first hydrogen bond originating from selection 1
and current[-1]
is the final hydrogen bond ending in selection 2. The
output sums up all the information in the current frame and is a dictionary
with a user-defined key and the value is the weight of the corresponding key.
During the analysis phase, the function analysis iterates through all the water
bridges and modify the output in-place. At the end of the analysis, the keys
from all the frames are collected and the corresponding values will be summed
up and returned.
def analysis(current, output, u):
r'''This function defines how the type of water bridge should be
specified.
Parameters
----------
current : list
The current water bridge being analysed is a list of hydrogen bonds
from selection 1 to selection 2.
output : dict
A dictionary which is modified in-place where the key is the type
of the water bridge and the value is the weight of this type of
water bridge.
u : MDAnalysis.universe
The current Universe for looking up atoms.'''
# decompose the first hydrogen bond.
sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle =
current[0]
# decompose the last hydrogen bond.
atom1, heavy_atom1, sele2_index, sele2_heavy_index, dist, angle =
current[-1]
# expand the atom index to the resname, resid, atom names
sele1 = u.atoms[sele1_index]
sele2 = u.atoms[sele2_index]
(s1_resname, s1_resid, s1_name) = (sele1.resname, sele1.resid,
sele1.name)
(s2_resname, s2_resid, s2_name) = (sele2.resname, sele2.resid,
sele2.name)
# if the residue name is ASP and the atom name is OD2 or OD1,
# the atom name is changed to OD
if s2_resname == 'ASP' and (s2_name == 'OD1' or s2_name == 'OD2'):
s2_name = 'OD'
# setting up the key which defines this type of water bridge.
key = (s1_resname, s1_resid, s1_name, s2_resname, s2_resid, s2_name)
# The number of this type of water bridge is incremented by 1.
output[key] += 1
w.count_by_type(analysis_func=analysis)
Returns
[(('ARG', 1, 'O', 'ASP', 3, 'OD'), 1.0),]
Note that the result is arranged in the format of
(key, the proportion of time)
. When no custom analysis function is supplied
, the key is expanded and is formatted as
[('ARG', 1, 'O', 'ASP', 3, 'OD', 1.0),]
Some people might only interested in contacts between residues and pay no
attention to the details regarding the atom name. However, since multiple water
bridges can exist between two residues, which sometimes can give a result such
that the water bridge between two residues exists 300% of the time. Though this
might be a desirable result for some people, others might want the water bridge
between two residues to be only counted once per frame. This can also be
achieved by supplying an analysis function to
count_by_type()
.
def analysis(current, output, u):
'''This function defines how the type of water bridge should be specified
.
Parameters
----------
current : list
The current water bridge being analysed is a list of hydrogen bonds
from selection 1 to selection 2.
output : dict
A dictionary which is modified in-place where the key is the type
of the water bridge and the value is the weight of this type of
water bridge.
u : MDAnalysis.universe
The current Universe for looking up atoms.
'''
# decompose the first hydrogen bond.
sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle =
current[0]
# decompose the last hydrogen bond.
atom1, heavy_atom1, sele2_index, sele2_heavy_index, dist, angle =
current[-1]
# expand the atom index to the resname, resid, atom names
sele1 = u.atoms[sele1_index]
sele2 = u.atoms[sele2_index]
(s1_resname, s1_resid, s1_name) = (sele1.resname, sele1.resid,
sele1.name)
(s2_resname, s2_resid, s2_name) = (sele2.resname, sele2.resid,
sele2.name)
# s1_name and s2_name are not included in the key
key = (s1_resname, s1_resid, s2_resname, s2_resid)
# Each residue is only counted once per frame
output[key] = 1
w.count_by_type(analysis_func=analysis)
Returns
[(('ARG', 1, 'ASP', 3), 1.0),]
On the other hand, other people may insist that the first order and
second-order water bridges shouldn’t be mixed together, which can also be
achieved by supplying an analysis function to
count_by_type()
.
def analysis(current, output, u):
'''This function defines how the type of water bridge should be specified
.
Parameters
----------
current : list
The current water bridge being analysed is a list of hydrogen bonds
from selection 1 to selection 2.
output : dict
A dictionary which is modified in-place where the key is the type
of the water bridge and the value is the weight of this type of
water bridge.
u : MDAnalysis.universe
The current Universe for looking up atoms.
'''
# decompose the first hydrogen bond.
sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle =
current[0]
# decompose the last hydrogen bond.
atom1, heavy_atom1, sele2_index, sele2_heavy_index, dist, angle =
current[-1]
# expand the atom index to the resname, resid, atom names
sele1 = u.atoms[sele1_index]
sele2 = u.atoms[sele2_index]
(s1_resname, s1_resid, s1_name) = (sele1.resname, sele1.resid,
sele1.name)
(s2_resname, s2_resid, s2_name) = (sele2.resname, sele2.resid,
sele2.name)
# order of the current water bridge is computed
order_of_water_bridge = len(current) - 1
# and is included in the key
key = (s1_resname, s1_resid, s2_resname, s2_resid, order_of_water_bridge)
# The number of this type of water bridge is incremented by 1.
output[key] += 1
w.count_by_type(analysis_func=analysis)
The extra number 1 precede the 1.0 indicate that this is a first order water bridge
[(('ARG', 1, 'ASP', 3, 1), 1.0),]
Some people might not be interested in the interactions related to arginine.
The undesirable interactions can be discarded by supplying an analysis function
to count_by_type()
.
def analysis(current, output, u):
'''This function defines how the type of water bridge should be
specified.
Parameters
----------
current : list
The current water bridge being analysed is a list of hydrogen bonds
from selection 1 to selection 2.
output : dict
A dictionary which is modified in-place where the key is the type
of the water bridge and the value is the number of this type of
water bridge.
u : MDAnalysis.universe
The current Universe for looking up atoms.
'''
# decompose the first hydrogen bond.
sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle =
current[0]
# decompose the last hydrogen bond.
atom1, heavy_atom1, sele2_index, sele2_heavy_index, dist, angle =
current[-1]
# expand the atom index to the resname, resid, atom names
sele1 = u.atoms[sele1_index]
sele2 = u.atoms[sele2_index]
(s1_resname, s1_resid, s1_name) = (sele1.resname, sele1.resid,
sele1.name)
(s2_resname, s2_resid, s2_name) = (sele2.resname, sele2.resid,
sele2.name)
if not s1_resname == 'ARG':
key = (s1_resname, s1_resid, s2_resname, s2_resid)
output[key] += 1
w.count_by_type(analysis_func=analysis)
Returns nothing in this case
[,]
Additional keywords can be supplied to the analysis function by passing through
count_by_type()
.
def analysis(current, output, **kwargs):
...
w.count_by_type(analysis_func=analysis, **kwargs)
4.3.3.7. Use count_by_time
¶
count_by_type()
aggregates data across frames, which
might be desirable in some cases but not the others.
count_by_time()
provides additional functionality
for aggregating results for each frame.
The default behaviour of count_by_time()
is counting
the number of water bridges from selection 1 to selection 2 for each frame.
Take the previous ASP, ARG salt bridge for example:
w.count_by_time()
As one water bridge is found in both frames, the method returns
[(1.0, 1), (2.0, 1), ]
Similar to count_by_type()
The behaviour of count_by_time()
can also be
modified by supplying an analysis function.
Suppose we want to count
- the number of water molecules involved in bridging selection 1 to selection 2.
- only if water bridge terminates in atom name OD1 of ASP.
- only when water bridge is joined by less than two water.
The analysis function can be written as:
def analysis(current, output, u, **kwargs):
'''This function defines how the counting of water bridge should be
specified.
Parameters
----------
current : list
The current water bridge being analysed is a list of hydrogen bonds
from selection 1 to selection 2.
output : dict
A dictionary which is modified in-place where the key is the type
of the water bridge and the value is the number of this type of
water bridge.
u : MDAnalysis.universe
The current Universe for looking up atoms.
'''
# decompose the first hydrogen bond.
sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle =
current[0]
# decompose the last hydrogen bond.
atom1, heavy_atom1, sele2_index, sele2_heavy_index, dist, angle =
current[-1]
# expand the atom index to the resname, resid, atom names
sele1 = u.atoms[sele1_index]
sele2 = u.atoms[sele2_index]
(s1_resname, s1_resid, s1_name) =
(sele1.resname, sele1.resid, sele1.name)
(s2_resname, s2_resid, s2_name) =
(sele2.resname, sele2.resid, sele2.name)
# only the residue name is ASP and the atom name is OD1,
if s2_resname == 'ASP' and s2_name == 'OD1':
# only if the order of water bridge is less than 2
if len(current) -1 < 2:
# extract all water molecules involved in the water bridge
# extract the first water from selection 1
s1_index, to_index, (s1_resname, s1_resid, s1_name),
(to_resname, to_resid, to_name), dist, angle = current[0]
key = (to_resname, to_resid)
output[key] = 1
# extract all the waters between selection 1 and selection 2
for hbond in current[1:-1]:
# decompose the hydrogen bond.
from_index, to_index, (from_resname, from_resid, from_name),
(to_resname, to_resid, to_name), dist, angle = hbond
# add first water
key1 = (from_resname, from_resid)
output[key1] = 1
# add second water
key2 = (to_resname, to_resid)
output[key2] = 1
# extract the last water to selection 2
from_index, s2_index, (from_resname, from_resid, from_name),
(s2_resname, s2_resid, s2_name), dist, angle = current[-1]
key = (from_resname, from_resid)
output[key] = 1
w.count_by_time(analysis_func=analysis)
Returns
[(1.0, 1), (2.0, 0),]
4.3.3.8. Classes¶
-
class
MDAnalysis.analysis.hydrogenbonds.wbridge_analysis.
WaterBridgeAnalysis
(universe, selection1='protein', selection2='not resname SOL', water_selection='resname SOL', order=1, selection1_type='both', update_selection=False, update_water_selection=True, filter_first=True, distance_type='hydrogen', distance=3.0, angle=120.0, forcefield='CHARMM27', donors=None, acceptors=None, output_format='sele1_sele2', debug=None, pbc=False, **kwargs)[source]¶ Perform a water bridge analysis
The analysis of the trajectory is performed with the
WaterBridgeAnalysis.run()
method. The result is stored inWaterBridgeAnalysis.results.timeseries
. Seerun()
for the format.New in version 0.17.0.
Set up the calculation of water bridges between two selections in a universe.
The timeseries is accessible as the attribute
WaterBridgeAnalysis.results.timeseries
.If no hydrogen bonds are detected or if the initial check fails, look at the log output (enable with
MDAnalysis.start_logging()
and set verbose=True
). It is likely that the default names for donors and acceptors are not suitable (especially for non-standard ligands). In this case, either change the forcefield or use customized donors and/or acceptors.Parameters: - universe (Universe) – Universe object
- selection1 (str (optional)) – Selection string for first selection [‘protein’]
- selection2 (str (optional)) – Selection string for second selection [‘not resname SOL’] This string selects everything except water where water is assumed to have a residue name as SOL.
- water_selection (str (optional)) –
Selection string for bridging water selection [‘resname SOL’] The default selection assumes that the water molecules have residue name “SOL”. Change it to the appropriate selection for your specific force field.
However, in theory, this selection can be anything which forms a hydrogen bond with selection 1 and selection 2.
- order (int (optional)) –
The maximum number of water bridges linking both selections. if the order is set to 3, then all the residues linked with less than three water molecules will be detected. [1]
Computation of high order water bridges can be very time-consuming. Think carefully before running the calculation, do you really want to compute the 20th order water bridge between domain A and domain B or you just want to know the third order water bridge between two residues.
- selection1_type ({"donor", "acceptor", "both"} (optional)) – Selection 1 can be ‘donor’, ‘acceptor’ or ‘both’. Note that the value for selection1_type automatically determines how selection2 handles donors and acceptors: If selection1 contains ‘both’ then selection2 will also contain ‘both’. If selection1 is set to ‘donor’ then selection2 is ‘acceptor’ (and vice versa). [‘both’].
- update_selection (bool (optional)) – Update selection 1 and 2 at each frame. Setting to
True
if the selection is not static. Selections are filtered first to speed up performance. Thus, setting toTrue
is recommended if contact surface between selection 1 and selection 2 is constantly changing. [False
] - update_water_selection (bool (optional)) –
Update selection of water at each frame. Setting to
False
is only recommended when the total amount of water molecules in the simulation are small and when water molecules remain static across the simulation.However, in normal simulations, only a tiny proportion of water is engaged in the formation of water bridge. It is recommended to update the water selection and set keyword filter_first to
True
so as to filter out water not residing between the two selections. [True
] - filter_first (bool (optional)) – Filter the water selection to only include water within 4 Å +
order * (2 Å + distance) away from both selection 1 and
selection 2.
Selection 1 and selection 2 are both filtered to only include atoms
with the same distance away from the other selection. [
True
] - distance (float (optional)) – Distance cutoff for hydrogen bonds; only interactions with a H-A distance <= distance (and the appropriate D-H-A angle, see angle) are recorded. (Note: distance_type can change this to the D-A distance.) [3.0]
- angle (float (optional)) – Angle cutoff for hydrogen bonds; an ideal H-bond has an angle of 180º. A hydrogen bond is only recorded if the D-H-A angle is >= angle. The default of 120º also finds fairly non-specific hydrogen interactions and possibly better value is 150º. [120.0]
- forcefield ({"CHARMM27", "GLYCAM06", "other"} (optional)) – Name of the forcefield used. Switches between different
DEFAULT_DONORS
andDEFAULT_ACCEPTORS
values. [“CHARMM27”] - donors (sequence (optional)) – Extra H donor atom types (in addition to those in
DEFAULT_DONORS
), must be a sequence. - acceptors (sequence (optional)) – Extra H acceptor atom types (in addition to those in
DEFAULT_ACCEPTORS
), must be a sequence. - distance_type ({"hydrogen", "heavy"} (optional)) – Measure hydrogen bond lengths between donor and acceptor heavy atoms (“heavy”) or between donor hydrogen and acceptor heavy atom (“hydrogen”). If using “heavy” then one should set the distance cutoff to a higher value such as 3.5 Å. [“hydrogen”]
- output_format ({"sele1_sele2", "donor_acceptor"} (optional)) – Setting the output format for timeseries and table. If set to “sele1_sele2”, for each hydrogen bond, the one close to selection 1 will be placed before selection 2. If set to “donor_acceptor”, the donor will be placed before acceptor. “sele1_sele2”]
- debug (bool (optional)) – If set to
True
enables per-frame debug logging. This is disabled by default because it generates a very large amount of output in the log file. (Note that a logger must have been started to see the output, e.g. usingMDAnalysis.start_logging()
.) - verbose (bool (optional)) – Toggle progress output. (Can also be given as keyword argument to
run()
.)
Notes
In order to speed up processing, atoms are filtered by a coarse distance criterion before a detailed hydrogen bonding analysis is performed (filter_first =
True
).If selection 1 and selection 2 are very mobile during the simulation and the contact surface is constantly changing (i.e. residues are moving farther than 4 Å + order * (2 Å + distance)), you might consider setting the update_selection keywords to
True
to ensure correctness.-
timesteps
¶ List of the times of each timestep. This can be used together with
timeseries
to find the specific time point of a water bridge existence.
-
results.
network
¶ Network representation of the water network.
New in version 2.0.0.
-
network
¶ Alias to the
results.network
attribute.Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use
results.network
instead.
-
table
¶ Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please generate the table with
generate_table()
instead.
-
results.
timeseries
¶ List of hydrogen bonds between the selection 1 or selection 2 and the bridging waters, for each frame.
New in version 2.0.0.
-
timeseries
¶ Alias to the
results.timeseries
attribute.Deprecated since version 2.0.0: Will be removed in MDAnalysis 3.0.0. Please use
results.timeseries
instead.
-
DEFAULT_ACCEPTORS
= {'CHARMM27': ('SG', 'SD', 'OD1', 'OG1', 'OW', 'OH', 'OH2', 'ND1', 'OC2', 'NE2', 'OC1', 'O', 'OG', 'OE2', 'OD2', 'OE1'), 'GLYCAM06': ('NT', 'OH', 'OW', 'OY', 'OS', 'O', 'SM', 'N', 'O2'), 'other': ()}¶ default atom names that are treated as hydrogen acceptors (see Default heavy atom names for CHARMM27 force field.); use the keyword acceptors to add a list of additional acceptor names.
-
DEFAULT_DONORS
= {'CHARMM27': ('SG', 'NZ', 'OG1', 'NE1', 'OW', 'OH', 'OH2', 'ND1', 'NE2', 'NH1', 'ND2', 'NE', 'NH2', 'OG', 'N'), 'GLYCAM06': ('NT', 'OH', 'OW', 'N3', 'N'), 'other': ()}¶ default heavy atom names whose hydrogens are treated as donors (see Default heavy atom names for CHARMM27 force field.); use the keyword donors to add a list of additional donor names.
-
count_by_time
(analysis_func=None, **kwargs)[source]¶ Counts the number of water bridges per timestep.
The counting behaviour can be adjusted by supplying analysis_func. See Use count_by_time for details.
Returns: counts – Returns a time series N(t)
whereN
is the total number of observed water bridges at timet
.Return type: list
-
count_by_type
(analysis_func=None, **kwargs)[source]¶ Counts the frequency of water bridge of a specific type.
If one atom A from selection 1 is linked to atom B from selection 2 through one or more bridging waters, an entity will be created and the proportion of time that this linkage exists in the whole simulation will be calculated.
The identification of a specific type of water bridge can be modified by supplying the analysis_func function. See Use count_by_type for detail.
Returns: counts – Returns a list
containing atom indices for A and B, residue names, residue numbers, atom names (for both A and B) and the fraction of the total time during which the water bridge was detected. This method returns None if methodWaterBridgeAnalysis.run()
was not executed first.Return type: list
-
generate_table
(output_format=None)[source]¶ Generate a normalised table of the results.
Parameters: output_format ({'sele1_sele2', 'donor_acceptor'}) – The output format of the table can be changed a fashion similar to WaterBridgeAnalysis.results.timeseries
by changing the labels of the columns of the participating atoms.Returns: - table (numpy.recarray) – A “tidy” table with one hydrogen bond per row, labeled according to output_format and containing information of atom_1, atom_2, distance, and angle.
- .. versionchanged:: 2.0.0 – Return the generated table (as well as storing it as
table
). - .. deprecated:: 2.0.0 – In release 3.0.0,
generate_table()
will _only_ return the table and no longer store it intable
.
-
r_cov
= {'N': 1.31, 'O': 1.31, 'P': 1.58, 'S': 1.55}¶ A
collections.defaultdict
of covalent radii of common donors (used in :meth`_get_bonded_hydrogens_list` to check if a hydrogen is sufficiently close to its donor heavy atom). Values are stored for N, O, P, and S. Any other heavy atoms are assumed to have hydrogens covalently bound at a maximum distance of 1.5 Å.
-
timesteps_by_type
(analysis_func=None, **kwargs)[source]¶ Frames during which each water bridges existed, sorted by each water bridges.
Processes
WaterBridgeAnalysis.results.network
and returns alist
containing atom indices, residue names, residue numbers (from selection 1 and selection 2) and each timestep at which the water bridge was detected.Similar to
count_by_type()
andcount_by_time()
, the behavior can be adjusted by supplying an analysis_func.Returns: data Return type: list