Skip to content

Commit

Permalink
Fixed AtomGroup.center's behavior when using compounds in conjugation…
Browse files Browse the repository at this point in the history 
… with unwrapping (PR #2992)

fixes #2984

(cherry picked from commit beb5232)
  • Loading branch information
@mnmelo @orbeckst
mnmelo authored and orbeckst committed Oct 24, 2020
1 parent 73cd1e6 commit 49cb0a8
Show file tree
Hide file tree
Showing 3 changed files with 56 additions and 25 deletions.
4 changes: 3 additions & 1 deletion package/CHANGELOG
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ The rules for this file:
------------------------------------------------------------------------------
??/??/20 richardjgowers, IAlibay, orbeckst, tylerjereddy, jbarnoud,
yuxuanzhuang, lilyminium, VOD555, p-j-smith, bieniekmateusz,
calcraven, ianmkenney, rcrehuet
calcraven, ianmkenney, rcrehuet, manuel.nuno.melo

* 1.0.1

Expand All @@ -25,6 +25,8 @@ Fixes
* Testsuite does not use any more matplotlib.use('agg') (#2191)
* The methods provided by topology attributes now appear in the
documentation (Issue #1845)
* AtomGroup.center now works correctly for compounds + unwrapping
(Issue #2984)
* In ChainReader, read_frame does not trigger change of iterating position.
(Issue #2723, PR #2815)
* empty_atomgroup.select_atoms('name *') now returns an empty
Expand Down
10 changes: 8 additions & 2 deletions package/MDAnalysis/core/groups.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -851,12 +851,18 @@ def center(self, weights, pbc=False, compound='group', unwrap=False):

# Sort positions and weights by compound index and promote to dtype if
# required:
sort_indices = np.argsort(compound_indices)

# are we already sorted? argsorting and fancy-indexing can be expensive
if np.any(np.diff(compound_indices) < 0):
sort_indices = np.argsort(compound_indices)
else:
sort_indices = slice(None)
compound_indices = compound_indices[sort_indices]

# Unwrap Atoms
if unwrap:
coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
coords = atoms.unwrap(compound=comp, reference=None,
inplace=False)[sort_indices]
else:
coords = atoms.positions[sort_indices]
if weights is None:
Expand Down
67 changes: 45 additions & 22 deletions testsuite/MDAnalysisTests/core/test_atomgroup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -993,6 +993,12 @@ def ag(self):
group.wrap(inplace=True)
return group

@pytest.fixture()
def unordered_ag(self, ag):
ndx = np.arange(len(ag))
np.random.shuffle(ndx)
return ag[ndx]

@pytest.fixture()
def ref_noUnwrap_residues(self):
return {
Expand All @@ -1015,12 +1021,12 @@ def ref_Unwrap_residues(self):
'COG': np.array([[21.356, 41.685, 40.501],
[44.577, 43.312, 79.039],
[ 2.204, 27.722, 54.023]], dtype=np.float32),
'COM': np.array([[20.815, 42.013, 39.802],
[44.918, 43.282, 79.325],
[2.045, 28.243, 54.127]], dtype=np.float32),
'MOI': np.array([[16747.486, -1330.489, 2938.243],
[-1330.489, 19315.253, 3306.212],
[ 2938.243, 3306.212, 8990.481]]),
'COM': np.array([[21.286, 41.664, 40.465],
[44.528, 43.426, 78.671],
[ 2.111, 27.871, 53.767]], dtype=np.float32),
'MOI': np.array([[16687.941, -1330.617, 2925.883],
[-1330.617, 19256.178, 3354.832],
[ 2925.883, 3354.832, 8989.946]]),
'Asph': 0.2969491080,
}

Expand Down Expand Up @@ -1060,6 +1066,12 @@ def test_UnWrapFlag_residues(self, ag, ref_Unwrap_residues):
assert_almost_equal(ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec)
assert_almost_equal(ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec)

def test_UnWrapFlag_residues_unordered(self, unordered_ag, ref_Unwrap_residues):
assert_almost_equal(unordered_ag.center_of_geometry(unwrap=True, compound='residues'), ref_Unwrap_residues['COG'], self.prec)
assert_almost_equal(unordered_ag.center_of_mass(unwrap=True, compound='residues'), ref_Unwrap_residues['COM'], self.prec)
assert_almost_equal(unordered_ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec)
assert_almost_equal(unordered_ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec)

def test_default(self, ref_noUnwrap):
u = UnWrapUniverse(is_triclinic=False)
group = u.atoms[31:39] # molecules 11
Expand Down Expand Up @@ -1283,22 +1295,27 @@ def test_center_of_mass_compounds(self, ag, name, compound):

@pytest.mark.parametrize('name, compound', (('resids', 'residues'),
('segids', 'segments')))
def test_center_of_geometry_compounds_pbc(self, ag, name, compound):
@pytest.mark.parametrize('unwrap', (True, False))
def test_center_of_geometry_compounds_pbc(self, ag, name, compound,
unwrap):
ag.dimensions = [50, 50, 50, 90, 90, 90]
ref = [a.center_of_geometry() for a in ag.groupby(name).values()]
ref = [a.center_of_geometry(unwrap=unwrap)
for a in ag.groupby(name).values()]
ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32),
ag.dimensions)
cog = ag.center_of_geometry(pbc=True, compound=compound)
ag.dimensions)
cog = ag.center_of_geometry(pbc=True, compound=compound, unwrap=unwrap)
assert_almost_equal(cog, ref, decimal=5)

@pytest.mark.parametrize('name, compound', (('resids', 'residues'),
('segids', 'segments')))
def test_center_of_mass_compounds_pbc(self, ag, name, compound):
@pytest.mark.parametrize('unwrap', (True, False))
def test_center_of_mass_compounds_pbc(self, ag, name, compound, unwrap):
ag.dimensions = [50, 50, 50, 90, 90, 90]
ref = [a.center_of_mass() for a in ag.groupby(name).values()]
ref = [a.center_of_mass(unwrap=unwrap)
for a in ag.groupby(name).values()]
ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32),
ag.dimensions)
com = ag.center_of_mass(pbc=True, compound=compound)
ag.dimensions)
com = ag.center_of_mass(pbc=True, compound=compound, unwrap=unwrap)
assert_almost_equal(com, ref, decimal=5)

