1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 3# 4# MDAnalysis --- https://www.mdanalysis.org 5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors 6# (see the file AUTHORS for the full list of names) 7# 8# Released under the GNU Public Licence, v2 or any higher version 9# 10# Please cite your use of MDAnalysis in published work: 11# 12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, 13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. 14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics 15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th 16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. 17# 18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. 19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. 20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 21# 22 23# Water Bridge Analysis 24r"""Water Bridge analysis --- :mod:`MDAnalysis.analysis.hbonds.WaterBridgeAnalysis` 25=============================================================================== 26 27:Author: Zhiyi Wu 28:Year: 2017 29:Copyright: GNU Public License v3 30:Maintainer: Zhiyi Wu <zhiyi.wu@gtc.ox.ac.uk>, `@xiki-tempula`_ on GitHub 31 32 33.. _`@xiki-tempula`: https://github.com/xiki-tempula 34 35 36Given a :class:`~MDAnalysis.core.universe.Universe` (simulation 37trajectory with 1 or more frames) measure all water bridges for each 38frame between selections 1 and 2. 39Water bridge is defined as a bridging water which simultaneously forms 40two hydrogen bonds with atoms from both selection 1 and selection 2. 41 42A water bridge can form between two hydrogen bond acceptors. 43 44e.g. -CO\ :sub:`2`\ :sup:`-`:···H−O−H···:\ :sup:`-`\ O\ :sub:`2`\ C- 45 46A water bridge can also form between two hydrogen bond donors. 47 48e.g. -NH···:O:···HN- (where O is the oxygen of a bridging water) 49 50A hydrogen bond acceptor and another hydrogen bond donor can be bridged by a 51water. 52 53e.g. -CO\ :sub:`2`\ :sup:`-`:···H−O:···HN- (where H−O is part of **H−O**\ −H) 54 55The :class:`WaterBridgeAnalysis` class is modeled after the \ 56:class:`~MDAnalysis.analysis.hbonds.hbond_analysis.HydrogenBondAnalysis`. 57 58The following keyword arguments are important to control the behavior of the 59water bridge analysis: 60 61 - *water_selection* (``resname SOL``): the selection string for the bridging 62 water 63 - donor-acceptor *distance* (Å): 3.0 64 - Angle *cutoff* (degrees): 120.0 65 - *forcefield* to switch between default values for different force fields 66 - *donors* and *acceptors* atom types (to add additional atom names) 67 68.. _wb_Analysis_Output: 69 70Output 71------ 72 73The results are a list of hydrogen bonds between the selection 1 or selection 2 74and the bridging water. 75 76Each list is formated similar to the \ :attr:`HydrogenBondAnalysis.timeseries 77<MDAnalysis.analysis.hbonds.hbond_analysis.HydrogenBondAnalysis.timeseries>` 78and contains 79 80 - the **identities** of donor and acceptor heavy-atoms, 81 - the **distance** between the heavy atom acceptor atom and the hydrogen atom 82 - the **angle** donor-hydrogen-acceptor angle (180º is linear). 83 84Water bridge data are returned per frame, which is stored in \ 85:attr:`WaterBridgeAnalysis.timeseries` (In the following description, ``#`` 86indicates comments that are not part of the output.):: 87 88 results = [ 89 [ # frame 1 90 # hbonds linking the selection 1 and selection 2 to the bridging 91 # water 1 92 [ # hbond 1 from selection 1 to the bridging water 1 93 <donor index (0-based)>, 94 <acceptor index (0-based)>, <donor string>, <acceptor string>, 95 <distance>, <angle> 96 ], 97 [ # hbond 2 from selection 1 to the bridging water 1 98 <donor index (0-based)>, 99 <acceptor index (0-based)>, <donor string>, <acceptor string>, 100 <distance>, <angle> 101 ], 102 [ # hbond 1 from selection 2 to the bridging water 1 103 <donor index (0-based)>, 104 <acceptor index (0-based)>, <donor string>, <acceptor string>, 105 <distance>, <angle> 106 ], 107 [ # hbond 2 from selection 2 to the bridging water 1 108 <donor index (0-based)>, 109 <acceptor index (0-based)>, <donor string>, <acceptor string>, 110 <distance>, <angle> 111 ], 112 113 # hbonds linking the selection 1 and selection 2 to the bridging 114 # water 2 115 [ # hbond 1 from selection 1 to the bridging water 2 116 <donor index (0-based)>, 117 <acceptor index (0-based)>, <donor string>, <acceptor string>, 118 <distance>, <angle> 119 ], 120 [ # hbond 1 from selection 2 to the bridging water 2 121 <donor index (0-based)>, 122 <acceptor index (0-based)>, <donor string>, <acceptor string>, 123 <distance>, <angle> 124 ], 125 .... 126 ], 127 [ # frame 2 128 [ ... ], [ ... ], ... 129 ], 130 ... 131 ] 132 133Using the :meth:`WaterBridgeAnalysis.generate_table` method one can reformat 134the results as a flat "normalised" table that is easier to import into a 135database or dataframe for further processing. 136:meth:`WaterBridgeAnalysis.save_table` saves the table to a pickled file. The 137table itself is a :class:`numpy.recarray`. 138 139Detection of water bridges 140--------------------------- 141Water bridges are recorded if a bridging water simultaneously forms two 142hydrogen bonds with selection 1 and selection 2. 143 144Hydrogen bonds are detected as is described in \ 145:class:`~MDAnalysis.analysis.hbonds.hbond_analysis.HydrogenBondAnalysis`, see \ 146:ref:`Detection-of-hydrogen-bonds`. 147 148The lists of donor and acceptor names can be extended by providing lists of 149atom names in the `donors` and `acceptors` keywords to 150:class:`WaterBridgeAnalysis`. If the lists are entirely inappropriate 151(e.g. when analysing simulations done with a force field that uses very 152different atom names) then one should either use the value "other" for 153`forcefield` to set no default values, or derive a new class and set the 154default list oneself:: 155 156 class WaterBridgeAnalysis_OtherFF(WaterBridgeAnalysis): 157 DEFAULT_DONORS = {"OtherFF": tuple(set([...]))} 158 DEFAULT_ACCEPTORS = {"OtherFF": tuple(set([...]))} 159 160Then simply use the new class instead of the parent class and call it with 161`forcefield` = "OtherFF". Please also consider contributing the list of heavy 162atom names to MDAnalysis. 163 164How to perform WaterBridgeAnalysis 165----------------------------------- 166 167All water bridges between arginine and aspartic acid can be analysed with :: 168 169 import MDAnalysis 170 import MDAnalysis.analysis.hbonds 171 172 u = MDAnalysis.Universe('topology', 'trajectory') 173 w = MDAnalysis.analysis.hbonds.WaterBridgeAnalysis(u, 'resname ARG', 'resname ASP') 174 w.run() 175 176The results are stored as the attribute 177:attr:`WaterBridgeAnalysis.timeseries`; see :ref:`wb_Analysis_Output` for the 178format. 179 180An example of using the :attr:`~WaterBridgeAnalysis.timeseries` would be 181detecting the percentage of time a certain water bridge exits. 182 183Trajectory :code:`u` has two frames, where the first frame contains a water 184bridge from the oxygen of the first arginine to the oxygen of the third 185aspartate. No water bridge is detected in the second frame. :: 186 187 print(w.timeseries) 188 189prints out (the comments are not part of the data structure but are added here 190for clarity): :: 191 192 [ # frame 1 193 # A water bridge SOL2 links O from ARG1 and ASP3 194 [[0,1,'ARG1:O', 'SOL2:HW1',3.0,180], 195 [2,3,'SOL2:HW2','ASP3:O', 3.0,180], 196 ], 197 # frame 2 198 # No water bridge detected 199 [] 200 ] 201 202To calculate the percentage, we can iterate through :code:`w.timeseries`. :: 203 204 water_bridge_presence = [] 205 for frame in w.timeseries: 206 if frame: 207 water_bridge_presence.append(True) 208 else: 209 water_bridge_presence.append(False) 210 p_bridge = float(sum(water_bridge_presence))/len(water_bridge_presence) 211 print("Fraction of time with water bridge present: {}".format(p_bridge)) 212 213In the example above, :code:`p_bridge` would become 0.5, i.e., for 50% of the 214trajectory a water bridge was detected between the selected residues. 215 216Alternatively, :meth:`~WaterBridgeAnalysis.count_by_type` can also be used to 217generate the frequence of all water bridges in the simulation. :: 218 219 w.count_by_type() 220 221Returns :: 222 223 [(0, 3, 'ARG', 1, 'O', 'ASP', 3, 'O', 0.5)] 224 225For further data analysis, it is convenient to process the 226:attr:`~WaterBridgeAnalysis.timeseries` data into a normalized table with the 227:meth:`~WaterBridgeAnalysis.generate_table` method, which creates a new data 228structure :attr:`WaterBridgeAnalysis.table` that contains one row for each 229observation of a hydrogen bond:: 230 231 w.generate_table() 232 233This table can then be easily turned into, e.g., a `pandas.DataFrame`_, and 234further analyzed:: 235 236 import pandas as pd 237 df = pd.DataFrame.from_records(w.table) 238 239 240.. _pandas.DataFrame: http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html 241 242Classes 243------- 244 245.. autoclass:: WaterBridgeAnalysis 246 :members: 247 248 .. attribute:: timesteps 249 250 List of the times of each timestep. This can be used together with 251 :attr:`~WaterBridgeAnalysis.timeseries` to find the specific time point 252 of a water bridge existence, or see :attr:`~WaterBridgeAnalysis.table`. 253 254 .. attribute:: table 255 256 A normalised table of the data in 257 :attr:`WaterBridgeAnalysis.timeseries`, generated by 258 :meth:`WaterBridgeAnalysis.generate_table`. It is a 259 :class:`numpy.recarray` with the following columns: 260 261 0. "time" 262 1. "donor_index" 263 2. "acceptor_index" 264 3. "donor_resnm" 265 4. "donor_resid" 266 5. "donor_atom" 267 6. "acceptor_resnm" 268 7. "acceptor_resid" 269 8. "acceptor_atom" 270 9. "distance" 271 10. "angle" 272 273 It takes up more space than :attr:`~WaterBridgeAnalysis.timeseries` but 274 it is easier to analyze and to import into databases or dataframes. 275 276 277 .. rubric:: Example 278 279 For example, to create a `pandas.DataFrame`_ from ``h.table``:: 280 281 import pandas as pd 282 df = pd.DataFrame.from_records(w.table) 283""" 284from __future__ import absolute_import, division 285import six 286 287from collections import defaultdict 288import logging 289import warnings 290 291import numpy as np 292 293from .hbond_analysis import HydrogenBondAnalysis 294from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch 295from MDAnalysis.lib.log import ProgressMeter 296from MDAnalysis.lib import distances 297from MDAnalysis import SelectionWarning 298 299logger = logging.getLogger('MDAnalysis.analysis.wbridges') 300 301class WaterBridgeAnalysis(HydrogenBondAnalysis): 302 """Perform a water bridge analysis 303 304 The analysis of the trajectory is performed with the 305 :meth:`WaterBridgeAnalysis.run` method. The result is stored in 306 :attr:`WaterBridgeAnalysis.timeseries`. See 307 :meth:`~WaterBridgeAnalysis.run` for the format. 308 309 :class:`WaterBridgeAnalysis` uses the same default atom names as 310 :class:`~MDAnalysis.analysis.hbonds.hbond_analysis.HydrogenBondAnalysis`, 311 see :ref:`Default atom names for hydrogen bonding analysis` 312 313 314 .. versionadded:: 0.17.0 315 """ 316 def __init__(self, universe, selection1='protein', 317 selection2='not resname SOL', water_selection='resname SOL', 318 selection1_type='both', update_selection1=False, 319 update_selection2=False, update_water_selection=True, 320 filter_first=True, distance_type='hydrogen', distance=3.0, 321 angle=120.0, forcefield='CHARMM27', donors=None, 322 acceptors=None, debug=None, verbose=False, **kwargs): 323 """Set up the calculation of water bridges between two selections in a 324 universe. 325 326 The timeseries is accessible as the attribute 327 :attr:`WaterBridgeAnalysis.timeseries`. 328 329 If no hydrogen bonds are detected or if the initial check fails, look 330 at the log output (enable with :func:`MDAnalysis.start_logging` and set 331 `verbose` ``=True``). It is likely that the default names for donors 332 and acceptors are not suitable (especially for non-standard 333 ligands). In this case, either change the `forcefield` or use 334 customized `donors` and/or `acceptors`. 335 336 Parameters 337 ---------- 338 universe : Universe 339 Universe object 340 selection1 : str (optional) 341 Selection string for first selection ['protein'] 342 selection2 : str (optional) 343 Selection string for second selection ['not resname SOL'] 344 This string selects everything except water where water is assumed 345 to have a residue name as SOL. 346 water_selection : str (optional) 347 Selection string for bridging water selection ['resname SOL'] 348 The default selection assumes that the water molecules have residue 349 name "SOL". Change it to the appropriate selection for your 350 specific force field. 351 352 However, in theory this selection can be anything which forms 353 hydrogen bond with selection 1 and selection 2. 354 selection1_type : {"donor", "acceptor", "both"} (optional) 355 Selection 1 can be 'donor', 'acceptor' or 'both'. Note that the 356 value for `selection1_type` automatically determines how 357 `selection2` handles donors and acceptors: If `selection1` contains 358 'both' then `selection2` will also contain 'both'. If `selection1` 359 is set to 'donor' then `selection2` is 'acceptor' (and vice versa). 360 ['both']. 361 update_selection1 : bool (optional) 362 Update selection 1 at each frame. Setting to ``True`` if the 363 selection is not static. Selection 1 is filtered first to speed up 364 performance. Thus, setting to ``True`` is recommended if contact 365 surface between selection 1 and selection 2 is constantly 366 changing. [``False``] 367 update_selection2 : bool (optional) 368 Similiar to *update_selection1* but is acted upon selection 2. 369 [``False``] 370 update_water_selection : bool (optional) 371 Update selection of water at each frame. Setting to ``False`` is 372 **only** recommended when the total amount of water molecules in the 373 simulation are small and when water molecules remain static across 374 the simulation. 375 376 However, in normal simulations, only a tiny proportion of water is 377 engaged in the formation of water bridge. It is recommended to 378 update the water selection and set keyword `filter_first` to 379 ``True`` so as to filter out water not residing between the two 380 selections. [``True``] 381 filter_first : bool (optional) 382 Filter the water selection to only include water within 3 * 383 `distance` away from `both` selection 1 and selection 2. 384 Selection 1 and selection 2 are both filtered to only include atoms 385 3 * `distance` away from the other selection. [``True``] 386 distance : float (optional) 387 Distance cutoff for hydrogen bonds; only interactions with a H-A 388 distance <= `distance` (and the appropriate D-H-A angle, see 389 `angle`) are recorded. (Note: `distance_type` can change this to 390 the D-A distance.) [3.0] 391 angle : float (optional) 392 Angle cutoff for hydrogen bonds; an ideal H-bond has an angle of 393 180º. A hydrogen bond is only recorded if the D-H-A angle is 394 >= `angle`. The default of 120º also finds fairly non-specific 395 hydrogen interactions and a possibly better value is 150º. [120.0] 396 forcefield : {"CHARMM27", "GLYCAM06", "other"} (optional) 397 Name of the forcefield used. Switches between different 398 :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS` and 399 :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS` values. 400 ["CHARMM27"] 401 donors : sequence (optional) 402 Extra H donor atom types (in addition to those in 403 :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS`), must be a sequence. 404 acceptors : sequence (optional) 405 Extra H acceptor atom types (in addition to those in 406 :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS`), must be a 407 sequence. 408 distance_type : {"hydrogen", "heavy"} (optional) 409 Measure hydrogen bond lengths between donor and acceptor heavy 410 attoms ("heavy") or between donor hydrogen and acceptor heavy 411 atom ("hydrogen"). If using "heavy" then one should set the 412 *distance* cutoff to a higher value such as 3.5 Å. ["hydrogen"] 413 debug : bool (optional) 414 If set to ``True`` enables per-frame debug logging. This is disabled 415 by default because it generates a very large amount of output in 416 the log file. (Note that a logger must have been started to see 417 the output, e.g. using :func:`MDAnalysis.start_logging`.) 418 verbose : bool (optional) 419 Toggle progress output. (Can also be given as keyword argument to 420 :meth:`run`.) 421 422 Notes 423 ----- 424 In order to speed up processing, atoms are filtered by a coarse 425 distance criterion before a detailed hydrogen bonding analysis is 426 performed (`filter_first` = ``True``). 427 428 If selection 1 and selection 2 are very mobile during the simulation 429 and the contact surface is constantly changing (i.e. residues are 430 moving farther than 3 x `distance`), you might consider setting the 431 `update_selection1` and `update_selection2` keywords to ``True`` to 432 ensure correctness. 433 """ 434 self.water_selection = water_selection 435 self.update_water_selection = update_water_selection 436 super(WaterBridgeAnalysis, self).__init__( 437 universe=universe, selection1=selection1, selection2=selection2, 438 selection1_type=selection1_type, 439 update_selection1=update_selection1, 440 update_selection2=update_selection2, filter_first=filter_first, 441 distance_type=distance_type, distance=distance, angle=angle, 442 forcefield=forcefield, donors=donors, acceptors=acceptors, 443 debug=debug, verbose=verbose, **kwargs) 444 self._update_water_selection() 445 446 def _log_parameters(self): 447 """Log important parameters to the logfile.""" 448 logger.info("WBridge analysis: selection1 = %r (update: %r)", 449 self.selection1, self.update_selection1) 450 logger.info("WBridge analysis: selection2 = %r (update: %r)", 451 self.selection2, self.update_selection2) 452 logger.info("WBridge analysis: water selection = %r (update: %r)", 453 self.water_selection, self.update_water_selection) 454 logger.info("WBridge analysis: criterion: donor %s atom and acceptor \ 455 atom distance <= %.3f A", self.distance_type, 456 self.distance) 457 logger.info("WBridge analysis: criterion: angle D-H-A >= %.3f degrees", 458 self.angle) 459 logger.info("WBridge analysis: force field %s to guess donor and \ 460 acceptor names", self.forcefield) 461 logger.info("WBridge analysis: bonded hydrogen detection algorithm: \ 462 %r", self.detect_hydrogens) 463 464 def _update_selection_1(self): 465 self._s1 = self.u.select_atoms(self.selection1) 466 self._s2 = self.u.select_atoms(self.selection2) 467 if self.filter_first and self._s1: 468 self.logger_debug('Size of selection 1 before filtering:' 469 ' {} atoms'.format(len(self._s1))) 470 ns_selection_1 = AtomNeighborSearch(self._s1) 471 self._s1 = ns_selection_1.search(self._s2, 3. * self.distance) 472 self.logger_debug("Size of selection 1: {0} atoms".format(len(self._s1))) 473 self._s1_donors = {} 474 self._s1_donors_h = {} 475 self._s1_acceptors = {} 476 if self.selection1_type in ('donor', 'both'): 477 self._s1_donors = self._s1.select_atoms( 478 'name {0}'.format(' '.join(self.donors))) 479 self._s1_donors_h = {} 480 for i, d in enumerate(self._s1_donors): 481 tmp = self._get_bonded_hydrogens(d) 482 if tmp: 483 self._s1_donors_h[i] = tmp 484 self.logger_debug("Selection 1 donors: {0}".format(len(self._s1_donors))) 485 self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_donors_h))) 486 if self.selection1_type in ('acceptor', 'both'): 487 self._s1_acceptors = self._s1.select_atoms( 488 'name {0}'.format(' '.join(self.acceptors))) 489 self.logger_debug("Selection 1 acceptors: {0}".format(len(self._s1_acceptors))) 490 491 def _sanity_check(self, selection, htype): 492 """sanity check the selections 1 and 2 493 494 *selection* is 1 or 2, *htype* is "donors" or "acceptors" 495 """ 496 assert selection in (1, 2) 497 assert htype in ("donors", "acceptors") 498 atoms = getattr(self, "_s{0}_{1}".format(selection, htype)) 499 update = getattr(self, "update_selection{0}".format(selection)) 500 if not atoms: 501 errmsg = "No {1} found in selection {0}. " \ 502 "You might have to specify a custom '{1}' keyword.".format( 503 selection, htype) 504 warnings.warn(errmsg, category=SelectionWarning) 505 logger.warning(errmsg) 506 507 def _update_water_selection(self): 508 self._water = self.u.select_atoms(self.water_selection) 509 self.logger_debug('Size of water selection before filtering:' 510 ' {} atoms'.format(len(self._water))) 511 if self.filter_first: 512 ns_water_selection = AtomNeighborSearch(self._water) 513 self._water = ns_water_selection.search(self._s1, 514 3. * self.distance) 515 self._water = ns_water_selection.search(self._s2, 516 3. * self.distance) 517 518 self.logger_debug("Size of water selection: {0} atoms".format(len(self._water))) 519 self._water_donors = {} 520 self._water_donors_h = {} 521 self._water_acceptors = {} 522 if not self._water: 523 logger.warning("Water selection '{0}' did not select any atoms." 524 .format(str(self.water_selection)[:80])) 525 else: 526 self._water_donors = self._water.select_atoms( 527 'name {0}'.format(' '.join(self.donors))) 528 self._water_donors_h = {} 529 for i, d in enumerate(self._water_donors): 530 tmp = self._get_bonded_hydrogens(d) 531 if tmp: 532 self._water_donors_h[i] = tmp 533 self.logger_debug("Water donors: {0}".format(len(self._water_donors))) 534 self.logger_debug("Water donor hydrogens: {0}".format(len(self._water_donors_h))) 535 self._water_acceptors = self._water.select_atoms( 536 'name {0}'.format(' '.join(self.acceptors))) 537 self.logger_debug("Water acceptors: {0}".format(len(self._water_acceptors))) 538 539 def run(self, start=None, stop=None, step=None, verbose=None, debug=None): 540 """Analyze trajectory and produce timeseries. 541 542 Stores the water bridge data per frame as 543 :attr:`WaterBridgeAnalysis.timeseries` (see there for output 544 format). 545 546 Parameters 547 ---------- 548 start : int (optional) 549 starting frame-index for analysis, ``None`` is the first one, 0. 550 `start` and `stop` are 0-based frame indices and are used to slice 551 the trajectory (if supported) [``None``] 552 stop : int (optional) 553 last trajectory frame for analysis, ``None`` is the last one 554 [``None``] 555 step : int (optional) 556 read every `step` between `start` (included) and `stop` (excluded), 557 ``None`` selects 1. [``None``] 558 verbose : bool (optional) 559 toggle progress meter output 560 :class:`~MDAnalysis.lib.log.ProgressMeter` [``True``] 561 debug : bool (optional) 562 enable detailed logging of debugging information; this can create 563 *very big* log files so it is disabled (``False``) by default; 564 setting `debug` toggles the debug status for 565 :class:`WaterBridgeAnalysis`, namely the value of 566 :attr:`WaterBridgeAnalysis.debug`. 567 568 See Also 569 -------- 570 :meth:`WaterBridgeAnalysis.generate_table` : 571 processing the data into a different format. 572 """ 573 self._setup_frames(self.u.trajectory, start, stop, step) 574 575 logger.info("WBridge analysis: starting") 576 logger.debug("WBridge analysis: donors %r", self.donors) 577 logger.debug("WBridge analysis: acceptors %r", self.acceptors) 578 logger.debug("WBridge analysis: water bridge %r", self.water_selection) 579 580 if debug is not None and debug != self.debug: 581 self.debug = debug 582 logger.debug("Toggling debug to %r", self.debug) 583 if not self.debug: 584 logger.debug("WBridge analysis: For full step-by-step debugging output use debug=True") 585 586 self._timeseries = [] 587 self.timesteps = [] 588 self._water_network = [] 589 590 if verbose is None: 591 verbose = self._verbose 592 pm = ProgressMeter(self.n_frames, 593 format="WBridge frame {current_step:5d}: {step:5d}/{numsteps} [{percentage:5.1f}%]\r", 594 verbose=verbose) 595 596 logger.info("Starting analysis (frame index start=%d stop=%d, step=%d)", 597 self.start, self.stop, self.step) 598 599 for progress, ts in enumerate(self.u.trajectory[self.start:self.stop:self.step]): 600 # all bonds for this timestep 601 602 # dict of tuples (atom.index, atom.index) for quick check if 603 # we already have the bond (to avoid duplicates) 604 605 frame = ts.frame 606 timestep = ts.time 607 self.timesteps.append(timestep) 608 pm.echo(progress, current_step=frame) 609 self.logger_debug("Analyzing frame %(frame)d, timestep %(timestep)f ps", vars()) 610 if self.update_selection1: 611 self._update_selection_1() 612 if self.update_selection2: 613 self._update_selection_2() 614 if self.update_water_selection: 615 self._update_water_selection() 616 617 s1_frame_results_dict = defaultdict(list) 618 if (self.selection1_type in ('donor', 'both') and 619 self._water_acceptors): 620 621 self.logger_debug("Selection 1 Donors <-> Water Acceptors") 622 ns_acceptors = AtomNeighborSearch(self._water_acceptors) 623 for i, donor_h_set in self._s1_donors_h.items(): 624 d = self._s1_donors[i] 625 for h in donor_h_set: 626 res = ns_acceptors.search(h, self.distance) 627 for a in res: 628 donor_atom = h if self.distance_type != 'heavy' else d 629 dist = distances.calc_bonds(donor_atom.position, 630 a.position) 631 if dist <= self.distance: 632 angle = distances.calc_angles(d.position, h.position, 633 a.position) 634 angle = np.rad2deg(angle) 635 if angle >= self.angle: 636 self.logger_debug( 637 "S1-D: {0!s} <-> W-A: {1!s} {2:f} A, {3:f} DEG"\ 638 .format(h.index, a.index, dist, angle)) 639 s1_frame_results_dict[(a.resname, a.resid)].append( 640 (h.index, a.index, 641 (h.resname, h.resid, h.name), 642 (a.resname, a.resid, a.name), 643 dist, angle)) 644 645 if (self.selection1_type in ('acceptor', 'both') and 646 self._s1_acceptors): 647 648 self.logger_debug("Selection 1 Acceptors <-> Water Donors") 649 ns_acceptors = AtomNeighborSearch(self._s1_acceptors) 650 for i, donor_h_set in self._water_donors_h.items(): 651 d = self._water_donors[i] 652 for h in donor_h_set: 653 res = ns_acceptors.search(h, self.distance) 654 for a in res: 655 donor_atom = h if self.distance_type != 'heavy' else d 656 dist = distances.calc_bonds(donor_atom.position, 657 a.position) 658 if dist <= self.distance: 659 angle = distances.calc_angles(d.position, h.position, 660 a.position) 661 angle = np.rad2deg(angle) 662 if angle >= self.angle: 663 self.logger_debug( 664 "S1-A: {0!s} <-> W-D: {1!s} {2:f} A, {3:f} DEG"\ 665 .format(a.index, h.index, dist, angle)) 666 s1_frame_results_dict[(h.resname, h.resid)].append( 667 (h.index, a.index, 668 (h.resname, h.resid, h.name), 669 (a.resname, a.resid, a.name), 670 dist, angle)) 671 672 # Narrow down the water selection 673 selection_resn_id = list(s1_frame_results_dict.keys()) 674 if not selection_resn_id: 675 self._timeseries.append([]) 676 continue 677 selection_resn_id = ['(resname {} and resid {})'.format( 678 resname, resid) for resname, resid in selection_resn_id] 679 water_bridges = self._water.select_atoms(' or '.join(selection_resn_id)) 680 self.logger_debug("Size of water bridge selection: {0} atoms".format(len(water_bridges))) 681 if not water_bridges: 682 logger.warning("No water forming hydrogen bonding with selection 1.") 683 water_bridges_donors = water_bridges.select_atoms( 684 'name {0}'.format(' '.join(self.donors))) 685 water_bridges_donors_h = {} 686 for i, d in enumerate(water_bridges_donors): 687 tmp = self._get_bonded_hydrogens(d) 688 if tmp: 689 water_bridges_donors_h[i] = tmp 690 self.logger_debug("water bridge donors: {0}".format(len(water_bridges_donors))) 691 self.logger_debug("water bridge donor hydrogens: {0}".format(len(water_bridges_donors_h))) 692 water_bridges_acceptors = water_bridges.select_atoms( 693 'name {0}'.format(' '.join(self.acceptors))) 694 self.logger_debug("water bridge: {0}".format(len(water_bridges_acceptors))) 695 696 # Finding the hydrogen bonds between water bridge and selection 2 697 s2_frame_results_dict = defaultdict(list) 698 if self._s2_acceptors: 699 self.logger_debug("Water bridge Donors <-> Selection 2 Acceptors") 700 ns_acceptors = AtomNeighborSearch(self._s2_acceptors) 701 for i, donor_h_set in water_bridges_donors_h.items(): 702 d = water_bridges_donors[i] 703 for h in donor_h_set: 704 res = ns_acceptors.search(h, self.distance) 705 for a in res: 706 donor_atom = h if self.distance_type != 'heavy' else d 707 dist = distances.calc_bonds(donor_atom.position, 708 a.position) 709 if dist <= self.distance: 710 angle = distances.calc_angles(d.position, h.position, 711 a.position) 712 angle = np.rad2deg(angle) 713 if angle >= self.angle: 714 self.logger_debug( 715 "WB-D: {0!s} <-> S2-A: {1!s} {2:f} A, {3:f} DEG"\ 716 .format(h.index, a.index, dist, angle)) 717 s2_frame_results_dict[(h.resname, h.resid)].append( 718 (h.index, a.index, 719 (h.resname, h.resid, h.name), 720 (a.resname, a.resid, a.name), 721 dist, angle)) 722 723 if water_bridges_acceptors: 724 self.logger_debug("Selection 2 Donors <-> Selection 2 Acceptors") 725 ns_acceptors = AtomNeighborSearch(water_bridges_acceptors) 726 for i, donor_h_set in self._s2_donors_h.items(): 727 d = self._s2_donors[i] 728 for h in donor_h_set: 729 res = ns_acceptors.search(h, self.distance) 730 for a in res: 731 donor_atom = h if self.distance_type != 'heavy' else d 732 dist = distances.calc_bonds(donor_atom.position, 733 a.position) 734 if dist <= self.distance: 735 angle = distances.calc_angles(d.position, h.position, 736 a.position) 737 angle = np.rad2deg(angle) 738 if angle >= self.angle: 739 self.logger_debug( 740 "WB-A: {0!s} <-> S2-D: {1!s} {2:f} A, {3:f} DEG"\ 741 .format(a.index, h.index, dist, angle)) 742 s2_frame_results_dict[(a.resname, a.resid)].append( 743 (h.index, a.index, 744 (h.resname, h.resid, h.name), 745 (a.resname, a.resid, a.name), 746 dist, angle)) 747 748 # Generate the water network 749 water_network = {} 750 for key in s2_frame_results_dict: 751 s1_frame_results = set(s1_frame_results_dict[key]) 752 s2_frame_results = set(s2_frame_results_dict[key]) 753 if len(s1_frame_results.union(s2_frame_results)) > 1: 754 # Thus if selection 1 and selection 2 are the same and both 755 # only form a single hydrogen bond with a water, this entry 756 # won't be included. 757 water_network[key] = [s1_frame_results, 758 s2_frame_results.difference(s1_frame_results)] 759 # Generate frame_results 760 frame_results = [] 761 for s1_frame_results, s2_frame_results in water_network.values(): 762 frame_results.extend(list(s1_frame_results)) 763 frame_results.extend(list(s2_frame_results)) 764 765 self._timeseries.append(frame_results) 766 self._water_network.append(water_network) 767 768 769 logger.info("WBridge analysis: complete; timeseries %s.timeseries", 770 self.__class__.__name__) 771 @staticmethod 772 def _reformat_hb(hb, atomformat="{0[0]!s}{0[1]!s}:{0[2]!s}"): 773 """ 774 .. deprecated:: 1.0 775 This is a compatibility layer so the water bridge 776 _timeseries to timeserie conversion can perform as it does in the 777 hydrogen bond analysis. 778 779 For more information about the function of _reformat_hb in hydrogen 780 bond analysis, please refer to the _reformat_hb in the hydrogen 781 bond analysis module. 782 783 As the _timeseries to timeserie conversion will be deprecated in 1.0 784 this function will automatically lose its value. 785 """ 786 787 return (list(hb[:2]) 788 + [atomformat.format(hb[2]), atomformat.format(hb[3])] 789 + list(hb[4:])) 790 791 @property 792 def timeseries(self): 793 r'''Time series of water bridges. 794 795 Due to the intrinsic complexity of water bridge problem, where several 796 atoms :math:`n` can be linked by the same water. A simple ``atom 1 - 797 water bridge - atom 2`` mapping will create a very large list 798 :math:`n!/(2(n-2)!)` due to the rule of combination. 799 800 Thus, the output is arranged based on each individual bridging water 801 allowing the later reconstruction of the water network. The hydrogen 802 bonds from selection 1 to the **first** bridging water is followed by 803 the hydrogen bonds from selection 2 to the **same** bridging water. 804 After that, hydrogen bonds from selection 1 to the **second** bridging 805 water is succeeded by hydrogen bonds from selection 2 to the **same 806 second** bridging water. An example of the output is given in 807 :ref:`wb_Analysis_Output`. 808 809 Note 810 ---- 811 To find an acceptor atom in :attr:`Universe.atoms` by 812 *index* one would use ``u.atoms[acceptor_index]``. 813 814 The :attr:`timeseries` is a managed attribute and it is generated 815 from the underlying data in :attr:`_timeseries` every time the 816 attribute is accessed. It is therefore costly to call and if 817 :attr:`timeseries` is needed repeatedly it is recommended that you 818 assign to a variable:: 819 820 w = WaterBridgeAnalysis(u) 821 w.run() 822 timeseries = w.timeseries 823 824 See Also 825 -------- 826 :attr:`table` : structured array of the data 827 ''' 828 return super(WaterBridgeAnalysis, self).timeseries 829 830 def generate_table(self): 831 """Generate a normalised table of the results. 832 833 The table is stored as a :class:`numpy.recarray` in the 834 attribute :attr:`~WaterBridgeAnalysis.table`. 835 836 See Also 837 -------- 838 WaterBridgeAnalysis.table 839 """ 840 super(WaterBridgeAnalysis, self).generate_table() 841 842 def count_by_type(self): 843 """Counts the frequency of water bridge of a specific type. 844 845 If one atom *A* from *selection 1* is linked to atom *B* from 846 *selection 2* through a bridging water, an entity will be created and 847 the proportion of time that this linkage exists in the whole simulation 848 will be calculated. 849 850 Returns a :class:`numpy.recarray` containing atom indices for *A* and 851 *B*, residue names, residue numbers, atom names (for both A and B) and 852 the fraction of the total time during which the water bridge was 853 detected. This method returns None if method 854 :meth:`WaterBridgeAnalysis.run` was not executed first. 855 856 Returns 857 ------- 858 counts : numpy.recarray 859 Each row of the array contains data to define a unique water 860 bridge together with the frequency (fraction of the total time) 861 that it has been observed. 862 """ 863 if not self._has_timeseries(): 864 return 865 866 wbridges = defaultdict(int) 867 for wframe in self._water_network: 868 pairs = set([]) 869 for water, (selection1, selection2) in wframe.items(): 870 for (donor_index, acceptor_index, donor, acceptor, distance, 871 angle) in selection1: 872 if donor[:2] == water: 873 sele1 = acceptor 874 sele1_index = acceptor_index 875 else: 876 sele1 = donor 877 sele1_index = donor_index 878 sele1_resnm, sele1_resid, sele1_atom = sele1 879 880 for (donor_index, acceptor_index, donor, acceptor, 881 distance, angle) in selection2: 882 if donor[:2] == water: 883 sele2 = acceptor 884 sele2_index = acceptor_index 885 else: 886 sele2 = donor 887 sele2_index = donor_index 888 sele2_resnm, sele2_resid, sele2_atom = sele2 889 890 key = (sele1_index, sele2_index) 891 if key in pairs: 892 continue 893 else: 894 pairs.add(key) 895 wb_key = (sele1_index, sele2_index, sele1_resnm, 896 sele1_resid, sele1_atom, sele2_resnm, sele2_resid, 897 sele2_atom) 898 wbridges[wb_key] += 1 899 900 # build empty output table 901 dtype = [ 902 ("sele1_index", int), ("sele2_index", int), ('sele1_resnm', 'U4'), 903 ('sele1_resid', int), ('sele1_atom', 'U4'), ('sele2_resnm', 'U4'), 904 ('sele2_resid', int), ('sele2_atom', 'U4'), 905 ('frequency', float) 906 ] 907 out = np.empty((len(wbridges),), dtype=dtype) 908 909 tsteps = float(len(self.timesteps)) 910 for cursor, (key, count) in enumerate(six.iteritems(wbridges)): 911 out[cursor] = key + (count / tsteps,) 912 913 r = out.view(np.recarray) 914 return r 915