From ff13d12cffd4f5866011d492fe85e842a9ce3c7d Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Tue, 30 Dec 2025 16:17:22 +0530 Subject: [PATCH 1/6] Add cylindrical selection test using PeriodicKDTree --- .../core/test_atomselections.py | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 214b97e1cee..7cb960334ca 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -33,6 +33,7 @@ from MDAnalysis import SelectionError, SelectionWarning from MDAnalysis.core.selection import Parser from MDAnalysis.lib.distances import distance_array +from MDAnalysis.lib.pkdtree import PeriodicKDTree from MDAnalysis.tests.datafiles import ( DCD, GRO, @@ -1774,3 +1775,46 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name + +@pytest.fixture +def universe(): + #Fixture providing a small test Universe. + return mda.Universe(PSF, DCD) + +def test_cylindrical_layer_selection_periodic_kdtree(universe): + #Test cylindrical layer selection using PeriodicKDTree. + # Center point + center = universe.atoms[0].position + inner_radius = 5.0 + outer_radius = 10.0 + zmin, zmax = center[2] - 15.0, center[2] + 15.0 + + # MDAnalysis selection for cylindrical layer + sel = universe.select_atoms( + f"cylayer {inner_radius} {outer_radius} {zmin} {zmax} " + f"point {center[0]} {center[1]} {center[2]}" + ) + + # KDTree-based selection + kdtree = PeriodicKDTree(outer_radius, universe.dimensions) + kdtree.set_coords(universe.atoms.positions) + indices = kdtree.search(center, outer_radius) + + # Filter by cylindrical layer (radius and z-range) + pos = universe.atoms.positions[indices] + dxy = np.sqrt((pos[:, 0] - center[0])**2 + (pos[:, 1] - center[1])**2) + + mask = ( + (inner_radius <= dxy) & + (dxy <= outer_radius) & + (zmin <= pos[:, 2]) & + (pos[:, 2] <= zmax) + ) + + expected = indices[mask] + + # Assert KDTree selection matches MDAnalysis selection + assert_array_equal( + np.sort(expected), + np.sort(sel.indices) + ) \ No newline at end of file From f58bd52dfcd78a1d9eb21ec9a82262e7b777ec42 Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Thu, 1 Jan 2026 12:07:52 +0530 Subject: [PATCH 2/6] Format cylindrical selection test with pinned Black version --- .../core/test_atomselections.py | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 7cb960334ca..89ceacae468 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -1776,13 +1776,15 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name + @pytest.fixture def universe(): - #Fixture providing a small test Universe. + # Fixture providing a small test Universe. return mda.Universe(PSF, DCD) + def test_cylindrical_layer_selection_periodic_kdtree(universe): - #Test cylindrical layer selection using PeriodicKDTree. + # Test cylindrical layer selection using PeriodicKDTree. # Center point center = universe.atoms[0].position inner_radius = 5.0 @@ -1802,19 +1804,16 @@ def test_cylindrical_layer_selection_periodic_kdtree(universe): # Filter by cylindrical layer (radius and z-range) pos = universe.atoms.positions[indices] - dxy = np.sqrt((pos[:, 0] - center[0])**2 + (pos[:, 1] - center[1])**2) + dxy = np.sqrt((pos[:, 0] - center[0]) ** 2 + (pos[:, 1] - center[1]) ** 2) mask = ( - (inner_radius <= dxy) & - (dxy <= outer_radius) & - (zmin <= pos[:, 2]) & - (pos[:, 2] <= zmax) + (inner_radius <= dxy) + & (dxy <= outer_radius) + & (zmin <= pos[:, 2]) + & (pos[:, 2] <= zmax) ) expected = indices[mask] # Assert KDTree selection matches MDAnalysis selection - assert_array_equal( - np.sort(expected), - np.sort(sel.indices) - ) \ No newline at end of file + assert_array_equal(np.sort(expected), np.sort(sel.indices)) From d05ed3d06591fd1e89644e6e6ad3b2735d968a50 Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Fri, 2 Jan 2026 19:56:18 +0530 Subject: [PATCH 3/6] Add sanity test for cylayer selection using COG --- .../core/test_atomselections.py | 47 +++++-------------- 1 file changed, 11 insertions(+), 36 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 89ceacae468..a6f560d1c75 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -57,7 +57,7 @@ waterPSF, ) from numpy.lib import NumpyVersion -from numpy.testing import assert_equal +from numpy.testing import assert_array_equal, assert_equal from MDAnalysisTests import make_Universe @@ -1776,44 +1776,19 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name - -@pytest.fixture +@pytest.fixture(scope="module") def universe(): - # Fixture providing a small test Universe. - return mda.Universe(PSF, DCD) - + u = mda.Universe(PSF, DCD) + u.dimensions = np.array([100.0, 100.0, 100.0, 90.0, 90.0, 90.0]) + return u -def test_cylindrical_layer_selection_periodic_kdtree(universe): - # Test cylindrical layer selection using PeriodicKDTree. - # Center point - center = universe.atoms[0].position - inner_radius = 5.0 - outer_radius = 10.0 - zmin, zmax = center[2] - 15.0, center[2] + 15.0 +def test_cylayer_selection_parses_correctly(universe): + universe.dimensions = np.array([100, 100, 100, 90, 90, 90]) - # MDAnalysis selection for cylindrical layer sel = universe.select_atoms( - f"cylayer {inner_radius} {outer_radius} {zmin} {zmax} " - f"point {center[0]} {center[1]} {center[2]}" + "cylayer 5 10 15 -15 name CA" ) - # KDTree-based selection - kdtree = PeriodicKDTree(outer_radius, universe.dimensions) - kdtree.set_coords(universe.atoms.positions) - indices = kdtree.search(center, outer_radius) - - # Filter by cylindrical layer (radius and z-range) - pos = universe.atoms.positions[indices] - dxy = np.sqrt((pos[:, 0] - center[0]) ** 2 + (pos[:, 1] - center[1]) ** 2) - - mask = ( - (inner_radius <= dxy) - & (dxy <= outer_radius) - & (zmin <= pos[:, 2]) - & (pos[:, 2] <= zmax) - ) - - expected = indices[mask] - - # Assert KDTree selection matches MDAnalysis selection - assert_array_equal(np.sort(expected), np.sort(sel.indices)) + # Basic sanity checks + assert len(sel) > 0 + assert sel.n_atoms <= universe.atoms.n_atoms \ No newline at end of file From f61e50e4b816130d0346ea38ea11147a53685b53 Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Fri, 2 Jan 2026 22:16:34 +0530 Subject: [PATCH 4/6] Format cylayer selection test for Black --- .../core/test_atomselections.py | 85 +++++-------------- 1 file changed, 23 insertions(+), 62 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index a6f560d1c75..97f81ccf357 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -128,9 +128,7 @@ def test_resid_single(self, universe): def test_resid_range(self, universe): sel = universe.select_atoms("resid 100:105") assert_equal(sel.n_atoms, 89) - assert_equal( - sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"] - ) + assert_equal(sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"]) def test_selgroup(self, universe): sel = universe.select_atoms("not resid 100") @@ -157,23 +155,15 @@ def test_resnum_range(self, universe): sel = universe.select_atoms("resnum 100:105") assert_equal(sel.n_atoms, 89) assert_equal(sel.residues.resids, range(100, 106)) - assert_equal( - sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"] - ) + assert_equal(sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"]) def test_resname(self, universe): sel = universe.select_atoms("resname LEU") - assert_equal( - sel.n_atoms, 304, "Failed to find all 'resname LEU' atoms." - ) - assert_equal( - sel.n_residues, 16, "Failed to find all 'resname LEU' residues." - ) + assert_equal(sel.n_atoms, 304, "Failed to find all 'resname LEU' atoms.") + assert_equal(sel.n_residues, 16, "Failed to find all 'resname LEU' residues.") assert_equal( sorted(sel.indices), - sorted( - universe.select_atoms("segid 4AKE and resname LEU").indices - ), + sorted(universe.select_atoms("segid 4AKE and resname LEU").indices), "selected 'resname LEU' atoms are not the same as auto-generated s4AKE.LEU", ) @@ -187,9 +177,7 @@ def test_atom(self, universe): assert_equal(sel.resnames, ["GLY"]) assert_equal( sel.positions, - np.array( - [[20.38685226, -3.44224262, -5.92158318]], dtype=np.float32 - ), + np.array([[20.38685226, -3.44224262, -5.92158318]], dtype=np.float32), ) def test_atom_empty(self, universe): @@ -331,10 +319,7 @@ def test_same_resname(self, universe): assert_equal( len(sel), 331, - ( - "Found a wrong number of atoms with same resname as " - "resids 10 or 11" - ), + ("Found a wrong number of atoms with same resname as " "resids 10 or 11"), ) # fmt: off target_resids = np.array( @@ -429,9 +414,7 @@ def test_no_space_around_parentheses(self, universe): def test_concatenated_selection(self, universe): E151 = universe.select_atoms("segid 4AKE").select_atoms("resid 151") # note that this is not quite phi... HN should be C of prec. residue - phi151 = E151.atoms.select_atoms( - "name HN", "name N", "name CA", "name CB" - ) + phi151 = E151.atoms.select_atoms("name HN", "name N", "name CA", "name CB") assert_equal(len(phi151), 4) assert_equal( phi151[0].name, @@ -443,9 +426,7 @@ def test_global(self, universe): """Test the `global` modifier keyword (Issue 268)""" ag = universe.select_atoms("resname LYS and name NZ") # Lys amines within 4 angstrom of the backbone. - ag1 = universe.select_atoms( - "resname LYS and name NZ and around 4 backbone" - ) + ag1 = universe.select_atoms("resname LYS and name NZ and around 4 backbone") ag2 = ag.select_atoms("around 4 global backbone") assert_equal(ag2.indices, ag1.indices) @@ -496,12 +477,8 @@ def universe(self): def test_protein(self, universe): # must include non-standard residues sel = universe.select_atoms("protein or resname HAO or resname ORT") - assert_equal( - sel.n_atoms, universe.atoms.n_atoms, "failed to select peptide" - ) - assert_equal( - sel.n_residues, 6, "failed to select all peptide residues" - ) + assert_equal(sel.n_atoms, universe.atoms.n_atoms, "failed to select peptide") + assert_equal(sel.n_residues, 6, "failed to select all peptide residues") def test_resid_single(self, universe): sel = universe.select_atoms("resid 12") @@ -591,9 +568,7 @@ def test_same_fragment(self, universe, selstr): # This test comes here because it's a system with solvent, # and thus multiple fragments. sel = universe.select_atoms(selstr) - errmsg = ( - "Found a wrong number of atoms " "on the same fragment as id 1" - ) + errmsg = "Found a wrong number of atoms " "on the same fragment as id 1" assert_equal(len(sel), 3341, errmsg) errmsg = ( "Found a differ set of atoms when using the 'same " @@ -696,9 +671,7 @@ def test_passing_rdkit_kwargs_to_converter(self): def test_passing_max_matches_to_converter(self, u2): with pytest.warns(UserWarning, match="Your smarts-based") as wsmg: sel = u2.select_atoms("smarts C", smarts_kwargs=dict(maxMatches=2)) - sel2 = u2.select_atoms( - "smarts C", smarts_kwargs=dict(maxMatches=1000) - ) + sel2 = u2.select_atoms("smarts C", smarts_kwargs=dict(maxMatches=1000)) assert sel.n_atoms == 2 assert sel2.n_atoms == 3 @@ -909,9 +882,7 @@ class TestOrthogonalDistanceSelections(BaseDistanceSelection): def u(self): return mda.Universe(TRZ_psf, TRZ) - @pytest.mark.parametrize( - "meth, periodic", [("distmat", True), ("distmat", False)] - ) + @pytest.mark.parametrize("meth, periodic", [("distmat", True), ("distmat", False)]) def test_cyzone(self, u, meth, periodic): sel = Parser.parse("cyzone 5 4 -4 resid 2", u.atoms) sel.periodic = periodic @@ -1179,9 +1150,7 @@ def test_flip(self, prop, ag, op): # reference group, doing things forwards ref = ag[func(getattr(ag, self.plurals[prop]), 1.5)] - selstr = "prop 1.5 {op} {prop}".format( - op=self.opposites[op], prop=prop - ) + selstr = "prop 1.5 {op} {prop}".format(op=self.opposites[op], prop=prop) sel = ag.select_atoms(selstr) assert_equal(set(ref.indices), set(sel.indices)) @@ -1212,9 +1181,7 @@ class TestSelectionErrors(object): @staticmethod @pytest.fixture() def universe(): - return make_Universe( - ("names", "masses", "resids", "resnames", "resnums") - ) + return make_Universe(("names", "masses", "resids", "resnames", "resnums")) @pytest.mark.parametrize( "selstr", @@ -1317,9 +1284,7 @@ def test_string_selections(self, ref, sel, universe): ], ) def test_range_selections(self, seltype, ref, sel, universe): - self._check_sels( - ref.format(typ=seltype), sel.format(typ=seltype), universe - ) + self._check_sels(ref.format(typ=seltype), sel.format(typ=seltype), universe) class TestICodeSelection(object): @@ -1663,9 +1628,7 @@ def test_bool_sel_error(): def test_error_selection_for_strange_dtype(): with pytest.raises(ValueError, match="No base class defined for dtype"): - MDAnalysis.core.selection.gen_selection_class( - "star", "stars", dict, "atom" - ) + MDAnalysis.core.selection.gen_selection_class("star", "stars", dict, "atom") @pytest.mark.parametrize( @@ -1732,9 +1695,7 @@ def test_chirality(smi, chirality): assert u.atoms[0].chirality == "" assert u.atoms[1].chirality == chirality - assert_equal( - u.atoms[:3].chiralities, np.array(["", chirality, ""], dtype="U1") - ) + assert_equal(u.atoms[:3].chiralities, np.array(["", chirality, ""], dtype="U1")) @pytest.mark.parametrize( @@ -1776,19 +1737,19 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name + @pytest.fixture(scope="module") def universe(): u = mda.Universe(PSF, DCD) u.dimensions = np.array([100.0, 100.0, 100.0, 90.0, 90.0, 90.0]) return u + def test_cylayer_selection_parses_correctly(universe): universe.dimensions = np.array([100, 100, 100, 90, 90, 90]) - sel = universe.select_atoms( - "cylayer 5 10 15 -15 name CA" - ) + sel = universe.select_atoms("cylayer 5 10 15 -15 name CA") # Basic sanity checks assert len(sel) > 0 - assert sel.n_atoms <= universe.atoms.n_atoms \ No newline at end of file + assert sel.n_atoms <= universe.atoms.n_atoms From b7c3bb6ef963ddea23a2962a8ac3b12be063bcb6 Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Sun, 4 Jan 2026 09:45:19 +0530 Subject: [PATCH 5/6] Address review feedback for cylayer selection test --- .../MDAnalysisTests/core/test_atomselections.py | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index 97f81ccf357..c9758eadd1d 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -33,7 +33,6 @@ from MDAnalysis import SelectionError, SelectionWarning from MDAnalysis.core.selection import Parser from MDAnalysis.lib.distances import distance_array -from MDAnalysis.lib.pkdtree import PeriodicKDTree from MDAnalysis.tests.datafiles import ( DCD, GRO, @@ -57,7 +56,7 @@ waterPSF, ) from numpy.lib import NumpyVersion -from numpy.testing import assert_array_equal, assert_equal +from numpy.testing import assert_equal from MDAnalysisTests import make_Universe @@ -1737,19 +1736,11 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name - -@pytest.fixture(scope="module") -def universe(): +def test_cylayer_selection_parses_correctly(): u = mda.Universe(PSF, DCD) u.dimensions = np.array([100.0, 100.0, 100.0, 90.0, 90.0, 90.0]) - return u - - -def test_cylayer_selection_parses_correctly(universe): - universe.dimensions = np.array([100, 100, 100, 90, 90, 90]) - sel = universe.select_atoms("cylayer 5 10 15 -15 name CA") + sel = u.select_atoms("cylayer 5 10 15 -15 name CA") - # Basic sanity checks assert len(sel) > 0 - assert sel.n_atoms <= universe.atoms.n_atoms + assert sel.n_atoms <= u.atoms.n_atoms From 12d7d6ed707df629988c0cd7b21a674c2fa9d352 Mon Sep 17 00:00:00 2001 From: rishinamdev17 Date: Sat, 17 Jan 2026 14:54:02 +0530 Subject: [PATCH 6/6] Test cylayer selection against PeriodicKDTree reference --- .../core/test_atomselections.py | 105 ++++++++++++++---- 1 file changed, 83 insertions(+), 22 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_atomselections.py b/testsuite/MDAnalysisTests/core/test_atomselections.py index c9758eadd1d..c00282fdcfc 100644 --- a/testsuite/MDAnalysisTests/core/test_atomselections.py +++ b/testsuite/MDAnalysisTests/core/test_atomselections.py @@ -33,6 +33,7 @@ from MDAnalysis import SelectionError, SelectionWarning from MDAnalysis.core.selection import Parser from MDAnalysis.lib.distances import distance_array +from MDAnalysis.lib.pkdtree import PeriodicKDTree from MDAnalysis.tests.datafiles import ( DCD, GRO, @@ -57,6 +58,7 @@ ) from numpy.lib import NumpyVersion from numpy.testing import assert_equal +from numpy.testing import assert_array_equal from MDAnalysisTests import make_Universe @@ -127,7 +129,9 @@ def test_resid_single(self, universe): def test_resid_range(self, universe): sel = universe.select_atoms("resid 100:105") assert_equal(sel.n_atoms, 89) - assert_equal(sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"]) + assert_equal( + sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"] + ) def test_selgroup(self, universe): sel = universe.select_atoms("not resid 100") @@ -154,15 +158,23 @@ def test_resnum_range(self, universe): sel = universe.select_atoms("resnum 100:105") assert_equal(sel.n_atoms, 89) assert_equal(sel.residues.resids, range(100, 106)) - assert_equal(sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"]) + assert_equal( + sel.residues.resnames, ["GLY", "ILE", "ASN", "VAL", "ASP", "TYR"] + ) def test_resname(self, universe): sel = universe.select_atoms("resname LEU") - assert_equal(sel.n_atoms, 304, "Failed to find all 'resname LEU' atoms.") - assert_equal(sel.n_residues, 16, "Failed to find all 'resname LEU' residues.") + assert_equal( + sel.n_atoms, 304, "Failed to find all 'resname LEU' atoms." + ) + assert_equal( + sel.n_residues, 16, "Failed to find all 'resname LEU' residues." + ) assert_equal( sorted(sel.indices), - sorted(universe.select_atoms("segid 4AKE and resname LEU").indices), + sorted( + universe.select_atoms("segid 4AKE and resname LEU").indices + ), "selected 'resname LEU' atoms are not the same as auto-generated s4AKE.LEU", ) @@ -176,7 +188,9 @@ def test_atom(self, universe): assert_equal(sel.resnames, ["GLY"]) assert_equal( sel.positions, - np.array([[20.38685226, -3.44224262, -5.92158318]], dtype=np.float32), + np.array( + [[20.38685226, -3.44224262, -5.92158318]], dtype=np.float32 + ), ) def test_atom_empty(self, universe): @@ -318,7 +332,10 @@ def test_same_resname(self, universe): assert_equal( len(sel), 331, - ("Found a wrong number of atoms with same resname as " "resids 10 or 11"), + ( + "Found a wrong number of atoms with same resname as " + "resids 10 or 11" + ), ) # fmt: off target_resids = np.array( @@ -413,7 +430,9 @@ def test_no_space_around_parentheses(self, universe): def test_concatenated_selection(self, universe): E151 = universe.select_atoms("segid 4AKE").select_atoms("resid 151") # note that this is not quite phi... HN should be C of prec. residue - phi151 = E151.atoms.select_atoms("name HN", "name N", "name CA", "name CB") + phi151 = E151.atoms.select_atoms( + "name HN", "name N", "name CA", "name CB" + ) assert_equal(len(phi151), 4) assert_equal( phi151[0].name, @@ -425,7 +444,9 @@ def test_global(self, universe): """Test the `global` modifier keyword (Issue 268)""" ag = universe.select_atoms("resname LYS and name NZ") # Lys amines within 4 angstrom of the backbone. - ag1 = universe.select_atoms("resname LYS and name NZ and around 4 backbone") + ag1 = universe.select_atoms( + "resname LYS and name NZ and around 4 backbone" + ) ag2 = ag.select_atoms("around 4 global backbone") assert_equal(ag2.indices, ag1.indices) @@ -476,8 +497,12 @@ def universe(self): def test_protein(self, universe): # must include non-standard residues sel = universe.select_atoms("protein or resname HAO or resname ORT") - assert_equal(sel.n_atoms, universe.atoms.n_atoms, "failed to select peptide") - assert_equal(sel.n_residues, 6, "failed to select all peptide residues") + assert_equal( + sel.n_atoms, universe.atoms.n_atoms, "failed to select peptide" + ) + assert_equal( + sel.n_residues, 6, "failed to select all peptide residues" + ) def test_resid_single(self, universe): sel = universe.select_atoms("resid 12") @@ -567,7 +592,9 @@ def test_same_fragment(self, universe, selstr): # This test comes here because it's a system with solvent, # and thus multiple fragments. sel = universe.select_atoms(selstr) - errmsg = "Found a wrong number of atoms " "on the same fragment as id 1" + errmsg = ( + "Found a wrong number of atoms " "on the same fragment as id 1" + ) assert_equal(len(sel), 3341, errmsg) errmsg = ( "Found a differ set of atoms when using the 'same " @@ -670,7 +697,9 @@ def test_passing_rdkit_kwargs_to_converter(self): def test_passing_max_matches_to_converter(self, u2): with pytest.warns(UserWarning, match="Your smarts-based") as wsmg: sel = u2.select_atoms("smarts C", smarts_kwargs=dict(maxMatches=2)) - sel2 = u2.select_atoms("smarts C", smarts_kwargs=dict(maxMatches=1000)) + sel2 = u2.select_atoms( + "smarts C", smarts_kwargs=dict(maxMatches=1000) + ) assert sel.n_atoms == 2 assert sel2.n_atoms == 3 @@ -881,7 +910,9 @@ class TestOrthogonalDistanceSelections(BaseDistanceSelection): def u(self): return mda.Universe(TRZ_psf, TRZ) - @pytest.mark.parametrize("meth, periodic", [("distmat", True), ("distmat", False)]) + @pytest.mark.parametrize( + "meth, periodic", [("distmat", True), ("distmat", False)] + ) def test_cyzone(self, u, meth, periodic): sel = Parser.parse("cyzone 5 4 -4 resid 2", u.atoms) sel.periodic = periodic @@ -1149,7 +1180,9 @@ def test_flip(self, prop, ag, op): # reference group, doing things forwards ref = ag[func(getattr(ag, self.plurals[prop]), 1.5)] - selstr = "prop 1.5 {op} {prop}".format(op=self.opposites[op], prop=prop) + selstr = "prop 1.5 {op} {prop}".format( + op=self.opposites[op], prop=prop + ) sel = ag.select_atoms(selstr) assert_equal(set(ref.indices), set(sel.indices)) @@ -1180,7 +1213,9 @@ class TestSelectionErrors(object): @staticmethod @pytest.fixture() def universe(): - return make_Universe(("names", "masses", "resids", "resnames", "resnums")) + return make_Universe( + ("names", "masses", "resids", "resnames", "resnums") + ) @pytest.mark.parametrize( "selstr", @@ -1283,7 +1318,9 @@ def test_string_selections(self, ref, sel, universe): ], ) def test_range_selections(self, seltype, ref, sel, universe): - self._check_sels(ref.format(typ=seltype), sel.format(typ=seltype), universe) + self._check_sels( + ref.format(typ=seltype), sel.format(typ=seltype), universe + ) class TestICodeSelection(object): @@ -1627,7 +1664,9 @@ def test_bool_sel_error(): def test_error_selection_for_strange_dtype(): with pytest.raises(ValueError, match="No base class defined for dtype"): - MDAnalysis.core.selection.gen_selection_class("star", "stars", dict, "atom") + MDAnalysis.core.selection.gen_selection_class( + "star", "stars", dict, "atom" + ) @pytest.mark.parametrize( @@ -1694,7 +1733,9 @@ def test_chirality(smi, chirality): assert u.atoms[0].chirality == "" assert u.atoms[1].chirality == chirality - assert_equal(u.atoms[:3].chiralities, np.array(["", chirality, ""], dtype="U1")) + assert_equal( + u.atoms[:3].chiralities, np.array(["", chirality, ""], dtype="U1") + ) @pytest.mark.parametrize( @@ -1736,11 +1777,31 @@ def test_formal_charge_selection(sel, size, name): assert len(ag) == size assert ag.atoms[0].name == name -def test_cylayer_selection_parses_correctly(): + +def test_cylayer_selection_matches_periodic_kdtree(): u = mda.Universe(PSF, DCD) u.dimensions = np.array([100.0, 100.0, 100.0, 90.0, 90.0, 90.0]) sel = u.select_atoms("cylayer 5 10 15 -15 name CA") - assert len(sel) > 0 - assert sel.n_atoms <= u.atoms.n_atoms + ref = u.select_atoms("name CA") + center = ref.center_of_geometry() + + inner, outer, zmax, zmin = 5.0, 10.0, 15.0, -15.0 + search_radius = np.sqrt(outer**2 + zmax**2) + + kdtree = PeriodicKDTree(u.dimensions) + kdtree.set_coords(u.atoms.positions, cutoff=search_radius) + + idx = kdtree.search(center, search_radius) + + pos = u.atoms.positions[idx] + dxy = np.sqrt((pos[:, 0] - center[0]) ** 2 + (pos[:, 1] - center[1]) ** 2) + dz = pos[:, 2] - center[2] + + mask = (inner <= dxy) & (dxy <= outer) & (zmin <= dz) & (dz <= zmax) + + expected = np.sort(idx[mask]) + result = np.sort(sel.indices) + + assert_array_equal(result, expected)