@pytest.mark.parametrize('name, compound', (('molnums', 'molecules'),
Expand All @@ -1319,24 +1336,30 @@ def test_center_of_mass_compounds_special(self, ag_molfrg,

@pytest.mark.parametrize('name, compound', (('molnums', 'molecules'),
('fragindices', 'fragments')))
@pytest.mark.parametrize('unwrap', (True, False))
def test_center_of_geometry_compounds_special_pbc(self, ag_molfrg,
name, compound):
name, compound, unwrap):
ag_molfrg.dimensions = [50, 50, 50, 90, 90, 90]
ref = [a.center_of_geometry() for a in ag_molfrg.groupby(name).values()]
ref = [a.center_of_geometry(unwrap=unwrap)
for a in ag_molfrg.groupby(name).values()]
ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32),
ag_molfrg.dimensions)
cog = ag_molfrg.center_of_geometry(pbc=True, compound=compound)
ag_molfrg.dimensions)
cog = ag_molfrg.center_of_geometry(pbc=True, compound=compound,
unwrap=unwrap)
assert_almost_equal(cog, ref, decimal=5)

@pytest.mark.parametrize('name, compound', (('molnums', 'molecules'),
('fragindices', 'fragments')))
@pytest.mark.parametrize('unwrap', (True, False))
def test_center_of_mass_compounds_special_pbc(self, ag_molfrg,
name, compound):
name, compound, unwrap):
ag_molfrg.dimensions = [50, 50, 50, 90, 90, 90]
ref = [a.center_of_mass() for a in ag_molfrg.groupby(name).values()]
ref = [a.center_of_mass(unwrap=unwrap)
for a in ag_molfrg.groupby(name).values()]
ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32),
ag_molfrg.dimensions)
com = ag_molfrg.center_of_mass(pbc=True, compound=compound)
ag_molfrg.dimensions)
com = ag_molfrg.center_of_mass(pbc=True, compound=compound,
unwrap=unwrap)
assert_almost_equal(com, ref, decimal=5)

def test_center_wrong_compound(self, ag):
Expand Down

0 comments on commit 49cb0a8

Please sign in to comment.