clease - pypi Package Compare versions

Comparing version 1.0.3 to 1.0.4

clease/calculator/calc_cache.py

+import numpy as np
+from ase.calculators.calculator import Calculator, PropertyNotImplementedError
+
+
+class CleaseCacheCalculator(Calculator):
+    """This is a simple clease calculator cache object,
+    similar to the ASE singlepoint calculator, but less restrictive"""
+
+    name = "clease_cache"
+
+    def __init__(self, **results):
+        super().__init__()
+        self.results = {}
+
+        self.results.update(**results)
+
+    def __str__(self) -> str:
+        tokens = []
+        for key, val in sorted(self.results.items()):
+            if np.isscalar(val):
+                txt = f"{key}={val}"
+            else:
+                txt = f"{key}=..."
+            tokens.append(txt)
+        return "{}({})".format(self.__class__.__name__, ", ".join(tokens))
+
+    def get_property(self, name, atoms=None, allow_calculation=True):
+        if name not in self.results:
+            raise PropertyNotImplementedError("The property '{name}' is not available.")
+
+        result = self.results[name]
+        if isinstance(result, np.ndarray):
+            result = result.copy()
+        return result
+
+    def check_state(self, *args, **kwargs):
+        """Will never fail the check_state check"""
+        # pylint: disable=unused-argument
+        return []

clease.egg-info/PKG-INFO

 Metadata-Version: 2.1
 Name: clease
-Version: 1.0.3
+Version: 1.0.4
 Summary: CLuster Expansion in Atomistic Simulation Environment
@@ -5,0 +5,0 @@ Home-page: https://gitlab.com/computationalmaterials/clease/

clease.egg-info/requires.txt

@@ -16,21 +16,21 @@ ase>=3.22
 [all]
-pip
-pytest-cov
-ipython
-black>=22.1.0
-mock
-twine
+build
 pytest
-cython
-sphinx_rtd_theme
 pyclean>=2.0.0
-pytest-benchmark[histogram]>=3.4.1
-pre-commit
 sphinx
 pylint
-clease-gui
+black>=22.1.0
+ipython
+cython
+twine
+pip
+pytest-mock
 clang-format>=14.0.3
-build
+pytest-cov
+pre-commit
 tox>=3.24.0
-pytest-mock
+sphinx_rtd_theme
+clease-gui
+pytest-benchmark[histogram]>=3.4.1
+mock
 
@@ -37,0 +37,0 @@ [dev]

clease.egg-info/SOURCES.txt

@@ -37,2 +37,3 @@ LICENSE.md
 clease/calculator/__init__.py
+clease/calculator/calc_cache.py
 clease/calculator/clease.py
@@ -39,0 +40,0 @@ clease/calculator/clease_vol_dep.py

clease/_version.txt

@@ -1,1 +0,1 @@
-1.0.3
+1.0.4

clease/calculator/__init__.py

@@ -0,1 +1,2 @@
+from .calc_cache import CleaseCacheCalculator
 from .clease import Clease, MovedIgnoredAtomError
@@ -7,2 +8,3 @@ from .clease_vol_dep import CleaseVolDep
     "Clease",
+    "CleaseCacheCalculator",
     "MovedIgnoredAtomError",
@@ -9,0 +11,0 @@ "attach_calculator",

clease/cluster/cluster_fingerprint.py

@@ -57,5 +57,5 @@ from typing import Any, Iterable
             return False
-        # We disable rtol, as it is otherwise messing with the numerical
-        # accuracy we're trying to achieve.
-        return np.allclose(self.fp, other.fp, atol=self.tol, rtol=0)
+        # Functionally equivalent to np.allclose(self.fp, other.fp, atol=self.tol, rtol=0)
+        # but ~10x faster.
+        return (np.abs(self.fp - other.fp) < self.tol).all()
 
@@ -62,0 +62,0 @@ def __getitem__(self, i):

clease/cluster/cluster_generator.py

 from itertools import product
 from typing import List, Tuple, Dict, Set, Iterator, Iterable, Union, Sequence
+from math import sqrt
 import numpy as np
@@ -47,25 +48,9 @@ from ase import Atoms
 
