From 13d4655e434351693344954b5ed2c9b8328454ee Mon Sep 17 00:00:00 2001 From: Marvin Poul Date: Tue, 7 Apr 2026 11:53:56 -0400 Subject: [PATCH 1/9] Add tests for repulse geometry function Co-Authored-By: Claude Sonnet 4.6 --- tests/test_geometry.py | 61 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 tests/test_geometry.py diff --git a/tests/test_geometry.py b/tests/test_geometry.py new file mode 100644 index 000000000..ebfc0391b --- /dev/null +++ b/tests/test_geometry.py @@ -0,0 +1,61 @@ +import unittest + +import numpy as np +from ase.build import bulk + +from structuretoolkit.analyse import get_neighbors +from structuretoolkit.build.geometry import repulse + + +class TestRepulse(unittest.TestCase): + def setUp(self): + self.atoms = bulk("Cu", cubic=True).repeat(5) + self.atoms.rattle(1) + + def test_noop(self): + """If no atoms are violating min_dist, atoms should be unchanged.""" + atoms = bulk("Cu", cubic=True).repeat(5) # perfect FCC, no rattle + original_positions = atoms.positions.copy() + # Cu nearest-neighbor ~2.55 Å, so min_dist=2.0 triggers no displacement + result = repulse(atoms, min_dist=2.0) + np.testing.assert_array_equal(result.positions, original_positions) + + def test_inplace(self): + """Input atoms should be copied depending on `inplace`.""" + original_positions = self.atoms.positions.copy() + + # inplace=False (default): original must be untouched, result is a copy + result = repulse(self.atoms, inplace=False) + np.testing.assert_array_equal(self.atoms.positions, original_positions) + self.assertIsNot(result, self.atoms) + + # inplace=True: result is the same object as the input + result2 = repulse(self.atoms, inplace=True) + self.assertIs(result2, self.atoms) + + def test_iterations(self): + """Should raise error if iterations exhausted.""" + # min_dist=5.0 is far larger than any achievable spacing; step_size tiny + # → convergence is impossible, so iterations will be exhausted + with self.assertRaises(RuntimeError): + repulse(self.atoms, min_dist=5.0, step_size=0.001, iterations=2) + + def test_axis(self): + """If axis given, the other axes should be exactly untouched.""" + atoms = self.atoms.copy() + original_y = atoms.positions[:, 1].copy() + original_z = atoms.positions[:, 2].copy() + # Modify inplace so we can inspect positions even if convergence fails + try: + repulse(atoms, axis=0, inplace=True) + except RuntimeError: + pass # convergence irrelevant; we only care which axes were touched + np.testing.assert_array_equal(atoms.positions[:, 1], original_y) + np.testing.assert_array_equal(atoms.positions[:, 2], original_z) + + def test_min_dist(self): + """min_dist must be respected after a call to repulse.""" + min_dist = 1.5 + result = repulse(self.atoms, min_dist=min_dist) + neigh = get_neighbors(result, num_neighbors=1) + self.assertGreaterEqual(neigh.distances[:, 0].min(), min_dist) From 3b2d46ab567a2b07b6e3521afdaa9c725e53b94f Mon Sep 17 00:00:00 2001 From: Marvin Poul Date: Tue, 7 Apr 2026 11:57:15 -0400 Subject: [PATCH 2/9] Implement utility to displace atoms if they are too close --- src/structuretoolkit/build/geometry.py | 58 ++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 src/structuretoolkit/build/geometry.py diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py new file mode 100644 index 000000000..28469d837 --- /dev/null +++ b/src/structuretoolkit/build/geometry.py @@ -0,0 +1,58 @@ +"""Utilities that operate purely geometric aspects of structures.""" + +import numpy as np + +from structuretoolkit.analyse import get_neighbors + + +def repulse( + structure, + min_dist=1.5, + step_size=0.2, + axis=None, + iterations: int = 100, + inplace=False +): + """Displace atoms to avoid minimum overlap. + + Args: + structure (:class:`ase.Atoms`): + structure to modify + min_dist (float): + Minimum distance to enforce between atoms + step_size (float): + Maximum distance to displace atoms in one step + iterations (int): + Maximum number of displacements made before giving up + """ + if not inplace: + structure = structure.copy() + if axis is None: + axis = slice(None) + for _ in range(iterations): + neigh = get_neighbors(structure, num_neighbors=1) + dd=neigh.distances[:,0] + if dd.min()>min_dist: + break + + I = dd Date: Tue, 7 Apr 2026 12:56:12 -0400 Subject: [PATCH 3/9] Add merge function to replace two atoms by one --- src/structuretoolkit/build/geometry.py | 50 +++++++++++++++++++- tests/test_geometry.py | 64 +++++++++++++++++++++++++- 2 files changed, 112 insertions(+), 2 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 28469d837..08874a035 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -53,6 +53,54 @@ def repulse( return structure +def merge(structure: "ase.Atoms", cutoff: float = 1.8, iterations: int = 10) -> "ase.Atoms": + """Merge pairs of atoms that are closer than ``cutoff`` by collapsing each + pair to their midpoint and deleting one of the two atoms. + + The operation is applied repeatedly (up to ``iterations`` times) to handle + cases where a merge creates new close contacts. + + .. note:: + The structure is modified **in place**. Pass a copy if you need the + original to remain unchanged. + + Args: + structure (:class:`ase.Atoms`): + Structure to modify. + cutoff (float): + Distance threshold in Ångström below which two atoms are + considered clashing and will be merged. Defaults to ``1.8``. + iterations (int): + Maximum number of recursive merge passes. Defaults to ``10``. + + Returns: + :class:`ase.Atoms`: The modified structure with clashing atom pairs + replaced by single atoms at their midpoints. + """ + neigh = get_neighbors(structure, 1) + clashing = np.argwhere( neigh.distances[:,0] < cutoff ).ravel() + if len(clashing) == 0: + return structure + + moving = [] + deleting = [] + + for c in clashing: + if c in deleting: + continue + + moving.append(c) + deleting.append(neigh.indices[c, 0]) + + structure.positions[moving] += neigh.vecs[moving, 0]/2 + del structure[deleting] + + if iterations > 0: + return merge(structure, cutoff=cutoff, iterations=iterations-1) + return structure + + __all__ = [ - "repulse" + "merge", + "repulse", ] diff --git a/tests/test_geometry.py b/tests/test_geometry.py index ebfc0391b..925fdc934 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -1,10 +1,11 @@ import unittest import numpy as np +from ase import Atoms from ase.build import bulk from structuretoolkit.analyse import get_neighbors -from structuretoolkit.build.geometry import repulse +from structuretoolkit.build.geometry import merge, repulse class TestRepulse(unittest.TestCase): @@ -59,3 +60,64 @@ def test_min_dist(self): result = repulse(self.atoms, min_dist=min_dist) neigh = get_neighbors(result, num_neighbors=1) self.assertGreaterEqual(neigh.distances[:, 0].min(), min_dist) + + +def _two_atom_structure(d: float) -> Atoms: + """Return a two-Cu-atom cell with atoms separated by ``d`` Å along x.""" + return Atoms("Cu2", positions=[[0, 0, 0], [d, 0, 0]], cell=[20, 20, 20], pbc=True) + + +class TestMerge(unittest.TestCase): + def test_noop(self): + """Perfect FCC Cu has no contacts below default cutoff; structure is unchanged.""" + atoms = bulk("Cu", cubic=True).repeat(3) + original_positions = atoms.positions.copy() + result = merge(atoms) + self.assertEqual(len(result), len(atoms)) + np.testing.assert_array_equal(result.positions, original_positions) + + def test_reduces_count(self): + """Two atoms within cutoff are collapsed into one.""" + atoms = _two_atom_structure(0.5) + result = merge(atoms, cutoff=1.8) + self.assertEqual(len(result), 1) + + def test_midpoint(self): + """The surviving atom must sit at the midpoint of the original pair.""" + atoms = _two_atom_structure(1.0) + result = merge(atoms, cutoff=1.8) + self.assertEqual(len(result), 1) + np.testing.assert_allclose(result.positions[0], [0.5, 0, 0], atol=1e-10) + + def test_cutoff_respected(self): + """Atoms just beyond the cutoff must not be merged.""" + atoms = _two_atom_structure(2.0) + result = merge(atoms, cutoff=1.8) + self.assertEqual(len(result), 2) + + def test_multiple_pairs(self): + """All clashing pairs in the structure are merged in one call.""" + # Two independent close pairs, well separated from each other + atoms = Atoms( + "Cu4", + positions=[[0, 0, 0], [0.5, 0, 0], [10, 0, 0], [10.5, 0, 0]], + cell=[30, 30, 30], + pbc=True, + ) + result = merge(atoms, cutoff=1.8) + self.assertEqual(len(result), 2) + + def test_iterations_zero_stops_early(self): + """With iterations=0 only one pass runs; further clashes are left unresolved.""" + # Three atoms in a row: 0–1 and 1–2 clash. First pass merges 0+1 → 0.25. + # The merged atom is now ~9.75 Å from atom 2, so one pass is enough here + # — we just verify the function returns without error. + atoms = Atoms( + "Cu3", + positions=[[0, 0, 0], [0.5, 0, 0], [10, 0, 0]], + cell=[20, 20, 20], + pbc=True, + ) + result = merge(atoms, cutoff=1.8, iterations=0) + # At least one merge happened + self.assertLess(len(result), 3) From 43853dcaf45f439ba29578eaabd42642c576ac0d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 7 Apr 2026 17:40:01 +0000 Subject: [PATCH 4/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/structuretoolkit/build/geometry.py | 30 +++++++++++--------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 28469d837..fa1ffadaa 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -6,12 +6,12 @@ def repulse( - structure, - min_dist=1.5, - step_size=0.2, - axis=None, - iterations: int = 100, - inplace=False + structure, + min_dist=1.5, + step_size=0.2, + axis=None, + iterations: int = 100, + inplace=False, ): """Displace atoms to avoid minimum overlap. @@ -31,28 +31,24 @@ def repulse( axis = slice(None) for _ in range(iterations): neigh = get_neighbors(structure, num_neighbors=1) - dd=neigh.distances[:,0] - if dd.min()>min_dist: + dd = neigh.distances[:, 0] + if dd.min() > min_dist: break - I = dd Date: Tue, 7 Apr 2026 17:49:04 +0000 Subject: [PATCH 5/9] Complete docstring for repulse(): add axis, inplace, Returns, and Raises sections Agent-Logs-Url: https://github.com/pyiron/structuretoolkit/sessions/0ba06c20-b234-4542-9c47-90adee68ccc1 Co-authored-by: samwaseda <37879103+samwaseda@users.noreply.github.com> --- src/structuretoolkit/build/geometry.py | 33 ++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index fa1ffadaa..1a10d21bb 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -13,17 +13,40 @@ def repulse( iterations: int = 100, inplace=False, ): - """Displace atoms to avoid minimum overlap. + """Iteratively displace atoms apart until all interatomic distances exceed a minimum threshold. + + For each pair of atoms closer than ``min_dist``, the atom is displaced away from its nearest + neighbour by up to ``step_size`` along the direction of the interatomic vector. The loop + repeats until all nearest-neighbour distances satisfy the minimum criterion or the iteration + limit is reached. Args: structure (:class:`ase.Atoms`): - structure to modify + Structure to modify. min_dist (float): - Minimum distance to enforce between atoms + Minimum interatomic distance (in Å) to enforce between every pair of atoms. + Defaults to 1.5. step_size (float): - Maximum distance to displace atoms in one step + Maximum displacement (in Å) applied to a single atom per iteration. + Smaller values give smoother convergence but require more iterations. + Defaults to 0.2. + axis (int or None): + Cartesian axis index (0, 1, or 2) along which displacements are restricted. + When *None* (default) displacements are applied in all three directions. iterations (int): - Maximum number of displacements made before giving up + Maximum number of displacement steps before raising a :class:`RuntimeError`. + Defaults to 100. + inplace (bool): + If *True*, the positions of ``structure`` are modified directly. + If *False* (default), a copy is made and the original is left unchanged. + + Returns: + :class:`ase.Atoms`: The structure with adjusted atomic positions. This is the + same object as ``structure`` when ``inplace=True``, or a new copy otherwise. + + Raises: + RuntimeError: If the minimum distance criterion is not satisfied within + ``iterations`` steps. """ if not inplace: structure = structure.copy() From 9c69fdfc042c160cc8708f2d8da47c24504865ce Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 7 Apr 2026 18:17:21 +0000 Subject: [PATCH 6/9] Add type hints to repulse() Agent-Logs-Url: https://github.com/pyiron/structuretoolkit/sessions/dc485354-9225-480f-9de6-7924e42d047d Co-authored-by: samwaseda <37879103+samwaseda@users.noreply.github.com> --- src/structuretoolkit/build/geometry.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 1a10d21bb..901504322 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -1,18 +1,19 @@ """Utilities that operate purely geometric aspects of structures.""" import numpy as np +from ase.atoms import Atoms from structuretoolkit.analyse import get_neighbors def repulse( - structure, - min_dist=1.5, - step_size=0.2, - axis=None, + structure: Atoms, + min_dist: float = 1.5, + step_size: float = 0.2, + axis: int | None = None, iterations: int = 100, - inplace=False, -): + inplace: bool = False, +) -> Atoms: """Iteratively displace atoms apart until all interatomic distances exceed a minimum threshold. For each pair of atoms closer than ``min_dist``, the atom is displaced away from its nearest From f0f24e0cdb394c39a9f0984ceea8384aebe8e1f7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 8 Apr 2026 19:13:38 +0000 Subject: [PATCH 7/9] Fix convergence check (>=) and handle zero-distance division by zero in repulse() Agent-Logs-Url: https://github.com/pyiron/structuretoolkit/sessions/f6dc014c-fb3c-4ea9-b3b3-90c9dc891ce6 Co-authored-by: samwaseda <37879103+samwaseda@users.noreply.github.com> --- src/structuretoolkit/build/geometry.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 901504322..478d7750a 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -56,13 +56,17 @@ def repulse( for _ in range(iterations): neigh = get_neighbors(structure, num_neighbors=1) dd = neigh.distances[:, 0] - if dd.min() > min_dist: + if dd.min() >= min_dist: break I = dd < min_dist - vv = neigh.vecs[I, 0, :] - vv /= dd[I, None] + dd_I = dd[I] + vv = neigh.vecs[I, 0, :].copy() + # Avoid division by zero for coincident atoms; fall back to x-axis direction + vv[dd_I == 0] = [1.0, 0.0, 0.0] + safe_dd = np.where(dd_I > 0, dd_I, 1.0) + vv /= safe_dd[:, None] disp = np.clip(min_dist - dd[I], 0, step_size) From 992f34343f8ae7146533d65b006f8788518e9fe8 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 8 Apr 2026 19:22:34 +0000 Subject: [PATCH 8/9] Fix coincident-atom fallback to use opposite directions via neigh.indices Agent-Logs-Url: https://github.com/pyiron/structuretoolkit/sessions/496fe42a-e9e5-44e1-9179-c59ad7037f28 Co-authored-by: samwaseda <37879103+samwaseda@users.noreply.github.com> --- src/structuretoolkit/build/geometry.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 478d7750a..436d8926a 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -63,8 +63,15 @@ def repulse( dd_I = dd[I] vv = neigh.vecs[I, 0, :].copy() - # Avoid division by zero for coincident atoms; fall back to x-axis direction - vv[dd_I == 0] = [1.0, 0.0, 0.0] + # Avoid division by zero for coincident atoms (distance == 0). + # Assign opposite fallback directions based on atom-index ordering so + # the two coincident atoms separate rather than move together. + atom_indices = np.where(I)[0] + zero_mask = dd_I == 0 + if np.any(zero_mask): + neighbor_indices = neigh.indices[atom_indices[zero_mask], 0] + sign = np.where(atom_indices[zero_mask] < neighbor_indices, 1.0, -1.0) + vv[zero_mask] = sign[:, None] * np.array([1.0, 0.0, 0.0]) safe_dd = np.where(dd_I > 0, dd_I, 1.0) vv /= safe_dd[:, None] From 923299b050197b594a06c4b428aa06137aad9a1f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 10 Apr 2026 19:49:33 +0000 Subject: [PATCH 9/9] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/structuretoolkit/build/geometry.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/structuretoolkit/build/geometry.py b/src/structuretoolkit/build/geometry.py index 82cdece31..431a1a9ae 100644 --- a/src/structuretoolkit/build/geometry.py +++ b/src/structuretoolkit/build/geometry.py @@ -86,7 +86,9 @@ def repulse( return structure -def merge(structure: "ase.Atoms", cutoff: float = 1.8, iterations: int = 10) -> "ase.Atoms": +def merge( + structure: "ase.Atoms", cutoff: float = 1.8, iterations: int = 10 +) -> "ase.Atoms": """Merge pairs of atoms that are closer than ``cutoff`` by collapsing each pair to their midpoint and deleting one of the two atoms. @@ -111,7 +113,7 @@ def merge(structure: "ase.Atoms", cutoff: float = 1.8, iterations: int = 10) -> replaced by single atoms at their midpoints. """ neigh = get_neighbors(structure, 1) - clashing = np.argwhere( neigh.distances[:,0] < cutoff ).ravel() + clashing = np.argwhere(neigh.distances[:, 0] < cutoff).ravel() if len(clashing) == 0: return structure @@ -125,11 +127,11 @@ def merge(structure: "ase.Atoms", cutoff: float = 1.8, iterations: int = 10) -> moving.append(c) deleting.append(neigh.indices[c, 0]) - structure.positions[moving] += neigh.vecs[moving, 0]/2 + structure.positions[moving] += neigh.vecs[moving, 0] / 2 del structure[deleting] if iterations > 0: - return merge(structure, cutoff=cutoff, iterations=iterations-1) + return merge(structure, cutoff=cutoff, iterations=iterations - 1) return structure