Merge branch 'develop' into gsoc-1

.github/workflows/gh-ci.yaml

@@ -127,7 +127,7 @@ jobs:
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         file: coverage.xml
-        fail_ci_if_error: True
+        fail_ci_if_error: False
         verbose: True
 
 
@@ -259,22 +259,22 @@ jobs:
 
     - name: check_package_build
       run: |
-        DISTRIBUTION=$(ls -t1 dist/MDAnalysis-*.tar.gz | head -n1)
+        DISTRIBUTION=$(ls -t1 dist/mdanalysis-*.tar.gz | head -n1)
         test -n "${DISTRIBUTION}" || { echo "no distribution dist/MDAnalysis-*.tar.gz found"; exit 1; }
         twine check $DISTRIBUTION
 
     - name: check_testsuite_build
       run: |
-        DISTRIBUTION=$(ls -t1 dist/MDAnalysisTests-*.tar.gz | head -n1)
+        DISTRIBUTION=$(ls -t1 dist/mdanalysistests-*.tar.gz | head -n1)
         test -n "${DISTRIBUTION}" || { echo "no distribution dist/MDAnalysisTests-*.tar.gz found"; exit 1; }
         twine check $DISTRIBUTION
 
     - name: install sdist
       working-directory: ./dist
       run: |
         ls -a .
-        python -m pip install MDAnalysis-*.tar.gz
-        python -m pip install MDAnalysisTests-*.tar.gz
+        python -m pip install mdanalysis-*.tar.gz
+        python -m pip install mdanalysistests-*.tar.gz
 
     - name: run tests
       working-directory: ./dist

README.rst

@@ -102,19 +102,18 @@ Citation
 When using MDAnalysis in published work, please cite the following
 two papers:
 
-*   R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy,
-    M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski,
-    S. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis:
+*   R\. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy,
+    M\. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski,
+    S\. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis:
     A Python package for the rapid analysis of molecular
     dynamics simulations. In S. Benthall and S. Rostrup,
     editors, Proceedings of the 15th Python in Science
     Conference, pages 102-109, Austin, TX, 2016. SciPy.
-    doi:`10.25080/Majora-629e541a-00e`_    
-
+    doi: `10.25080/Majora-629e541a-00e`_    
 *   N. Michaud-Agrawal, E. J. Denning, T. B. Woolf,
     and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular
     Dynamics Simulations. *J. Comput. Chem.* **32** (2011), 2319--2327.
-    doi:`10.1002/jcc.21787`_
+    doi: `10.1002/jcc.21787`_
 
 For citations of included algorithms and sub-modules please see the references_.
 
@@ -194,4 +193,4 @@ For citations of included algorithms and sub-modules please see the references_.
 
 .. |discussions| image:: https://img.shields.io/github/discussions/MDAnalysis/MDAnalysis
    :alt: GitHub Discussions
-   :target: https://github.com/MDAnalysis/mdanalysis/discussions
+   :target: https://github.com/MDAnalysis/mdanalysis/discussions

package/AUTHORS

@@ -236,8 +236,13 @@ Chronological list of authors
   2024
   - Aditya Keshari
   - Philipp Stärk
+  - Kai Niklas Spauszus
+  - Sampurna Mukherjee
+  - Leon Wehrhan
+  - Valerij Talagayev
   - Luna Morrow
 
+
 External code
 -------------
 

package/CHANGELOG

@@ -37,12 +37,16 @@ Deprecations
 
 ??/??/?? IAlibay, HeetVekariya, marinegor, lilyminium, RMeli,
          ljwoods2, aditya292002, pstaerk, PicoCentauri, BFedder,
-         tyler.je.reddy
+         tyler.je.reddy, SampurnaM, leonwehrhan, kainszs, orionarcher
 
  * 2.8.0
 
 Fixes