-    def eucledian_distance_vec(
-        self, x1: Union[np.ndarray, FourVector], x2: Union[np.ndarray, FourVector]
-    ) -> np.ndarray:
-        """
-        Calculate the difference between two vectors. Either they are in FourVector representation.
-        in which case they are translated into cartesian coordiantes, or they are assumed to be
-        NumPy arrays in cartesian coordinates.
+    def _as_euc(self, x: Union[np.ndarray, FourVector]) -> np.ndarray:
+        """Helper function to translate vector to NumPy array format"""
+        if isinstance(x, FourVector):
+            return self._fv_to_cart(x)
+        # Assume it's already in cartesian coordinates
+        return x
 
-        :param x1: First vector
-        :param x2: Second vector
-        """
-
-        def as_euc(x: Union[np.ndarray, FourVector]) -> np.ndarray:
-            """Helper function to translate vector to NumPy array format"""
-            if isinstance(x, FourVector):
-                return self._fv_to_cart(x)
-            # Assume it's already in cartesian coordinates
-            return x
-
-        euc1 = as_euc(x1)
-        euc2 = as_euc(x2)
-        return euc2 - euc1
-
     def many_to_four_vector(
@@ -136,3 +121,5 @@ self, cartesian: np.ndarray, sublattices: Sequence[int]
     def eucledian_distance(
-        self, x1: Union[np.ndarray, FourVector], x2: Union[np.ndarray, FourVector]
+        self,
+        x1: Union[np.ndarray, FourVector],
+        x2: Union[np.ndarray, FourVector],
     ) -> float:
@@ -150,4 +137,6 @@ """
 
-        d = self.eucledian_distance_vec(x1, x2)
-        return np.linalg.norm(d, ord=2)
+        # Find the distance vector
+        euc1 = self._as_euc(x1)
+        euc2 = self._as_euc(x2)
+        return _1d_array_norm(euc2 - euc1)
 
@@ -164,4 +153,4 @@ @property
 
-                diag = cellT.dot(w)
-                length: float = np.linalg.norm(diag, ord=2)
+                diag: np.ndarray = cellT.dot(w)
+                length: float = _1d_array_norm(diag)
                 yield length
@@ -296,24 +285,37 @@
 
-    def _get_internal_distances(self, figure: Figure) -> List[List[float]]:
+    def _get_internal_distances(self, figure: Figure, sort: bool = False) -> np.ndarray:
         """
-        Return all the internal distances of a figure
+        Return all the internal distances of a figure.
+
+        If ``sort`` is True, then distances are sorted from
+        largest->smallest distance for each component.
+
+        Returns a 2D array with the internal distances. If the array is *unsorted*,
+            then entry ``(i, j)`` is the eucledian distance between figures ``i`` and ``j``.
         """
-        dists = []
-        for x0 in figure.components:
-            tmp_dist = []
-            for x1 in figure.components:
-                dist = self.eucledian_distance(x0, x1)
-                # float() to ensure we're not working on a NumPy floating number
-                dist = round(float(dist), 6)
-                tmp_dist.append(dist)
-            dists.append(sorted(tmp_dist, reverse=True))
-        return dists
+        comp = figure.components
 
+        # Convert each FourVector component to cartesian coordinates
+        x = np.array([self._fv_to_cart(c) for c in comp], dtype=float)
+        # Build the distance array, such that elements (i, j) is the
+        # Euclidean distance betweens elements i and j.
+        arr = np.linalg.norm(x[:, np.newaxis] - x, axis=-1, ord=2)
+        # Round to 6 digits for numerical stability.
+        arr = np.round(arr, 6)
+        if sort:
+            # fliplr for reversing the sort. Sorted such that each row has largest->smallest
+            # distance (i.e. like reverse=True does in Python's sort() method).
+            arr = np.fliplr(np.sort(arr, axis=1))
+        return arr
+
     def _order_by_internal_distances(self, figure: Figure) -> Figure:
         """Order the Figure by internal distances. Returns a new instance of the Figure."""
-        dists = self._get_internal_distances(figure)
+        dists = self._get_internal_distances(figure, sort=True)
         fvs = figure.components
-        zipped = sorted(list(zip(dists, fvs)), reverse=True)
-        return Figure((x[1] for x in zipped))
 
+        # Sort according to the figure with the largest internal distances, lexographically.
+        # Largest internal distances first.
+        order = np.lexsort(dists.T)[::-1]
+        return Figure((fvs[i] for i in order))
+
     def to_atom_index(self, cluster: Cluster, lut: Dict[FourVector, int]) -> List[List[int]]:
@@ -334,3 +336,3 @@ """
         """Find the equivalent sites of a figure."""
-        dists = self._get_internal_distances(figure)
+        dists = self._get_internal_distances(figure, sort=True)
         equiv_sites = []
@@ -357,5 +359,4 @@ for i, d1 in enumerate(dists):
         """Return the maximum distance of a figure."""
-        internal_dists = self._get_internal_distances(figure)
-        max_dists = [x[0] for x in internal_dists]
-        return max(max_dists)
+        internal_dists = self._get_internal_distances(figure, sort=False)
+        return internal_dists.max()
 
@@ -455,1 +456,12 @@
     return ClusterFingerprint(fp=inner)
+
+
+def _1d_array_norm(arr: np.ndarray) -> float:
+    """Calculate the norm of a 1-d vector.
+    Assumes the vector is 1-dimensional and real, and will not be checked.
+
+    More efficient norm than ``np.linalg.norm`` for small vectors,
+    since we already know the dimensionality, all numbers are real,
+    and that the result is scalar.
+    """
+    return sqrt(arr.dot(arr))

clease/cluster/cluster_manager.py

@@ -1,2 +0,2 @@
-from typing import Sequence, Set, Dict, List, Iterator, Tuple
+from typing import Sequence, Set, Dict, List, Iterator, Tuple, Callable
 import logging
@@ -66,12 +66,15 @@ from itertools import product
         # If no background atoms are present, this will do nothing.
-        delete = [atom.index for atom in atoms if self.is_background_atom(atom)]
-        delete.sort(reverse=True)
-        for i in delete:
-            del atoms[i]
+        mask = self.is_background_atoms(atoms)
+        indices = np.arange(len(atoms))[mask]
+        del atoms[indices]
         return atoms
 
-    def is_background_atom(self, atom: ase.Atom) -> bool:
-        """Check whether an atom is a background atom."""
-        return atom.symbol in self.background_syms
+    def is_background_atoms(self, atoms: ase.Atoms) -> np.ndarray:
+        """Find every single background atom in an entire atoms object.
 
+        Returns a NumPy array of booleans indicating if a site is a background site."""
+        # Convert background syms to a list, as otherwise numpy treats the
+        # set as an entire element.
+        return np.isin(atoms.symbols, list(self.background_syms))
+
     def __eq__(self, other):
@@ -248,2 +251,17 @@ return self.clusters == other.clusters and self.generator == other.generator
 
+    def make_all_four_vectors(self, atoms: ase.Atoms) -> Tuple[np.ndarray, List[FourVector]]:
+        """Construct the FourVector to every site which is not a background site.
+
+        Returns an array of indices and a list with the corresponding FourVectors.
+        Indices for background sites are skipped.
+        """
+        # Find every non-background site
+        mask = ~self.is_background_atoms(atoms)
+        # Array corresponding to the atomic index of sites we calculate the FVs for.
+        indices = np.arange(len(atoms))[mask]
+        tags = atoms.get_tags()[mask]
+        pos = atoms.positions[mask]
+        four_vectors = self.generator.many_to_four_vector(pos, tags)
+        return indices, four_vectors
+
     def create_four_vector_lut(self, template: ase.Atoms) -> Dict[FourVector, int]:
@@ -259,9 +277,4 @@ """
         """
-        lut = {}
-        for atom in template:
-            if self.is_background_atom(atom):
-                # No need to make a lookup for a background atom
-                continue
-            vec = self.generator.to_four_vector(atom.position, atom.tag)
-            lut[vec] = atom.index
+        indices, four_vectors = self.make_all_four_vectors(template)
+        lut = {fv: idx for idx, fv in zip(indices, four_vectors)}
         return lut
@@ -299,7 +312,43 @@
         """
-        cell = template.get_cell()
-
         # Get the unique four-vectors which are present in all of our clusters.
         unique = self.unique_four_vectors()
 
+        wrap_fnc = self._prepare_wrap_fnc(template, unique)
+
+        lut = self.create_four_vector_lut(template)
+
+        # Map the un-translated unique 4-vectors to their index
+        unique_indx_lut = self.fourvec_to_indx(template, unique)
+        # call int() to convert from NumPy integer to python integer
+        unique_index = [int(unique_indx_lut[u]) for u in unique]
+
+        # The inverse of the LUT gives us the index of the template to its FourVector
+        # Background atoms will be missing in this table.
+        idx_to_fv = {idx: fv for fv, idx in lut.items()}
+
+        def _make_site_mapping(index: int) -> Dict[int, int]:
+            """Helper function to calculate the translation mapping for each
+            atomic site."""
+            try:
+                vec = idx_to_fv[index]
+            except KeyError:
+                # Background site
+                return {}
+
+            # Calculate the four vectors wrapped back into the supercell
+            four_vecs = wrap_fnc(vec)
+
+            # Map the translated four-vectors to the unique index
+            return {uidx: lut[fv] for uidx, fv in zip(unique_index, four_vecs)}
+
+        # Calculate the mapping for each site in the template.
+        trans_mat = [_make_site_mapping(i) for i in range(len(template))]
+        return TransMatrix(trans_mat)
+
+    def _prepare_wrap_fnc(
+        self, template: ase.Atoms, unique: Sequence[FourVector]
+    ) -> Callable[[FourVector], Iterator[FourVector]]:
+        """Prepare the wrapping function which will map a given translation FourVector
+        into the unique FourVectors with that site as a given reference.
+        """
         # Only set "trivial_supercell" to True, if we allow that path.
@@ -314,10 +363,19 @@ # Setting self._allow_trivial_path to False is only for testing purposes
         if trivial_supercell:
-            nx, ny, nz = tools.get_repetition(self.prim, template)
+            rep_arr = tools.get_repetition(self.prim, template)
+            # Pre-calculate the XYZ arrays of the unique FourVectors,
+            # as well as the sublattices. We utilize the NumPy vectorization
+            # for calculating the shift + modulo operation.
+            unique_xyz = np.array([u.xyz_array for u in unique], dtype=int)
+            sublattices = [u.sublattice for u in unique]
             wrap_fnc = functools.partial(
-                self._wrap_four_vectors_trivial, unique=unique, nx=nx, ny=ny, nz=nz
+                self._wrap_four_vectors_trivial,
+                unique_xyz=unique_xyz,
+                sublattices=sublattices,
+                rep_arr=rep_arr,
             )
-            logger.debug("Trivial supercell with repetition: (%d, %d, %d)", nx, ny, nz)
+            logger.debug("Trivial supercell with repetition: (%d, %d, %d)", *rep_arr)
         else:
             # Choose the generalized pathway.
             # Pre-calculate the inverse of the cell for faster wrapping of the positions.
+            cell = template.get_cell()
             cell_T_inv = np.linalg.inv(cell.T)
@@ -328,44 +386,18 @@ wrap_fnc = functools.partial(
             logger.debug("Non-trivial supercell, will wrap using cartesian coordinates")
+        return wrap_fnc
 
-        lut = self.create_four_vector_lut(template)
-
-        # Map the un-translated unique 4-vectors to their index
-        unique_indx_lut = self.fourvec_to_indx(template, unique)
-        # call int() to convert from NumPy integer to python integer
-        unique_index = [int(unique_indx_lut[u]) for u in unique]
-
-        def _make_site_mapping(atom: ase.Atom) -> Dict[int, int]:
-            """Helper function to calculate the translation mapping for each
-            atomic site."""
-            if self.is_background_atom(atom):
-                # This atom is considered a background, it has no mapping
-                return {}
-
-            # Translate the atom into its four-vector representation
-            vec = self.generator.to_four_vector(atom.position, sublattice=atom.tag)
-
-            # Calculate the four vectors wrapped back into the supercell
-            four_vecs = wrap_fnc(vec)
-
-            # Get the index of the translated four-vector
-            indices = [lut[fv] for fv in four_vecs]
-
-            return dict(zip(unique_index, indices))
-
-        # Calculate the mapping for each site in the template.
-        trans_mat = list(map(_make_site_mapping, template))
-        return TransMatrix(trans_mat)
-
     def _wrap_four_vectors_trivial(
         self,
         translation_vector: FourVector,
-        unique: Sequence[FourVector],
-        nx: int,
-        ny: int,
-        nz: int,
+        unique_xyz: np.ndarray,
+        sublattices: List[int],
+        rep_arr: np.ndarray,
     ) -> Iterator[FourVector]:
         """Wrap FourVectors using the trivial shift+modulo operation"""
         # pylint: disable=no-self-use
-        # Create as a generator, no need to assign this into a new list.
-        return (u.shift_xyz_and_modulo(translation_vector, nx, ny, nz) for u in unique)
+        # We use the .tolist() method, faster to iterate the Python list than the NumPy array
+        # for building the subsequent FourVectors.
+        translated = np.mod(unique_xyz + translation_vector.xyz_array, rep_arr).tolist()
+        for xyz, s in zip(translated, sublattices):
+            yield FourVector(*xyz, s)
 
@@ -372,0 +404,0 @@ def _wrap_four_vectors_general(

clease/cluster/cluster.py

@@ -131,4 +131,5 @@ from typing import Sequence, Dict, Any, List
         target_figure = self._order_equiv_sites(target_figure)
+        ref_row = trans_matrix[ref_indx]
         for figure in self.indices:
-            translated_figure = [trans_matrix[ref_indx][x] for x in figure]
+            translated_figure = [ref_row[x] for x in figure]
             translated_figure = self._order_equiv_sites(translated_figure)
@@ -143,3 +144,3 @@ if translated_figure == target_figure:
 
-    def get_all_figure_keys(self):
+    def get_all_figure_keys(self) -> List[str]:
         return [self.get_figure_key(fig) for fig in self.indices]
@@ -146,0 +147,0 @@

clease/datastructures/four_vector.py

@@ -92,27 +92,2 @@ from __future__ import annotations
 
-    def shift_xyz_and_modulo(self, other: FourVector, nx: int, ny: int, nz: int) -> FourVector:
-        """
-        Shift the four vector, similarly to "shift_xyz". Additionally, apply a modulo
-        of (nx, ny, nz) to (ix, iy, iz).
-        Useful for wrapping a four-vector back into a supercell, which is a (nx, ny, nz)
-        repetition of the primitive. The sublattice remains the same.
-
-        Return a new FourVector instance.
-
-        Example:
-
-        >>> from clease.datastructures.four_vector import FourVector
-        >>> a = FourVector(-1, -1, 3, 0)
-        >>> b = FourVector(0, 0, 0, 0)
-        >>> a.shift_xyz_and_modulo(b, 2, 2, 1)
-        FourVector(ix=1, iy=1, iz=0, sublattice=0)
-        """
-        if not isinstance(other, FourVector):
-            raise NotImplementedError(f"Shift must be by another FourVector, got {other!r}")
-
-        ix = (self.ix + other.ix) % nx
-        iy = (self.iy + other.iy) % ny
-        iz = (self.iz + other.iz) % nz
-        return FourVector(ix, iy, iz, self.sublattice)
-
     def _validate(self) -> None:
@@ -119,0 +94,0 @@ """Method for explicitly checking the validity of the FourVector.

clease/montecarlo/__init__.py

@@ -9,2 +9,3 @@ # pylint: disable=undefined-variable
 from . import base
+from . import observers
 from .bias_potential import BiasPotential
@@ -28,2 +29,3 @@ from .montecarlo import Montecarlo
     "swap_move_index_tracker",
+    "observers",
 )
