random structure generation

bsym/configuration.py

@@ -185,7 +185,7 @@ def tolist(self) -> list[int]:
         Returns:
             (List)
         """
-        return list(self.vector)
+        return self.vector.tolist() # type: ignore[no-any-return]
 
     def pprint(self) -> None:
         print(" ".join([str(e) for e in self.tolist()]))

bsym/configuration_space.py

@@ -121,6 +121,59 @@ def unique_configurations(self,
             )  
         return self.enumerate_configurations(generator, verbose=verbose)
 
+    def random_unique_configurations(
+        self,
+        site_distribution: dict[int, int],
+        n: int,
+        sampling: str = 'degeneracy_weighted',
+        seed: int | None = None,
+    ) -> list[Configuration]:
+        """Generate n random symmetry-inequivalent configurations.
+
+        Args:
+            site_distribution: Dictionary mapping species labels to counts.
+            n: Number of unique configurations to generate.
+            sampling: Sampling method. Either 'degeneracy_weighted' (default) or
+                'uniform'. 'degeneracy_weighted' samples configurations with
+                probability proportional to their degeneracy. 'uniform' samples
+                uniformly over equivalence classes.
+            seed: Random seed for reproducibility.
+
+        Returns:
+            List of n unique Configuration objects with count attributes set.
+
+        Raises:
+            ValueError: If sampling is not 'degeneracy_weighted' or 'uniform'.
+        """
+        if sampling not in ('degeneracy_weighted', 'uniform'):
+            raise ValueError(
+                f"sampling must be 'degeneracy_weighted' or 'uniform', got '{sampling}'"
+            )
+
+        rng = np.random.default_rng(seed)
+        seen: set[bytes] = set()
+        unique_configs: list[Configuration] = []
+
+        while len(unique_configs) < n:
+            config = self._generate_random_configuration(site_distribution, rng)
+            config_hash = config.as_bytes()
+
+            if config_hash in seen:
+                continue
+
+            equivalents = config.get_byte_equivalents(self.symmetry_group)
+            degeneracy = len(equivalents)
+
+            if sampling == 'uniform':
+                if rng.random() >= 1.0 / degeneracy:
+                    continue
+
+            seen.update(equivalents)
+            config.count = degeneracy
+            unique_configs.append(config)
+
+        return unique_configs
+
     def unique_colourings(self, colours, verbose=False):
         """
         Find the symmetry inequivalent colourings for a given number of 'colours'.