- * Add support for TPR files produced by Gromacs 2024.1
+ * Fix doctest errors of analysis/pca.py related to rounding issues 
+   (Issue #3925, PR #4377)
+ * Convert openmm Quantity to raw value for KE and PE in OpenMMSimulationReader.
+ * Atomname methods can handle empty groups (Issue #2879, PR #4529)
+ * Add support for TPR files produced by Gromacs 2024.1 (PR #4523)
  * Remove mutable data from ``progressbar_kwargs`` argument in ``AnalysisBase.run()``
    (PR #4459)
  * Fix ChainReader `__repr__()` method when sub-reader is MemoryReader 
@@ -53,12 +57,17 @@ Fixes
    `DistanceMatrix.run(frames=...)` (PR #4433)
  * Fix doctest errors of lib/picklable_file_io.py (Issue #3925, PR #4371)
  * Fix deploy action to use the correct version of the pypi upload action.
+ * Fix groups.py doctests using sphinx directives (Issue #3925, PR #4374)
 
 Enhancements
+ * Added a tqdm progress bar for `MDAnalysis.analysis.pca.PCA.transform()`
+   (PR #4531)
  * Improved performance of PDBWriter (Issue #2785, PR #4472)
  * Added parsing of arbitrary columns of the LAMMPS dump parser. (Issue #3504)
  * Documented the r0 attribute in the `Contacts` class and added the 
    `n_initial_contacts` attribute, with documentation. (Issue #2604, PR #4415)
+ * Implement average structures with iterative algorithm from 
+   DOI 10.1021/acs.jpcb.7b11988. (Issue #2039, PR #4524)
 
 Changes
  * As per SPEC0 the minimum supported Python version has been raised

package/MDAnalysis/analysis/align.py

@@ -80,14 +80,14 @@
 
 Note that in this example translations have not been removed. In order
 to look at the pure rotation one needs to superimpose the centres of
-mass (or geometry) first:
+mass (or geometry) first::
 
    >>> rmsd(mobile.select_atoms('name CA').positions, ref.select_atoms('name CA').positions, center=True)
    21.892591663632704
 
 This has only done a translational superposition. If you want to also do a
 rotational superposition use the superposition keyword. This will calculate a
-minimized RMSD between the reference and mobile structure.
+minimized RMSD between the reference and mobile structure::
 
    >>> rmsd(mobile.select_atoms('name CA').positions, ref.select_atoms('name CA').positions, 
    ...      superposition=True)
@@ -174,6 +174,7 @@
 .. autoclass:: AlignTraj
 .. autoclass:: AverageStructure
 .. autofunction:: rotation_matrix
+.. autofunction:: iterative_average
 
 
 Helper functions
@@ -212,6 +213,8 @@
 from MDAnalysis.coordinates.memory import MemoryReader
 from MDAnalysis.lib.util import get_weights
 from MDAnalysis.lib.util import deprecate   # remove 3.0
+from MDAnalysis.lib.log import ProgressBar
+from ..due import due, Doi
 
 from .base import AnalysisBase
 
@@ -541,6 +544,125 @@ def alignto(mobile, reference, select=None, weights=None,
     return old_rmsd, new_rmsd
 
 
+@due.dcite(
+        Doi("10.1021/acs.jpcb.7b11988"),
+        description="Iterative Calculation of Opimal Reference",
+        path="MDAnalysis.analysis.align.iterative_average"
+)
+def iterative_average(
+    mobile, reference=None, select='all', weights=None, niter=100,
+    eps=1e-6, verbose=False, **kwargs
+):
+    """Iteratively calculate an optimal reference that is also the average
+    structure after an RMSD alignment.
+
+    The optimal reference is defined as average
+    structure of a trajectory, with the optimal reference used as input.
+    This function computes the optimal reference by using a starting
+    reference for the average structure, which is used as the reference
+    to calculate the average structure again. This is repeated until the
+    reference structure has converged. :footcite:p:`Linke2018`
+
+    Parameters
+    ----------
+    mobile : mda.Universe
+        Universe containing trajectory to be fitted to reference.
+    reference : mda.Universe (optional)
+        Universe containing the initial reference structure.
+    select : str or tuple or dict (optional)
+        Atom selection for fitting a substructue. Default is set to all.
+        Can be tuple or dict to define different selection strings for
+        mobile and target.
+    weights : str, array_like (optional)
+        Weights that can be used. If `None` use equal weights, if `'mass'`
+        use masses of ref as weights or give an array of arbitrary weights.
+    niter : int (optional)
+        Maximum number of iterations.
+    eps : float (optional)
+        RMSD distance at which reference and average are assumed to be
+        equal.
+    verbose : bool (optional)
+        Verbosity.
+    **kwargs : dict (optional)
+        AverageStructure kwargs.
+
+    Returns
+    -------
+    avg_struc : AverageStructure
+        AverageStructure result from the last iteration.
+
+    Example
+    -------
+    `iterative_average` can be used to obtain a :class:`MDAnalysis.Universe`
+    with the optimal reference structure.
+
+    ::
+
+        import MDAnalysis as mda
+        from MDAnalysis.analysis import align
+        from MDAnalysisTests.datafiles import PSF, DCD
+
+        u = mda.Universe(PSF, DCD)
+        av = align.iterative_average(u, u, verbose=True)
+
+        averaged_universe = av.results.universe
+
+    References
+    ----------
+
+    .. footbibliography::
+
+    .. versionadded:: 2.8.0
+    """
+    if not reference:
+        reference = mobile
+
+    select = rms.process_selection(select)
+    ref = mda.Merge(reference.select_atoms(*select['reference']))
+    sel_mobile = select['mobile'][0]
+
+    weights = get_weights(ref.atoms, weights)
+
+    drmsd = np.inf
+    for i in ProgressBar(range(niter)):
+        # found a converged structure
+        if drmsd < eps:
+            break
+
+        avg_struc = AverageStructure(
+            mobile, reference=ref, select={
+                'mobile': sel_mobile, 'reference': 'all'
+                },
+            weights=weights, **kwargs
+        ).run()
+        drmsd = rms.rmsd(ref.atoms.positions, avg_struc.results.positions,
+                         weights=weights)
+        ref = avg_struc.results.universe
+
+        if verbose:
+            logger.debug(
+                f"iterative_average(): i = {i}, "
+                f"rmsd-change = {drmsd:.5f}, "
+                f"ave-rmsd = {avg_struc.results.rmsd:.5f}"
+            )
+
+    else:
+        errmsg = (
+            "iterative_average(): Did not converge in "
+            f"{niter} iterations to DRMSD < {eps}. "
+            f"Final average RMSD = {avg_struc.results.rmsd:.5f}"
+        )
+        logger.error(errmsg)
+        raise RuntimeError(errmsg)
+
+    logger.info(
+        f"iterative_average(): Converged to DRMSD < {eps}. "
+        f"Final average RMSD = {avg_struc.results.rmsd:.5f}"
+    )
+
+    return avg_struc
+
+
 class AlignTraj(AnalysisBase):
     """RMS-align trajectory to a reference structure using a selection.
 

package/MDAnalysis/analysis/helix_analysis.py

@@ -121,7 +121,7 @@ def vector_of_best_fit(coordinates):
     """
     centered = coordinates - coordinates.mean(axis=0)
     Mt_M = np.matmul(centered.T, centered)
-    u, s, vh = np.linalg.linalg.svd(Mt_M)
+    u, s, vh = np.linalg.svd(Mt_M)
     vector = vh[0]
 
     # does vector face first local helix origin?

package/MDAnalysis/analysis/pca.py

@@ -118,6 +118,7 @@
 
 import numpy as np
 import scipy.integrate
+from tqdm.auto import tqdm
 
 from MDAnalysis import Universe
 from MDAnalysis.analysis.align import _fit_to
@@ -355,12 +356,12 @@ def n_components(self, n):
                 n = len(self._variance)
             self.results.variance = self._variance[:n]
             self.results.cumulated_variance = (np.cumsum(self._variance) /
-                                       np.sum(self._variance))[:n]
+                                               np.sum(self._variance))[:n]
             self.results.p_components = self._p_components[:, :n]
         self._n_components = n
 
     def transform(self, atomgroup, n_components=None, start=None, stop=None,
-                  step=None):
+                  step=None, verbose=False):
         """Apply the dimensionality reduction on a trajectory
 
         Parameters
@@ -382,6 +383,11 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None,
             Include every `step` frames in the PCA transform. If set to
             ``None`` (the default) then every frame is analyzed (i.e., same as
             ``step=1``).
+        verbose : bool, optional
+            ``verbose = True`` option displays a progress bar for the
+            iterations of transform. ``verbose = False`` disables the
+            progress bar, just returns the pca_space array when the
+            calculations are finished.
 
         Returns
         -------
@@ -391,14 +397,17 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None,
         .. versionchanged:: 0.19.0
            Transform now requires that :meth:`run` has been called before,
            otherwise a :exc:`ValueError` is raised.
+        .. versionchanged:: 2.8.0
+           Transform now has shows a tqdm progressbar, which can be toggled
+           on with ``verbose = True``, or off with ``verbose = False``
         """
         if not self._calculated:
             raise ValueError('Call run() on the PCA before using transform')
 
         if isinstance(atomgroup, Universe):
             atomgroup = atomgroup.atoms
 
-        if(self._n_atoms != atomgroup.n_atoms):
+        if self._n_atoms != atomgroup.n_atoms:
             raise ValueError('PCA has been fit for'
                              '{} atoms. Your atomgroup'
                              'has {} atoms'.format(self._n_atoms,
@@ -415,7 +424,10 @@ def transform(self, atomgroup, n_components=None, start=None, stop=None,
 
         dot = np.zeros((n_frames, dim))
 
-        for i, ts in enumerate(traj[start:stop:step]):
+        for i, ts in tqdm(enumerate(traj[start:stop:step]), disable=not verbose,
+                          total=len(traj[start:stop:step])
+                          ):
+
             xyz = atomgroup.positions.ravel() - self._xmean
             dot[i] = np.dot(xyz, self._p_components[:, :dim])
 
@@ -561,7 +573,7 @@ def project_single_frame(self, components=None, group=None, anchor=None):
             for res in group.residues:
                 # n_common is the number of pca atoms in a residue
                 n_common = pca_res_counts[np.where(
-                           pca_res_indices == res.resindex)][0]
+                            pca_res_indices == res.resindex)][0]
                 non_pca_atoms = np.append(non_pca_atoms,
                                           res.atoms.n_atoms - n_common)
             # index_extrapolate records the anchor number for each non-PCA atom
@@ -638,10 +650,10 @@ def rmsip(self, other, n_components=None):
             >>> first_interval = pca.PCA(u, select="backbone").run(start=0, stop=25)
             >>> second_interval = pca.PCA(u, select="backbone").run(start=25, stop=50)
             >>> last_interval = pca.PCA(u, select="backbone").run(start=75)
-            >>> first_interval.rmsip(second_interval, n_components=3)
-            0.38147609331128324
-            >>> first_interval.rmsip(last_interval, n_components=3)
-            0.17478244043688052
+            >>> round(first_interval.rmsip(second_interval, n_components=3), 6)
+            0.381476
+            >>> round(first_interval.rmsip(last_interval, n_components=3), 6)
+            0.174782
 
 
         See also
@@ -814,14 +826,14 @@ def rmsip(a, b, n_components=None):
         >>> first_interval = pca.PCA(u, select="backbone").run(start=0, stop=25)
         >>> second_interval = pca.PCA(u, select="backbone").run(start=25, stop=50)
         >>> last_interval = pca.PCA(u, select="backbone").run(start=75)
-        >>> pca.rmsip(first_interval.results.p_components.T,
+        >>> round(pca.rmsip(first_interval.results.p_components.T,
         ...           second_interval.results.p_components.T,
-        ...           n_components=3)
-        0.38147609331128324
-        >>> pca.rmsip(first_interval.results.p_components.T,
+        ...           n_components=3), 6)
+        0.381476
+        >>> round(pca.rmsip(first_interval.results.p_components.T,
         ...           last_interval.results.p_components.T,
-        ...           n_components=3)
-        0.17478244043688052
+        ...           n_components=3), 6)
+        0.174782
 
 
     .. versionadded:: 1.0.0