@@ -30,0 +32,0 @@ __all__ = ADDITIONAL + (

clease/montecarlo/montecarlo.py

@@ -150,2 +150,9 @@ """Monte Carlo method for ase."""
 
+    def iter_observers(self) -> Iterator[MCObserver]:
+        """Directly iterate the attached observers without also getting
+        information about the interval."""
+        # Remember that the observer list contains tuples of (interval, observer)
+        for _, obs in self.observers:
+            yield obs
+
     def initialize_run(self):
@@ -152,0 +159,0 @@ """Prepare MC object for a new run."""

clease/montecarlo/observers/lowest_energy_observer.py

 import logging
 import numpy as np
-from ase.calculators.singlepoint import SinglePointCalculator
+import ase
+from clease.datastructures import MCStep
+from clease.calculator import CleaseCacheCalculator
 from .mc_observer import MCObserver
@@ -15,2 +17,7 @@
 
+    track_cf: bool
+        Whether to keep a copy of the correlation functions for the
+        emin structure. If enabled, this will be stored in the ``lowest_energy_cf``
+        variable. Defaults to False.
+
     verbose: bool
@@ -22,16 +29,37 @@ If `True`, progress messages will be printed
 
-    def __init__(self, atoms, verbose=False):
+    def __init__(self, atoms: ase.Atoms, track_cf: bool = False, verbose: bool = False):
         super().__init__()
         self.atoms = atoms
-        self.lowest_energy = np.inf
+        self.track_cf = track_cf
+        self.verbose = verbose
+        # Silence pylint
         self.lowest_energy_cf = None
-        self.emin_atoms = None
-        self.verbose = verbose
+        self.reset()
 
-    def reset(self):
+    def reset(self) -> None:
+        self.lowest_energy_cf = None
+        self.emin_atoms: ase.Atoms = self.atoms.copy()
+        # Keep a reference directly to the calculator cache
+        self._calc_cache = CleaseCacheCalculator()
+        self.emin_atoms.calc = self._calc_cache
         self.lowest_energy = np.inf
-        self.lowest_energy_cf = None
-        self.emin_atoms = None
 
     @property
+    def lowest_energy(self) -> float:
+        return self.emin_results["energy"]
+
+    @lowest_energy.setter
+    def lowest_energy(self, value: float) -> None:
+        """The currently best observed energy"""
+        self.emin_results["energy"] = value
+        # Keep a copy without needing to look up the calculator
+        # all the time, for quicker lookup.
+        self._lowest_energy_cache = value
+
+    @property
+    def emin_results(self) -> dict:
+        """The results dictionary of the lowest energy atoms"""
+        return self._calc_cache.results
+
+    @property
     def calc(self):
@@ -42,5 +70,6 @@ return self.atoms.calc
     def energy(self):
+        """The energy of the current atoms object (not the emin energy)"""
         return self.calc.results["energy"]
 
-    def __call__(self, system_changes):
+    def observe_step(self, mc_step: MCStep) -> None:
         """
@@ -50,23 +79,23 @@ Check if the current state has lower energy and store the current
 
-        system_changes: list
-            System changes. See doc-string of
-            `clease.montecarlo.observers.MCObserver`
+        mc_step: MCStep
+             Instance of MCStep with information on the latest step.
         """
-
-        if self.emin_atoms is None or self.energy < self.lowest_energy:
-            self._update_emin()
-            dE = self.energy - self.lowest_energy
-            msg = "Found new low energy structure. "
-            msg += f"New energy: {self.lowest_energy} eV. Change: {dE} eV"
-            logger.info(msg)
+        # We do the comparison even if the move was rejected,
+        # as no energy has been recorded for the first step, yet.
+        # The internal energy cache is np.inf when nothing has been recorded.
+        if mc_step.energy < self._lowest_energy_cache:
+            self._update_emin(mc_step)
+            dE = mc_step.energy - self.lowest_energy
+            msg = "Found new low energy structure. New energy: %.3f eV. Change: %.3f eV"
+            logger.info(msg, self.lowest_energy, dE)
             if self.verbose:
-                print(msg)
+                print(msg % (self.lowest_energy, dE))
 
-    def _update_emin(self):
-        self.lowest_energy_cf = self.calc.get_cf()
-        self.lowest_energy = self.energy
+    def _update_emin(self, mc_step: MCStep) -> None:
+        if self.track_cf:
+            self.lowest_energy_cf = self.calc.get_cf()
+        self.lowest_energy = mc_step.energy
 
-        # Store emin atoms, and attach a cache calculator
-        self.emin_atoms = self.atoms.copy()
-        calc_cache = SinglePointCalculator(self.emin_atoms, energy=self.energy)
-        self.emin_atoms.calc = calc_cache
+        # Avoid copying the atoms object multiple times. We can only ever have
+        # changed the symbols, so copy those over from the atoms object.
+        self.emin_atoms.numbers = self.atoms.numbers

clease/montecarlo/observers/move_obs.py

 from __future__ import annotations
 from typing import Iterator, List
 import ase
-from ase.calculators.singlepoint import SinglePointCalculator
 from clease.datastructures import MCStep
+from clease.calculator import CleaseCacheCalculator
 from .mc_observer import MCObserver
@@ -67,4 +67,4 @@
     # Attach a calculator with the energies
-    calc = SinglePointCalculator(atoms, energy=step.energy)
+    calc = CleaseCacheCalculator(energy=step.energy)
     atoms.calc = calc
     return atoms

clease/montecarlo/sgc_montecarlo.py

@@ -39,7 +39,8 @@ from typing import Sequence, Dict
             generator = RandomFlip(symbols, atoms)
+
+        self.averager = SGCObserver(atoms.calc)
         super().__init__(atoms, temp, generator=generator)
+
         self.symbols = symbols
 
-        self.averager = SGCObserver(self.atoms.calc)
-
         self.chem_pots = []
@@ -54,3 +55,3 @@ self.chem_pot_names = []
         has_attached_obs = False
-        for obs in self.observers:
+        for obs in self.iter_observers():
             if obs.name == "SGCObserver":
@@ -188,2 +189,8 @@ has_attached_obs = True
 
+    def reset_averagers(self) -> None:
+        """Reset the energy averagers, including the internal SGC Observer"""
+        super().reset_averagers()
+        # Also remember to reset the internal SGC averager
+        self.averager.reset()
+
     def get_thermodynamic_quantities(self, reset_eci=True):
@@ -190,0 +197,0 @@ """Compute thermodynamic quantities.

clease/regression/sequential_cluster_ridge.py

@@ -29,5 +29,8 @@ import logging
         Number of alpha values
+
+    verbose: bool
+        Print information about fit after completion
     """
 
-    def __init__(self, min_alpha=1e-10, max_alpha=10.0, num_alpha=20):
+    def __init__(self, min_alpha=1e-10, max_alpha=10.0, num_alpha=20, verbose: bool = False):
         super().__init__()
@@ -37,2 +40,3 @@ self.min_alpha = min_alpha
         self.num_alpha = num_alpha
+        self.verbose = verbose
 
@@ -97,3 +101,4 @@ @staticmethod
         res[: len(best_coeff)] = best_coeff
-        self._print_summary(cvs, coeffs)
+        if self.verbose:
+            self._print_summary(cvs, coeffs)
         return res

clease/tools.py

@@ -504,3 +504,3 @@ # pylint: disable=too-many-lines
 def list2str(array):
-    return "-".join(str(x) for x in array)
+    return "-".join([str(x) for x in array])
 
@@ -507,0 +507,0 @@

cxx/include/basis_function.hpp

@@ -28,2 +28,12 @@ #ifndef BASIS_FUNCTION_H
 
+    /* For preparing a vector of basis functions that maps a decoration number
+    to the old and new basis functions simultaneously.
+    This vector contains copies of the underlying basis functions.
+    The pair<double, double> is constructed as follows:
+    first -> Basis function of the *new* symbol
+    second -> Basis function of the *old* symbol
+    */
+    std::vector<std::pair<double, double>> prepare_bfs_new_old(unsigned int new_id,
+                                                               unsigned int old_id) const;
+
     /** Stream operator */
@@ -30,0 +40,0 @@ friend std::ostream &operator<<(std::ostream &out, const BasisFunction &bf);

cxx/include/ce_updater.hpp

@@ -50,2 +50,12 @@ #ifndef CE_UPDATER_H
 };
