Source code for clease.montecarlo.montecarlo

"""Monte Carlo method for ase."""
from typing import Dict, Union, Iterator, Any
from datetime import datetime
import time
import logging
import random
import math
from collections import Counter
from ase import Atoms
from ase.units import kB
from clease.datastructures import SystemChanges, MCStep
from .mc_evaluator import CEMCEvaluator, MCEvaluator
from .base import BaseMC
from .averager import Averager
from .bias_potential import BiasPotential
from .observers import MCObserver
from .trial_move_generator import TrialMoveGenerator, RandomSwap

logger = logging.getLogger(__name__)

# pylint: disable=too-many-instance-attributes
[docs]class Montecarlo(BaseMC): """Class for running Monte Carlo at a fixed composition (canonical). For more information, also see the documentation of the parent class :class:`~clease.montecarlo.base.BaseMC`. Args: system (Union[ase.Atoms, MCEvaluator]): Either an ASE Atoms object with an attached calculator, or a pre-initialized :class:`~clease.montecarlo.mc_evaluator.MCEvaluator` object. temp (float): Temperature of Monte Carlo simulation in Kelvin generator (TrialMoveGenerator, optional): A :class:`~clease.montecarlo.trial_move_generator.TrialMoveGenerator` object that produces trial moves. Defaults to None. """ NAME = "MonteCarlo" def __init__( self, system: Union[Atoms, MCEvaluator], temp: float, generator: TrialMoveGenerator = None, ): # We cannot cause an energy calculation trigger in init, # so we defer these quantities until needed. self.current_energy = None self.new_energy = None self.mean_energy = None self.energy_squared = None super().__init__(system, temp) if generator is None: self.generator = _make_default_swap_generator(self.evaluator) else: self.generator = generator # List of observers that will be called every n-th step # similar to the ones used in the optimization routines self.observers = [] self.bias_potentials = [] self.current_step = 0 self.num_accepted = 0 self.status_every_sec = 30 self.trial_move = [] self.quit = False def _on_temp_change(self) -> None: """Reset the energy averagers after a change in temperature""" self.reset_averagers() self._reset_internal_counters() def update_current_energy(self) -> None: self.current_energy = self.evaluator.get_energy() self.current_energy += sum( bias.calculate_from_scratch(self.atoms) for bias in self.bias_potentials ) logger.debug("Updating current energy to %s", self.current_energy) self.evaluator.reset() def _initialize_energies(self) -> None: """Initialize the energy quantities, causes a trigger of an energy evaluation.""" self.update_current_energy() self.new_energy = self.current_energy # Create new averager objects, only the first time # (since they are initialized as None in __init__) if self.mean_energy is None: self.mean_energy = Averager(ref_value=self.current_energy) if self.energy_squared is None: self.energy_squared = Averager(ref_value=self.current_energy)
[docs] def reset(self) -> None: """Reset all member variables to their original values.""" logger.debug("Resetting.") for _, obs in self.observers: obs.reset() self.evaluator.reset() self._reset_internal_counters() self.reset_averagers()
[docs] def reset_averagers(self) -> None: """Reset the energy averagers.""" # Averagers are initialized to None in the constructor. for averager in (self.mean_energy, self.energy_squared): if averager is not None: averager.clear()
def _reset_internal_counters(self) -> None: """Reset the step counters which are used internally""" self.current_step = 0 self.num_accepted = 0
[docs] def add_bias(self, potential: BiasPotential): """Add a new bias potential. Parameters: potential: Potential to be added """ if not isinstance(potential, BiasPotential): raise TypeError("potential has to be of type BiasPotential") self.bias_potentials.append(potential)
[docs] def attach(self, obs: MCObserver, interval: int = 1): """Attach observers to be called on a given MC step interval. Parameters: obs: MCObserver Observer to be added interval: int How often the observer should be called """ if not obs.interval_ok(interval): name = type(obs).__name__ raise ValueError(f"Invalid interval for {name}. Check docstring of the observer.") self.observers.append((interval, obs))
[docs] 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
[docs] def initialize_run(self): """Prepare MC object for a new run.""" logger.debug("Initializing run") self.generator.initialize(self.atoms) # Ensure the evaluator is properly synchronized. self.evaluator.synchronize() # Initialize/update relevant energy quantities self._initialize_energies() # Reset the internal step counters self._reset_internal_counters()
[docs] def run(self, steps: int = 100, call_observers: bool = True) -> None: """Run Monte Carlo simulation. Parameters: steps: int Number of steps in the MC simulation call_observers: bool Should the observers be called during this run? Can be turned off for running burn-ins. The energy averagers will still be updated, even if this flag is disabled. Defaults to True. """ # Construct the iterator, make the preparations for starting the run mc_iter = self.irun(steps, call_observers=call_observers) start = time.perf_counter() prev = self.current_step info_enabled = logger.isEnabledFor(logging.INFO) # Do we emit INFO logs? for _ in mc_iter: # We only want to do this calculation if logging is enabled for INFO. if info_enabled and time.perf_counter() - start > self.status_every_sec: ms_per_step = 1000.0 * self.status_every_sec / (self.current_step - prev) "%d of %d steps. %.2f ms per step. Acceptance rate: %.2f", self.current_step, steps, ms_per_step, self.current_accept_rate, ) prev = self.current_step start = time.perf_counter()
[docs] def irun(self, steps: int, call_observers: bool = True) -> Iterator[MCStep]: """Run Monte Carlo simulation as an iterator. Can be used to inspect the MC after each step, for example, to print the energy every 5 steps, one could do: >>> mc = Montecarlo(...) # doctest: +SKIP >>> for mc_step in mc.irun(500): # doctest: +SKIP ... if mc_step.step % 5 == 0: ... print(f"Current energy: {} eV") The iterator yields individual instances of :class:`~clease.datastructures.mc_step.MCStep` for each step which is taken. Parameters: steps: int Number of steps in the MC simulation call_observers: bool Should the observers be called during this run? Can be turned off for running burn-ins. The energy averagers will still be updated, even if this flag is disabled. Defaults to True. """"Starting MC run with %d steps.", steps) # We first ensure the MC object is initialized, and then we # construct the iterator. # This ensures everything is prepared before the first step # is consumed from the iterator. self.initialize_run() # Now we create the iterator return self._irun(steps, call_observers=call_observers)
def _irun(self, steps: int, call_observers: bool = True) -> Iterator[MCStep]: """Create the MC iterator""" # Offset range by 1, so that we start with current step = 1 for _ in range(steps): step = self._mc_step(call_observers=call_observers) en = self.mean_energy += en # E * E is slightly faster than E ** 2 self.energy_squared += en * en if self.quit: # Breaking out will cause the iterator to end. logger.warning("Quit signal detected. Breaking the MC loop.") break # Yield, allow inspection before continuing. yield step else: # Loop was not broken, we reached the max number of steps"Reached maximum number of steps (%d mc steps)", steps) logger.debug("Reached end of MC iterator.") @property def meta_info(self) -> Dict[str, str]: """Return dict with meta info.""" # Get the timestamp with millisecond precision timestamp ="milliseconds") meta_info = { "timestamp": timestamp, } return meta_info @property def current_accept_rate(self) -> float: """Return the current accept rate as a value between 0 and 1.""" if self.current_step == 0: # No steps have been taken yet. return 0.0 return self.num_accepted / self.current_step
[docs] def get_thermodynamic_quantities(self) -> Dict[str, Any]: """Compute thermodynamic quantities.""" quantities = {} mean_energy = self.mean_energy.mean quantities["energy"] = mean_energy mean_sq = self.energy_squared.mean quantities["heat_capacity"] = (mean_sq - mean_energy**2) / (kB * self.temperature**2) quantities["energy_var"] = mean_sq - mean_energy**2 quantities["temperature"] = self.temperature quantities["accept_rate"] = self.current_accept_rate quantities["n_mc_steps"] = self.current_step at_count = self.count_atoms() for key, value in at_count.items(): name = f"{key}_conc" conc = value / len(self.atoms) quantities[name] = conc # Add some more info that can be useful quantities.update(self.meta_info) # Add information from observers quantities.update(self._get_obs_averages()) return quantities
def _get_obs_averages(self) -> Dict[str, Any]: """Get average measurements from observers""" obs_avgs = {} for obs in self.iter_observers(): obs_avgs.update(obs.get_averages()) return obs_avgs def _calculate_step(self, system_changes: SystemChanges): """Calculate energies given a step, and decide if we accept the step. Returns boolean if system changes are accepted. Parameters: system_changes: list List with system changes """ self.evaluator.apply_system_changes(system_changes) # Evaluate the energy quantity after applying the changes to the system. self.new_energy = self.evaluator.get_energy(applied_changes=system_changes) # NOTE: As this is called after calculate, the changes has # already been introduced to the system for bias in self.bias_potentials: self.new_energy += bias(system_changes) accept = self._do_accept(self.current_energy, self.new_energy) if accept: # Changes accepted, finalize evaluator. self.evaluator.keep_system_changes(system_changes) else: # Undo changes self.evaluator.undo_system_changes(system_changes) return accept def _do_accept(self, current_energy: float, new_energy: float) -> bool: """Decide if we accept a state, based on the energies. Return a bool on whether the move was accepted. Parameters: :param current_energy: Energy of the current configuration :param new_energy: Energy of the new configuration. """ # Standard Metropolis acceptance criteria if new_energy < current_energy: return True energy_diff = new_energy - current_energy probability = math.exp(-energy_diff / self.kT) return random.random() <= probability def _move_accepted(self, system_changes: SystemChanges) -> None: self.num_accepted += 1 self.generator.on_move_accepted(system_changes) self.current_energy = self.new_energy def _move_rejected(self, system_changes: SystemChanges) -> None: self.generator.on_move_rejected(system_changes)
[docs] def count_atoms(self) -> Dict[str, int]: """Count the number of each element.""" return dict(Counter(self.atoms.symbols))
def _mc_step(self, call_observers: bool = True) -> MCStep: """Make one Monte Carlo step by swithing two atoms.""" self.current_step += 1 system_changes = self.generator.get_trial_move() self.trial_move = system_changes # Calculate step, and whether we accept it move_accepted = self._calculate_step(system_changes) updater = self._move_accepted if move_accepted else self._move_rejected updater(system_changes) step = MCStep(self.current_step, self.current_energy, move_accepted, system_changes) if call_observers: # Execute all observers self.execute_observers(step) return step def execute_observers(self, last_step: MCStep): for interval, obs in self.observers: if self.current_step % interval == 0: obs.observe_step(last_step)
def _make_default_swap_generator(evaluator: MCEvaluator) -> RandomSwap: """Construct the default RandomSwap trial move generator. If the evaluator object is not a Cluster Expansion evaluator, then any swaps are allowed. If the evaluator is a Cluster Expansion MC evaluator object, then only non-background sites are considered for random swaps. This does *not* constrain swaps between sublattices. """ atoms = evaluator.atoms if not isinstance(evaluator, CEMCEvaluator): # Running a non-CE Evaluator return RandomSwap(atoms) # We're a CE evaluator, constrain the background sites non_bkg_indices = evaluator.settings.non_background_indices return RandomSwap(atoms, indices=non_bkg_indices)