LinearDensity now working with residues/segments

package/CHANGELOG

@@ -18,6 +18,8 @@ The rules for this file:
  * 2.2.0
 
 Fixes
+  * Fixed LinearDensity not working with grouping='residues' or 'segments'
+    (Issue #3571, PR #3572)
   * When converting an AtomGroup to an RDKit molecule (PR #3044):
     - the atoms are now in the same order
     - the output of `atom.GetMonomerInfo().GetName()` now follows the

package/MDAnalysis/analysis/lineardensity.py

@@ -43,8 +43,9 @@ class LinearDensity(AnalysisBase):
     select : AtomGroup
           any atomgroup
     grouping : str {'atoms', 'residues', 'segments', 'fragments'}
-          Density profiles will be computed on the center of geometry
-          of a selected group of atoms
+          Density profiles will be computed either on the atom positions (in
+          the case of 'atoms') or on the center of mass of the specified
+          grouping unit ('residues', 'segments', or 'fragments').
     binsize : float
           Bin width in Angstrom used to build linear density
           histograms. Defines the resolution of the resulting density
@@ -79,6 +80,8 @@ class LinearDensity(AnalysisBase):
        ldens.run()
        print(ldens.results.x.pos)
 
+       # alternatively, specify grouping/binsize like this
+       ldens = LinearDensity(selection, grouping='residues', binsize=1.0)
 
     .. versionadded:: 0.14.0
 
@@ -96,6 +99,12 @@ class LinearDensity(AnalysisBase):
        Results are now instances of
        :class:`~MDAnalysis.core.analysis.Results` allowing access
        via key and attribute.
+
+    .. versionchanged:: 2.2.0
+       Fixed a bug that caused LinearDensity to fail if grouping="residues"
+       or grouping="segments" were set.
+       Residues, segments, and fragments will be analysed based on their centre
+       of mass, not centre of geometry as previously stated.
     """
 
     def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs):
@@ -141,18 +150,19 @@ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs):
         self.totalmass = None
 
     def _prepare(self):
-        # group must be a local variable, otherwise there will be
-        # issues with parallelization
-        group = getattr(self._ags[0], self.grouping)
-
         # Get masses and charges for the selection
-        try:  # in case it's not an atom
-            self.masses = np.array([elem.total_mass() for elem in group])
-            self.charges = np.array([elem.total_charge() for elem in group])
-        except AttributeError:  # much much faster for atoms
+        if self.grouping == "atoms":
             self.masses = self._ags[0].masses
             self.charges = self._ags[0].charges
 
+        elif self.grouping in ["residues", "segments", "fragments"]:
+            self.masses = self._ags[0].total_mass(compound=self.grouping)
+            self.charges = self._ags[0].total_charge(compound=self.grouping)
+
+        else:
+            raise AttributeError(
+                f"{self.grouping} is not a valid value for grouping.")
+
         self.totalmass = np.sum(self.masses)
 
     def _single_frame(self):
@@ -163,8 +173,8 @@ def _single_frame(self):
         if self.grouping == 'atoms':
             positions = self._ags[0].positions  # faster for atoms
         else:
-            # COM for res/frag/etc
-            positions = np.array([elem.centroid() for elem in self.group])
+            # Centre of mass for residues, segments, fragments
+            positions = self._ags[0].center_of_mass(compound=self.grouping)
 
         for dim in ['x', 'y', 'z']:
             idx = self.results[dim]['dim']
@@ -199,10 +209,20 @@ def _conclude(self):
             for key in ['pos', 'pos_std', 'char', 'char_std']:
                 self.results[dim][key] /= self.n_frames
             # Compute standard deviation for the error
-            self.results[dim]['pos_std'] = np.sqrt(self.results[dim][
-                'pos_std'] - np.square(self.results[dim]['pos']))
-            self.results[dim]['char_std'] = np.sqrt(self.results[dim][
-                'char_std'] - np.square(self.results[dim]['char']))
+            # For certain tests in testsuite, floating point imprecision
+            # can lead to negative radicands of tiny magnitude (yielding nan).
+            # radicand_pos and radicand_char are therefore calculated first and
+            # and negative values set to 0 before the square root
+            # is calculated.
+            radicand_pos = self.results[dim][
+                'pos_std'] - np.square(self.results[dim]['pos'])
+            radicand_pos[radicand_pos < 0] = 0
+            self.results[dim]['pos_std'] = np.sqrt(radicand_pos)
+
+            radicand_char = self.results[dim][
+                'char_std'] - np.square(self.results[dim]['char'])
+            radicand_char[radicand_char < 0] = 0
+            self.results[dim]['char_std'] = np.sqrt(radicand_char)
 
         for dim in ['x', 'y', 'z']:
             norm = k * self.results[dim]['slice_volume']

testsuite/MDAnalysisTests/analysis/test_lineardensity.py

@@ -22,18 +22,91 @@
 #
 import MDAnalysis as mda
 import numpy as np
+import pytest
 
 from MDAnalysisTests.datafiles import waterPSF, waterDCD
 from MDAnalysis.analysis.lineardensity import LinearDensity
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_allclose
 
 
-def test_serial():
+def test_invalid_grouping():
+    """Invalid groupings raise AttributeError"""
     universe = mda.Universe(waterPSF, waterDCD)
     sel_string = 'all'
     selection = universe.select_atoms(sel_string)
+    with pytest.raises(AttributeError):
+        # centroid is attribute of AtomGroup, but not valid here
+        ld = LinearDensity(selection, grouping="centroid", binsize=5)
+        ld.run()
 
-    xpos = np.array([0., 0., 0., 0.0072334, 0.00473299, 0.,
-                          0., 0., 0., 0.])
-    ld = LinearDensity(selection, binsize=5).run()
-    assert_almost_equal(xpos, ld.results['x']['pos'])
+
+# test data for grouping='atoms'
+expected_masses_atoms = np.array([15.9994, 1.008, 1.008, 15.9994, 1.008, 1.008,
+                                  15.9994, 1.008, 1.008, 15.9994, 1.008, 1.008,
+                                  15.9994, 1.008, 1.008])
+expected_charges_atoms = np.array([-0.834, 0.417, 0.417, -0.834, 0.417,
+                                   0.417, -0.834, 0.417, 0.417, -0.834,
+                                   0.417, 0.417, -0.834, 0.417, 0.417])
+expected_xpos_atoms = np.array([0., 0., 0., 0.0072334, 0.00473299, 0.,
+                                0., 0., 0., 0.])
+expected_xchar_atoms = np.array([0., 0., 0., 2.2158751e-05, -2.2158751e-05,
+                                 0., 0., 0., 0., 0.])
+
+# test data for grouping='residues'
+expected_masses_residues = np.array([18.0154, 18.0154, 18.0154, 18.0154,
+                                     18.0154])
+expected_charges_residues = np.array([0, 0, 0, 0, 0])
+expected_xpos_residues = np.array([0., 0., 0., 0.00717983, 0.00478656,
+                                   0., 0., 0., 0., 0.])
+expected_xchar_residues = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
+
+# test data for grouping='segments'
+expected_masses_segments = np.array([90.0770])
+expected_charges_segments = np.array([0])
+expected_xpos_segments = np.array([0., 0., 0., 0.01196639, 0.,
+                                   0., 0., 0., 0., 0.])
+expected_xchar_segments = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
+
+# test data for grouping='fragments'
+expected_masses_fragments = np.array([18.0154, 18.0154, 18.0154, 18.0154,
+                                      18.0154])
+expected_charges_fragments = np.array([0, 0, 0, 0, 0])
+expected_xpos_fragments = np.array([0., 0., 0., 0.00717983, 0.00478656,
+                                   0., 0., 0., 0., 0.])
+expected_xchar_fragments = np.array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
+
+
+@pytest.mark.parametrize("grouping, expected_masses, expected_charges,\
+                         expected_xpos, expected_xchar", [
+                         ("atoms",
+                          expected_masses_atoms,
+                          expected_charges_atoms,
+                          expected_xpos_atoms,
+                          expected_xchar_atoms),
+                         ("residues",
+                          expected_masses_residues,
+                          expected_charges_residues,
+                          expected_xpos_residues,
+                          expected_xchar_residues),
+                         ("segments",
+                          expected_masses_segments,
+                          expected_charges_segments,
+                          expected_xpos_segments,
+                          expected_xchar_segments),
+                         ("fragments",
+                          expected_masses_fragments,
+                          expected_charges_fragments,
+                          expected_xpos_fragments,
+                          expected_xchar_fragments)
+                         ])
+def test_lineardensity(grouping, expected_masses, expected_charges,
+                       expected_xpos, expected_xchar):
+    universe = mda.Universe(waterPSF, waterDCD)
+    sel_string = 'all'
+    selection = universe.select_atoms(sel_string)
+    ld = LinearDensity(selection, grouping, binsize=5).run()
+    assert_allclose(ld.masses, expected_masses)
+    assert_allclose(ld.charges, expected_charges)
+    # rtol changed here due to floating point imprecision
+    assert_allclose(ld.results['x']['pos'], expected_xpos, rtol=1e-06)
+    assert_allclose(ld.results['x']['char'], expected_xchar)