diff --git a/package/AUTHORS b/package/AUTHORS index 69cebb5e184..cbe9a5c6343 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -148,6 +148,7 @@ Chronological list of authors - William Glass - Marcello Sega - Edis Jakupovic + - Nicholas Craven External code ------------- diff --git a/package/CHANGELOG b/package/CHANGELOG index f9c3aca686d..0261e85555c 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,11 +14,14 @@ The rules for this file: ------------------------------------------------------------------------------ ??/??/?? tylerjereddy, richardjgowers, IAlibay, hmacdope, orbeckst, cbouy, - lilyminium, daveminh, jbarnoud, yuxuanzhuang, VOD555, ianmkenney + lilyminium, daveminh, jbarnoud, yuxuanzhuang, VOD555, ianmkenney, + calcraven， xiki-tempula * 2.0.0 Fixes + * Fixed reading in masses and charges from a hoomdxml file + (Issue #2888, PR #2889) * Bond attribute is no longer set if PDB file contains no CONECT records (Issue #2832) * ResidueAttrs now have Atom as a target class by default; ICodes and @@ -40,7 +43,9 @@ Fixes was thrown when finding D-H pairs via the topology if `hydrogens` was an empty AtomGroup (Issue #2848) * Fixed the DMSParser, allowing the creation of multiple segids sharing - residues with identical resids (Issue #1387, PR #2872) + residues with identical resids (Issue #1387, PR #2872) + * H5MD files are now pickleable with H5PYPicklable (Issue #2890, PR #2894) + * Fixed Janin analysis residue filtering (include CYSH) (Issue #2898) Enhancements * Refactored analysis.helanal into analysis.helix_analysis @@ -64,6 +69,9 @@ Enhancements * Dead code removed from the TPR parser and increased test coverage (PR #2840) * TPR parser exposes the elements topology attribute (PR #2858, see Issue #2553) * Added `H5MDReader` to coordinate readers (Issue #762, PR #2787) + * Added new kwargs `select_remove` and `select_protein` to + analysis.dihedrals.Janin analysis to give user more fine grained control + over selections (PR #2899) * Added an RDKit converter that works for any input with all hydrogens explicit in the topology (Issue #2468, PR #2775) @@ -89,6 +97,7 @@ Changes * Changes the minimal NumPy version to 1.16.0 (Issue #2827, PR #2831) * Sets the minimal RDKit version for CI to 2020.03.1 (Issue #2827, PR #2831) * Removes deprecated waterdynamics.HydrogenBondLifetimes (PR #2842) + * Make NeighborSearch return empty atomgroup, residue, segments instead of list (Issue #2892, PR #2907) Deprecations diff --git a/package/MDAnalysis/analysis/dihedrals.py b/package/MDAnalysis/analysis/dihedrals.py index 7ad6fa68d5d..cbbf0c19bfa 100644 --- a/package/MDAnalysis/analysis/dihedrals.py +++ b/package/MDAnalysis/analysis/dihedrals.py @@ -21,7 +21,7 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # r"""Dihedral angles analysis --- :mod:`MDAnalysis.analysis.dihedrals` -=========================================================================== +================================================================= :Author: Henry Mull :Year: 2018 @@ -35,7 +35,7 @@ A list of time steps that contain angles of interest is generated and can be easily plotted if desired. For the :class:`~MDAnalysis.analysis.dihedrals.Ramachandran` and :class:`~MDAnalysis.analysis.dihedrals.Janin` classes, basic plots can be -generated using the method :meth:`Ramachandran.plot()` or :meth:`Janin.plot()`. +generated using the method :meth:`Ramachandran.plot` or :meth:`Janin.plot`. These plots are best used as references, but they also allow for user customization. @@ -72,9 +72,11 @@ ~~~~~~~~~~~~~~~~~~~~~ The :class:`~MDAnalysis.analysis.dihedrals.Ramachandran` class allows for the -quick calculation of phi and psi angles. Unlike the :class:`~MDanalysis.analysis.dihedrals.Dihedral` -class which takes a list of `atomgroups`, this class only needs a list of -residues or atoms from those residues. The previous example can repeated with:: +quick calculation of classical Ramachandran plots [Ramachandran1963]_ in the +backbone :math:`phi` and :math:`psi` angles. Unlike the +:class:`~MDanalysis.analysis.dihedrals.Dihedral` class which takes a list of +`atomgroups`, this class only needs a list of residues or atoms from those +residues. The previous example can repeated with:: u = mda.Universe(GRO, XTC) r = u.select_atoms("resid 5-10") @@ -85,7 +87,8 @@ class which takes a list of `atomgroups`, this class only needs a list of import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=plt.figaspect(1)) - R.plot(ax=ax, color='k', marker='s') + R.plot(ax=ax, color='k', marker='o', ref=True) + fig.tight_layout() as shown in the example :ref:`Ramachandran plot figure `. @@ -99,28 +102,29 @@ class which takes a list of `atomgroups`, this class only needs a list of trajectory (XTC). The contours in the background are the "allowed region" and the "marginally allowed" regions. -The Janin class works in the same way, only needing a list of residues; see the -:ref:`Janin plot figure ` as an example. To plot the data -yourself, the angles can be accessed using :attr:`Ramachandran.angles` or -:attr:`Janin.angles`. - -Reference plots can be added to the axes for both the Ramachandran and Janin -classes using the kwarg ``ref=True``. The Ramachandran reference data -(:data:`~MDAnalysis.analysis.data.filenames.Rama_ref`) and Janin reference data -(:data:`~MDAnalysis.analysis.data.filenames.Janin_ref`) were made using data -obtained from a large selection of 500 PDB files, and were analyzed using these -classes. The allowed and marginally allowed regions of the Ramachandran reference -plt have cutoffs set to include 90% and 99% of the data points, and the Janin -reference plot has cutoffs for 90% and 98% of the data points. The list of PDB -files used for the referece plots was taken from [Lovell2003]_ and information -about general Janin regions was taken from [Janin1978]_. +To plot the data yourself, the angles can be accessed using +:attr:`Ramachandran.angles`. .. Note:: - These classes are prone to errors if the topology contains duplicate or missing - atoms (e.g. atoms with `altloc` or incomplete residues). If the topology has as - an `altloc` attribute, you must specify only one `altloc` for the atoms with - more than one (``"protein and not altloc B"``). + The Ramachandran analysis is prone to errors if the topology contains + duplicate or missing atoms (e.g. atoms with `altloc` or incomplete + residues). If the topology has as an `altloc` attribute, you must specify + only one `altloc` for the atoms with more than one (``"protein and not + altloc B"``). + + +Janin analysis +~~~~~~~~~~~~~~ + +Janin plots [Janin1978]_ for side chain conformations (:math:`\chi_1` and +:math:`chi_2` angles) can be created with the +:class:`~MDAnalysis.analysis.dihedrals.Janin` class. It works in the same way, +only needing a list of residues; see the :ref:`Janin plot figure +` as an example. + +The data for the angles can be accessed in the attribute +:attr:`Janin.angles`. .. _figure-janin: @@ -128,10 +132,39 @@ class which takes a list of `atomgroups`, this class only needs a list of :scale: 50 % :alt: Janin plot - Janin plot for residues 5 to 10 of AdK, sampled from the AdK test trajectory + Janin plot for all residues of AdK, sampled from the AdK test trajectory (XTC). The contours in the background are the "allowed region" and the "marginally allowed" regions for all possible residues. +.. Note:: + + The Janin analysis is prone to errors if the topology contains duplicate or + missing atoms (e.g. atoms with `altloc` or incomplete residues). If the + topology has as an `altloc` attribute, you must specify only one `altloc` + for the atoms with more than one (``"protein and not altloc B"``). + + Furthermore, many residues do not have a :math:`\chi_2` dihedral and if the + selections of residues is not carefully filtered to only include those + residues with *both* sidechain dihedrals then a :exc:`ValueError` with the + message *Too many or too few atoms selected* is raised. + + +Reference plots +~~~~~~~~~~~~~~~ + +Reference plots can be added to the axes for both the Ramachandran and Janin +classes using the kwarg ``ref=True``. The Ramachandran reference data +(:data:`~MDAnalysis.analysis.data.filenames.Rama_ref`) and Janin reference data +(:data:`~MDAnalysis.analysis.data.filenames.Janin_ref`) were made using data +obtained from a large selection of 500 PDB files, and were analyzed using these +classes [Mull2018]_. The allowed and marginally allowed regions of the +Ramachandran reference plot have cutoffs set to include 90% and 99% of the data +points, and the Janin reference plot has cutoffs for 90% and 98% of the data +points. The list of PDB files used for the reference plots was taken from +[Lovell2003]_ and information about general Janin regions was taken from +[Janin1978]_. + + Analysis Classes ---------------- @@ -168,6 +201,19 @@ class which takes a list of `atomgroups`, this class only needs a list of References ---------- + +.. [Ramachandran1963] G. Ramachandran, C. Ramakrishnan, and + V. Sasisekharan. (1963) Stereochemistry of polypeptide chain + configurations. *Journal of Molecular Biology*, 7(1):95 – 99. doi: + `10.1016/S0022-2836(63)80023-6 + `_ + +.. [Janin1978] Joël Janin, Shoshanna Wodak, Michael Levitt, and Bernard + Maigret. (1978). "Conformation of amino acid side-chains in + proteins". *Journal of Molecular Biology* 125(3): 357-386. doi: + `10.1016/0022-2836(78)90408-4 + `_ + .. [Lovell2003] Simon C. Lovell, Ian W. Davis, W. Bryan Arendall III, Paul I. W. de Bakker, J. Michael Word, Michael G. Prisant, Jane S. Richardson, and David C. Richardson (2003). "Structure validation by @@ -175,10 +221,10 @@ class which takes a list of `atomgroups`, this class only needs a list of :math:`C_{\beta}` deviation". *Proteins* 50(3): 437-450. doi: `10.1002/prot.10286 `_ -.. [Janin1978] Joël Janin, Shoshanna Wodak, Michael Levitt, and Bernard - Maigret. (1978). "Conformation of amino acid side-chains in - proteins". *Journal of Molecular Biology* 125(3): 357-386. doi: - `10.1016/0022-2836(78)90408-4 `_ +.. [Mull2018] H. Mull and O. Beckstein. SPIDAL Summer REU 2018 Dihedral + Analysis in MDAnalysis. Technical report, Arizona State University, 8 + 2018. doi: `10.6084/m9.figshare.6957296 + `_ """ import numpy as np @@ -244,7 +290,8 @@ def _conclude(self): self.angles = np.rad2deg(np.array(self.angles)) class Ramachandran(AnalysisBase): - """Calculate :math:`\phi` and :math:`\psi` dihedral angles of selected residues. + r"""Calculate :math:`\phi` and :math:`\psi` dihedral angles of selected + residues. :math:`\phi` and :math:`\psi` angles will be calculated for each residue corresponding to `atomgroup` for each time step in the trajectory. A @@ -263,25 +310,25 @@ class Ramachandran(AnalysisBase): ca_name: str (optional) name for the alpha-carbon atom check_protein: bool (optional) - whether to raise an error if the provided atomgroup is not a + whether to raise an error if the provided atomgroup is not a subset of protein atoms Example ------- - For standard proteins, the default arguments will suffice to run a + For standard proteins, the default arguments will suffice to run a Ramachandran analysis:: r = Ramachandran(u.select_atoms('protein')).run() - For proteins with non-standard residues, or for calculating dihedral - angles for other linear polymers, you can switch off the protein checking - and provide your own atom names in place of the typical peptide backbone + For proteins with non-standard residues, or for calculating dihedral + angles for other linear polymers, you can switch off the protein checking + and provide your own atom names in place of the typical peptide backbone atoms:: r = Ramachandran(u.atoms, c_name='CX', n_name='NT', ca_name='S', check_protein=False).run() - The above analysis will calculate angles from a "phi" selection of + The above analysis will calculate angles from a "phi" selection of CX'-NT-S-CX and "psi" selections of NT-S-CX-NT'. Raises @@ -292,13 +339,14 @@ class Ramachandran(AnalysisBase): Note ---- - If ``check_protein`` is ``True`` and the residue selection is beyond - the scope of the protein and, then an error will be raised. + If ``check_protein`` is ``True`` and the residue selection is beyond + the scope of the protein and, then an error will be raised. If the residue selection includes the first or last residue, then a warning will be raised and they will be removed from the list of residues, but the analysis will still run. If a :math:`\phi` or :math:`\psi` selection cannot be made, that residue will be removed from the analysis. + .. versionchanged:: 1.0.0 added c_name, n_name, ca_name, and check_protein keyword arguments """ @@ -337,9 +385,9 @@ def __init__(self, atomgroup, c_name='C', n_name='N', ca_name='CA', # find n, c, ca keep_prev = [sum(r.atoms.names==c_name)==1 for r in prev] rnames = [n_name, c_name, ca_name] - keep_res = [all(sum(r.atoms.names==n)==1 for n in rnames) + keep_res = [all(sum(r.atoms.names == n) == 1 for n in rnames) for r in residues] - keep_next = [sum(r.atoms.names==n_name)==1 for r in nxt] + keep_next = [sum(r.atoms.names == n_name) == 1 for r in nxt] # alright we'll keep these keep = np.array(keep_prev) & np.array(keep_res) & np.array(keep_next) @@ -372,8 +420,10 @@ def _conclude(self): self.angles = np.rad2deg(np.array(self.angles)) def plot(self, ax=None, ref=False, **kwargs): - """Plots data into standard ramachandran plot. Each time step in - :attr:`Ramachandran.angles` is plotted onto the same graph. + """Plots data into standard Ramachandran plot. + + Each time step in :attr:`Ramachandran.angles` is plotted onto + the same graph. Parameters ---------- @@ -385,6 +435,9 @@ def plot(self, ax=None, ref=False, **kwargs): Adds a general Ramachandran plot which shows allowed and marginally allowed regions + kwargs : optional + All other kwargs are passed to :func:`matplotlib.pyplot.scatter`. + Returns ------- ax : :class:`matplotlib.axes.Axes` @@ -397,8 +450,13 @@ def plot(self, ax=None, ref=False, **kwargs): ax.axhline(0, color='k', lw=1) ax.axvline(0, color='k', lw=1) ax.set(xticks=range(-180, 181, 60), yticks=range(-180, 181, 60), - xlabel=r"$\phi$ (deg)", ylabel=r"$\psi$ (deg)") - if ref == True: + xlabel=r"$\phi$", ylabel=r"$\psi$") + degree_formatter = plt.matplotlib.ticker.StrMethodFormatter( + r"{x:g}$\degree$") + ax.xaxis.set_major_formatter(degree_formatter) + ax.yaxis.set_major_formatter(degree_formatter) + + if ref: X, Y = np.meshgrid(np.arange(-180, 180, 4), np.arange(-180, 180, 4)) levels = [1, 17, 15000] @@ -410,7 +468,8 @@ def plot(self, ax=None, ref=False, **kwargs): class Janin(Ramachandran): - """Calculate :math:`\chi_1` and :math:`\chi_2` dihedral angles of selected residues. + r"""Calculate :math:`\chi_1` and :math:`\chi_2` dihedral angles of + selected residues. :math:`\chi_1` and :math:`\chi_2` angles will be calculated for each residue corresponding to `atomgroup` for each time step in the trajectory. A @@ -420,46 +479,66 @@ class Janin(Ramachandran): Note ---- If the residue selection is beyond the scope of the protein, then an error - will be raised. If the residue selection includes the residues ALA, CYS, - GLY, PRO, SER, THR, or VAL, then a warning will be raised and they will be - removed from the list of residues, but the analysis will still run. Some - topologies have altloc attribues which can add duplicate atoms to the - selection and must be removed. + will be raised. If the residue selection includes the residues ALA, CYS*, + GLY, PRO, SER, THR, or VAL (the default of the `select_remove` keyword + argument) then a warning will be raised and they will be removed from the + list of residues, but the analysis will still run. Some topologies have + altloc attribues which can add duplicate atoms to the selection and must be + removed. """ - def __init__(self, atomgroup, **kwargs): - """Parameters + def __init__(self, atomgroup, + select_remove="resname ALA CYS* GLY PRO SER THR VAL", + select_protein="protein", + **kwargs): + r"""Parameters ---------- atomgroup : AtomGroup or ResidueGroup atoms for residues for which :math:`\chi_1` and :math:`\chi_2` are calculated + select_remove : str + selection string to remove residues that do not have :math:`chi_2` + angles + + select_protein : str + selection string to subselect protein-only residues from + `atomgroup` to check that only amino acids are selected; if you + have non-standard amino acids then adjust this selection to include + them + Raises ------ ValueError - If the selection of residues is not contained within the protein + if the final selection of residues is not contained within the + protein (as determined by + ``atomgroup.select_atoms(select_protein)``) ValueError - If not enough or too many atoms are found for a residue in the - selection, usually due to missing atoms or alternative locations + if not enough or too many atoms are found for a residue in the + selection, usually due to missing atoms or alternative locations, + or due to non-standard residues + + .. versionchanged:: 2.0.0 + `select_remove` and `select_protein` keywords were added """ super(Ramachandran, self).__init__( atomgroup.universe.trajectory, **kwargs) self.atomgroup = atomgroup residues = atomgroup.residues - protein = atomgroup.universe.select_atoms("protein").residues - remove = residues.atoms.select_atoms("resname ALA CYS GLY PRO SER" - " THR VAL").residues + protein = atomgroup.select_atoms(select_protein).residues + remove = residues.atoms.select_atoms(select_remove).residues if not residues.issubset(protein): raise ValueError("Found atoms outside of protein. Only atoms " - "inside of a 'protein' selection can be used to " - "calculate dihedrals.") + "inside of a protein " + f"(select_protein='{select_protein}') can be " + "used to calculate dihedrals.") elif len(remove) != 0: - warnings.warn("All ALA, CYS, GLY, PRO, SER, THR, and VAL residues" - " have been removed from the selection.") + warnings.warn(f"All residues selected with '{select_remove}' " + "have been removed from the selection.") residues = residues.difference(remove) self.ag1 = residues.atoms.select_atoms("name N") @@ -480,8 +559,10 @@ def _conclude(self): self.angles = (np.rad2deg(np.array(self.angles)) + 360) % 360 def plot(self, ax=None, ref=False, **kwargs): - """Plots data into standard Janin plot. Each time step in - :attr:`Janin.angles` is plotted onto the same graph. + """Plots data into standard Janin plot. + + Each time step in :attr:`Janin.angles` is plotted onto the + same graph. Parameters ---------- @@ -493,6 +574,9 @@ def plot(self, ax=None, ref=False, **kwargs): Adds a general Janin plot which shows allowed and marginally allowed regions + kwargs : optional + All other kwargs are passed to :func:`matplotlib.pyplot.scatter`. + Returns ------- ax : :class:`matplotlib.axes.Axes` @@ -505,8 +589,13 @@ def plot(self, ax=None, ref=False, **kwargs): ax.axhline(180, color='k', lw=1) ax.axvline(180, color='k', lw=1) ax.set(xticks=range(0, 361, 60), yticks=range(0, 361, 60), - xlabel=r"$\chi1$ (deg)", ylabel=r"$\chi2$ (deg)") - if ref == True: + xlabel=r"$\chi_1$", ylabel=r"$\chi_2$") + degree_formatter = plt.matplotlib.ticker.StrMethodFormatter( + r"{x:g}$\degree$") + ax.xaxis.set_major_formatter(degree_formatter) + ax.yaxis.set_major_formatter(degree_formatter) + + if ref: X, Y = np.meshgrid(np.arange(0, 360, 6), np.arange(0, 360, 6)) levels = [1, 6, 600] colors = ['#A1D4FF', '#35A1FF'] diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py index 8947f9ef7ca..2dc2cf1a06f 100644 --- a/package/MDAnalysis/analysis/msd.py +++ b/package/MDAnalysis/analysis/msd.py @@ -263,7 +263,7 @@ class EinsteinMSD(AnalysisBase): Dimensionality :math:`d` of the MSD. timeseries : :class:`numpy.ndarray` The averaged MSD over all the particles with respect to lag-time. - msd_per_particle : :class:`numpy.ndarray` + msds_by_particle : :class:`numpy.ndarray` The MSD of each individual particle with respect to lag-time. ag : :class:`AtomGroup` The :class:`AtomGroup` resulting from your selection diff --git a/package/MDAnalysis/coordinates/H5MD.py b/package/MDAnalysis/coordinates/H5MD.py index bf33eec169b..a49594282d5 100644 --- a/package/MDAnalysis/coordinates/H5MD.py +++ b/package/MDAnalysis/coordinates/H5MD.py @@ -31,7 +31,7 @@ interface of the HDF5 library to improve parallel reads and writes. The HDF5 library and `h5py`_ must be installed; otherwise, H5MD files -cannot be read by MDAnalysis. If `h5py`_ is not installed, a +cannot be read by MDAnalysis. If `h5py`_ is not installed, a :exc:`RuntimeError` is raised. Units @@ -94,10 +94,10 @@ Building parallel h5py and HDF5 on Linux ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -Building a working parallel HDF5/h5py/mpi4py environment can be +Building a working parallel HDF5/h5py/mpi4py environment can be challenging and is often specific to your local computing resources, e.g., the supercomputer that you're running on typically already has -its preferred MPI installation. As a starting point we provide +its preferred MPI installation. As a starting point we provide instructions that worked in a specific, fairly generic environment. These instructions successfully built parallel HDF5/h5py @@ -133,8 +133,8 @@ CC=mpicc HDF5_DIR=$HDF5_PATH python setup.py build python setup.py install -If you have questions or want to share how you managed to build -parallel hdf5/h5py/mpi4py please let everyone know on the +If you have questions or want to share how you managed to build +parallel hdf5/h5py/mpi4py please let everyone know on the `MDAnalysis forums`_. .. _`H5MD`: https://nongnu.org/h5md/index.html @@ -180,6 +180,9 @@ .. automethod:: H5MDReader._reopen +.. autoclass:: H5PYPicklable + :members: + """ import numpy as np @@ -190,6 +193,15 @@ import h5py except ImportError: HAS_H5PY = False + + # Allow building documentation even if h5py is not installed + import types + + class MockH5pyFile: + pass + h5py = types.ModuleType("h5py") + h5py.File = MockH5pyFile + else: HAS_H5PY = True @@ -228,23 +240,23 @@ class H5MDReader(base.ReaderBase): See `h5md documentation `_ for a detailed overview of the H5MD file format. - - The reader attempts to convert units in the trajectory file to - the standard MDAnalysis units (:mod:`MDAnalysis.units`) if + + The reader attempts to convert units in the trajectory file to + the standard MDAnalysis units (:mod:`MDAnalysis.units`) if `convert_units` is set to ``True``. - + Additional data in the *observables* group of the H5MD file are loaded into the :attr:`Timestep.data` dictionary. - + Only 3D-periodic boxes or no periodicity are supported; for no periodicity, :attr:`Timestep.dimensions` will return ``None``. - + Although H5MD can store varying numbers of particles per time step as produced by, e.g., GCMC simulations, MDAnalysis can currently only process a fixed number of particles per step. If the number of particles changes a :exc:`ValueError` is raised. - The :class:`H5MDReader` reads .h5md files with the following + The :class:`H5MDReader` reads .h5md files with the following HDF5 hierarchy: .. code-block:: text @@ -277,21 +289,21 @@ class H5MDReader(base.ReaderBase): \-- (position) \-- [step] , gives frame \-- [time] , gives time - +-- units + +-- unit \-- [value] , gives numpy arrary of positions with shape (n_atoms, 3) +-- unit \-- (velocity) \-- [step] , gives frame \-- [time] , gives time - +-- units + +-- unit \-- [value] , gives numpy arrary of velocities with shape (n_atoms, 3) +-- unit \-- (force) \-- [step] , gives frame \-- [time] , gives time - +-- units + +-- unit \-- [value] , gives numpy arrary of forces with shape (n_atoms, 3) +-- unit @@ -557,13 +569,15 @@ def open_trajectory(self): if self._comm is not None: # can only pass comm argument to h5py.File if driver='mpio' assert self._driver == 'mpio' - self._file = h5py.File(self.filename, 'r', # pragma: no cover - driver=self._driver, - comm=self._comm) + self._file = H5PYPicklable(name=self.filename, # pragma: no cover + mode='r', + driver=self._driver, + comm=self._comm) elif self._driver is not None: - self._file = h5py.File(self.filename, 'r', driver=self._driver) + self._file = H5PYPicklable(name=self.filename, mode='r', + driver=self._driver) else: - self._file = h5py.File(self.filename, 'r') + self._file = H5PYPicklable(name=self.filename, mode='r') # pulls first key out of 'particles' # allows for arbitrary name of group1 in 'particles' self._particle_group = self._file['particles'][ @@ -727,3 +741,88 @@ def has_forces(self): @has_forces.setter def has_forces(self, value: bool): self._has['force'] = value + + def __getstate__(self): + state = self.__dict__.copy() + del state['_particle_group'] + return state + + def __setstate__(self, state): + self.__dict__ = state + self._particle_group = self._file['particles'][ + list(self._file['particles'])[0]] + self[self.ts.frame] + + +class H5PYPicklable(h5py.File): + """H5PY file object (read-only) that can be pickled. + + This class provides a file-like object (as returned by + :class:`h5py.File`) that, + unlike standard Python file objects, + can be pickled. Only read mode is supported. + + When the file is pickled, filename, mode, driver, and comm of + :class:`h5py.File` in the file are saved. On unpickling, the file + is opened by filename, mode, driver. This means that for a successful + unpickle, the original file still has to be accessible with its filename. + + Parameters + ---------- + filename : str or file-like + a filename given a text or byte string. + driver : str (optional) + H5PY file driver used to open H5MD file + + Example + ------- + :: + + f = H5PYPicklable('filename', 'r') + print(f['particles/trajectory/position/value'][0]) + f.close() + + can also be used as context manager:: + + with H5PYPicklable('filename', 'r'): + print(f['particles/trajectory/position/value'][0]) + + Note + ---- + Pickling of an `h5py.File` opened with `driver="mpio"` and an MPI + communicator is currently not supported + + See Also + --------- + :class:`MDAnalysis.lib.picklable_file_io.FileIOPicklable` + :class:`MDAnalysis.lib.picklable_file_io.BufferIOPicklable` + :class:`MDAnalysis.lib.picklable_file_io.TextIOPicklable` + :class:`MDAnalysis.lib.picklable_file_io.GzipPicklable` + :class:`MDAnalysis.lib.picklable_file_io.BZ2Picklable` + + + .. versionadded:: 2.0.0 + """ + + def __getstate__(self): + driver = self.driver + # Current issues: Need a way to retrieve MPI communicator object + # from self and pickle MPI.Comm object. Parallel driver is excluded + # from test because h5py calls for an MPI configuration when driver is + # 'mpio', so this will need to be patched in the test function. + if driver == 'mpio': # pragma: no cover + raise TypeError("Parallel pickling of `h5py.File` with" # pragma: no cover + " 'mpio' driver is currently not supported.") + + return {'name': self.filename, + 'mode': self.mode, + 'driver': driver} + + def __setstate__(self, state): + self.__init__(name=state['name'], + mode=state['mode'], + driver=state['driver']) + + def __getnewargs__(self): + """Override the h5py getnewargs to skip its error message""" + return () diff --git a/package/MDAnalysis/coordinates/chemfiles.py b/package/MDAnalysis/coordinates/chemfiles.py index 8538c8fc27c..bd80f674790 100644 --- a/package/MDAnalysis/coordinates/chemfiles.py +++ b/package/MDAnalysis/coordinates/chemfiles.py @@ -50,11 +50,11 @@ HAS_CHEMFILES = False # Allow building documentation even if chemfiles is not installed - import imp + import types class MockTrajectory: pass - chemfiles = imp.new_module("chemfiles") + chemfiles = types.ModuleType("chemfiles") chemfiles.Trajectory = MockTrajectory else: HAS_CHEMFILES = True diff --git a/package/MDAnalysis/lib/NeighborSearch.py b/package/MDAnalysis/lib/NeighborSearch.py index 2bee44a3135..f4693047465 100644 --- a/package/MDAnalysis/lib/NeighborSearch.py +++ b/package/MDAnalysis/lib/NeighborSearch.py @@ -32,18 +32,16 @@ from MDAnalysis.lib.distances import capped_distance from MDAnalysis.lib.util import unique_int_1d -from MDAnalysis.core.groups import AtomGroup, Atom - class AtomNeighborSearch(object): """This class can be used to find all atoms/residues/segments within the radius of a given query position. - For the neighbor search, this class uses the BioPython KDTree and its - wrapper PeriodicKDTree for non-periodic and periodic systems, respectively. + For the neighbor search, this class is a wrapper around + :class:`~MDAnalysis.lib.distances.capped_distance`. """ - def __init__(self, atom_group, box=None, bucket_size=10): + def __init__(self, atom_group, box=None): """ Parameters @@ -55,16 +53,10 @@ def __init__(self, atom_group, box=None, bucket_size=10): :attr:`MDAnalysis.trajectory.base.Timestep.dimensions` when periodic boundary conditions should be taken into account for the calculation of contacts. - bucket_size : int - Number of entries in leafs of the KDTree. If you suffer poor - performance you can play around with this number. Increasing the - `bucket_size` will speed up the construction of the KDTree but - slow down the search. """ self.atom_group = atom_group self._u = atom_group.universe self._box = box - #self.kdtree = PeriodicKDTree(box=box, leafsize=bucket_size) def search(self, atoms, radius, level='A'): """ @@ -73,21 +65,37 @@ def search(self, atoms, radius, level='A'): Parameters ---------- - atoms : AtomGroup, MDAnalysis.core.groups.Atom - list of atoms + atoms : AtomGroup, MDAnalysis.core.groups.AtomGroup + AtomGroup object radius : float Radius for search in Angstrom. level : str char (A, R, S). Return atoms(A), residues(R) or segments(S) within *radius* of *atoms*. + + Returns + ------- + AtomGroup : :class:`~MDAnalysis.core.groups.AtomGroup` + When ``level='A'``, AtomGroup is being returned. + ResidueGroup : :class:`~MDAnalysis.core.groups.ResidueGroup` + When ``level='R'``, ResidueGroup is being returned. + SegmentGroup : :class:`~MDAnalysis.core.groups.SegmentGroup` + When ``level='S'``, SegmentGroup is being returned. + + + .. versionchanged:: 2.0.0 + Now returns :class:`AtomGroup` (when empty this is now an empty + :class:`AtomGroup` instead of an empty list), :class:`ResidueGroup`, + or a :class:`SegmentGroup` """ unique_idx = [] - if isinstance(atoms, Atom): - positions = atoms.position.reshape(1, 3) - else: - positions = atoms.positions - - pairs = capped_distance(positions, self.atom_group.positions, + try: + # For atom groups, take the positions attribute + position = atoms.positions + except AttributeError: + # For atom, take the position attribute + position = atoms.position + pairs = capped_distance(position, self.atom_group.positions, radius, box=self._box, return_distances=False) if pairs.size > 0: @@ -106,15 +114,12 @@ def _index2level(self, indices, level): char (A, R, S). Return atoms(A), residues(R) or segments(S) within *radius* of *atoms*. """ - n_atom_list = self.atom_group[indices] + atomgroup = self.atom_group[indices] if level == 'A': - if not n_atom_list: - return [] - else: - return n_atom_list + return atomgroup elif level == 'R': - return list({a.residue for a in n_atom_list}) + return atomgroup.residues elif level == 'S': - return list(set([a.segment for a in n_atom_list])) + return atomgroup.segments else: raise NotImplementedError('{0}: level not implemented'.format(level)) diff --git a/package/MDAnalysis/lib/log.py b/package/MDAnalysis/lib/log.py index a063878bd57..d3a6782aa12 100644 --- a/package/MDAnalysis/lib/log.py +++ b/package/MDAnalysis/lib/log.py @@ -176,57 +176,161 @@ def emit(self, record): pass -def echo(s='', replace=False, newline=True): - r"""Simple string output that immediately prints to the console. +class ProgressBar(tqdm): + r"""Display a visual progress bar and time estimate. + + The :class:`ProgressBar` decorates an iterable object, returning an + iterator which acts exactly like the original iterable, but prints a + dynamically updating progressbar every time a value is requested. See the + example below for how to use it when iterating over the frames of a + trajectory. - The string `s` is modified according to the keyword arguments and then - printed to :const:`sys.stderr`, which is immediately flushed. Parameters ---------- - s : str - The string to output. - replace : bool - If ``True``, the string will be printed on top of the current line. In - practice, ``\r`` is added at the beginning of the string. - newline : bool - If ``True``, a newline is added at the end of the string. - - """ - if replace: - s = '\r' + s - if newline: - end='\n' - else: - end='' - print(s, file=sys.stderr, end=end) - sys.stderr.flush() - - -def _new_format(template, variables): - """Format a string that follows the {}-based syntax. - """ - return template.format(**variables) - + iterable : iterable, optional + Iterable to decorate with a progressbar. + Leave blank to manually manage the updates. + verbose : bool, optional + If ``True`` (the default) then show the progress bar, *unless* the + `disable` keyword is set to ``True`` (`disable` takes precedence over + `verbose`). If `verbose` is set to ``None`` then the progress bar is + displayed (like ``True``), *unless* this is a non-TTY output device + (see `disable`). + desc : str, optional + Prefix for the progressbar. + total : int or float, optional + The number of expected iterations. If unspecified, + ``len(iterable)`` is used if possible. If ``float("inf")`` or as a last + resort, only basic progress statistics are displayed + (no ETA, no progressbar). + leave : bool, optional + If [default: ``True``], keeps all traces of the progressbar + upon termination of iteration. + If ``None``, will leave only if `position` is 0. + file : :class:`io.TextIOWrapper` or :class:`io.StringIO`, optional + Specifies where to output the progress messages (default: + :data:`sys.stderr`). Uses :meth:`file.write` and :meth:`file.flush` + methods. For encoding, see `write_bytes`. + ncols : int, optional + The width of the entire output message. If specified, + dynamically resizes the progressbar to stay within this bound. + If unspecified, attempts to use environment width. The + fallback is a meter width of 10 and no limit for the counter and + statistics. If 0, will not print any meter (only stats). + mininterval : float, optional + Minimum progress display update interval [default: 0.1] seconds. + maxinterval : float, optional + Maximum progress display update interval [default: 10] seconds. + Automatically adjusts `miniters` to correspond to `mininterval` + after long display update lag. Only works if `dynamic_miniters` + or monitor thread is enabled. + miniters : int or float, optional + Minimum progress display update interval, in iterations. + If 0 and `dynamic_miniters`, will automatically adjust to equal + `mininterval` (more CPU efficient, good for tight loops). + If > 0, will skip display of specified number of iterations. + Tweak this and `mininterval` to get very efficient loops. + If your progress is erratic with both fast and slow iterations + (network, skipping items, etc) you should set miniters=1. + ascii : bool or str, optional + If unspecified or ``False``, use unicode (smooth blocks) to fill + the meter. The fallback is to use ASCII characters " 123456789#". + disable : bool, optional + Whether to disable the entire progressbar wrapper + [default: ``False``]. If set to None, disable on non-TTY. + unit : str, optional + String that will be used to define the unit of each iteration + [default: it]. + unit_scale : bool or int or float, optional + If 1 or True, the number of iterations will be reduced/scaled + automatically and a metric prefix following the + International System of Units standard will be added + (kilo, mega, etc.) [default: ``False``]. If any other non-zero + number, will scale `total` and `n`. + dynamic_ncols : bool, optional + If set, constantly alters `ncols` and `nrows` to the + environment (allowing for window resizes) [default: ``False``]. + smoothing : float, optional + Exponential moving average smoothing factor for speed estimates + (ignored in GUI mode). Ranges from 0 (average speed) to 1 + (current/instantaneous speed) [default: 0.3]. + bar_format : str, optional + Specify a custom bar string formatting. May impact performance. + [default: ``'{l_bar}{bar}{r_bar}'``], where ``l_bar='{desc}: + {percentage:3.0f}%|'`` and ``r_bar='| {n_fmt}/{total_fmt} + [{elapsed}<{remaining}, {rate_fmt}{postfix}]'`` + + Possible vars: l_bar, bar, r_bar, n, n_fmt, total, total_fmt, + percentage, elapsed, elapsed_s, ncols, nrows, desc, unit, + rate, rate_fmt, rate_noinv, rate_noinv_fmt, + rate_inv, rate_inv_fmt, postfix, unit_divisor, + remaining, remaining_s. + + Note that a trailing ": " is automatically removed after {desc} + if the latter is empty. + initial : int or float, optional + The initial counter value. Useful when restarting a progress bar + [default: 0]. If using :class:`float`, consider specifying ``{n:.3f}`` + or similar in `bar_format`, or specifying `unit_scale`. + position : int, optional + Specify the line offset to print this bar (starting from 0) + Automatic if unspecified. + Useful to manage multiple bars at once (e.g., from threads). + postfix : dict or \*, optional + Specify additional stats to display at the end of the bar. + Calls ``set_postfix(**postfix)`` if possible (:class:`dict`). + unit_divisor : float, optional + [default: 1000], ignored unless `unit_scale` is ``True``. + write_bytes : bool, optional + If (default: ``None``) and `file` is unspecified, + bytes will be written in Python 2. If `True` will also write + bytes. In all other cases will default to unicode. + lock_args : tuple, optional + Passed to `refresh` for intermediate output + (initialisation, iterating, and updating). + nrows : int, optional + The screen height. If specified, hides nested bars outside this + bound. If unspecified, attempts to use environment height. + The fallback is 20. + + Returns + ------- + out : decorated iterator. + + Example + ------- + To get a progress bar when analyzing a trajectory:: + + from MDAnalysis.lib.log import ProgressBar + + ... + + for ts in ProgressBar(u.trajectory): + # perform analysis + + + will produce something similar to :: + + 30%|███████████ | 3/10 [00:13<00:30, 4.42s/it] + + in a terminal or Jupyter notebook. + + + See Also + -------- + The :class:`ProgressBar` is derived from :class:`tqdm.auto.tqdm`; see the + `tqdm documentation`_ for further details on how to use it. + + + + .. _`tqdm documentation`: https://tqdm.github.io/ -def _legacy_format(template, variables): - """Format a string that follows the %-based syntax. """ - return template % variables - - -def _guess_string_format(template): - """Guess if the template follow {}-based or %-based syntax. - """ - match = re.search(r'%\((step|numsteps|percentage)\)', template) - if match is None: - return _new_format - else: - return _legacy_format - - -class ProgressBar(tqdm): def __init__(self, *args, **kwargs): + """""" + # ^^^^ keep the empty doc string to avoid Sphinx doc errors with the + # original doc string from tqdm.auto.tqdm verbose = kwargs.pop('verbose', True) # disable: Whether to disable the entire progressbar wrapper [default: False]. # If set to None, disable on non-TTY. diff --git a/package/MDAnalysis/topology/HoomdXMLParser.py b/package/MDAnalysis/topology/HoomdXMLParser.py index 5a0f5da1aae..b8e7baa0613 100644 --- a/package/MDAnalysis/topology/HoomdXMLParser.py +++ b/package/MDAnalysis/topology/HoomdXMLParser.py @@ -126,7 +126,6 @@ def parse(self, **kwargs): pass else: attrs[attrname] = attr(np.array(vals, dtype=dtype)) - for attrname, attr, in ( ('bond', Bonds), ('angle', Angles), @@ -142,10 +141,10 @@ def parse(self, **kwargs): vals = [] attrs[attrname] = attr(vals) - if not 'masses' in attrs: - attrs['masses'] = Masses(np.zeros(natoms)) - if not 'charges' in attrs: - attrs['charges'] = Charges(np.zeros(natoms, dtype=np.float32)) + if 'mass' not in attrs: + attrs['mass'] = Masses(np.zeros(natoms)) + if 'charge' not in attrs: + attrs['charge'] = Charges(np.zeros(natoms, dtype=np.float32)) attrs = list(attrs.values()) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 4e27f81012d..57871292ff3 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -40,35 +40,39 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): ''' - Rotates the trajectory by a given angle on a given axis. The axis is defined by + Rotates the trajectory by a given angle on a given axis. The axis is defined by the user, combining the direction vector and a point. This point can be the center - of geometry or the center of mass of a user defined AtomGroup, or an array defining + of geometry or the center of mass of a user defined AtomGroup, or an array defining custom coordinates. - + Examples -------- - - e.g. rotate the coordinates by 90 degrees on a axis formed by the [0,0,1] vector and + + e.g. rotate the coordinates by 90 degrees on a axis formed by the [0,0,1] vector and the center of geometry of a given AtomGroup: - + .. code-block:: python - + + from MDAnalysis import transformations + ts = u.trajectory.ts angle = 90 - ag = u.atoms() + ag = u.atoms d = [0,0,1] - rotated = MDAnalysis.transformations.rotate(angle, direction=d, ag=ag)(ts) - + rotated = transformations.rotate.rotateby(angle, direction=d, ag=ag)(ts) + e.g. rotate the coordinates by a custom axis: - + .. code-block:: python + from MDAnalysis import transformations + ts = u.trajectory.ts angle = 90 p = [1,2,3] d = [0,0,1] - rotated = MDAnalysis.transformations.rotate(angle, direction=d, point=point)(ts) - + rotated = transformations.rotate.rotateby(angle, direction=d, point=p)(ts) + Parameters ---------- angle: float @@ -97,7 +101,7 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): Returns ------- MDAnalysis.coordinates.base.Timestep - + Warning ------- Wrapping/unwrapping the trajectory or performing PBC corrections may not be possible @@ -132,7 +136,7 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): center_method = partial(atoms.center, weights, pbc=wrap) else: raise ValueError('A point or an AtomGroup must be specified') - + def wrapped(ts): if point is None: position = center_method() @@ -143,8 +147,8 @@ def wrapped(ts): translation = matrix[:3, 3] ts.positions= np.dot(ts.positions, rotation) ts.positions += translation - + return ts - + return wrapped - + diff --git a/package/doc/sphinx/source/documentation_pages/coordinates/pickle_readers.rst b/package/doc/sphinx/source/documentation_pages/coordinates/pickle_readers.rst index 0866590297f..fbcc881ab0c 100644 --- a/package/doc/sphinx/source/documentation_pages/coordinates/pickle_readers.rst +++ b/package/doc/sphinx/source/documentation_pages/coordinates/pickle_readers.rst @@ -1,4 +1,4 @@ -.. Contains the formatted docstrings for the serialization of universe located +.. Contains the formatted docstrings for the serialization of universe located .. mainly in 'MDAnalysis/libs/pickle_file_io.py' .. _serialization: @@ -12,7 +12,7 @@ and what developers should do to serialize a new reader. To make sure every Trajectory reader can be successfully serialized, we implement picklable I/O classes (see :ref:`implemented-fileio`). -When the file is pickled, filename and other necessary attributes of the open +When the file is pickled, filename and other necessary attributes of the open file handle are saved. On unpickling, the file is opened by filename. This means that for a successful unpickle, the original file still has to be accessible with its filename. To retain the current frame of the trajectory, @@ -25,12 +25,12 @@ How to serialize a new reader File Access ^^^^^^^^^^^ -If the new reader uses :func:`util.anyopen()` +If the new reader uses :func:`util.anyopen()` (e.g. :class:`MDAnalysis.coordinates.PDB.PDBReader`), the reading handler can be pickled without modification. If the new reader uses I/O classes from other package (e.g. :class:`MDAnalysis.coordinates.GSD.GSDReader`), -and cannot be pickled natively, create a new picklable class inherited from +and cannot be pickled natively, create a new picklable class inherited from the file class in that package (e.g. :class:`MDAnalysis.coordinates.GSD.GSDPicklable`), adding :func:`__getstate__`, @@ -40,10 +40,10 @@ to allow file handler serialization. To seek or not to seek ^^^^^^^^^^^^^^^^^^^^^^ -Some I/O classes support :func:`seek` and :func:`tell` functions to allow the file +Some I/O classes support :func:`seek` and :func:`tell` functions to allow the file to be pickled with an offset. It is normally not needed for MDAnalysis with random access. But if error occurs during testing, find a way to make the offset work. -Maybe this I/O class supports frame indexing? Maybe the file handler inside this I/O +Maybe this I/O class supports frame indexing? Maybe the file handler inside this I/O class supports offset? For example, in :class:`MDAnalysis.coordinates.TRZ.TRZReader`, @@ -99,3 +99,4 @@ Currently implemented picklable IO Formats * :class:`MDAnalysis.coordinates.GSD.GSDPicklable` * :class:`MDAnalysis.coordinates.TRJ.NCDFPicklable` * :class:`MDAnalysis.coordinates.chemfiles.ChemfilesPicklable` +* :class:`MDAnalysis.coordinates.H5MD.H5PYPicklable` diff --git a/package/doc/sphinx/source/images/janin_demo_plot.png b/package/doc/sphinx/source/images/janin_demo_plot.png index f64cf64265d..b531fbc306a 100644 Binary files a/package/doc/sphinx/source/images/janin_demo_plot.png and b/package/doc/sphinx/source/images/janin_demo_plot.png differ diff --git a/package/doc/sphinx/source/images/rama_demo_plot.png b/package/doc/sphinx/source/images/rama_demo_plot.png index c81ef10e328..5007bd7a314 100644 Binary files a/package/doc/sphinx/source/images/rama_demo_plot.png and b/package/doc/sphinx/source/images/rama_demo_plot.png differ diff --git a/testsuite/MDAnalysisTests/analysis/test_dihedrals.py b/testsuite/MDAnalysisTests/analysis/test_dihedrals.py index d05963a3e71..07419830286 100644 --- a/testsuite/MDAnalysisTests/analysis/test_dihedrals.py +++ b/testsuite/MDAnalysisTests/analysis/test_dihedrals.py @@ -26,9 +26,10 @@ import pytest import MDAnalysis as mda -from MDAnalysisTests.datafiles import (GRO, XTC, DihedralArray, DihedralsArray, - RamaArray, GLYRamaArray, JaninArray, - LYSJaninArray, PDB_rama, PDB_janin) +from MDAnalysisTests.datafiles import (GRO, XTC, TPR, DihedralArray, + DihedralsArray, RamaArray, GLYRamaArray, + JaninArray, LYSJaninArray, PDB_rama, + PDB_janin) from MDAnalysis.analysis.dihedrals import Dihedral, Ramachandran, Janin @@ -106,7 +107,7 @@ def test_outside_protein_length(self, universe): with pytest.raises(ValueError): rama = Ramachandran(universe.select_atoms("resid 220"), check_protein=True).run() - + def test_outside_protein_unchecked(self, universe): rama = Ramachandran(universe.select_atoms("resid 220"), check_protein=False).run() @@ -133,15 +134,26 @@ class TestJanin(object): def universe(self): return mda.Universe(GRO, XTC) + @pytest.fixture() + def universe_tpr(self): + return mda.Universe(TPR, XTC) + @pytest.fixture() def janin_ref_array(self): return np.load(JaninArray) def test_janin(self, universe, janin_ref_array): - janin = Janin(universe.select_atoms("protein")).run() + self._test_janin(universe, janin_ref_array) + + def test_janin_tpr(self, universe_tpr, janin_ref_array): + """Test that CYSH are filtered (#2898)""" + self._test_janin(universe_tpr, janin_ref_array) + + def _test_janin(self, u, ref_array): + janin = Janin(u.select_atoms("protein")).run() # Test precision lowered to account for platform differences with osx - assert_almost_equal(janin.angles, janin_ref_array, 3, + assert_almost_equal(janin.angles, ref_array, 3, err_msg="error: dihedral angles should " "match test values") diff --git a/testsuite/MDAnalysisTests/lib/test_neighborsearch.py b/testsuite/MDAnalysisTests/lib/test_neighborsearch.py index e4b912fc6b3..a6f06b3d6e2 100644 --- a/testsuite/MDAnalysisTests/lib/test_neighborsearch.py +++ b/testsuite/MDAnalysisTests/lib/test_neighborsearch.py @@ -49,3 +49,14 @@ def test_search(universe): ns_res = ns.search(universe.atoms[20], 20) pns_res = pns.search(universe.atoms[20], 20) assert_equal(ns_res, pns_res) + + +def test_zero(universe): + """Check if empty atomgroup, residue, segments are returned""" + ns = NeighborSearch.AtomNeighborSearch(universe.atoms[:10]) + ns_res = ns.search(universe.atoms[20], 0.1, level='A') + assert ns_res == universe.atoms[[]] + ns_res = ns.search(universe.atoms[20], 0.1, level='R') + assert ns_res == universe.atoms[[]].residues + ns_res = ns.search(universe.atoms[20], 0.1, level='S') + assert ns_res == universe.atoms[[]].segments diff --git a/testsuite/MDAnalysisTests/topology/test_hoomdxml.py b/testsuite/MDAnalysisTests/topology/test_hoomdxml.py index cd48bef4886..018470647f4 100644 --- a/testsuite/MDAnalysisTests/topology/test_hoomdxml.py +++ b/testsuite/MDAnalysisTests/topology/test_hoomdxml.py @@ -20,6 +20,8 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # +from numpy.testing import assert_almost_equal + import MDAnalysis as mda from MDAnalysisTests.topology.base import ParserBase @@ -52,7 +54,7 @@ def test_angles(self, top): def test_dihedrals(self, top): assert len(top.dihedrals.values) == 576 assert isinstance(top.dihedrals.values[0], tuple) - + def test_impropers(self, top): assert len(top.impropers.values) == 0 @@ -73,3 +75,12 @@ def test_dihedrals_identity(self, top): for b in ((0, 1, 2, 3), (1, 2, 3, 4), (2, 3, 4, 5), (3, 4, 5, 6)): assert (b in vals) or (b[::-1] in vals) assert ((0, 250, 350, 450) not in vals) + + def test_read_masses(self, top): + assert_almost_equal(top.masses.values, 1.0) + + def test_read_charges(self, top): + # note: the example topology file contains 0 for all charges which + # is the same as the default so this test does not fully test + # reading of charges from the file (#2888) + assert_almost_equal(top.charges.values, 0.0) diff --git a/testsuite/MDAnalysisTests/utils/test_pickleio.py b/testsuite/MDAnalysisTests/utils/test_pickleio.py index b005f478c12..bbc7a5a5f23 100644 --- a/testsuite/MDAnalysisTests/utils/test_pickleio.py +++ b/testsuite/MDAnalysisTests/utils/test_pickleio.py @@ -50,6 +50,10 @@ from MDAnalysis.coordinates.chemfiles import ( ChemfilesPicklable ) +from MDAnalysis.coordinates.H5MD import HAS_H5PY +from MDAnalysis.coordinates.H5MD import ( + H5PYPicklable +) from MDAnalysis.tests.datafiles import ( PDB, @@ -58,7 +62,9 @@ MMTF_gz, GMS_ASYMOPT, GSD, - NCDF + NCDF, + TPR_xvf, + H5MD_xvf ) @@ -200,3 +206,19 @@ def test_Chemfiles_with_write_mode(tmpdir): with pytest.raises(ValueError, match=r"Only read mode"): chemfiles_io = ChemfilesPicklable(tmpdir.mkdir("xyz").join('t.xyz'), mode='w') + + +@pytest.mark.skipif(not HAS_H5PY, reason="h5py not installed") +def test_H5MD_pickle(): + h5md_io = H5PYPicklable(H5MD_xvf, 'r') + h5md_io_pickled = pickle.loads(pickle.dumps(h5md_io)) + assert_equal(h5md_io['particles/trajectory/position/value'][0], + h5md_io_pickled['particles/trajectory/position/value'][0]) + + +@pytest.mark.skipif(not HAS_H5PY, reason="h5py not installed") +def test_H5MD_pickle_with_driver(): + h5md_io = H5PYPicklable(H5MD_xvf, 'r', driver='core') + h5md_io_pickled = pickle.loads(pickle.dumps(h5md_io)) + assert_equal(h5md_io['particles/trajectory/position/value'][0], + h5md_io_pickled['particles/trajectory/position/value'][0])