Add isolayer selection method

package/CHANGELOG

@@ -13,7 +13,8 @@ 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
 
@@ -25,6 +26,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).
+  * Added isolayer selection method (Issue #3845)
 
 
 

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,

package/MDAnalysis/core/selection.py

@@ -356,6 +356,48 @@ 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)
+        inner = np.array(pairs_inner).T[1]
+        outer = np.array(pairs_outer).T[1]
+
+        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

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)
+
     @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)