@@ -221,6 +274,38 @@ def unique_configurations_by_composition(self,
             print(f"  Total unique configurations: {sum(len(configs) for configs in results.values())}")
 
         return results
+
+    def _generate_random_configuration(
+        self,
+        site_distribution: dict[int, int],
+        rng: np.random.Generator,
+    ) -> Configuration:
+        """Generate a random configuration with the given site distribution.
+
+        Args:
+            site_distribution: Dictionary mapping species labels to counts.
+            rng: Random number generator.
+
+        Returns:
+            A random Configuration with the specified distribution.
+        """
+        n_sites = sum(site_distribution.values())
+        config = np.empty(n_sites, dtype=int)
+        available_indices = np.arange(n_sites)
+
+        # Process all but the last species
+        species_list = list(site_distribution.items())
+        for species, count in species_list[:-1]:
+            selected = _select_random_indices(available_indices, count, rng)
+            config[selected] = species
+            # Remove selected indices from available
+            available_indices = np.setdiff1d(available_indices, selected)
+
+        # Last species gets remaining indices
+        last_species, _ = species_list[-1]
+        config[available_indices] = last_species
+
+        return Configuration(config)
 
 def apply_species_mapping(config, mapping_vector):
     """
@@ -273,3 +358,20 @@ def permutation_as_config_number(p):
         tot *= 10
         tot += num
     return tot
+
+def _select_random_indices(
+    available_indices: np.ndarray,
+    count: int,
+    rng: np.random.Generator,
+) -> np.ndarray:
+    """Select count random indices from available_indices.
+
+    Args:
+        available_indices: Array of indices to select from.
+        count: Number of indices to select.
+        rng: Random number generator.
+
+    Returns:
+        Array of selected indices.
+    """
+    return rng.choice(available_indices, size=count, replace=False)

bsym/interface/pymatgen.py

@@ -496,4 +496,63 @@ def unique_structure_substitutions_by_composition(
 
         results[composition_tuple] = structures
 
-    return results
+    return results
+
+def random_unique_structure_substitutions(
+    structure,
+    to_substitute,
+    site_distribution,
+    n,
+    sampling='degeneracy_weighted',
+    seed=None,
+    atol=1e-5,
+):
+    """
+    Generate n random symmetry-unique structures by substituting sites in a pymatgen structure.
+
+    Args:
+        structure (pymatgen.Structure): The parent structure.
+        to_substitute (str): Atom label for the sites to be substituted.
+        site_distribution (dict): Dictionary mapping species to counts, e.g. {'O': 8, 'F': 8}.
+        n (int): Number of unique structures to generate.
+        sampling (str): Sampling method. Either 'degeneracy_weighted' (default) or 'uniform'.
+            'degeneracy_weighted' samples configurations with probability proportional
+            to their degeneracy. 'uniform' samples uniformly over equivalence classes.
+        seed (int, optional): Random seed for reproducibility.
+        atol (float): Tolerance factor for coordinate mapping. Default=1e-5.
+
+    Returns:
+        list[Structure]: A list of n unique Structure objects. Each has a
+            `number_of_equivalent_configurations` attribute.
+    """
+    site_substitution_index = list(structure.indices_from_symbol(to_substitute))
+
+    config_space = configuration_space_from_structure(
+        structure,
+        subset=site_substitution_index,
+        atol=atol
+    )
+
+    numeric_site_distribution, numeric_site_mapping = parse_site_distribution(
+        site_distribution
+    )
+
+    configurations = config_space.random_unique_configurations(
+        site_distribution=numeric_site_distribution,
+        n=n,
+        sampling=sampling,
+        seed=seed,
+    )
+
+    unique_structures = []
+    for config in configurations:
+        species_for_sites = [numeric_site_mapping[i] for i in config.tolist()]
+        new_structure = new_structure_from_substitution(
+            structure,
+            site_substitution_index,
+            species_for_sites
+        )
+        new_structure.number_of_equivalent_configurations = config.count
+        unique_structures.append(new_structure)
+
+    return unique_structures

bsym/symmetry_group.py

@@ -14,17 +14,17 @@ class SymmetryGroup:
 
     e.g.::
 
-        SymmetryGroup( symmetry_operations=[ s1, s2, s3 ] )
+        SymmetryGroup( symmetry_operations=[s1, s2, s3])
 
     where `s1`, `s2`, and `s3` are :any:`SymmetryOperation` objects.
 
     :any:`SymmetryGroup` objects can also be created from files using the class methods::
 
-        SymmetryGroup.read_from_file( filename )
+        SymmetryGroup.read_from_file(filename)
 
     and::
 
-        SymmetryGroup.read_from_file_with_labels( filename )
+        SymmetryGroup.read_from_file_with_labels(filename)
     """
 
     class_str = 'SymmetryGroup'
@@ -74,8 +74,8 @@ def unique_index_mappings(self) -> NDArray[np.int_]:
         return self._unique_mappings
 
     def operate_on(self,
-                   configuration: Configuration,
-                   minimal_set: bool=False) -> list[Configuration]:
+        configuration: Configuration,
+        minimal_set: bool=False) -> list[Configuration]:
         """
         Returns a list of Configurations generated by applying every symmetry operation in this symmetry group.
 

bsym/version.py

@@ -1 +1 @@
-__version__ = "2.0.0"
+__version__ = "2.1.0"

docs/source/user_guide/index.rst

@@ -12,3 +12,4 @@ These guides show you how to solve common crystallographic problems using the ``
    fixed_composition
    varying_composition
    multi_level_disorder
+   random_sampling