+/* Cache structure which contains objects that are useful for the updating of the spin product */
+struct SpinProductCache {
+    int ref_indx;        // The current reference index
+    int original_index;  // The original "untranslated" index.
+    /* A cache of the relevant translation matrix row
+    corresponding to the ref_indx */
+    int *trans_matrix_row;
+    // The basis functions for the new and old site symbols
+    std::vector<std::pair<double, double>> bfs_new_old;
+};
 
@@ -85,8 +95,14 @@ class CEUpdater {
     double spin_product_one_atom(int ref_indx, const Cluster &indx_list,
-                                 const std::vector<int> &dec, int ref_id);
+                                 const std::vector<int> &dec, int ref_id) const;
     // Calculate the change in spin product going from old_symb_id to new_symb_id
-    double spin_product_one_atom_delta(int ref_indx, const Cluster &indx_list,
-                                       const std::vector<int> &dec, int old_symb_id,
-                                       int new_symb_id);
+    double spin_product_one_atom_delta(const SpinProductCache &sp_cache, const Cluster &indx_list,
+                                       const equiv_deco_t &equiv_deco) const;
 
+    SpinProductCache build_sp_cache(const SymbolChange &symb_change, unsigned int old_symb_id,
+                                    unsigned int new_symb_id) const;
+
+    /* Locate the original index for the "untranslated" site, given a
+    reference index, i.e. the inverse operation of the translation matrix.*/
+    int get_original_index(int ref_indx) const;
+
     /**
@@ -93,0 +109,0 @@ Calculates the new energy given a set of system changes

cxx/include/cluster_list.hpp

 #ifndef CLUSTER_LIST_H
 #define CLUSTER_LIST_H
-#include <map>
 #include <set>
 #include <string>
+#include <unordered_map>
 #include <vector>
@@ -30,6 +30,6 @@
    private:
-    std::map<std::string, std::set<unsigned int> > symm_group_by_name;
-    std::map<std::string, std::vector<Cluster> > clusters;
+    std::unordered_map<std::string, std::set<unsigned int>> symm_group_by_name;
+    std::unordered_map<std::string, std::vector<Cluster>> clusters;
 };
 
 #endif

cxx/include/cluster.hpp

@@ -6,6 +6,6 @@ #ifndef CLUSTER_H
 #include <iostream>
-#include <map>
 #include <memory>
 #include <set>
 #include <string>
+#include <unordered_map>
 #include <vector>
@@ -15,3 +15,3 @@
 typedef std::vector<std::vector<int>> equiv_deco_t;
-typedef std::map<std::string, equiv_deco_t> all_equiv_deco_t;
+typedef std::unordered_map<std::string, equiv_deco_t> all_equiv_deco_t;
 
@@ -18,0 +18,0 @@ class Cluster {

cxx/include/row_sparse_struct_matrix.hpp

@@ -65,2 +65,10 @@ #ifndef ROW_SPARSE_STRUCT_MATRIX_H
 
+    const unsigned int get_num_non_zero() const {
+        return this->num_non_zero;
+    };
+
+    int *get_allowed_lookup_values() const {
+        return this->allowed_lookup_values;
+    };
+
    private:
@@ -67,0 +75,0 @@ int *allowed_lookup_values{nullptr};

cxx/src/basis_function.cpp

@@ -52,2 +52,15 @@ #include "basis_function.hpp"
 
+std::vector<std::pair<double, double>> BasisFunction::prepare_bfs_new_old(
+    unsigned int new_symb_id, unsigned int old_symb_id) const {
+    std::vector<std::pair<double, double>> prepped;
+    prepped.reserve(this->num_bf_values);
+
+    for (unsigned int dec_num = 0; dec_num < this->num_bf_values; dec_num++) {
+        double new_bf = this->get(dec_num, new_symb_id);
+        double old_bf = this->get(dec_num, old_symb_id);
+        prepped.emplace_back(std::make_pair(new_bf, old_bf));
+    }
+    return prepped;
+}
+
 void BasisFunction::swap(const BasisFunction &other) {
@@ -54,0 +67,0 @@ this->symb_ptr = other.symb_ptr;

cxx/src/cluster_list.cpp

 #include "cluster_list.hpp"
 
-using namespace std;
-
 void ClusterList::clear() {
@@ -15,7 +13,8 @@ clusters.clear();
 
-bool ClusterList::is_in_symm_group(const string &name, unsigned int symm_group) const {
-    return symm_group_by_name.at(name).find(symm_group) != symm_group_by_name.at(name).end();
+bool ClusterList::is_in_symm_group(const std::string &name, unsigned int symm_group) const {
+    auto group = symm_group_by_name.at(name);
+    return group.find(symm_group) != group.end();
 }
 
-const Cluster &ClusterList::get(const string &name, unsigned int symm_group) const {
+const Cluster &ClusterList::get(const std::string &name, unsigned int symm_group) const {
     for (const Cluster &cluster : clusters.at(name)) {
@@ -27,3 +26,3 @@ if (cluster.symm_group == symm_group) {
 
-    throw runtime_error("Did not find cluster in the requested symmetry group!");
+    throw std::runtime_error("Did not find cluster in the requested symmetry group!");
 }
@@ -43,3 +42,3 @@
 
-void ClusterList::unique_indices(set<int> &unique_indx) const {
+void ClusterList::unique_indices(std::set<int> &unique_indx) const {
     for (auto iter = clusters.begin(); iter != clusters.end(); ++iter) {
@@ -46,0 +45,0 @@ for (const Cluster &cluster : iter->second) {

cxx/src/row_sparse_struct_matrix.cpp

@@ -28,2 +28,3 @@ #include "row_sparse_struct_matrix.hpp"
     allowed_lookup_values = new int[num_non_zero];
+    this->num_rows = n_rows;
 }
@@ -30,0 +31,0 @@

cxx/src/with_numpy/ce_updater.cpp

@@ -230,3 +230,3 @@ #include "ce_updater.hpp"
 double CEUpdater::spin_product_one_atom(int ref_indx, const Cluster &cluster,
-                                        const vector<int> &dec, int ref_id) {
+                                        const vector<int> &dec, int ref_id) const {
     double sp = 0.0;
@@ -263,5 +263,20 @@
 
-double CEUpdater::spin_product_one_atom_delta(int ref_indx, const Cluster &cluster,
-                                              const vector<int> &dec, int old_symb_id,
-                                              int new_symb_id) {
+int CEUpdater::get_original_index(int ref_indx) const {
+    int *trans_matrix_row = trans_matrix.get_row(ref_indx);
+    int *allowed_lu = trans_matrix.get_allowed_lookup_values();
+    for (unsigned int j = 0; j < trans_matrix.get_num_non_zero(); j++) {
+        int col = allowed_lu[j];
+        int indx = trans_matrix.lookup_in_row(trans_matrix_row, col);
+        if (indx == ref_indx) {
+            return col;
+        }
+    }
+    std::stringstream err;
+    err << "Did not locate original index for ref index: " << ref_indx;
+    throw std::runtime_error(err.str());
+}
+
+double CEUpdater::spin_product_one_atom_delta(const SpinProductCache &sp_cache,
+                                              const Cluster &cluster,
+                                              const equiv_deco_t &equiv_deco) const {
     // Keep track of the change in spin-product
@@ -276,29 +291,36 @@ double sp_delta = 0.0;
 
-    // Cache the relevant row from the trans matrix.
-    int *trans_matrix_row = trans_matrix.get_row(ref_indx);
+    int *tm_row = sp_cache.trans_matrix_row;
 
-    for (unsigned int i = 0; i < num_indx; i++) {
-        // Calculate the spin product for both new and the old (ref)
-        double sp_temp_new = 1.0;
-        double sp_temp_ref = 1.0;
+    for (const vector<int> &dec : equiv_deco) {
+        // Iterate each site in the cluster
+        for (unsigned int i = 0; i < num_indx; i++) {
+            // Calculate the spin product for both new and the old (ref)
+            double sp_temp_new = 1.0;
+            double sp_temp_ref = 1.0;
+            // The constant term to the spin product from the sites which didn't change.
+            double sp_const = 1.0;
 
-        // Use pointer arithmetics in the inner most loop
-        const int *indices = &indx_list[i][0];
+            // Use pointer arithmetics in the inner most loop
+            const int *indices = &indx_list[i][0];
 
-        for (unsigned int j = 0; j < n_memb; j++) {
-            double bf_new, bf_ref;
-            int trans_index = trans_matrix.lookup_in_row(trans_matrix_row, indices[j]);
-            int dec_j = dec[j];
-            if (trans_index == ref_indx) {
-                bf_new = basis_functions.get(dec_j, new_symb_id);
-                bf_ref = basis_functions.get(dec_j, old_symb_id);
-            } else {
-                double bf = basis_functions.get(dec_j, symbols_with_id->id(trans_index));
-                bf_new = bf;
-                bf_ref = bf;
+            for (unsigned int j = 0; j < n_memb; j++) {
+                const int site_index = indices[j];
+                const int dec_j = dec[j];
+                if (site_index == sp_cache.original_index) {
+                    // This site is the reference index.
+                    // Look up the BF values directly for the new and old symbols
+                    const std::pair<double, double> &ref_bf = sp_cache.bfs_new_old[dec_j];
+                    sp_temp_new *= ref_bf.first;
+                    sp_temp_ref *= ref_bf.second;
+                } else {
+                    // Look up the symbol of the non-reference site, which hasn't changed.
+                    const int trans_index = trans_matrix.lookup_in_row(tm_row, site_index);
+                    sp_const *= basis_functions.get(dec_j, symbols_with_id->id(trans_index));
+                }
             }
-            sp_temp_new *= bf_new;
-            sp_temp_ref *= bf_ref;
+            // The change in spin-product is the difference in SP between the site(s) which changed,
+            // multiplied by the constant SP from the other un-changed sites (since these contribute
+            // equally before and after the change).
+            sp_delta += (sp_temp_new - sp_temp_ref) * sp_const * dup_factors[i];
         }
-        sp_delta += (sp_temp_new - sp_temp_ref) * dup_factors[i];
     }
@@ -322,2 +344,19 @@ return sp_delta;
 
+SpinProductCache CEUpdater::build_sp_cache(const SymbolChange &symb_change,
+                                           unsigned int old_symb_id,
+                                           unsigned int new_symb_id) const {
+    int ref_indx = symb_change.indx;
+    // Look up the untranslated index of the reference index.
+    int orig_indx = this->get_original_index(ref_indx);
+    // Cache the relevant row from the trans matrix.
+    int *trans_matrix_row = this->trans_matrix.get_row(ref_indx);
+
+    // Cache the basis functions for each decoration number specifically
+    // for the old and new symbols for faster lookup.
+    auto bfs_new_old = this->basis_functions.prepare_bfs_new_old(new_symb_id, old_symb_id);
+
+    SpinProductCache sp_cache = {ref_indx, orig_indx, trans_matrix_row, bfs_new_old};
+    return sp_cache;
+}
+
 void CEUpdater::update_cf(SymbolChange &symb_change) {
@@ -347,2 +386,9 @@ if (symb_change.old_symb == symb_change.new_symb) {
 
+    /* The following prepares a range of properties which will be
+     useful for all of the clusters, so we don't compute more
+     than we have to inside the main ECI loop (and especially inside the
+     spin-product loop). */
+
+    SpinProductCache sp_cache = this->build_sp_cache(symb_change, old_symb_id, new_symb_id);
+
     if (atoms != nullptr) {
@@ -371,8 +417,8 @@ set_symbol_in_atoms(atoms, symb_change.indx, symb_change.new_symb);
             unsigned int dec = parsed.dec_num;
-            next_cf[i] = current_cf[i] + (basis_functions.get(dec, new_symb_id) -
-                                          basis_functions.get(dec, old_symb_id)) /
-                                             num_non_bkg_sites;
+            const std::pair<double, double> &bf_ref = sp_cache.bfs_new_old[dec];
+            next_cf[i] = current_cf[i] + (bf_ref.first - bf_ref.second) / num_non_bkg_sites;
             continue;
         }
 
