From e54ff8903c7e5c59ba01c6ced0d0ba24b8e414d9 Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Thu, 9 Mar 2017 09:56:29 -0700 Subject: [PATCH 1/9] analysis.density: fixed top level Example section in docs Do not use Examples as a heading UNLESS inside a function/class doc because the NumPy reST parser changes it to a rubric heading. This breaks document structure. --- package/MDAnalysis/analysis/rms.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/analysis/rms.py b/package/MDAnalysis/analysis/rms.py index 14fb5ce9f1e..9b9ded629ba 100644 --- a/package/MDAnalysis/analysis/rms.py +++ b/package/MDAnalysis/analysis/rms.py @@ -51,8 +51,8 @@ :mod:`MDAnalysis.lib.qcprot` implements the fast RMSD algorithm. -Examples --------- +Example applications +-------------------- Calculating RMSD for multiple domains ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 3da846f8ede9b0d3e8128702c14ccc139455b8ad Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Thu, 23 Mar 2017 03:18:32 -0700 Subject: [PATCH 2/9] analysis.psa: doc fixes - converted all docs to numpy style - added additional references - see also: @tylerjereddy 's scipy.spatial.distance.directed_hausdorff() --- package/MDAnalysis/analysis/psa.py | 928 ++++++++++++++++++----------- 1 file changed, 575 insertions(+), 353 deletions(-) diff --git a/package/MDAnalysis/analysis/psa.py b/package/MDAnalysis/analysis/psa.py index 4850c24cf74..84510728952 100644 --- a/package/MDAnalysis/analysis/psa.py +++ b/package/MDAnalysis/analysis/psa.py @@ -251,50 +251,67 @@ def get_path_metric_func(name): def sqnorm(v, axis=None): """Compute the sum of squares of elements along specified axes. - :Arguments: - *v* - :class:`numpy.ndarray` of coordinates - *axes* - None or int or tuple of ints, optional - Axes or axes along which a sum is performed. The default (*axes* = - ``None``) performs a sum over all the dimensions of the input array. - The value of *axes* may be negative, in which case it counts from the - last axis to the zeroth axis. + Parameters + ---------- + v : numpy.ndarray + coordinates + axes : None / int / tuple (optional) + Axes or axes along which a sum is performed. The default + (*axes* = ``None``) performs a sum over all the dimensions of + the input array. The value of *axes* may be negative, in + which case it counts from the last axis to the zeroth axis. + + Returns + ------- + float + the sum of the squares of the elements of `v` along `axes` - :Returns: - float, the sum of the squares of the elements of *v* along *axes* """ return np.sum(v*v, axis=axis) def get_msd_matrix(P, Q, axis=None): - """Generate the matrix of pairwise mean-squared deviations (MSDs) between - all pairs of points in *P* and *Q*, each pair having a point from *P* and a - point from *Q*. + r"""Generate the matrix of pairwise mean-squared deviations between paths. - *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time + The MSDs between all pairs of points in `P` and `Q` are + calculated, each pair having a point from `P` and a point from + `Q`. + + `P` (`Q`) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). The pairwise MSD + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). The pairwise MSD matrix has dimensions :math:`N_p` by :math:`N_q`. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path - :Returns: - :class:`numpy.ndarray` of pairwise MSDs between points in *P* and points - in *Q* + Returns + ------- + msd_matrix : numpy.ndarray + matrix of pairwise MSDs between points in `P` and points + in `Q` + + Notes + ----- + We calculate the MSD matrix + + .. math:: + M_{ij} = ||p_i - q_j||^2 + + where :math:`p_i \in P` and :math:`q_j \in Q`. """ return np.asarray([sqnorm(p - Q, axis=axis) for p in P]) def get_coord_axes(path): - """Return the number of atoms and the axes (*axis*) corresponding to atoms + """Return the number of atoms and the axes corresponding to atoms and coordinates for a given path. - The *path* is assumed to be a :class:`numpy.ndarray` where the 0th axis + The `path` is assumed to be a :class:`numpy.ndarray` where the 0th axis corresponds to a frame (a snapshot of coordinates). The :math:`3N` (Cartesian) coordinates are assumed to be either: @@ -304,13 +321,15 @@ def get_coord_axes(path): number and the 2nd axis contains the *x*,*y*,*z* coordinates of each atom. - :Arguments: - *path* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + path : numpy.ndarray + representing a path - :Returns: - ``(int, (int, ...))``, the number of atoms and the axes - containing coordinates + Returns + ------- + (int, (int, ...)) + the number of atoms and the axes containing coordinates """ path_dimensions = len(path.shape) @@ -330,23 +349,29 @@ def get_coord_axes(path): def hausdorff(P, Q): - r"""Calculate the Hausdorff distance between two paths. + r"""Calculate the symmetric Hausdorff distance between two paths. *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). *P* (*Q*) has + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). *P* (*Q*) has either shape |3Dp| (|3Dq|), or |2Dp| (|2Dq|) in flattened form. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path - :Returns: - float, the Hausdorff distance between paths *P* and *Q* + Returns + ------- + float + the Hausdorff distance between paths `P` and `Q` + + Example + ------- + Calculate the Hausdorff distance between two halves of a trajectory: - Example:: >>> from MDAnalysis.tests.datafiles import PSF, DCD >>> u = Universe(PSF,DCD) >>> mid = len(u.trajectory)/2 @@ -361,6 +386,34 @@ def hausdorff(P, Q): 4.7786639840135905 >>> hausdorff(P,Q[::-1]) # hausdorff distance w/ reversed 2nd trajectory 4.7786639840135905 + + Note that reversing the path does not change the Hausdorff distance. + + Notes + ----- + The Hausdorff distance is calculated in a brute force manner from the + distance matrix without further optimizations, essentially following + [Huttenlocher1993]_. + + .. SeeAlso:: + + :func:`scipy.spatial.distance.directed_hausdorff` is an optimized + implementation of the early break algorithm of [Taha2015]_; note that + one still has to calculate the *symmetric* Hausdorff distance as + `max(directed_hausdorff(P, Q)[0], directed_hausdorff(Q, P)[0])`. + + + References + ---------- + .. [Huttenlocher1993] D. P. Huttenlocher, G. A. Klanderman, and + W. J. Rucklidge. Comparing images using the Hausdorff distance. IEEE + Transactions on Pattern Analysis and Machine Intelligence, + 15(9):850–863, 1993. + + .. [Taha2015] A. A. Taha and A. Hanbury. An efficient algorithm for + calculating the exact Hausdorff distance. IEEE Transactions On Pattern + Analysis And Machine Intelligence, 37:2153-63, 2015. + """ N, axis = get_coord_axes(P) d = get_msd_matrix(P, Q, axis=axis) @@ -373,23 +426,28 @@ def hausdorff_wavg(P, Q): *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). *P* (*Q*) has + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). *P* (*Q*) has either shape |3Dp| (|3Dq|), or |2Dp| (|2Dq|) in flattened form. The nearest neighbor distances for *P* (to *Q*) and those of *Q* (to *P*) are averaged individually to get the average nearest neighbor distance for *P* and likewise for *Q*. These averages are then summed and divided by 2 to get a measure that gives equal weight to *P* and *Q*. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path - :Returns: - float, the weighted average Hausdorff distance between paths *P* and *Q* + Returns + ------- + float + the weighted average Hausdorff distance between paths `P` and `Q` + + Example + ------- - Example:: >>> from MDAnalysis import Universe >>> from MDAnalysis.tests.datafiles import PSF, DCD >>> u = Universe(PSF,DCD) @@ -405,6 +463,13 @@ def hausdorff_wavg(P, Q): 2.5669644353703447 >>> hausdorff_wavg(P,Q[::-1]) # weighted avg hausdorff dist w/ Q reversed 2.5669644353703447 + + Notes + ----- + The weighted average Hausdorff distance is not a true metric (it does not + obey the triangle inequality); see [Seyler2015]_ for further details. + + """ N, axis = get_coord_axes(P) d = get_msd_matrix(P, Q, axis=axis) @@ -417,23 +482,28 @@ def hausdorff_avg(P, Q): *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). *P* (*Q*) has + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). *P* (*Q*) has either shape |3Dp| (|3Dq|), or |2Dp| (|2Dq|) in flattened form. The nearest neighbor distances for *P* (to *Q*) and those of *Q* (to *P*) are all averaged together to get a mean nearest neighbor distance. This measure biases the average toward the path that has more snapshots, whereas weighted average Hausdorff gives equal weight to both paths. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path - :Returns: - float, the average Hausdorff distance between paths *P* and *Q* + Returns + ------- + float + the average Hausdorff distance between paths `P` and `Q` + + Example + ------- - Example:: >>> from MDAnalysis.tests.datafiles import PSF, DCD >>> u = Universe(PSF,DCD) >>> mid = len(u.trajectory)/2 @@ -448,6 +518,13 @@ def hausdorff_avg(P, Q): 2.5669646575869005 >>> hausdorff_avg(P,Q[::-1]) # hausdorff distance w/ reversed 2nd trajectory 2.5669646575869005 + + + Notes + ----- + The average Hausdorff distance is not a true metric (it does not obey the + triangle inequality); see [Seyler2015]_ for further details. + """ N, axis = get_coord_axes(P) d = get_msd_matrix(P, Q, axis=axis) @@ -456,54 +533,75 @@ def hausdorff_avg(P, Q): def hausdorff_neighbors(P, Q): - r"""Calculate the Hausdorff distance between two paths. + r"""Find the Hausdorff neighbors of two paths. *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). *P* (*Q*) has + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). *P* (*Q*) has either shape |3Dp| (|3Dq|), or |2Dp| (|2Dq|) in flattened form. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path + + Returns + ------- + dict + dictionary of two pairs of numpy arrays, the first pair (key + "frames") containing the indices of (Hausdorff) nearest + neighbors for `P` and `Q`, respectively, the second (key + "distances") containing (corresponding) nearest neighbor + distances for `P` and `Q`, respectively + + Notes + ----- + Hausdorff neighbors are those points on the two paths that are separated by + the Hausdorff distance. They are the farthest nearest neighbors and are + maximally different in the sense of the Hausdorff distance [Seyler2015]_. + + .. SeeAlso:: + + :func:`scipy.spatial.distance.directed_hausdorff` can also provide the + Hausdorff neighbors. - :Returns: - dictionary of two pairs of numpy arrays, the first pair containing the - indices of (Hausdorff) nearest neighbors for *P* and *Q*, respectively, the - second containing (corresponding) nearest neighbor distances for *P* and - *Q*, respectively """ N, axis = get_coord_axes(P) d = get_msd_matrix(P, Q, axis=axis) - nearest_neighbors = \ - {'frames' : \ - (np.argmin(d, axis=1), np.argmin(d, axis=0)), \ - 'distances' : \ - ((np.amin(d,axis=1)/N)**0.5, (np.amin(d, axis=0)/N)**0.5), \ - } + nearest_neighbors = { + 'frames' : (np.argmin(d, axis=1), np.argmin(d, axis=0)), + 'distances' : ((np.amin(d,axis=1)/N)**0.5, (np.amin(d, axis=0)/N)**0.5) + } return nearest_neighbors def discrete_frechet(P, Q): - r"""Calculate the discrete Frechet distance between two paths. + r"""Calculate the discrete Fréchet distance between two paths. *P* (*Q*) is a :class:`numpy.ndarray` of :math:`N_p` (:math:`N_q`) time steps, :math:`N` atoms, and :math:`3N` coordinates (e.g., - :meth:`MDAnalysis.core.groups.AtomGroup.coordinates`). *P* (*Q*) has + :attr:`MDAnalysis.core.groups.AtomGroup.positions`). *P* (*Q*) has either shape |3Dp| (|3Dq|), or :|2Dp| (|2Dq|) in flattened form. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path + Parameters + ---------- + P : numpy.ndarray + the points in the first path + Q : numpy.ndarray + the points in the second path - :Returns: - float, the discrete Frechet distance between paths *P* and *Q* + Returns + ------- + float + the discrete Fréchet distance between paths *P* and *Q* + + Example + ------- + Calculate the discrete Fréchet distance between two halves of a + trajectory. - Example:: >>> u = Universe(PSF,DCD) >>> mid = len(u.trajectory)/2 >>> ca = u.select_atoms('name CA') @@ -517,7 +615,43 @@ def discrete_frechet(P, Q): 4.7786639840135905 >>> discrete_frechet(P,Q[::-1]) # frechet distance w/ 2nd trj reversed 2nd 6.8429011177113832 + + Note that reversing the direction increased the Fréchet distance: + it is sensitive to the direction of the path. + + Notes + ----- + + The discrete Fréchet metric is an approximation to the continuous Fréchet + metric [Frechet1906]_ [Alt1995]_. The calculation of the continuous + Fréchet distance is implemented with the dynamic programming algorithm of + [EiterMannila1994]_ [EiterMannila1997]_. + + + References + ---------- + + .. [Frechet1906] M. Fréchet. Sur quelques points du calcul + fonctionnel. Rend. Circ. Mat. Palermo, 22(1):1–72, Dec. 1906. + + .. [Alt1995] H. Alt and M. Godau. Computing the Fréchet distance between + two polygonal curves. Int J Comput Geometry & Applications, + 5(01n02):75–91, 1995. doi: `10.1142/S0218195995000064`_ + + .. _`10.1142/S0218195995000064`: http://doi.org/10.1142/S0218195995000064 + + .. [EiterMannila1994] T. Eiter and H. Mannila. Computing discrete Fréchet + distance. Technical Report CD-TR 94/64, Christian Doppler Laboratory for + Expert Systems, Technische Universität Wien, Wien, 1994. + + .. [EiterMannila1997] T. Eiter and H. Mannila. Distance measures for point + sets and their computation. Acta Informatica, 34:109–133, 1997. doi: `10.1007/s002360050075`_. + + + .. _10.1007/s002360050075: http://doi.org/10.1007/s002360050075 + """ + N, axis = get_coord_axes(P) Np, Nq = len(P), len(Q) d = get_msd_matrix(P, Q, axis=axis) @@ -563,19 +697,24 @@ def dist_mat_to_vec(N, i, j): should be greater than *i+1*, corresponding to the upper triangle of the distance matrix. - :Arguments: - *N* - int, size of the distance matrix (of shape *N*-by-*N*) - *i* - int, row index (starting at 0) of the distance matrix - *j* - int, column index (starting at 0) of the distance matrix + Parameters + ---------- + N : int + size of the distance matrix (of shape *N*-by-*N*) + i : int + row index (starting at 0) of the distance matrix + j : int + column index (starting at 0) of the distance matrix + + Returns + ------- + int + index (of the matrix element) in the corresponding distance vector - :Returns: - int, index (of the matrix element) in the corresponding distance vector """ - if not (isinstance(N, numbers.Integral) or isinstance(i, numbers.Integral) or isinstance(j, numbers.Integral)): + if not (isinstance(N, numbers.Integral) or isinstance(i, numbers.Integral) + or isinstance(j, numbers.Integral)): err_str = "N, i, j all must be of type int" raise ValueError(err_str) @@ -602,8 +741,10 @@ def dist_mat_to_vec(N, i, j): class Path(object): - """Pre-process a :class:`MDAnalysis.Universe` object: (1) fit the - trajectory to a reference structure, (2) convert fitted time series to a + """Represent a path based on a :class:`~MDAnalysis.core.universe.Universe`. + + Pre-process a :class:`Universe` object: (1) fit the trajectory to a + reference structure, (2) convert fitted time series to a :class:`numpy.ndarray` representation of :attr:`Path.path`. The analysis is performed with :meth:`PSAnalysis.run` and stores the result @@ -612,20 +753,20 @@ class Path(object): alignment of the original trajectories to a reference structure. .. versionadded:: 0.9.1 + """ def __init__(self, universe, reference, ref_select='name CA', path_select='all', ref_frame=0): """Setting up trajectory alignment and fitted path generation. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`MDAnalysis.Universe` object containing a trajectory - *reference* - reference structure; :class:`MDAnalysis.Universe` object; if - ``None`` then *traj* is used (uses the current time step of the - object) [``None``] - *ref_select* + reference : Universe + reference structure (uses `ref_frame` from the trajectory) + ref_select : str or dict or tuple (optional) The selection to operate on for rms fitting; can be one of: 1. any valid selection string for @@ -641,12 +782,12 @@ def __init__(self, universe, reference, ref_select='name CA', can also each be a list of selection strings (to generate an AtomGroup with defined atom order as described under :ref:`ordered-selections-label`). - *path_select* + ref_frame : int + frame index to select the coordinate frame from + `ref_select.trajectory` + path_select : selection_string atom selection composing coordinates of (fitted) path; if ``None`` - then *path_select* is set to *ref_select* [``None``] - - .. _ClustalW: http://www.clustal.org/ - .. _STAMP: http://www.compbio.dundee.ac.uk/manuals/stamp.4.2/ + then `path_select` is set to `ref_select` [``None``] """ self.u_original = universe @@ -666,26 +807,32 @@ def __init__(self, universe, reference, ref_select='name CA', def fit_to_reference(self, filename=None, prefix='', postfix='_fit', rmsdfile=None, targetdir=os.path.curdir, mass_weighted=False, tol_mass=0.1): - """Align each trajectory frame to the reference structure with - :func:`MDAnalysis.analysis.align.AlignTraj`. + """Align each trajectory frame to the reference structure - :Arguments: - *filename* + Parameters + ---------- + filename : str (optional) file name for the RMS-fitted trajectory or pdb; defaults to the original trajectory filename (from :attr:`Path.u_original`) with *prefix* prepended - *prefix* + prefix : str (optional) prefix for auto-generating the new output filename - *rmsdfile* + rmsdfile : str (optional) file name for writing the RMSD time series [``None``] - *mass_weighted* - do a mass-weighted RMSD fit - *tol_mass* + mass_weighted : bool (optional) + do a mass-weighted RMSD fit, default is ``False`` + tol_mass : float (optional) Reject match if the atomic masses for matched atoms differ by more - than *tol_mass* [0.1] + than `tol_mass` [0.1] - :Returns: - :class:`MDAnalysis.Universe` object containing a fitted trajectory + Returns + ------- + Universe + :class:`MDAnalysis.Universe` object containing a fitted trajectory + + Notes + ----- + Uses :class:`MDAnalysis.analysis.align.AlignTraj` for the fitting. """ head, tail = os.path.split(self.trj_name) oldname, ext = os.path.splitext(tail) @@ -714,24 +861,28 @@ def to_path(self, fitted=False, select=None, flat=False): (if *flat* is ``True``) :class:`numpy.ndarray` representation of the fitted trajectory (with dimensions |3D| or |2D|, respectively). - :Arguments: - *fitted* + Parameters + ---------- + fitted : bool (optional) construct a :attr:`Path.path` from the :attr:`Path.u_fitted` trajectory; if ``False`` then :attr:`Path.path` is generated with the trajectory from :attr:`Path.u_original` [``False``] - *select* + select : str (optional) the selection for constructing the coordinates of each frame in :attr:`Path.path`; if ``None`` then :attr:`Path.path_select` is used, else it is overridden by *select* [``None``] - *flat* + flat : bool (optional) represent :attr:`Path.path` as a 2D (|2D|) :class:`numpy.ndarray`; if ``False`` then :attr:`Path.path` is a 3D (|3D|) :class:`numpy.ndarray` [``False``] - :Returns: - :class:`numpy.ndarray` representing a time series of atomic positions - of an :class:`MDAnalysis.core.groups.AtomGroup` selection from - :attr:`Path.u_fitted.trajectory` + Returns + ------- + numpy.ndarray + representing a time series of atomic positions of an + :class:`MDAnalysis.core.groups.AtomGroup` selection from + :attr:`Path.u_fitted.trajectory` + """ select = select if select is not None else self.path_select if fitted: @@ -756,8 +907,10 @@ def to_path(self, fitted=False, select=None, flat=False): def run(self, align=False, filename=None, postfix='_fit', rmsdfile=None, targetdir=os.path.curdir, mass_weighted=False, tol_mass=0.1, flat=False): - r"""Generate a path from a trajectory and reference structure, aligning - to a reference structure if specified. + r"""Generate a path from a trajectory and reference structure. + + As part of the path generation, the trajectory can be superimposed + ("aligned") to a reference structure if specified. This is a convenience method to generate a fitted trajectory from an inputted universe (:attr:`Path.u_original`) and reference structure @@ -771,32 +924,35 @@ def run(self, align=False, filename=None, postfix='_fit', rmsdfile=None, ``*`` operator, as in ``MDAnalysis.Universe(*(top_name, newtraj_name))``. - :Arguments: - *align* + Parameters + ---------- + align : bool (optional) Align trajectory to atom selection :attr:`Path.ref_select` of :attr:`Path.u_reference`. If ``True``, a universe containing an aligned trajectory is produced with :meth:`Path.fit_to_reference` [``False``] - *filename* + filename : str (optional) filename for the RMS-fitted trajectory or pdb; defaults to the original trajectory filename (from :attr:`Path.u_original`) with *prefix* prepended - *prefix* + postfix : str (optional) prefix for auto-generating the new output filename - *rmsdfile* + rmsdfile : str (optional) file name for writing the RMSD time series [``None``] - *mass_weighted* + mass_weighted : bool (optional) do a mass-weighted RMSD fit - *tol_mass* + tol_mass : float (optional) Reject match if the atomic masses for matched atoms differ by more than *tol_mass* [0.1] - *flat* + flat : bool (optional) represent :attr:`Path.path` with 2D (|2D|) :class:`numpy.ndarray`; if ``False`` then :attr:`Path.path` is a 3D (|3D|) :class:`numpy.ndarray` [``False``] - :Returns: - A tuple of the topology name and new trajectory name. + Returns + ------- + topology_trajectory : tuple + A tuple of the topology name and new trajectory name. """ if align: self.u_fitted = self.fit_to_reference( \ @@ -812,13 +968,15 @@ def get_num_atoms(self): Must run :meth:`Path.to_path` prior to calling this method. - :Returns: - ``int``, the number of atoms in the :class:`Path` + Returns + ------- + int + the number of atoms in the :class:`Path` + """ if self.natoms is None: - err_str = "No path data; do 'Path.to_path()' first." - raise ValueError(err_str) + raise ValueError("No path data; do 'Path.to_path()' first.") return self.natoms @@ -869,14 +1027,15 @@ def __init__(self, npaths, i, j): distances are actually computed), or a single ID (index) in the corresponding distance vector. - :Arguments: - *npaths* - int, total number of paths in :class:`PSA` used to generate *this* - :class:`PSAPair` - *i* - int, row index (starting at 0) of the distance matrix - *j* - int, column index (starting at 0) of the distance matrix + Parameters + ---------- + npaths : int + total number of paths in :class:`PSA` used to generate *this* + :class:`PSAPair` + i : int + row index (starting at 0) of the distance matrix + j : int + column index (starting at 0) of the distance matrix """ self.npaths = npaths self.matrix_idx = (i,j) @@ -898,14 +1057,17 @@ def _dvec_idx(self, i, j): index *j* should be greater than *i+1*, corresponding to the upper triangle of the distance matrix. - :Arguments: - *i* - int, row index (starting at 0) of the distance matrix - *j* - int, column index (starting at 0) of the distance matrix - - :Returns: - int, (matrix element) index in the corresponding distance vector + Parameters + ---------- + i : int + row index (starting at 0) of the distance matrix + j : int + column index (starting at 0) of the distance matrix + + Returns + ------- + int + (matrix element) index in the corresponding distance vector """ return (self.npaths*i) + j - (i+2)*(i+1)/2 @@ -923,13 +1085,15 @@ def compute_nearest_neighbors(self, P,Q, N=None): :math:`i^\text{th}` path and *Q* is the :math:`j^\text{th}` path among the set of *N* total paths in the comparison. - :Arguments: - *P* - :class:`numpy.ndarray` representing a path - *Q* - :class:`numpy.ndarray` representing a path - *N* - int, size of the distance matrix (of shape *N*-by-*N*) [``None``] + Parameters + ---------- + P : numpy.ndarray + representing a path + Q : numpy.ndarray + representing a path + N : int + size of the distance matrix (of shape *N*-by-*N*) [``None``] + """ hn = hausdorff_neighbors(P, Q) self.nearest_neighbors['frames'] = hn['frames'] @@ -974,21 +1138,25 @@ def get_nearest_neighbors(self, frames=True, distances=True): calling :meth:`compute_nearest_neighbors`. At least one of *frames* or *distances* must be ``True``, or else a ``NoDataError`` is raised. - :Arguments: - *frames* - boolean; if ``True``, return nearest neighbor frame indices + Parameters + ---------- + frames : bool + if ``True``, return nearest neighbor frame indices [``True``] - *distances* - boolean; if ``True``, return nearest neighbor distances [``True``] + distances : bool + if ``True``, return nearest neighbor distances [``True``] + + Returns + ------- + dict or tuple + If both *frames* and *distances* are ``True``, return the entire + dictionary (:attr:`nearest_neighbors`); if only *frames* is + ``True``, return a pair of :class:`numpy.ndarray` containing the + indices of the frames (for the pair of paths) of the nearest + neighbors; if only *distances* is ``True``, return a pair of + :class:`numpy.ndarray` of the nearest neighbor distances (for the + pair of paths). - :Returns: - If both *frames* and *distances* are ``True``, return the entire - dictionary (:attr:`nearest_neighbors`); if only *frames* is - ``True``, return a pair of :class:`numpy.ndarray` containing the - indices of the frames (for the pair of paths) of the nearest - neighbors; if only *distances* is ``True``, return a pair of - :class:`numpy.ndarray` of the nearest neighbor distances (for the pair - of paths). """ if self.nearest_neighbors['distances'] is None: err_str = "Nearest neighbors have not been calculated yet;" \ @@ -1018,19 +1186,22 @@ def get_hausdorff_pair(self, frames=True, distance=True): :meth:`find_hausdorff_pair`. At least one of *frames* or *distance* must be ``True``, or else a ``NoDataError`` is raised. - :Arguments: - *frames* - boolean; if ``True``, return the indices of the frames + Parameters + ---------- + frames : bool + if ``True``, return the indices of the frames of the Hausdorff pair [``True``] - *distances* - boolean; if ``True``, return Hausdorff distance [``True``] - - :Returns: - If both *frames* and *distance* are ``True``, return the entire - dictionary (:attr:`hausdorff_pair`); if only *frames* is - ``True``, return a pair of ``int`` containing the indices of the - frames (one index per path) of the Hausdorff pair; if only *distance* - is ``True``, return the Hausdorff distance for this path pair. + distances : bool + if ``True``, return Hausdorff distance [``True``] + + Returns + ------- + dict or tuple + If both *frames* and *distance* are ``True``, return the entire + dictionary (:attr:`hausdorff_pair`); if only *frames* is + ``True``, return a pair of ``int`` containing the indices of the + frames (one index per path) of the Hausdorff pair; if only *distance* + is ``True``, return the Hausdorff distance for this path pair. """ if self.hausdorff_pair['distance'] is None: err_str = "Hausdorff pair has not been calculated yet;" \ @@ -1069,15 +1240,16 @@ def __init__(self, universes, reference=None, ref_select='name CA', The mutual similarity between all unique pairs of trajectories are computed using a selected path metric. - :Arguments: - *universes* + Parameters + ---------- + universes : list a list of universes (:class:`MDAnalysis.Universe` object), each containing a trajectory - *reference* + reference : Universe reference coordinates; :class:`MDAnalysis.Universe` object; if - ``None`` the first time step of the first item in *trajs* is used + ``None`` the first time step of the first item in `universes` is used [``None``] - *ref_select* + ref_select : str or dict or tuple The selection to operate on; can be one of: 1. any valid selection string for @@ -1093,22 +1265,28 @@ def __init__(self, universes, reference=None, ref_select='name CA', can also each be a list of selection strings (to generate an AtomGroup with defined atom order as described under :ref:`ordered-selections-label`). - *mass_weighted* + mass_weighted : bool do a mass-weighted RMSD fit [``False``] - *tol_mass* + tol_mass : float Reject match if the atomic masses for matched atoms differ by more than *tol_mass* [0.1] - *ref_frame* + ref_frame : int frame index to select frame from *reference* [0] - *path_select* + path_select : str atom selection composing coordinates of (fitted) path; if ``None`` then *path_select* is set to *ref_select* [``None``] - *targetdir* - output files are saved there [.] - *labels* + targetdir : str + output files are saved there; if ``None`` then "./psadata" is + created and used [.] + labels : list list of strings, names of trajectories to be analyzed (:class:`MDAnalysis.Universe`); if ``None``, defaults to trajectory names [``None``] + + + .. _ClustalW: http://www.clustal.org/ + .. _STAMP: http://www.compbio.dundee.ac.uk/manuals/stamp.4.2/ + """ self.universes = universes self.u_reference = self.universes[0] if reference is None else reference @@ -1177,35 +1355,37 @@ def __init__(self, universes, reference=None, ref_select='name CA', def generate_paths(self, **kwargs): """Generate paths, aligning each to reference structure if necessary. - :Keywords: - *align* + Parameters + ---------- + align : bool Align trajectories to atom selection :attr:`PSAnalysis.ref_select` of :attr:`PSAnalysis.u_reference` [``False``] - *filename* + filename : str strings representing base filename for fitted trajectories and paths [``None``] - *infix* + infix : str additional tag string that is inserted into the output filename of the fitted trajectory files [''] - *mass_weighted* + mass_weighted : bool do a mass-weighted RMSD fit - *tol_mass* + tol_mass : float Reject match if the atomic masses for matched atoms differ by more than *tol_mass* - *ref_frame* + ref_frame : int frame index to select frame from *reference* - *flat* + flat : bool represent :attr:`Path.path` as a 2D (|2D|) :class:`numpy.ndarray`; if ``False`` then :attr:`Path.path` is a 3D (|3D|) :class:`numpy.ndarray` [``False``] - *save* - boolean; if ``True``, pickle list of names for fitted trajectories + save : bool + if ``True``, pickle list of names for fitted trajectories [``True``] - *store* - boolean; if ``True`` then writes each path (:class:`numpy.ndarray`) + store : bool + if ``True`` then writes each path (:class:`numpy.ndarray`) in :attr:`PSAnalysis.paths` to compressed npz (numpy) files [``False``] + The fitted trajectories are written to new files in the "/trj_fit" subdirectory in :attr:`PSAnalysis.targetdir` named "filename(*trajectory*)XXX*infix*_psa", where "XXX" is a number between @@ -1214,6 +1394,7 @@ def generate_paths(self, **kwargs): format in the "/paths" subdirectory in :attr:`PSAnalysis.targetdir` for persistence and can be accessed as the attribute :attr:`PSAnalysis.paths`. + """ align = kwargs.pop('align', False) filename = kwargs.pop('filename', 'fitted') @@ -1258,18 +1439,24 @@ def run(self, **kwargs): A number of parameters can be changed from the defaults. The result is stored as the array :attr:`PSAnalysis.D`. - :Keywords: - *metric* + Parameters + ---------- + metric : str or callable selection string specifying the path metric to measure pairwise - distances among :attr:`PSAnalysis.paths` [``'hausdorff'``] - *start*, *stop*, *step* - start and stop frame index with step size: analyze + distances among :attr:`PSAnalysis.paths` or a callable with the + same call signature as :func:`hausdorff` + [``'hausdorff'``] + start : int + `start` and `stop` frame index with `step` size: analyze ``trajectory[start:stop:step]`` [``None``] - *store* - boolean; if ``True`` then writes :attr:`PSAnalysis.D` to text and + stop : int + step : int + store : bool + if ``True`` then writes :attr:`PSAnalysis.D` to text and compressed npz (numpy) files [``True``] - *filename* + filename : str string, filename to save :attr:`PSAnalysis.D` + """ metric = kwargs.pop('metric', 'hausdorff') start = kwargs.pop('start', None) @@ -1310,15 +1497,18 @@ def run_pairs_analysis(self, **kwargs): of Hausdorff pairs analysis that is available for each pair of path, including nearest neighbors lists and the Hausdorff pairs. - :Keywords: - *start*, *stop*, *step* - start and stop frame index with step size: analyze + Parameters + ---------- + start : int + `start` and `stop` frame index with `step` size: analyze ``trajectory[start:stop:step]`` [``None``] - *neighbors* - boolean; if ``True``, then stores dictionary of nearest neighbor + stop : int + step : int + neighbors : bool + if ``True``, then stores dictionary of nearest neighbor frames/distances in :attr:`PSAnalysis.NN` [``False``] - *hausdorff_pairs* - boolean; if ``True``, then stores dictionary of Hausdorff pair + hausdorff_pairs : bool + if ``True``, then stores dictionary of Hausdorff pair frames/distances in :attr:`PSAnalysis.HP` [``False``] """ start = kwargs.pop('start', None) @@ -1350,13 +1540,19 @@ def save_result(self, filename=None): """Save distance matrix :attr:`PSAnalysis.D` to a numpy compressed npz file and text file. - :Arguments: - *filename* - string, specifies filename [``None``] - The data are saved with :func:`numpy.savez_compressed` and :func:`numpy.savetxt` in the directory specified by :attr:`PSAnalysis.targetdir`. + + Parameters + ---------- + filename : str + specifies filename [``None``] + + Returns + ------- + filename : str + """ filename = filename or 'psa_distances' head = self.targetdir + self.datadirs['distance_matrices'] @@ -1373,12 +1569,18 @@ def save_result(self, filename=None): def save_paths(self, filename=None): """Save fitted :attr:`PSAnalysis.paths` to numpy compressed npz files. - :Arguments: - *filename* - string, specifies filename [``None``] - The data are saved with :func:`numpy.savez_compressed` in the directory specified by :attr:`PSAnalysis.targetdir`. + + Parameters + ---------- + filename : str + specifies filename [``None``] + + Returns + ------- + filename : str + """ filename = filename or 'path_psa' head = self.targetdir + self.datadirs['paths'] @@ -1413,28 +1615,33 @@ def load(self): def plot(self, filename=None, linkage='ward', count_sort=False, distance_sort=False, figsize=4.5, labelsize=12): - """Plot a clustered distance matrix using method *linkage* along with - the corresponding dendrogram. Rows (and columns) are identified using - the list of strings specified by :attr:`PSAnalysis.labels`. + """Plot a clustered distance matrix. - :Arguments: - *filename* - string, save figure to *filename* [``None``] - *linkage* - string, name of linkage criterion for clustering [``'ward'``] - *count_sort* - boolean, see scipy.cluster.hierarchy.dendrogram [``False``] - *distance_sort* - boolean, see scipy.cluster.hierarchy.dendrogram [``False``] - *figsize* + Usese method *linkage* and plots the corresponding dendrogram. Rows + (and columns) are identified using the list of strings specified by + :attr:`PSAnalysis.labels`. + + If `filename` is supplied then the figure is also written to file (the + suffix determines the file type, e.g. pdf, png, eps, ...). All other + keyword arguments are passed on to :func:`matplotlib.pyplot.imshow`. + + + Parameters + ---------- + filename : str + save figure to *filename* [``None``] + linkage : str + name of linkage criterion for clustering [``'ward'``] + count_sort : bool + see :func:`scipy.cluster.hierarchy.dendrogram` [``False``] + distance_sort : bool + see :func:`scipy.cluster.hierarchy.dendrogram` [``False``] + figsize : float set the vertical size of plot in inches [``4.5``] - *labelsize* + labelsize : float set the font size for colorbar labels; font size for path labels on dendrogram default to 3 points smaller [``12``] - If *filename* is supplied then the figure is also written to file (the - suffix determines the file type, e.g. pdf, png, eps, ...). All other - keyword arguments are passed on to :func:`pylab.imshow`. """ from matplotlib.pyplot import figure, colorbar, cm, savefig, clf @@ -1511,29 +1718,33 @@ def plot(self, filename=None, linkage='ward', count_sort=False, def plot_annotated_heatmap(self, filename=None, linkage='ward', \ count_sort=False, distance_sort=False, \ figsize=8, annot_size=6.5): - """Plot a clustered distance matrix using method *linkage* with - annotated distances in the matrix. Rows (and columns) are identified - using the list of strings specified by :attr:`PSAnalysis.labels`. + """Plot a clustered distance matrix. - :Arguments: - *filename* - string, save figure to *filename* [``None``] - *count_sort* - boolean, see scipy.cluster.hierarchy.dendrogram [``False``] - *distance_sort* - boolean, see scipy.cluster.hierarchy.dendrogram [``False``] - *linkage* - string, name of linkage criterion for clustering [``'ward'``] - *figsize* - set the vertical size of plot in inches [``8``] - *annot_size* - float, font size of annotation labels on heat map [``6.5``] - - If *filename* is supplied then the figure is also written to file (the + Uses method `linkage` and plots annotated distances in the matrix. Rows + (and columns) are identified using the list of strings specified by + :attr:`PSAnalysis.labels`. + + If `filename` is supplied then the figure is also written to file (the suffix determines the file type, e.g. pdf, png, eps, ...). All other - keyword arguments are passed on to :func:`pylab.imshow`. + keyword arguments are passed on to :func:`matplotlib.pyplot.imshow`. + + Parameters + ---------- + filename : str + save figure to *filename* [``None``] + linkage : str + name of linkage criterion for clustering [``'ward'``] + count_sort : bool + see :func:`scipy.cluster.hierarchy.dendrogram` [``False``] + distance_sort : bool + see :func:`scipy.cluster.hierarchy.dendrogram` [``False``] + figsize : float + set the vertical size of plot in inches [``4.5``] + annot_size : float + font size of annotation labels on heat map [``6.5``] + """ - from matplotlib.pylab import figure, colorbar, cm, savefig, clf + from matplotlib.pyplot import figure, colorbar, cm, savefig, clf try: import seaborn.apionly as sns @@ -1604,30 +1815,34 @@ def plot_nearest_neighbors(self, filename=None, idx=0, \ multiplot=False, aspect_ratio=1.75, \ labelsize=12): """Plot nearest neighbor distances as a function of normalized frame - number (mapped to the interval *[0,1]*). + number. - :Arguments: - *filename* - string, save figure to *filename* [``None``] - *idx* - integer, index of path (pair) comparison to plot [``0``] - *labels* - (string, string), pair of names to label nearest neighbor distance + The frame number is mapped to the interval *[0, 1]*. + + If `filename` is supplied then the figure is also written to file (the + suffix determines the file type, e.g. pdf, png, eps, ...). All other + keyword arguments are passed on to :func:`matplotlib.pyplot.imshow`. + + Parameters + ---------- + filename : str + save figure to *filename* [``None``] + idx : int + index of path (pair) comparison to plot [``0``] + labels : (str, str) + pair of names to label nearest neighbor distance curves [``('Path 1', 'Path 2')``] - *figsize* - float, set the vertical size of plot in inches [``4.5``] - *multiplot* - boolean, set to ``True`` to enable plotting multiple nearest + figsize : float + set the vertical size of plot in inches [``4.5``] + multiplot : bool + set to ``True`` to enable plotting multiple nearest neighbor distances on the same figure [``False``] - *aspect_ratio* - float, set the ratio of width to height of the plot [``1.75``] - *labelsize* + aspect_ratio : float + set the ratio of width to height of the plot [``1.75``] + labelsize : float set the font size for colorbar labels; font size for path labels on dendrogram default to 3 points smaller [``12``] - If *filename* is supplied then the figure is also written to file (the - suffix determines the file type, e.g. pdf, png, eps, ...). All other - keyword arguments are passed on to :func:`pylab.imshow`. """ from matplotlib.pyplot import figure, savefig, tight_layout, clf, show try: @@ -1693,14 +1908,15 @@ def cluster(self, distArray, method='ward', count_sort=False, \ color_threshold=4): """Cluster trajectories and optionally plot the dendrogram. - :Arguments: - *method* - string, name of linkage criterion for clustering [``'ward'``] - *no_plot* - boolean, if ``True``, do not render the dendrogram [``False``] - *no_labels* - boolean, if ``True`` then do not label dendrogram [``True``] - *color_threshold* + Parameters + ---------- + method : str + name of linkage criterion for clustering [``'ward'``] + no_plot : bool + if ``True``, do not render the dendrogram [``False``] + no_labels : bool + if ``True`` then do not label dendrogram [``True``] + color_threshold : float For brevity, let t be the color_threshold. Colors all the descendent links below a cluster node k the same color if k is the first node below the cut threshold t. All links connecting @@ -1709,9 +1925,12 @@ def cluster(self, distArray, method='ward', count_sort=False, \ colored blue. If color_threshold is None or ‘default’, corresponding with MATLAB(TM) behavior, the threshold is set to 0.7*max(Z[:,2]). [``4``]] - :Returns: - list of indices representing the row-wise order of the objects - after clustering + + Returns + ------- + list + list of indices representing the row-wise order of the objects + after clustering """ import matplotlib from scipy.cluster.hierarchy import linkage, dendrogram @@ -1728,7 +1947,9 @@ def cluster(self, distArray, method='ward', count_sort=False, \ def _get_plot_obj_locs(self): """Find and return coordinates for dendrogram, heat map, and colorbar. - :Returns: + Returns + ------- + tuple tuple of coordinates for placing the dendrogram, heat map, and colorbar in the plot. """ @@ -1759,8 +1980,10 @@ def get_num_atoms(self): Must run :meth:`PSAnalysis.generate_paths` prior to calling this method. - :Returns: - int, the number of atoms in :class:`PSA`'s :class:`Path`s' + Returns + ------- + int + the number of atoms in :class:`PSA`'s :class:`Path`s' """ if self.natoms is None: err_str = "No path data; do 'PSAnalysis.generate_paths()' first." @@ -1775,8 +1998,10 @@ def get_num_paths(self): Must run :meth:`PSAnalysis.generate_paths` prior to calling this method. - :Returns: - int, the number of paths in :class:`PSA` + Returns + ------- + int + the number of paths in :class:`PSA` """ if self.npaths is None: err_str = "No path data; do 'PSAnalysis.generate_paths()' first." @@ -1791,9 +2016,11 @@ def get_paths(self): Must run :meth:`PSAnalysis.generate_paths` prior to calling this method. - :Returns: - list of :class:`numpy.ndarray` representations of paths in - :class:`PSA` + Returns + ------- + list + list of :class:`numpy.ndarray` representations of paths in + :class:`PSA` """ if self.paths is None: err_str = "No path data; do 'PSAnalysis.generate_paths()' first." @@ -1808,13 +2035,15 @@ def get_pairwise_distances(self, vectorform=False): Must run :meth:`PSAnalysis.run` with ``store=True`` prior to calling this method. - :Arguments: - *vectorform* - boolean, if ``True``, return the distance vector instead [``False``] + Parameters + ---------- + vectorform : bool + if ``True``, return the distance vector instead [``False``] - :Returns: - :class:`numpy.ndarray` representation of the distance matrix (or - vector) + Returns + ------- + numpy.ndarray + representation of the distance matrix (or vector) """ if self.D is None: @@ -1829,15 +2058,16 @@ def get_pairwise_distances(self, vectorform=False): @property def psa_pairs(self): - """Get the list of :class:`PSAPair` instances for each pair of paths. - - :meth:`psa_pairs` is a list of :class:`PSAPair` whose elements are - pairs of paths that have been compared using - :meth:`PSAnalysis.run_pairs_analysis`. Each :class:`PSAPair` - contains nearest neighbor and Hausdorff pair information specific to a - pair of paths. The nearest neighbor frames and distances for a - :class:`PSAPair` can be accessed in the nearest neighbor dictionary - using the keys 'frames' and 'distances', respectively. E.g., + """The list of :class:`PSAPair` instances for each pair of paths. + + :attr:`psa_pairs` is a list of all :class:`PSAPair` objects (in + distance vector order). The elements of a :class:`PSAPair` are pairs of + paths that have been compared using + :meth:`PSAnalysis.run_pairs_analysis`. Each :class:`PSAPair` contains + nearest neighbor and Hausdorff pair information specific to a pair of + paths. The nearest neighbor frames and distances for a :class:`PSAPair` + can be accessed in the nearest neighbor dictionary using the keys + 'frames' and 'distances', respectively. E.g., :attr:`PSAPair.nearest_neighbors['distances']` returns a *pair* of :class:`numpy.ndarray` corresponding to the nearest neighbor distances for each path. Similarly, Hausdorff pair information can be accessed @@ -1848,9 +2078,6 @@ def psa_pairs(self): Must run :meth:`PSAnalysis.run_pairs_analysis` prior to calling this method. - :Returns: - list of all :class:`PSAPair` objects (in distance vector order) - """ if self._psa_pairs is None: err_str = "No nearest neighbors data; do" \ @@ -1861,21 +2088,19 @@ def psa_pairs(self): @property def hausdorff_pairs(self): - """Get the Hausdorff pair for each (unique) pairs of paths. + """The Hausdorff pair for each (unique) pairs of paths. - This method returns a list of Hausdorff pair information, where each - element is a dictionary containing the pair of frames and the - (Hausdorff) distance between a pair of paths. See - :meth:`PSAnalysis.psa_pairs` and :attr:`PSAPair.hausdorff_pair` for - more information about accessing Hausdorff pair data. + This attribute contains a list of Hausdorff pair information (in + distance vector order), where each element is a dictionary containing + the pair of frames and the (Hausdorff) distance between a pair of + paths. See :meth:`PSAnalysis.psa_pairs` and + :attr:`PSAPair.hausdorff_pair` for more information about accessing + Hausdorff pair data. .. note:: Must run :meth:`PSAnalysis.run_pairs_analysis` with ``hausdorff_pairs=True`` prior to calling this method. - :Returns: - list of all Hausdorff pairs (in distance vector order) - """ if self._HP is None: err_str = "No Hausdorff pairs data; do " \ @@ -1886,21 +2111,18 @@ def hausdorff_pairs(self): @property def nearest_neighbors(self): - """Get the nearest neighbors for each (unique) pair of paths. + """The nearest neighbors for each (unique) pair of paths. - This method returns a list of nearest neighbor information, where each - element is a dictionary containing the nearest neighbor frames and - distances between a pair of paths. See :meth:`PSAnalysis.psa_pairs` - and :attr:`PSAPair.nearest_neighbors` for more information about - accessing nearest neighbor data. + This attribute contains a list of nearest neighbor information (in + distance vector order), where each element is a dictionary containing + the nearest neighbor frames and distances between a pair of paths. See + :meth:`PSAnalysis.psa_pairs` and :attr:`PSAPair.nearest_neighbors` for + more information about accessing nearest neighbor data. .. note:: Must run :meth:`PSAnalysis.run_pairs_analysis` with ``neighbors=True`` prior to calling this method. - :Returns: - list of all nearest neighbors (in distance vector order) - """ if self._NN is None: err_str = "No nearest neighbors data; do" \ From 0bfbc9c37ae3013bd056c336643f421232565b50 Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Thu, 23 Mar 2017 03:21:10 -0700 Subject: [PATCH 3/9] analysis: numpyfying/fixing docs to numpy style in multiple modules - analysis.align: formatting fixes - analysis.contacts: formatting fixes - analysis.diffusionmap: formatting fixes and section headers (cannot use 'Examples' as a normal section header because it is rewritten by sphinx.ext.napoleon as a rubric) - analysis.distances: numpyfied - analysis.hbonds.hbond_analysis: numpyfied --- package/MDAnalysis/analysis/align.py | 103 +++++++++--------- package/MDAnalysis/analysis/contacts.py | 25 ++--- package/MDAnalysis/analysis/density.py | 28 ++--- package/MDAnalysis/analysis/diffusionmap.py | 44 ++++---- package/MDAnalysis/analysis/distances.py | 14 ++- .../analysis/hbonds/hbond_analysis.py | 50 +++++---- 6 files changed, 136 insertions(+), 128 deletions(-) diff --git a/package/MDAnalysis/analysis/align.py b/package/MDAnalysis/analysis/align.py index debf25ff47a..dcb8f0e74e7 100644 --- a/package/MDAnalysis/analysis/align.py +++ b/package/MDAnalysis/analysis/align.py @@ -174,6 +174,7 @@ .. autofunction:: _fit_to .. autofunction:: fasta2select +.. autofunction:: sequence_alignment .. autofunction:: get_matching_atoms """ @@ -281,9 +282,9 @@ def _fit_to(mobile_coordinates, ref_coordinates, mobile_atoms, Parameters ---------- - mobile_coordinates : array + mobile_coordinates : ndarray Coordinates of atoms to be aligned - ref_coordinates : array + ref_coordinates : ndarray Coordinates of atoms to be fit against mobile_atoms : AtomGroup Atoms to be translated @@ -291,7 +292,7 @@ def _fit_to(mobile_coordinates, ref_coordinates, mobile_atoms, array of xyz coordinate of mobile center of mass ref_com: ndarray array of xyz coordinate of reference center of mass - weights : numpy array, optional + weights : numpy array (optional) Array to be used for weighted rmsd Returns @@ -352,7 +353,7 @@ def alignto(mobile, reference, select="all", mass_weighted=None, weights=None, reference : Universe or AtomGroup reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole :class:`~MDAnalysis.core.universe.Universe` - select: string or dict, optional + select: string or dict (optional) 1. any valid selection string for :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` that produces identical selections in `mobile` and `reference`; or @@ -364,15 +365,15 @@ def alignto(mobile, reference, select="all", mass_weighted=None, weights=None, When using 2. or 3. with *sel1* and *sel2* then these selections can also each be a list of selection strings (to generate a AtomGroup with defined atom order as described under :ref:`ordered-selections-label`). - mass_weighted : boolean, optional (deprecated) + mass_weighted : boolean (optional) (deprecated) ``True`` uses the masses :meth:`reference.masses` as weights for the RMSD fit. - weights : str/array_like, optional + weights : str/array_like (optional) weights to be used for fit. Can be either 'mass' or an array_like - tol_mass: float, optional + tol_mass: float (optional) Reject match if the atomic masses for matched atoms differ by more than *tol_mass*, default [0.1] - strict: boolean, optional + strict: boolean (optional) ``True`` Will raise :exc:`SelectionError` if a single atom does not match between the two selections. @@ -380,7 +381,7 @@ def alignto(mobile, reference, select="all", mass_weighted=None, weights=None, Will try to prepare a matching selection by dropping residues with non-matching atoms. See :func:`get_matching_atoms` for details. - subselection : string, optional + subselection : string (optional) Apply the transformation only to this selection. ``None`` [default] @@ -399,26 +400,27 @@ def alignto(mobile, reference, select="all", mass_weighted=None, weights=None, new_rmsd : float RMSD after spatial alignment - See Also - -------- - AlignTraj: More efficient method for RMSD-fitting trajectories. + .. SeeAlso:: + AlignTraj: More efficient method for RMSD-fitting trajectories. + .. _ClustalW: http://www.clustal.org/ .. _STAMP: http://www.compbio.dundee.ac.uk/manuals/stamp.4.2/ + .. versionchanged:: 0.8 Added check that the two groups describe the same atoms including the new *tol_mass* keyword. .. versionchanged:: 0.10.0 Uses :func:`get_matching_atoms` to work with incomplete selections - and new *strict* keyword. The new default is to be lenient whereas - the old behavior was the equivalent of *strict* = ``True``. + and new `strict` keyword. The new default is to be lenient whereas + the old behavior was the equivalent of ``strict = True``. .. versionchanged:: 0.16.0 new general 'weights' kwarg replace mass_weights, deprecated 'mass_weights' .. deprecated:: 0.16.0 - Instead of ``mass_weighted=True`` use new ``weights='mass'` + Instead of ``mass_weighted=True`` use new ``weights='mass'`` """ if select in ('all', None): @@ -510,36 +512,37 @@ def __init__(self, mobile, reference, select='all', filename=None, Universe containing trajectory to be fitted to reference reference : Universe Universe containing trajectory frame to be used as reference - select : string, optional + select : string (optional) Set as default to all, is used for Universe.select_atoms to choose subdomain to be fitted against - filename : string, optional + filename : string (optional) Provide a filename for results to be written to - prefix : string, optional + prefix : string (optional) Provide a string to prepend to filename for results to be written to - mass_weighted : boolean, optional (deprecated) - Boolean, if true will rmsd will be mass-weighted corresponding to - the atoms selected in the reference trajectory + mass_weighted : boolean (optional, deprecated) + Boolean, if ``True`` will rmsd will be mass-weighted corresponding + to the atoms selected in the reference trajectory. However, this is + deprecated: rather use ``weights="mass"``. weights : str/array_like (optional) Can either be 'mass' to use reference masses as weights or an arbitrary array - tol_mass : float, optional + tol_mass : float (optional) Tolerance given to `get_matching_atoms` to find appropriate atoms - strict : boolean, optional + strict : boolean (optional) Force `get_matching_atoms` to fail if atoms can't be found using exact methods - force : boolean, optional + force : boolean (optional) Force overwrite of filename for rmsd-fitting - verbose : boolean, optional + verbose : boolean (optional) Set logger to show more information - start : int, optional + start : int (optional) First frame of trajectory to analyse, Default: 0 - stop : int, optional + stop : int (optional) Last frame of trajectory to analyse, Default: -1 - step : int, optional + step : int (optional) Step between frames to analyse, Default: 1 - in_memory : boolean, optional + in_memory : boolean (optional) *Permanently* switch `mobile` to an in-memory trajectory so that alignment can be done in-place, which can improve performance substantially in some cases. In this case, no file @@ -562,6 +565,7 @@ def __init__(self, mobile, reference, select='all', filename=None, new general 'weights' kwarg replace mass_weights, deprecated 'mass_weights' .. deprecated:: 0.16.0 Instead of ``mass_weighted=True`` use new ``weights='mass'` + """ select = rms.process_selection(select) self.ref_atoms = reference.select_atoms(*select['reference']) @@ -868,15 +872,15 @@ def sequence_alignment(mobile, reference, match_score=2, mismatch_penalty=-1, Atom group to be aligned reference : AtomGroup Atom group to be aligned against - match_score : float, optional, default 2 + match_score : float (optional), default 2 score for matching residues, default 2 - mismatch_penalty : float, optional, default -1 + mismatch_penalty : float (optional), default -1 penalty for residues that do not match , default : -1 - gap_penalty : float, optional, default -2 + gap_penalty : float (optional), default -2 penalty for opening a gap; the high default value creates compact alignments for highly identical sequences but might not be suitable for sequences with low identity, default : -2 - gapextension_penalty : float, optional, default -0.1 + gapextension_penalty : float (optional), default -0.1 penalty for extending a gap, default: -0.1 Returns @@ -943,31 +947,31 @@ def fasta2select(fastafilename, is_aligned=False, fastafilename : str, path to filename FASTA file with first sequence as reference and second the one to be aligned (ORDER IS IMPORTANT!) - is_aligned : boolean, optional - ``False`` : [default] + is_aligned : boolean (optional) + ``False`` (default) run clustalw for sequence alignment; ``True`` use the alignment in the file (e.g. from STAMP) [``False``] - ref_offset : int, optional + ref_offset : int (optional) add this number to the column number in the FASTA file to get the original residue number, default: 0 - target_offset : int, optional + target_offset : int (optional) add this number to the column number in the FASTA file to get the original residue number, default: 0 - ref_resids : str, optional + ref_resids : str (optional) sequence of resids as they appear in the reference structure - target_resids : str, optional + target_resids : str (optional) sequence of resids as they appear in the target - alnfilename : str, optional + alnfilename : str (optional) filename of ClustalW alignment (clustal format) that is produced by *clustalw* when *is_aligned* = ``False``. default ``None`` uses the name and path of *fastafilename* and subsititutes the suffix with '.aln'. - treefilename: str, optional + treefilename: str (optional) filename of ClustalW guide tree (Newick format); if default ``None`` the the filename is generated from *alnfilename* with the suffix '.dnd' instead of '.aln' - clustalw : str, optional + clustalw : str (optional) path to the ClustalW (or ClustalW2) binary; only needed for `is_aligned` = ``False``, default: "ClustalW2" @@ -975,13 +979,14 @@ def fasta2select(fastafilename, is_aligned=False, ------- select_dict : dict dictionary with 'reference' and 'mobile' selection string - that can be used immediately in :func:`rms_fit_trj` as + that can be used immediately in :class:`AlignTraj` as ``select=select_dict``. - See Also - -------- - :func:`sequence_alignment`, which does not require external - programs. + + .. SeeAlso:: + :func:`sequence_alignment`, which does not require external + programs. + .. _ClustalW: http://www.clustal.org/ .. _STAMP: http://www.compbio.dundee.ac.uk/manuals/stamp.4.2/ @@ -1158,10 +1163,10 @@ def get_matching_atoms(ag1, ag2, tol_mass=0.1, strict=False): ag2 : AtomGroup Second :class:`~MDAnalysis.core.groups.AtomGroup` instance that is compared - tol_mass : float, optional + tol_mass : float (optional) Reject if the atomic masses for matched atoms differ by more than `tol_mass` [0.1] - strict : boolean, optional + strict : boolean (optional) ``True`` Will raise :exc:`SelectionError` if a single atom does not match between the two selections. diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index 1f533e54f15..2722d583a32 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -21,28 +21,27 @@ # """ -================================================================ Native contacts analysis --- :mod:`MDAnalysis.analysis.contacts` ================================================================ - -Analysis of native contacts *Q* over a trajectory. Native contacts of a -conformation are contacts that exist in a reference structure and in the -conformation. Contacts in the reference structure are always defined as being -closer then a distance `radius`. The fraction of native contacts for a -conformation can be calculated in different ways. This module supports 3 -different metrics liseted below, as wel as custom metrics. +This module contains classes to analyze native contacts *Q* over a +trajectory. Native contacts of a conformation are contacts that exist +in a reference structure and in the conformation. Contacts in the +reference structure are always defined as being closer then a distance +`radius`. The fraction of native contacts for a conformation can be +calculated in different ways. This module supports 3 different metrics +listed below, as well as custom metrics. 1. *Hard Cut*: To count as a contact the atoms *i* and *j* have to be at least as close as in the reference structure. 2. *Soft Cut*: The atom pair *i* and *j* is assigned based on a soft potential - that is 1 for if the distance is 0, 1/2 if the distance is the same as in + that is 1 if the distance is 0, 1/2 if the distance is the same as in the reference and 0 for large distances. For the exact definition of the potential and parameters have a look at function :func:`soft_cut_q`. 3. *Radius Cut*: To count as a contact the atoms *i* and *j* cannot be further - apart then some distance `radius`. + apart than some distance `radius`. The "fraction of native contacts" *Q(t)* is a number between 0 and 1 and calculated as the total number of native contacts for a given time frame @@ -139,9 +138,9 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The :class:`Contacts` class has been designed to be extensible for your own -analysis. As an example we will analysis when the acidic and basic groups of -are in contact which each other, this means that at least one of the contacts -formed in the reference is closer then 2.5 Å. +analysis. As an example we will analyze when the acidic and basic groups of AdK +are in contact which each other; this means that at least one of the contacts +formed in the reference is closer than 2.5 Å. For this we define a new function to determine if any contact is closer than 2.5 Å; this function must implement the API prescribed by :class:`Contacts`:: diff --git a/package/MDAnalysis/analysis/density.py b/package/MDAnalysis/analysis/density.py index 76cef9855de..febe09df275 100644 --- a/package/MDAnalysis/analysis/density.py +++ b/package/MDAnalysis/analysis/density.py @@ -382,7 +382,7 @@ def convert_length(self, unit='Angstrom'): Parameters ---------- - unit : str, optional + unit : str (optional) unit that the grid should be converted to: one of "Angstrom", "nm" @@ -406,7 +406,7 @@ def convert_density(self, unit='Angstrom'): Parameters ---------- - unit : str, optional + unit : str (optional) The target unit that the density should be converted to. `unit` can be one of the following: @@ -475,43 +475,43 @@ def density_from_Universe(universe, delta=1.0, atomselection='name OH2', ---------- universe : MDAnalysis.Universe :class:`MDAnalysis.Universe` object with a trajectory - atomselection : str, optional + atomselection : str (optional) selection string (MDAnalysis syntax) for the species to be analyzed ["name OH2"] - delta : float, optional + delta : float (optional) bin size for the density grid in Angstroem (same in x,y,z) [1.0] - start : int, optional - stop : int, optional - step : int, optional + start : int (optional) + stop : int (optional) + step : int (optional) Slice the trajectory as ``trajectory[start:stop:step]``; default is to read the whole trajectory. metadata : dict. optional `dict` of additional data to be saved with the object; the meta data are passed through as they are. - padding : float, optional + padding : float (optional) increase histogram dimensions by padding (on top of initial box size) in Angstroem [2.0] - soluteselection : str, optional + soluteselection : str (optional) MDAnalysis selection for the solute, e.g. "protein" [``None``] - cutoff : float, optional + cutoff : float (optional) With `cutoff`, select " NOT WITHIN OF " (Special routines that are faster than the standard ``AROUND`` selection); any value that evaluates to ``False`` (such as the default 0) disables this special selection. - update_selection : bool, optional + update_selection : bool (optional) Should the selection of atoms be updated for every step? [``False``] - ``True``: atom selection is updated for each frame, can be slow - ``False``: atoms are only selected at the beginning - verbose : bool, optional + verbose : bool (optional) Print status update to the screen for every *interval* frame? [``True``] - ``False``: no status updates when a new frame is processed - ``True``: status update every frame (including number of atoms processed, which is interesting with ``update_selection=True``) - interval : int, optional + interval : int (optional) Show status update every `interval` frame [1] - parameters : dict, optional + parameters : dict (optional) `dict` with some special parameters for :class:`Density` (see docs) Returns diff --git a/package/MDAnalysis/analysis/diffusionmap.py b/package/MDAnalysis/analysis/diffusionmap.py index 27439901982..ae2bc8d1495 100644 --- a/package/MDAnalysis/analysis/diffusionmap.py +++ b/package/MDAnalysis/analysis/diffusionmap.py @@ -20,8 +20,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -""" -Diffusion map --- :mod:`MDAnalysis.analysis.diffusionmap` +"""Diffusion map --- :mod:`MDAnalysis.analysis.diffusionmap` ===================================================================== :Authors: Eugen Hruska, John Detlefs @@ -139,35 +138,30 @@ ---------- If you use this Dimension Reduction method in a publication, please -cite: - -.. [Lafon1] -Coifman, Ronald R., Lafon, Stephane Diffusion maps. Appl. Comput. Harmon. -Anal. 21, 5–30 (2006). +cite [Lafon1]_. -For more information --------------------- +If you choose the default metric, this module uses the fast QCP algorithm +[Theobald2005]_ to calculate the root mean square distance (RMSD) between two +coordinate sets (as implemented in +:func:`MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix`). When using this +module in published work please cite [Theobald2005]_. -.. [deLaPorte1] -J. de la Porte, B. M. Herbst, W. Hereman, S. J. van der Walt. -An Introduction to Diffusion Maps. -.. [Clementi1] -Rohrdanz, M. A, Zheng, W, Maggioni, M, & Clementi, C. -Determination of reaction coordinates via locally scaled -diffusion map. J. Chem. Phys. 134, 124116 (2011). +.. [Lafon1] Coifman, Ronald R., Lafon, Stephane. Diffusion + maps. Appl. Comput. Harmon. Anal. 21, 5–30 (2006). -.. [Ferguson1] -Ferguson, A. L.; Panagiotopoulos, A. Z.; Kevrekidis, I. G. -Debenedetti, P. G. Nonlinear dimensionality reduction in molecular simulation: -The diffusion map approach Chem. Phys. Lett. 509, 1−11 (2011) +.. [deLaPorte1] J. de la Porte, B. M. Herbst, W. Hereman, S. J. van der Walt. + An Introduction to Diffusion Maps. In: The 19th Symposium of the + Pattern Recognition Association of South Africa (2008). -.. If you choose the default metric, this module uses the fast QCP algorithm -[Theobald2005]_ to calculate the root mean square distance (RMSD) between -two coordinate sets (as implemented -in :func:`MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix`). -When using this module in published work please cite [Theobald2005]_. +.. [Clementi1] Rohrdanz, M. A, Zheng, W, Maggioni, M, & Clementi, C. + Determination of reaction coordinates via locally scaled diffusion + map. J. Chem. Phys. 134, 124116 (2011). +.. [Ferguson1] Ferguson, A. L.; Panagiotopoulos, A. Z.; Kevrekidis, I. G. + Debenedetti, P. G. Nonlinear dimensionality reduction in molecular + simulation: The diffusion map approach Chem. Phys. Lett. 509, 1−11 + (2011) """ from six.moves import range diff --git a/package/MDAnalysis/analysis/distances.py b/package/MDAnalysis/analysis/distances.py index 30ba85309a7..fb478ac4257 100644 --- a/package/MDAnalysis/analysis/distances.py +++ b/package/MDAnalysis/analysis/distances.py @@ -95,6 +95,10 @@ def contact_matrix(coord, cutoff=15.0, returntype="numpy", box=None): :mod:`scipy.sparse` is require for using *sparse* matrices; if it cannot be imported then an `ImportError` is raised. + See Also + -------- + :mod:`MDAnalysis.analysis.contacts` for native contact analysis + .. versionchanged:: 0.11.0 Keyword *suppress_progmet* and *progress_meter_freq* were removed. @@ -179,10 +183,11 @@ def between(group, A, B, distance): ---------- group : AtomGroup Find members of `group` that are between `A` and `B` - A, B : AtomGroups - `A` and `B` are :class:`~MDAnalysis.core.groups.AtomGroup` - instances. Works best if `group` is bigger than either `A` or - `B`. + A : AtomGroup + B : AtomGroup + `A` and `B` are the groups of atoms between which atoms in + `group` are searched for. The function works is more + efficient if `group` is bigger than either `A` or `B`. distance : float maximum distance for an atom to be counted as in the vicinity of `A` or `B` @@ -195,6 +200,7 @@ def between(group, A, B, distance): .. versionadded: 0.7.5 + """ ns_group = AtomNeighborSearch(group) resA = set(ns_group.search(A, distance)) diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py index 8976cd2a7cf..45f39b6b75d 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py @@ -214,8 +214,8 @@ class HydrogenBondAnalysis_OtherFF(HydrogenBondAnalysis): .. _`10.1002/prot.340090204`: http://dx.doi.org/10.1002/prot.340090204 -Example -------- +How to perform HydrogenBondAnalysis +----------------------------------- All protein-water hydrogen bonds can be analysed with :: @@ -464,61 +464,62 @@ def __init__(self, universe, selection1='protein', selection2='all', selection1_ you can improve performance by setting the *update_selection* keywords to ``False``. - :Arguments: - *universe* + Parameters + ---------- + universe : universe Universe object - *selection1* + selection1 : str (optional) Selection string for first selection ['protein'] - *selection2* + selection2 : str (optional) Selection string for second selection ['all'] - *selection1_type* + selection1_type : str (optional) Selection 1 can be 'donor', 'acceptor' or 'both'. Note that the value for *selection1_type* automatically determines how *selection2* handles donors and acceptors: If *selection1* contains 'both' then *selection2* will also contain *both*. If *selection1* is set to 'donor' then *selection2* is 'acceptor' (and vice versa). ['both']. - *update_selection1* + update_selection1 : bool (optional) Update selection 1 at each frame? [``False``] - *update_selection2* + update_selection2 : bool (optional) Update selection 2 at each frame? [``False``] - *filter_first* + filter_first : bool (optional) Filter selection 2 first to only atoms 3*distance away [``True``] - *distance* + distance : float (optional) Distance cutoff for hydrogen bonds; only interactions with a H-A distance <= *distance* (and the appropriate D-H-A angle, see *angle*) are recorded. (Note: *distance_type* can change this to the D-A distance.) [3.0] - *angle* + angle : float (optional) Angle cutoff for hydrogen bonds; an ideal H-bond has an angle of 180º. A hydrogen bond is only recorded if the D-H-A angle is >= *angle*. The default of 120º also finds fairly non-specific hydrogen interactions and a possibly better value is 150º. [120.0] - *forcefield* + forcefield : {"CHARMM27", "GLYCAM06", "other"} (optional) Name of the forcefield used. Switches between different :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS` and :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS` values. Available values: "CHARMM27", "GLYCAM06", "other" ["CHARMM27"] - *donors* + donors : sequence (optional) Extra H donor atom types (in addition to those in :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS`), must be a sequence. - *acceptors* + acceptors : sequence (optional) Extra H acceptor atom types (in addition to those in :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS`), must be a sequence. - *start* + start : int (optional) starting frame-index for analysis, ``None`` is the first one, 0. *start* and *stop* are 0-based frame indices and are used to slice the trajectory (if supported) [``None``] - *stop* + stop : int (optional) last trajectory frame for analysis, ``None`` is the last one [``None``] - *step* + step : int (optional) read every *step* between *start* and *stop*, ``None`` selects 1. Note that not all trajectory reader from 1 [``None``] - *debug* + debug : bool (optional) If set to ``True`` enables per-frame debug logging. This is disabled by default because it generates a very large amount of output in the log file. (Note that a logger must have been started to see the output, e.g. using :func:`MDAnalysis.start_logging`.) - *detect_hydrogens* + detect_hydrogens : {"distance", "heuristic"} (optional) Determine the algorithm to find hydrogens connected to donor atoms. Can be "distance" (default; finds all hydrogens in the donor's residue within a cutoff of the donor) or "heuristic" @@ -526,14 +527,17 @@ def __init__(self, universe, selection1='protein', selection2='all', selection1_ always give the correct answer but "heuristic" is faster, especially when the donor list is updated each for each frame. ["distance"] - *distance_type* + distance_type : {"hydrogen", "heavy"} (optional) Measure hydrogen bond lengths between donor and acceptor heavy attoms ("heavy") or between donor hydrogen and acceptor heavy atom ("hydrogen"). If using "heavy" then one should set the *distance* cutoff to a higher value such as 3.5 Å. ["hydrogen"] - :Raises: :exc:`SelectionError` is raised for each static selection without - the required donors and/or acceptors. + Raises + ------ + :exc:`SelectionError` + is raised for each static selection without the required + donors and/or acceptors. .. versionchanged:: 0.7.6 New *verbose* keyword (and per-frame debug logging disabled by From 69b6abe8c453eec062af750c7625fdfdd76bfb7c Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Thu, 30 Mar 2017 02:02:01 -0700 Subject: [PATCH 4/9] analysis.helanal: doc update and code cleanup - numpified docs - removed kwargs start and end for resid selection from def helanal_main() and helanal_trajectory() because this can be easily done inside the selection string and neither start nor end are used further in the code. helanal_trajectory() uses the resid of the first and last residue extensively for reporting so we now get these resids from the selection itself. --- package/CHANGELOG | 6 +- package/MDAnalysis/analysis/helanal.py | 160 +++++++++++++------------ 2 files changed, 85 insertions(+), 81 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index aab5164bf18..f78a51d39ca 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -137,7 +137,7 @@ Changes parameters anymore (Issue #1025) * atoms.rotate can be given center of rotation (Issue #1022) * XTCFile and TRRFile only raise IOError now on error. - * Deprecate usage of MDAnalaysis.core.AtomGroup + * Deprecate usage of MDAnalysis.core.AtomGroup * get_writer_for() returns NullWriter when filename=None instead of raising TypeError and filename is now a required arg instead of kwarg * moved coordinates.base.ChainReader to coordinates.chain.ChainReader @@ -145,8 +145,8 @@ Changes * totaltime now considers the first frame to be at time 0 (Issue #1137) * analysis.rms.RMSF now conforms to standard analysis API (PR #1136) * analysis.rms/align function keyword `mass_weighted` is replaced - by generic `weights='mass'`. This new keyword allows to pass it an - arbitrary weights array as well. (PR 1136) + by generic `weights='mass'`. This new keyword allows to pass it an + arbitrary weights array as well. (PR #1136) * Renamed various base classes for clarity. Iobase -> IOBase, Reader -> ReaderBase, SingleFrameReader -> SingleFrameReaderBase, Writer -> WriterBase, TopologyReader -> TopologyReaderBase, diff --git a/package/MDAnalysis/analysis/helanal.py b/package/MDAnalysis/analysis/helanal.py index f9f86bd8c03..1c3903e371e 100644 --- a/package/MDAnalysis/analysis/helanal.py +++ b/package/MDAnalysis/analysis/helanal.py @@ -158,7 +158,7 @@ def mean_abs_dev(a, mean_a=None): return np.mean(np.fabs(a - mean_a)) def helanal_trajectory(universe, selection="name CA", - start=None, end=None, begin=None, finish=None, + begin=None, finish=None, matrix_filename="bending_matrix.dat", origin_pdbfile="origin.pdb", summary_filename="summary.txt", @@ -169,70 +169,72 @@ def helanal_trajectory(universe, selection="name CA", twist_filename="unit_twist.xvg", prefix="helanal_", ref_axis=None, verbose=None, quiet=None): - """Perform HELANAL_ helix analysis on all frames in *universe*. + """Perform HELANAL helix analysis on all frames in `universe`. + + Parameters + ---------- + universe : Universe + selection : str (optional) + selection string that selects Calpha atoms [``"name CA"``] + begin : float (optional) + start analysing for time (ps) >= *begin*; ``None`` starts from the + beginning [``None``] + finish : float (optional) + stop analysis for time (ps) =< *finish*; ``None`` goes to the + end of the trajectory [``None``] + matrix_filename : str (optional) + Output file- bending matrix [``"bending_matrix.dat"``] + origin_pdbfile : str (optional) + Output file- origin pdb file [``"origin.pdb"``] + summary_filename : str (optional) + Output file- all of the basic data [``"summary.txt"``] + screw_filename : str (optional) + Output file- local tilts of individual residues from 2 to n-1 + [``"screw.xvg"``] + tilt_filename : str (optional) + Output file- tilt of line of best fit applied to origin axes + [``"local_tilt.xvg"``] + bend_filename : str (optional) + Output file- local bend angles between successive local helix axes + [``"local_bend.xvg"``] + twist_filename : str (optional) + Output file- local unit twist between successive helix turns + [``"unit_twist.xvg"``] + prefix : str (optional) + Prefix to add to all output file names; set to ``None`` to disable + [``"helanal__"``] + ref_axis : array_like (optional) + Calculate tilt angle relative to the axis; if ``None`` then ``[0,0,1]`` + is chosen [``None``] + verbose : bool (optional) + Toggle diagnostic outputs. [``True``] + + Raises + ------ + FinishTimeException + If the specified finish time precedes the specified start time or + current time stamp of trajectory object. .. Note:: Only a single helix is analyzed. Use the selection to specify the helix, e.g. with "name CA and resid 1:20" or use start=1, stop=20. - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` - - :Keywords: - *selection* - selection string that selects Calpha atoms [``"name CA"``] - *start* - start residue resid - *end* - end residue resid - *begin* - start analysing for time (ps) >= *begin*; ``None`` starts from the - beginning [``None``] - *finish* - stop analysis for time (ps) =< *finish*; ``None`` goes to the - end of the trajectory [``None``] - *matrix_filename* - Output file- bending matrix [``"bending_matrix.dat"``] - *origin_pdbfile* - Output file- origin pdb file [``"origin.pdb"``] - *summary_filename* - Output file- all of the basic data [``"summary.txt"``] - *screw_filename* - Output file- local tilts of individual residues from 2 to n-1 - [``"screw.xvg"``] - *tilt_filename* - Output file- tilt of line of best fit applied to origin axes - [``"local_tilt.xvg"``] - *bend_filename* - Output file- local bend angles between successive local helix axes - [``"local_bend.xvg"``] - *twist_filename* - Output file- local unit twist between successive helix turns - [``"unit_twist.xvg"``] - *prefix* - Prefix to add to all output file names; set to ``None`` to disable - [``"helanal__"``] - *ref_axis* - Calculate tilt angle relative to the axis; if ``None`` then ``[0,0,1]`` - is chosen [``None``] - *verbose* - Toggle diagnostic outputs. [``True``] - - :Raises: - FinishTimeException - If the specified finish time precedes the specified start time or - current time stamp of trajectory object. .. versionchanged:: 0.13.0 - New *quiet* keyword to silence frame progress output and most of the + New `quiet` keyword to silence frame progress output and most of the output that used to be printed to stdout is now logged to the logger *MDAnalysis.analysis.helanal* (at logelevel *INFO*). - .. deprecated:: 0.16 - The *quiet* keyword argument is deprecated in favor of the New - *verbose* one. + .. versionchanged:: 0.16.0 + Removed the `start` and `end` keywords for selecting residues because this can + be accomplished more transparently with `selection`. The first and last resid + are directly obtained from the selection. + + .. deprecated:: 0.16.0 + The `quiet` keyword argument is deprecated in favor of the new + `verbose` one. + """ verbose = _set_verbose(verbose, quiet, default=True) if ref_axis is None: @@ -242,13 +244,8 @@ def helanal_trajectory(universe, selection="name CA", # two atoms ref_axis = np.asarray(ref_axis) - if not (start is None and end is None): - if start is None: - start = universe.atoms[0].resid - if end is None: - end = universe.atoms[-1].resid - selection += " and resid {start:d}:{end:d}".format(**vars()) ca = universe.select_atoms(selection) + start, end = ca.resids[[0, -1]] trajectory = universe.trajectory if finish is not None: @@ -483,8 +480,8 @@ def stats(some_list): return [list_mean, list_sd, list_abdev] -def helanal_main(pdbfile, selection="name CA", start=None, end=None, ref_axis=None): - """Simple HELANAL run on a single frame PDB/GRO. +def helanal_main(pdbfile, selection="name CA", ref_axis=None): + """Simple HELANAL_ run on a single frame PDB/GRO. Computed data are returned as a dict and also logged at level INFO to the logger *MDAnalysis.analysis.helanal*. A simple way to enable a logger is to @@ -493,21 +490,30 @@ def helanal_main(pdbfile, selection="name CA", start=None, end=None, ref_axis=No .. Note:: Only a single helix is analyzed. Use the selection to specify the helix, e.g. with "name CA and resid 1:20". + Parameters + ---------- + pdbfile : str + filename of the single-frame input file + selection : str (optional) + selection string, default is "name CA" to select all C-alpha atoms. + ref_axis : array_like (optional) + Calculate tilt angle relative to the axis; if ``None`` then ``[0,0,1]`` + is chosen [``None``] + Returns ------- - - Dict with keys for - * Height: mean, stdev, abs dev - * Twist: mean, stdev, abs dev - * Residues/turn: mean, stdev, abs dev - * Local bending angles: array for computed angles (per residue) - * Unit twist angles: array for computed angles (per residue) - * Best fit tilt - * Rotation angles: local screw angles (per residue) + result : dict + The `result` contains keys + * Height: mean, stdev, abs dev + * Twist: mean, stdev, abs dev + * Residues/turn: mean, stdev, abs dev + * Local bending angles: array for computed angles (per residue) + * Unit twist angles: array for computed angles (per residue) + * Best fit tilt + * Rotation angles: local screw angles (per residue) Example ------- - Analyze helix 8 in AdK (PDB 4AKE); the standard logger is started and writes output to the file ``MDAnalysis.log``:: @@ -518,15 +524,13 @@ def helanal_main(pdbfile, selection="name CA", start=None, end=None, ref_axis=No .. versionchanged:: 0.13.0 All output is returned as a dict and logged to the logger *MDAnalysis.analysis.helanal* instead of being printed to stdout. + + .. versionchanged:: 0.16.0 + Removed the `start` and `end` keywords for selecting residues because this can + be accomplished more transparently with `selection`. """ universe = MDAnalysis.Universe(pdbfile) - if not (start is None and end is None): - if start is None: - start = universe.atoms[0].resid - if end is None: - end = universe.atoms[-1].resid - selection += " and resid {start:d}:{end:d}".format(**vars()) ca = universe.select_atoms(selection) logger.info("Analysing %d/%d residues", ca.n_atoms, universe.atoms.n_residues) From d396244148f5529bdd9b1dfc64a7f01f3ea33558 Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Thu, 30 Mar 2017 02:39:59 -0700 Subject: [PATCH 5/9] analysis.leaflet: numpified docs and added test for optimize_cutoff() --- package/MDAnalysis/analysis/leaflet.py | 150 ++++++++++-------- .../MDAnalysisTests/analysis/test_leaflet.py | 10 +- 2 files changed, 97 insertions(+), 63 deletions(-) diff --git a/package/MDAnalysis/analysis/leaflet.py b/package/MDAnalysis/analysis/leaflet.py index ae8b3c8bc4a..af957f3a5ac 100644 --- a/package/MDAnalysis/analysis/leaflet.py +++ b/package/MDAnalysis/analysis/leaflet.py @@ -1,5 +1,5 @@ # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.mdanalysis.org # Copyright (c) 2006-2016 The MDAnalysis Development Team and contributors @@ -25,14 +25,10 @@ Leaflet identification --- :mod:`MDAnalysis.analysis.leaflet` ============================================================== -:Author: Oliver Beckstein -:Year: 2010 -:Copyright: GNU Public License v3 - -Algorithm: - 1. build a graph of all phosphate distances < cutoff - 2. identify the largest connected subgraphs - 3. analyse first and second largest graph, which correspond to the leaflets +This module implements the *LeafletFinder* algorithm, described in +[Michaud-Agrawal2011]_. It can identify the lipids in a bilayer of +arbitrary shape and topology, including planar and undulating bilayers +under periodic boundary conditions or vesicles. One can use this information to identify @@ -51,6 +47,20 @@ .. MDAnalysisCookbook_: https://github.com/MDAnalysis/MDAnalysisCookbook/tree/master/examples + +Algorithm +--------- + +1. build a graph of all phosphate distances < cutoff +2. identify the largest connected subgraphs +3. analyse first and second largest graph, which correspond to the leaflets + +For further details see [Michaud-Agrawal2011]_. + + +Classes and Functions +--------------------- + .. autoclass:: LeafletFinder :members: @@ -73,6 +83,31 @@ class LeafletFinder(object): """Identify atoms in the same leaflet of a lipid bilayer. + This class implements the *LeafletFinder* algorithm [Michaud-Agrawal2011]_. + + Parameters + ---------- + universe : Universe or str + :class:`MDAnalysis.Universe` or a file name (e.g., in PDB or + GRO format) + selection : AtomGroup or str + A AtomGroup instance or a + :meth:`Universe.select_atoms` selection string + for atoms that define the lipid head groups, e.g. + universe.atoms.PO4 or "name PO4" or "name P*" + cutoff : float (optional) + head group-defining atoms within a distance of `cutoff` + Angstroms are deemed to be in the same leaflet [15.0] + pbc : bool (optional) + take periodic boundary conditions into account [``False``] + sparse : bool (optional) + ``None``: use fastest possible routine; ``True``: use slow + sparse matrix implementation (for large systems); ``False``: + use fast :func:`~MDAnalysis.lib.distances.distance_array` + implementation [``None``]. + + Example + ------- The components of the graph are stored in the list :attr:`LeafletFinder.components`; the atoms in each component are numbered consecutively, starting at 0. To obtain the atoms in the input structure @@ -93,29 +128,6 @@ class LeafletFinder(object): """ def __init__(self, universe, selectionstring, cutoff=15.0, pbc=False, sparse=None): - """Initialize from a *universe* or pdb file. - - :Arguments: - *universe* - :class:`MDAnalysis.Universe` or a PDB file name. - *selection* - :class:`MDAnalysis.core.groups.AtomGroup` or a - :meth:`MDAnalysis.Universe.select_atoms` selection string - for atoms that define the lipid head groups, e.g. - universe.atoms.PO4 or "name PO4" or "name P*" - :Keywords: - *cutoff* - head group-defining atoms within a distance of *cutoff* - Angstroms are deemed to be in the same leaflet [15.0] - *pbc* - take periodic boundary conditions into account (only works - for orthorhombic boxes) [``False``] - *sparse* - ``None``: use fastest possible routine; ``True``: use slow - sparse matrix implementation (for large systems); ``False``: - use fast :func:`~MDAnalysis.analysis.distances.distance_array` - implementation [``None``]. - """ universe = core.universe.as_Universe(universe) self.universe = universe self.selectionstring = selectionstring @@ -218,48 +230,62 @@ def write_selection(self, filename, **kwargs): options. """ sw = selections.get_writer(filename, kwargs.pop('format', None)) - writer = sw( - filename, mode=kwargs.pop('mode', 'wa'), - preamble="leaflets based on selection={selectionstring!r} cutoff={cutoff:f}\n".format(**vars(self)), - **kwargs) - for i, ag in enumerate(self.groups_iter()): - name = "leaflet_{0:d}".format((i + 1)) - writer.write(ag, name=name) + with sw(filename, mode=kwargs.pop('mode', 'w'), + preamble="leaflets based on selection={selectionstring!r} cutoff={cutoff:f}\n".format( + **vars(self)), + **kwargs) as writer: + for i, ag in enumerate(self.groups_iter()): + name = "leaflet_{0:d}".format((i + 1)) + writer.write(ag, name=name) def __repr__(self): - return "".format(self.selectionstring, self.cutoff, self.selection.n_atoms, - len(self.components)) + return "".format( + self.selectionstring, self.cutoff, self.selection.n_atoms, + len(self.components)) def optimize_cutoff(universe, selection, dmin=10.0, dmax=20.0, step=0.5, max_imbalance=0.2, **kwargs): - """Find cutoff that minimizes number of disconnected groups. + r"""Find cutoff that minimizes number of disconnected groups. Applies heuristics to find best groups: 1. at least two groups (assumes that there are at least 2 leaflets) 2. reject any solutions for which: - `|N0 - N1|/|N0 + N1| > *max_imbalance*` - - Ni = number of lipids in group i. This heuristic picks groups with - balanced numbers of lipids. - - :Arguments: - *universe* - :class:`MDAnalysis.Universe` instance - *selection* - selection string as used for :class:`LeafletFinder` - *dmin*, *dmax*, *step* - scan cutoffs from *dmin* to *dmax* at stepsize*step (in Angstroms) - *max_imbalance* - tuning parameter for the balancing heuristic (2) [0.2] - *kwargs* - other arguments for :class:`LeafletFinder` - - :Returns: ``(cutoff,N)`` optimum cutoff and number of groups found - :Raises: can die in various ways if really no appropriate number of groups - can be found; needs to be made more robust + .. math:: + + \frac{|N_0 - N_1|}{|N_0 + N_1|} > \mathrm{max_imbalance} + + with :math:`N_i` being the number of lipids in group + :math:`i`. This heuristic picks groups with balanced numbers of + lipids. + + Parameters + ---------- + universe : Universe + :class:`MDAnalysis.Universe` instance + selection : AtomGroup or str + AtomGroup or selection string as used for :class:`LeafletFinder` + dmin : float (optional) + dmax : float (optional) + step : float (optional) + scan cutoffs from `dmin` to `dmax` at stepsize `step` (in Angstroms) + max_imbalance : float (optional) + tuning parameter for the balancing heuristic [0.2] + kwargs : other keyword arguments + other arguments for :class:`LeafletFinder` + + Returns + ------- + (cutoff, N) + optimum cutoff and number of groups found + + + .. Note:: This function can die in various ways if really no + appropriate number of groups can be found; it ought to be + made more robust. + """ kwargs.pop('cutoff', None) # not used, so we filter it _sizes = [] diff --git a/testsuite/MDAnalysisTests/analysis/test_leaflet.py b/testsuite/MDAnalysisTests/analysis/test_leaflet.py index 25dcb507811..c1eb79092a6 100644 --- a/testsuite/MDAnalysisTests/analysis/test_leaflet.py +++ b/testsuite/MDAnalysisTests/analysis/test_leaflet.py @@ -22,7 +22,7 @@ import MDAnalysis from MDAnalysisTests import module_not_found -from numpy.testing import TestCase, assert_equal, dec +from numpy.testing import TestCase, assert_equal, assert_almost_equal, dec import numpy as np from MDAnalysisTests.datafiles import Martini_membrane_gro @@ -37,6 +37,7 @@ def setUp(self): def tearDown(self): del self.universe + del self.lipid_heads del self.lipid_head_string def test_leaflet_finder(self): @@ -58,3 +59,10 @@ def test_string_vs_atomgroup_proper(self): groups_string = lfls_string.groups() assert_equal(groups_string[0].indices, groups_ag[0].indices) assert_equal(groups_string[1].indices, groups_ag[1].indices) + + def test_optimize_cutoff(self): + from MDAnalysis.analysis.leaflet import optimize_cutoff + cutoff, N = optimize_cutoff(self.universe, self.lipid_heads, pbc=True) + assert_equal(N, 2) + assert_almost_equal(cutoff, 10.5, decimal=4) + From ff5f7d26a1f3d7e649d2cf51138c5dfcee56e359 Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Fri, 31 Mar 2017 02:57:39 -0700 Subject: [PATCH 6/9] analysis.legacy.x3dna: numpified docs and use OrderedDict --- package/MDAnalysis/analysis/legacy/x3dna.py | 237 +++++++++++++------- 1 file changed, 155 insertions(+), 82 deletions(-) diff --git a/package/MDAnalysis/analysis/legacy/x3dna.py b/package/MDAnalysis/analysis/legacy/x3dna.py index 68b9741a3a5..dcd3df4f99b 100644 --- a/package/MDAnalysis/analysis/legacy/x3dna.py +++ b/package/MDAnalysis/analysis/legacy/x3dna.py @@ -121,7 +121,7 @@ Utilities --------- -.. autoclass:: ApplicationError +.. autoexception:: ApplicationError """ from __future__ import print_function @@ -137,6 +137,7 @@ import subprocess import tempfile import textwrap +from collections import OrderedDict from MDAnalysis import ApplicationError from MDAnalysis.lib.util import which, realpath, asiterable @@ -147,13 +148,28 @@ def mean_std_from_x3dnaPickle(profile): - """Get mean and standard deviation of helicoidal parameters from a saved *profile*. + """Get mean and standard deviation of helicoidal parameters from a saved `profile`. - The *profile* should have been saved with :meth:`BaseX3DNA.save`. Then + The `profile` should have been saved with :meth:`BaseX3DNA.save`. Then load it with :: profile = cPickle.load(open("x3dna.pickle")) h_mean, h_std = mean_std_from_x3dnaPickle(profile) + + Arguments + --------- + profile : dict + A :attr:`X3DNA.profiles` dict with results from the + :class:`X3DNA` analysis. + + Returns + ------- + (list, list) + The tuple contains two lists with the means and the standard deviations + for the helicoidal parameters. The order for both lists is ``[shear, + stretch, stagger, buckle, propeller, opening, shift, slide, rise, tilt, + roll, twist]``. + """ if profile.x3dna_param is False: bp_shear, bp_stretch, bp_stagger, bp_rise, bp_shift, bp_slide, bp_buckle, bp_prop, bp_open, bp_tilt, bp_roll,\ @@ -217,8 +233,30 @@ def mean_std_from_x3dnaPickle(profile): class BaseX3DNA(object): """Baseclass for X3DNA_ analysis, providing plotting and utility functions. + When methods return helicoidal basepair parameter as lists, then the order + is always + + ====== ============== + index parameter + ====== ============== + 0 shear + 1 stretch + 2 stagger + 3 buckle + 4 propeller + 5 opening + 6 shift + 7 slide + 8 rise + 9 tilt + 10 roll + 11 twist + ====== ============== + + for each nucleic acid pair. .. _X3DNA: http://x3dna.org + """ def save(self, filename="x3dna.pickle"): @@ -235,9 +273,15 @@ def save(self, filename="x3dna.pickle"): cPickle.dump(self.profiles, open(filename, "wb"), cPickle.HIGHEST_PROTOCOL) def mean_std(self): - """H.mean_std() returns the mean and standard deviation of base parameters (order of base parameters:Shear, - Stretch,Stagger,Buckle,Propeller,Opening,Shift,Slide,Rise,Tilt,Roll,Twist) for each nucleic acid pair. - + """Returns the mean and standard deviation of base parameters. + + Returns + ------- + (list, list) + The tuple contains two lists with the means and the standard deviations + for the helicoidal parameters. The order for both lists is ``[shear, + stretch, stagger, buckle, propeller, opening, shift, slide, rise, tilt, + roll, twist]``. """ bp_shear, bp_stretch, bp_stagger, bp_rise, bp_shift, bp_slide, bp_buckle, bp_prop, bp_open, bp_tilt, bp_roll,\ @@ -275,8 +319,14 @@ def mean_std(self): return na_avg, na_std def mean(self): - """H.mean() returns the mean value for the base parameters (order of base parameters:Shear,Stretch,Stagger, - Buckle,Propeller,Opening,Shift,Slide,Rise,Tilt,Roll,Twist) for each nucleic acid pair. + """Returns the mean value for the base parameters. + + Returns + ------- + list + The list contains the means for the helicoidal parameters. The + order is ``[shear, stretch, stagger, buckle, propeller, opening, + shift, slide, rise, tilt, roll, twist]``. """ bp_shear, bp_stretch, bp_stagger, bp_rise, bp_shift, bp_slide, bp_buckle, bp_prop, bp_open, bp_tilt, bp_roll,\ @@ -310,11 +360,16 @@ def mean(self): return na_avg def std(self): - """H.std() returns the standard deviation of base parameters (order of base parameters:Shear,Stretch,Stagger, - Buckle,Propeller,Opening,Shift,Slide,Rise,Tilt,Roll,Twist) for each nucleic acid pair. + """Returns the standard deviation for the base parameters. - """ + Returns + ------- + list + The list contains the standard deviations for the helicoidal + parameters. The order is ``[shear, stretch, stagger, buckle, + propeller, opening, shift, slide, rise, tilt, roll, twist]``. + """ bp_shear, bp_stretch, bp_stagger, bp_rise, bp_shift, bp_slide, bp_buckle, bp_prop, bp_open, bp_tilt, bp_roll,\ bp_twist = [], [], [], [], [], [], [], [], [], [], [], [] for i in range(len(self.profiles)): @@ -346,11 +401,16 @@ def std(self): return na_std def plot(self, **kwargs): - """Plot base parameter profiles in a 1D graph - (plotting of the mean and standard deviation parameter value for each individual base pair). + """Plot time-averaged base parameters for each basse pair in a 1D graph. + + One plot is produced for each parameter. It shows the the mean and + standard deviation for each individual base pair. Each plot is saved to + PNG file with name ".png". - :Keywords: - *Not available at this time* + Parameters + ---------- + ax : matplotlib.pyplot.Axes (optional) + Provide `ax` to have all plots plotted in the same axes. """ import matplotlib.pyplot as plt @@ -371,9 +431,10 @@ def plot(self, **kwargs): ax.figure.clf() def sorted_profiles_iter(self): - """Return an iterator over profiles sorted by frame/order parameter *q*. + """Return an iterator over profiles sorted by frame/order parameter. - The iterator produces tuples ``(q, profile)``. + The iterator produces tuples ``(q, profile)``. Typically, `q` is the + frame number. """ if self.profiles is None: raise StopIteration @@ -398,6 +459,26 @@ class X3DNA(BaseX3DNA): The class also provides some simple plotting functions of the collected data such as :meth:`X3DNA.plot` or :meth:`X3DNA.plot3D`. + When methods return helicoidal basepair parameter as lists, then the order + is always + + ====== ============== + index parameter + ====== ============== + 0 shear + 1 stretch + 2 stagger + 3 buckle + 4 propeller + 5 opening + 6 shift + 7 slide + 8 rise + 9 tilt + 10 roll + 11 twist + ====== ============== + .. versionadded:: 0.8 .. _`X3DNA docs`: http://forum.x3dna.org/ @@ -406,33 +487,29 @@ class X3DNA(BaseX3DNA): def __init__(self, filename, **kwargs): """Set up parameters to run X3DNA_ on PDB *filename*. - :Arguments: - - *filename* - - The *filename* is used as input for X3DNA in the "xdna_ensemble" - command. It specifies the name of a PDB coordinate file - to be used. This must be in Brookhaven protein databank format - or something closely approximating this. Both ATOM and HETATM - records are read. Note that if water molecules or ions are - present in the channel these can be ignored on read by the use - of the *ignore_residues* keyword. - - :Keywords: - - *executable* - + Parameters + ---------- + filename : str + The `filename` is used as input for X3DNA in the + :program:`xdna_ensemble` command. It specifies the name of a + PDB coordinate file to be used. This must be in Brookhaven + protein databank format or something closely approximating + this. + executable : str (optional) Path to the :program:`xdna_ensemble` executable directories (e.g. ``/opt/x3dna/2.1 and /opt/x3dna/2.1/bin``) must be set and then added to export in bashrc file. See X3DNA documentation for set-up instructions. + x3dna_param : bool (optional) + Determines whether base step or base pair parameters will be + calculated. If ``True`` (default) then stacked *base step* + parameters will be analyzed. If ``False`` then stacked *base + pair* parameters will be analyzed. + logfile : str (optional) + Write output from X3DNA to `logfile` (default: "bp_step.par") - *x3dna_param* - Determines whether base step or base pair parameters will be - calculated. If True then stacked base step parameters will be - analyzed [Default is True]. If False then stacked base pair - parameters will be analyzed. + .. SeeAlso:: :class:`X3DNAtraj` """ # list of temporary files, to be cleaned up on __del__ self.tempfiles = [ @@ -469,7 +546,7 @@ def __init__(self, filename, **kwargs): logger.error("Executable %(program)r not found, should have been %(path)r.", vars()) # results - self.profiles = {} + self.profiles = OrderedDict() def run(self, **kwargs): """Run X3DNA on the input file.""" @@ -503,33 +580,29 @@ def run(self, **kwargs): def collect(self, **kwargs): """Parse the output from a X3DNA run into numpy recarrays. - Can deal with outputs containing multiple frames. Output format: + Can deal with outputs containing multiple frames. The method saves the result as :attr:`X3DNA.profiles`, a dictionary indexed by the frame number. Each entry is a :class:`np.recarray`. - If the keyword *outdir* is supplied (e.g. ".") then each profile is + If the keyword `outdir` is supplied (e.g. ".") then each profile is saved to a gzipped data file. - :Keywords: - *run* + Parameters + ---------- + run : str, int (optional identifier, free form [1] - *outdir* - save output data under *outdir*/*run* if set to any other + outdir : str (optional) + save output data under `outdir`/`run` if set to any other value but ``None`` [``None``] """ - # Shear Stretch Stagger Buckle Prop-Tw Opening Shift Slide Rise Tilt - # Roll Twist - #0123456789.0123456789.0123456789.0123456789.0123456789.0123456789. + # Shear Stretch Stagger Buckle Prop-Tw Opening Shift Slide Rise Tilt Roll Twist + #0123456789.0123456789.0123456789.0123456789.0123456789.0123456789.123456789.123456789.123456789.123456789.123456789.123456789.123456789. # 11 22 33 44 - #123456789.123456789.123456789.123456789.123456789.123456789.123456789. - #1 13 - #T-A -0.033 -0.176 0.158 -12.177 -8.979 1.440 0.000 0.000 0.000 0.000 - # 0.000 0.000 - #C-G -0.529 0.122 -0.002 -7.983 -10.083 -0.091 -0.911 1.375 3.213 -0.766 - # -4.065 41.492 + #T-A -0.033 -0.176 0.158 -12.177 -8.979 1.440 0.000 0.000 0.000 0.000 0.000 0.000 + #C-G -0.529 0.122 -0.002 -7.983 -10.083 -0.091 -0.911 1.375 3.213 -0.766 -4.065 41.492 # only parse bp_step.par x3dna_output = kwargs.pop("x3dnaout", self.logfile) @@ -549,8 +622,7 @@ def collect(self, **kwargs): logger.info("Found %d input files based on glob pattern %s", length, self.filename) # one recarray for each frame, indexed by frame number - # TODO: use ordered dict - self.profiles = {} + self.profiles = OrderedDict() logger.info("Run %s: Reading %d X3DNA profiles from %r", run, length, x3dna_output) x3dna_profile_no = 0 @@ -651,31 +723,39 @@ class X3DNAtraj(BaseX3DNA): trajectory formats understood by MDAnalysis. This class can take any universe and feed it to X3DNA. By default it sequentially creates a PDB for each frame and runs X3DNA on the frame. + """ def __init__(self, universe, **kwargs): """Set up the class. - :Arguments: - - *universe* + Parameters + ---------- + universe : Universe The input trajectory as part of a - :class:`~MDAnalysis.core.universe.Universe`. trajectory is + :class:`~MDAnalysis.core.universe.Universe`; the trajectory is converted to a sequence of PDB files and X3DNA is run on each - individual file. (Use the *start*, *stop*, and *step* keywords + individual file. (Use the `start`, `stop`, and `step` keywords to slice the trajectory.) - - *start*, *stop*, *step* + selection : str (optional) + MDAnalysis selection string (default: "nucleic") to select the + atoms that should be analyzed. + start : int (optional) + stop : int (optional) + step : int (optional) frame indices to slice the trajectory as - ``universe.trajectory[start, stop, step]`` - - *x3dna_param* + ``universe.trajectory[start, stop, step]``; by default, the whole + trajectory is analyzed. + x3dna_param : bool (optional) indicates whether stacked bases or stacked base-pairs will be - analyzed. True is bases or False is stacked base-pairs [Default - is True]. + analyzed. ``True`` is bases and ``False`` is stacked base-pairs + [Default is ``True``]. + kwargs : keyword arguments (optional) + All other keywords are passed on to :class:`X3DNA` (see there + for description). + - *kwargs* - All other keywords are passed on to :class:`X3DNA` (see there for description). + .. SeeAlso:: :class:`X3DNA` """ self.universe = universe @@ -693,8 +773,9 @@ def __init__(self, universe, **kwargs): def run(self, **kwargs): """Run X3DNA on the whole trajectory and collect profiles. - Keyword arguments *start*, *stop*, and *step* can be used to only - analyse part of the trajectory. + Keyword arguments `start`, `stop`, and `step` can be used to only + analyse part of the trajectory. The defaults are the values provided to + the class constructor. """ start = kwargs.pop('start', self.start) stop = kwargs.pop('stop', self.stop) @@ -703,12 +784,10 @@ def run(self, **kwargs): x3dna_kw = self.x3dna_kwargs.copy() x3dna_kw.update(kwargs) - profiles = {} + profiles = OrderedDict() nucleic = self.universe.select_atoms(self.selection) for ts in self.universe.trajectory[start:stop:step]: - print(ts.frame) - print(nucleic.atoms) logger.info("X3DNA analysis frame %4d ", ts.frame) fd, pdbfile = tempfile.mkstemp(suffix=".pdb") os.close(fd) @@ -716,12 +795,6 @@ def run(self, **kwargs): os.system("find_pair {0!s} 355d.bps".format(pdbfile)) try: nucleic.write(pdbfile) - ''' - if int(ts.frame) in [int(start+1), 1]: - #print start - os.system("find_pair %s 355d.bps" %(pdbfile)) - #print "complete" - ''' x3dna_profiles = self.run_x3dna(pdbfile, **x3dna_kw) finally: try: @@ -737,7 +810,7 @@ def run(self, **kwargs): self.profiles = profiles def run_x3dna(self, pdbfile, **kwargs): - """Run X3DNA on a single PDB file *pdbfile*.""" + """Run X3DNA on a single PDB file `pdbfile`.""" kwargs['x3dna_param'] = self.x3dna_param H = X3DNA(pdbfile, **kwargs) H.run() From 056f0976100e02c82a354a76be79f0669549926e Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Fri, 31 Mar 2017 03:17:15 -0700 Subject: [PATCH 7/9] analysis.nuclinfo: numpyfied docs --- package/MDAnalysis/analysis/nuclinfo.py | 413 +++++++++++++++--------- 1 file changed, 264 insertions(+), 149 deletions(-) diff --git a/package/MDAnalysis/analysis/nuclinfo.py b/package/MDAnalysis/analysis/nuclinfo.py index 6847743462f..9d22733534b 100644 --- a/package/MDAnalysis/analysis/nuclinfo.py +++ b/package/MDAnalysis/analysis/nuclinfo.py @@ -20,9 +20,6 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -# MDAnalysis -- nucleic acid analysis -# Copyright (c) 2011 Elizabeth Denning - """ Nucleic acid analysis --- :mod:`MDAnalysis.analysis.nuclinfo` ============================================================= @@ -113,24 +110,34 @@ def wc_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): - """Watson-Crick basepair distance for residue *i* with residue *bp*. + """Watson-Crick basepair distance for residue `i` with residue `bp`. The distance of the nitrogen atoms in a Watson-Crick hydrogen bond is computed. - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *seg1* - segment id for first base - *i* - resid of the first base - *seg2* - segment id for second base - *bp* - resid of the second base + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + i : int + resid of the first base + bp : int + resid of the second base + seg1 : str (optional) + segment id for first base ["SYSTEM"] + seg2 : str (optional) + segment id for second base ["SYSTEM"] + + Returns + ------- + float + Watson-Crick base pair distance + + + Notes + ----- + If failure occurs be sure to check the segment identification. - .. NOTE:: If failure occurs be sure to check the segment identification. .. versionadded:: 0.7.6 """ @@ -146,24 +153,33 @@ def wc_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): def minor_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): - """Minor-Groove basepair distance for residue *i* with residue *bp*. + """Minor-Groove basepair distance for residue `i` with residue `bp`. The distance of the nitrogen and oxygen atoms in a Minor-groove hydrogen bond is computed. - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *seg1* - segment id for first base - *i* - resid of the first base - *seg2* - segment id for second base - *bp* - resid of the second base + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + i : int + resid of the first base + bp : int + resid of the second base + seg1 : str (optional) + segment id for first base ["SYSTEM"] + seg2 : str (optional) + segment id for second base ["SYSTEM"] + + Returns + ------- + float + Minor groove base pair distance + + Notes + ----- + If failure occurs be sure to check the segment identification. - .. NOTE:: If failure occurs be sure to check the segment identification. .. versionadded:: 0.7.6 """ @@ -179,25 +195,34 @@ def minor_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): def major_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): - """Major-Groove basepair distance for residue *i* with residue *bp*. + """Major-Groove basepair distance for residue `i` with residue `bp`. The distance of the nitrogen and oxygen atoms in a Major-groove hydrogen bond is computed. - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the - trajectory *i* - *seg1* - segment id for first base - *i* - resid of the first base - *seg2* - segment id for second base - *bp* - resid of the second base - - .. NOTE:: If failure occurs be sure to check the segment identification + + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + i : int + resid of the first base + bp : int + resid of the second base + seg1 : str (optional) + segment id for first base ["SYSTEM"] + seg2 : str (optional) + segment id for second base ["SYSTEM"] + + Returns + ------- + float + Major groove base pair distance + + Notes + ----- + If failure occurs be sure to check the segment identification. + .. versionadded:: 0.7.6 """ @@ -219,18 +244,25 @@ def major_pair(universe, i, bp, seg1="SYSTEM", seg2="SYSTEM"): def phase_cp(universe, seg, i): - """Pseudo-angle describing the phase of the ribose pucker for residue *i* using the CP method. + """Pseudo-angle describing the phase of the ribose pucker for residue `i` using the CP method. The angle is computed by the positions of atoms in the ribose ring. - :Arguments: - *universe* + + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + float + phase angle in degrees - *segid* - segment identity of resid - *i* - resid of the base .. versionadded:: 0.7.6 """ @@ -287,17 +319,24 @@ def phase_cp(universe, seg, i): def phase_as(universe, seg, i): - """Pseudo-angle describing the phase of the ribose pucker for residue *i* using the AS method + """Pseudo-angle describing the phase of the ribose pucker for residue `i` using the AS method The angle is computed by the position vector of atoms in the ribose ring. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segment identity of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + float + phase angle in degrees + .. versionadded:: 0.7.6 """ @@ -353,19 +392,32 @@ def phase_as(universe, seg, i): def tors(universe, seg, i): - """Backbone dihedrals includes alpha, beta, gamma, delta, epsilon, zeta, chi + """Calculation of nucleic backbone dihedral angles. + + The dihedral angles are alpha, beta, gamma, delta, epsilon, zeta, chi. - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position of atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + [alpha, beta, gamma, delta, epsilon, zeta, chi] : list of floats + torsion angles in degrees + + Notes + ----- + If failure occurs be sure to check the segment identification. + - .. NOTE:: If failure occurs be sure to check the segment identification + .. versionadded:: 0.7.6 """ a = universe.select_atoms(" atom {0!s} {1!s} O3\' ".format(seg, i - 1), @@ -433,15 +485,22 @@ def tors(universe, seg, i): def tors_alpha(universe, seg, i): """alpha backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + alpha : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -458,15 +517,22 @@ def tors_alpha(universe, seg, i): def tors_beta(universe, seg, i): """beta backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + beta : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -483,15 +549,22 @@ def tors_beta(universe, seg, i): def tors_gamma(universe, seg, i): """ Gamma backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + gamma : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -508,15 +581,22 @@ def tors_gamma(universe, seg, i): def tors_delta(universe, seg, i): """delta backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + delta : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -533,15 +613,22 @@ def tors_delta(universe, seg, i): def tors_eps(universe, seg, i): """Epsilon backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + epsilon : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -558,15 +645,22 @@ def tors_eps(universe, seg, i): def tors_zeta(universe, seg, i): """Zeta backbone dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. - :Arguments: - *universe* + Parameters + ---------- + universe : Universe :class:`~MDAnalysis.core.universe.Universe` containing the trajectory - *segid* - segid of resid - *i* - resid of the base + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + zeta : float + torsion angle in degrees + .. versionadded:: 0.7.6 """ @@ -583,16 +677,22 @@ def tors_zeta(universe, seg, i): def tors_chi(universe, seg, i): """chi nucleic acid dihedral - The dihedral is computed based on position atoms for resid *i*. + The dihedral is computed based on position atoms for resid `i`. + + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + chi : float + torsion angle in degrees - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the - trajectory - *segid* - segid of resid - *i* - resid of the base .. versionadded:: 0.7.6 """ @@ -614,19 +714,27 @@ def tors_chi(universe, seg, i): def hydroxyl(universe, seg, i): """2-hydroxyl dihedral. Useful only for RNA calculations. - .. Note:: This dihedral calculation will only work if using atom names as - documented by charmm force field parameters. + .. Note:: This dihedral calculation will only work if using atom + names as documented by charmm force field parameters, + namely "C1', C2', O2', H2'". + + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the trajectory + seg : str + segment id for base + i : int + resid of the first base + + Returns + ------- + hydroxyl_angle : float + torsion angle in degrees - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the - trajectory - *segid* - segid of resid - *i* - resid of the base .. versionadded:: 0.7.6 + """ h = universe.select_atoms(" atom {0!s} {1!s} C1\' ".format(seg, i), " atom {0!s} {1!s} C2\' ".format(seg, i), @@ -646,27 +754,34 @@ def pseudo_dihe_baseflip(universe, bp1, bp2, i, seg1="SYSTEM", seg2="SYSTEM", seg3="SYSTEM"): """pseudo dihedral for flipped bases. Useful only for nucleic acid base flipping - The dihedral is computed based on position atoms for resid *i* - - .. Note:: This dihedral calculation will only work if using atom names as - documented by charmm force field parameters. - - :Arguments: - *universe* - :class:`~MDAnalysis.core.universe.Universe` containing the - trajectory - *segid1* - segid of resid base pairing with bp2 - *bp1* - resid that base pairs with bp2 - *segid2* - segid same as that of segid of flipping resid - *bp2* - resid below the base that flips - *segid3* - segid of resid that flips - *i* - resid of the base that flips + The dihedral is computed based on position atoms for resid `i` + + .. Note:: This dihedral calculation will only work if using atom names as + documented by charmm force field parameters. + + Parameters + ---------- + universe : Universe + :class:`~MDAnalysis.core.universe.Universe` containing the + trajectory + bp1 : int + resid that base pairs with `bp2` + bp2 : int + resid below the base that flips + i : int + resid of the base that flips + segid1 : str (optional) + segid of resid base pairing with `bp2` + segid2 : str (optional) + segid, same as that of segid of flipping resid `i` + segid3 : str (optional) + segid of resid `i` that flips + + Returns + ------- + float + pseudo dihedral angle in degrees + .. versionadded:: 0.8.0 """ From 2e2017bd7d6796bc6f8cd5741abf3ac357b3c551 Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Fri, 31 Mar 2017 03:44:59 -0700 Subject: [PATCH 8/9] analysis.rdf: fixed numpy reST --- package/MDAnalysis/analysis/rdf.py | 47 ++++++++++++++++-------------- 1 file changed, 25 insertions(+), 22 deletions(-) diff --git a/package/MDAnalysis/analysis/rdf.py b/package/MDAnalysis/analysis/rdf.py index 6c78571ec4e..e7e22e327f5 100644 --- a/package/MDAnalysis/analysis/rdf.py +++ b/package/MDAnalysis/analysis/rdf.py @@ -24,12 +24,15 @@ Radial Distribution Functions --- :mod:`MDAnalysis.analysis.rdf` ================================================================ -Tools for calculating pair distribution functions +Tools for calculating pair distribution functions ("radial +distribution functions" or "RDF"). + +.. Not Implemented yet: +.. - Structure factor? +.. - Coordination number -# TODO - - Structure factor? - - Coordination number """ + import numpy as np from ..lib.util import blocks_of @@ -44,44 +47,44 @@ class InterRDF(AnalysisBase): Arguments --------- - g1 + g1 : AtomGroup First AtomGroup - g2 + g2 : AtomGroup Second AtomGroup - - Keywords - -------- - nbins + nbins : int (optional) Number of bins in the histogram [75] - range + range : tuple or list (optional) The size of the RDF [0.0, 15.0] - exclusion_block + exclusion_block : tuple (optional) A tuple representing the tile to exclude from the distance array. [None] - start - The frame to start at [0] - stop - The frame to end at [-1] - step - The step size through the trajectory in frames [0] + start : int (optional) + The frame to start at (default is first) + stop : int (optional) + The frame to end at (default is last) + step : int (optional) + The step size through the trajectory in frames (default is + every frame) Example ------- - First create the InterRDF object, by supplying two AtomGroups - then use the `run` method + First create the :class:`InterRDF` object, by supplying two + AtomGroups then use the :meth:`run` method :: rdf = InterRDF(ag1, ag2) rdf.run() - Results are available through the .bins and .rdf attributes + Results are available through the :attr:`bins` and :attr:`rdf` + attributes:: plt.plot(rdf.bins, rdf.rdf) The `exclusion_block` keyword allows the masking of pairs from within the same molecule. For example, if there are 7 of each - atom in each molecule, the exclusion mask (7, 7) can be used. + atom in each molecule, the exclusion mask `(7, 7)` can be used. .. versionadded:: 0.13.0 + """ def __init__(self, g1, g2, nbins=75, range=(0.0, 15.0), exclusion_block=None, From bec81f282cf0eaa2b48d24499ff2a5d120ce2bed Mon Sep 17 00:00:00 2001 From: Oliver Beckstein Date: Fri, 31 Mar 2017 04:06:44 -0700 Subject: [PATCH 9/9] analysis.waterdynamics: numpyfied docs (TODO: PEP8) --- package/MDAnalysis/analysis/waterdynamics.py | 398 +++++++++++-------- 1 file changed, 231 insertions(+), 167 deletions(-) diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index 6008f56430d..8fb1e7a423c 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -1,5 +1,5 @@ # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # # MDAnalysis --- http://www.mdanalysis.org # Copyright (c) 2006-2016 The MDAnalysis Development Team and contributors @@ -20,8 +20,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -""" -Water dynamics analysis --- :mod:`MDAnalysis.analysis.waterdynamics` +"""Water dynamics analysis --- :mod:`MDAnalysis.analysis.waterdynamics` ======================================================================= :Author: Alejandro Bernardin @@ -30,51 +29,66 @@ .. 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]_. +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). +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. +.. [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. +.. [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 +.. [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. +.. [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. +.. [Brodka1994] Aleksander Brodka (1994). Diffusion in restricted volume, + Molecular Physics, 1994, Vol. 82, No. 5, 1075-1078. -.. [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. +.. [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 +.. [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: -Examples --------- +Example use of the analysis classes +----------------------------------- HydrogenBondLifetimes ~~~~~~~~~~~~~~~~~~~~~ -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 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:: +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 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 HydrogenBondLifetimes as HBL @@ -90,7 +104,7 @@ for HBLc, HBLi in HBL_analysis.timeseries: print("{time} {HBLc} {time} {HBLi}".format(time=time, HBLc=HBLc, HBLi=HBLi)) time += 1 - + #we can also plot our data plt.figure(1,figsize=(18, 6)) @@ -110,15 +124,22 @@ plt.show() -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 lifetimes should decay faster than 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 lifetimes should decay faster than +intermittent. WaterOrientationalRelaxation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -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:: +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 WaterOrientationalRelaxation as WOR @@ -158,20 +179,24 @@ plt.title('WOR dip') plt.plot(range(0,time),[column[2] for column in WOR_analysis.timeseries]) - plt.show() + plt.show() -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. +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. AngularDistribution ~~~~~~~~~~~~~~~~~~~ -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):: +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 AngularDistribution as AD @@ -214,13 +239,17 @@ plt.show() -where P(cos(theta)) is the angular distribution or angular probabilities. +where `P(cos(theta))` is the angular distribution or angular probabilities. + 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:: + +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 MeanSquareDisplacement as MSD @@ -235,7 +264,7 @@ for msd in MSD_analysis.timeseries: print("{time} {msd}".format(time=time, msd=msd)) time += 1 - + #Plot plt.xlabel('time') plt.ylabel('MSD') @@ -243,12 +272,17 @@ plt.plot(range(0,time),MSD_analysis.timeseries) plt.show() + +.. _SP-examples: + 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:: + +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 SurvivalProbability as SP @@ -270,17 +304,19 @@ plt.title('Survival Probability') plt.plot(range(0,time),MSD_analysis.timeseries) plt.show() - -.. Output + +.. _Output: 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):: + +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 @@ -294,8 +330,9 @@ WaterOrientationalRelaxation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Water orientational relaxation (WOR) data is returned per window timestep, which is stored in -:attr:`WaterOrientationalRelaxation.timeseries`:: + +Water orientational relaxation (WOR) data is returned per window timestep, +which is stored in :attr:`WaterOrientationalRelaxation.timeseries`:: results = [ [ # time t0 @@ -309,8 +346,10 @@ AngularDistribution ~~~~~~~~~~~~~~~~~~~ + Angular distribution (AD) data is returned per vector, which is stored in -:attr:`AngularDistribution.graph`. In fact, AngularDistribution returns a histogram:: +:attr:`AngularDistribution.graph`. In fact, AngularDistribution returns a +histogram:: results = [ [ # OH vector values @@ -327,8 +366,10 @@ 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`:: + +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 @@ -337,8 +378,10 @@ 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`:: + +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 @@ -370,7 +413,6 @@ :members: :inherited-members: - """ from __future__ import print_function from six.moves import range, zip_longest @@ -383,47 +425,53 @@ class HydrogenBondLifetimes(object): - r""" - 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 - Relaxation" (HBPR). In the continuos case we have: + r"""Hydrogen bond lifetime analysis + + 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 continuous 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: + 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. + 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. + + + Parameters + ---------- + universe : Universe + Universe object + selection1 : str + Selection string for first selection [‘byres name OH2’]. + It could be any selection available in MDAnalysis, not just water. + selection2 : str + Selection string to analize its HBL against selection1 + t0 : int + frame where analysis begins + tf : int + frame where analysis ends + dtmax : int + Maximum dt size, `dtmax` < `tf` or it will crash. + nproc : int + Number of processors to use, by default is 1. - .. versionadded:: 0.11.0 - - :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. + .. versionadded:: 0.11.0 """ - def __init__(self,universe ,selection1 ,selection2, t0 , tf , dtmax, nproc = 1): + def __init__(self, universe, selection1, selection2, t0, tf, dtmax, nproc=1): self.universe = universe self.selection1 = selection1 self.selection2 = selection2 @@ -618,10 +666,12 @@ def run(self, **kwargs): self.timeseries = self._getGraphics(h_list.timeseries, self.t0, self.tf, self.dtmax) class WaterOrientationalRelaxation(object): - r""" - Function to evaluate the Water Orientational Relaxation proposed by Yu-ling Yeh - and Chung-Yuan Mou [Yeh1999_]. WaterOrientationalRelaxation indicates "how fast" water molecules are rotating - or changing direction. This is a time correlation function given by: + r"""Water orientation relaxation analysis + + Function to evaluate the Water Orientational Relaxation proposed by Yu-ling + Yeh 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:: C_{\hat u}(\tau)=\langle \mathit{P}_2[\mathbf{\hat{u}}(t_0)\cdot\mathbf{\hat{u}}(t_0+\tau)]\rangle @@ -629,23 +679,26 @@ class WaterOrientationalRelaxation(object): 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* - 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. + Parameters + ---------- + universe : Universe + Universe object + selection : str + Selection string for water [‘byres name OH2’]. + t0 : int + frame where analysis begins + tf : int + frame where analysis ends + dtmax : int + Maximum dt size, `dtmax` < `tf` or it will crash. + + + .. versionadded:: 0.11.0 """ - def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): + def __init__(self, universe, selection, t0, tf, dtmax, nproc=1): self.universe = universe self.selection = selection self.t0 = t0 @@ -766,7 +819,7 @@ def _sameMolecTandDT(self,selection,t0d,tf): sort = sorted(list(a.intersection(b))) return sort - def _selection_serial(self,universe,selection_str): + def _selection_serial(self, universe, selection_str): selection = [] for ts in universe.trajectory: selection.append(universe.select_atoms(selection_str)) @@ -776,7 +829,7 @@ def _selection_serial(self,universe,selection_str): # Second Legendre polynomial lg2 = lambda self,x : (3*x*x - 1)/2 - def run(self,**kwargs): + def run(self, **kwargs): """ Analyze trajectory and produce timeseries """ @@ -796,31 +849,33 @@ def run(self,**kwargs): class AngularDistribution(object): - 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). + r"""Angular distribution function analysis - .. versionadded:: 0.11.0 + 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/np/reference/generated/np.histogram.html + Parameters + ---------- + universe : Universe + Universe object + selection : str + Selection string to evaluate its angular distribution [‘byres name OH2’] + bins : int (optional) + Number of bins to create the histogram by means of :func:`numpy.histogram [40] + axis : {'x', 'y', 'z'} (optional) + Axis to create angle with the vector (HH, OH or dipole) and calculate + cosine theta ['z']. + .. versionadded:: 0.11.0 """ - def __init__(self,universe,selection_str,bins=40,nproc=1,axis="z"): + def __init__(self, universe, selection_str, bins=40, nproc=1, axis="z"): self.universe = universe self.selection_str = selection_str self.bins = bins @@ -932,9 +987,10 @@ def _selection_serial(self,universe,selection_str): class MeanSquareDisplacement(object): - r""" - Function to evaluate the Mean Square Displacement (MSD_). The MSD gives the average distance that - particles travels. The MSD is given by: + r"""Mean square displacement analysis + + 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 @@ -945,23 +1001,25 @@ class MeanSquareDisplacement(object): .. _MSD: http://en.wikipedia.org/wiki/Mean_squared_displacement - .. versionadded:: 0.11.0 - :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. + Parameters + ---------- + universe : Universe + Universe object + selection : str + Selection string for water [‘byres name OH2’]. + t0 : int + frame where analysis begins + tf : int + frame where analysis ends + dtmax : int + Maximum dt size, `dtmax` < `tf` or it will crash. + + .. versionadded:: 0.11.0 """ - def __init__(self,universe,selection,t0,tf,dtmax,nproc=1): + def __init__(self, universe, selection, t0, tf, dtmax, nproc=1): self.universe = universe self.selection = selection self.t0 = t0 @@ -1069,30 +1127,36 @@ def run(self,**kwargs): class SurvivalProbability(object): - 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: + r"""Survival probability analysis + + 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: .. 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. + 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* - 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 + Parameters + ---------- + universe : Universe + Universe object + selection : str + Selection string; any selection is allowed. With this selection you + define the region/zone where to analyze, e.g.: "selection_a" and "zone" + (see `SP-examples`_ ) + t0 : int + frame where analysis begins + tf : int + frame where analysis ends + dtmax : int + Maximum dt size, `dtmax` < `tf` or it will crash. + + + .. versionadded:: 0.11.0 """