From 35ffd18abe9147900418a9bbbec93c36394d404f Mon Sep 17 00:00:00 2001 From: Alejandro Bernardin Date: Tue, 2 Jun 2015 14:41:28 -0300 Subject: [PATCH 1/4] New module to analize water dynamics. This module have five new analysis methods: Hydrogen Bond Lifetime (HBL) Water Orientation Relaxation (WOR) Angular Distribution (AD) Mean Square Displacement (MSD) Surivival Probability (SP) --- package/MDAnalysis/analysis/waterdynamics.py | 1057 +++++++++++++++++ .../analysis/waterdynamics.rst | 2 + .../documentation_pages/analysis_modules.rst | 1 + testsuite/MDAnalysisTests/data/watdyn.dcd | Bin 0 -> 2876 bytes testsuite/MDAnalysisTests/data/watdyn.psf | 54 + testsuite/MDAnalysisTests/datafiles.py | 4 + testsuite/MDAnalysisTests/test_analysis.py | 38 +- 7 files changed, 1155 insertions(+), 1 deletion(-) create mode 100644 package/MDAnalysis/analysis/waterdynamics.py create mode 100644 package/doc/sphinx/source/documentation_pages/analysis/waterdynamics.rst create mode 100644 testsuite/MDAnalysisTests/data/watdyn.dcd create mode 100644 testsuite/MDAnalysisTests/data/watdyn.psf diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py new file mode 100644 index 00000000000..ef0d021688b --- /dev/null +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -0,0 +1,1057 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; encoding: utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- http://mdanalysis.googlecode.com +# Copyright (c) 2006-2012 Naveen Michaud-Agrawal, +# Elizabeth J. Denning, Oliver Beckstein, +# and contributors (see website for details) +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and +# O. Beckstein. MDAnalysis: A Toolkit for the Analysis of +# Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319--2327, +# doi:10.1002/jcc.21787 +# + +""" +Water dynamics analysis --- :mod:`MDAnalysis.analysis.waterdynamics` +======================================================================= + +:Author: Alejandro Bernardin +:Year: 2014-2015 +:Copyright: GNU Public License v2 + +This module provides functions to analize water dynamics trajectories and water interactions with other molecules. +The functions in this module are: water orientational relaxation (WOR) [Yeh1999]_, hydrogen bond lifetimes (HBL) [Rapaport1983]_, +angular distribution (AD) [Grigera1995]_, mean square displacement (MSD) [Brodka1994]_ and survival probability (SP) [Liu2004]_. + +For more information about this type of analysis please refer to [Araya-Secchi2014]_ (water in a protein cavity) and [Milischuk2011]_ (water in a nanopore). + +.. rubric:: References + +.. [Rapaport1983] D.C. Rapaport (1983): Hydrogen bonds in water, Molecular Physics: An International + Journal at the Interface Between Chemistry and Physics, 50:5, 1151-1162. + +.. [Yeh1999] Yu-ling Yeh and Chung-Yuan Mou (1999). + Orientational Relaxation Dynamics of Liquid Water Studied by Molecular Dynamics + Simulation, J. Phys. Chem. B 1999, 103, 3699-3705. + +.. [Grigera1995] Raul Grigera, Susana G. Kalko and Jorge Fischbarg (1995). Wall-Water Interface. + A Molecular Dynamics Study, Langmuir 1996,12,154-158 + +.. [Liu2004] Pu Liu, Edward Harder, and B. J. Berne (2004).On the Calculation of Diffusion Coefficients + in Confined Fluids and Interfaces with an Application to the Liquid-Vapor Interface of + Water, J. Phys. Chem. B 2004, 108, 6595-6602. + +.. [Brodka1994] Aleksander Brodka (1994). Diffusion in restricted volume, Molecular Physics, 1994, Vol. + 82, No. 5, 1075-1078. + +.. [Araya-Secchi2014] Raul Araya-Secchi, Seung-gu Kang, Tien Huynh, Alejandro Bernardin, + Agustin D. Martinez, Juan C. Saez, Tomas Perez-Acle, Ruhong Zhou. Characterization of a + novel water pocket inside the human Cx26 hemi-channel structure. +.. [Milischuk2011] Anatoli A. Milischuk and Branka M. Ladanyi. Structure and dynamics of water confined + in silica nanopores. J. Chem. Phys. 135, 174709 (2011); doi: 10.1063/1.3657408 + + + +.. examples + +Examples +-------- + +HBL +~~~~ + +Analyzing hydrogen bond lifetimes :class:`HBL`, both continuos and intermittent. In this case we are analyzing +how residue 38 interact with a water sphere of radius 6.0 centered on the geometric center of protein and +residue 42. If the HBL are very stable, we can assume that residue 38 is hydrophilic, on the other +hand, if the HBL are very unstable, we can assume that residue 38 is hydrophobic:: + + import MDAnalysis + from MDAnalysis.analysis.waterdynamics import HBL + + u = MDAnalysis.Universe(pdb, trajectory) + selection1 = "byres name OH2 and sphzone 6.0 protein and resid 42" + selection2 = "resid 38" + HBL_analysis = HBL(universe, selection1, selection2, 0, 2000, 30) + HBL_analysis.run() + i=0 + #now we print the data ready to graph. The first two columns are the HBLc vs t graph and + #the second two columns are the HBLi vs t graph + for HBLc, HBLi in HBL_analysis.timeseries: + print i +" "+ HBLc +" "+ i +" "+ HBLi + i+=1 + +where HBLc is the value for the continuos hydrogen bond lifetime and HBLi is the value for the intermittent +hydrogen bond lifetime, t0 = 0, tf = 2000 and dtmax = 30. In this way we create 30 windows timestep +(30 values in x axis). The continuos hydrogen bond lifetime should decay faster than intermittent. + +WOR +~~~ + +Analyzing water orientational relaxation :class:`WOR`. In this case we are analyzing "how fast" water molecules are rotating/changin direction. If WOR is very stable we can assume that water molecules are rotating/changing direction very slow, on the other hand, if WOR decay very fast, we can assume that water molecules are rotating/changing direction very fast:: + + import MDAnalysis + from MDAnalysis.analysis.waterdynamics import WOR + + u = MDAnalysis.Universe(pdb, trajectory) + selection = "byres name OH2 and sphzone 6.0 protein and resid 42" + WOR_analysis = WOR(universe, selection, 0, 1000, 20) + WOR_analysis.run() + i=0 + #now we print the data ready to graph. The first two columns are WOR_OH vs t graph, + #the second two columns are WOR_HH vs t graph and the third two columns are WOR_dip vs t graph + for WOR_OH, WOR_HH, WOR_dip in WOR_analysis.timeseries: + print i +" "+ WOR_OH +" "+ i +" "+ WOR_HH +" "+ i +" "+ WOR_dip + +where t0 = 0, tf = 1000 and dtmax = 20. In this way we create 20 windows timesteps (20 values in the x axis), +the first window is created with 1000 timestep average (1000/1), the second window is created with 500 +timestep average(1000/2), the third window is created with 333 timestep average (1000/3) and so on. + +AD +~~ + +Analyzing angular distribution :class:`AD` for OH vector, HH vector and dipole vector. It returns +a line histogram with vector orientation preference. A straight line in the output graphic means no preferential +orientation in water molecules. In this case we are analyzing if water molecules have some orientational +preference, in this way we can see if water molecules are under an electric field or if they are interacting +with something (residue, protein, etc):: + + import MDAnalysis + from MDAnalysis.analysis.waterdynamics import AD + + u = MDAnalysis.Universe(pdb, trajectory) + selection = "byres name OH2 and sphzone 6.0 (protein and (resid 42 or resid 26) )" + bins = 30 + AD_analysis = AD(universe,selection,bins) + AD_analysis.run() + #now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector , + #the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns + #are P(cos(theta)) vs cos(theta) for dipole vector + for i in range(bins): + print AD_analysis.graph[0][i] +" "+ AD_analysis.graph[1][i] +" "+ AD_analysis.graph[2][i] + +where P(cos(theta)) is the angular distribution or angular probabilities. + +MSD +~~~ +Analyzing mean square displacement :class:`MSD` for water molecules. In this case we are analyzing the average distance +that water molecules travels inside protein in XYZ direction (cylindric zone of radius 11[nm], Zmax 4.0[nm] and Zmin -8.0[nm]). A strong +rise mean a fast movement of water molecules, a weak rise mean slow movement of particles:: + + import MDAnalysis + from MDAnalysis.analysis.waterdynamics import MSD + + u = MDAnalysis.Universe(pdb, trajectory) + selection = "byres name OH2 and cyzone 11.0 4.0 -8.0 protein" + MSD_analysis = MSD(universe, selection, 0, 1000, 20) + MSD_analysis.run() + #now we print data ready to graph. The graph + #represents MSD vs t + i=0 + for msd in MSD_analysis.timeseries: + print i +" "+ msd + i += 1 + + +SP +~~ +Analyzing survival probability :class:`SP` for water molecules. In this case we are analyzing how long water +molecules remain in a sphere of radius 12.3 centered in the geometrical center of resid 42, 26, 34 and 80. +A slow decay of SP means a long permanence time of water molecules in the zone, on the +other hand, a fast decay means a short permanence time:: + + import MDAnalysis + from MDAnalysis.analysis.waterdynamics import SP + + u = MDAnalysis.Universe(pdb, trajectory) + selection = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26 or resid 34 or resid 80) " + SP_analysis = SP(universe, selection, 0, 100, 20) + SP_analysis.run() + #now we print data ready to graph. The graph + #represents SP vs t + i=0 + for sp in SP_analysis.timeseries: + print i +" "+ sp + i += 1 + +.. Output + +Output +------ + +HBL +~~~ +Hydrogen bond lifetime data is returned per window timestep, which is stored in +:attr:`HBL.timeseries` (in all the following descriptions, # indicates comments that are not part of the output):: + + results = [ + [ # time t0 + , + ], + [ # time t1 + , + ], + ... + ] + +WOR +~~~ +Water orientational relaxation data is returned per window timestep, which is stored in +:attr:`WOR.timeseries`:: + + results = [ + [ # time t0 + , , + ], + [ # time t1 + , , + ], + ... + ] + +AD +~~~ +Angular distribution data is returned per vector, which is stored in +:attr:`AD.graph`. In fact, AD returns a histogram:: + + results = [ + [ # OH vector values + # the values are order in this way: + , , ... + ], + [ # HH vector values + , , ... + ], + [ # dip vector values + , , ... + ], + ] + +MSD +~~~ +Mean Square Displacement data is returned in a list, which each element represents a MSD value in its respective +window timestep. Data is stored in :attr:`MSD.timeseries`:: + + results = [ + #MSD values orders by window timestep + , , ... + ] + +SP +~~ +Survival Probability data is returned in a list, which each element represents a MSD value in its respective +window timestep. Data is stored in :attr:`SP.timeseries`:: + + results = [ + # SP values order by window timestep + , , ... + ] + + + +Classes +-------- + +.. autoclass:: HBL + :members: + :inherited-members: + +.. autoclass:: WOR + :members: + :inherited-members: + +.. autoclass:: AD + :members: + :inherited-members: + +.. autoclass:: MSD + :members: + :inherited-members: + +.. autoclass:: SP + :members: + :inherited-members: + + +""" +from MDAnalysis import * +import MDAnalysis.analysis.hbonds +import MDAnalysis.core.Selection +import numpy +import argparse +import sys +import multiprocessing +import time +import itertools + +class HBL: + r""" + This is a autocorrelation function that gives the "Hydrogen Bond Lifetimes" proposed by D.C. Rapaport [Rapaport1983]_. From this + function we can obtain the continuos and intermittent behavior of hydrogen bonds in time. A + fast decay in these parameters indicate a fast change in HBs connectivity. A slow decay + indicate very stables hydrogen bonds, like in ice. The HBL is also know as "Hydrogen Bond Population + Relaxation" (HBPR). In the continuos case we have: + + .. math:: + C_{HB}^c(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h'_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)} + + where :math:`h'_{ij}(t_0+\tau)=1` if there is a H-bond between a pair :math:`ij` during time interval + :math:`t_0+\tau` (continuos) and :math:`h'_{ij}(t_0+\tau)=0` otherwise. In the intermittent case + we have: + + .. math:: + C_{HB}^i(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)} + + where :math:`h_{ij}(t_0+\tau)=1` if there is a H-bond between a pair :math:`ij` at time + :math:`t_0+\tau` (intermittent) and :math:`h_{ij}(t_0+\tau)=0` otherwise. + + :Arguments: + *universe* + Universe object + *selection1* + Selection string for first selection [‘byres name OH2’] + It could be any selection available in MDAnalysis, not just water. + *selection2* + Selection string to analize its HBL against selection1 + *t0* + Time where analysis begin + *tf* + Time where analysis end + *dtmax* + Maximum dt size, dtmax < tf or it will crash. + *nproc* + Number of processors to use, by default is 1. + + """ + def __init__(self,universe ,selection1 ,selection2, t0 , tf , dtmax, nproc = 1): + self.universe = universe + self.selection1 = selection1 + self.selection2 = selection2 + self.t0 = t0 + self.tf = tf - 1 + self.dtmax = dtmax + self.nproc = nproc + self.timeseries = None + + def _getC_i(self,HBP,t0,t): + """ + This function give the intermitent Hydrogen Bond Lifetime + C_i = / between t0 and t + """ + C_i = 0 + for i in range(len(HBP[t0])): + for j in range(len(HBP[t])): + if (HBP[t0][i][0] == HBP[t][j][0] and HBP[t0][i][1] == HBP[t][j][1]): + C_i += 1 + break + if len(HBP[t0]) == 0 : + return 0.0 + else: + return float(C_i)/len(HBP[t0]) + + def _getC_c(self,HBP,t0,t): + """ + This function give the continous Hydrogen Bond Lifetime + C_c = / between t0 and t + """ + C_c = 0 + dt = 1 + begt0 = t0 + HBP_cp = HBP + HBP_t0 = HBP[t0] + newHBP = [] + if t0==t: + return 1.0 + while t0+dt <= t: + for i in range(len(HBP_t0)): + for j in range(len(HBP_cp[t0+dt])): + if (HBP_t0[i][0] == HBP_cp[t0+dt][j][0] and HBP_t0[i][1] == HBP_cp[t0+dt][j][1]): + newHBP.append(HBP_t0[i]) + break + C_c = len(newHBP) + t0 += dt + HBP_t0 = newHBP + newHBP = [] + if len(HBP[begt0]) == 0 : + return 0 + else: + return C_c/float(len(HBP[begt0])) + + def _intervC_c(self,HBP,t0,tf,dt): + """ + This function gets all the data for the h(t0)h(t0+dt)', where + t0 = 1,2,3,...,tf. This function give us one point of the final graphic + HBL vs t + """ + a = 0 + count = 0 + for i in range(len(HBP)): + if (t0+dt <= tf): + if t0 == t0+dt: + b = self._getC_c(HBP,t0,t0) + break + b = self._getC_c(HBP,t0,t0+dt) + t0 += dt + a += b + count += 1 + if count == 0: + return 1.0 + return a/count + + def _intervC_i(self,HBP,t0,tf,dt): + """ + This function gets all the data for the h(t0)h(t0+dt), where + t0 = 1,2,3,...,tf. This function give us a point of the final graphic + HBL vs t + """ + a = 0 + count = 0 + for i in range(len(HBP)): + if (t0+dt <= tf ): + b = self._getC_i(HBP,t0,t0+dt) + t0 += dt + a += b + count += 1 + return a/count + + def _finalGraphGetC_i(self,HBP,t0,tf,maxdt): + """ + This function gets the final data of the C_i graph. + """ + output = [] + for dt in range(maxdt): + a = self._intervC_i(HBP,t0,tf,dt) + output.append(a) + return output + + def _finalGraphGetC_c(self,HBP,t0,tf,maxdt): + """ + This function gets the final data of the C_c graph. + """ + output = [] + for dt in range(maxdt): + a = self._intervC_c(HBP,t0,tf,dt) + output.append(a) + return output + + def _getGraphics(self,HBP,t0,tf,maxdt): + """ + Function that join all the results into a graphics. + """ + a = [] + cont = self._finalGraphGetC_c(HBP,t0,tf,maxdt) + inte = self._finalGraphGetC_i(HBP,t0,tf,maxdt) + for i in range(len(cont)): + fix = [cont[i],inte[i]] + a.append(fix) + return a + + def _HBA(self,ts , conn , universe , selAtom1 , selAtom2 ): + """ + Main function for calculate C_i and C_c in parallel. + """ + finalGetResidue1 = selAtom1 + finalGetResidue2 = selAtom2 + frame = ts.frame + h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(universe, finalGetResidue1 , finalGetResidue2 , distance=3.5, angle=120.0, start=frame-1 , stop=frame) + while True: + try: + h.run() + break + except: + print "error" + print "trying again" + sys.stdout.flush() + sys.stdout.flush() + conn.send( h.timeseries[0] ) + conn.close() + + + def run(self,**kwargs): + """ + Analyze trajectory and produce timeseries + """ + h_list = [] + i = 0 + if (self.nproc > 1): + while i < len(self.universe.trajectory): + jobs = [] + k=i + for j in range(self.nproc): + #start + print "ts=",i+1 + if i >= len(self.universe.trajectory): + break + conn_parent, conn_child = multiprocessing.Pipe(False) + while True: + try: + #new thread + jobs.append(( multiprocessing.Process(target = self._HBA ,args = (self.universe.trajectory[i],conn_child,self.universe, self.selection1 , self.selection2,) ) , conn_parent) ) + break + except: + print "error in jobs.append" + jobs[j][0].start() + i = i + 1 + + for j in range(self.nproc): + if k >= len(self.universe.trajectory): + break + rec01 = jobs[j][1] + received = rec01.recv() + h_list.append(received) + jobs[j][0].join() + k += 1 + self.timeseries = self._getGraphics( h_list , 0 , self.tf-1 , self.dtmax ) + else: + h_list = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(self.universe, self.selection1 , self.selection2,distance=3.5, angle=120.0) + h_list.run() + self.timeseries = self._getGraphics( h_list.timeseries , self.t0 , self.tf , self.dtmax ) + +class WOR: + r""" + Function to evaluate the Water Orientational Relaxation proposed by Yu-ling Yeh + and Chung-Yuan Mou [Yeh1999_]. WOR indicates "how fast" water molecules are rotating + or changing direction. This is a time correlation function given by: + + .. math:: + C_{\hat u}(\tau)=\langle \mathit{P}_2[\mathbf{\hat{u}}(t_0)\cdot\mathbf{\hat{u}}(t_0+\tau)]\rangle + + where :math:`P_2=(3x^2-1)/2` is the second-order Legendre polynomial and :math:`\hat{u}` is + a unit vector along HH, OH or dipole vector. + + + :Arguments: + *universe* + Universe object + *selection* + Selection string, only models with 3 atoms molecules are allowed (TIP3, TIP3P, etc) + *t0* + Time where analysis begin + *tf* + Time where analysis end + *dtmax* + Maximum dt size window, dtmax < tf or it will crash. + + """ + + def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): + self.universe = universe + self.selection = selection + self.t0 = t0 + self.tf = tf + self.dtmax= dtmax + self.nproc = nproc + self.timeseries = None + + def _repitedIndex(self,selection,dt,totalFrames): + """ + Indicate the comparation between all the t+dt. + The results is a list of list with all the repited index per frame (or time). + Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... + Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... + Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... + """ + rep=[] + for i in range(int(round( (totalFrames-1)/float(dt) ) ) ): + if ( dt*i+dt < totalFrames ): + rep.append(self._sameMolecTandDT(selection,dt*i,(dt*i)+dt)) + return rep + + def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): + """ + Give one point to promediate and get one point of the graphic C_vect vs t + Ex: t0=1 and tau=1 so calculate the t0-tau=1-2 intervale. + Ex: t0=5 and tau=3 so calcultate the t0-tau=5-8 intervale. + i = come from getMeanOneCOHPoint (named j) (int) + """ + valOH = 0 + valHH = 0 + valdip= 0 + n = 0 + for j in range(len(repInd[i])/3): + begj = 3*j + universe.trajectory[t0] + Ot0 = repInd[i][begj] + H1t0 = repInd[i][begj+1] + H2t0 = repInd[i][begj+2] + OHVector0 = H1t0.position - Ot0.position + HHVector0 = H1t0.position-H2t0.position + dipVector0 = ((H1t0.position + H2t0.position)*0.5)-Ot0.position + + universe.trajectory[t0+dt] + Otp = repInd[i][begj] + H1tp = repInd[i][begj+1] + H2tp = repInd[i][begj+2] + + OHVectorp = H1tp.position- Otp.position + HHVectorp = H1tp.position - H2tp.position + dipVectorp = ((H1tp.position + H2tp.position)*0.5)-Otp.position + + normOHVector0 = numpy.linalg.norm(OHVector0) + normOHVectorp = numpy.linalg.norm(OHVectorp) + normHHVector0 = numpy.linalg.norm(HHVector0) + normHHVectorp = numpy.linalg.norm(HHVectorp) + normdipVector0 = numpy.linalg.norm(dipVector0) + normdipVectorp = numpy.linalg.norm(dipVectorp) + + unitOHVector0 = [OHVector0[0]/normOHVector0,OHVector0[1]/normOHVector0,OHVector0[2]/normOHVector0] + unitOHVectorp = [OHVectorp[0]/normOHVectorp,OHVectorp[1]/normOHVectorp,OHVectorp[2]/normOHVectorp] + unitHHVector0 = [HHVector0[0]/normHHVector0,HHVector0[1]/normHHVector0,HHVector0[2]/normHHVector0] + unitHHVectorp = [HHVectorp[0]/normHHVectorp,HHVectorp[1]/normHHVectorp,HHVectorp[2]/normHHVectorp] + unitdipVector0 = [dipVector0[0]/normdipVector0,dipVector0[1]/normdipVector0,dipVector0[2]/normdipVector0] + unitdipVectorp = [dipVectorp[0]/normdipVectorp,dipVectorp[1]/normdipVectorp,dipVectorp[2]/normdipVectorp] + + valOH += self.lg2(numpy.dot(unitOHVector0,unitOHVectorp)) + valHH += self.lg2(numpy.dot(unitHHVector0,unitHHVectorp)) + valdip += self.lg2(numpy.dot(unitdipVector0,unitdipVectorp)) + n += 1 + valOH = valOH/n + valHH = valHH/n + valdip = valdip/n + return (valOH,valHH,valdip) + + def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): + """ + This function get one point of the graphic C_OH vs t. It uses the + _getOneDeltaCOHPoint() function to calculate the average. + + """ + repInd = self._repitedIndex(selection1,dt,totalFrames) + sumsdt = 0 + n = 0.0 + sumDeltaOH = 0.0 + sumDeltaHH = 0.0 + sumDeltadip = 0.0 + valOHList = [] + valHHList = [] + valdipList = [] + + for j in range(totalFrames/dt-1): + a = self._getOneDeltaCOHPoint(universe,repInd,j,sumsdt,dt) + sumDeltaOH += a[0] + sumDeltaHH += a[1] + sumDeltadip += a[2] + valOHList.append(a[0]) + valHHList.append(a[1]) + valdipList.append(a[2]) + sumsdt += dt + n += 1 + return (sumDeltaOH/n,sumDeltaHH/n,sumDeltadip/n) + + def _sameMolecTandDT(self,selection,t0d,tf): + """ + Compare the molecules in the t0d selection and the t0d+dt selection and + select only the particles that are repited in both frame. This is to consider + only the molecules that remains in the selection after the dt time has elapsed. + The result is a list with the indexs of the atoms. + """ + a = set(selection[t0d]) + b = set(selection[tf]) + sort = sorted(list(a.intersection(b))) + return sort + + def _selection_serial(self,universe,selection_str): + selection = [] + for ts in universe.trajectory: + selection.append(universe.selectAtoms(selection_str)) + print ts.frame + return selection + + # Second Legendre polynomial + lg2 = lambda self,x : (3*x*x - 1)/2 + + def run(self,**kwargs): + """ + Analyze trajectory and produce timeseries + """ + + #All the selection to an array, this way is faster than selecting later. + if self.nproc==1: + selection_out = self._selection_serial(self.universe,self.selection) + else: + #selection_out = self._selection_parallel(self.universe,self.selection,self.nproc) + #parallel selection to be implemented + selection_out = self._selection_serial(self.universe,self.selection) + self.timeseries = [] + for dt in list(range(1,self.dtmax+1)): + output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + self.timeseries.append(output) + + + +class AD: + r""" + The angular distribution function (AD) is defined as the distribution + probability of the cosine of the :math:`\theta` angle formed by the OH vector, HH vector + or dipolar vector of water molecules and a vector :math:`\hat n` parallel to chosen axis + (z is the default value). The cosine is define as :math:`\cos \theta = \hat u \cdot \hat n`, where :math:`\hat u` is OH, HH or dipole vector. + It creates a histogram and returns a list of lists, see Output_. The AD is also know as Angular Probability (AP). + + :Arguments: + *universe* + Universe object + *selection* + Selection string to evaluate its angular distribution [‘byres name OH2’] + *bins* + Number of bins to create the histogram by means of numpy.histogram_ [40] + *axis* + Axis to create angle with the vector (HH, OH or dipole) and calculate cosine theta ['z']. Options: 'x', + 'y', 'z' + + .. _numpy.histogram: http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html + + + """ + def __init__(self,universe,selection_str,bins=40,nproc=1,axis="z"): + self.universe = universe + self.selection_str = selection_str + self.bins = bins + self.nproc = nproc + self.axis = axis + self.graph = None + + def _getCosTheta(self,universe,selection,axis): + valOH = [] + valHH = [] + valdip= [] + i = 0 + while i <= (len(selection)-1): + universe.trajectory[i] + line = selection[ i ].positions + j=0 + while j < len(line): + Ot0 = numpy.array(line[j]) + H1t0 = numpy.array(line[j+1]) + H2t0 = numpy.array(line[j+2]) + + OHVector0 = H1t0 - Ot0 + HHVector0 = H1t0 - H2t0 + dipVector0 = (H1t0 + H2t0)*0.5 - Ot0 + + normOHVector0 = numpy.linalg.norm(OHVector0) + normHHVector0 = numpy.linalg.norm(HHVector0) + normdipVector0 = numpy.linalg.norm(dipVector0) + + unitOHVector0 = [OHVector0[0]/normOHVector0,OHVector0[1]/normOHVector0,OHVector0[2]/normOHVector0] + unitHHVector0 = [HHVector0[0]/normHHVector0,HHVector0[1]/normHHVector0,HHVector0[2]/normHHVector0] + unitdipVector0 = [dipVector0[0]/normdipVector0,dipVector0[1]/normdipVector0,dipVector0[2]/normdipVector0] + + if axis == "z": + valOH.append(numpy.dot(unitOHVector0,[0,0,1])) + valHH.append(numpy.dot(unitHHVector0,[0,0,1])) + valdip.append(numpy.dot(unitdipVector0,[0,0,1])) + + elif axis == "x": + valOH.append(numpy.dot(unitOHVector0,[1,0,0])) + valHH.append(numpy.dot(unitHHVector0,[1,0,0])) + valdip.append(numpy.dot(unitdipVector0,[1,0,0])) + + elif axis == "y": + valOH.append(numpy.dot(unitOHVector0,[0,1,0])) + valHH.append(numpy.dot(unitHHVector0,[0,1,0])) + valdip.append(numpy.dot(unitdipVector0,[0,1,0])) + + + j += 3 + i += 1 + return (valOH,valHH,valdip) + + def _getHistogram(self,universe,selection,bins,axis): + """ + This function gets a normalized histogram of the cos(theta) values. It return a list of list. + """ + a = self._getCosTheta(universe,selection,axis) + cosThetaOH = a[0] + cosThetaHH = a[1] + cosThetadip = a[2] + lencosThetaOH = len(cosThetaOH) + lencosThetaHH = len(cosThetaHH) + lencosThetadip = len(cosThetadip) + histInterval = bins + histcosThetaOH = numpy.histogram(cosThetaOH,histInterval, normed = True) + histcosThetaHH = numpy.histogram(cosThetaHH,histInterval, normed = True) + histcosThetadip = numpy.histogram(cosThetadip,histInterval, normed = True) + + return (histcosThetaOH,histcosThetaHH,histcosThetadip) + + def _hist2column(self,aList): + """ + This function transform from the histogram format + to a column format. + """ + a = [] + for x in itertools.izip_longest(*aList, fillvalue="."): + a.append(" ".join(str(i) for i in x)) + return a + + def run(self,**kwargs): + """ + Function to evaluate the angular distribution of cos(theta) + """ + + if self.nproc ==1: + selection = self._selection_serial(self.universe,self.selection_str) + else: + #not implemented yet + #selection = self._selection_parallel(self.universe,self.selection_str,self.nproc) + selection = self._selection_serial(self.universe,self.selection_str) + + self.graph = [] + output=self._getHistogram(self.universe,selection,self.bins,self.axis) + #this is to format the exit of the file + #maybe this output could be improved + listOH = [list(output[0][1]),list(output[0][0])] + listHH = [list(output[1][1]),list(output[1][0])] + listdip = [list(output[2][1]),list(output[2][0])] + + self.graph.append(self._hist2column(listOH)) + self.graph.append(self._hist2column(listHH)) + self.graph.append(self._hist2column(listdip)) + + def _selection_serial(self,universe,selection_str): + selection = [] + for ts in universe.trajectory: + selection.append(universe.selectAtoms(selection_str)) + print ts.frame + return selection + + +class MSD: + r""" + Function to evaluate the Mean Square Displacement (MSD_). The MSD gives the average distance that + particles travels. The MSD is given by: + + .. math:: + \langle\Delta r(t)^2\rangle = 2nDt + + where :math:`r(t)` is the position of particle in time :math:`t`, :math:`\Delta r(t)` is the displacement + after time lag :math:`t`, :math:`n` is the dimensionality, in this case :math:`n=3`, :math:`D` is the diffusion + coefficient and :math:`t` is the time. + + .. _MSD: http://en.wikipedia.org/wiki/Mean_squared_displacement + + :Arguments: + *universe* + Universe object + *selection* + Selection string + *t0* + Time where analysis begin + *tf* + Time where analysis end + *dtmax* + Maximum dt size window, dtmax < tf or it will crash. + + """ + + def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): + self.universe = universe + self.selection = selection + self.t0 = t0 + self.tf = tf + self.dtmax= dtmax + self.nproc = nproc + self.timeseries = None + + def _repitedIndex(self,selection,dt,totalFrames): + """ + Indicate the comparation between all the t+dt. + The results is a list of list with all the repited index per frame (or time). + Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... + Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... + Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... + """ + rep=[] + for i in range(int(round( (totalFrames-1)/float(dt) ) ) ): + if ( dt*i+dt < totalFrames ): + rep.append(self._sameMolecTandDT(selection,dt*i,(dt*i)+dt)) + return rep + + def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): + """ + Give one point to promediate and get one point of the grapic C_vect vs t + Ex: t0=1 and dt=1 so calculate the t0-dt=1-2 intervale. + Ex: t0=5 and dt=3 so calcultate the t0-dt=5-8 intervale + i = come from getMeanOneCOHPoint (named j) (int) + """ + valO = 0 + n = 0 + for j in range(len(repInd[i])/3): + begj = 3*j + universe.trajectory[t0] + #Plus zero is to avoid 0to be equal to 0tp + Ot0 = repInd[i][begj].position + 0 + + universe.trajectory[t0+dt] + #Plus zero is to avoid 0to be equal to 0tp + Otp = repInd[i][begj].position + 0 + + #position oxygen + OVector = Ot0 - Otp + #here it is the difference with waterdynamics.WOR + valO += numpy.dot(OVector, OVector) + n += 1 + valO = valO/n + return (valO) + + def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): + """ + This function get one point of the graphic C_OH vs t. It's uses the + _getOneDeltaCOHPoint() function to calculate the average. + + """ + repInd = self._repitedIndex(selection1,dt,totalFrames) + sumsdt = 0 + n = 0.0 + sumDeltaO = 0.0 + valOList = [] + + for j in range(totalFrames/dt-1): + a = self._getOneDeltaCOHPoint(universe,repInd,j,sumsdt,dt) + print "a=",a + sumDeltaO += a + valOList.append(a) + sumsdt += dt + n += 1 + return (sumDeltaO/n) + + def _sameMolecTandDT(self,selection,t0d,tf): + """ + Compare the molecules in the t0d selection and the t0d+dt selection and + select only the particles that are repited in both frame. This is to consider + only the molecules that remains in the selection after the dt time has elapsed. + The result is a list with the indexs of the atoms. + """ + a = set(selection[t0d]) + b = set(selection[tf]) + sort = sorted(list(a.intersection(b))) + return sort + + def _selection_serial(self,universe,selection_str): + selection = [] + for ts in universe.trajectory: + selection.append(universe.selectAtoms(selection_str)) + print ts.frame + return selection + + def run(self,**kwargs): + """ + Analyze trajectory and produce timeseries + """ + + #All the selection to an array, this way is faster than selecting later. + if self.nproc==1: + selection_out = self._selection_serial(self.universe,self.selection) + else: + #parallel not yet implemented + #selection = selection_parallel(universe,selection_str,nproc) + selection_out = self._selection_serial(self.universe,self.selection) + self.timeseries = [] + for dt in list(range(1,self.dtmax+1)): + output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + self.timeseries.append(output) + + +class SP: + r""" + Function to evaluate the Survival Probability (SP). The SP gives the probability + for a group of particles to remain in certain region. The SP is given by: + SP = equation + + .. math:: + P(\tau) = \frac1T \sum_{t=1}^T \frac{N(t,t+\tau)}{N(t)} + + :Arguments: + *universe* + Universe object + *selection* + Selection string, any selection is allowed, with this selection you define the region/zone where + to analize, i.e.: "selection_a" and "zone" (see SP examples_ ) + *t0* + Time where analysis begin + *tf* + Time where analysis end + *dtmax* + Maximum dt size window, dtmax < tf or it will crash + + """ + + def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): + self.universe = universe + self.selection = selection + self.t0 = t0 + self.tf = tf + self.dtmax= dtmax + self.nproc = nproc + self.timeseries = None + + def _getOneDeltaCOHPoint(self,selection, totalFrames, t0, tau): + """ + Give one point to promediate and get one point of the graphic C_vect vs t + Ex: t0=1 and tau=1 so calculate the t0-tau=1-2 intervale. + Ex: t0=5 and tau=3 so calcultate the t0-tau=5-8 intervale. + """ + Ntau = self._NumPart_tau(selection, totalFrames, t0, tau) + Nt = float(self._NumPart(selection,t0)) + + return Ntau/Nt + + def _getMeanOneCOHPoint(self,universe,selection1,selection_str,wint,totalFrames): + """ + This function get one point of the graphic P(t) vs t. It uses the + _getOneDeltaCOHPoint() function to calculate the average. + + """ + n = 0.0 + sumDeltaP = 0.0 + valPList = [] + for frame in range(totalFrames-wint): + a = self._getOneDeltaCOHPoint(selection1,totalFrames ,frame, wint) + sumDeltaP += a + n += 1 + + return sumDeltaP/n + + def _NumPart_tau(self,selection, totalFrames, t0,tau): + """ + Compare the molecules in t0 selection and t0+tau selection and + select only the particles that remaing from t0 to t0+tau. It returns + the number of remaining particles. + """ + a = set(selection[t0]) + i=0 + while (t0+i) < t0+tau and (t0+i) < totalFrames: + b = set(selection[t0+i]) + a = a.intersection(b) + i += 1 + return len(a) + + def _NumPart(self, selection, t): + return len(selection[t]) + + def _selection_serial(self,universe,selection_str): + selection = [] + for ts in universe.trajectory: + selection.append(universe.selectAtoms(selection_str)) + print ts.frame + return selection + + def run(self,**kwargs): + """ + Analyze trajectory and produce timeseries + """ + + #All the selection to an array, this way is faster than selecting later. + if self.nproc==1: + selection_out = self._selection_serial(self.universe,self.selection) + else: + #selection = selection_parallel(universe,selection_str,nproc) + #parallel selection to be implemented + selection_out = self._selection_serial(self.universe,self.selection) + self.timeseries = [] + for dt in list(range(1,self.dtmax+1)): + output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + self.timeseries.append(output) diff --git a/package/doc/sphinx/source/documentation_pages/analysis/waterdynamics.rst b/package/doc/sphinx/source/documentation_pages/analysis/waterdynamics.rst new file mode 100644 index 00000000000..2db6600b282 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/waterdynamics.rst @@ -0,0 +1,2 @@ +.. automodule:: MDAnalysis.analysis.waterdynamics + diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index 1dbe857305f..2b9a555672e 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -88,5 +88,6 @@ Volumetric analysis :maxdepth: 1 analysis/density + analysis/waterdynamics diff --git a/testsuite/MDAnalysisTests/data/watdyn.dcd b/testsuite/MDAnalysisTests/data/watdyn.dcd new file mode 100644 index 0000000000000000000000000000000000000000..8eec31cd81c38aaf3eaa1089dddf3664e1800bca GIT binary patch literal 2876 zcmbuB3s98T6~{r11nUckh-l5KqeevpvDRSN12j8^Ph7cbIyF9o{g8Jrz# z{KKrwrR$|4I~NDIIoOL9triLE_M%^D@ABCT!PX*+`Ex8kpJO@izq9?p)xpJD0Z zHNnftax@toCG{9TS|`xs!O zMngpv_LUP&im9PNJ@s^6q@vOODw_U+)2*p0Iugw&HQ*Ah-=w6Vz6;cQj?vp317W~2 zW6=?tfKHtRc;pHUt|Xy&X$l;pQ}B4M99OHeFmrbiyte0|B)S++`ir4=@GK6_E!Hwf zPj6rgEeyD{rG>e+?i9~GG)4l;c;y72Aa1@kmG5m2=VT}^mg9WEAw)_tVG&sf54#h{ z6P03TZ7D)8RwAYF{W<7vxyB~?ney3Fp0GPbYi;h@c}m9YZBpuo`SF--F??wH7QUcq zE9+k|i5)9!RpzhSz#AkRxnArnErYcuIlMp8z>BJ}xLidK+871Rt*0&vHHEL^q<4yu z{%|#Y?5L(rp%%8rE2y%F)3?J1&Y(3p9?NVJ@y1667^Go$L=uFP3z@sn6RMMN`MXS1 z{Js$OgHNEzzZmXcl_8<70%47%+H-KO@EY6JY{UajJYd%*@3i?WVX@@(x^CsvPAUJz zH-dZLcIRtWhO)qqomg;Ci^6@!dhSpW&Rt(k(=u3D%3DwX|fmn)JAudbSgtZmp*SjtUZes-iCMfinpGIvykXVsUIvJnCMi3bmMsjC*o4 zJ(gp`ffVeYkOlXi0({eygB=fwu`sv{X&q+~=31s@@G|-uODY@BeI0*hE3cc0W4Wc| z&@NqeZEFxe)F$O^X&ZUP$ZafP`C>NCD__A!`tWDVqxd%|Q?(3?3mI|~)p)g}3c+tw zG^|tj*Zv0T?NHJDRW)Q5rJ)O`BjXoDUB@+)VWpt&9aUub$AS061%~l39~F=7vI7Xw zq(bLdBK(czxI!st?@z<`bF#4Xasm44a}W|+3g1O#=-pp|UC-aw!VjK*W+&TCxqF1b zU}%Q;%G0Hig@ev3JEKDR(e5=|x^^GmGh!S2U4Q{Q+qOa}&G6>IFSl?qoUUcSqLru` zsYXa%73OuS=;sHVMm}qxZn=tF&(~1ReGUDXRY%#s(NKw(h7Nlv$a{;54!<5agSPBA z7Pq_5RFz)h+LVUw?lO(vuSYbaioOhmG!*9eT^ZR2rGwV#1a(w17Wx%F2d{V+{ zenmc6%RqcZfzj0*y?blm<)Egq`iw$DF45>n72R1@OOKlCDCuu?r230c6S^8o6Z+xN zP);(%z-wWaV;nN(Cc`CSKYTBwBD_Qf@4zIa)~Dg~5m^YZ%|v#10XACXz(!tz`NF+0 zC!i8n=f7uglkTwedAhuH&p1B)`hBsTpNFJPlU(Db7s0PT^5j3}WOCc~wM@^zkTu-Q zR!n`cns-Q-abwpRS_a;M3V0@uA{ry}PJ1x4&q)9c%H)L2_b zE<@|-k*SiZgr3m+lm>h+?DdI*_3>oX+QtAisR(^2!m2c#n?Nu1ecU&FmCeuH`0)ct*r66E}!5$kq;_+EcT!0BQct3Dz-L_0DQy+%-wd;7GFqv3a$Kyh_OJt|PxGoX)8l>6oFK^3pj~ foKn;Ho;tD$Q&WaP1MNGiq(1+PbZ?sKUkv^U5*Gv% literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/data/watdyn.psf b/testsuite/MDAnalysisTests/data/watdyn.psf new file mode 100644 index 00000000000..75301b17a68 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/watdyn.psf @@ -0,0 +1,54 @@ +PSF + + 3 !NTITLE + REMARKS original generated structure x-plor psf file + REMARKS topology toppar_water_ions.top + REMARKS segment WAT { first NONE; last NONE; auto none } + + 15 !NATOM + 1 WAT 5 TIP3 OH2 OT -0.834000 15.9994 0 + 2 WAT 5 TIP3 H1 HT 0.417000 1.0080 0 + 3 WAT 5 TIP3 H2 HT 0.417000 1.0080 0 + 4 WAT 7 TIP3 OH2 OT -0.834000 15.9994 0 + 5 WAT 7 TIP3 H1 HT 0.417000 1.0080 0 + 6 WAT 7 TIP3 H2 HT 0.417000 1.0080 0 + 7 WAT 8 TIP3 OH2 OT -0.834000 15.9994 0 + 8 WAT 8 TIP3 H1 HT 0.417000 1.0080 0 + 9 WAT 8 TIP3 H2 HT 0.417000 1.0080 0 + 10 WAT 15 TIP3 OH2 OT -0.834000 15.9994 0 + 11 WAT 15 TIP3 H1 HT 0.417000 1.0080 0 + 12 WAT 15 TIP3 H2 HT 0.417000 1.0080 0 + 13 WAT 21 TIP3 OH2 OT -0.834000 15.9994 0 + 14 WAT 21 TIP3 H1 HT 0.417000 1.0080 0 + 15 WAT 21 TIP3 H2 HT 0.417000 1.0080 0 + + 15 !NBOND: bonds + 1 2 1 3 2 3 4 5 + 4 6 5 6 7 8 7 9 + 8 9 10 11 10 12 11 12 + 13 14 13 15 14 15 + + 5 !NTHETA: angles + 2 1 3 5 4 6 8 7 9 + 11 10 12 14 13 15 + + 0 !NPHI: dihedrals + + + 0 !NIMPHI: impropers + + + 0 !NDON: donors + + + 0 !NACC: acceptors + + + 0 !NNB + + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 + + 1 0 !NGRP + 0 0 0 + diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 1905b4f9448..9b37544044d 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -78,6 +78,7 @@ "GMS_SYMOPT", # GAMESS D4h optimization "GMS_ASYMSURF", # GAMESS C1 surface "two_water_gro", "two_water_gro_nonames", # for bond guessing, 2 water molecules, one with weird names + "waterPSF","waterDCD", ] from pkg_resources import resource_filename @@ -208,3 +209,6 @@ two_water_gro = resource_filename(__name__, "data/two_water_gro.gro") two_water_gro_nonames = resource_filename(__name__, "data/two_water_gro_nonames.gro") + +waterPSF = resource_filename(__name__, 'data/watdyn.psf') +waterDCD = resource_filename(__name__, 'data/watdyn.dcd') diff --git a/testsuite/MDAnalysisTests/test_analysis.py b/testsuite/MDAnalysisTests/test_analysis.py index e2456bb27a2..dd983a99a29 100644 --- a/testsuite/MDAnalysisTests/test_analysis.py +++ b/testsuite/MDAnalysisTests/test_analysis.py @@ -19,6 +19,7 @@ import MDAnalysis.analysis.align import MDAnalysis.analysis.hbonds import MDAnalysis.analysis.helanal +import MDAnalysis.analysis.waterdynamics from MDAnalysis import SelectionError, FinishTimeException from numpy.testing import * @@ -28,7 +29,7 @@ import errno import tempfile -from MDAnalysis.tests.datafiles import PSF, DCD, FASTA, PDB_helix, PDB_HOLE, XTC_HOLE, GRO, XTC +from MDAnalysis.tests.datafiles import PSF, DCD, FASTA, PDB_helix, PDB_HOLE, XTC_HOLE, GRO, XTC, waterDCD, waterPSF from . import executable_not_found_runtime @@ -359,3 +360,38 @@ def test_xtc_striding(self): MDAnalysis.analysis.helanal.helanal_trajectory(u, selection=sel, finish=5) except IndexError: self.fail("IndexError consistent with Issue 188.") + + +class TestWaterdynamics(TestCase): + def setUp(self): + self.universe = MDAnalysis.Universe(waterPSF, waterDCD) + self.selection1 = "byres name OH2" + self.selection2 = self.selection1 + + def test_HBL(self): + hbl = MDAnalysis.analysis.waterdynamics.HBL(self.universe, self.selection1, self.selection2, 0, 5, 3) + hbl.run() + assert_equal(round(hbl.timeseries[2][1],5), 0.75) + + def test_WOR(self): + wor = MDAnalysis.analysis.waterdynamics.WOR(self.universe, self.selection1, 0, 5, 2) + wor.run() + assert_equal(round(wor.timeseries[1][2],5), 0.35887) + + def test_AD(self): + ad = MDAnalysis.analysis.waterdynamics.AD(self.universe,self.selection1,40) + ad.run() + #debo arreglar la salida del AD + #must fix AD output + assert_equal(str(ad.graph[0][39]), str("0.951172947884 0.48313682125") ) + + def test_MSD(self): + msd = MDAnalysis.analysis.waterdynamics.MSD(self.universe, self.selection1, 0, 10, 2) + msd.run() + assert_equal(round(msd.timeseries[1],5), 0.03984) + + def test_SP(self): + sp = MDAnalysis.analysis.waterdynamics.SP(self.universe, self.selection1, 0, 6, 3) + sp.run() + assert_equal(round(sp.timeseries[1],5), 1.0) + From 79038a727bdcc2f6d8d72cddb7f270335a88208c Mon Sep 17 00:00:00 2001 From: Alejandro Bernardin Date: Tue, 9 Jun 2015 15:15:21 -0300 Subject: [PATCH 2/4] Added corrected suggestions in #295. Documentation fixed and updated. --- package/MDAnalysis/analysis/waterdynamics.py | 284 +++++++++---------- testsuite/MDAnalysisTests/test_analysis.py | 22 +- 2 files changed, 149 insertions(+), 157 deletions(-) diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index ef0d021688b..4d0a7d183d7 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -1,18 +1,17 @@ -# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; encoding: utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- http://www.MDAnalysis.org +# Copyright (c) 2006-2015 Naveen Michaud-Agrawal, Elizabeth J. Denning, Oliver Beckstein +# and contributors (see AUTHORS for the full list) # -# MDAnalysis --- http://mdanalysis.googlecode.com -# Copyright (c) 2006-2012 Naveen Michaud-Agrawal, -# Elizabeth J. Denning, Oliver Beckstein, -# and contributors (see website for details) # Released under the GNU Public Licence, v2 or any higher version # # Please cite your use of MDAnalysis in published work: # -# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and -# O. Beckstein. MDAnalysis: A Toolkit for the Analysis of -# Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319--2327, -# doi:10.1002/jcc.21787 +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # """ @@ -21,7 +20,7 @@ :Author: Alejandro Bernardin :Year: 2014-2015 -:Copyright: GNU Public License v2 +:Copyright: GNU Public License v3 This module provides functions to analize water dynamics trajectories and water interactions with other molecules. The functions in this module are: water orientational relaxation (WOR) [Yeh1999]_, hydrogen bond lifetimes (HBL) [Rapaport1983]_, @@ -48,9 +47,9 @@ .. [Brodka1994] Aleksander Brodka (1994). Diffusion in restricted volume, Molecular Physics, 1994, Vol. 82, No. 5, 1075-1078. -.. [Araya-Secchi2014] Raul Araya-Secchi, Seung-gu Kang, Tien Huynh, Alejandro Bernardin, - Agustin D. Martinez, Juan C. Saez, Tomas Perez-Acle, Ruhong Zhou. Characterization of a - novel water pocket inside the human Cx26 hemi-channel structure. +.. [Araya-Secchi2014] Araya-Secchi, R., Tomas Perez-Acle, Seung-gu Kang, Tien Huynh, Alejandro Bernardin, Yerko Escalona, Jose-Antonio Garate, Agustin D. Martinez, + Isaac E. Garcia, Juan C. Saez, Ruhong Zhou (2014). Characterization of a novel water pocket inside the human Cx26 hemichannel structure. Biophysical journal, 107(3), 599-612. + .. [Milischuk2011] Anatoli A. Milischuk and Branka M. Ladanyi. Structure and dynamics of water confined in silica nanopores. J. Chem. Phys. 135, 174709 (2011); doi: 10.1063/1.3657408 @@ -61,16 +60,16 @@ Examples -------- -HBL -~~~~ +Hydrogen bond lifetimes +~~~~~~~~~~~~~~~~~~~~~~~ -Analyzing hydrogen bond lifetimes :class:`HBL`, both continuos and intermittent. In this case we are analyzing +Analyzing hydrogen bond lifetimes (HBL) :class:`HydrogenBondLifetimes`, both continuos and intermittent. In this case we are analyzing how residue 38 interact with a water sphere of radius 6.0 centered on the geometric center of protein and -residue 42. If the HBL are very stable, we can assume that residue 38 is hydrophilic, on the other -hand, if the HBL are very unstable, we can assume that residue 38 is hydrophobic:: +residue 42. If the hydrogen bond lifetimes are very stable, we can assume that residue 38 is hydrophilic, on the other +hand, if the are very unstable, we can assume that residue 38 is hydrophobic:: import MDAnalysis - from MDAnalysis.analysis.waterdynamics import HBL + from MDAnalysis.analysis.waterdynamics import HydrogenBondLifetimes as HBL u = MDAnalysis.Universe(pdb, trajectory) selection1 = "byres name OH2 and sphzone 6.0 protein and resid 42" @@ -84,17 +83,17 @@ print i +" "+ HBLc +" "+ i +" "+ HBLi i+=1 -where HBLc is the value for the continuos hydrogen bond lifetime and HBLi is the value for the intermittent +where HBLc is the value for the continuos hydrogen bond lifetimes and HBLi is the value for the intermittent hydrogen bond lifetime, t0 = 0, tf = 2000 and dtmax = 30. In this way we create 30 windows timestep -(30 values in x axis). The continuos hydrogen bond lifetime should decay faster than intermittent. +(30 values in x axis). The continuos hydrogen bond lifetimes should decay faster than intermittent. -WOR -~~~ +WaterOrientationalRelaxation +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Analyzing water orientational relaxation :class:`WOR`. In this case we are analyzing "how fast" water molecules are rotating/changin direction. If WOR is very stable we can assume that water molecules are rotating/changing direction very slow, on the other hand, if WOR decay very fast, we can assume that water molecules are rotating/changing direction very fast:: +Analyzing water orientational relaxation (WOR) :class:`WaterOrientationalRelaxation`. In this case we are analyzing "how fast" water molecules are rotating/changing direction. If WOR is very stable we can assume that water molecules are rotating/changing direction very slow, on the other hand, if WOR decay very fast, we can assume that water molecules are rotating/changing direction very fast:: import MDAnalysis - from MDAnalysis.analysis.waterdynamics import WOR + from MDAnalysis.analysis.waterdynamics import WaterOrientationalRelaxation as WOR u = MDAnalysis.Universe(pdb, trajectory) selection = "byres name OH2 and sphzone 6.0 protein and resid 42" @@ -110,17 +109,17 @@ the first window is created with 1000 timestep average (1000/1), the second window is created with 500 timestep average(1000/2), the third window is created with 333 timestep average (1000/3) and so on. -AD -~~ +AngularDistribution +~~~~~~~~~~~~~~~~~~~ -Analyzing angular distribution :class:`AD` for OH vector, HH vector and dipole vector. It returns +Analyzing angular distribution (AD) :class:`AngularDistribution` for OH vector, HH vector and dipole vector. It returns a line histogram with vector orientation preference. A straight line in the output graphic means no preferential orientation in water molecules. In this case we are analyzing if water molecules have some orientational preference, in this way we can see if water molecules are under an electric field or if they are interacting with something (residue, protein, etc):: import MDAnalysis - from MDAnalysis.analysis.waterdynamics import AD + from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD u = MDAnalysis.Universe(pdb, trajectory) selection = "byres name OH2 and sphzone 6.0 (protein and (resid 42 or resid 26) )" @@ -135,14 +134,14 @@ where P(cos(theta)) is the angular distribution or angular probabilities. -MSD -~~~ -Analyzing mean square displacement :class:`MSD` for water molecules. In this case we are analyzing the average distance +MeanSquareDisplacement +~~~~~~~~~~~~~~~~~~~~~~ +Analyzing mean square displacement (MSD):class:`MeanSquareDisplacement` for water molecules. In this case we are analyzing the average distance that water molecules travels inside protein in XYZ direction (cylindric zone of radius 11[nm], Zmax 4.0[nm] and Zmin -8.0[nm]). A strong rise mean a fast movement of water molecules, a weak rise mean slow movement of particles:: import MDAnalysis - from MDAnalysis.analysis.waterdynamics import MSD + from MDAnalysis.analysis.waterdynamics import MeanSquareDisplacement as MSD u = MDAnalysis.Universe(pdb, trajectory) selection = "byres name OH2 and cyzone 11.0 4.0 -8.0 protein" @@ -156,15 +155,15 @@ i += 1 -SP -~~ -Analyzing survival probability :class:`SP` for water molecules. In this case we are analyzing how long water +SurvivalProbability +~~~~~~~~~~~~~~~~~~~ +Analyzing survival probability (SP) :class:`SurvivalProbability` for water molecules. In this case we are analyzing how long water molecules remain in a sphere of radius 12.3 centered in the geometrical center of resid 42, 26, 34 and 80. A slow decay of SP means a long permanence time of water molecules in the zone, on the other hand, a fast decay means a short permanence time:: import MDAnalysis - from MDAnalysis.analysis.waterdynamics import SP + from MDAnalysis.analysis.waterdynamics import SurvivalProbability as SP u = MDAnalysis.Universe(pdb, trajectory) selection = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26 or resid 34 or resid 80) " @@ -182,10 +181,10 @@ Output ------ -HBL -~~~ -Hydrogen bond lifetime data is returned per window timestep, which is stored in -:attr:`HBL.timeseries` (in all the following descriptions, # indicates comments that are not part of the output):: +HydrogenBondLifetimes +~~~~~~~~~~~~~~~~~~~~~ +Hydrogen bond lifetimes (HBL) data is returned per window timestep, which is stored in +:attr:`HydrogenBondLifetimes.timeseries` (in all the following descriptions, # indicates comments that are not part of the output):: results = [ [ # time t0 @@ -197,10 +196,10 @@ ... ] -WOR -~~~ -Water orientational relaxation data is returned per window timestep, which is stored in -:attr:`WOR.timeseries`:: +WaterOrientationalRelaxation +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Water orientational relaxation (WOR) data is returned per window timestep, which is stored in +:attr:`WaterOrientationalRelaxation.timeseries`:: results = [ [ # time t0 @@ -212,10 +211,10 @@ ... ] -AD -~~~ -Angular distribution data is returned per vector, which is stored in -:attr:`AD.graph`. In fact, AD returns a histogram:: +AngularDistribution +~~~~~~~~~~~~~~~~~~~ +Angular distribution (AD) data is returned per vector, which is stored in +:attr:`AngularDistribution.graph`. In fact, AngularDistribution returns a histogram:: results = [ [ # OH vector values @@ -230,20 +229,20 @@ ], ] -MSD -~~~ -Mean Square Displacement data is returned in a list, which each element represents a MSD value in its respective -window timestep. Data is stored in :attr:`MSD.timeseries`:: +MeanSquareDisplacement +~~~~~~~~~~~~~~~~~~~~~~ +Mean Square Displacement (MSD) data is returned in a list, which each element represents a MSD value in its respective +window timestep. Data is stored in :attr:`MeanSquareDisplacement.timeseries`:: results = [ #MSD values orders by window timestep , , ... ] -SP -~~ -Survival Probability data is returned in a list, which each element represents a MSD value in its respective -window timestep. Data is stored in :attr:`SP.timeseries`:: +SurvivalProbability +~~~~~~~~~~~~~~~~~~~ +Survival Probability (SP) data is returned in a list, which each element represents a SP value in its respective +window timestep. Data is stored in :attr:`SurvivalProbability.timeseries`:: results = [ # SP values order by window timestep @@ -255,41 +254,37 @@ Classes -------- -.. autoclass:: HBL +.. autoclass:: HydrogenBondLifetimes :members: :inherited-members: -.. autoclass:: WOR +.. autoclass:: WaterOrientationalRelaxation :members: :inherited-members: -.. autoclass:: AD +.. autoclass:: AngularDistribution :members: :inherited-members: -.. autoclass:: MSD +.. autoclass:: MeanSquareDisplacement :members: :inherited-members: -.. autoclass:: SP +.. autoclass:: SurvivalProbability :members: :inherited-members: """ -from MDAnalysis import * + import MDAnalysis.analysis.hbonds -import MDAnalysis.core.Selection import numpy -import argparse -import sys import multiprocessing -import time import itertools -class HBL: +class HydrogenBondLifetimes: r""" - This is a autocorrelation function that gives the "Hydrogen Bond Lifetimes" proposed by D.C. Rapaport [Rapaport1983]_. From this + This is a autocorrelation function that gives the "Hydrogen Bond Lifetimes" (HBL) proposed by D.C. Rapaport [Rapaport1983]_. From this function we can obtain the continuos and intermittent behavior of hydrogen bonds in time. A fast decay in these parameters indicate a fast change in HBs connectivity. A slow decay indicate very stables hydrogen bonds, like in ice. The HBL is also know as "Hydrogen Bond Population @@ -510,10 +505,10 @@ def run(self,**kwargs): h_list.run() self.timeseries = self._getGraphics( h_list.timeseries , self.t0 , self.tf , self.dtmax ) -class WOR: +class WaterOrientationalRelaxation: r""" Function to evaluate the Water Orientational Relaxation proposed by Yu-ling Yeh - and Chung-Yuan Mou [Yeh1999_]. WOR indicates "how fast" water molecules are rotating + and Chung-Yuan Mou [Yeh1999_]. WaterOrientationalRelaxation indicates "how fast" water molecules are rotating or changing direction. This is a time correlation function given by: .. math:: @@ -546,10 +541,10 @@ def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): self.nproc = nproc self.timeseries = None - def _repitedIndex(self,selection,dt,totalFrames): + def _repeatedIndex(self,selection,dt,totalFrames): """ Indicate the comparation between all the t+dt. - The results is a list of list with all the repited index per frame (or time). + The results is a list of list with all the repeated index per frame (or time). Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... @@ -560,12 +555,12 @@ def _repitedIndex(self,selection,dt,totalFrames): rep.append(self._sameMolecTandDT(selection,dt*i,(dt*i)+dt)) return rep - def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): + def _getOneDeltaPoint(self,universe, repInd, i ,t0, dt): """ Give one point to promediate and get one point of the graphic C_vect vs t Ex: t0=1 and tau=1 so calculate the t0-tau=1-2 intervale. Ex: t0=5 and tau=3 so calcultate the t0-tau=5-8 intervale. - i = come from getMeanOneCOHPoint (named j) (int) + i = come from getMeanOnePoint (named j) (int) """ valOH = 0 valHH = 0 @@ -613,13 +608,13 @@ def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): valdip = valdip/n return (valOH,valHH,valdip) - def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): + def _getMeanOnePoint(self,universe,selection1,selection_str,dt,totalFrames): """ This function get one point of the graphic C_OH vs t. It uses the - _getOneDeltaCOHPoint() function to calculate the average. + _getOneDeltaPoint() function to calculate the average. """ - repInd = self._repitedIndex(selection1,dt,totalFrames) + repInd = self._repeatedIndex(selection1,dt,totalFrames) sumsdt = 0 n = 0.0 sumDeltaOH = 0.0 @@ -630,7 +625,7 @@ def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): valdipList = [] for j in range(totalFrames/dt-1): - a = self._getOneDeltaCOHPoint(universe,repInd,j,sumsdt,dt) + a = self._getOneDeltaPoint(universe,repInd,j,sumsdt,dt) sumDeltaOH += a[0] sumDeltaHH += a[1] sumDeltadip += a[2] @@ -644,7 +639,7 @@ def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): def _sameMolecTandDT(self,selection,t0d,tf): """ Compare the molecules in the t0d selection and the t0d+dt selection and - select only the particles that are repited in both frame. This is to consider + select only the particles that are repeated in both frame. This is to consider only the molecules that remains in the selection after the dt time has elapsed. The result is a list with the indexs of the atoms. """ @@ -677,12 +672,12 @@ def run(self,**kwargs): selection_out = self._selection_serial(self.universe,self.selection) self.timeseries = [] for dt in list(range(1,self.dtmax+1)): - output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + output = self._getMeanOnePoint(self.universe,selection_out,self.selection,dt,self.tf) self.timeseries.append(output) -class AD: +class AngularDistribution: r""" The angular distribution function (AD) is defined as the distribution probability of the cosine of the :math:`\theta` angle formed by the OH vector, HH vector @@ -717,45 +712,42 @@ def _getCosTheta(self,universe,selection,axis): valOH = [] valHH = [] valdip= [] + i = 0 while i <= (len(selection)-1): universe.trajectory[i] line = selection[ i ].positions - j=0 - while j < len(line): - Ot0 = numpy.array(line[j]) - H1t0 = numpy.array(line[j+1]) - H2t0 = numpy.array(line[j+2]) - - OHVector0 = H1t0 - Ot0 - HHVector0 = H1t0 - H2t0 - dipVector0 = (H1t0 + H2t0)*0.5 - Ot0 - normOHVector0 = numpy.linalg.norm(OHVector0) - normHHVector0 = numpy.linalg.norm(HHVector0) - normdipVector0 = numpy.linalg.norm(dipVector0) - - unitOHVector0 = [OHVector0[0]/normOHVector0,OHVector0[1]/normOHVector0,OHVector0[2]/normOHVector0] - unitHHVector0 = [HHVector0[0]/normHHVector0,HHVector0[1]/normHHVector0,HHVector0[2]/normHHVector0] - unitdipVector0 = [dipVector0[0]/normdipVector0,dipVector0[1]/normdipVector0,dipVector0[2]/normdipVector0] - + Ot0 = line[::3] + H1t0 = line[1::3] + H2t0 = line[2::3] + + OHVector0 = H1t0 - Ot0 + HHVector0 = H1t0 - H2t0 + dipVector0 = (H1t0 + H2t0)*0.5 - Ot0 + + unitOHVector0 = OHVector0/numpy.linalg.norm(OHVector0, axis = 1)[:,None] + unitHHVector0 = HHVector0/numpy.linalg.norm(HHVector0, axis = 1)[:,None] + unitdipVector0 = dipVector0/numpy.linalg.norm(dipVector0, axis = 1)[:,None] + + j=0 + while j < len(line)/3: if axis == "z": - valOH.append(numpy.dot(unitOHVector0,[0,0,1])) - valHH.append(numpy.dot(unitHHVector0,[0,0,1])) - valdip.append(numpy.dot(unitdipVector0,[0,0,1])) + valOH.append(unitOHVector0[j][2]) + valHH.append(unitHHVector0[j][2]) + valdip.append(unitdipVector0[j][2]) elif axis == "x": - valOH.append(numpy.dot(unitOHVector0,[1,0,0])) - valHH.append(numpy.dot(unitHHVector0,[1,0,0])) - valdip.append(numpy.dot(unitdipVector0,[1,0,0])) + valOH.append(unitOHVector0[j][0]) + valHH.append(unitHHVector0[j][0]) + valdip.append(unitdipVector0[j][0]) elif axis == "y": - valOH.append(numpy.dot(unitOHVector0,[0,1,0])) - valHH.append(numpy.dot(unitHHVector0,[0,1,0])) - valdip.append(numpy.dot(unitdipVector0,[0,1,0])) - - - j += 3 + valOH.append(unitOHVector0[j][1]) + valHH.append(unitHHVector0[j][1]) + valdip.append(unitdipVector0[j][1]) + + j += 1 i += 1 return (valOH,valHH,valdip) @@ -819,7 +811,7 @@ def _selection_serial(self,universe,selection_str): return selection -class MSD: +class MeanSquareDisplacement: r""" Function to evaluate the Mean Square Displacement (MSD_). The MSD gives the average distance that particles travels. The MSD is given by: @@ -856,10 +848,10 @@ def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): self.nproc = nproc self.timeseries = None - def _repitedIndex(self,selection,dt,totalFrames): + def _repeatedIndex(self,selection,dt,totalFrames): """ Indicate the comparation between all the t+dt. - The results is a list of list with all the repited index per frame (or time). + The results is a list of list with all the repeated index per frame (or time). Ex: dt=1, so compare frames (1,2),(2,3),(3,4)... Ex: dt=2, so compare frames (1,3),(3,5),(5,7)... Ex: dt=3, so compare frames (1,4),(4,7),(7,10)... @@ -870,12 +862,12 @@ def _repitedIndex(self,selection,dt,totalFrames): rep.append(self._sameMolecTandDT(selection,dt*i,(dt*i)+dt)) return rep - def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): + def _getOneDeltaPoint(self,universe, repInd, i ,t0, dt): """ Give one point to promediate and get one point of the grapic C_vect vs t Ex: t0=1 and dt=1 so calculate the t0-dt=1-2 intervale. Ex: t0=5 and dt=3 so calcultate the t0-dt=5-8 intervale - i = come from getMeanOneCOHPoint (named j) (int) + i = come from getMeanOnePoint (named j) (int) """ valO = 0 n = 0 @@ -891,26 +883,26 @@ def _getOneDeltaCOHPoint(self,universe, repInd, i ,t0, dt): #position oxygen OVector = Ot0 - Otp - #here it is the difference with waterdynamics.WOR + #here it is the difference with waterdynamics.WaterOrientationalRelaxation valO += numpy.dot(OVector, OVector) n += 1 valO = valO/n return (valO) - def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): + def _getMeanOnePoint(self,universe,selection1,selection_str,dt,totalFrames): """ This function get one point of the graphic C_OH vs t. It's uses the - _getOneDeltaCOHPoint() function to calculate the average. + _getOneDeltaPoint() function to calculate the average. """ - repInd = self._repitedIndex(selection1,dt,totalFrames) + repInd = self._repeatedIndex(selection1,dt,totalFrames) sumsdt = 0 n = 0.0 sumDeltaO = 0.0 valOList = [] for j in range(totalFrames/dt-1): - a = self._getOneDeltaCOHPoint(universe,repInd,j,sumsdt,dt) + a = self._getOneDeltaPoint(universe,repInd,j,sumsdt,dt) print "a=",a sumDeltaO += a valOList.append(a) @@ -921,7 +913,7 @@ def _getMeanOneCOHPoint(self,universe,selection1,selection_str,dt,totalFrames): def _sameMolecTandDT(self,selection,t0d,tf): """ Compare the molecules in the t0d selection and the t0d+dt selection and - select only the particles that are repited in both frame. This is to consider + select only the particles that are repeated in both frame. This is to consider only the molecules that remains in the selection after the dt time has elapsed. The result is a list with the indexs of the atoms. """ @@ -951,31 +943,34 @@ def run(self,**kwargs): selection_out = self._selection_serial(self.universe,self.selection) self.timeseries = [] for dt in list(range(1,self.dtmax+1)): - output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + output = self._getMeanOnePoint(self.universe,selection_out,self.selection,dt,self.tf) self.timeseries.append(output) -class SP: +class SurvivalProbability: r""" Function to evaluate the Survival Probability (SP). The SP gives the probability for a group of particles to remain in certain region. The SP is given by: - SP = equation .. math:: P(\tau) = \frac1T \sum_{t=1}^T \frac{N(t,t+\tau)}{N(t)} + where :math:`T` is the maximum time of simulation, :math:`\tau` is the timestep and + :math:`N` the number of particles in certain time. + + :Arguments: - *universe* - Universe object - *selection* - Selection string, any selection is allowed, with this selection you define the region/zone where - to analize, i.e.: "selection_a" and "zone" (see SP examples_ ) - *t0* - Time where analysis begin - *tf* - Time where analysis end - *dtmax* - Maximum dt size window, dtmax < tf or it will crash + *universe* + Universe object + *selection* + Selection string, any selection is allowed, with this selection you define the region/zone where + to analize, i.e.: "selection_a" and "zone" (see SP examples_ ) + *t0* + Time where analysis begin + *tf* + Time where analysis end + *dtmax* + Maximum dt size window, dtmax < tf or it will crash """ @@ -988,7 +983,7 @@ def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): self.nproc = nproc self.timeseries = None - def _getOneDeltaCOHPoint(self,selection, totalFrames, t0, tau): + def _getOneDeltaPoint(self,selection, totalFrames, t0, tau): """ Give one point to promediate and get one point of the graphic C_vect vs t Ex: t0=1 and tau=1 so calculate the t0-tau=1-2 intervale. @@ -999,17 +994,16 @@ def _getOneDeltaCOHPoint(self,selection, totalFrames, t0, tau): return Ntau/Nt - def _getMeanOneCOHPoint(self,universe,selection1,selection_str,wint,totalFrames): + def _getMeanOnePoint(self,universe,selection1,selection_str,wint,totalFrames): """ This function get one point of the graphic P(t) vs t. It uses the - _getOneDeltaCOHPoint() function to calculate the average. + _getOneDeltaPoint() function to calculate the average. """ n = 0.0 sumDeltaP = 0.0 - valPList = [] for frame in range(totalFrames-wint): - a = self._getOneDeltaCOHPoint(selection1,totalFrames ,frame, wint) + a = self._getOneDeltaPoint(selection1,totalFrames ,frame, wint) sumDeltaP += a n += 1 @@ -1046,12 +1040,12 @@ def run(self,**kwargs): #All the selection to an array, this way is faster than selecting later. if self.nproc==1: - selection_out = self._selection_serial(self.universe,self.selection) + selection_out = self._selection_serial(self.universe, self.selection) else: #selection = selection_parallel(universe,selection_str,nproc) #parallel selection to be implemented - selection_out = self._selection_serial(self.universe,self.selection) + selection_out = self._selection_serial(self.universe, self.selection) self.timeseries = [] for dt in list(range(1,self.dtmax+1)): - output = self._getMeanOneCOHPoint(self.universe,selection_out,self.selection,dt,self.tf) + output = self._getMeanOnePoint(self.universe, selection_out, self.selection, dt, self.tf) self.timeseries.append(output) diff --git a/testsuite/MDAnalysisTests/test_analysis.py b/testsuite/MDAnalysisTests/test_analysis.py index dd983a99a29..5a4684924c6 100644 --- a/testsuite/MDAnalysisTests/test_analysis.py +++ b/testsuite/MDAnalysisTests/test_analysis.py @@ -368,30 +368,28 @@ def setUp(self): self.selection1 = "byres name OH2" self.selection2 = self.selection1 - def test_HBL(self): - hbl = MDAnalysis.analysis.waterdynamics.HBL(self.universe, self.selection1, self.selection2, 0, 5, 3) + def test_HydrogenBondLifetimes(self): + hbl = MDAnalysis.analysis.waterdynamics.HydrogenBondLifetimes(self.universe, self.selection1, self.selection2, 0, 5, 3) hbl.run() assert_equal(round(hbl.timeseries[2][1],5), 0.75) - def test_WOR(self): - wor = MDAnalysis.analysis.waterdynamics.WOR(self.universe, self.selection1, 0, 5, 2) + def test_WaterOrientationalRelaxation(self): + wor = MDAnalysis.analysis.waterdynamics.WaterOrientationalRelaxation(self.universe, self.selection1, 0, 5, 2) wor.run() assert_equal(round(wor.timeseries[1][2],5), 0.35887) - def test_AD(self): - ad = MDAnalysis.analysis.waterdynamics.AD(self.universe,self.selection1,40) + def test_AngularDistribution(self): + ad = MDAnalysis.analysis.waterdynamics.AngularDistribution(self.universe,self.selection1,40) ad.run() - #debo arreglar la salida del AD - #must fix AD output assert_equal(str(ad.graph[0][39]), str("0.951172947884 0.48313682125") ) - def test_MSD(self): - msd = MDAnalysis.analysis.waterdynamics.MSD(self.universe, self.selection1, 0, 10, 2) + def test_MeanSquareDisplacement(self): + msd = MDAnalysis.analysis.waterdynamics.MeanSquareDisplacement(self.universe, self.selection1, 0, 10, 2) msd.run() assert_equal(round(msd.timeseries[1],5), 0.03984) - def test_SP(self): - sp = MDAnalysis.analysis.waterdynamics.SP(self.universe, self.selection1, 0, 6, 3) + def test_SurvivalProbability(self): + sp = MDAnalysis.analysis.waterdynamics.SurvivalProbability(self.universe, self.selection1, 0, 6, 3) sp.run() assert_equal(round(sp.timeseries[1],5), 1.0) From b27bcda3227074378c1c26f18e6b586f21d7c0bd Mon Sep 17 00:00:00 2001 From: Alejandro Bernardin Date: Tue, 9 Jun 2015 16:37:18 -0300 Subject: [PATCH 3/4] Solved merge conflict. --- testsuite/MDAnalysisTests/datafiles.py | 6 ++---- testsuite/MDAnalysisTests/test_analysis.py | 10 ---------- 2 files changed, 2 insertions(+), 14 deletions(-) diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index d63b9567be0..d904448ae1d 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -79,7 +79,7 @@ "GMS_ASYMSURF", # GAMESS C1 surface "two_water_gro", "two_water_gro_nonames", # for bond guessing, 2 water molecules, one with weird names "DLP_CONFIG", "DLP_CONFIG_order", "DLP_CONFIG_minimal", # dl_poly 4 config file - "DLP_HISTORY", "DLP_HISTORY_order", "DLP_HISTORY_minimal" # dl_poly 4 history file + "DLP_HISTORY", "DLP_HISTORY_order", "DLP_HISTORY_minimal", # dl_poly 4 history file "waterPSF","waterDCD", ] @@ -212,14 +212,12 @@ two_water_gro = resource_filename(__name__, "data/two_water_gro.gro") two_water_gro_nonames = resource_filename(__name__, "data/two_water_gro_nonames.gro") -<<<<<<< HEAD DLP_CONFIG = resource_filename(__name__, "data/dlpoly/CONFIG") DLP_CONFIG_order = resource_filename(__name__, "data/dlpoly/CONFIG_order") DLP_CONFIG_minimal = resource_filename(__name__, "data/dlpoly/CONFIG_minimal") DLP_HISTORY = resource_filename(__name__, "data/dlpoly/HISTORY") DLP_HISTORY_order = resource_filename(__name__, "data/dlpoly/HISTORY_order") DLP_HISTORY_minimal = resource_filename(__name__, "data/dlpoly/HISTORY_minimal") -======= + waterPSF = resource_filename(__name__, 'data/watdyn.psf') waterDCD = resource_filename(__name__, 'data/watdyn.dcd') ->>>>>>> feature-watdyn diff --git a/testsuite/MDAnalysisTests/test_analysis.py b/testsuite/MDAnalysisTests/test_analysis.py index f361741326c..8bcc87ed597 100644 --- a/testsuite/MDAnalysisTests/test_analysis.py +++ b/testsuite/MDAnalysisTests/test_analysis.py @@ -356,8 +356,6 @@ def test_xtc_striding(self): """Check for sustained resolution of Issue 188.""" u = self.universe sel = self.selection -<<<<<<< HEAD - assert_raises(FinishTimeException, MDAnalysis.analysis.helanal.helanal_trajectory, u, selection=sel, finish=5) @@ -366,13 +364,6 @@ def test_xtc_striding(self): # MDAnalysis.analysis.helanal.helanal_trajectory(u, selection=sel, finish=5) # except IndexError: # self.fail("IndexError consistent with Issue 188.") -======= - with assert_raises(FinishTimeException): - try: - MDAnalysis.analysis.helanal.helanal_trajectory(u, selection=sel, finish=5) - except IndexError: - self.fail("IndexError consistent with Issue 188.") - class TestWaterdynamics(TestCase): def setUp(self): @@ -405,4 +396,3 @@ def test_SurvivalProbability(self): sp.run() assert_equal(round(sp.timeseries[1],5), 1.0) ->>>>>>> feature-watdyn From eacddb5edf728cca77feebfe63776fb930a37b99 Mon Sep 17 00:00:00 2001 From: Alejandro Bernardin Date: Wed, 10 Jun 2015 16:21:48 -0300 Subject: [PATCH 4/4] Added version information and CHANGELOG update. --- package/CHANGELOG | 14 ++++++++++++++ package/MDAnalysis/analysis/waterdynamics.py | 11 +++++++++++ 2 files changed, 25 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index c5b4ee79607..3e528c207cb 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -10,6 +10,20 @@ The rules for this file: use tabs but use spaces for formatting ------------------------------------------------------------------------------ +??/??/15 alejob + + * 0.11.0 + + Enhancements + * Waterdynamics analysis module added, including five analysis classes: Hydrogen Bond + Lifetimes, Water Orientational Relaxation, Angular Distribution, Mean Square + Displacement and Survival Probability. Documentation and test are included. + + Changes + + Fixes + + ??/??/15 richardjgowers * 0.10.1 diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index 4d0a7d183d7..9510b3d1d8e 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -22,6 +22,8 @@ :Year: 2014-2015 :Copyright: GNU Public License v3 +.. versionadded:: 0.11.0 + This module provides functions to analize water dynamics trajectories and water interactions with other molecules. The functions in this module are: water orientational relaxation (WOR) [Yeh1999]_, hydrogen bond lifetimes (HBL) [Rapaport1983]_, angular distribution (AD) [Grigera1995]_, mean square displacement (MSD) [Brodka1994]_ and survival probability (SP) [Liu2004]_. @@ -87,6 +89,7 @@ hydrogen bond lifetime, t0 = 0, tf = 2000 and dtmax = 30. In this way we create 30 windows timestep (30 values in x axis). The continuos hydrogen bond lifetimes should decay faster than intermittent. + WaterOrientationalRelaxation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -302,6 +305,8 @@ class HydrogenBondLifetimes: where :math:`h_{ij}(t_0+\tau)=1` if there is a H-bond between a pair :math:`ij` at time :math:`t_0+\tau` (intermittent) and :math:`h_{ij}(t_0+\tau)=0` otherwise. + + .. versionadded:: 0.11.0 :Arguments: *universe* @@ -517,6 +522,7 @@ class WaterOrientationalRelaxation: where :math:`P_2=(3x^2-1)/2` is the second-order Legendre polynomial and :math:`\hat{u}` is a unit vector along HH, OH or dipole vector. + .. versionadded:: 0.11.0 :Arguments: *universe* @@ -685,6 +691,8 @@ class AngularDistribution: (z is the default value). The cosine is define as :math:`\cos \theta = \hat u \cdot \hat n`, where :math:`\hat u` is OH, HH or dipole vector. It creates a histogram and returns a list of lists, see Output_. The AD is also know as Angular Probability (AP). + .. versionadded:: 0.11.0 + :Arguments: *universe* Universe object @@ -825,6 +833,8 @@ class MeanSquareDisplacement: .. _MSD: http://en.wikipedia.org/wiki/Mean_squared_displacement + .. versionadded:: 0.11.0 + :Arguments: *universe* Universe object @@ -958,6 +968,7 @@ class SurvivalProbability: where :math:`T` is the maximum time of simulation, :math:`\tau` is the timestep and :math:`N` the number of particles in certain time. + .. versionadded:: 0.11.0 :Arguments: *universe*