From 4819464618de18bc0460eb11bb6dfba132193f48 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Mon, 17 Jun 2019 13:15:32 +0530 Subject: [PATCH 01/46] 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/46] 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/46] 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/46] 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/46] 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/46] 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 1783fa0b7cf6fc50fbf8580cd599fafaf35a8274 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 07/46] Add unwrap to group.radius_of_gyration --- package/MDAnalysis/core/topologyattrs.py | 10 ++++++++- .../MDAnalysisTests/core/test_atomgroup.py | 22 +++++++++++++------ 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 566d70788c6..94c398cf4d9 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -953,15 +953,23 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) masses = atomgroup.masses + unwrap = kwargs.pop('unwrap', False) + compound = kwargs.pop('compound', 'group') + + 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() - com = atomgroup.center_of_mass(pbc=pbc) if pbc: recenteredpos = atomgroup.pack_into_box(inplace=False) - com + elif unwrap: + recenteredpos = atomgroup.unwrap(compound=compound) - com else: recenteredpos = atomgroup.positions - com diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c5070a01022..3e8c93b5f35 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -853,11 +853,12 @@ def ref_noUnwrap_residues(self): [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 ]], dtype=np.float32), 'MOI': np.array([ [7333.79167791, -211.8997285, -721.50785456], [-211.8997285, 7059.07470427, -91.32156884], - [-721.50785456, -91.32156884, 6509.31735029]]), + [-721.50785456, -91.32156884, 6509.31735029]], dtype=np.float32), + 'ROG': 27.713009, } @@ -870,9 +871,10 @@ def ref_Unwrap_residues(self): '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]]) + 'MOI': np.array([[16747.486, -1330.489, 2938.243], + [-1330.489, 19315.253, 3306.212], + [2938.243, 3306.212, 8990.481]], dtype=np.float32), + 'ROG': 40.686541, } @pytest.fixture() @@ -883,7 +885,8 @@ def ref_noUnwrap(self): 'MOI': np.array([ [0.0, 0.0, 0.0], [0.0, 98.6542, 0.0], - [0.0, 0.0, 98.65421327]]), + [0.0, 0.0, 98.65421327]], dtype=np.float32), + 'ROG': 0.9602093, } @pytest.fixture() @@ -894,18 +897,21 @@ def ref_Unwrap(self): 'MOI': np.array([ [0.0, 0.0, 0.0], [0.0, 132.673, 0.0], - [0.0, 0.0, 132.673]]), + [0.0, 0.0, 132.673]], dtype=np.float32), + 'ROG': 1.1135230, } 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.radius_of_gyration(compound='residues'), ref_noUnwrap_residues['ROG'], 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.radius_of_gyration(unwrap=True, compound='residues'), ref_Unwrap_residues['ROG'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -916,6 +922,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.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -925,6 +932,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.radius_of_gyration(unwrap=True), ref_Unwrap['ROG'], self.prec) class TestPBCFlag(object): From 09c8aaefc2a132358d4dc4db3e3c035d5acab07c Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 17 Jul 2018 13:03:01 +0100 Subject: [PATCH 08/46] added transformations based on alignto ; added translations fitting transformation --- .../MDAnalysis/transformations/__init__.py | 1 + package/MDAnalysis/transformations/fit.py | 197 ++++++++++++++++++ .../transformations/test_fit.py | 140 +++++++++++++ 3 files changed, 338 insertions(+) create mode 100644 package/MDAnalysis/transformations/fit.py create mode 100644 testsuite/MDAnalysisTests/transformations/test_fit.py diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index b1776527c59..e74ce842444 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -96,3 +96,4 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager +from .fit import alignto, fit_translation diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py new file mode 100644 index 00000000000..03cc88f1c11 --- /dev/null +++ b/package/MDAnalysis/transformations/fit.py @@ -0,0 +1,197 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +"""\ +Fitting transformations --- :mod:`MDAnalysis.transformations.fit` +================================================================= + +Translate and rotates the coordinates of a given trajectory to align +a given AtomGroup to a reference structure. + +""" +from __future__ import absolute_import + +import numpy as np +from functools import partial + +from ..analysis import align + +def alignto(ag, reference, select='all', weights=None, + subselection=None, tol_mass=0.1): + + """Perform a spatial superposition by minimizing the RMSD. + + Spatially align the group of atoms `ag` to `reference` by + doing a RMSD fit on `select` atoms. + + A detailed description on the way the fitting is performed can be found in + :func:`MDAnalysis.analysis.align.alignto` + + Example + ------- + + ..code-block::python + + ag = u.select_atoms("protein") + ref = mda.Universe("reference.pdb") + transform = MDAnalysis.transformations.alignto(ag, ref, wheights="mass") + u.trajectory.add_transformations(transform) + + Parameters + ---------- + ag : Universe or AtomGroup + structure to be aligned, a + :class:`~MDAnalysis.core.groups.AtomGroup` or a whole + :class:`~MDAnalysis.core.universe.Universe` + reference : Universe or AtomGroup + reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` + or a whole :class:`~MDAnalysis.core.universe.Universe` + select : str or dict or tuple, optional + The selection to operate on; can be one of: + + 1. any valid selection string for + :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` that + produces identical selections in `ag` and `reference`; or + + 2. a dictionary ``{'ag': sel1, 'reference': sel2}`` where *sel1* + and *sel2* are valid selection strings that are applied to + `ag` and `reference` respectively (the + :func:`MDAnalysis.analysis.align.fasta2select` function returns such + a dictionary based on a ClustalW_ or STAMP_ sequence alignment); or + + 3. a tuple ``(sel1, sel2)`` + + When using 2. or 3. with *sel1* and *sel2* then these selection strings + are applied to `ag` and `reference` respectively and should + generate *groups of equivalent atoms*. *sel1* and *sel2* can each also + be a *list of selection strings* to generate a + :class:`~MDAnalysis.core.groups.AtomGroup` with defined atom order as + described under :ref:`ordered-selections-label`). + weights : {"mass", ``None``} or array_like, optional + choose weights. With ``"mass"`` uses masses as weights; with ``None`` + weigh each atom equally. If a float array of the same length as + `ag` is provided, use each element of the `array_like` as a + weight for the corresponding atom in `mobile`. + tol_mass: float, optional + Reject match if the atomic masses for matched atoms differ by more than + `tol_mass`, default [0.1] + subselection : str or AtomGroup or None (optional) + Apply the transformation only to this selection. + + ``None`` [default] + Apply to ``mobile.universe.atoms`` (i.e., all atoms in the + context of the selection from `ag` such as the rest of a + protein, ligands and the surrounding water) + *selection-string* + Apply to ``mobile.select_atoms(selection-string)`` + :class:`~MDAnalysis.core.groups.AtomGroup` + Apply to the arbitrary group of atoms + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + + """ + + def wrapped(ts): + align.alignto(ag, reference, select, weights, subselection, tol_mass) + + return ts + + return wrapped + + +def fit_translation(ag, reference, plane=None, center_of="geometry"): + """ + Translates a given AtomGroup so that its center of geometry/mass matches + the respective center of the given reference. A plane can be given by the + user using the option `plane`, and will result in the removal of + the translation motions of the AtomGroup over that particular plane. + + Example + ------- + + ..code-block::python + + ag = u.select_atoms("protein") + ref = mda.Universe("reference.pdb") + transform = MDAnalysis.transformations.fit_translation(ag, ref, center_of="mass") + u.trajectory.add_transformations(transform) + + Parameters + ---------- + ag : Universe or AtomGroup + structure to translate, a + :class:`~MDAnalysis.core.groups.AtomGroup` or a whole + :class:`~MDAnalysis.core.universe.Universe` + reference : Universe or AtomGroup + reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole + :class:`~MDAnalysis.core.universe.Universe` + plane: str, optional + used to define the plane on which the translations will be removed. Defined as a + string of the plane. Suported planes are yz, xz and xy planes. + center_of: str, optional + used to choose the method of centering on the given atom group. Can be 'geometry' + or 'mass'. Default is "geometry". + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + """ + + if plane is not None: + if plane not in ('xy', 'yz', 'xz'): + raise ValueError('{} is not a valid plane'.format(plane)) + try: + if center_of == 'geometry': + ref_center = partial(reference.atoms.center_of_geometry) + ag_center = partial(ag.atoms.center_of_geometry) + elif center_of == 'mass': + ref_center = partial(reference.atoms.center_of_mass) + ag_center = partial(ag.atoms.center_of_mass) + else: + raise ValueError('{} is not a valid argument for "center_of"'.format(center_of)) + except AttributeError: + if center_of == 'mass': + raise AttributeError('Either {} or {} is not an Universe/AtomGroup object with masses'.format(ag, reference)) + else: + raise ValueError('Either {} or {} is not an Universe/AtomGroup object'.format(ag, reference)) + + reference = np.asarray(ref_center(), np.float32) + + def wrapped(ts): + center = np.asarray(ag_center(), np.float32) + if plane == 'yz': + vector = np.asarray([0, reference[1] - center[1], reference[2] - center[2]], np.float32) + if plane == 'xz': + vector = np.asarray([reference[0] - center[0], 0, reference[2] - center[2]], np.float32) + if plane == 'xy': + vector = np.asarray([reference[0] - center[0], reference[1] - center[1], 0], np.float32) + else: + vector = reference - center + ts.positions += vector + + return ts + + return wrapped + + \ No newline at end of file diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py new file mode 100644 index 00000000000..6b2718b8cd0 --- /dev/null +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -0,0 +1,140 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +from __future__ import absolute_import + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +from MDAnalysisTests import make_Universe +from MDAnalysis.transformations.fit import alignto, fit_translation + + +@pytest.fixture() +def test_universe(): + # make a test universe + test = make_Universe(('masses', ), trajectory=True) + ref = make_Universe(('masses', ), trajectory=True) + ref.trajectory.ts.positions += np.asarray([10, 10, 10], np.float32) + return test, ref + + +def test_alignto_coords(test_universe): + # when aligning the two universes do their coordinates become similar? + test_u = test_universe[0] + ref_u = test_universe[1] + alignto(test_u, ref_u, weights="mass")(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + + +def test_alignto_transformations_api(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + transform = alignto(test_u, ref_u, weights="mass") + test_u.trajectory.add_transformations(transform) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + + +def test_fit_translation_bad_ag(test_universe): + ts = test_universe[0].trajectory.ts + test_u = test_universe[0] + ref_u = test_universe[1] + # what happens if something other than an AtomGroup or Universe is given? + bad_ag = 1 + with pytest.raises(ValueError): + fit_translation(bad_ag, ref_u)(ts) + + +def test_fit_translation_bad_center(test_universe): + ts = test_universe[0].trajectory.ts + test_u = test_universe[0] + ref_u = test_universe[1] + # what happens if a bad string for center is given? + bad_center = "123" + with pytest.raises(ValueError): + fit_translation(test_u, ref_u, center_of = bad_center)(ts) + + +def test_fit_translation_bad_plane(test_universe): + ts = test_universe[0].trajectory.ts + test_u = test_universe[0] + ref_u = test_universe[1] + # what happens if a bad string for center is given? + bad_plane = "123" + with pytest.raises(ValueError): + fit_translation(test_u, ref_u, plane=bad_plane)(ts) + + +def test_fit_translation_no_masses(test_universe): + ts = test_universe[0].trajectory.ts + test_u = test_universe[0] + # create a universe without masses + ref_u = make_Universe() + # what happens Universe without masses is given? + with pytest.raises(AttributeError): + fit_translation(test_u, ref_u, center_of="mass")(ts) + + +def test_fit_translation_no_options(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + fit_translation(test_u, ref_u)(test_u.trajectory.ts) + # what happens when no options are passed? + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + + +def test_fit_translation_com(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + fit_translation(test_u, ref_u, center_of="mass")(test_u.trajectory.ts) + # what happens when the center o mass is used? + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + + +def test_fit_translation_plane(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + # translate the test universe on the x and y coordinates only + fit_translation(test_u, ref_u, plane="xy")(test_u.trajectory.ts) + # the reference is 10 angstrom in the z coordinate above the test universe + shiftz = np.asanyarray([0, 0, -10], np.float32) + ref_coordinates = ref_u.trajectory.ts.positions + shiftz + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_coordinates, decimal=6) + + +def test_fit_translation_all_options(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + # translate the test universe on the x and y coordinates only + fit_translation(test_u, ref_u, plane="xy", center_of="mass")(test_u.trajectory.ts) + # the reference is 10 angstrom in the z coordinate above the test universe + shiftz = np.asanyarray([0, 0, -10], np.float32) + ref_coordinates = ref_u.trajectory.ts.positions + shiftz + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_coordinates, decimal=6) + + +def test_fit_translation_transformations_api(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + transform = fit_translation(test_u, ref_u) + test_u.trajectory.add_transformations(transform) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) From ce3ac2ba55c741c3130a11b20ff2436718725621 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 17 Jul 2018 23:12:35 +0100 Subject: [PATCH 09/46] added rotation+translation fitting with plane support --- package/MDAnalysis/transformations/fit.py | 136 +++++++++++------- .../transformations/test_fit.py | 9 +- 2 files changed, 93 insertions(+), 52 deletions(-) diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py index 03cc88f1c11..74f037471d8 100644 --- a/package/MDAnalysis/transformations/fit.py +++ b/package/MDAnalysis/transformations/fit.py @@ -33,10 +33,11 @@ import numpy as np from functools import partial -from ..analysis import align +from ..analysis import align +from ..lib.util import get_weights +from ..lib.transformations import euler_from_matrix, euler_matrix -def alignto(ag, reference, select='all', weights=None, - subselection=None, tol_mass=0.1): +def alignto(ag, reference, weights=None): """Perform a spatial superposition by minimizing the RMSD. @@ -65,46 +66,11 @@ def alignto(ag, reference, select='all', weights=None, reference : Universe or AtomGroup reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole :class:`~MDAnalysis.core.universe.Universe` - select : str or dict or tuple, optional - The selection to operate on; can be one of: - - 1. any valid selection string for - :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` that - produces identical selections in `ag` and `reference`; or - - 2. a dictionary ``{'ag': sel1, 'reference': sel2}`` where *sel1* - and *sel2* are valid selection strings that are applied to - `ag` and `reference` respectively (the - :func:`MDAnalysis.analysis.align.fasta2select` function returns such - a dictionary based on a ClustalW_ or STAMP_ sequence alignment); or - - 3. a tuple ``(sel1, sel2)`` - - When using 2. or 3. with *sel1* and *sel2* then these selection strings - are applied to `ag` and `reference` respectively and should - generate *groups of equivalent atoms*. *sel1* and *sel2* can each also - be a *list of selection strings* to generate a - :class:`~MDAnalysis.core.groups.AtomGroup` with defined atom order as - described under :ref:`ordered-selections-label`). weights : {"mass", ``None``} or array_like, optional choose weights. With ``"mass"`` uses masses as weights; with ``None`` weigh each atom equally. If a float array of the same length as `ag` is provided, use each element of the `array_like` as a weight for the corresponding atom in `mobile`. - tol_mass: float, optional - Reject match if the atomic masses for matched atoms differ by more than - `tol_mass`, default [0.1] - subselection : str or AtomGroup or None (optional) - Apply the transformation only to this selection. - - ``None`` [default] - Apply to ``mobile.universe.atoms`` (i.e., all atoms in the - context of the selection from `ag` such as the rest of a - protein, ligands and the surrounding water) - *selection-string* - Apply to ``mobile.select_atoms(selection-string)`` - :class:`~MDAnalysis.core.groups.AtomGroup` - Apply to the arbitrary group of atoms Returns ------- @@ -113,7 +79,7 @@ def alignto(ag, reference, select='all', weights=None, """ def wrapped(ts): - align.alignto(ag, reference, select, weights, subselection, tol_mass) + align.alignto(ag, reference, weights=weights) return ts @@ -121,20 +87,22 @@ def wrapped(ts): def fit_translation(ag, reference, plane=None, center_of="geometry"): - """ - Translates a given AtomGroup so that its center of geometry/mass matches + + """Translates a given AtomGroup so that its center of geometry/mass matches the respective center of the given reference. A plane can be given by the user using the option `plane`, and will result in the removal of the translation motions of the AtomGroup over that particular plane. Example ------- + Removing the translations of a given AtomGroup `ag` on the XY plane by fitting + its center of mass to the center of mass of a reference `ref`: ..code-block::python ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") - transform = MDAnalysis.transformations.fit_translation(ag, ref, center_of="mass") + transform = MDAnalysis.transformations.fit_translation(ag, ref, plane="xy", center_of="mass") u.trajectory.add_transformations(transform) Parameters @@ -161,6 +129,8 @@ def fit_translation(ag, reference, plane=None, center_of="geometry"): if plane is not None: if plane not in ('xy', 'yz', 'xz'): raise ValueError('{} is not a valid plane'.format(plane)) + axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} + plane = axes[plane] try: if center_of == 'geometry': ref_center = partial(reference.atoms.center_of_geometry) @@ -180,18 +150,82 @@ def fit_translation(ag, reference, plane=None, center_of="geometry"): def wrapped(ts): center = np.asarray(ag_center(), np.float32) - if plane == 'yz': - vector = np.asarray([0, reference[1] - center[1], reference[2] - center[2]], np.float32) - if plane == 'xz': - vector = np.asarray([reference[0] - center[0], 0, reference[2] - center[2]], np.float32) - if plane == 'xy': - vector = np.asarray([reference[0] - center[0], reference[1] - center[1], 0], np.float32) - else: - vector = reference - center + vector = reference - center + if plane: + vector[plane] = 0 ts.positions += vector return ts return wrapped - \ No newline at end of file + +def fit_rot_trans(ag, reference, plane=None, weights=None): + """Perform a spatial superposition by minimizing the RMSD. + + Spatially align the group of atoms `ag` to `reference` by doing a RMSD + fit on `select` atoms. + + This fit works as a way o remove translations and rotations of a given + AtomGroup in a trajectory. A plane can be given using the flag `plane` + so that only translations and rotations in that particular plane are + removed. This is useful for protein-membrane systems to where the membrane + must remain in the same orientation. + + Example + ------- + Removing the translations and rotations of a given AtomGroup `ag` on the XY plane + by fitting it to a reference `ref`, using the masses as weights for the RMSD fit. + + ..code-block::python + + ag = u.select_atoms("protein") + ref = mda.Universe("reference.pdb") + transform = MDAnalysis.transformations.fit_rot_trans(ag, ref, plane="xy", weights="mass") + u.trajectory.add_transformations(transform) + + Parameter + --------- + ag : Universe or AtomGroup + structure to translate, a + :class:`~MDAnalysis.core.groups.AtomGroup` or a whole + :class:`~MDAnalysis.core.universe.Universe` + reference : Universe or AtomGroup + reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole + :class:`~MDAnalysis.core.universe.Universe` + plane: str, optional + used to define the plane on which the translations will be removed. Defined as a + string of the plane. Suported planes are "yz", "xz" and "xy" planes. + weights : {"mass", ``None``} or array_like, optional + choose weights. With ``"mass"`` uses masses as weights; with ``None`` + weigh each atom equally. If a float array of the same length as + `ag` is provided, use each element of the `array_like` as a + weight for the corresponding atom in `mobile`. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + """ + if plane is not None: + if plane not in ('xy', 'yz', 'xz'): + raise ValueError('{} is not a valid plane'.format(plane)) + axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} + plane = axes[plane] + reference, ag = align.get_matching_atoms(reference.atoms, ag.atoms) + weights = align.get_weights(reference.atoms, weights) + ref_com = reference.atoms.center(weights) + ref_coordinates = reference.atoms.positions - ref_com + def wrapped(ts): + mobile_com = ag.atoms.center(weights) + mobile_coordinates = ts.positions - mobile_com + rotation, dump = align.rotation_matrix(mobile_coordinates, ref_coordinates, weights=weights) + if plane: + euler_angs = euler_from_matrix(rotation, axes='sxyz') + euler_angs[plane] = 0 + rotation = euler_matrix(euler_angs[0], euler_angs[1], euler_angs[2], axes='sxyz') + ts.positions = ts.positions + (ref_com - mobile_com) + ts.positions = np.dot(ts.positions, rotation) + + return ts + + return wrapped diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py index 6b2718b8cd0..155e848e7a2 100644 --- a/testsuite/MDAnalysisTests/transformations/test_fit.py +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -26,7 +26,7 @@ from numpy.testing import assert_array_almost_equal from MDAnalysisTests import make_Universe -from MDAnalysis.transformations.fit import alignto, fit_translation +from MDAnalysis.transformations.fit import alignto, fit_translation, fit_rot_trans @pytest.fixture() @@ -138,3 +138,10 @@ def test_fit_translation_transformations_api(test_universe): transform = fit_translation(test_u, ref_u) test_u.trajectory.add_transformations(transform) assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + + +def test_fit_rot_trans_coords(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + fit_rot_trans(test_u, ref_u, weights="mass")(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) From 5a702d3caaa83ada1e5b6baee5d3628bf083387a Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 17 Jul 2018 23:14:37 +0100 Subject: [PATCH 10/46] updated __init__ --- package/MDAnalysis/transformations/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index e74ce842444..ba21c332d1f 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -96,4 +96,4 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager -from .fit import alignto, fit_translation +from .fit import alignto, fit_translation, fit_rot_trans From a0390ac89149aa87cb8529b55cb0f8545189660b Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Wed, 18 Jul 2018 20:23:24 +0100 Subject: [PATCH 11/46] corrected some errors in fit_rot_trans ; added more tests --- package/MDAnalysis/transformations/fit.py | 26 +++-- .../transformations/test_fit.py | 99 +++++++++++++++++-- 2 files changed, 108 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py index 74f037471d8..b5e2e29c11a 100644 --- a/package/MDAnalysis/transformations/fit.py +++ b/package/MDAnalysis/transformations/fit.py @@ -211,20 +211,30 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): raise ValueError('{} is not a valid plane'.format(plane)) axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} plane = axes[plane] - reference, ag = align.get_matching_atoms(reference.atoms, ag.atoms) - weights = align.get_weights(reference.atoms, weights) - ref_com = reference.atoms.center(weights) - ref_coordinates = reference.atoms.positions - ref_com + try: + if ag.atoms.n_residues != reference.atoms.n_residues: + raise ValueError("{} and {} have mismatched number of residues".format(ag,reference)) + except AttributeError: + raise AttributeError("{} or {} is not valid Universe/AtomGroup".format(ag,reference)) + ref, mobile = align.get_matching_atoms(reference.atoms, ag.atoms) + try: + weights = align.get_weights(ref.atoms, weights=weights) + except (ValueError, TypeError): + raise TypeError("weights must be {'mass', None} or an iterable of the " + "same size as the atomgroup.") + ref_com = ref.center(weights) + ref_coordinates = ref.atoms.positions - ref_com def wrapped(ts): - mobile_com = ag.atoms.center(weights) - mobile_coordinates = ts.positions - mobile_com + mobile_com = mobile.atoms.center(weights) + mobile_coordinates = mobile.atoms.positions - mobile_com rotation, dump = align.rotation_matrix(mobile_coordinates, ref_coordinates, weights=weights) if plane: euler_angs = euler_from_matrix(rotation, axes='sxyz') euler_angs[plane] = 0 rotation = euler_matrix(euler_angs[0], euler_angs[1], euler_angs[2], axes='sxyz') - ts.positions = ts.positions + (ref_com - mobile_com) - ts.positions = np.dot(ts.positions, rotation) + ts.positions = ts.positions - mobile_com + ts.positions = np.dot(ts.positions, rotation.T) + ts.positions = ts.positions + ref_com return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py index 155e848e7a2..1a2ed37ef5b 100644 --- a/testsuite/MDAnalysisTests/transformations/test_fit.py +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -27,6 +27,7 @@ from MDAnalysisTests import make_Universe from MDAnalysis.transformations.fit import alignto, fit_translation, fit_rot_trans +from MDAnalysis.lib.transformations import rotation_matrix @pytest.fixture() @@ -34,7 +35,7 @@ def test_universe(): # make a test universe test = make_Universe(('masses', ), trajectory=True) ref = make_Universe(('masses', ), trajectory=True) - ref.trajectory.ts.positions += np.asarray([10, 10, 10], np.float32) + ref.atoms.positions += np.asarray([10, 10, 10], np.float32) return test, ref @@ -42,16 +43,26 @@ def test_alignto_coords(test_universe): # when aligning the two universes do their coordinates become similar? test_u = test_universe[0] ref_u = test_universe[1] - alignto(test_u, ref_u, weights="mass")(test_u.trajectory.ts) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + ref_com = ref_u.atoms.center(None) + ref_u.trajectory.ts.positions -= ref_com + R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_u.trajectory.ts.positions += ref_com + alignto(test_u, ref_u, weights=None)(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) def test_alignto_transformations_api(test_universe): test_u = test_universe[0] ref_u = test_universe[1] - transform = alignto(test_u, ref_u, weights="mass") + ref_com = ref_u.atoms.center(None) + ref_u.trajectory.ts.positions -= ref_com + R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_u.trajectory.ts.positions += ref_com + transform = alignto(test_u, ref_u, weights=None) test_u.trajectory.add_transformations(transform) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) def test_fit_translation_bad_ag(test_universe): @@ -138,10 +149,80 @@ def test_fit_translation_transformations_api(test_universe): transform = fit_translation(test_u, ref_u) test_u.trajectory.add_transformations(transform) assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) - -def test_fit_rot_trans_coords(test_universe): + +def test_fit_rot_trans_bad_universe(test_universe): test_u = test_universe[0] + bad_u = 1 + ref_u= bad_u + with pytest.raises(AttributeError): + fit_rot_trans(test_u, ref_u)(test_u.trajectory.ts) + + +def test_fit_rot_trans_shorter_universe(test_universe): ref_u = test_universe[1] - fit_rot_trans(test_u, ref_u, weights="mass")(test_u.trajectory.ts) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) + bad_u =test_universe[0].atoms[0:5] + test_u= bad_u + with pytest.raises(ValueError): + fit_rot_trans(test_u, ref_u)(test_u.trajectory.ts) + + +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_fit_rot_trans_bad_weights(test_universe, weights): + test_u = test_universe[0] + ref_u = test_universe[1] + bad_weights = weights + with pytest.raises(TypeError): + fit_rot_trans(test_u, ref_u, weights=bad_weights)(test_u.trajectory.ts) + + +@pytest.mark.parametrize('plane', ( + " ", + "totallynotaplane", + "xyz", + 123456789, + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_fit_rot_trans_bad_plane(test_universe, plane): + test_u = test_universe[0] + ref_u = test_universe[1] + bad_plane = plane + with pytest.raises(ValueError): + fit_rot_trans(test_u, ref_u, plane=bad_plane)(test_u.trajectory.ts) + + +def test_fit_rot_trans_no_options(test_universe): + test_u = test_universe[0] + ref_u = test_universe[1] + ref_com = ref_u.atoms.center(None) + ref_u.trajectory.ts.positions -= ref_com + R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_u.trajectory.ts.positions += ref_com + fit_rot_trans(test_u, ref_u)(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) + + +def test_fit_rot_trans_plane(test_universe): + # the reference is rotated in the x axis so removing the translations and rotations + # in the yz plane should return the same as the fitting without specifying a plane + test_u = test_universe[0] + ref_u = test_universe[1] + ref_com = ref_u.atoms.center(None) + ref_u.trajectory.ts.positions -= ref_com + R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_u.trajectory.ts.positions += ref_com + fit_rot_trans(test_u, ref_u, plane="yz")(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file From 9f7791e6a1809e9b5f8f8059cab4482a4d723be2 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 24 Jul 2018 22:55:58 +0100 Subject: [PATCH 12/46] corrected plane rot+trans fitting; removed alignto --- .../MDAnalysis/transformations/__init__.py | 2 +- package/MDAnalysis/transformations/fit.py | 131 ++++-------- .../transformations/test_fit.py | 201 ++++++++++-------- 3 files changed, 156 insertions(+), 178 deletions(-) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index ba21c332d1f..31c7241889c 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -96,4 +96,4 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager -from .fit import alignto, fit_translation, fit_rot_trans +from .fit import fit_translation, fit_rot_trans diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py index b5e2e29c11a..3fe8595fe25 100644 --- a/package/MDAnalysis/transformations/fit.py +++ b/package/MDAnalysis/transformations/fit.py @@ -37,56 +37,8 @@ from ..lib.util import get_weights from ..lib.transformations import euler_from_matrix, euler_matrix -def alignto(ag, reference, weights=None): - - """Perform a spatial superposition by minimizing the RMSD. - - Spatially align the group of atoms `ag` to `reference` by - doing a RMSD fit on `select` atoms. - - A detailed description on the way the fitting is performed can be found in - :func:`MDAnalysis.analysis.align.alignto` - - Example - ------- - - ..code-block::python - - ag = u.select_atoms("protein") - ref = mda.Universe("reference.pdb") - transform = MDAnalysis.transformations.alignto(ag, ref, wheights="mass") - u.trajectory.add_transformations(transform) - - Parameters - ---------- - ag : Universe or AtomGroup - structure to be aligned, a - :class:`~MDAnalysis.core.groups.AtomGroup` or a whole - :class:`~MDAnalysis.core.universe.Universe` - reference : Universe or AtomGroup - reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` - or a whole :class:`~MDAnalysis.core.universe.Universe` - weights : {"mass", ``None``} or array_like, optional - choose weights. With ``"mass"`` uses masses as weights; with ``None`` - weigh each atom equally. If a float array of the same length as - `ag` is provided, use each element of the `array_like` as a - weight for the corresponding atom in `mobile`. - - Returns - ------- - MDAnalysis.coordinates.base.Timestep - - """ - - def wrapped(ts): - align.alignto(ag, reference, weights=weights) - - return ts - - return wrapped - -def fit_translation(ag, reference, plane=None, center_of="geometry"): +def fit_translation(ag, reference, plane=None, weights=None): """Translates a given AtomGroup so that its center of geometry/mass matches the respective center of the given reference. A plane can be given by the @@ -102,7 +54,7 @@ def fit_translation(ag, reference, plane=None, center_of="geometry"): ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") - transform = MDAnalysis.transformations.fit_translation(ag, ref, plane="xy", center_of="mass") + transform = mda.transformations.fit_translation(ag, ref, plane="xy", center_of="mass") u.trajectory.add_transformations(transform) Parameters @@ -117,9 +69,11 @@ def fit_translation(ag, reference, plane=None, center_of="geometry"): plane: str, optional used to define the plane on which the translations will be removed. Defined as a string of the plane. Suported planes are yz, xz and xy planes. - center_of: str, optional - used to choose the method of centering on the given atom group. Can be 'geometry' - or 'mass'. Default is "geometry". + weights : {"mass", ``None``} or array_like, optional + choose weights. With ``"mass"`` uses masses as weights; with ``None`` + weigh each atom equally. If a float array of the same length as + `ag` is provided, use each element of the `array_like` as a + weight for the corresponding atom in `ag`. Returns ------- @@ -127,31 +81,29 @@ def fit_translation(ag, reference, plane=None, center_of="geometry"): """ if plane is not None: - if plane not in ('xy', 'yz', 'xz'): - raise ValueError('{} is not a valid plane'.format(plane)) axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} - plane = axes[plane] + try: + plane = axes[plane] + except (TypeError, KeyError): + raise ValueError('{} is not a valid plane'.format(plane)) try: - if center_of == 'geometry': - ref_center = partial(reference.atoms.center_of_geometry) - ag_center = partial(ag.atoms.center_of_geometry) - elif center_of == 'mass': - ref_center = partial(reference.atoms.center_of_mass) - ag_center = partial(ag.atoms.center_of_mass) - else: - raise ValueError('{} is not a valid argument for "center_of"'.format(center_of)) + if ag.atoms.n_residues != reference.atoms.n_residues: + raise ValueError("{} and {} have mismatched number of residues".format(ag,reference)) except AttributeError: - if center_of == 'mass': - raise AttributeError('Either {} or {} is not an Universe/AtomGroup object with masses'.format(ag, reference)) - else: - raise ValueError('Either {} or {} is not an Universe/AtomGroup object'.format(ag, reference)) + raise AttributeError("{} or {} is not valid Universe/AtomGroup".format(ag,reference)) + ref, mobile = align.get_matching_atoms(reference.atoms, ag.atoms) + try: + weights = align.get_weights(ref.atoms, weights=weights) + except (ValueError, TypeError): + raise ValueError("weights must be {'mass', None} or an iterable of the " + "same size as the atomgroup.") - reference = np.asarray(ref_center(), np.float32) + ref_com = np.asarray(ref.center(weights), np.float32) def wrapped(ts): - center = np.asarray(ag_center(), np.float32) - vector = reference - center - if plane: + mobile_com = np.asarray(mobile.atoms.center(weights), np.float32) + vector = ref_com - mobile_com + if plane is not None: vector[plane] = 0 ts.positions += vector @@ -181,7 +133,7 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") - transform = MDAnalysis.transformations.fit_rot_trans(ag, ref, plane="xy", weights="mass") + transform = mda.transformations.fit_rot_trans(ag, ref, plane="xy", weights="mass") u.trajectory.add_transformations(transform) Parameter @@ -194,23 +146,24 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole :class:`~MDAnalysis.core.universe.Universe` plane: str, optional - used to define the plane on which the translations will be removed. Defined as a - string of the plane. Suported planes are "yz", "xz" and "xy" planes. + used to define the plane on which the rotations and translations will be removed. + Defined as a string of the plane. Suported planes are "yz", "xz" and "xy" planes. weights : {"mass", ``None``} or array_like, optional choose weights. With ``"mass"`` uses masses as weights; with ``None`` weigh each atom equally. If a float array of the same length as `ag` is provided, use each element of the `array_like` as a - weight for the corresponding atom in `mobile`. + weight for the corresponding atom in `ag`. Returns ------- MDAnalysis.coordinates.base.Timestep """ if plane is not None: - if plane not in ('xy', 'yz', 'xz'): - raise ValueError('{} is not a valid plane'.format(plane)) axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} - plane = axes[plane] + try: + plane = axes[plane] + except (TypeError, KeyError): + raise ValueError('{} is not a valid plane'.format(plane)) try: if ag.atoms.n_residues != reference.atoms.n_residues: raise ValueError("{} and {} have mismatched number of residues".format(ag,reference)) @@ -220,21 +173,27 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): try: weights = align.get_weights(ref.atoms, weights=weights) except (ValueError, TypeError): - raise TypeError("weights must be {'mass', None} or an iterable of the " + raise ValueError("weights must be {'mass', None} or an iterable of the " "same size as the atomgroup.") ref_com = ref.center(weights) - ref_coordinates = ref.atoms.positions - ref_com + ref_coordinates = ref.atoms.positions - ref_com + def wrapped(ts): mobile_com = mobile.atoms.center(weights) mobile_coordinates = mobile.atoms.positions - mobile_com rotation, dump = align.rotation_matrix(mobile_coordinates, ref_coordinates, weights=weights) - if plane: - euler_angs = euler_from_matrix(rotation, axes='sxyz') - euler_angs[plane] = 0 - rotation = euler_matrix(euler_angs[0], euler_angs[1], euler_angs[2], axes='sxyz') + vector = ref_com + if plane is not None: + matrix = np.r_[rotation, np.zeros(3).reshape(1,3)] + matrix = np.c_[matrix, np.zeros(4)] + euler_angs = np.asarray(euler_from_matrix(matrix, axes='sxyz'), np.float32) + for i in range(0, euler_angs.size): + euler_angs[i] = ( euler_angs[plane] if i == plane else 0) + rotation = euler_matrix(euler_angs[0], euler_angs[1], euler_angs[2], axes='sxyz')[:3, :3] + vector[plane] = mobile_com[plane] ts.positions = ts.positions - mobile_com ts.positions = np.dot(ts.positions, rotation.T) - ts.positions = ts.positions + ref_com + ts.positions = ts.positions + vector return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py index 1a2ed37ef5b..f5bd0913a08 100644 --- a/testsuite/MDAnalysisTests/transformations/test_fit.py +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -26,12 +26,12 @@ from numpy.testing import assert_array_almost_equal from MDAnalysisTests import make_Universe -from MDAnalysis.transformations.fit import alignto, fit_translation, fit_rot_trans +from MDAnalysis.transformations.fit import fit_translation, fit_rot_trans from MDAnalysis.lib.transformations import rotation_matrix @pytest.fixture() -def test_universe(): +def fit_universe(): # make a test universe test = make_Universe(('masses', ), trajectory=True) ref = make_Universe(('masses', ), trajectory=True) @@ -39,91 +39,89 @@ def test_universe(): return test, ref -def test_alignto_coords(test_universe): - # when aligning the two universes do their coordinates become similar? - test_u = test_universe[0] - ref_u = test_universe[1] - ref_com = ref_u.atoms.center(None) - ref_u.trajectory.ts.positions -= ref_com - R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] - ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) - ref_u.trajectory.ts.positions += ref_com - alignto(test_u, ref_u, weights=None)(test_u.trajectory.ts) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) - - -def test_alignto_transformations_api(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] - ref_com = ref_u.atoms.center(None) - ref_u.trajectory.ts.positions -= ref_com - R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] - ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) - ref_u.trajectory.ts.positions += ref_com - transform = alignto(test_u, ref_u, weights=None) - test_u.trajectory.add_transformations(transform) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) - - -def test_fit_translation_bad_ag(test_universe): - ts = test_universe[0].trajectory.ts - test_u = test_universe[0] - ref_u = test_universe[1] +@pytest.mark.parametrize('universe', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotanag', + 1) +) +def test_fit_translation_bad_ag(fit_universe, universe): + ts = fit_universe[0].trajectory.ts + test_u = fit_universe[0] + ref_u = fit_universe[1] # what happens if something other than an AtomGroup or Universe is given? - bad_ag = 1 - with pytest.raises(ValueError): - fit_translation(bad_ag, ref_u)(ts) + with pytest.raises(AttributeError): + fit_translation(universe, ref_u)(ts) -def test_fit_translation_bad_center(test_universe): - ts = test_universe[0].trajectory.ts - test_u = test_universe[0] - ref_u = test_universe[1] +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_fit_translation_bad_weights(fit_universe, weights): + ts = fit_universe[0].trajectory.ts + test_u = fit_universe[0] + ref_u = fit_universe[1] # what happens if a bad string for center is given? - bad_center = "123" with pytest.raises(ValueError): - fit_translation(test_u, ref_u, center_of = bad_center)(ts) + fit_translation(test_u, ref_u, weights=weights)(ts) -def test_fit_translation_bad_plane(test_universe): - ts = test_universe[0].trajectory.ts - test_u = test_universe[0] - ref_u = test_universe[1] +@pytest.mark.parametrize('plane', ( + 1, + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + "xyz", + "notaplane") +) +def test_fit_translation_bad_plane(fit_universe, plane): + ts = fit_universe[0].trajectory.ts + test_u = fit_universe[0] + ref_u = fit_universe[1] # what happens if a bad string for center is given? - bad_plane = "123" with pytest.raises(ValueError): - fit_translation(test_u, ref_u, plane=bad_plane)(ts) + fit_translation(test_u, ref_u, plane=plane)(ts) -def test_fit_translation_no_masses(test_universe): - ts = test_universe[0].trajectory.ts - test_u = test_universe[0] +def test_fit_translation_no_masses(fit_universe): + ts = fit_universe[0].trajectory.ts + test_u = fit_universe[0] # create a universe without masses ref_u = make_Universe() # what happens Universe without masses is given? with pytest.raises(AttributeError): - fit_translation(test_u, ref_u, center_of="mass")(ts) + fit_translation(test_u, ref_u, weights="mass")(ts) -def test_fit_translation_no_options(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_translation_no_options(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] fit_translation(test_u, ref_u)(test_u.trajectory.ts) # what happens when no options are passed? assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) -def test_fit_translation_com(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] - fit_translation(test_u, ref_u, center_of="mass")(test_u.trajectory.ts) +def test_fit_translation_com(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] + fit_translation(test_u, ref_u, weights="mass")(test_u.trajectory.ts) # what happens when the center o mass is used? assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) -def test_fit_translation_plane(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_translation_plane(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] # translate the test universe on the x and y coordinates only fit_translation(test_u, ref_u, plane="xy")(test_u.trajectory.ts) # the reference is 10 angstrom in the z coordinate above the test universe @@ -132,36 +130,45 @@ def test_fit_translation_plane(test_universe): assert_array_almost_equal(test_u.trajectory.ts.positions, ref_coordinates, decimal=6) -def test_fit_translation_all_options(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_translation_all_options(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] # translate the test universe on the x and y coordinates only - fit_translation(test_u, ref_u, plane="xy", center_of="mass")(test_u.trajectory.ts) + fit_translation(test_u, ref_u, plane="xy", weights="mass")(test_u.trajectory.ts) # the reference is 10 angstrom in the z coordinate above the test universe shiftz = np.asanyarray([0, 0, -10], np.float32) ref_coordinates = ref_u.trajectory.ts.positions + shiftz assert_array_almost_equal(test_u.trajectory.ts.positions, ref_coordinates, decimal=6) -def test_fit_translation_transformations_api(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_translation_transformations_api(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] transform = fit_translation(test_u, ref_u) test_u.trajectory.add_transformations(transform) assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) -def test_fit_rot_trans_bad_universe(test_universe): - test_u = test_universe[0] - bad_u = 1 - ref_u= bad_u +@pytest.mark.parametrize('universe', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotanag', + 1) +) +def test_fit_rot_trans_bad_universe(fit_universe, universe): + test_u = fit_universe[0] + ref_u= universe with pytest.raises(AttributeError): fit_rot_trans(test_u, ref_u)(test_u.trajectory.ts) -def test_fit_rot_trans_shorter_universe(test_universe): - ref_u = test_universe[1] - bad_u =test_universe[0].atoms[0:5] +def test_fit_rot_trans_shorter_universe(fit_universe): + ref_u = fit_universe[1] + bad_u =fit_universe[0].atoms[0:5] test_u= bad_u with pytest.raises(ValueError): fit_rot_trans(test_u, ref_u)(test_u.trajectory.ts) @@ -176,11 +183,11 @@ def test_fit_rot_trans_shorter_universe(test_universe): np.array([0, 1, 2, 3, 4]), np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) ) -def test_fit_rot_trans_bad_weights(test_universe, weights): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_rot_trans_bad_weights(fit_universe, weights): + test_u = fit_universe[0] + ref_u = fit_universe[1] bad_weights = weights - with pytest.raises(TypeError): + with pytest.raises(ValueError): fit_rot_trans(test_u, ref_u, weights=bad_weights)(test_u.trajectory.ts) @@ -194,17 +201,16 @@ def test_fit_rot_trans_bad_weights(test_universe, weights): np.array([0, 1, 2, 3, 4]), np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) ) -def test_fit_rot_trans_bad_plane(test_universe, plane): - test_u = test_universe[0] - ref_u = test_universe[1] - bad_plane = plane +def test_fit_rot_trans_bad_plane(fit_universe, plane): + test_u = fit_universe[0] + ref_u = fit_universe[1] with pytest.raises(ValueError): - fit_rot_trans(test_u, ref_u, plane=bad_plane)(test_u.trajectory.ts) + fit_rot_trans(test_u, ref_u, plane=plane)(test_u.trajectory.ts) -def test_fit_rot_trans_no_options(test_universe): - test_u = test_universe[0] - ref_u = test_universe[1] +def test_fit_rot_trans_no_options(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] ref_com = ref_u.atoms.center(None) ref_u.trajectory.ts.positions -= ref_com R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] @@ -214,15 +220,28 @@ def test_fit_rot_trans_no_options(test_universe): assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) -def test_fit_rot_trans_plane(test_universe): +def test_fit_rot_trans_plane(fit_universe): # the reference is rotated in the x axis so removing the translations and rotations # in the yz plane should return the same as the fitting without specifying a plane - test_u = test_universe[0] - ref_u = test_universe[1] + test_u = fit_universe[0] + ref_u = fit_universe[1] ref_com = ref_u.atoms.center(None) ref_u.trajectory.ts.positions -= ref_com R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) ref_u.trajectory.ts.positions += ref_com fit_rot_trans(test_u, ref_u, plane="yz")(test_u.trajectory.ts) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) + + +def test_fit_rot_trans_transformations_api(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1] + ref_com = ref_u.atoms.center(None) + ref_u.trajectory.ts.positions -= ref_com + R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_u.trajectory.ts.positions += ref_com + transform = fit_rot_trans(test_u, ref_u) + test_u.trajectory.add_transformations(transform) +# assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file From 408243a27274f80667af85fb49b4d8a7ed0ec880 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Thu, 26 Jul 2018 18:06:45 +0100 Subject: [PATCH 13/46] added plane tests --- .../transformations/test_fit.py | 39 ++++++++++++++----- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py index f5bd0913a08..2f9b3008c33 100644 --- a/testsuite/MDAnalysisTests/transformations/test_fit.py +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -119,13 +119,21 @@ def test_fit_translation_com(fit_universe): assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) -def test_fit_translation_plane(fit_universe): +@pytest.mark.parametrize('plane', ( + "yz", + "xz", + "xy") +) +def test_fit_translation_plane(fit_universe, plane): test_u = fit_universe[0] ref_u = fit_universe[1] - # translate the test universe on the x and y coordinates only - fit_translation(test_u, ref_u, plane="xy")(test_u.trajectory.ts) - # the reference is 10 angstrom in the z coordinate above the test universe - shiftz = np.asanyarray([0, 0, -10], np.float32) + axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} + idx = axes[plane] + # translate the test universe on the plane coordinates only + fit_translation(test_u, ref_u, plane=plane)(test_u.trajectory.ts) + # the reference is 10 angstrom in all coordinates above the test universe + shiftz = np.asanyarray([0, 0, 0], np.float32) + shiftz[idx] = -10 ref_coordinates = ref_u.trajectory.ts.positions + shiftz assert_array_almost_equal(test_u.trajectory.ts.positions, ref_coordinates, decimal=6) @@ -220,18 +228,29 @@ def test_fit_rot_trans_no_options(fit_universe): assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) -def test_fit_rot_trans_plane(fit_universe): +@pytest.mark.parametrize('plane', ( + "yz", + "xz", + "xy") +) +def test_fit_rot_trans_plane(fit_universe, plane): # the reference is rotated in the x axis so removing the translations and rotations # in the yz plane should return the same as the fitting without specifying a plane test_u = fit_universe[0] ref_u = fit_universe[1] ref_com = ref_u.atoms.center(None) + mobile_com = test_u.atoms.center(None) + axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} + idx = axes[plane] + rotaxis = np.asarray([0,0,0]) + rotaxis[idx]=1 ref_u.trajectory.ts.positions -= ref_com - R = rotation_matrix(np.pi/3, [1,0,0])[:3,:3] + R = rotation_matrix(np.pi/3, rotaxis)[:3,:3] ref_u.trajectory.ts.positions = np.dot(ref_u.trajectory.ts.positions, R) + ref_com[idx] = mobile_com[idx] ref_u.trajectory.ts.positions += ref_com - fit_rot_trans(test_u, ref_u, plane="yz")(test_u.trajectory.ts) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) + fit_rot_trans(test_u, ref_u, plane=plane)(test_u.trajectory.ts) + assert_array_almost_equal(test_u.trajectory.ts.positions[:,idx], ref_u.trajectory.ts.positions[:,idx], decimal=3) def test_fit_rot_trans_transformations_api(fit_universe): @@ -244,4 +263,4 @@ def test_fit_rot_trans_transformations_api(fit_universe): ref_u.trajectory.ts.positions += ref_com transform = fit_rot_trans(test_u, ref_u) test_u.trajectory.add_transformations(transform) -# assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file From e5fef8e56da92605cca9516ff85735b23454f3cc Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Sat, 28 Jul 2018 17:57:27 +0100 Subject: [PATCH 14/46] improved docs --- package/MDAnalysis/transformations/fit.py | 26 ++++++++++++------- .../trajectory_transformations.rst | 1 + .../transformations/fit.rst | 1 + 3 files changed, 18 insertions(+), 10 deletions(-) create mode 100644 package/doc/sphinx/source/documentation_pages/transformations/fit.rst diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py index 3fe8595fe25..6f0a326041c 100644 --- a/package/MDAnalysis/transformations/fit.py +++ b/package/MDAnalysis/transformations/fit.py @@ -24,9 +24,13 @@ Fitting transformations --- :mod:`MDAnalysis.transformations.fit` ================================================================= -Translate and rotates the coordinates of a given trajectory to align +Translate and/or rotates the coordinates of a given trajectory to align a given AtomGroup to a reference structure. - + +.. autofunction:: fit_translation + +.. autofunction:: fit_rot_trans + """ from __future__ import absolute_import @@ -50,11 +54,12 @@ def fit_translation(ag, reference, plane=None, weights=None): Removing the translations of a given AtomGroup `ag` on the XY plane by fitting its center of mass to the center of mass of a reference `ref`: - ..code-block::python + .. code-block:: python ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") - transform = mda.transformations.fit_translation(ag, ref, plane="xy", center_of="mass") + transform = mda.transformations.fit_translation(ag, ref, plane="xy", + center_of="mass") u.trajectory.add_transformations(transform) Parameters @@ -127,17 +132,18 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): Example ------- Removing the translations and rotations of a given AtomGroup `ag` on the XY plane - by fitting it to a reference `ref`, using the masses as weights for the RMSD fit. + by fitting it to a reference `ref`, using the masses as weights for the RMSD fit: - ..code-block::python + .. code-block:: python ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") - transform = mda.transformations.fit_rot_trans(ag, ref, plane="xy", weights="mass") + transform = mda.transformations.fit_rot_trans(ag, ref, plane="xy", + weights="mass") u.trajectory.add_transformations(transform) - - Parameter - --------- + + Parameters + ---------- ag : Universe or AtomGroup structure to translate, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole diff --git a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst index a5f29f4a371..051f7b726fa 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -126,3 +126,4 @@ e.g. giving a workflow as a keyword argument when defining the universe: ./transformations/translate ./transformations/rotate ./transformations/positionaveraging + ./transformations/fit diff --git a/package/doc/sphinx/source/documentation_pages/transformations/fit.rst b/package/doc/sphinx/source/documentation_pages/transformations/fit.rst new file mode 100644 index 00000000000..a5fe6293a83 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/transformations/fit.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.transformations.fit \ No newline at end of file From cdc1205cf779bf095bfef71b36efde182b0b1219 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Sun, 28 Apr 2019 16:01:12 -0700 Subject: [PATCH 15/46] MAINT: reviewer adjusts PR 1991 * test_fit.py now uses proper __future__ imports, which should alleviate code style / linting issues for old style division * test_fit.py adjusted to bump coverage of fit.py to 100 % --- testsuite/MDAnalysisTests/transformations/test_fit.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_fit.py b/testsuite/MDAnalysisTests/transformations/test_fit.py index 2f9b3008c33..84b7f71b954 100644 --- a/testsuite/MDAnalysisTests/transformations/test_fit.py +++ b/testsuite/MDAnalysisTests/transformations/test_fit.py @@ -19,7 +19,7 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -from __future__ import absolute_import +from __future__ import division, print_function, absolute_import import numpy as np import pytest @@ -110,6 +110,11 @@ def test_fit_translation_no_options(fit_universe): # what happens when no options are passed? assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=6) +def test_fit_translation_residue_mismatch(fit_universe): + test_u = fit_universe[0] + ref_u = fit_universe[1].residues[:-1].atoms + with pytest.raises(ValueError, match='number of residues'): + fit_translation(test_u, ref_u)(test_u.trajectory.ts) def test_fit_translation_com(fit_universe): test_u = fit_universe[0] @@ -263,4 +268,4 @@ def test_fit_rot_trans_transformations_api(fit_universe): ref_u.trajectory.ts.positions += ref_com transform = fit_rot_trans(test_u, ref_u) test_u.trajectory.add_transformations(transform) - assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) \ No newline at end of file + assert_array_almost_equal(test_u.trajectory.ts.positions, ref_u.trajectory.ts.positions, decimal=3) From df3684c921c538e5c4df3bb42fda3ea9780f3f44 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Tue, 11 Jun 2019 19:35:38 -0700 Subject: [PATCH 16/46] MAINT: final reviewer adjustments in PR 1991 * just minor DOC adjustments noticed in a final review --- package/MDAnalysis/transformations/fit.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/transformations/fit.py b/package/MDAnalysis/transformations/fit.py index 6f0a326041c..63ee6213544 100644 --- a/package/MDAnalysis/transformations/fit.py +++ b/package/MDAnalysis/transformations/fit.py @@ -59,7 +59,7 @@ def fit_translation(ag, reference, plane=None, weights=None): ag = u.select_atoms("protein") ref = mda.Universe("reference.pdb") transform = mda.transformations.fit_translation(ag, ref, plane="xy", - center_of="mass") + weights="mass") u.trajectory.add_transformations(transform) Parameters @@ -121,9 +121,9 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): """Perform a spatial superposition by minimizing the RMSD. Spatially align the group of atoms `ag` to `reference` by doing a RMSD - fit on `select` atoms. + fit. - This fit works as a way o remove translations and rotations of a given + This fit works as a way to remove translations and rotations of a given AtomGroup in a trajectory. A plane can be given using the flag `plane` so that only translations and rotations in that particular plane are removed. This is useful for protein-membrane systems to where the membrane @@ -145,7 +145,7 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): Parameters ---------- ag : Universe or AtomGroup - structure to translate, a + structure to translate and rotate, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole :class:`~MDAnalysis.core.universe.Universe` reference : Universe or AtomGroup @@ -153,7 +153,7 @@ def fit_rot_trans(ag, reference, plane=None, weights=None): :class:`~MDAnalysis.core.universe.Universe` plane: str, optional used to define the plane on which the rotations and translations will be removed. - Defined as a string of the plane. Suported planes are "yz", "xz" and "xy" planes. + Defined as a string of the plane. Supported planes are "yz", "xz" and "xy" planes. weights : {"mass", ``None``} or array_like, optional choose weights. With ``"mass"`` uses masses as weights; with ``None`` weigh each atom equally. If a float array of the same length as From 7eaf796bc180016b6970603d668203c4fc013ff7 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 17/46] Add unwrap to group.radius_of_gyration --- package/MDAnalysis/core/topologyattrs.py | 10 +++++++++- testsuite/MDAnalysisTests/core/test_atomgroup.py | 8 ++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 22a71302b30..c70b02a4644 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -953,15 +953,23 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) masses = atomgroup.masses + unwrap = kwargs.pop('unwrap', False) + compound = kwargs.pop('compound', 'group') + + 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() - com = atomgroup.center_of_mass(pbc=pbc) if pbc: recenteredpos = atomgroup.pack_into_box(inplace=False) - com + elif unwrap: + recenteredpos = atomgroup.unwrap(compound=compound) - com else: recenteredpos = atomgroup.positions - com diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c17c3013233..72074c51666 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -859,6 +859,7 @@ def ref_noUnwrap_residues(self): [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), 'Asph': 0.02060121, + 'ROG': 27.713009, } @pytest.fixture() @@ -874,6 +875,7 @@ def ref_Unwrap_residues(self): [-1330.489, 19315.253, 3306.212], [ 2938.243, 3306.212, 8990.481]]), 'Asph': 0.2969491080, + 'ROG': 40.686541, } @pytest.fixture() @@ -886,6 +888,7 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, + 'ROG': 0.9602093, } @pytest.fixture() @@ -898,6 +901,7 @@ def ref_Unwrap(self): [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), 'Asph': 1.0, + 'ROG': 1.1135230, } def test_default_residues(self, ag, ref_noUnwrap_residues): @@ -905,12 +909,14 @@ def test_default_residues(self, ag, ref_noUnwrap_residues): 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) + assert_almost_equal(ag.radius_of_gyration(compound='residues'), ref_noUnwrap_residues['ROG'], 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) + assert_almost_equal(ag.radius_of_gyration(unwrap=True, compound='residues'), ref_Unwrap_residues['ROG'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -922,6 +928,7 @@ def test_default(self, ref_noUnwrap): 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) + assert_almost_equal(group.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -932,6 +939,7 @@ def test_UnWrapFlag(self, ref_Unwrap): 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) + assert_almost_equal(group.radius_of_gyration(unwrap=True), ref_Unwrap['ROG'], self.prec) class TestPBCFlag(object): From 94cac246b99de97a8bff6520b888c38dfaf898e5 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:02:25 +0530 Subject: [PATCH 18/46] Adds docstring --- package/MDAnalysis/core/topologyattrs.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index c70b02a4644..2c9a5aed342 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -992,6 +992,10 @@ def shape_parameter(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ---- From f72ca1834795a4bbb9c18bd4d5ca4c064fd1b94d Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:08:43 +0530 Subject: [PATCH 19/46] Adds CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 83a6dbe1228..b25b1c3cbc6 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -61,6 +61,7 @@ Enhancements * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) * added unwrap keyword to asphericity (PR #2290) + * added unwrap keyword to radius_of_gyration (PR #2294) Changes * added official support for Python 3.7 (PR #1963) From 9fe2c6b5d05eaad0cd3d732d52888459e281872f Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 11 Jul 2019 03:01:08 +0530 Subject: [PATCH 20/46] Adds decorator --- package/MDAnalysis/core/topologyattrs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 2c9a5aed342..61a0d6a97ef 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -937,6 +937,7 @@ def moment_of_inertia(group, **kwargs): ('moment_of_inertia', moment_of_inertia)) @warn_if_not_unique + @check_pbc_and_unwrap def radius_of_gyration(group, **kwargs): """Radius of gyration. @@ -953,8 +954,7 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword - .. versionchanged:: 0.20.0 Added *unwrap* parameter - + .. versionchanged:: 0.20.0 Added *unwrap* and *compound* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) From 6a791d48f96370f133145ab323bc5172cc849ed3 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Mon, 17 Jun 2019 13:15:32 +0530 Subject: [PATCH 21/46] Adds unwrap to center --- package/MDAnalysis/core/groups.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index ef80087d794..f8a4cb7429d 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -750,6 +750,8 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': + if unwrap: + coords = atoms.unwrap(inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) elif unwrap: @@ -787,6 +789,8 @@ 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) @@ -797,6 +801,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=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: From 9c18ce45e3455f042fbe26c01a91bf72564dce86 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 18 Jun 2019 23:26:28 +0530 Subject: [PATCH 22/46] Adds tests --- package/MDAnalysis/core/groups.py | 3 +-- testsuite/MDAnalysisTests/core/test_atomgroup.py | 1 - testsuite/MDAnalysisTests/core/util.py | 11 +++++++++++ 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index f8a4cb7429d..6e5756cb2f4 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -751,7 +751,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) elif unwrap: @@ -801,7 +801,6 @@ def center(self, weights, pbc=None, compound='group', unwrap=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 c17c3013233..8ba0216b95d 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -539,7 +539,6 @@ def test_center_unwrap(self, level, compound, is_triclinic): # get the expected results center = group.center(weights=None, pbc=False, compound=compound, unwrap=True) - ref_center = u.center(compound=compound) assert_almost_equal(ref_center, center, decimal=4) diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 653c743d19c..a54d38984b9 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -601,6 +601,17 @@ def center(self, compound): for base in range(23, 47, 8): loc_center = np.mean(relpos[base:base + 8, :], axis=0) center_pos[pos,:] = loc_center + 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": From 93d5ce8cffdcbb14b55ee47924aead25571fefbe Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 10:41:37 +0530 Subject: [PATCH 23/46] Review changes --- package/MDAnalysis/core/groups.py | 5 +++-- testsuite/MDAnalysisTests/core/util.py | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 6e5756cb2f4..78b3c8f67bd 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -750,6 +750,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: @@ -789,8 +792,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/util.py b/testsuite/MDAnalysisTests/core/util.py index a54d38984b9..42d95d457c0 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -606,12 +606,12 @@ def center(self, compound): 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 918d5a08c1ff11eacdc44e6cc19bf7502bb6f4aa Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 18:26:59 +0530 Subject: [PATCH 24/46] add unwrap to 'cog' --- package/MDAnalysis/core/groups.py | 5 +++-- testsuite/MDAnalysisTests/core/test_atomgroup.py | 5 ++++- testsuite/MDAnalysisTests/core/util.py | 10 +--------- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 78b3c8f67bd..0b338132ab7 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -753,14 +753,15 @@ 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) 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 diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 8ba0216b95d..e42c63bb38b 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -838,6 +838,7 @@ class TestUnwrapFlag(object): prec = 3 @pytest.fixture() +<<<<<<< HEAD def ag(self): universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues[0:3] @@ -885,7 +886,8 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, - } + + } @pytest.fixture() def ref_Unwrap(self): @@ -922,6 +924,7 @@ def test_default(self, ref_noUnwrap): 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) group = u.atoms[31:39] # molecules 11 diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 42d95d457c0..56b4c958250 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -604,15 +604,7 @@ def center(self, compound): 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_center = np.mean(relpos[base:base + 4, :], axis=0) - center_pos[pos,:] = loc_center - pos+=1 - for base in range(23, 47, 8): - loc_center = np.mean(relpos[base:base + 8, :], axis=0) - center_pos[pos,:] = loc_center - pos+=1 + if compound == "group": center_pos = center_pos[11] From 3ebbe238338eba86bcec94de6010326eb3a17756 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 4 Jul 2019 03:04:34 +0530 Subject: [PATCH 25/46] 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 e844781e34929e6afafcee7d132f37a30045b9ef Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 26/46] Add unwrap to group.radius_of_gyration --- package/MDAnalysis/core/topologyattrs.py | 10 +++++++++- testsuite/MDAnalysisTests/core/test_atomgroup.py | 12 +++++++++--- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 22a71302b30..c70b02a4644 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -953,15 +953,23 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) masses = atomgroup.masses + unwrap = kwargs.pop('unwrap', False) + compound = kwargs.pop('compound', 'group') + + 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() - com = atomgroup.center_of_mass(pbc=pbc) if pbc: recenteredpos = atomgroup.pack_into_box(inplace=False) - com + elif unwrap: + recenteredpos = atomgroup.unwrap(compound=compound) - com else: recenteredpos = atomgroup.positions - com diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index e42c63bb38b..915da14e639 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -859,6 +859,7 @@ def ref_noUnwrap_residues(self): [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), 'Asph': 0.02060121, + 'ROG': 27.713009, } @pytest.fixture() @@ -874,6 +875,7 @@ def ref_Unwrap_residues(self): [-1330.489, 19315.253, 3306.212], [ 2938.243, 3306.212, 8990.481]]), 'Asph': 0.2969491080, + 'ROG': 40.686541, } @pytest.fixture() @@ -886,8 +888,8 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, - - } + 'ROG': 0.9602093, + } @pytest.fixture() def ref_Unwrap(self): @@ -899,6 +901,7 @@ def ref_Unwrap(self): [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), 'Asph': 1.0, + 'ROG': 1.1135230, } def test_default_residues(self, ag, ref_noUnwrap_residues): @@ -906,12 +909,14 @@ def test_default_residues(self, ag, ref_noUnwrap_residues): 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) + assert_almost_equal(ag.radius_of_gyration(compound='residues'), ref_noUnwrap_residues['ROG'], 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) + assert_almost_equal(ag.radius_of_gyration(unwrap=True, compound='residues'), ref_Unwrap_residues['ROG'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -923,7 +928,7 @@ def test_default(self, ref_noUnwrap): 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) - + assert_almost_equal(group.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -934,6 +939,7 @@ def test_UnWrapFlag(self, ref_Unwrap): 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) + assert_almost_equal(group.radius_of_gyration(unwrap=True), ref_Unwrap['ROG'], self.prec) class TestPBCFlag(object): From c37d2ac924f09254ea873995cbdfe7e7e38b1526 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:02:25 +0530 Subject: [PATCH 27/46] Adds docstring --- package/MDAnalysis/core/topologyattrs.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index c70b02a4644..2c9a5aed342 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -992,6 +992,10 @@ def shape_parameter(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ---- From c5ca308bdfdfed7cc1da27df95c40ed19028fa87 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:08:43 +0530 Subject: [PATCH 28/46] Adds CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 83a6dbe1228..b25b1c3cbc6 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -61,6 +61,7 @@ Enhancements * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) * added unwrap keyword to asphericity (PR #2290) + * added unwrap keyword to radius_of_gyration (PR #2294) Changes * added official support for Python 3.7 (PR #1963) From 2e17578719a7cce1620c7e9de20abd80d14783b8 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 11 Jul 2019 03:01:08 +0530 Subject: [PATCH 29/46] Adds decorator --- package/MDAnalysis/core/topologyattrs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 2c9a5aed342..61a0d6a97ef 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -937,6 +937,7 @@ def moment_of_inertia(group, **kwargs): ('moment_of_inertia', moment_of_inertia)) @warn_if_not_unique + @check_pbc_and_unwrap def radius_of_gyration(group, **kwargs): """Radius of gyration. @@ -953,8 +954,7 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword - .. versionchanged:: 0.20.0 Added *unwrap* parameter - + .. versionchanged:: 0.20.0 Added *unwrap* and *compound* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) From d66ab5f097340fee6f8dcee37516755362c40c07 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Mon, 17 Jun 2019 13:15:32 +0530 Subject: [PATCH 30/46] Adds unwrap to center --- package/MDAnalysis/core/groups.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index ef80087d794..123352335d8 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -690,8 +690,11 @@ 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, optional - If ``True``, compounds will be unwrapped before computing their centers. + 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 ------- @@ -750,6 +753,8 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': + if unwrap: + coords = atoms.unwrap(inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) elif unwrap: @@ -787,6 +792,8 @@ 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) From 95a7c0fcdf340b035461327a051e72c88833ac8e Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 18 Jun 2019 23:26:28 +0530 Subject: [PATCH 31/46] Adds tests --- package/MDAnalysis/core/groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 123352335d8..153e19f404d 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -754,7 +754,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) elif unwrap: From ca5268f1c3658d8bb23c6208e1f83aaf4ee07ea0 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 10:41:37 +0530 Subject: [PATCH 32/46] Review changes --- package/MDAnalysis/core/groups.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 153e19f404d..78b3c8f67bd 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -690,11 +690,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 ------- @@ -753,6 +750,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: @@ -792,8 +792,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) From d2948bbf209be2f7e519e1993c7454dd770293d5 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 18:26:59 +0530 Subject: [PATCH 33/46] add unwrap to 'cog' --- package/MDAnalysis/core/groups.py | 5 +++-- testsuite/MDAnalysisTests/core/test_atomgroup.py | 3 +++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 78b3c8f67bd..0b338132ab7 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -753,14 +753,15 @@ 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) 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 diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c17c3013233..66c3c45a3ec 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -839,6 +839,7 @@ class TestUnwrapFlag(object): prec = 3 @pytest.fixture() +<<<<<<< HEAD def ag(self): universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues[0:3] @@ -923,9 +924,11 @@ def test_default(self, ref_noUnwrap): 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) group = u.atoms[31:39] # molecules 11 + # Changing masses for center_of_mass 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) From f789e3464592b5cf4099a5325c70db9299faa268 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 4 Jul 2019 03:04:34 +0530 Subject: [PATCH 34/46] 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 3afb0a51fabc9ed7c0abba5d91d82c95be913580 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 35/46] Add unwrap to group.radius_of_gyration --- package/MDAnalysis/core/topologyattrs.py | 10 +++++++++- testsuite/MDAnalysisTests/core/test_atomgroup.py | 12 +++++++++--- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 22a71302b30..c70b02a4644 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -953,15 +953,23 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) masses = atomgroup.masses + unwrap = kwargs.pop('unwrap', False) + compound = kwargs.pop('compound', 'group') + + 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() - com = atomgroup.center_of_mass(pbc=pbc) if pbc: recenteredpos = atomgroup.pack_into_box(inplace=False) - com + elif unwrap: + recenteredpos = atomgroup.unwrap(compound=compound) - com else: recenteredpos = atomgroup.positions - com diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 66c3c45a3ec..2395ff2d19d 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -839,7 +839,6 @@ class TestUnwrapFlag(object): prec = 3 @pytest.fixture() -<<<<<<< HEAD def ag(self): universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues[0:3] @@ -860,6 +859,8 @@ def ref_noUnwrap_residues(self): [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), 'Asph': 0.02060121, + 'ROG': 27.713009, + } @pytest.fixture() @@ -875,6 +876,7 @@ def ref_Unwrap_residues(self): [-1330.489, 19315.253, 3306.212], [ 2938.243, 3306.212, 8990.481]]), 'Asph': 0.2969491080, + 'ROG': 40.686541, } @pytest.fixture() @@ -887,6 +889,7 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, + 'ROG': 0.9602093, } @pytest.fixture() @@ -899,6 +902,7 @@ def ref_Unwrap(self): [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), 'Asph': 1.0, + 'ROG': 1.1135230, } def test_default_residues(self, ag, ref_noUnwrap_residues): @@ -906,12 +910,14 @@ def test_default_residues(self, ag, ref_noUnwrap_residues): 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) + assert_almost_equal(ag.radius_of_gyration(compound='residues'), ref_noUnwrap_residues['ROG'], 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) + assert_almost_equal(ag.radius_of_gyration(unwrap=True, compound='residues'), ref_Unwrap_residues['ROG'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -923,7 +929,7 @@ def test_default(self, ref_noUnwrap): 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) - + assert_almost_equal(group.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -935,7 +941,7 @@ def test_UnWrapFlag(self, ref_Unwrap): 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) - + assert_almost_equal(group.radius_of_gyration(unwrap=True), ref_Unwrap['ROG'], self.prec) class TestPBCFlag(object): From 9d8e61812dc8b6fc039025606d1b8d395fa4370e Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 36/46] Add unwrap to group.radius_of_gyration --- testsuite/MDAnalysisTests/core/test_atomgroup.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 2395ff2d19d..dd84bc41707 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -860,7 +860,10 @@ def ref_noUnwrap_residues(self): [-721.50785456, -91.32156884, 6509.31735029]]), 'Asph': 0.02060121, 'ROG': 27.713009, +<<<<<<< HEAD +======= +>>>>>>> Add unwrap to group.radius_of_gyration } @pytest.fixture() @@ -943,6 +946,7 @@ def test_UnWrapFlag(self, ref_Unwrap): assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) assert_almost_equal(group.radius_of_gyration(unwrap=True), ref_Unwrap['ROG'], self.prec) + class TestPBCFlag(object): prec = 3 From 0d149500012b036472d098eef264f283b45bb24a Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:02:25 +0530 Subject: [PATCH 37/46] Adds docstring --- package/MDAnalysis/core/topologyattrs.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index c70b02a4644..2c9a5aed342 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -992,6 +992,10 @@ def shape_parameter(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ---- From bab1bf74c469afa0cfde1f73bf366bf8a6406ea2 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:08:43 +0530 Subject: [PATCH 38/46] Adds CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 83a6dbe1228..b25b1c3cbc6 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -61,6 +61,7 @@ Enhancements * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) * added unwrap keyword to asphericity (PR #2290) + * added unwrap keyword to radius_of_gyration (PR #2294) Changes * added official support for Python 3.7 (PR #1963) From a287e5db3df59b58880642ec1ecf2282cfd8ba98 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 11 Jul 2019 03:01:08 +0530 Subject: [PATCH 39/46] Adds decorator --- package/MDAnalysis/core/topologyattrs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 2c9a5aed342..61a0d6a97ef 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -937,6 +937,7 @@ def moment_of_inertia(group, **kwargs): ('moment_of_inertia', moment_of_inertia)) @warn_if_not_unique + @check_pbc_and_unwrap def radius_of_gyration(group, **kwargs): """Radius of gyration. @@ -953,8 +954,7 @@ def radius_of_gyration(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword - .. versionchanged:: 0.20.0 Added *unwrap* parameter - + .. versionchanged:: 0.20.0 Added *unwrap* and *compound* parameter """ atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) From 9f51e9aeb1feb04b0344c5fc11e7396f977d2989 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Mon, 17 Jun 2019 13:15:32 +0530 Subject: [PATCH 40/46] Adds unwrap to center --- package/MDAnalysis/core/groups.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index ef80087d794..f8a4cb7429d 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -750,6 +750,8 @@ def center(self, weights, pbc=None, compound='group', unwrap=False): comp = compound.lower() if comp == 'group': + if unwrap: + coords = atoms.unwrap(inplace=False) if pbc: coords = atoms.pack_into_box(inplace=False) elif unwrap: @@ -787,6 +789,8 @@ 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) @@ -797,6 +801,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=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: From 7c96dedffed5f9bdab969688227f41a1ee171a8e Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 18 Jun 2019 23:26:28 +0530 Subject: [PATCH 41/46] Adds tests --- package/MDAnalysis/core/groups.py | 3 +-- testsuite/MDAnalysisTests/core/test_atomgroup.py | 1 - testsuite/MDAnalysisTests/core/util.py | 11 +++++++++++ 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index f8a4cb7429d..6e5756cb2f4 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -751,7 +751,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) elif unwrap: @@ -801,7 +801,6 @@ def center(self, weights, pbc=None, compound='group', unwrap=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 dd84bc41707..af6cdba2a15 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -539,7 +539,6 @@ def test_center_unwrap(self, level, compound, is_triclinic): # get the expected results center = group.center(weights=None, pbc=False, compound=compound, unwrap=True) - ref_center = u.center(compound=compound) assert_almost_equal(ref_center, center, decimal=4) diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 653c743d19c..a54d38984b9 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -601,6 +601,17 @@ def center(self, compound): for base in range(23, 47, 8): loc_center = np.mean(relpos[base:base + 8, :], axis=0) center_pos[pos,:] = loc_center + 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": From 2fd92d82d510185195eeb9da8adb47fc388d0174 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 10:41:37 +0530 Subject: [PATCH 42/46] Review changes --- package/MDAnalysis/core/groups.py | 5 +++-- testsuite/MDAnalysisTests/core/util.py | 8 ++++---- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 6e5756cb2f4..78b3c8f67bd 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -750,6 +750,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: @@ -789,8 +792,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/util.py b/testsuite/MDAnalysisTests/core/util.py index a54d38984b9..42d95d457c0 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -606,12 +606,12 @@ def center(self, compound): 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 8145606bda41ce7d1bdd80b15437eab7f7fb97a5 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 19 Jun 2019 18:26:59 +0530 Subject: [PATCH 43/46] add unwrap to 'cog' --- package/MDAnalysis/core/groups.py | 5 +++-- testsuite/MDAnalysisTests/core/test_atomgroup.py | 11 +++++++---- testsuite/MDAnalysisTests/core/util.py | 10 +--------- 3 files changed, 11 insertions(+), 15 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 78b3c8f67bd..0b338132ab7 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -753,14 +753,15 @@ 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) 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 diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index af6cdba2a15..078f8c3024d 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -838,6 +838,7 @@ class TestUnwrapFlag(object): prec = 3 @pytest.fixture() +<<<<<<< HEAD def ag(self): universe = mda.Universe(TRZ_psf, TRZ) group = universe.residues[0:3] @@ -859,10 +860,6 @@ def ref_noUnwrap_residues(self): [-721.50785456, -91.32156884, 6509.31735029]]), 'Asph': 0.02060121, 'ROG': 27.713009, -<<<<<<< HEAD - -======= ->>>>>>> Add unwrap to group.radius_of_gyration } @pytest.fixture() @@ -891,8 +888,13 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, +<<<<<<< HEAD 'ROG': 0.9602093, } +======= + + } +>>>>>>> add unwrap to 'cog' @pytest.fixture() def ref_Unwrap(self): @@ -933,6 +935,7 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) assert_almost_equal(group.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) + def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) group = u.atoms[31:39] # molecules 11 diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 42d95d457c0..56b4c958250 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -604,15 +604,7 @@ def center(self, compound): 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_center = np.mean(relpos[base:base + 4, :], axis=0) - center_pos[pos,:] = loc_center - pos+=1 - for base in range(23, 47, 8): - loc_center = np.mean(relpos[base:base + 8, :], axis=0) - center_pos[pos,:] = loc_center - pos+=1 + if compound == "group": center_pos = center_pos[11] From 6200e3212d66feecaff9f5cb8d8bec7f86761411 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Thu, 4 Jul 2019 03:04:34 +0530 Subject: [PATCH 44/46] 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 0edc9a7bd945bd458dd521e98e42fcb92d72fc23 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Tue, 9 Jul 2019 15:52:05 +0530 Subject: [PATCH 45/46] Add unwrap to group.radius_of_gyration --- testsuite/MDAnalysisTests/core/test_atomgroup.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 078f8c3024d..f7f4412f903 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -888,13 +888,8 @@ def ref_noUnwrap(self): [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), 'Asph': 1.0, -<<<<<<< HEAD 'ROG': 0.9602093, } -======= - - } ->>>>>>> add unwrap to 'cog' @pytest.fixture() def ref_Unwrap(self): @@ -935,7 +930,6 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) assert_almost_equal(group.radius_of_gyration(), ref_noUnwrap['ROG'], self.prec) - def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) group = u.atoms[31:39] # molecules 11 From 68e5329700c99617e295d3d6598c0e22a8a2ca62 Mon Sep 17 00:00:00 2001 From: Ninad Bhat Date: Wed, 10 Jul 2019 21:02:25 +0530 Subject: [PATCH 46/46] Adds docstring --- package/MDAnalysis/core/topologyattrs.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 61a0d6a97ef..2ba01f3bccf 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -946,6 +946,10 @@ def radius_of_gyration(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ----