Check each atom pair against its own distance tolerance

pyxtal/crystal.py

@@ -353,6 +353,9 @@ def _set_ion_wyckoffs(self, numIon, specie, cell, wyks):
 
         """
         numIon_added = 0
+        # The like-like tolerance. It is the right one for the checks that stay
+        # within one orbit of one species -- `short_distances` and `merge` --
+        # while `check_wp` looks the pair tolerance up per pair.
         tol = self.tol_matrix.get_tol(specie, specie)
         wyckoff_sites_tmp = []
 
@@ -431,11 +434,39 @@ def _set_ion_wyckoffs(self, numIon, specie, cell, wyks):
         return None
 
     def check_wp(self, wyckoff_sites_tmp, wyks, cell, new_site, tol):
+        """
+        Check a candidate Wyckoff site against the sites already placed.
+
+        Every pair is checked at its own tolerance, `tol_matrix[A][B]`, rather
+        than at the tolerance of the species being placed. The two coincide
+        only for pairs of like species: with `Tol_matrix(prototype="atomic")`
+        the pair tolerance is `f * (r_A + r_B)`, so using `f * 2 * r_new` for
+        every pair permits overlaps whenever the species being placed is the
+        smaller of the two, and rejects legal structures whenever it is the
+        larger.
+
+        Args:
+            wyckoff_sites_tmp: sites already placed for the current species
+            wyks: sites already placed for the preceding species
+            cell: 3x3 matrix of lattice vectors
+            new_site: the candidate `atom_site`
+            tol: fallback tolerance, used for a pair the `Tol_matrix` has no
+                value for (an element with no tabulated radius)
+
+        Returns:
+            True if every pair clears its own tolerance
+        """
         # Check current WP against existing WP's
         if new_site is None:
             return False
 
-        return all(new_site.check_with_ws2(ws, cell, tol) for ws in wyckoff_sites_tmp + wyks)
+        for ws in wyckoff_sites_tmp + wyks:
+            pair_tol = self.tol_matrix.get_tol(new_site.specie, ws.specie)
+            if pair_tol is None:
+                pair_tol = tol
+            if not new_site.check_with_ws2(ws, cell, pair_tol):
+                return False
+        return True
 
     def _check_consistency(self, site, numIon):
         num = 0

tests/test_crystal.py

@@ -3,12 +3,17 @@
 import os
 import unittest
 
+import numpy as np
 import pymatgen.analysis.structure_matcher as sm
+from ase.neighborlist import neighbor_list
 from pymatgen.core import Structure
 
 from pyxtal import pyxtal
+from pyxtal.crystal import random_crystal
 from pyxtal.lattice import Lattice
 from pyxtal.symmetry import Hall, Wyckoff_position
+from pyxtal.tolerance import Tol_matrix
+from pyxtal.wyckoff_site import atom_site
 
 
 def resource_filename(package_name, resource_path):
@@ -147,6 +152,83 @@ def test_from_tabular(self):
                                                 N_grids=100)
         assert(len(reps)==8)
 
+class TestDistanceTolerance(unittest.TestCase):
+    """Every pair of atoms must clear the `Tol_matrix` entry of *that pair*.
+
+    Regression test: `check_wp` used to be handed a single tolerance, the
+    like-like one of the species being placed, and applied it to every pair.
+    With `prototype="atomic"` the pair tolerance is `f * (r_A + r_B)`, so
+    `f * 2 * r_new` is too small whenever the species being placed is the
+    smaller of the two -- which lets the two overlap -- and too large whenever
+    it is the larger, which rejects legal structures.
+    """
+
+    @staticmethod
+    def worst_pair_ratio(struc, tm):
+        """`min(d / tol(pair))` over pairs of distinct atoms, images included.
+
+        Below 1.0 means at least one pair is closer than it was allowed to be.
+
+        An atom against its own periodic image is excluded: those are governed
+        by the cell, not by `check_wp`, and are not checked at all for an orbit
+        of multiplicity 1 (`short_distances` has no pair to look at), so a
+        lattice vector shorter than the like-like tolerance survives
+        generation. That is a separate gap from the one this class covers.
+        """
+        atoms = struc.to_ase()
+        numbers = atoms.numbers
+        elements = sorted({int(n) for n in numbers})
+        cutoff = max(tm.get_tol(a, b) for a in elements for b in elements)
+        first, second, dist = neighbor_list("ijd", atoms, cutoff)
+        distinct = first != second
+        first, second, dist = first[distinct], second[distinct], dist[distinct]
+        if len(dist) == 0:
+            return np.inf
+        tols = np.array([tm.get_tol(int(numbers[a]), int(numbers[b]))
+                         for a, b in zip(first, second)])
+        return float(np.min(dist / tols))
+
+    def test_check_wp_rejects_a_contact_below_the_pair_tolerance(self):
+        # Cs-O has to clear 2.04 A, O-O only 0.91 A. A Cs-O contact of 1.5 A
+        # sits between the two: legal under O's like-like tolerance, illegal
+        # under the pair's.
+        tm = Tol_matrix(prototype="atomic", factor=1.3)
+        assert tm.get_tol("O", "O") < 1.5 < tm.get_tol("Cs", "O")
+
+        wp = Wyckoff_position.from_group_and_letter(1, "1a")
+        cell = np.eye(3) * 10.0
+        placed = atom_site(wp, [0.0, 0.0, 0.0], "Cs")
+        candidate = atom_site(wp, [0.15, 0.0, 0.0], "O")
+
+        class _Stub:
+            tol_matrix = tm
+
+        accepted = random_crystal.check_wp(
+            _Stub(), [], [placed], cell, candidate, tm.get_tol("O", "O")
+        )
+        assert not accepted
+
+    def test_generated_structures_honour_the_pair_tolerance(self):
+        # Cs and O differ by 3.5x in covalent radius, and every site is a
+        # general position, so a wrong tolerance shows up as a real overlap.
+        # Both species orders are checked: the sites are placed one species at
+        # a time, so a tolerance taken from the species being placed makes the
+        # outcome depend on the order they are given in.
+        tm = Tol_matrix(prototype="atomic", factor=1.3)
+        for species, num_ions in ((["Cs", "O"], [1, 3]), (["O", "Cs"], [3, 1])):
+            sites = [["1a"] * n for n in num_ions]
+            for seed in range(10):
+                struc = pyxtal()
+                struc.from_random(3, 1, species, num_ions, sites=sites,
+                                  tm=tm, random_state=seed)
+                assert struc.valid
+                ratio = self.worst_pair_ratio(struc, tm)
+                assert ratio >= 1.0, (
+                    f"{species}, seed {seed}: closest contact is {ratio:.3f} "
+                    f"of the tolerance it was generated under"
+                )
+
+
 class TestAtomic2D(unittest.TestCase):
     def test_single_specie(self):
         struc = pyxtal()