Changed LinearDensity for Issue MDAnalysis#2508

LinearDensity now works with updating AtomGroups. A new test was
added to verify this.
Result variables "pos" and "char" as well as their "std" were
renamed to "mass_density", "charge_density" and "stddev", respectively.
The bin edges of np.histogram are now saved for xyz as "hist_bin_edges"

package/CHANGELOG

@@ -35,12 +35,18 @@ Fixes
   * Updated the badge in README of MDAnalysisTest to point to Github CI Actions
 
 Enhancements
+  * LinearDensity now works with updating AtomGroups (Issue #2508, PR #    )
   * Improves the RDKitConverter's accuracy (PR #3044): AtomGroups containing
     monatomic ion charges or edge cases with nitrogen, sulfur, phosphorus and
     conjugated systems should now have correctly assigned bond orders and
     charges.
 
 Changes
+  * LinearDensity's results container was changed. (Issue #2508, PR #    )
+    - "pos" is now "mass_density",
+    - "char" is now "charge_density"
+    - "std" entries are now "stddev"
+    - "hist_bin_edges" added to each dimension xyz.
   * ITPParser no longer adds an angle for water molecules that have the SETTLE
     constraint applied to them. This mirrors the way the constraint is handled
     in AMBER. (Issue #2417, PR #3564)

package/MDAnalysis/analysis/lineardensity.py

@@ -57,16 +57,19 @@ class LinearDensity(AnalysisBase):
     ----------
     results.x.dim : int
            index of the [xyz] axes
-    results.x.pos : numpy.ndarray
+    results.x.mass_density : numpy.ndarray
            mass density in [xyz] direction
-    results.x.pos_std : numpy.ndarray
+    results.x.mass_density_stddev : numpy.ndarray
            standard deviation of the mass density in [xyz] direction
-    results.x.char : numpy.ndarray
+    results.x.charge_density : numpy.ndarray
            charge density in [xyz] direction
-    results.x.char_std : numpy.ndarray
+    results.x.charge_density_stddev : numpy.ndarray
            standard deviation of the charge density in [xyz] direction
     results.x.slice_volume : float
            volume of bin in [xyz] direction
+    results.x.hist_bin_edges : numpy.ndarray
+           edges of histogram bins for mass/charge densities, useful for, e.g.,
+           plotting of histogram data.
 
     Example
     -------
@@ -78,7 +81,7 @@ class LinearDensity(AnalysisBase):
 
        ldens = LinearDensity(selection)
        ldens.run()
-       print(ldens.results.x.pos)
+       print(ldens.results.x.mass_density)
 
 
     Alternatively, other types of grouping can be selected using the
@@ -115,6 +118,15 @@ class LinearDensity(AnalysisBase):
        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.
+       LinearDensity now works with updating atom groups.
+       Changed names of result containers
+         :attr:`results.x.pos` -> :attr:`results.x.mass_density`
+         :attr:`results.x.pos_std` -> :attr:`results.x.mass_density_stddev`
+         :attr:`results.x.char` -> :attr:`results.x.charge_density`
+         :attr:`results.x.char_std` -> :attr:`results.x.charge_density_stddev`
+       Added new result container :attr:`results.x.hist_bin_edges`. It contains
+       the bin edges of the histrogram bins for calculated densities and can be
+       used for easier plotting of histogram data.
     """
 
     def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs):
@@ -145,7 +157,8 @@ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs):
         self.nbins = bins.max()
         slices_vol = self.volume / bins
 
-        self.keys = ['pos', 'pos_std', 'char', 'char_std']
+        self.keys = ['mass_density', 'mass_density_stddev',
+                     'charge_density', 'charge_density_stddev']
 
         # Initialize results array with zeros
         for dim in self.results:
@@ -154,12 +167,12 @@ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs):
             for key in self.keys:
                 self.results[dim][key] = np.zeros(self.nbins)
 
-        # Variables later defined in _prepare() method
+        # Variables later defined in _single_frame() method
         self.masses = None
         self.charges = None
         self.totalmass = None
 
-    def _prepare(self):
+    def _single_frame(self):
         # Get masses and charges for the selection
         if self.grouping == "atoms":
             self.masses = self._ags[0].masses
@@ -175,7 +188,6 @@ def _prepare(self):
 
         self.totalmass = np.sum(self.masses)
 
-    def _single_frame(self):
         self.group = getattr(self._ags[0], self.grouping)
         self._ags[0].wrap(compound=self.grouping)
 
@@ -189,8 +201,8 @@ def _single_frame(self):
         for dim in ['x', 'y', 'z']:
             idx = self.results[dim]['dim']
 
-            key = 'pos'
-            key_std = 'pos_std'
+            key = 'mass_density'
+            key_std = 'mass_density_stddev'
             # histogram for positions weighted on masses
             hist, _ = np.histogram(positions[:, idx],
                                    weights=self.masses,
@@ -200,39 +212,42 @@ def _single_frame(self):
             self.results[dim][key] += hist
             self.results[dim][key_std] += np.square(hist)
 
-            key = 'char'
-            key_std = 'char_std'
+            key = 'charge_density'
+            key_std = 'charge_density_stddev'
             # histogram for positions weighted on charges
-            hist, _ = np.histogram(positions[:, idx],
-                                   weights=self.charges,
-                                   bins=self.nbins,
-                                   range=(0.0, max(self.dimensions)))
+            hist, bin_edges = np.histogram(positions[:, idx],
+                                           weights=self.charges,
+                                           bins=self.nbins,
+                                           range=(0.0, max(self.dimensions)))
 
             self.results[dim][key] += hist
             self.results[dim][key_std] += np.square(hist)
+            self.results[dim]['hist_bin_edges'] = bin_edges
 
     def _conclude(self):
         k = 6.022e-1  # divide by avodagro and convert from A3 to cm3
 
         # Average results over the number of configurations
         for dim in ['x', 'y', 'z']:
-            for key in ['pos', 'pos_std', 'char', 'char_std']:
+            for key in ['mass_density', 'mass_density_stddev',
+                        'charge_density', 'charge_density_stddev']:
                 self.results[dim][key] /= self.n_frames
             # Compute standard deviation for the error
             # 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
+            # radicand_mass and radicand_charge are therefore calculated first
             # 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_mass = self.results[dim]['mass_density_stddev'] - \
+                np.square(self.results[dim]['mass_density'])
+            radicand_mass[radicand_mass < 0] = 0
+            self.results[dim]['mass_density_stddev'] = np.sqrt(radicand_mass)
 
-            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)
+            radicand_charge = self.results[dim]['charge_density_stddev'] - \
+                np.square(self.results[dim]['charge_density'])
+            radicand_charge[radicand_charge < 0] = 0
+            self.results[dim]['charge_density_stddev'] = \
+                np.sqrt(radicand_charge)
 
         for dim in ['x', 'y', 'z']:
             norm = k * self.results[dim]['slice_volume']
@@ -243,12 +258,12 @@ def _add_other_results(self, other):
         # For parallel analysis
         results = self.results
         for dim in ['x', 'y', 'z']:
-            key = 'pos'
-            key_std = 'pos_std'
+            key = 'mass_density'
+            key_std = 'mass_density_stddev'
             results[dim][key] += other[dim][key]
             results[dim][key_std] += other[dim][key_std]
 
-            key = 'char'
-            key_std = 'char_std'
+            key = 'charge_density'
+            key_std = 'charge_density_stddev'
             results[dim][key] += other[dim][key]
             results[dim][key_std] += other[dim][key_std]

testsuite/MDAnalysisTests/analysis/test_lineardensity.py

@@ -27,6 +27,7 @@
 from MDAnalysisTests.datafiles import waterPSF, waterDCD
 from MDAnalysis.analysis.lineardensity import LinearDensity
 from numpy.testing import assert_allclose
+from MDAnalysis.core._get_readers import get_reader_for
 
 
 def test_invalid_grouping():
@@ -47,53 +48,93 @@ def test_invalid_grouping():
 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.,
+expected_xmass_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.])
+expected_xcharge_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,
+expected_xmass_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.])
+expected_xcharge_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.,
+expected_xmass_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.])
+expected_xcharge_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.])
+expected_xmass_fragments = np.array([0., 0., 0., 0.00717983, 0.00478656,
+                                     0., 0., 0., 0., 0.])
+expected_xcharge_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", [
+@pytest.mark.parametrize("grouping, expected_masses, expected_charges, expected_xmass, expected_xcharge", [
     ("atoms", expected_masses_atoms, expected_charges_atoms,
-     expected_xpos_atoms, expected_xchar_atoms),
+     expected_xmass_atoms, expected_xcharge_atoms),
     ("residues", expected_masses_residues, expected_charges_residues,
-     expected_xpos_residues, expected_xchar_residues),
+     expected_xmass_residues, expected_xcharge_residues),
     ("segments", expected_masses_segments, expected_charges_segments,
-     expected_xpos_segments, expected_xchar_segments),
+     expected_xmass_segments, expected_xcharge_segments),
     ("fragments", expected_masses_fragments, expected_charges_fragments,
-     expected_xpos_fragments, expected_xchar_fragments)
+     expected_xmass_fragments, expected_xcharge_fragments)
 ])
 def test_lineardensity(grouping, expected_masses, expected_charges,
-                       expected_xpos, expected_xchar):
+                       expected_xmass, expected_xcharge):
     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)
+    assert_allclose(ld.results.x.mass_density, expected_xmass, rtol=1e-06)
+    assert_allclose(ld.results.x.charge_density, expected_xcharge)
+
+
+def make_Universe_updating_atomgroup():
+    """Generate a universe for testing whether LinearDensity works with
+       updating atom groups."""
+    n_atoms = 3
+    u = mda.Universe.empty(n_atoms=n_atoms,
+                           n_residues=n_atoms,
+                           n_segments=n_atoms,
+                           atom_resindex=np.arange(n_atoms),
+                           residue_segindex=np.arange(n_atoms))
+
+    for attr in ["charges", "masses"]:
+        u.add_TopologyAttr(attr, values=np.ones(n_atoms))
+
+    coords = np.array([
+                      [[1., 1., 1.], [1., 2., 1.], [2., 1., 1.]],
+                      [[1., 1., 2.], [1., 2., 1.], [2., 1., 1.]],
+                      [[1., 1., 3.], [1., 2., 1.], [2., 1., 1.]],
+                      [[1., 1., 4.], [1., 2., 1.], [2., 1., 1.]],
+                      [[1., 1., 5.], [1., 2., 1.], [2., 1., 1.]]
+                      ])
+
+    u.trajectory = get_reader_for(coords)(coords, order='fac', n_atoms=n_atoms)
+
+    for ts in u.trajectory:
+        ts.dimensions = np.array([2, 2, 6, 90, 90, 90])
+
+    return u
+
+
+def test_updating_atomgroup():
+    expected_z_pos = np.array([0., 0.91331783, 0.08302889, 0., 0., 0.])
+    u = make_Universe_updating_atomgroup()
+    selection = u.select_atoms("prop z < 3", updating=True)
+    ld = LinearDensity(selection, binsize=1).run()
+    assert_allclose(ld.results.z.mass_density, expected_z_pos)
+    # Test whether histogram bins are saved correctly.
+    expected_bin_edges = np.arange(0, 7)
+    assert_allclose(ld.results.x.hist_bin_edges, expected_bin_edges)