Source code for clease.montecarlo.sgc_montecarlo

from typing import Sequence, Dict, Any
import numpy as np
from ase import Atoms
from ase.units import kB
from clease.calculator import Clease
from clease.settings import ClusterExpansionSettings
from .observers import SGCObserver
from .montecarlo import Montecarlo
from .trial_move_generator import TrialMoveGenerator, RandomFlip

class InvalidChemicalPotentialError(Exception):

[docs]class SGCMonteCarlo(Montecarlo): """ Class for running Monte Carlo in the Semi-Grand Canonical Ensebmle (i.e., fixed number of atoms, but varying composition) See the docstring of :class:`clease.montecarlo.Montecarlo` :param atoms: Atoms object (with CLEASE calculator attached!) :param temp: Temperature in kelvin :param symbols: Possible symbols to be used in swaps :param generator: Generator that produces trial moves """ def __init__( self, atoms: Atoms, temp: float, symbols: Sequence[str] = (), generator: TrialMoveGenerator = None, observe_singlets: bool = False, ): if not isinstance(atoms, Atoms): raise ValueError(f"atoms must be an Atoms object, got {atoms!r}") if not isinstance(atoms.calc, Clease): raise ValueError( f"Atoms must have a Clease calculator object attached, got {atoms.calc!r}." ) if generator is None: if len(symbols) <= 1: raise ValueError("At least 2 symbols have to be specified") # Only select indices which are not considered background. non_bkg = atoms.calc.settings.non_background_indices generator = RandomFlip(symbols, atoms, indices=non_bkg) self.averager = SGCObserver(atoms.calc, observe_singlets=observe_singlets) super().__init__(atoms, temp, generator=generator) self.symbols = symbols self.chem_pots = [] self.chem_pot_names = [] self.has_attached_avg = False = "SGCMonteCarlo" self._chemical_potential = None self.chem_pot_in_eci = False has_attached_obs = False for obs in self.iter_observers(): if == "SGCObserver": has_attached_obs = True self.averager = obs break if not has_attached_obs: self.attach(self.averager) @property def observe_singlets(self) -> bool: return self.averager.observe_singlets def _check_symbols(self): """ Override because there are no restriction on the symbols here """
[docs] def reset(self): """ Reset the simulation object """ super().reset() self.averager.reset()
@property def calc(self) -> Clease: return self.atoms.calc @property def settings(self) -> ClusterExpansionSettings: return self.calc.settings @property def chemical_potential(self): return self._chemical_potential @chemical_potential.setter def chemical_potential(self, chem_pot: Dict[str, float]): eci = self.calc.eci if any(key not in eci for key in chem_pot): msg = "A chemical potential not being trackted is added. Make " msg += "sure that all the following keys are in the ECIs before " msg += "they are passed to the calculator: " msg += f"{list(chem_pot.keys())}\n" msg += "(Add them with a zero ECI value if they are not supposed " msg += "to be included.)" raise InvalidChemicalPotentialError(msg) self._chemical_potential = chem_pot if self.chem_pot_in_eci: self._reset_eci_to_original(self.calc.eci) self._include_chemical_potential_in_eci(chem_pot, self.calc.eci) def _include_chemical_potential_in_eci(self, chem_pot: Dict[str, float], eci: Dict[str, float]): """ Including the chemical potentials in the ECIs Parameters: chem_pot: dict Chemical potentials eci: dict Original ECIs """ self.chem_pots = [] self.chem_pot_names = [] keys = list(chem_pot.keys()) keys.sort() for key in keys: self.chem_pots.append(chem_pot[key]) self.chem_pot_names.append(key) current_eci = eci.get(key, 0.0) eci[key] = current_eci - chem_pot[key] calc = self.calc calc.update_eci(eci) self.chem_pot_in_eci = True self.current_energy = calc.calculate(None, None, None) return eci def _reset_eci_to_original(self, eci_with_chem_pot: Dict[str, float]): """ Resets the ECIs to their original value :parma dict eci_with_chem_pot: ECIs with chemical potential included """ for name, val in zip(self.chem_pot_names, self.chem_pots): eci_with_chem_pot[name] += val calc = self.calc calc.update_eci(eci_with_chem_pot) self.chem_pot_in_eci = False self.current_energy = calc.calculate(None, None, None) return eci_with_chem_pot
[docs] def reset_eci(self): """Return the ECIs.""" if self.chem_pot_in_eci: self._reset_eci_to_original(self.calc.eci)
[docs] def run(self, steps: int = 10, call_observers: bool = True, chem_pot: Dict[str, float] = None): """ Run Monte Carlo simulation. See :py:meth:`` Parameters: chem_pot: dict Chemical potentials. The keys should correspond to one of the singlet terms. A typical form of this is {"c1_0":-1.0,c1_1_1.0} """ if chem_pot is None and self.chemical_potential is None: ex_chem_pot = {"c1_1": -0.1, "c1_2": 0.05} raise ValueError( f"No chemical potentials given. Has to be a dictionary of the form {ex_chem_pot}" ) if chem_pot is not None: self.chemical_potential = chem_pot self.reset() super().run(steps=steps, call_observers=call_observers)
[docs] def singlet2composition(self, avg_singlets: Dict[str, float]): """Convert singlets to composition.""" bf = self.settings.basis_functions matrix = np.zeros((len(self.symbols), len(self.symbols))) index = {s: i for i, s in enumerate(self.symbols)} for i, b in enumerate(bf): for s, col in index.items(): matrix[i, col] = b[s] matrix[-1, :] = 1.0 rhs = np.zeros(len(self.symbols)) rhs[:-1] = avg_singlets rhs[-1] = 1.0 x = np.linalg.solve(matrix, rhs) res = {} for s, i in index.items(): name = f"{s}_conc" res[name] = x[i] return res
[docs] 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()
[docs] def get_thermodynamic_quantities(self, reset_eci: bool = False) -> Dict[str, Any]: """Compute thermodynamic quantities. Parameters: reset_eci: bool If True, the chemical potential will be removed from the ECIs. """ # Note - in order to correctly measure averagers from the SGC observer, # we need some information from the SGC MC object. So we directly extract the averages # from that observer here. avg_obs = self.averager # Alias N = self.averager.counter averages = {} averages["energy"] = averages["sgc_energy"] = averages["sgc_heat_capacity"] = avg_obs.energy_sq.mean -**2 averages["sgc_heat_capacity"] /= kB * self.temperature**2 averages["temperature"] = self.temperature averages["n_mc_steps"] = self.averager.counter averages["accept_rate"] = self.current_accept_rate # Singlets are more expensive to measure than the other quantities, # so only measure them upon request. if self.observe_singlets: # Add singlets and chemical potential to the dictionary # pylint: disable=consider-using-enumerate singlets = avg_obs.singlets / N singlets_sq = avg_obs.quantities["singlets_sq"] / N averages["singlet_energy"] = natoms = len(self.atoms) for i, chem_pot in enumerate(self.chem_pots): averages["singlet_energy"] += chem_pot * singlets[i] * natoms for i in range(len(singlets)): name = f"singlet_{self.chem_pot_names[i]}" averages[name] = singlets[i] name = f"var_singlet_{self.chem_pot_names[i]}" averages[name] = singlets_sq[i] - singlets[i] ** 2 # Measure concentration from the singlets try: avg_conc = self.singlet2composition(singlets) averages.update(avg_conc) # pylint: disable=broad-except except Exception as exc: print("Could not find average singlets!") print(exc) # Always measure the chemical potentials. for chem_pot_name, chem_pot in zip(self.chem_pot_names, self.chem_pots): key = f"mu_{chem_pot_name}" averages[key] = chem_pot averages.update(self.meta_info) # Add information from observers averages.update(self._get_obs_averages()) if reset_eci: self._reset_eci_to_original(self.atoms.calc.eci) return averages