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 dynamics analysis --- :mod:`MDAnalysis.analysis.waterdynamics` 24======================================================================= 25 26:Author: Alejandro Bernardin 27:Year: 2014-2015 28:Copyright: GNU Public License v3 29 30.. versionadded:: 0.11.0 31 32This module provides functions to analize water dynamics trajectories and water 33interactions with other molecules. The functions in this module are: water 34orientational relaxation (WOR) [Yeh1999]_, hydrogen bond lifetimes (HBL) 35[Rapaport1983]_, angular distribution (AD) [Grigera1995]_, mean square 36displacement (MSD) [Brodka1994]_ and survival probability (SP) [Liu2004]_. 37 38For more information about this type of analysis please refer to 39[Araya-Secchi2014]_ (water in a protein cavity) and [Milischuk2011]_ (water in 40a nanopore). 41 42.. rubric:: References 43 44.. [Rapaport1983] D.C. Rapaport (1983): Hydrogen bonds in water, Molecular 45 Physics: An International Journal at the Interface Between 46 Chemistry and Physics, 50:5, 1151-1162. 47 48.. [Yeh1999] Yu-ling Yeh and Chung-Yuan Mou (1999). Orientational Relaxation 49 Dynamics of Liquid Water Studied by Molecular Dynamics Simulation, 50 J. Phys. Chem. B 1999, 103, 3699-3705. 51 52.. [Grigera1995] Raul Grigera, Susana G. Kalko and Jorge Fischbarg 53 (1995). Wall-Water Interface. A Molecular Dynamics Study, 54 Langmuir 1996,12,154-158 55 56.. [Liu2004] Pu Liu, Edward Harder, and B. J. Berne (2004).On the Calculation 57 of Diffusion Coefficients in Confined Fluids and Interfaces with 58 an Application to the Liquid-Vapor Interface of Water, 59 J. Phys. Chem. B 2004, 108, 6595-6602. 60 61.. [Brodka1994] Aleksander Brodka (1994). Diffusion in restricted volume, 62 Molecular Physics, 1994, Vol. 82, No. 5, 1075-1078. 63 64.. [Araya-Secchi2014] Araya-Secchi, R., Tomas Perez-Acle, Seung-gu Kang, Tien 65 Huynh, Alejandro Bernardin, Yerko Escalona, Jose-Antonio 66 Garate, Agustin D. Martinez, Isaac E. Garcia, Juan 67 C. Saez, Ruhong Zhou (2014). Characterization of a novel 68 water pocket inside the human Cx26 hemichannel 69 structure. Biophysical journal, 107(3), 599-612. 70 71.. [Milischuk2011] Anatoli A. Milischuk and Branka M. Ladanyi. Structure and 72 dynamics of water confined in silica 73 nanopores. J. Chem. Phys. 135, 174709 (2011); doi: 74 10.1063/1.3657408 75 76 77Example use of the analysis classes 78----------------------------------- 79 80HydrogenBondLifetimes 81~~~~~~~~~~~~~~~~~~~~~ 82 83Analyzing hydrogen bond lifetimes (HBL) :class:`HydrogenBondLifetimes`, both 84continuos and intermittent. In this case we are analyzing how residue 38 85interact with a water sphere of radius 6.0 centered on the geometric center of 86protein and residue 42. If the hydrogen bond lifetimes are very stable, we can 87assume that residue 38 is hydrophilic, on the other hand, if the are very 88unstable, we can assume that residue 38 is hydrophobic:: 89 90 import MDAnalysis 91 from MDAnalysis.analysis.waterdynamics import HydrogenBondLifetimes as HBL 92 93 u = MDAnalysis.Universe(pdb, trajectory) 94 selection1 = "byres name OH2 and sphzone 6.0 protein and resid 42" 95 selection2 = "resid 38" 96 HBL_analysis = HBL(universe, selection1, selection2, 0, 2000, 30) 97 HBL_analysis.run() 98 time = 0 99 #now we print the data ready to plot. The first two columns are the HBLc vs t 100 #plot and the second two columns are the HBLi vs t graph 101 for HBLc, HBLi in HBL_analysis.timeseries: 102 print("{time} {HBLc} {time} {HBLi}".format(time=time, HBLc=HBLc, HBLi=HBLi)) 103 time += 1 104 105 #we can also plot our data 106 plt.figure(1,figsize=(18, 6)) 107 108 #HBL continuos 109 plt.subplot(121) 110 plt.xlabel('time') 111 plt.ylabel('HBLc') 112 plt.title('HBL Continuos') 113 plt.plot(range(0,time),[column[0] for column in HBL_analysis.timeseries]) 114 115 #HBL intermitent 116 plt.subplot(122) 117 plt.xlabel('time') 118 plt.ylabel('HBLi') 119 plt.title('HBL Intermitent') 120 plt.plot(range(0,time),[column[1] for column in HBL_analysis.timeseries]) 121 122 plt.show() 123 124where HBLc is the value for the continuos hydrogen bond lifetimes and HBLi is 125the value for the intermittent hydrogen bond lifetime, t0 = 0, tf = 2000 and 126dtmax = 30. In this way we create 30 windows timestep (30 values in x 127axis). The continuos hydrogen bond lifetimes should decay faster than 128intermittent. 129 130 131WaterOrientationalRelaxation 132~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 133 134Analyzing water orientational relaxation (WOR) 135:class:`WaterOrientationalRelaxation`. In this case we are analyzing "how fast" 136water molecules are rotating/changing direction. If WOR is very stable we can 137assume that water molecules are rotating/changing direction very slow, on the 138other hand, if WOR decay very fast, we can assume that water molecules are 139rotating/changing direction very fast:: 140 141 import MDAnalysis 142 from MDAnalysis.analysis.waterdynamics import WaterOrientationalRelaxation as WOR 143 144 u = MDAnalysis.Universe(pdb, trajectory) 145 selection = "byres name OH2 and sphzone 6.0 protein and resid 42" 146 WOR_analysis = WOR(universe, selection, 0, 1000, 20) 147 WOR_analysis.run() 148 time = 0 149 #now we print the data ready to plot. The first two columns are WOR_OH vs t plot, 150 #the second two columns are WOR_HH vs t graph and the third two columns are WOR_dip vs t graph 151 for WOR_OH, WOR_HH, WOR_dip in WOR_analysis.timeseries: 152 print("{time} {WOR_OH} {time} {WOR_HH} {time} {WOR_dip}".format(time=time, WOR_OH=WOR_OH, WOR_HH=WOR_HH,WOR_dip=WOR_dip)) 153 time += 1 154 155 #now, if we want, we can plot our data 156 plt.figure(1,figsize=(18, 6)) 157 158 #WOR OH 159 plt.subplot(131) 160 plt.xlabel('time') 161 plt.ylabel('WOR') 162 plt.title('WOR OH') 163 plt.plot(range(0,time),[column[0] for column in WOR_analysis.timeseries]) 164 165 #WOR HH 166 plt.subplot(132) 167 plt.xlabel('time') 168 plt.ylabel('WOR') 169 plt.title('WOR HH') 170 plt.plot(range(0,time),[column[1] for column in WOR_analysis.timeseries]) 171 172 #WOR dip 173 plt.subplot(133) 174 plt.xlabel('time') 175 plt.ylabel('WOR') 176 plt.title('WOR dip') 177 plt.plot(range(0,time),[column[2] for column in WOR_analysis.timeseries]) 178 179 plt.show() 180 181where t0 = 0, tf = 1000 and dtmax = 20. In this way we create 20 windows 182timesteps (20 values in the x axis), the first window is created with 1000 183timestep average (1000/1), the second window is created with 500 timestep 184average(1000/2), the third window is created with 333 timestep average (1000/3) 185and so on. 186 187AngularDistribution 188~~~~~~~~~~~~~~~~~~~ 189 190Analyzing angular distribution (AD) :class:`AngularDistribution` for OH vector, 191HH vector and dipole vector. It returns a line histogram with vector 192orientation preference. A straight line in the output plot means no 193preferential orientation in water molecules. In this case we are analyzing if 194water molecules have some orientational preference, in this way we can see if 195water molecules are under an electric field or if they are interacting with 196something (residue, protein, etc):: 197 198 import MDAnalysis 199 from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD 200 201 u = MDAnalysis.Universe(pdb, trajectory) 202 selection = "byres name OH2 and sphzone 6.0 (protein and (resid 42 or resid 26) )" 203 bins = 30 204 AD_analysis = AD(universe,selection,bins) 205 AD_analysis.run() 206 #now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector , 207 #the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns 208 #are P(cos(theta)) vs cos(theta) for dipole vector 209 for bin in range(bins): 210 print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin])) 211 212 #and if we want to graph our results 213 plt.figure(1,figsize=(18, 6)) 214 215 #AD OH 216 plt.subplot(131) 217 plt.xlabel('cos theta') 218 plt.ylabel('P(cos theta)') 219 plt.title('PDF cos theta for OH') 220 plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]]) 221 222 #AD HH 223 plt.subplot(132) 224 plt.xlabel('cos theta') 225 plt.ylabel('P(cos theta)') 226 plt.title('PDF cos theta for HH') 227 plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]]) 228 229 #AD dip 230 plt.subplot(133) 231 plt.xlabel('cos theta') 232 plt.ylabel('P(cos theta)') 233 plt.title('PDF cos theta for dipole') 234 plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]]) 235 236 plt.show() 237 238 239where `P(cos(theta))` is the angular distribution or angular probabilities. 240 241 242MeanSquareDisplacement 243~~~~~~~~~~~~~~~~~~~~~~ 244 245Analyzing mean square displacement (MSD) :class:`MeanSquareDisplacement` for 246water molecules. In this case we are analyzing the average distance that water 247molecules travels inside protein in XYZ direction (cylindric zone of radius 24811[nm], Zmax 4.0[nm] and Zmin -8.0[nm]). A strong rise mean a fast movement of 249water molecules, a weak rise mean slow movement of particles:: 250 251 import MDAnalysis 252 from MDAnalysis.analysis.waterdynamics import MeanSquareDisplacement as MSD 253 254 u = MDAnalysis.Universe(pdb, trajectory) 255 selection = "byres name OH2 and cyzone 11.0 4.0 -8.0 protein" 256 MSD_analysis = MSD(universe, selection, 0, 1000, 20) 257 MSD_analysis.run() 258 #now we print data ready to graph. The graph 259 #represents MSD vs t 260 time = 0 261 for msd in MSD_analysis.timeseries: 262 print("{time} {msd}".format(time=time, msd=msd)) 263 time += 1 264 265 #Plot 266 plt.xlabel('time') 267 plt.ylabel('MSD') 268 plt.title('MSD') 269 plt.plot(range(0,time),MSD_analysis.timeseries) 270 plt.show() 271 272 273.. _SP-examples: 274 275SurvivalProbability 276~~~~~~~~~~~~~~~~~~~ 277 278Analyzing survival probability (SP) :class:`SurvivalProbability` for water 279molecules. In this case we are analyzing how long water molecules remain in a 280sphere of radius 12.3 centered in the geometrical center of resid 42, 26, 34 281and 80. A slow decay of SP means a long permanence time of water molecules in 282the zone, on the other hand, a fast decay means a short permanence time:: 283 284 import MDAnalysis 285 from MDAnalysis.analysis.waterdynamics import SurvivalProbability as SP 286 import matplotlib.pyplot as plt 287 288 universe = MDAnalysis.Universe(pdb, trajectory) 289 selection = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26 or resid 34 or resid 80) " 290 sp = SP(universe, selection, verbose=True) 291 sp.run(start=0, stop=100, tau_max=20) 292 tau_timeseries = sp.tau_timeseries 293 sp_timeseries = sp.sp_timeseries 294 295 # print in console 296 for tau, sp in zip(tau_timeseries, sp_timeseries): 297 print("{time} {sp}".format(time=tau, sp=sp)) 298 299 # plot 300 plt.xlabel('Time') 301 plt.ylabel('SP') 302 plt.title('Survival Probability') 303 plt.plot(taus, sp_timeseries) 304 plt.show() 305 306 307.. _Output: 308 309Output 310------ 311 312HydrogenBondLifetimes 313~~~~~~~~~~~~~~~~~~~~~ 314 315Hydrogen bond lifetimes (HBL) data is returned per window timestep, which is 316stored in :attr:`HydrogenBondLifetimes.timeseries` (in all the following 317descriptions, # indicates comments that are not part of the output):: 318 319 results = [ 320 [ # time t0 321 <HBL_c>, <HBL_i> 322 ], 323 [ # time t1 324 <HBL_c>, <HBL_i> 325 ], 326 ... 327 ] 328 329WaterOrientationalRelaxation 330~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 331 332Water orientational relaxation (WOR) data is returned per window timestep, 333which is stored in :attr:`WaterOrientationalRelaxation.timeseries`:: 334 335 results = [ 336 [ # time t0 337 <WOR_OH>, <WOR_HH>, <WOR_dip> 338 ], 339 [ # time t1 340 <WOR_OH>, <WOR_HH>, <WOR_dip> 341 ], 342 ... 343 ] 344 345AngularDistribution 346~~~~~~~~~~~~~~~~~~~ 347 348Angular distribution (AD) data is returned per vector, which is stored in 349:attr:`AngularDistribution.graph`. In fact, AngularDistribution returns a 350histogram:: 351 352 results = [ 353 [ # OH vector values 354 # the values are order in this way: <x_axis y_axis> 355 <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ... 356 ], 357 [ # HH vector values 358 <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ... 359 ], 360 [ # dip vector values 361 <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ... 362 ], 363 ] 364 365MeanSquareDisplacement 366~~~~~~~~~~~~~~~~~~~~~~ 367 368Mean Square Displacement (MSD) data is returned in a list, which each element 369represents a MSD value in its respective window timestep. Data is stored in 370:attr:`MeanSquareDisplacement.timeseries`:: 371 372 results = [ 373 #MSD values orders by window timestep 374 <MSD_t0>, <MSD_t1>, ... 375 ] 376 377SurvivalProbability 378~~~~~~~~~~~~~~~~~~~ 379 380Survival Probability (SP) computes two lists: a list of taus (:attr:`SurvivalProbability.tau_timeseries`) and a list of their corresponding mean survival 381probabilities (:attr:`SurvivalProbability.sp_timeseries`). Additionally, a list :attr:`SurvivalProbability.sp_timeseries_data` is provided which contains 382a list of SPs for each tau, which can be used to compute their distribution, etc. 383 384 results = [ tau1, tau2, ..., tau_n ], [ sp_tau1, sp_tau2, ..., sp_tau_n] 385 386Additionally, for each 387 388Classes 389-------- 390 391.. autoclass:: HydrogenBondLifetimes 392 :members: 393 :inherited-members: 394 395.. autoclass:: WaterOrientationalRelaxation 396 :members: 397 :inherited-members: 398 399.. autoclass:: AngularDistribution 400 :members: 401 :inherited-members: 402 403.. autoclass:: MeanSquareDisplacement 404 :members: 405 :inherited-members: 406 407.. autoclass:: SurvivalProbability 408 :members: 409 :inherited-members: 410 411""" 412from __future__ import print_function, division, absolute_import 413 414import warnings 415 416from six.moves import range, zip_longest 417 418import numpy as np 419import multiprocessing 420 421import MDAnalysis.analysis.hbonds 422from MDAnalysis.lib.log import ProgressMeter 423 424 425class HydrogenBondLifetimes(object): 426 r"""Hydrogen bond lifetime analysis 427 428 This is a autocorrelation function that gives the "Hydrogen Bond Lifetimes" 429 (HBL) proposed by D.C. Rapaport [Rapaport1983]_. From this function we can 430 obtain the continuous and intermittent behavior of hydrogen bonds in 431 time. A fast decay in these parameters indicate a fast change in HBs 432 connectivity. A slow decay indicate very stables hydrogen bonds, like in 433 ice. The HBL is also know as "Hydrogen Bond Population Relaxation" 434 (HBPR). In the continuos case we have: 435 436 .. math:: 437 C_{HB}^c(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h'_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)} 438 439 where :math:`h'_{ij}(t_0+\tau)=1` if there is a H-bond between a pair 440 :math:`ij` during time interval :math:`t_0+\tau` (continuos) and 441 :math:`h'_{ij}(t_0+\tau)=0` otherwise. In the intermittent case we have: 442 443 .. math:: 444 C_{HB}^i(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)} 445 446 where :math:`h_{ij}(t_0+\tau)=1` if there is a H-bond between a pair 447 :math:`ij` at time :math:`t_0+\tau` (intermittent) and 448 :math:`h_{ij}(t_0+\tau)=0` otherwise. 449 450 451 Parameters 452 ---------- 453 universe : Universe 454 Universe object 455 selection1 : str 456 Selection string for first selection [‘byres name OH2’]. 457 It could be any selection available in MDAnalysis, not just water. 458 selection2 : str 459 Selection string to analize its HBL against selection1 460 t0 : int 461 frame where analysis begins 462 tf : int 463 frame where analysis ends 464 dtmax : int 465 Maximum dt size, `dtmax` < `tf` or it will crash. 466 nproc : int 467 Number of processors to use, by default is 1. 468 469 470 .. versionadded:: 0.11.0 471 """ 472 473 def __init__(self, universe, selection1, selection2, t0, tf, dtmax, 474 nproc=1): 475 self.universe = universe 476 self.selection1 = selection1 477 self.selection2 = selection2 478 self.t0 = t0 479 self.tf = tf - 1 480 self.dtmax = dtmax 481 self.nproc = nproc 482 self.timeseries = None 483 484 def _getC_i(self, HBP, t0, t): 485 """ 486 This function give the intermitent Hydrogen Bond Lifetime 487 C_i = <h(t0)h(t)>/<h(t0)> between t0 and t 488 """ 489 C_i = 0 490 for i in range(len(HBP[t0])): 491 for j in range(len(HBP[t])): 492 if (HBP[t0][i][0] == HBP[t][j][0] and HBP[t0][i][1] == HBP[t][j][1]): 493 C_i += 1 494 break 495 if len(HBP[t0]) == 0: 496 return 0.0 497 else: 498 return float(C_i) / len(HBP[t0]) 499 500 def _getC_c(self, HBP, t0, t): 501 """ 502 This function give the continous Hydrogen Bond Lifetime 503 C_c = <h(t0)h'(t)>/<h(t0)> between t0 and t 504 """ 505 C_c = 0 506 dt = 1 507 begt0 = t0 508 HBP_cp = HBP 509 HBP_t0 = HBP[t0] 510 newHBP = [] 511 if t0 == t: 512 return 1.0 513 while t0 + dt <= t: 514 for i in range(len(HBP_t0)): 515 for j in range(len(HBP_cp[t0 + dt])): 516 if (HBP_t0[i][0] == HBP_cp[t0 + dt][j][0] and 517 HBP_t0[i][1] == HBP_cp[t0 + dt][j][1]): 518 newHBP.append(HBP_t0[i]) 519 break 520 C_c = len(newHBP) 521 t0 += dt 522 HBP_t0 = newHBP 523 newHBP = [] 524 if len(HBP[begt0]) == 0: 525 return 0 526 else: 527 return C_c / float(len(HBP[begt0])) 528 529 def _intervC_c(self, HBP, t0, tf, dt): 530 """ 531 This function gets all the data for the h(t0)h(t0+dt)', where 532 t0 = 1,2,3,...,tf. This function give us one point of the final plot 533 HBL vs t 534 """ 535 a = 0 536 count = 0 537 for i in range(len(HBP)): 538 if (t0 + dt <= tf): 539 if t0 == t0 + dt: 540 b = self._getC_c(HBP, t0, t0) 541 break 542 b = self._getC_c(HBP, t0, t0 + dt) 543 t0 += dt 544 a += b 545 count += 1 546 if count == 0: 547 return 1.0 548 return a / count 549 550 def _intervC_i(self, HBP, t0, tf, dt): 551 """ 552 This function gets all the data for the h(t0)h(t0+dt), where 553 t0 = 1,2,3,...,tf. This function give us a point of the final plot 554 HBL vs t 555 """ 556 a = 0 557 count = 0 558 for i in range(len(HBP)): 559 if (t0 + dt <= tf): 560 b = self._getC_i(HBP, t0, t0 + dt) 561 t0 += dt 562 a += b 563 count += 1 564 return a / count 565 566 def _finalGraphGetC_i(self, HBP, t0, tf, maxdt): 567 """ 568 This function gets the final data of the C_i graph. 569 """ 570 output = [] 571 for dt in range(maxdt): 572 a = self._intervC_i(HBP, t0, tf, dt) 573 output.append(a) 574 return output 575 576 def _finalGraphGetC_c(self, HBP, t0, tf, maxdt): 577 """ 578 This function gets the final data of the C_c graph. 579 """ 580 output = [] 581 for dt in range(maxdt): 582 a = self._intervC_c(HBP, t0, tf, dt) 583 output.append(a) 584 return output 585 586 def _getGraphics(self, HBP, t0, tf, maxdt): 587 """ 588 Function that join all the results into a plot. 589 """ 590 a = [] 591 cont = self._finalGraphGetC_c(HBP, t0, tf, maxdt) 592 inte = self._finalGraphGetC_i(HBP, t0, tf, maxdt) 593 for i in range(len(cont)): 594 fix = [cont[i], inte[i]] 595 a.append(fix) 596 return a 597 598 def _HBA(self, ts, conn, universe, selAtom1, selAtom2, 599 verbose=False): 600 """ 601 Main function for calculate C_i and C_c in parallel. 602 """ 603 finalGetResidue1 = selAtom1 604 finalGetResidue2 = selAtom2 605 frame = ts.frame 606 h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(universe, 607 finalGetResidue1, 608 finalGetResidue2, 609 distance=3.5, 610 angle=120.0, 611 start=frame - 1, 612 stop=frame) 613 while True: 614 try: 615 h.run(verbose=verbose) 616 break 617 except: 618 print("error") 619 print("trying again") 620 sys.stdout.flush() 621 sys.stdout.flush() 622 conn.send(h.timeseries[0]) 623 conn.close() 624 625 def run(self, **kwargs): 626 """Analyze trajectory and produce timeseries""" 627 h_list = [] 628 i = 0 629 if (self.nproc > 1): 630 while i < len(self.universe.trajectory): 631 jobs = [] 632 k = i 633 for j in range(self.nproc): 634 # start 635 print("ts=", i + 1) 636 if i >= len(self.universe.trajectory): 637 break 638 conn_parent, conn_child = multiprocessing.Pipe(False) 639 while True: 640 try: 641 # new thread 642 jobs.append( 643 (multiprocessing.Process( 644 target=self._HBA, 645 args=(self.universe.trajectory[i], 646 conn_child, self.universe, 647 self.selection1, self.selection2,)), 648 conn_parent)) 649 break 650 except: 651 print("error in jobs.append") 652 jobs[j][0].start() 653 i = i + 1 654 655 for j in range(self.nproc): 656 if k >= len(self.universe.trajectory): 657 break 658 rec01 = jobs[j][1] 659 received = rec01.recv() 660 h_list.append(received) 661 jobs[j][0].join() 662 k += 1 663 self.timeseries = self._getGraphics( 664 h_list, 0, self.tf - 1, self.dtmax) 665 else: 666 h_list = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(self.universe, 667 self.selection1, 668 self.selection2, 669 distance=3.5, 670 angle=120.0) 671 h_list.run(**kwargs) 672 self.timeseries = self._getGraphics( 673 h_list.timeseries, self.t0, self.tf, self.dtmax) 674 675 676class WaterOrientationalRelaxation(object): 677 r"""Water orientation relaxation analysis 678 679 Function to evaluate the Water Orientational Relaxation proposed by Yu-ling 680 Yeh and Chung-Yuan Mou [Yeh1999_]. WaterOrientationalRelaxation indicates 681 "how fast" water molecules are rotating or changing direction. This is a 682 time correlation function given by: 683 684 .. math:: 685 C_{\hat u}(\tau)=\langle \mathit{P}_2[\mathbf{\hat{u}}(t_0)\cdot\mathbf{\hat{u}}(t_0+\tau)]\rangle 686 687 where :math:`P_2=(3x^2-1)/2` is the second-order Legendre polynomial and :math:`\hat{u}` is 688 a unit vector along HH, OH or dipole vector. 689 690 691 Parameters 692 ---------- 693 universe : Universe 694 Universe object 695 selection : str 696 Selection string for water [‘byres name OH2’]. 697 t0 : int 698 frame where analysis begins 699 tf : int 700 frame where analysis ends 701 dtmax : int 702 Maximum dt size, `dtmax` < `tf` or it will crash. 703 704 705 .. versionadded:: 0.11.0 706 707 """ 708 709 def __init__(self, universe, selection, t0, tf, dtmax, nproc=1): 710 self.universe = universe 711 self.selection = selection 712 self.t0 = t0 713 self.tf = tf 714 self.dtmax = dtmax 715 self.nproc = nproc 716 self.timeseries = None 717 718 def _repeatedIndex(self, selection, dt, totalFrames): 719 """ 720 Indicates the comparation between all the t+dt. 721 The results is a list of list with all the repeated index per frame 722 (or time). 723 Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... 724 Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... 725 Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... 726 """ 727 rep = [] 728 for i in range(int(round((totalFrames - 1) / float(dt)))): 729 if (dt * i + dt < totalFrames): 730 rep.append(self._sameMolecTandDT( 731 selection, dt * i, (dt * i) + dt)) 732 return rep 733 734 def _getOneDeltaPoint(self, universe, repInd, i, t0, dt): 735 """ 736 Gives one point to calculate the mean and gets one point of the plot 737 C_vect vs t. 738 Ex: t0=1 and tau=1 so calculate the t0-tau=1-2 intervale. 739 Ex: t0=5 and tau=3 so calcultate the t0-tau=5-8 intervale. 740 i = come from getMeanOnePoint (named j) (int) 741 """ 742 valOH = 0 743 valHH = 0 744 valdip = 0 745 n = 0 746 for j in range(len(repInd[i]) // 3): 747 begj = 3 * j 748 universe.trajectory[t0] 749 Ot0 = repInd[i][begj] 750 H1t0 = repInd[i][begj + 1] 751 H2t0 = repInd[i][begj + 2] 752 OHVector0 = H1t0.position - Ot0.position 753 HHVector0 = H1t0.position - H2t0.position 754 dipVector0 = ((H1t0.position + H2t0.position) * 0.5) - Ot0.position 755 756 universe.trajectory[t0 + dt] 757 Otp = repInd[i][begj] 758 H1tp = repInd[i][begj + 1] 759 H2tp = repInd[i][begj + 2] 760 761 OHVectorp = H1tp.position - Otp.position 762 HHVectorp = H1tp.position - H2tp.position 763 dipVectorp = ((H1tp.position + H2tp.position) * 0.5) - Otp.position 764 765 normOHVector0 = np.linalg.norm(OHVector0) 766 normOHVectorp = np.linalg.norm(OHVectorp) 767 normHHVector0 = np.linalg.norm(HHVector0) 768 normHHVectorp = np.linalg.norm(HHVectorp) 769 normdipVector0 = np.linalg.norm(dipVector0) 770 normdipVectorp = np.linalg.norm(dipVectorp) 771 772 unitOHVector0 = [OHVector0[0] / normOHVector0, 773 OHVector0[1] / normOHVector0, 774 OHVector0[2] / normOHVector0] 775 unitOHVectorp = [OHVectorp[0] / normOHVectorp, 776 OHVectorp[1] / normOHVectorp, 777 OHVectorp[2] / normOHVectorp] 778 unitHHVector0 = [HHVector0[0] / normHHVector0, 779 HHVector0[1] / normHHVector0, 780 HHVector0[2] / normHHVector0] 781 unitHHVectorp = [HHVectorp[0] / normHHVectorp, 782 HHVectorp[1] / normHHVectorp, 783 HHVectorp[2] / normHHVectorp] 784 unitdipVector0 = [dipVector0[0] / normdipVector0, 785 dipVector0[1] / normdipVector0, 786 dipVector0[2] / normdipVector0] 787 unitdipVectorp = [dipVectorp[0] / normdipVectorp, 788 dipVectorp[1] / normdipVectorp, 789 dipVectorp[2] / normdipVectorp] 790 791 valOH += self.lg2(np.dot(unitOHVector0, unitOHVectorp)) 792 valHH += self.lg2(np.dot(unitHHVector0, unitHHVectorp)) 793 valdip += self.lg2(np.dot(unitdipVector0, unitdipVectorp)) 794 n += 1 795 return (valOH/n, valHH/n, valdip/n) if n > 0 else (0, 0, 0) 796 797 798 def _getMeanOnePoint(self, universe, selection1, selection_str, dt, 799 totalFrames): 800 """ 801 This function gets one point of the plot C_vec vs t. It uses the 802 _getOneDeltaPoint() function to calculate the average. 803 804 """ 805 repInd = self._repeatedIndex(selection1, dt, totalFrames) 806 sumsdt = 0 807 n = 0.0 808 sumDeltaOH = 0.0 809 sumDeltaHH = 0.0 810 sumDeltadip = 0.0 811 812 for j in range(totalFrames // dt - 1): 813 a = self._getOneDeltaPoint(universe, repInd, j, sumsdt, dt) 814 sumDeltaOH += a[0] 815 sumDeltaHH += a[1] 816 sumDeltadip += a[2] 817 sumsdt += dt 818 n += 1 819 820 # if no water molecules remain in selection, there is nothing to get 821 # the mean, so n = 0. 822 return (sumDeltaOH / n, sumDeltaHH / n, sumDeltadip / n) if n > 0 else (0, 0, 0) 823 824 def _sameMolecTandDT(self, selection, t0d, tf): 825 """ 826 Compare the molecules in the t0d selection and the t0d+dt selection and 827 select only the particles that are repeated in both frame. This is to 828 consider only the molecules that remains in the selection after the dt 829 time has elapsed. 830 The result is a list with the indexs of the atoms. 831 """ 832 a = set(selection[t0d]) 833 b = set(selection[tf]) 834 sort = sorted(list(a.intersection(b))) 835 return sort 836 837 def _selection_serial(self, universe, selection_str): 838 selection = [] 839 pm = ProgressMeter(universe.trajectory.n_frames, 840 interval=10, verbose=True) 841 for ts in universe.trajectory: 842 selection.append(universe.select_atoms(selection_str)) 843 pm.echo(ts.frame) 844 return selection 845 846 @staticmethod 847 def lg2(x): 848 """Second Legendre polynomial""" 849 return (3*x*x - 1)/2 850 851 def run(self, **kwargs): 852 """Analyze trajectory and produce timeseries""" 853 854 # All the selection to an array, this way is faster than selecting 855 # later. 856 if self.nproc == 1: 857 selection_out = self._selection_serial( 858 self.universe, self.selection) 859 else: 860 # selection_out = self._selection_parallel(self.universe, 861 # self.selection, self.nproc) 862 # parallel selection to be implemented 863 selection_out = self._selection_serial( 864 self.universe, self.selection) 865 self.timeseries = [] 866 for dt in list(range(1, self.dtmax + 1)): 867 output = self._getMeanOnePoint( 868 self.universe, selection_out, self.selection, dt, self.tf) 869 self.timeseries.append(output) 870 871 872class AngularDistribution(object): 873 r"""Angular distribution function analysis 874 875 The angular distribution function (AD) is defined as the distribution 876 probability of the cosine of the :math:`\theta` angle formed by the OH 877 vector, HH vector or dipolar vector of water molecules and a vector 878 :math:`\hat n` parallel to chosen axis (z is the default value). The cosine 879 is define as :math:`\cos \theta = \hat u \cdot \hat n`, where :math:`\hat 880 u` is OH, HH or dipole vector. It creates a histogram and returns a list 881 of lists, see Output_. The AD is also know as Angular Probability (AP). 882 883 884 Parameters 885 ---------- 886 universe : Universe 887 Universe object 888 selection : str 889 Selection string to evaluate its angular distribution ['byres name OH2'] 890 bins : int (optional) 891 Number of bins to create the histogram by means of :func:`numpy.histogram` 892 axis : {'x', 'y', 'z'} (optional) 893 Axis to create angle with the vector (HH, OH or dipole) and calculate 894 cosine theta ['z']. 895 896 897 .. versionadded:: 0.11.0 898 """ 899 900 def __init__(self, universe, selection_str, bins=40, nproc=1, axis="z"): 901 self.universe = universe 902 self.selection_str = selection_str 903 self.bins = bins 904 self.nproc = nproc 905 self.axis = axis 906 self.graph = None 907 908 def _getCosTheta(self, universe, selection, axis): 909 valOH = [] 910 valHH = [] 911 valdip = [] 912 913 i = 0 914 while i <= (len(selection) - 1): 915 universe.trajectory[i] 916 line = selection[i].positions 917 918 Ot0 = line[::3] 919 H1t0 = line[1::3] 920 H2t0 = line[2::3] 921 922 OHVector0 = H1t0 - Ot0 923 HHVector0 = H1t0 - H2t0 924 dipVector0 = (H1t0 + H2t0) * 0.5 - Ot0 925 926 unitOHVector0 = OHVector0 / \ 927 np.linalg.norm(OHVector0, axis=1)[:, None] 928 unitHHVector0 = HHVector0 / \ 929 np.linalg.norm(HHVector0, axis=1)[:, None] 930 unitdipVector0 = dipVector0 / \ 931 np.linalg.norm(dipVector0, axis=1)[:, None] 932 933 j = 0 934 while j < len(line) / 3: 935 if axis == "z": 936 valOH.append(unitOHVector0[j][2]) 937 valHH.append(unitHHVector0[j][2]) 938 valdip.append(unitdipVector0[j][2]) 939 940 elif axis == "x": 941 valOH.append(unitOHVector0[j][0]) 942 valHH.append(unitHHVector0[j][0]) 943 valdip.append(unitdipVector0[j][0]) 944 945 elif axis == "y": 946 valOH.append(unitOHVector0[j][1]) 947 valHH.append(unitHHVector0[j][1]) 948 valdip.append(unitdipVector0[j][1]) 949 950 j += 1 951 i += 1 952 return (valOH, valHH, valdip) 953 954 def _getHistogram(self, universe, selection, bins, axis): 955 """ 956 This function gets a normalized histogram of the cos(theta) values. It 957 return a list of list. 958 """ 959 a = self._getCosTheta(universe, selection, axis) 960 cosThetaOH = a[0] 961 cosThetaHH = a[1] 962 cosThetadip = a[2] 963 lencosThetaOH = len(cosThetaOH) 964 lencosThetaHH = len(cosThetaHH) 965 lencosThetadip = len(cosThetadip) 966 histInterval = bins 967 histcosThetaOH = np.histogram(cosThetaOH, histInterval, density=True) 968 histcosThetaHH = np.histogram(cosThetaHH, histInterval, density=True) 969 histcosThetadip = np.histogram(cosThetadip, histInterval, density=True) 970 971 return (histcosThetaOH, histcosThetaHH, histcosThetadip) 972 973 def _hist2column(self, aList): 974 """ 975 This function transform from the histogram format 976 to a column format. 977 """ 978 a = [] 979 for x in zip_longest(*aList, fillvalue="."): 980 a.append(" ".join(str(i) for i in x)) 981 return a 982 983 def run(self, **kwargs): 984 """Function to evaluate the angular distribution of cos(theta)""" 985 986 if self.nproc == 1: 987 selection = self._selection_serial( 988 self.universe, self.selection_str) 989 else: 990 # not implemented yet 991 # selection = self._selection_parallel(self.universe, 992 # self.selection_str,self.nproc) 993 selection = self._selection_serial( 994 self.universe, self.selection_str) 995 996 self.graph = [] 997 output = self._getHistogram( 998 self.universe, selection, self.bins, self.axis) 999 # this is to format the exit of the file 1000 # maybe this output could be improved 1001 listOH = [list(output[0][1]), list(output[0][0])] 1002 listHH = [list(output[1][1]), list(output[1][0])] 1003 listdip = [list(output[2][1]), list(output[2][0])] 1004 1005 self.graph.append(self._hist2column(listOH)) 1006 self.graph.append(self._hist2column(listHH)) 1007 self.graph.append(self._hist2column(listdip)) 1008 1009 def _selection_serial(self, universe, selection_str): 1010 selection = [] 1011 pm = ProgressMeter(universe.trajectory.n_frames, 1012 interval=10, verbose=True) 1013 for ts in universe.trajectory: 1014 selection.append(universe.select_atoms(selection_str)) 1015 pm.echo(ts.frame) 1016 return selection 1017 1018 1019class MeanSquareDisplacement(object): 1020 r"""Mean square displacement analysis 1021 1022 Function to evaluate the Mean Square Displacement (MSD_). The MSD gives the 1023 average distance that particles travels. The MSD is given by: 1024 1025 .. math:: 1026 \langle\Delta r(t)^2\rangle = 2nDt 1027 1028 where :math:`r(t)` is the position of particle in time :math:`t`, 1029 :math:`\Delta r(t)` is the displacement after time lag :math:`t`, 1030 :math:`n` is the dimensionality, in this case :math:`n=3`, 1031 :math:`D` is the diffusion coefficient and :math:`t` is the time. 1032 1033 .. _MSD: http://en.wikipedia.org/wiki/Mean_squared_displacement 1034 1035 1036 Parameters 1037 ---------- 1038 universe : Universe 1039 Universe object 1040 selection : str 1041 Selection string for water [‘byres name OH2’]. 1042 t0 : int 1043 frame where analysis begins 1044 tf : int 1045 frame where analysis ends 1046 dtmax : int 1047 Maximum dt size, `dtmax` < `tf` or it will crash. 1048 1049 1050 .. versionadded:: 0.11.0 1051 """ 1052 1053 def __init__(self, universe, selection, t0, tf, dtmax, nproc=1): 1054 self.universe = universe 1055 self.selection = selection 1056 self.t0 = t0 1057 self.tf = tf 1058 self.dtmax = dtmax 1059 self.nproc = nproc 1060 self.timeseries = None 1061 1062 def _repeatedIndex(self, selection, dt, totalFrames): 1063 """ 1064 Indicate the comparation between all the t+dt. 1065 The results is a list of list with all the repeated index per frame 1066 (or time). 1067 1068 - Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... 1069 - Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... 1070 - Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... 1071 """ 1072 rep = [] 1073 for i in range(int(round((totalFrames - 1) / float(dt)))): 1074 if (dt * i + dt < totalFrames): 1075 rep.append(self._sameMolecTandDT( 1076 selection, dt * i, (dt * i) + dt)) 1077 return rep 1078 1079 def _getOneDeltaPoint(self, universe, repInd, i, t0, dt): 1080 """ 1081 Gives one point to calculate the mean and gets one point of the plot 1082 C_vect vs t. 1083 1084 - Ex: t0=1 and dt=1 so calculate the t0-dt=1-2 interval. 1085 - Ex: t0=5 and dt=3 so calcultate the t0-dt=5-8 interva 1086 1087 i = come from getMeanOnePoint (named j) (int) 1088 """ 1089 valO = 0 1090 n = 0 1091 for j in range(len(repInd[i]) // 3): 1092 begj = 3 * j 1093 universe.trajectory[t0] 1094 # Plus zero is to avoid 0to be equal to 0tp 1095 Ot0 = repInd[i][begj].position + 0 1096 1097 universe.trajectory[t0 + dt] 1098 # Plus zero is to avoid 0to be equal to 0tp 1099 Otp = repInd[i][begj].position + 0 1100 1101 # position oxygen 1102 OVector = Ot0 - Otp 1103 # here it is the difference with 1104 # waterdynamics.WaterOrientationalRelaxation 1105 valO += np.dot(OVector, OVector) 1106 n += 1 1107 1108 # if no water molecules remain in selection, there is nothing to get 1109 # the mean, so n = 0. 1110 return valO/n if n > 0 else 0 1111 1112 def _getMeanOnePoint(self, universe, selection1, selection_str, dt, 1113 totalFrames): 1114 """ 1115 This function gets one point of the plot C_vec vs t. It's uses the 1116 _getOneDeltaPoint() function to calculate the average. 1117 1118 """ 1119 repInd = self._repeatedIndex(selection1, dt, totalFrames) 1120 sumsdt = 0 1121 n = 0.0 1122 sumDeltaO = 0.0 1123 valOList = [] 1124 1125 for j in range(totalFrames // dt - 1): 1126 a = self._getOneDeltaPoint(universe, repInd, j, sumsdt, dt) 1127 sumDeltaO += a 1128 valOList.append(a) 1129 sumsdt += dt 1130 n += 1 1131 1132 # if no water molecules remain in selection, there is nothing to get 1133 # the mean, so n = 0. 1134 return sumDeltaO/n if n > 0 else 0 1135 1136 def _sameMolecTandDT(self, selection, t0d, tf): 1137 """ 1138 Compare the molecules in the t0d selection and the t0d+dt selection and 1139 select only the particles that are repeated in both frame. This is to 1140 consider only the molecules that remains in the selection after the dt 1141 time has elapsed. The result is a list with the indexs of the atoms. 1142 """ 1143 a = set(selection[t0d]) 1144 b = set(selection[tf]) 1145 sort = sorted(list(a.intersection(b))) 1146 return sort 1147 1148 def _selection_serial(self, universe, selection_str): 1149 selection = [] 1150 pm = ProgressMeter(universe.trajectory.n_frames, 1151 interval=10, verbose=True) 1152 for ts in universe.trajectory: 1153 selection.append(universe.select_atoms(selection_str)) 1154 pm.echo(ts.frame) 1155 return selection 1156 1157 def run(self, **kwargs): 1158 """Analyze trajectory and produce timeseries""" 1159 1160 # All the selection to an array, this way is faster than selecting 1161 # later. 1162 if self.nproc == 1: 1163 selection_out = self._selection_serial( 1164 self.universe, self.selection) 1165 else: 1166 # parallel not yet implemented 1167 # selection = selection_parallel(universe, selection_str, nproc) 1168 selection_out = self._selection_serial( 1169 self.universe, self.selection) 1170 self.timeseries = [] 1171 for dt in list(range(1, self.dtmax + 1)): 1172 output = self._getMeanOnePoint( 1173 self.universe, selection_out, self.selection, dt, self.tf) 1174 self.timeseries.append(output) 1175 1176 1177class SurvivalProbability(object): 1178 r"""Survival probability analysis 1179 1180 Function to evaluate the Survival Probability (SP). The SP gives the 1181 probability for a group of particles to remain in certain region. The SP is 1182 given by: 1183 1184 .. math:: 1185 P(\tau) = \frac1T \sum_{t=1}^T \frac{N(t,t+\tau)}{N(t)} 1186 1187 where :math:`T` is the maximum time of simulation, :math:`\tau` is the 1188 timestep, :math:`N(t)` the number of particles at time t, and 1189 :math:`N(t, t+\tau)` is the number of particles at every frame from t to `\tau`. 1190 1191 1192 Parameters 1193 ---------- 1194 universe : Universe 1195 Universe object 1196 selection : str 1197 Selection string; any selection is allowed. With this selection you 1198 define the region/zone where to analyze, e.g.: "resname SOL and around 5 (resname LIPID)" 1199 and "resname ION and around 10 (resid 20)" (see `SP-examples`_ ) 1200 verbose : Boolean 1201 If True, prints progress and comments to the console. 1202 1203 1204 .. versionadded:: 0.11.0 1205 1206 """ 1207 1208 def __init__(self, universe, selection, t0=None, tf=None, dtmax=None, verbose=False): 1209 self.universe = universe 1210 self.selection = selection 1211 self.verbose = verbose 1212 1213 # backward compatibility 1214 self.start = self.stop = self.tau_max = None 1215 if t0 is not None: 1216 self.start = t0 1217 warnings.warn("t0 is deprecated, use run(start=t0) instead", category=DeprecationWarning) 1218 1219 if tf is not None: 1220 self.stop = tf 1221 warnings.warn("tf is deprecated, use run(stop=tf) instead", category=DeprecationWarning) 1222 1223 if dtmax is not None: 1224 self.tau_max = dtmax 1225 warnings.warn("dtmax is deprecated, use run(tau_max=dtmax) instead", category=DeprecationWarning) 1226 1227 def print(self, verbose, *args): 1228 if self.verbose: 1229 print(args) 1230 elif verbose: 1231 print(args) 1232 1233 def run(self, tau_max=20, start=0, stop=None, step=1, verbose=False): 1234 """ 1235 Computes and returns the survival probability timeseries 1236 1237 Parameters 1238 ---------- 1239 start : int 1240 Zero-based index of the first frame to be analysed 1241 stop : int 1242 Zero-based index of the last frame to be analysed (inclusive) 1243 step : int 1244 Jump every `step`'th frame 1245 tau_max : int 1246 Survival probability is calculated for the range :math:`1 <= \tau <= tau_max` 1247 verbose : Boolean 1248 Overwrite the constructor's verbosity 1249 1250 Returns 1251 ------- 1252 tau_timeseries : list 1253 tau from 1 to tau_max. Saved in the field tau_timeseries. 1254 sp_timeseries : list 1255 survival probability for each value of `tau`. Saved in the field sp_timeseries. 1256 """ 1257 1258 # backward compatibility (and priority) 1259 start = self.start if self.start is not None else start 1260 stop = self.stop if self.stop is not None else stop 1261 tau_max = self.tau_max if self.tau_max is not None else tau_max 1262 1263 # sanity checks 1264 if stop is not None and stop >= len(self.universe.trajectory): 1265 raise ValueError("\"stop\" must be smaller than the number of frames in the trajectory.") 1266 1267 if stop is None: 1268 stop = len(self.universe.trajectory) 1269 else: 1270 stop = stop + 1 1271 1272 if tau_max > (stop - start): 1273 raise ValueError("Too few frames selected for given tau_max.") 1274 1275 # load all frames to an array of sets 1276 selected_ids = [] 1277 for ts in self.universe.trajectory[start:stop]: 1278 self.print(verbose, "Loading frame:", ts) 1279 selected_ids.append(set(self.universe.select_atoms(self.selection).ids)) 1280 1281 tau_timeseries = np.arange(1, tau_max + 1) 1282 sp_timeseries_data = [[] for _ in range(tau_max)] 1283 1284 for t in range(0, len(selected_ids), step): 1285 Nt = len(selected_ids[t]) 1286 1287 if Nt == 0: 1288 self.print(verbose, 1289 "At frame {} the selection did not find any molecule. Moving on to the next frame".format(t)) 1290 continue 1291 1292 for tau in tau_timeseries: 1293 if t + tau >= len(selected_ids): 1294 break 1295 1296 # ids that survive from t to t + tau and at every frame in between 1297 Ntau = len(set.intersection(*selected_ids[t:t + tau + 1])) 1298 sp_timeseries_data[tau - 1].append(Ntau / float(Nt)) 1299 1300 # user can investigate the distribution and sample size 1301 self.sp_timeseries_data = sp_timeseries_data 1302 1303 self.tau_timeseries = tau_timeseries 1304 self.sp_timeseries = [np.mean(sp) for sp in sp_timeseries_data] 1305 return self 1306