From 389eb1ce7e134a4589db09a98f6692565229a869 Mon Sep 17 00:00:00 2001 From: BFedder Date: Thu, 24 Mar 2022 12:32:33 +0000 Subject: [PATCH 01/12] LinearDensity now working with residues/segments LinearDensity did not work with grouping='residues' and grouping='segments'. This is now fixed, and more tests have been added. Changes made: MDAnalysis/analysis/lineardensity.py: - in _prepare(self): - try/except statement replaced by if statements checking the grouping specified for LinearDensity. - masses and charges are obtained in a way appropriate for each case. - in _single_frame(self): - 'positions' variable assignment changed from elem.centroid() to elem.atoms.centroid() so it works for residues/segments - in _conclude(self): - Floating point arithmetic issues could lead to negative radicands. A check was implemented to set numbers which are close to zero to precisely zero before the square root (for standard deviation) is calculated. MDAnalysisTests/analysis/test_lineardensity.py: - renamed previous test function and added expected masses and charges - added test functions for the other groupings with their respective expected masses and charges --- package/CHANGELOG | 2 + package/MDAnalysis/analysis/lineardensity.py | 46 ++++++++++---- .../analysis/test_lineardensity.py | 60 ++++++++++++++++++- 3 files changed, 95 insertions(+), 13 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 1e9af4f7100..dca0dfdda26 100644 --- a/package/CHANGELOG +++ b/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 #3566 ) * Fixed DumpReader triclinic box dimensions reading. (Issue #3386, PR #3403) * Updated the badge in README of MDAnalysisTest to point to Github CI Actions diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 78ce2fe4433..ef419bf7e6e 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -96,6 +96,10 @@ 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. """ def __init__(self, select, grouping='atoms', binsize=0.25, **kwargs): @@ -146,13 +150,19 @@ def _prepare(self): 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 + if self.grouping in ["residues", "segments"]: + self.masses = np.array([elem.atoms.total_mass() for elem in group]) + self.charges = np.array( + [elem.atoms.total_charge() for elem in group]) + + if self.grouping == "fragments": + self.masses = np.array([elem.total_mass() for elem in group]) + self.charges = np.array([elem.total_charge() for elem in group]) + 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 geometry for residues, segments, fragments + positions = np.array([elem.atoms.centroid() for elem in self.group]) for dim in ['x', 'y', 'z']: idx = self.results[dim]['dim'] @@ -183,6 +193,7 @@ def _single_frame(self): key = 'char' key_std = 'char_std' # histogram for positions weighted on charges + hist, _ = np.histogram(positions[:, idx], weights=self.charges, bins=self.nbins, @@ -191,6 +202,7 @@ def _single_frame(self): self.results[dim][key] += hist self.results[dim][key_std] += np.square(hist) + def _conclude(self): k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 @@ -199,10 +211,24 @@ 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), + # rather than actual value of zero. radicand_pos and radicand_char + # are therefore tested for closeness to 0 first and changed where + # appropriate before the square root is calculated. + radicand_pos = self.results[dim][ + 'pos_std'] - np.square(self.results[dim]['pos']) + for i in range(len(radicand_pos)): + if np.isclose(radicand_pos, np.zeros(self.nbins))[i]: + radicand_pos[i] = 0 + self.results[dim]['pos_std'] = np.sqrt(radicand_pos) + + radicand_char = self.results[dim][ + 'char_std'] - np.square(self.results[dim]['char']) + for i in range(len(radicand_char)): + if np.isclose(radicand_char, np.zeros(self.nbins))[i]: + radicand_char[i] = 0 + self.results[dim]['char_std'] = np.sqrt(radicand_char) for dim in ['x', 'y', 'z']: norm = k * self.results[dim]['slice_volume'] diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 2fd2533046e..55e1149f691 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -28,12 +28,66 @@ from numpy.testing import assert_almost_equal -def test_serial(): + +def test_grouping_atoms(): + """For testing the case of grouping='atoms'""" universe = mda.Universe(waterPSF, waterDCD) sel_string = 'all' selection = universe.select_atoms(sel_string) + ld = LinearDensity(selection, grouping="atoms", binsize=5).run() + + expected_masses = 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 = 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]) + 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']) + + assert_almost_equal(ld.masses, expected_masses) + assert_almost_equal(ld.charges, expected_charges) + assert_almost_equal(ld.results['x']['pos'], xpos) + +def test_grouping_residues(): + """For testing the case of grouping='residues'""" + universe = mda.Universe(waterPSF, waterDCD) + sel_string = 'all' + selection = universe.select_atoms(sel_string) + ld = LinearDensity(selection, grouping="residues", binsize=5).run() + + expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) + expected_charges = np.array([0, 0, 0, 0, 0]) + + assert_almost_equal(ld.masses, expected_masses) + assert_almost_equal(ld.charges, expected_charges) + + +def test_grouping_segments(): + """For testing the case of grouping='segments'""" + universe = mda.Universe(waterPSF, waterDCD) + sel_string = 'all' + selection = universe.select_atoms(sel_string) + ld = LinearDensity(selection, grouping="segments", binsize=5).run() + + expected_masses = np.array([90.0770]) + expected_charges = np.array([0]) + assert_almost_equal(ld.masses, expected_masses) + assert_almost_equal(ld.charges, expected_charges) + +def test_grouping_fragments(): + """For testing the case of grouping='fragments'""" + universe = mda.Universe(waterPSF, waterDCD) + sel_string = 'all' + selection = universe.select_atoms(sel_string) + ld = LinearDensity(selection, grouping="fragments", binsize=5).run() + + expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) + expected_charges = np.array([0, 0, 0, 0, 0]) + + assert_almost_equal(ld.masses, expected_masses) + assert_almost_equal(ld.charges, expected_charges) + From d96245ec77a70316b3c0b61ff255ee28f8a8a79e Mon Sep 17 00:00:00 2001 From: BFedder Date: Thu, 24 Mar 2022 16:54:28 +0000 Subject: [PATCH 02/12] changes requested by IAlibay Changed the if/if/if for reading masses and charges to if/elif/else. Will now raise an AttributeError if invalid grouping is chosen. 'residues', 'segments' and 'fragments' now use the same logic (elem.atoms.total_mass()). --- package/MDAnalysis/analysis/lineardensity.py | 25 +++++--------------- 1 file changed, 6 insertions(+), 19 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index ef419bf7e6e..77c73ed30f6 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -23,7 +23,6 @@ """ Linear Density --- :mod:`MDAnalysis.analysis.lineardensity` =========================================================== - A tool to compute mass and charge density profiles along the three cartesian axes [xyz] of the simulation cell. Works only for orthorombic, fixed volume cells (thus for simulations in canonical NVT ensemble). @@ -37,7 +36,6 @@ class LinearDensity(AnalysisBase): """Linear density profile - Parameters ---------- select : AtomGroup @@ -51,7 +49,6 @@ class LinearDensity(AnalysisBase): profile (smaller --> higher resolution) verbose : bool, optional Show detailed progress of the calculation if set to ``True`` - Attributes ---------- results.x.dim : int @@ -66,37 +63,28 @@ class LinearDensity(AnalysisBase): standard deviation of the charge density in [xyz] direction results.x.slice_volume : float volume of bin in [xyz] direction - Example ------- First create a ``LinearDensity`` object by supplying a selection, then use the :meth:`run` method. Finally access the results stored in results, i.e. the mass density in the x direction. - .. code-block:: python - ldens = LinearDensity(selection) ldens.run() print(ldens.results.x.pos) - - .. versionadded:: 0.14.0 - .. versionchanged:: 1.0.0 Support for the ``start``, ``stop``, and ``step`` keywords has been removed. These should instead be passed to :meth:`LinearDensity.run`. The ``save()`` method was also removed, you can use ``np.savetxt()`` or ``np.save()`` on the :attr:`LinearDensity.results` dictionary contents instead. - .. versionchanged:: 1.0.0 Changed `selection` keyword to `select` - .. versionchanged:: 2.0.0 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. @@ -154,14 +142,14 @@ def _prepare(self): self.masses = self._ags[0].masses self.charges = self._ags[0].charges - if self.grouping in ["residues", "segments"]: + elif self.grouping in ["residues", "segments", "fragments"]: self.masses = np.array([elem.atoms.total_mass() for elem in group]) self.charges = np.array( [elem.atoms.total_charge() for elem in group]) - if self.grouping == "fragments": - self.masses = np.array([elem.total_mass() for elem in group]) - self.charges = np.array([elem.total_charge() for elem in group]) + else: + raise AttributeError( + f"{self.grouping} is not a valid value for grouping.") self.totalmass = np.sum(self.masses) @@ -174,7 +162,8 @@ def _single_frame(self): positions = self._ags[0].positions # faster for atoms else: # Centre of geometry for residues, segments, fragments - positions = np.array([elem.atoms.centroid() for elem in self.group]) + positions = np.array( + [elem.atoms.centroid() for elem in self.group]) for dim in ['x', 'y', 'z']: idx = self.results[dim]['dim'] @@ -193,7 +182,6 @@ def _single_frame(self): key = 'char' key_std = 'char_std' # histogram for positions weighted on charges - hist, _ = np.histogram(positions[:, idx], weights=self.charges, bins=self.nbins, @@ -202,7 +190,6 @@ def _single_frame(self): self.results[dim][key] += hist self.results[dim][key_std] += np.square(hist) - def _conclude(self): k = 6.022e-1 # divide by avodagro and convert from A3 to cm3 From 54768408187e873178c000fdaa0dbeb6d929d7e1 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Thu, 24 Mar 2022 17:43:20 +0000 Subject: [PATCH 03/12] Update test_lineardensity.py --- testsuite/MDAnalysisTests/analysis/test_lineardensity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 55e1149f691..f7a413d5328 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -28,7 +28,6 @@ from numpy.testing import assert_almost_equal - def test_grouping_atoms(): """For testing the case of grouping='atoms'""" universe = mda.Universe(waterPSF, waterDCD) @@ -52,6 +51,7 @@ def test_grouping_atoms(): assert_almost_equal(ld.charges, expected_charges) assert_almost_equal(ld.results['x']['pos'], xpos) + def test_grouping_residues(): """For testing the case of grouping='residues'""" universe = mda.Universe(waterPSF, waterDCD) @@ -78,6 +78,7 @@ def test_grouping_segments(): assert_almost_equal(ld.masses, expected_masses) assert_almost_equal(ld.charges, expected_charges) + def test_grouping_fragments(): """For testing the case of grouping='fragments'""" universe = mda.Universe(waterPSF, waterDCD) @@ -90,4 +91,3 @@ def test_grouping_fragments(): assert_almost_equal(ld.masses, expected_masses) assert_almost_equal(ld.charges, expected_charges) - From 8b23b5d506c7a72318eb40cd9efec3931440a4ed Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Thu, 24 Mar 2022 17:51:27 +0000 Subject: [PATCH 04/12] Update lineardensity.py --- package/MDAnalysis/analysis/lineardensity.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 77c73ed30f6..deea2892ccc 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -23,6 +23,7 @@ """ Linear Density --- :mod:`MDAnalysis.analysis.lineardensity` =========================================================== + A tool to compute mass and charge density profiles along the three cartesian axes [xyz] of the simulation cell. Works only for orthorombic, fixed volume cells (thus for simulations in canonical NVT ensemble). @@ -36,6 +37,7 @@ class LinearDensity(AnalysisBase): """Linear density profile + Parameters ---------- select : AtomGroup @@ -49,6 +51,7 @@ class LinearDensity(AnalysisBase): profile (smaller --> higher resolution) verbose : bool, optional Show detailed progress of the calculation if set to ``True`` + Attributes ---------- results.x.dim : int @@ -63,28 +66,37 @@ class LinearDensity(AnalysisBase): standard deviation of the charge density in [xyz] direction results.x.slice_volume : float volume of bin in [xyz] direction + Example ------- First create a ``LinearDensity`` object by supplying a selection, then use the :meth:`run` method. Finally access the results stored in results, i.e. the mass density in the x direction. + .. code-block:: python + ldens = LinearDensity(selection) ldens.run() print(ldens.results.x.pos) + + .. versionadded:: 0.14.0 + .. versionchanged:: 1.0.0 Support for the ``start``, ``stop``, and ``step`` keywords has been removed. These should instead be passed to :meth:`LinearDensity.run`. The ``save()`` method was also removed, you can use ``np.savetxt()`` or ``np.save()`` on the :attr:`LinearDensity.results` dictionary contents instead. + .. versionchanged:: 1.0.0 Changed `selection` keyword to `select` + .. versionchanged:: 2.0.0 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. From ba629537297420a203655e7e761867d42d6c9b89 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Thu, 24 Mar 2022 18:10:28 +0000 Subject: [PATCH 05/12] Update test_lineardensity.py Exchanged assert_almost_equal with assert_allclose --- .../analysis/test_lineardensity.py | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index f7a413d5328..716e8ee3672 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -25,7 +25,7 @@ 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_grouping_atoms(): @@ -47,9 +47,10 @@ def test_grouping_atoms(): xpos = np.array([0., 0., 0., 0.0072334, 0.00473299, 0., 0., 0., 0., 0.]) - assert_almost_equal(ld.masses, expected_masses) - assert_almost_equal(ld.charges, expected_charges) - assert_almost_equal(ld.results['x']['pos'], xpos) + 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'], xpos, rtol=1e-06) def test_grouping_residues(): @@ -62,8 +63,8 @@ def test_grouping_residues(): expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) expected_charges = np.array([0, 0, 0, 0, 0]) - assert_almost_equal(ld.masses, expected_masses) - assert_almost_equal(ld.charges, expected_charges) + assert_allclose(ld.masses, expected_masses) + assert_allclose(ld.charges, expected_charges) def test_grouping_segments(): @@ -75,8 +76,8 @@ def test_grouping_segments(): expected_masses = np.array([90.0770]) expected_charges = np.array([0]) - assert_almost_equal(ld.masses, expected_masses) - assert_almost_equal(ld.charges, expected_charges) + assert_allclose(ld.masses, expected_masses) + assert_allclose(ld.charges, expected_charges) def test_grouping_fragments(): @@ -89,5 +90,5 @@ def test_grouping_fragments(): expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) expected_charges = np.array([0, 0, 0, 0, 0]) - assert_almost_equal(ld.masses, expected_masses) - assert_almost_equal(ld.charges, expected_charges) + assert_allclose(ld.masses, expected_masses) + assert_allclose(ld.charges, expected_charges) From 5e51f5a152243632ad441d9cf65b9f5d383b87b6 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Fri, 25 Mar 2022 10:01:08 +0000 Subject: [PATCH 06/12] Update CHANGELOG --- package/CHANGELOG | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index dca0dfdda26..a8a25f6c6e4 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,7 +19,7 @@ The rules for this file: Fixes * Fixed LinearDensity not working with grouping='residues' or 'segments' - (Issue #3571, PR #3566 ) + (Issue #3571, PR #3572) * Fixed DumpReader triclinic box dimensions reading. (Issue #3386, PR #3403) * Updated the badge in README of MDAnalysisTest to point to Github CI Actions From 0aa743e20bae5ca04c66218716486243e999485a Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Thu, 31 Mar 2022 11:01:40 +0100 Subject: [PATCH 07/12] Update lineardensity.py Updated check of radicands in _conclude(self). --- package/MDAnalysis/analysis/lineardensity.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index deea2892ccc..cee490933bc 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -212,21 +212,18 @@ def _conclude(self): # Compute standard deviation for the error # For certain tests in testsuite, floating point imprecision # can lead to negative radicands of tiny magnitude (yielding nan), - # rather than actual value of zero. radicand_pos and radicand_char - # are therefore tested for closeness to 0 first and changed where - # appropriate before the square root is calculated. + # or tiny positive values, both of which should actually be zero. + # radicand_pos and radicand_char are therefore calculated first and + # and values that should be 0 are changed appropriately before + # the square root is calculated. radicand_pos = self.results[dim][ 'pos_std'] - np.square(self.results[dim]['pos']) - for i in range(len(radicand_pos)): - if np.isclose(radicand_pos, np.zeros(self.nbins))[i]: - radicand_pos[i] = 0 + radicand_pos[radicand_pos < 1e-9] = 0 self.results[dim]['pos_std'] = np.sqrt(radicand_pos) radicand_char = self.results[dim][ 'char_std'] - np.square(self.results[dim]['char']) - for i in range(len(radicand_char)): - if np.isclose(radicand_char, np.zeros(self.nbins))[i]: - radicand_char[i] = 0 + radicand_char[radicand_char < 1e-9] = 0 self.results[dim]['char_std'] = np.sqrt(radicand_char) for dim in ['x', 'y', 'z']: From 89725dd6c38bc3fb2d643b481df03f0d5a7db654 Mon Sep 17 00:00:00 2001 From: Bjarne Date: Mon, 4 Apr 2022 15:23:47 +0100 Subject: [PATCH 08/12] Addressing reviewer comments * Changed test of radicand to only convert negative numbers to 0 * Now uses CoM rather than CoG for density calculations * Changed Class docstring to reflect change to CoM * Now uses parameter for compound instead of list comprehension when calculating masses, charges, and CoM * refactored test_lineardensity.py, make use of pytest functionality better * Added test for case of invalid grouping selection --- package/MDAnalysis/analysis/lineardensity.py | 34 ++--- .../analysis/test_lineardensity.py | 130 ++++++++++-------- 2 files changed, 89 insertions(+), 75 deletions(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index cee490933bc..303ab87c8a6 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -28,7 +28,7 @@ cartesian axes [xyz] of the simulation cell. Works only for orthorombic, fixed volume cells (thus for simulations in canonical NVT ensemble). """ -import os.path as path +import os.path as pat import numpy as np @@ -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 @@ -100,6 +101,8 @@ class LinearDensity(AnalysisBase): .. 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): @@ -145,19 +148,14 @@ 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 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 = np.array([elem.atoms.total_mass() for elem in group]) - self.charges = np.array( - [elem.atoms.total_charge() for elem in group]) + self.masses = self._ags[0].total_mass(compound=self.grouping) + self.charges = self._ags[0].total_charge(compound=self.grouping) else: raise AttributeError( @@ -173,9 +171,8 @@ def _single_frame(self): if self.grouping == 'atoms': positions = self._ags[0].positions # faster for atoms else: - # Centre of geometry for residues, segments, fragments - positions = np.array( - [elem.atoms.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'] @@ -211,19 +208,18 @@ def _conclude(self): 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), - # or tiny positive values, both of which should actually be zero. + # can lead to negative radicands of tiny magnitude (yielding nan). # radicand_pos and radicand_char are therefore calculated first and - # and values that should be 0 are changed appropriately before - # the square root is calculated. + # 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 < 1e-9] = 0 + 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 < 1e-9] = 0 + radicand_char[radicand_char < 0] = 0 self.results[dim]['char_std'] = np.sqrt(radicand_char) for dim in ['x', 'y', 'z']: diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 716e8ee3672..1cb7b6f0e3e 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -22,73 +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_allclose -def test_grouping_atoms(): - """For testing the case of grouping='atoms'""" +def test_invalid_grouping(): + """Invalid groupings raise AttributeError""" universe = mda.Universe(waterPSF, waterDCD) sel_string = 'all' selection = universe.select_atoms(sel_string) - - ld = LinearDensity(selection, grouping="atoms", binsize=5).run() - - expected_masses = 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 = 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]) - - xpos = np.array([0., 0., 0., 0.0072334, 0.00473299, 0., - 0., 0., 0., 0.]) - - 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'], xpos, rtol=1e-06) - - -def test_grouping_residues(): - """For testing the case of grouping='residues'""" - universe = mda.Universe(waterPSF, waterDCD) - sel_string = 'all' - selection = universe.select_atoms(sel_string) - ld = LinearDensity(selection, grouping="residues", binsize=5).run() - - expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) - expected_charges = np.array([0, 0, 0, 0, 0]) - - assert_allclose(ld.masses, expected_masses) - assert_allclose(ld.charges, expected_charges) - - -def test_grouping_segments(): - """For testing the case of grouping='segments'""" + with pytest.raises(AttributeError): + # centroid is attribute of AtomGroup, but not valid here + ld = LinearDensity(selection, grouping="centroid", binsize=5) + ld.run() + + +# 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="segments", binsize=5).run() - - expected_masses = np.array([90.0770]) - expected_charges = np.array([0]) - assert_allclose(ld.masses, expected_masses) - assert_allclose(ld.charges, expected_charges) - - -def test_grouping_fragments(): - """For testing the case of grouping='fragments'""" - universe = mda.Universe(waterPSF, waterDCD) - sel_string = 'all' - selection = universe.select_atoms(sel_string) - ld = LinearDensity(selection, grouping="fragments", binsize=5).run() - - expected_masses = np.array([18.0154, 18.0154, 18.0154, 18.0154, 18.0154]) - expected_charges = np.array([0, 0, 0, 0, 0]) - + 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) From ea4e5611c3b36f840704987988b3035687d2b646 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Mon, 4 Apr 2022 15:33:40 +0100 Subject: [PATCH 09/12] Update lineardensity.py --- package/MDAnalysis/analysis/lineardensity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 303ab87c8a6..40238f0c5d8 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -28,7 +28,7 @@ cartesian axes [xyz] of the simulation cell. Works only for orthorombic, fixed volume cells (thus for simulations in canonical NVT ensemble). """ -import os.path as pat +import os.path as path import numpy as np From 9d2eb68ba3d039bd4cc2a6814f6aefbf53e7218f Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Tue, 5 Apr 2022 10:17:39 +0100 Subject: [PATCH 10/12] Update lineardensity.py Added second example to class docstring --- package/MDAnalysis/analysis/lineardensity.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index 40238f0c5d8..c5957bb039e 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -80,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 From 446de45d44cfd9dd5da85c83a6a128044370ad68 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Fri, 8 Apr 2022 18:49:42 +0100 Subject: [PATCH 11/12] Update package/MDAnalysis/analysis/lineardensity.py Co-authored-by: Irfan Alibay --- package/MDAnalysis/analysis/lineardensity.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/lineardensity.py b/package/MDAnalysis/analysis/lineardensity.py index c5957bb039e..a3bdaf28c77 100644 --- a/package/MDAnalysis/analysis/lineardensity.py +++ b/package/MDAnalysis/analysis/lineardensity.py @@ -80,8 +80,15 @@ class LinearDensity(AnalysisBase): ldens.run() print(ldens.results.x.pos) - # alternatively, specify grouping/binsize like this + + Alternatively, the other types of groupings can be selected using the ``grouping`` + keyword. For example to calculated the density based on the :class:`ResidueGroup`s + of the system: + + .. code-block:: python ldens = LinearDensity(selection, grouping='residues', binsize=1.0) + ldens.run() + .. versionadded:: 0.14.0 From 565bd4dd195029c2f3879e37d42ce9c709fab3d4 Mon Sep 17 00:00:00 2001 From: BFedder <80363742+BFedder@users.noreply.github.com> Date: Fri, 8 Apr 2022 18:50:03 +0100 Subject: [PATCH 12/12] Update testsuite/MDAnalysisTests/analysis/test_lineardensity.py Co-authored-by: Irfan Alibay --- .../analysis/test_lineardensity.py | 33 ++++++------------- 1 file changed, 10 insertions(+), 23 deletions(-) diff --git a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py index 1cb7b6f0e3e..00c62cc83aa 100644 --- a/testsuite/MDAnalysisTests/analysis/test_lineardensity.py +++ b/testsuite/MDAnalysisTests/analysis/test_lineardensity.py @@ -76,29 +76,16 @@ def test_invalid_grouping(): 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) - ]) +@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)