From 4819464618de18bc0460eb11bb6dfba132193f48 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Mon, 17 Jun 2019 13:15:32 +0530 Subject: [PATCH 01/10] Adds unwrap to center --- package/MDAnalysis/core/groups.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index e9dcfecc458..4ee7af04608 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -649,7 +649,7 @@ def isunique(self): return not np.count_nonzero(mask) @warn_if_not_unique - def center(self, weights, pbc=None, compound='group'): + def center(self, weights, pbc=None, compound='group', unwrap=False): """Weighted center of (compounds of) the group Computes the weighted center of :class:`Atoms` in the group. @@ -679,6 +679,11 @@ def center(self, weights, pbc=None, compound='group'): will be returned as an array of position vectors, i.e. a 2d array. Note that, in any case, *only* the positions of :class:`Atoms` *belonging to the group* will be taken into account. + unwrap : bool or None, optional + If ``True`` and `compound` is ``'group'``, the atoms will be unwrapped + before calculations. If ``True`` and `compound` is + ``'segments'`` or ``'residues'``, all molecules will be unwrapped before + calculation to keep the compounds intact. Returns ------- @@ -734,6 +739,8 @@ def center(self, weights, pbc=None, compound='group'): comp = compound.lower() if comp == 'group': + if unwrap: + coords = atoms.unwrap(inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) else: @@ -769,11 +776,19 @@ def center(self, weights, pbc=None, compound='group'): " one of 'group', 'residues', 'segments', " "'molecules', or 'fragments'.".format(compound)) + + # Sort positions and weights by compound index and promote to dtype if # required: sort_indices = np.argsort(compound_indices) compound_indices = compound_indices[sort_indices] - coords = atoms.positions[sort_indices] + + # Unwrap Atoms + if unwrap: + coords = atoms.unwrap(inplace=False) + else: + coords = atoms.positions[sort_indices] + if weights is None: coords = coords.astype(dtype, copy=False) else: From 300c955a98b331186f6ccb815e7b33f410a78290 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 18 Jun 2019 23:26:28 +0530 Subject: [PATCH 02/10] Adds tests --- package/MDAnalysis/core/groups.py | 5 +- .../MDAnalysisTests/core/test_atomgroup.py | 27 +++++++- testsuite/MDAnalysisTests/core/util.py | 65 +++++++++++++++++++ 3 files changed, 93 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 4ee7af04608..45fc061a2c6 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -740,7 +740,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': if unwrap: - coords = atoms.unwrap(inplace=False) + coords = atoms.unwrap(compound=comp, reference=None, inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) else: @@ -785,10 +785,9 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): # Unwrap Atoms if unwrap: - coords = atoms.unwrap(inplace=False) + coords = atoms.unwrap(compound=comp, reference=None, inplace=False) else: coords = atoms.positions[sort_indices] - if weights is None: coords = coords.astype(dtype, copy=False) else: diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 617cdb7b959..6da44f5b6d5 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -52,7 +52,7 @@ GRO ) from MDAnalysisTests import make_Universe, no_deprecated_call - +from MDAnalysisTests.core.util import UnWrapUniverse import pytest @@ -518,6 +518,31 @@ def test_center_wrong_shape(self, ag): with pytest.raises(ValueError): ag.center(weights) + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_center_unwrap(self, level, compound, is_triclinic): + u = UnWrapUniverse(is_triclinic=is_triclinic) + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + + # get the expected results + center = group.center(weights=None, compound=compound, unwrap=True) + + ref_center = u.center(compound=compound) + assert_almost_equal(ref_center, center, decimal=4) + class TestSplit(object): diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 664dbabbef9..ad2b911f38b 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -333,6 +333,7 @@ def __new__(cls, have_bonds=True, have_masses=True, have_molnums=True, # bind custom methods to universe: u.unwrapped_coords = cls.unwrapped_coords.__get__(u) u.wrapped_coords = cls.wrapped_coords.__get__(u) + u.center = cls.center.__get__(u) return u @@ -545,3 +546,67 @@ def wrapped_coords(self, compound, center): positions = relpos * np.array([a, a, a]) return positions.astype(np.float32) + + def center(self, compound): + """Returns centers which correspond to the unwrapped system. + + Parameters + ---------- + compound : {'atoms', 'group', 'segments', 'residues', 'molecules', \ + 'fragments'} + Which type of component is unwrapped. Note that for ``'group'``, + the result will only be correct *if the group is the entire system*. + + Note + ---- + This function assumes that all atom masses are equal. Therefore, the + returned coordinates for ``center='com'`` and ``center='cog'`` are + identical. + """ + + relpos = self.unwrapped_coords(compound, reference=None) + + comp = compound.lower() + if comp not in ['group', 'segments', 'residues', 'molecules', + 'fragments']: + raise ValueError("Unknown unwrap compound: {}".format(compound)) + + pos = 0 + + if compound=="residues": + center_pos = np.zeros((15, 3), dtype=np.float32) + else: + center_pos = np.zeros((12, 3), dtype=np.float32) + + for base in range(3): + loc_centre = relpos[base, :] + center_pos[pos,:] = loc_centre + pos+=1 + + for base in range(3, 15, 3): + loc_centre = np.mean(relpos[base:base + 3, :], axis=0) + center_pos[pos,:] = loc_centre + pos+=1 + + if compound=="residues": + for base in range(15, 47, 4): + loc_centre = np.mean(relpos[base:base + 4, :], axis=0) + center_pos[pos,:] = loc_centre + pos+=1 + else: + for base in range(15, 23, 4): + loc_centre = np.mean(relpos[base:base + 4, :], axis=0) + center_pos[pos,:] = loc_centre + pos+=1 + for base in range(23, 47, 8): + loc_centre = np.mean(relpos[base:base + 8, :], axis=0) + center_pos[pos,:] = loc_centre + pos+=1 + + if compound == "group": + center_pos = center_pos[11] + elif compound == "segments": + center_pos = center_pos[9:] + + return center_pos + From 6f21b3143c60199e95fc50bb7d95de878aeabf33 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 10:41:37 +0530 Subject: [PATCH 03/10] Review changes --- package/MDAnalysis/core/groups.py | 12 +++++------ .../MDAnalysisTests/core/test_atomgroup.py | 7 +++++++ testsuite/MDAnalysisTests/core/util.py | 20 +++++++++---------- 3 files changed, 22 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 45fc061a2c6..1e8ec7df630 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -679,11 +679,8 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): will be returned as an array of position vectors, i.e. a 2d array. Note that, in any case, *only* the positions of :class:`Atoms` *belonging to the group* will be taken into account. - unwrap : bool or None, optional - If ``True`` and `compound` is ``'group'``, the atoms will be unwrapped - before calculations. If ``True`` and `compound` is - ``'segments'`` or ``'residues'``, all molecules will be unwrapped before - calculation to keep the compounds intact. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. Returns ------- @@ -739,6 +736,9 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': + if unwrap and pbc: + raise ValueError("'unwrap' and 'pbc' cannot be true at the same time " + "for compound='group") if unwrap: coords = atoms.unwrap(compound=comp, reference=None, inplace=False) if pbc: @@ -776,8 +776,6 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): " one of 'group', 'residues', 'segments', " "'molecules', or 'fragments'.".format(compound)) - - # Sort positions and weights by compound index and promote to dtype if # required: sort_indices = np.argsort(compound_indices) diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 6da44f5b6d5..c0e9a1b58cc 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -543,6 +543,13 @@ def test_center_unwrap(self, level, compound, is_triclinic): ref_center = u.center(compound=compound) assert_almost_equal(ref_center, center, decimal=4) + def test_center_unwrap_pbc_true(self): + u = UnWrapUniverse(is_triclinic=False) + # select group appropriate for compound: + group = u.atoms[39:47] # molecule 12 + with pytest.raises(ValueError): + group.center(weights=None, compound="group", unwrap=True, pbc=True) + class TestSplit(object): diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index ad2b911f38b..653c743d19c 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -579,28 +579,28 @@ def center(self, compound): center_pos = np.zeros((12, 3), dtype=np.float32) for base in range(3): - loc_centre = relpos[base, :] - center_pos[pos,:] = loc_centre + loc_center = relpos[base, :] + center_pos[pos,:] = loc_center pos+=1 for base in range(3, 15, 3): - loc_centre = np.mean(relpos[base:base + 3, :], axis=0) - center_pos[pos,:] = loc_centre + loc_center = np.mean(relpos[base:base + 3, :], axis=0) + center_pos[pos,:] = loc_center pos+=1 if compound=="residues": for base in range(15, 47, 4): - loc_centre = np.mean(relpos[base:base + 4, :], axis=0) - center_pos[pos,:] = loc_centre + loc_center = np.mean(relpos[base:base + 4, :], axis=0) + center_pos[pos,:] = loc_center pos+=1 else: for base in range(15, 23, 4): - loc_centre = np.mean(relpos[base:base + 4, :], axis=0) - center_pos[pos,:] = loc_centre + loc_center = np.mean(relpos[base:base + 4, :], axis=0) + center_pos[pos,:] = loc_center pos+=1 for base in range(23, 47, 8): - loc_centre = np.mean(relpos[base:base + 8, :], axis=0) - center_pos[pos,:] = loc_centre + loc_center = np.mean(relpos[base:base + 8, :], axis=0) + center_pos[pos,:] = loc_center pos+=1 if compound == "group": From 55d0a84392e2e31e129e640e6a01d63aba00f11c Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 18:26:59 +0530 Subject: [PATCH 04/10] add unwrap to 'cog' --- package/MDAnalysis/core/groups.py | 11 +++++--- .../MDAnalysisTests/core/test_atomgroup.py | 26 +++++++++++++++++++ 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 1e8ec7df630..1b9c06a5955 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -739,12 +739,13 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): if unwrap and pbc: raise ValueError("'unwrap' and 'pbc' cannot be true at the same time " "for compound='group") - if unwrap: - coords = atoms.unwrap(compound=comp, reference=None, inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) else: coords = atoms.positions + if unwrap: + coords = atoms.unwrap(compound=comp, reference=None, inplace=False) + # If there's no atom, return its (empty) coordinates unchanged. if len(atoms) == 0: return coords @@ -816,7 +817,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): return centers @warn_if_not_unique - def center_of_geometry(self, pbc=None, compound='group'): + def center_of_geometry(self, pbc=None, compound='group', unwrap=False): """Center of geometry of (compounds of) the group. Computes the center of geometry (a.k.a. centroid) of @@ -840,6 +841,8 @@ def center_of_geometry(self, pbc=None, compound='group'): will be returned as an array of position vectors, i.e. a 2d array. Note that, in any case, *only* the positions of :class:`Atoms` *belonging to the group* will be taken into account. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. Returns ------- @@ -862,7 +865,7 @@ def center_of_geometry(self, pbc=None, compound='group'): .. versionchanged:: 0.20.0 Added ``'molecules'`` and ``'fragments'`` compounds """ - return self.center(None, pbc=pbc, compound=compound) + return self.center(None, pbc=pbc, compound=compound, unwrap=unwrap) centroid = center_of_geometry diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c0e9a1b58cc..cb4e211f35b 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -834,6 +834,32 @@ def test_chi1_nodep(self, PSFDCD): sel = PSFDCD.segments[0].residues[12].chi1_selection() # LYS +class TestUnwrapFlag(object): + + prec = 3 + + @pytest.fixture() + def ref_noUnwrap(self): + return { + 'COG': np.array([5.1, 7.5, 7. ], dtype=np.float32), + } + + @pytest.fixture() + def ref_Unwrap(self): + return { + 'COG': np.array([10.1, 7.5, 7. ], dtype=np.float32), + } + + def test_default(self, ref_noUnwrap): + u = UnWrapUniverse(is_triclinic=False) + group = u.atoms[31:39] # molecules 11 + assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec) + + def test_UnWrapFlag(self, ref_Unwrap): + u = UnWrapUniverse(is_triclinic=False) + group = u.atoms[31:39] # molecules 11 + assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec) + class TestPBCFlag(object): prec = 3 From 6ae780bea1b64b3542104c0521f1c19809047fc2 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 3 Jul 2019 22:13:26 +0530 Subject: [PATCH 05/10] check for pbc and unwrap --- package/CHANGELOG | 1 + package/MDAnalysis/core/groups.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 997e420616f..65684073954 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -57,6 +57,7 @@ Enhancements * survival probability additions: residues, intermittency, step with performance, (PR #2226) * added unwrap keyword to center (PR #2275) + * added unwrap keyword to center_of_geometry (PR #2279) Changes * added official support for Python 3.7 (PR #1963) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 129f2cb5de6..0b338132ab7 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -833,6 +833,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): return centers @warn_if_not_unique + @check_pbc_and_unwrap def center_of_geometry(self, pbc=None, compound='group', unwrap=False): """Center of geometry of (compounds of) the group. @@ -880,6 +881,7 @@ def center_of_geometry(self, pbc=None, compound='group', unwrap=False): .. versionchanged:: 0.19.0 Added `compound` parameter .. versionchanged:: 0.20.0 Added ``'molecules'`` and ``'fragments'`` compounds + .. versionchanged:: 0.20.0 Added `unwrap` parameter """ return self.center(None, pbc=pbc, compound=compound, unwrap=unwrap) From 6d33fe940b4069f1f930064893fd53fd221cd792 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 4 Jul 2019 03:04:34 +0530 Subject: [PATCH 06/10] Removes unnecessary lines --- package/MDAnalysis/core/groups.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 0b338132ab7..ef80087d794 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -750,18 +750,12 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': - if unwrap and pbc: - raise ValueError("'unwrap' and 'pbc' cannot be true at the same time " - "for compound='group") if pbc: coords = atoms.pack_into_box(inplace=False) elif unwrap: coords = atoms.unwrap(compound=comp, reference=None, inplace=False) else: coords = atoms.positions - if unwrap: - coords = atoms.unwrap(compound=comp, reference=None, inplace=False) - # If there's no atom, return its (empty) coordinates unchanged. if len(atoms) == 0: return coords From 017a6d252750e61e633c8c0fa75d84462f00be27 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 12:49:06 +0530 Subject: [PATCH 07/10] Adds unwrap to asphericity --- package/MDAnalysis/core/topologyattrs.py | 20 ++++++++--- .../MDAnalysisTests/core/test_atomgroup.py | 36 +++++++++++++++++++ 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index fe1d8b47919..c37e1b06abe 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1020,7 +1020,8 @@ def shape_parameter(group, **kwargs): ('shape_parameter', shape_parameter)) @warn_if_not_unique - def asphericity(group, pbc=None): + @check_pbc_and_unwrap + def asphericity(group, pbc=None, unwrap=None, compound='group'): """Asphericity. See [Dima2004b]_ for background information. @@ -1031,6 +1032,10 @@ def asphericity(group, pbc=None): If ``True``, move all atoms within the primary unit cell before calculation. If ``None`` use value defined in MDAnalysis.core.flags['use_pbc'] + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + Which type of component to keep together during unwrapping. Note ---- @@ -1051,6 +1056,7 @@ def asphericity(group, pbc=None): .. versionadded:: 0.7.7 .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms @@ -1058,12 +1064,16 @@ def asphericity(group, pbc=None): pbc = flags['use_pbc'] masses = atomgroup.masses + com = atomgroup.center_of_mass(pbc=pbc, unwrap=unwrap, compound=compound) + if compound is not 'group': + com = (com * group.masses[:, None]).sum(axis=0) / group.masses.sum() + if pbc: - recenteredpos = (atomgroup.pack_into_box(inplace=False) - - atomgroup.center_of_mass(pbc=True)) + recenteredpos = (atomgroup.pack_into_box(inplace=False) - com) + elif unwrap: + recenteredpos = (atomgroup.unwrap(inplace=False) - com) else: - recenteredpos = (atomgroup.positions - - atomgroup.center_of_mass(pbc=False)) + recenteredpos = (atomgroup.positions - com) tensor = np.zeros((3, 3)) for x in range(recenteredpos.shape[0]): diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 08983e9f3c9..cb9ceb4024b 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -838,11 +838,43 @@ class TestUnwrapFlag(object): prec = 3 + @pytest.fixture() + def ag(self): + universe = mda.Universe(TRZ_psf, TRZ) + group = universe.residues[0:3] + group.wrap(inplace=True) + return group + + @pytest.fixture() + def ref_noUnwrap_residues(self): + return { + 'COG': np.array([[21.356, 28.52, 36.762], + [32.062, 36.16, 27.679], + [27.071, 29.997, 28.506]], dtype=np.float32), + 'COM': np.array([[21.286, 28.407, 36.629], + [31.931, 35.814, 27.916], + [26.817, 29.41, 29.05]]), + 'Asph': 0.53054563878, + + } + + @pytest.fixture() + def ref_Unwrap_residues(self): + return { + 'COG': np.array([[21.356, 41.685, 40.501], + [44.577, 43.312, 79.039], + [2.204, 27.722, 54.023]], dtype=np.float32), + 'COM': np.array([[20.815, 42.013, 39.802], + [44.918, 43.282, 79.325], + [2.045, 28.243, 54.127]], dtype=np.float32), + 'Asph': 0.2969491080, + } @pytest.fixture() def ref_noUnwrap(self): return { 'COG': np.array([5.1, 7.5, 7. ], dtype=np.float32), 'COM': np.array([6.48785, 7.5, 7.0], dtype=np.float32), + 'Asph': 1.0, } @pytest.fixture() @@ -850,6 +882,7 @@ def ref_Unwrap(self): return { 'COG': np.array([10.1, 7.5, 7. ], dtype=np.float32), 'COM': np.array([6.8616, 7.5, 7.], dtype=np.float32), + 'Asph': 1.0, } def test_default(self, ref_noUnwrap): @@ -860,6 +893,7 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(), ref_noUnwrap['COM'], self.prec) + assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -868,6 +902,8 @@ def test_UnWrapFlag(self, ref_Unwrap): group.masses = [100.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(unwrap=True), ref_Unwrap['COM'], self.prec) + assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) + class TestPBCFlag(object): From 7c97ee02618e983707512a786e529c431c94be9f Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 16:24:58 +0530 Subject: [PATCH 08/10] Adds tests --- package/MDAnalysis/core/topologyattrs.py | 20 ++++++++++++++----- .../MDAnalysisTests/core/test_atomgroup.py | 16 +++++++++++++-- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 566d70788c6..c2ee24713b0 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1028,7 +1028,8 @@ def shape_parameter(group, **kwargs): ('shape_parameter', shape_parameter)) @warn_if_not_unique - def asphericity(group, pbc=None): + @check_pbc_and_unwrap + def asphericity(group, pbc=None, unwrap=None, compound='group'): """Asphericity. See [Dima2004b]_ for background information. @@ -1039,6 +1040,10 @@ def asphericity(group, pbc=None): If ``True``, move all atoms within the primary unit cell before calculation. If ``None`` use value defined in MDAnalysis.core.flags['use_pbc'] + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + Which type of component to keep together during unwrapping. Note ---- @@ -1059,6 +1064,7 @@ def asphericity(group, pbc=None): .. versionadded:: 0.7.7 .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms @@ -1066,12 +1072,16 @@ def asphericity(group, pbc=None): pbc = flags['use_pbc'] masses = atomgroup.masses + com = atomgroup.center_of_mass(pbc=pbc, unwrap=unwrap, compound=compound) + if compound is not 'group': + com = (com * group.masses[:, None]).sum(axis=0) / group.masses.sum() + if pbc: - recenteredpos = (atomgroup.pack_into_box(inplace=False) - - atomgroup.center_of_mass(pbc=True)) + recenteredpos = (atomgroup.pack_into_box(inplace=False) - com) + elif unwrap: + recenteredpos = (atomgroup.unwrap(inplace=False) - com) else: - recenteredpos = (atomgroup.positions - - atomgroup.center_of_mass(pbc=False)) + recenteredpos = (atomgroup.positions - com) tensor = np.zeros((3, 3)) for x in range(recenteredpos.shape[0]): diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c5070a01022..91401d6d54a 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -848,16 +848,18 @@ def ag(self): @pytest.fixture() def ref_noUnwrap_residues(self): return { + 'COG': np.array([[21.356, 28.52 , 36.762], [32.062, 36.16 , 27.679], [27.071, 29.997, 28.506]], dtype=np.float32), 'COM': np.array([[21.286, 28.407, 36.629], [31.931, 35.814, 27.916], - [26.817, 29.41 , 29.05 ]]), + [26.817, 29.41, 29.05]]), 'MOI': np.array([ [7333.79167791, -211.8997285, -721.50785456], [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), + 'Asph': 0.02060121, } @@ -866,13 +868,15 @@ def ref_Unwrap_residues(self): return { 'COG': np.array([[21.356, 41.685, 40.501], [44.577, 43.312, 79.039], + [ 2.204, 27.722, 54.023]], dtype=np.float32), 'COM': np.array([[20.815, 42.013, 39.802], [44.918, 43.282, 79.325], [2.045, 28.243, 54.127]], dtype=np.float32), 'MOI': np.array([[16747.486, -1330.489, 2938.243], [-1330.489, 19315.253, 3306.212], - [ 2938.243, 3306.212, 8990.481]]) + [ 2938.243, 3306.212, 8990.481]]), + 'Asph': 0.2969491080, } @pytest.fixture() @@ -884,6 +888,7 @@ def ref_noUnwrap(self): [0.0, 0.0, 0.0], [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), + 'Asph': 1.0, } @pytest.fixture() @@ -895,17 +900,21 @@ def ref_Unwrap(self): [0.0, 0.0, 0.0], [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), + 'Asph': 1.0, } def test_default_residues(self, ag, ref_noUnwrap_residues): assert_almost_equal(ag.center_of_geometry(compound='residues'), ref_noUnwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(compound='residues'), ref_noUnwrap_residues['COM'], self.prec) assert_almost_equal(ag.moment_of_inertia(compound='residues'), ref_noUnwrap_residues['MOI'], self.prec) + assert_almost_equal(ag.asphericity(compound='residues'), ref_noUnwrap_residues['Asph'], self.prec) def test_UnWrapFlag_residues(self, ag, ref_Unwrap_residues): assert_almost_equal(ag.center_of_geometry(unwrap=True, compound='residues'), ref_Unwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(unwrap=True, compound='residues'), ref_Unwrap_residues['COM'], self.prec) assert_almost_equal(ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec) + assert_almost_equal(ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec) + def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -916,6 +925,7 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(), ref_noUnwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(), ref_noUnwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -925,6 +935,8 @@ def test_UnWrapFlag(self, ref_Unwrap): assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(unwrap=True), ref_Unwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(unwrap=True), ref_Unwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) + class TestPBCFlag(object): From 0a5e05900af685aaa8819b85b26b9424c652e097 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 12:33:56 +0530 Subject: [PATCH 09/10] Add unwrap to asphericity --- package/CHANGELOG | 1 + package/MDAnalysis/core/topologyattrs.py | 20 ++++++++++++++----- .../MDAnalysisTests/core/test_atomgroup.py | 17 +++++++++++----- 3 files changed, 28 insertions(+), 10 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 1471e70456e..10509a1f56f 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -60,6 +60,7 @@ Enhancements * added unwrap keyword to center_of_geometry (PR #2279) * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) + * added unwrap keyword to as Changes * added official support for Python 3.7 (PR #1963) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 566d70788c6..22a71302b30 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1028,7 +1028,8 @@ def shape_parameter(group, **kwargs): ('shape_parameter', shape_parameter)) @warn_if_not_unique - def asphericity(group, pbc=None): + @check_pbc_and_unwrap + def asphericity(group, pbc=None, unwrap=None, compound='group'): """Asphericity. See [Dima2004b]_ for background information. @@ -1039,6 +1040,10 @@ def asphericity(group, pbc=None): If ``True``, move all atoms within the primary unit cell before calculation. If ``None`` use value defined in MDAnalysis.core.flags['use_pbc'] + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + Which type of component to keep together during unwrapping. Note ---- @@ -1059,6 +1064,7 @@ def asphericity(group, pbc=None): .. versionadded:: 0.7.7 .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* and *compound* parameter """ atomgroup = group.atoms @@ -1066,12 +1072,16 @@ def asphericity(group, pbc=None): pbc = flags['use_pbc'] masses = atomgroup.masses + com = atomgroup.center_of_mass(pbc=pbc, unwrap=unwrap, compound=compound) + if compound is not 'group': + com = (com * group.masses[:, None]).sum(axis=0) / group.masses.sum() + if pbc: - recenteredpos = (atomgroup.pack_into_box(inplace=False) - - atomgroup.center_of_mass(pbc=True)) + recenteredpos = (atomgroup.pack_into_box(inplace=False) - com) + elif unwrap: + recenteredpos = (atomgroup.unwrap(inplace=False) - com) else: - recenteredpos = (atomgroup.positions - - atomgroup.center_of_mass(pbc=False)) + recenteredpos = (atomgroup.positions - com) tensor = np.zeros((3, 3)) for x in range(recenteredpos.shape[0]): diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c5070a01022..3c67a145f61 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -848,17 +848,17 @@ def ag(self): @pytest.fixture() def ref_noUnwrap_residues(self): return { - 'COG': np.array([[21.356, 28.52 , 36.762], - [32.062, 36.16 , 27.679], + 'COG': np.array([[21.356, 28.52, 36.762], + [32.062, 36.16, 27.679], [27.071, 29.997, 28.506]], dtype=np.float32), 'COM': np.array([[21.286, 28.407, 36.629], [31.931, 35.814, 27.916], - [26.817, 29.41 , 29.05 ]]), + [26.817, 29.41, 29.05]]), 'MOI': np.array([ [7333.79167791, -211.8997285, -721.50785456], [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), - + 'Asph': 0.02060121, } @pytest.fixture() @@ -872,7 +872,8 @@ def ref_Unwrap_residues(self): [2.045, 28.243, 54.127]], dtype=np.float32), 'MOI': np.array([[16747.486, -1330.489, 2938.243], [-1330.489, 19315.253, 3306.212], - [ 2938.243, 3306.212, 8990.481]]) + [ 2938.243, 3306.212, 8990.481]]), + 'Asph': 0.2969491080, } @pytest.fixture() @@ -884,6 +885,7 @@ def ref_noUnwrap(self): [0.0, 0.0, 0.0], [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), + 'Asph': 1.0, } @pytest.fixture() @@ -895,17 +897,20 @@ def ref_Unwrap(self): [0.0, 0.0, 0.0], [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), + 'Asph': 1.0, } def test_default_residues(self, ag, ref_noUnwrap_residues): assert_almost_equal(ag.center_of_geometry(compound='residues'), ref_noUnwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(compound='residues'), ref_noUnwrap_residues['COM'], self.prec) assert_almost_equal(ag.moment_of_inertia(compound='residues'), ref_noUnwrap_residues['MOI'], self.prec) + assert_almost_equal(ag.asphericity(compound='residues'), ref_noUnwrap_residues['Asph'], self.prec) def test_UnWrapFlag_residues(self, ag, ref_Unwrap_residues): assert_almost_equal(ag.center_of_geometry(unwrap=True, compound='residues'), ref_Unwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(unwrap=True, compound='residues'), ref_Unwrap_residues['COM'], self.prec) assert_almost_equal(ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec) + assert_almost_equal(ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -916,6 +921,7 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(), ref_noUnwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(), ref_noUnwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -925,6 +931,7 @@ def test_UnWrapFlag(self, ref_Unwrap): assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(unwrap=True), ref_Unwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(unwrap=True), ref_Unwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) class TestPBCFlag(object): From a5d8c9792a8fbf5ea407d1465ad8c3b076138a0b Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 12:42:09 +0530 Subject: [PATCH 10/10] Add CHANGELOG --- package/CHANGELOG | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 10509a1f56f..83a6dbe1228 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -60,7 +60,7 @@ Enhancements * added unwrap keyword to center_of_geometry (PR #2279) * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) - * added unwrap keyword to as + * added unwrap keyword to asphericity (PR #2290) Changes * added official support for Python 3.7 (PR #1963)