Functionality for calculating energy changes

Description

For efficient Monte Carlo simulations it is pivotal that energy changes be evaluated locally. That is to say, when performing a trial step one should not compute the total energy of the entire system but rather the energy difference induced by the trial step (which affects only a certain part of the system).

The objective of this issue is to discuss a suitable strategy for how to implement this feature. (Note that the first implementation of the canonical ensemble (#163 (closed)) is scheduled to use the total energy approach, to be changed later.)

Background

Generally ensembles will carry out a trial move that induces a configurational change. The move is accepted or rejected depending on a condition associated with the ensemble (most often the Metropolis criterion), which involves computing the change in the respective observable.

Examples

• Canonical ensemble
  • trial move: select two sites with unalike species (decoration) and swap their decoration (A<-->B)
  • acceptance probability: \phi = \min\{ 1, \exp[-\Delta E/k_B T] \}, where \Delta E is the change of the total energy and T is the temperature (ensemble parameter)
• Semi-grand canonical (SGC) ensemble
  • trial move: select one site and change its decoration (A-->B or B-->A)
  • acceptance probability: \phi = \min\{ 1, \exp[-\Delta U/k_B T] \} with \Delta U = \Delta E + \Delta\mu\Delta c, where \Delta E is the change of the total energy, \Delta c is the change in the concentration, T is the temperature (ensemble parameter), and \Delta\mu is the difference in chemical potentials between A and B (ensemble parameter)
• Variance constrained semi-grand canonical (VCSGC) ensemble Phys. Rev. B 86, 134204 (2012)
  • trial move: select one site and change its decoration (A-->B or B-->A); identical to SGC ensemble
  • acceptance probability: \phi = \min\{ 1, \exp[-\Delta U/k_B T] \} with U=E+N\bar{\kappa}\left(c+\bar{\phi}/2\right)^2, where \Delta E is the change of the total energy, \Delta c is the change in the concentration, T is the temperature (ensemble parameter), and \bar{\phi} as well as \bar{\kappa} are ensemble parameters that control average and variance of the concentration, respectively
• Svendsen-Wang algorithm; strongly simplified here
  • trial move: identify clusters of sites with the same spin (or decoration)
  • acceptance probability: \phi = \min\{ 1, \exp[-\Delta E/k_B T] \}, where \Delta E is the change of the total energy and T is the temperature (ensemble parameter)

Each of these ensembles should be implemented for two or more components.

Notes

For any of these ensembles to be sampled efficiently one should compute \Delta E based on a local change in the spin vector \vec{\sigma}

\Delta E = f(\Delta \vec{\sigma}) = f(\vec{\sigma}_\text{trial} - \vec{\sigma}_\text{current})

rather than from the brute force total energy difference

\Delta E = \mathcal{H}(\vec{\sigma}_\text{trial}) - \mathcal{H}(\vec{\sigma}_\text{current}),

where \mathcal{H} is the Hamiltonian of the system.

Suggested implementation

The following structure suggestion has emerged from rather extensive IRL discussions involving @angqvist @williamz @erikfransson and @erhart.

icet / calculators / base_calculator
                   / cluster_expansion_calculator
                   / some_other_calculator

# icet/calculators/base_calculator.py

class BaseCalculator:
  def __init__(self, supercell, ...):
    # the following internal variables should be exposed as read-only
    # properties that can be queried by the ensemble to ensure that
    # the calculator is compatible with the trial moves of the ensemble
    self._displacements = None  # True for e.g., displacement MC
    self._swap_species = None   # True for e.g., SGC ensemble MC
  def calculate_total(self, configuration):
    ...
  def calculate_local_contribution(self, sites):
    """
    Calculate the sum of the contributions from clusters
    that involve any of the sites specified.

    Parameters
    ----------
    sites : list of ints
        sites for which to calculate the contribution to the property

    Returns
    -------
    float

    Example
    -------
    This function can be used to compute the energy change e.g., in the
    case of the canonical ensemble as follows::
        dE = calculate_local_contribution([index1, index2])
        # carry out trial move
        dE -= calculate_local_contribution([index1, index2])
        # accept/reject
    ...

# icet/calculators/cluster_expansion_calculator.py

class ClusterExpansionCalculator(BaseCalculator):
  def __init__(self, supercell, cluster_expansion):
    # type(supercell) = ase.Atoms
    # type(cluster_expansion) = icet.ClusterExpansion
    # call constructor of base class

# icet/ensembles/some_ensemble.py

class CanonicalEnsemble(BaseEnsemble):
  def __init__(self, calculator, ...):
    assert isinstance(calculator, BaseCalculator), \
      'calculator must be instance of BaseCalculator'
    ...
  def _trial_move(self,...):
    dE = calculate_local_contribution([index1, index2])
    # carry out trial move
    dE -= calculate_local_contribution([index1, index2])
    # accept/reject