+        // n-body
         int symm = trans_symm_group[symb_change.indx];
@@ -390,7 +436,4 @@
 
-        double delta_sp = 0.0;
-        for (const vector<int> &deco : equiv_deco) {
-            delta_sp += spin_product_one_atom_delta(symb_change.indx, cluster, deco, old_symb_id,
-                                                    new_symb_id);
-        }
+        // Calculate the change (delta) in spin product
+        double delta_sp = spin_product_one_atom_delta(sp_cache, cluster, equiv_deco);
 
@@ -894,10 +937,8 @@ delta_sp *= (static_cast<double>(size) / equiv_deco.size());
             const Cluster &cluster = clusters.get(prefix, symm);
-            unsigned int size = cluster.size;
-            assert(cluster_indices[0].size() == size);
 
             const equiv_deco_t &equiv_deco = cluster.get_equiv_deco(dec_str);
+            unsigned int ref_id = symbols_with_id->id(atom_no);
             double sp_temp = 0.0;
             for (const vector<int> &deco : equiv_deco) {
-                sp_temp +=
-                    spin_product_one_atom(atom_no, cluster, deco, symbols_with_id->id(atom_no));
+                sp_temp += spin_product_one_atom(atom_no, cluster, deco, ref_id);
             }
@@ -904,0 +945,0 @@ sp += sp_temp / equiv_deco.size();

