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 (MDAnalysis#2992)

unwrapping (closes MDAnalysis#2984).
  • Loading branch information
@mnmelo
mnmelo authored Oct 18, 2020
1 parent 26880f0 commit beb5232
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,11 +15,13 @@ The rules for this file:
------------------------------------------------------------------------------
??/??/?? tylerjereddy, richardjgowers, IAlibay, hmacdope, orbeckst, cbouy,
lilyminium, daveminh, jbarnoud, yuxuanzhuang, VOD555, ianmkenney,
calcraven，xiki-tempula, mieczyslaw
calcraven，xiki-tempula, mieczyslaw, manuel.nuno.melo

* 2.0.0

Fixes
* AtomGroup.center now works correctly for compounds + unwrapping
(Issue #2984)
* Avoid using deprecated array indexing in topology attributes
(Issue #2990, PR #2991)
* ParmedParser no longer guesses elements if they are not recognised, instead
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 @@ -830,12 +830,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 @@ -991,6 +991,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 @@ -1013,12 +1019,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 @@ -1058,6 +1064,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 @@ -1281,22 +1293,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 @@ -1317,24 +1334,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 beb5232

Please sign in to comment.