diff --git a/package/AUTHORS b/package/AUTHORS index 52a03eab250..c092d1ee731 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -198,6 +198,7 @@ Chronological list of authors - Tengyu Xie - Raymond Zhao - Haleema Khan + - Jennifer A Clark External code diff --git a/package/CHANGELOG b/package/CHANGELOG index 0cce35f9bd9..14caed72f2b 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -13,11 +13,13 @@ The rules for this file: * release numbers follow "Semantic Versioning" http://semver.org ------------------------------------------------------------------------------ -??/??/?? IAlibay, Luthaf, hmacdope, rafaelpap, jbarnoud, BFedder, aya9aladdin +??/??/?? IAlibay, Luthaf, hmacdope, rafaelpap, jbarnoud, BFedder, aya9aladdin, + jaclark5 * 2.4.0 Fixes + * Added isolayer selection method (Issue #3845) * Auxiliary; determination of representative frames: Removed undefined behaviour for cutoff values < -1 (PR # 3749) * Fixes distances.between not always returning AtomGroup (Issue #3794) @@ -25,7 +27,7 @@ Fixes * fixes guessed_attributes and read_attributes methods of the Topology class, as those methods were giving unexpected results for some topology attributes (e.g. bonds, angles) (PR #3779). - + Enhancements diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 097d4fa0b87..fc09831d850 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -3081,6 +3081,9 @@ def select_atoms(self, sel, *othersel, periodic=True, rtol=1e-05, sphlayer *inner radius* *outer radius* *selection* Similar to sphzone, but also excludes atoms that are within *inner radius* of the selection COG + isolayer *inner radius* *outer radius* *selection* + Similar to sphlayer, but will find layer around all reference + layer, creating an iso-surface. cyzone *externalRadius* *zMax* *zMin* *selection* selects all atoms within a cylindric zone centered in the center of geometry (COG) of a given selection, diff --git a/package/MDAnalysis/core/selection.py b/package/MDAnalysis/core/selection.py index 9692fc66395..d2c32d80b50 100644 --- a/package/MDAnalysis/core/selection.py +++ b/package/MDAnalysis/core/selection.py @@ -356,6 +356,46 @@ def _apply(self, group): return group[np.asarray(indices, dtype=np.int64)] +class IsoLayerSelection(Selection): + token = 'isolayer' + precedence = 1 + + def __init__(self, parser, tokens): + super().__init__(parser, tokens) + self.periodic = parser.periodic + self.inRadius = float(tokens.popleft()) + self.exRadius = float(tokens.popleft()) + self.sel = parser.parse_expression(self.precedence) + + @return_empty_on_apply + def _apply(self, group): + indices = [] + sel = self.sel.apply(group) + # All atoms in group that aren't in sel + sys = group[~np.in1d(group.indices, sel.indices)] + + if not sys or not sel: + return sys[[]] + + box = group.dimensions if self.periodic else None + pairs_outer = distances.capped_distance(sel.positions, sys.positions, + self.exRadius, box=box, + return_distances=False) + pairs_inner = distances.capped_distance(sel.positions, sys.positions, + self.inRadius, box=box, + return_distances=False) + + if pairs_outer.size > 0: + sys_ind_outer = np.sort(np.unique(pairs_outer[:,1])) + if pairs_inner.size > 0: + sys_ind_inner = np.sort(np.unique(pairs_inner[:,1])) + indices = sys_ind_outer[~np.in1d(sys_ind_outer, sys_ind_inner)] + else: + indices = sys_ind_outer + + return sys[np.asarray(indices, dtype=np.int64)] + + class SphericalZoneSelection(Selection): token = 'sphzone' precedence = 1 diff --git a/package/doc/sphinx/source/documentation_pages/selections.rst b/package/doc/sphinx/source/documentation_pages/selections.rst index 4f0233b7360..6ee7f26bc2d 100644 --- a/package/doc/sphinx/source/documentation_pages/selections.rst +++ b/package/doc/sphinx/source/documentation_pages/selections.rst @@ -245,6 +245,16 @@ sphzone *externalRadius* *selection* resid 130 or resid 80 ) )`` selects the center of geometry of protein, resid 130, resid 80 and creates a sphere of radius 6.0 around the COG. +isolayer *inner radius* *outer radius* *selection* + Similar to sphlayer, but will find layer around all referenced atoms. + For example, if the atom types for a polymer backbone were chosen, then + an isolayer parallel to the backbone will be generated. As another + example, if a set of ions were chosen as a reference to isolate the second + hydration layer, then they will all be included in the same group. + However, in the instance that a molecule is in the second hydration layer + of one ion and the first hydration layer of another, those atoms will not + be included. + cylayer *innerRadius* *externalRadius* *zMax* *zMin* *selection* selects all atoms within a cylindric layer centered in the center of geometry (COG) of a given selection, e.g. ``cylayer 5 10 10 -8 diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 90ff242c573..6915e6a95d1 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -195,6 +195,14 @@ def test_sphlayer(self, universe, selstr): sel = universe.select_atoms(selstr) assert_equal(len(sel), 66) + @pytest.mark.parametrize('selstr', [ + 'isolayer 4.0 6.0 bynum 1281', + 'isolayer 4.0 6.0 index 1280' + ]) + def test_isolayer(self, universe, selstr): + sel = universe.select_atoms(selstr) + assert_equal(len(sel), 66) + @pytest.mark.parametrize('selstr', [ 'sphzone 6.0 bynum 1281', 'sphzone 6.0 index 1280' @@ -689,6 +697,34 @@ def test_spherical_layer(self, u, periodic): assert ref == set(result.indices) + @pytest.mark.parametrize('periodic', (True, False)) + def test_isolayer(self, u, periodic): + + rmin, rmax = 2.4, 3 + if u.atoms.types[0].isdigit(): + r1 = u.select_atoms("resid 1 or resid 2 and type 7") + else: + r1 = u.select_atoms("resid 1 or resid 2 and type O") + + sel = Parser.parse("isolayer {} {} (index {} or index {})".format(rmin, rmax, *r1.ix), u.atoms) + if periodic: + sel.periodic = True + else: + sel.periodic = False + result = sel.apply(u.atoms) + + box = u.dimensions if periodic else None + cog1 = r1[0].position.reshape(1, 3) + d1 = distance_array(u.atoms.positions, cog1, box=box) + cog2 = r1[1].position.reshape(1, 3) + d2 = distance_array(u.atoms.positions, cog2, box=box) + + ref_inner = set(np.where((d1 < rmin) | (d2 < rmin))[0]) + ref_outer = set(np.where((d1 < rmax) | (d2 < rmax))[0]) + ref_outer -= ref_inner + + assert ref_outer == set(result.indices) and len(list(ref_outer)) > 0 + @pytest.mark.parametrize('periodic', (True, False)) def test_spherical_zone(self, u, periodic): sel = Parser.parse('sphzone 5.0 resid 1', u.atoms) @@ -806,6 +842,10 @@ def test_empty_sphlayer(self, u): empty = u.select_atoms('sphlayer 2.4 6.0 name NOT_A_NAME') assert len(empty) == 0 + def test_empty_isolayer(self, u): + empty = u.select_atoms('isolayer 2.4 6.0 name NOT_A_NAME') + assert len(empty) == 0 + def test_sphzone(self, u): r1 = u.select_atoms('resid 1') cog = r1.center_of_geometry().reshape(1, 3)