doc/source/releasenotes.rst

@@ -7,2 +7,14 @@ .. _releasenotes:
 
+1.0.4
+======
+* Added :class:`~clease.calculator.calc_cache.CleaseCacheCalculator`, as a primitive cache calculator
+  object with no cache validation.
+* Performance improvements to updating correlation functions.
+* Performance improvements to calculating the translation matrix, so the first
+  calculation of the clusters should be faster.
+* Performance improvements to the :class:`~clease.montecarlo.observers.lowest_energy_observer.LowestEnergyStructure`.
+  Correlation functions are also no longer tracked by default, but can be enabled with the ``track_cf`` key.
+* The default SGC observer in :class:`~clease.montecarlo.sgc_montecarlo.SGCMonteCarlo` should now be reset
+  automatically upon changing the temperature.
+
 1.0.3
@@ -9,0 +21,0 @@ ======

PKG-INFO

 Metadata-Version: 2.1
 Name: clease
-Version: 1.0.3
+Version: 1.0.4
 Summary: CLuster Expansion in Atomistic Simulation Environment
@@ -5,0 +5,0 @@ Home-page: https://gitlab.com/computationalmaterials/clease/

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

Improved metrics

Total package byte prevSize

1614335

0.62%

Number of package files

199

0.51%

Lines of code

16945

0.56%