Ensembles#

Canonical ensemble#

class mchammer.ensembles.CanonicalEnsemble(structure, calculator, temperature, user_tag=None, boltzmann_constant=8.617330337217213e-05, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allow one to simulate systems in the canonical ensemble (\(N_iVT\)), i.e. at constant temperature (\(T\)), number of atoms of each species (\(N_i\)), and volume (\(V\)).

The probability for a particular state in the canonical ensemble is proportional to the well-known Boltzmann factor,

\[\rho_{\text{C}} \propto \exp [ - E / k_B T ].\]

Since the concentrations or equivalently the number of atoms of each species is held fixed in the canonical ensemble, a trial step must conserve the concentrations. This is accomplished by randomly picking two unlike atoms and swapping their identities. The swap is accepted with probability

\[P = \min \{ 1, \, \exp [ - \Delta E / k_B T ] \},\]

where \(\Delta E\) is the change in potential energy caused by the swap.

The canonical ensemble provides an ideal framework for studying the properties of a system at a specific concentration. Properties such as potential energy or phenomena such as chemical ordering at a specific temperature can conveniently be studied by simulating at that temperature. The canonical ensemble is also a convenient tool for “optimizing” a system, i.e., finding its lowest energy chemical ordering. In practice, this is usually achieved by simulated annealing, i.e., the system is equilibrated at a high temperature, after which the temperature is continuously lowered until the acceptance probability is almost zero. In a well-behaved system, the chemical ordering at that point corresponds to a low-energy structure, possibly the global minimum at that particular concentration.

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • temperature (float) – Temperature \(T\) in appropriate units, commonly Kelvin.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e., units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random swap. This should be as long as the number of sublattices and should sum up to 1.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again.

Example

The following snippet illustrate how to carry out a simple Monte Carlo simulation in the canonical ensemble. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet ECIs are zero and only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

>>> # prepare initial configuration
>>> structure = prim.repeat(3)
>>> for k in range(5):
>>>     structure[k].symbol = 'Ag'

>>> # set up and run MC simulation
>>> calc = ClusterExpansionCalculator(structure, ce)
>>> mc = CanonicalEnsemble(structure=structure, calculator=calc,
...                        temperature=600,
...                        dc_filename='myrun_canonical.dc')
>>> mc.run(100)  # carry out 100 trial swaps
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run(number_of_trial_steps)#

Samples the ensemble for the given number of trial steps.

Observations are made whenever the trial step counter is a multiple of the respective interval, so the observation steps stay aligned when run() is called several times. The data container is written to file when the write period has elapsed and at the end of the run, if a file name was given.

Parameters:

number_of_trial_steps (int) – Number of trial steps to carry out.

Raises:

TypeError – If number_of_trial_steps is not an int.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Current temperature.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Canonical annealing#

class mchammer.ensembles.CanonicalAnnealing(structure, calculator, T_start, T_stop, n_steps, cooling_function='exponential', user_tag=None, boltzmann_constant=8.617330337217213e-05, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allow one to carry out simulated annealing in the canonical ensemble, i.e. the temperature is varied in pre-defined fashion while the composition is kept fixed. See CanonicalEnsemble for more information about the standard canonical ensemble.

The canonical annealing ensemble can be useful, for example, for finding ground states or generating low energy configurations.

The temperature control scheme is selected via the cooling_function keyword argument, while the initial and final temperature are set via the T_start and T_stop arguments. Several pre-defined temperature control schemes are available including 'linear' and 'exponential'. In the latter case the temperature varies logarithmically as a function of the MC step, emulating the exponential temperature dependence of the atomic exchange rate encountered in many materials. It is also possible to provide a user defined cooling function via the keyword argument. This function must comply with the following function header:

def cooling_function(step, T_start, T_stop, n_steps):
    T = ...  # compute temperature
    return T

Here step refers to the current MC trial step.

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (ClusterExpansionCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • T_start (float) – Temperature from which the annealing is started.

  • T_stop (float) – Final temperature for annealing.

  • n_steps (int) – Number of steps to take in the annealing simulation.

  • cooling_function (str) – to use the predefined cooling functions provide a string 'linear' or 'exponential', otherwise provide a function.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e., units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random swap. This should be as long as the number of sublattices and should sum up to 1.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again.

property T_start: float#

Starting temperature.

property T_stop: float#

Final temperature.

attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

property estimated_ground_state: Atoms#

Structure with lowest observed potential during run.

property estimated_ground_state_potential: float#

Lowest observed potential during run.

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property n_steps: int#

Number of steps to carry out.

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run()[source]#

Runs the annealing simulation.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Current temperature.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Semi-grand canonical ensemble#

class mchammer.ensembles.SemiGrandCanonicalEnsemble(structure, calculator, temperature, chemical_potentials, boltzmann_constant=8.617330337217213e-05, user_tag=None, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allow one to simulate systems in the semi-grand canonical (SGC) ensemble (\(N\Delta\mu_i VT\)), i.e. at constant temperature (\(T\)), total number of sites (\(N=\sum_i N_i\)), relative chemical potentials (\(\Delta\mu_i=\mu_i - \mu_1\), where \(i\) denotes the species), and volume (\(V\)).

The probability for a particular state in the SGC ensemble for a \(m\)-component system can be written

\[\rho_{\text{SGC}} \propto \exp\Big[ - \big( E - \sum_{i>1}^m \Delta\mu_i N_i \big) \big / k_B T \Big]\]

with the relative chemical potentials \(\Delta\mu_i = \mu_i - \mu_1\) and species counts \(N_i\). Unlike the canonical ensemble, the number of the respective species (or equivalently the concentrations) are allowed to vary in the SGC ensemble. A trial step thus consists of randomly picking an atom and changing its identity with probability

\[P = \min \Big\{ 1, \, \exp \big[ - \big( \Delta E - \sum_i \Delta \mu_i \Delta N_i \big) \big / k_B T \big] \Big\},\]

where \(\Delta E\) is the change in potential energy caused by the swap.

There exists a simple relation between the differences in chemical potential and the canonical free energy \(F\). In a binary system this relationship reads

\[\Delta \mu = \frac{1}{N} \frac{\partial F}{\partial c} ( N, V, T, \langle c \rangle).\]

Here \(c\) denotes concentration (\(c=N_i/N\)) and \(\langle c \rangle\) the average concentration observed in the simulation. By recording \(\langle c \rangle\) while gradually changing \(\Delta \mu\), one can thus in principle calculate the difference in canonical free energy between the pure phases (\(c=0\) or \(1\)) and any concentration by integrating \(\Delta \mu\) over that concentration range. In practice this requires that the recorded average concentration \(\langle c \rangle\) varies continuously with \(\Delta \mu\). This is not the case for materials with multiphase regions (such as miscibility gaps), because in such regions \(\Delta \mu\) maps to multiple concentrations. In a Monte Carlo simulation, this is typically manifested by discontinuous jumps in concentration. Such jumps mark the phase boundaries of a multiphase region and can thus be used to construct the phase diagram. To recover the free energy, however, such systems require sampling in other ensembles, such as the variance-constrained semi-grand canonical ensemble.

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • temperature (float) – Temperature \(T\) in appropriate units, commonly Kelvin.

  • chemical_potentials (dict[str, float]) – Chemical potential for each species \(\mu_i\). The key denotes the species, the value specifies the chemical potential in units that are consistent with the underlying cluster expansion.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e. units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random flip. This should be as long as the number of sublattices and should sum up to 1.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again.

Example

The following snippet illustrate how to carry out a simple Monte Carlo simulation in the semi-canonical ensemble. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet ECIs are zero and only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

>>> # set up and run MC simulation (T=600 K, delta_mu=0.8 eV/atom)
>>> structure = prim.repeat(3)
>>> calc = ClusterExpansionCalculator(structure, ce)
>>> mc = SemiGrandCanonicalEnsemble(structure=structure, calculator=calc,
...                                temperature=600,
...                                dc_filename='myrun_sgc.dc',
...                                chemical_potentials={'Ag': 0, 'Au': 0.8})
>>> mc.run(100)  # carry out 100 trial swaps
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property chemical_potentials: dict[int, float]#

Chemical potentials \(\mu_i\) (see parameters section above).

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run(number_of_trial_steps)#

Samples the ensemble for the given number of trial steps.

Observations are made whenever the trial step counter is a multiple of the respective interval, so the observation steps stay aligned when run() is called several times. The data container is written to file when the write period has elapsed and at the end of the run, if a file name was given.

Parameters:

number_of_trial_steps (int) – Number of trial steps to carry out.

Raises:

TypeError – If number_of_trial_steps is not an int.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Temperature \(T\) (see parameters section above).

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Semi-grand canonical annealing#

class mchammer.ensembles.SGCAnnealing(structure, calculator, T_start, T_stop, n_steps, chemical_potentials, cooling_function='exponential', boltzmann_constant=8.617330337217213e-05, user_tag=None, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allow one to carry out simulated annealing in the semi grand canonical ensemble, i.e., the temperature is varied in pre-defined fashion while the chemical potential is kept fixed. See mchammer.ensembles.SemiGrandCanonicalEnsemble for more information about the ensemble.

The temperature control scheme is selected via the cooling_function keyword argument, while the initial and final temperature are set via the T_start and T_stop arguments. Several pre-defined temperature control schemes are available including 'linear' and 'exponential'. In the latter case the temperature varies logarithmically as a function of the MC step, emulating the exponential temperature dependence of the atomic exchange rate encountered in many materials. It is also possible to provide a user defined cooling function via the keyword argument. This function must comply with the following function header:

def cooling_function(step, T_start, T_stop, n_steps):
    T = ...  # compute temperature
    return T

Here step refers to the current MC trial step.

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • chemical_potentials (dict[str, float]) – Chemical potential for each species \(\mu_i\). The key denotes the species, the value specifies the chemical potential in units that are consistent with the underlying cluster expansion.

  • calculator (ClusterExpansionCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • T_start (float) – Temperature from which the annealing is started.

  • T_stop (float) – Final temperature for annealing.

  • n_steps (int) – Number of steps to take in the annealing simulation.

  • cooling_function (str) – to use the predefined cooling functions provide a string 'linear' or 'exponential', otherwise provide a function.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e. units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random swap. This should be as long as the number of sublattices and should sum up to 1.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again.

attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property chemical_potentials: dict[int, float]#

Chemical potentials \(\mu_i\) (see parameters section above).

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

property estimated_ground_state: Atoms#

Structure with lowest observed potential during run.

property estimated_ground_state_potential: float#

Lowest observed potential during run.

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run()[source]#

Runs the annealing.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Current temperature.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Variance-constrained semi-grand canonical ensemble#

class mchammer.ensembles.VCSGCEnsemble(structure, calculator, temperature, phis, kappa, boltzmann_constant=8.617330337217213e-05, user_tag=None, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allow one to simulate systems in the variance-constrained semi-grand canonical (VCSGC) ensemble (\(N\phi\kappa VT\)), i.e. at constant temperature (\(T\)), total number of sites (\(N=\sum_i N_i\)), and two additional dimensionless parameters \(\phi\) and \(\kappa\), which constrain average and variance of the concentration, respectively.

The below examples treat the binary case, but the generalization of to ternaries and higher-order systems is straight-forward. The probability for a particular state in the VCSGC ensemble for a \(2\)-component system can be written

\[\rho_{\text{VCSGC}} \propto \exp\Big[ - E / k_B T + \kappa N ( c_1 + \phi_1 / 2 )^2 \Big],\]

where \(c_1\) represents the concentration of species 1, i.e., \(c_1=N_1/N\). (Please note that the quantities \(\kappa\) and \(\phi\) correspond, respectively, to \(\bar{\kappa}\) and \(\bar{\phi}\) in [SadErh12].) The \(\phi\) may refer to any of the two species. If \(\phi\) is specified for species A, an equivalent simulation can be carried out by specifying \(\phi_B\) as \(-2-\phi_A\). In general, simulations of \(N\)-component systems requires the specification of \(\phi\) for \(N-1\) elements.

Just like the semi-grand canonical ensemble, the VCSGC ensemble allows concentrations to change. A trial step consists of changing the identity of a randomly chosen atom and accepting the change with probability

\[P = \min \{ 1, \, \exp [ - \Delta E / k_B T + \kappa N \Delta c_1 (\phi_1 + \Delta c_1 + 2 c_1 ) ] \}.\]

Note that for a sufficiently large value of \(\kappa\), say 200, the probability density \(\rho_{\text{VCSGC}}\) is sharply peaked around \(c_1=-\phi_1 / 2\). In practice, this means that we can gradually change \(\phi_1\) from (using some margins) \(-2.1\) to \(0.1\) and take the system continuously from \(c_1 = 0\) to \(1\). The parameter \(\kappa\) constrains the fluctuations (or the variance) of the concentration at each value of \(\phi_1\), with higher values of \(\kappa\) meaning less fluctuations. Unlike the semi-grand canonical ensemble, one value of \(\phi_1\) maps to one and only one concentration also in multiphase regions. Since the derivative of the canonical free energy can be expressed in terms of parameters and observables of the VCSGC ensemble,

\[k_B T \kappa ( \phi_1 + 2 \langle c_1 \rangle ) = - \frac{1}{N} \frac{\partial F}{\partial c_1} (N, V, T, \langle c_1 \rangle ),\]

this ensemble allows for thermodynamic integration across multiphase regions. This means that we can construct phase diagrams by directly comparing the free energies of the different phases. This often makes the VCSGC ensemble more convenient than the semi-grand canonical ensemble when simulating materials with multiphase regions, such as alloys with miscibility gaps.

When using the VCSGC ensemble, please cite Phys. Rev. B 86, 134204 (2012) [SadErh12].

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • temperature (float) – Temperature \(T\) in appropriate units, commonly Kelvin.

  • phis (dict[str, float]) – Average constraint parameters \(\phi_i\). The key denotes the species. For a \(N\)-component sublattice, there should be \(N - 1\) different \(\phi_i\) (referred to as \(\bar{\phi}\) in [SadErh12]).

  • kappa (float) – parameter that constrains the variance of the concentration (referred to as \(\bar{\kappa}\) in [SadErh12])

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e., units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random flip. The list should be as long as the number of sublattices and should sum up to 1.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again.

Example

The following snippet illustrate how to carry out a simple Monte Carlo simulation in the variance-constrained semi-canonical ensemble. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet ECIs are zero and only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

>>> # set up and run MC simulation
>>> structure = prim.repeat(3)
>>> calc = ClusterExpansionCalculator(structure, ce)
>>> phi = 0.6
>>> mc = VCSGCEnsemble(structure=structure, calculator=calc,
...                   temperature=600,
...                   dc_filename='myrun_vcsgc.dc',
...                   phis={'Au': phi},
...                   kappa=200)
>>> mc.run(100)  # carry out 100 trial swaps
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property kappa: float#

Variance constraint parameter \(\bar{\kappa}\) (see parameters section above).

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property phis: dict[int, float]#

Average constraint parameters \(\phi_i\) for all species but one (referred to as \(\bar{\phi}\) in [SadErh12]).

property random_seed: int#

Seed used to initialize the random number generator.

run(number_of_trial_steps)#

Samples the ensemble for the given number of trial steps.

Observations are made whenever the trial step counter is a multiple of the respective interval, so the observation steps stay aligned when run() is called several times. The data container is written to file when the write period has elapsed and at the end of the run, if a file name was given.

Parameters:

number_of_trial_steps (int) – Number of trial steps to carry out.

Raises:

TypeError – If number_of_trial_steps is not an int.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Temperature \(T\) (see parameters section above).

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Hybrid ensemble#

class mchammer.ensembles.HybridEnsemble(structure, calculator, temperature, ensemble_specs, probabilities=None, boltzmann_constant=8.617330337217213e-05, user_tag=None, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, neighbor_sites_to_avoid=None)[source]#

Instances of this class allows one to combine multiple ensembles. In particular, a dictionary should be provided for each ensemble, which must include the type (ensemble) as well as the index of the sublattice (sublattice_index). In addition, it is possible to provide lists of allowed symbols (allowed_symbols) and site indices (allowed_sites) for the trial steps, where the allowed symbols must represent a subset of the elements that can occupy the sites on the specified sublattice. Note that additional arguments are required for the SGC and VCSGC ensembles, namely chemical potentials (chemical_potentials) for the former and constraint parameters (phis and kappa) for the latter. For more detailed information regarding the different ensembles, please see CanonicalEnsemble, SemiGrandCanonicalEnsemble, and VCSGCEnsemble.

This class is particularly useful for effectively sampling complex multi-component systems with several active sublattices, in which case different ensembles can be defined for each of the latter. The fact that it is possible to set the allowed chemical symbols means that one can vary the concentrations of a few selected species, with the help of a VCSGC or semi-grand canonical ensemble, while still allowing swaps between any two sites, using a canonical ensemble (see also the below example). It is advisable to carefully consider how to define the ensemble probabilities. By default the ensembles are weighted by the sizes of the corresponding sublattices, which should give suitable probabilities in most cases. As is shown in the example below, it might be prudent to provide different values if allowed symbols are provided as well as for cases where there are several ensembles that are active on different sublattices.

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • temperature (float) – Temperature \(T\) in appropriate units, commonly Kelvin.

  • ensemble_specs (list[dict]) –

    A list of dictionaries, which should contain the following items:

    • ensemble: Can be either 'vcsgc', 'semi-grand' or 'canonical', lowercase and uppercase letters or any combination thereof are accepted (required).

    • sublattice_index: Index for the sublattice of interest (required).

    • allowed_symbols: List of allowed chemical symbols. By default the symbols allowed on the sublattice by the cluster space.

    • allowed_sites: List of allowed sites. By default all sites of the sublattice.

    • chemical_potentials: Dictionary of chemical potentials for each species \(\mu_i\). The key denotes the species, the value specifies the chemical potential in units that are consistent with the underlying cluster expansion. Only applicable and required for SGC ensembles.

    • phis: Dictionary with average constraint parameters \(\phi_i\). The key denotes the species. For a \(N\)-component sublattice, there should be \(N - 1\) different \(\phi_i\) (referred to as \(\bar{\phi}\) in [SadErh12]). Only applicable and required for VCSGC ensembles.

    • kappa: Parameter that constrains the variance of the concentration (referred to as \(\bar{\kappa}\) in [SadErh12]). Only applicable and required for VCSGC ensembles.

  • probabilities (list[float] | None) – Probabilities for choosing a particular ensemble with the same length as ensemble specs. If left unspecified the probabilties are scaled to match the sizes of the associated sublattices.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e., units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • neighbor_sites_to_avoid (dict[int, list[int]] | None) – Sites that must not be occupied simultaneously, keyed by site index. A site counts as occupied when it is not held by a vacancy, and sites that are absent from the mapping are unconstrained. The mapping must be symmetric and the initial configuration must already satisfy it. The constraint applies to every trial step, whichever ensemble takes it, and is enforced by rejecting trial moves that would violate it, which keeps the acceptance criterion exact. A restrictive mapping can nevertheless leave the allowed configurations disconnected under the available trial moves, in which case only the reachable part of them is sampled. The mapping itself is not written to the data container, only a digest of it, so a restart has to be given the same mapping again. By default no constraint is applied.

Example

The following snippet illustrates how to carry out a simple Monte Carlo simulation using a combination of one canonical and one VCSGC ensemble. Specifically, the concentration of one species (Au) is kept constant while the others (Ag and Pd) are varied, while swaps are still allowed. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet ECIs are zero and only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3],
...                  chemical_symbols=['Ag', 'Au', 'Pd'])
>>> ce = ClusterExpansion(
...    cs, [0, 0, 0, 0.1, 0.1, 0.1, -0.02, -0.02, -0.02])

>>> # define structure object
>>> structure = prim.repeat(3)
>>> for i, atom in enumerate(structure):
>>>    if i % 2 == 0:
>>>        atom.symbol = 'Ag'
>>>    elif i % 3 == 0:
>>>        atom.symbol = 'Pd'

>>> # the default probabilities for this case would be [0.5, 0.5], but
>>> # since the VCSGC ensemble only performs flips on a subset of all
>>> # sites on the sublattice, namely those originally occupied by Ag
>>> # and Pd atoms, specific values will be provided
>>> weights = [len(structure),
...            len([s for s in structure.symbols if s != 'Au'])]
>>> norm = sum(weights)
>>> probabilities = [w / norm for w in weights]

>>> # set up and run MC simulation
>>> calc = ClusterExpansionCalculator(structure, ce)
>>> ensemble_specs = [
...    {'ensemble': 'canonical', 'sublattice_index': 0},
...    {'ensemble': 'vcsgc', 'sublattice_index': 0,
...     'phis': {'Ag': -0.2}, 'kappa': 200,
...     'allowed_symbols':['Ag', 'Pd']}]
>>> mc = HybridEnsemble(structure=structure, calculator=calc,
...                     ensemble_specs=ensemble_specs,
...                     temperature=600, probabilities=probabilities,
...                     dc_filename='myrun_hybrid.dc')
>>> mc.run(100)  # carry out 100 trial steps
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property probabilities: list[float]#

Ensemble propabilities.

property random_seed: int#

Seed used to initialize the random number generator.

run(number_of_trial_steps)#

Samples the ensemble for the given number of trial steps.

Observations are made whenever the trial step counter is a multiple of the respective interval, so the observation steps stay aligned when run() is called several times. The data container is written to file when the write period has elapsed and at the end of the run, if a file name was given.

Parameters:

number_of_trial_steps (int) – Number of trial steps to carry out.

Raises:

TypeError – If number_of_trial_steps is not an int.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Current temperature.

property trial_steps_per_ensemble: dict[str, int]#

Number of Monte Carlo trial steps for each ensemble.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Wang-Landau ensemble#

For analysis functions see here.

class mchammer.ensembles.WangLandauEnsemble(structure, calculator, energy_spacing, energy_limit_left=None, energy_limit_right=None, trial_move='swap', fill_factor_limit=1e-06, flatness_check_interval=None, flatness_limit=0.8, window_search_penalty=2.0, user_tag=None, dc_filename=None, random_seed=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None)[source]#

Instances of this class allow one to sample a system using the Wang-Landau (WL) algorithm, see Phys. Rev. Lett. 86, 2050 (2001) [WanLan01a]. The WL algorithm enables one to acquire the density of states (DOS) as a function of energy, from which one can readily calculate many thermodynamic observables as a function of temperature. To this end, the WL algorithm accumulates both the microcanonical entropy \(S(E)\) and a histogram \(H(E)\) on an energy grid with a predefined spacing (energy_spacing).

The algorithm is initialized as follows.

  1. Generate an initial configuration.

  2. Initialize counters for the microcanonical entropy \(S(E)\) and the histogram \(H(E)\) to zero.

  3. Set the fill factor \(f=1\).

It then proceeds as follows.

  1. Propose a new configuration (see trial_move).

  2. Accept or reject the new configuration with probability

    \[P = \min \{ 1, \, \exp [ S(E_\mathrm{cur}) - S(E_\mathrm{new}) ] \},\]

    where \(E_\mathrm{cur}\) and \(E_\mathrm{new}\) are the energies of the current and new configurations, respectively.

  3. Update the microcanonical entropy \(S(E)\leftarrow S(E) + f\) and histogram \(H(E) \leftarrow H(E) + 1\) where \(E\) is the energy of the system at the end of the move.

  4. Check the flatness of the histogram \(H(E)\). If \(H(E) > \chi \langle H(E)\rangle\,\forall E\) reset the histogram \(H(E) = 0\) and reduce the fill factor \(f \leftarrow f / 2\). The parameter \(\chi\) is set via flatness_limit.

  5. If \(f\) is smaller than fill_factor_limit terminate the loop, otherwise return to 1.

The microcanonical entropy \(S(E)\) and the histogram along with related information are written to the data container every time \(f\) is updated. Using the density \(\rho(E) = \exp S(E)\) one can then readily compute various thermodynamic quantities, including, e.g., the average energy:

\[\left<E\right> = \frac{\sum_E E \rho(E) \exp(-E / k_B T)}{ \sum_E \rho(E) \exp(-E / k_B T)}\]
Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Wang-Landau simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating potential changes.

  • trial_move (str) – One can choose between two different trial moves for generating new configurations. In a 'swap' move two sites are selected and their occupations are swapped. In a 'flip' move one site is selected and its occupation is flipped to a different species. While 'swap' moves conserve the concentrations of the species in the system, 'flip' moves allow one in principle to sample the full composition space. By default 'swap'.

  • energy_spacing (float) – Sets the bin size of the energy grid on which the microcanonical entropy \(S(E)\), and thus the density \(\exp S(E)\), is evaluated. The spacing should be small enough to capture the features of the density of states. Too small values will, however, render the convergence very tedious if not impossible.

  • energy_limit_left (float | None) – Sets the lower limit of the energy range within which the microcanonical entropy \(S(E)\) will be sampled. By default (None) no limit is imposed. Setting limits can be useful if only a part of the density of states is required.

  • energy_limit_right (float | None) – Sets the upper limit of the energy range within which the microcanonical entropy \(S(E)\) will be sampled. By default (None) no limit is imposed. Setting limits can be useful if only a part of the density of states is required.

  • fill_factor_limit (float) – If the fill_factor \(f\) falls below this value, the WL sampling loop is terminated. By default 1e-6.

  • flatness_check_interval (int | None) – For computational efficiency the flatness condition is only evaluated every flatness_check_interval-th trial step. By default (None) flatness_check_interval is set to 1000 times the number of sites in structure, i.e., 1000 Monte Carlo sweeps.

  • flatness_limit (float) – The histogram \(H(E)\) is deemed sufficiently flat if \(H(E) > \chi \left<H(E)\right>\,\forall E\). flatness_limit sets the parameter \(\chi\). By default 0.8.

  • window_search_penalty (float) – If energy_limit_left and/or energy_limit_right have been set, a modified acceptance probability, \(P=\min\{1,\,\exp[C_\mathrm{WSP}(d_\mathrm{new}- d_\mathrm{cur})]\}\), will be used until a configuration is found within the interval of interest. This parameter, specifically, corresponds to \(C_\mathrm{WSP}\), which controls how strongly moves that lead to an increase in the distance, i.e. difference in energy divided by the energy spacing, to the energy window (\(d_\mathrm{new}> d_\mathrm{cur}\)) should be penalized. A higher value leads to a lower acceptance probability for such moves. By default 2.0.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensemble specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random swap. The list must contain as many elements as there are sublattices and it needs to sum up to 1.

Example

The following snippet illustrates how to carry out a Wang-Landau simulation. For the purpose of demonstration, the parameters of the cluster expansion are set to obtain a two-dimensional square Ising model, one of the systems studied in the original work by Wang and Landau:

>>> from ase import Atoms
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator
>>> from mchammer.ensembles import WangLandauEnsemble

>>> # prepare cluster expansion
>>> prim = Atoms('Au', positions=[[0, 0, 0]], cell=[1, 1, 10], pbc=True)
>>> cs = ClusterSpace(prim, cutoffs=[1.1], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 2])

>>> # prepare initial configuration
>>> structure = prim.repeat((4, 4, 1))
>>> for k in range(8):
...     structure[k].symbol = 'Ag'

>>> # set up and run Wang-Landau simulation
>>> calculator = ClusterExpansionCalculator(structure, ce)
>>> mc = WangLandauEnsemble(structure=structure,
...                         calculator=calculator,
...                         energy_spacing=1,
...                         dc_filename='ising_2d_run.dc')
>>> # N.B.: in practice one requires many more steps
>>> mc.run(number_of_trial_steps=len(structure) * 100)
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property converged: bool | None#

Whether convergence has been achieved.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

property fill_factor: float#

Current value of the fill factor.

property fill_factor_history: dict[int, float]#

evolution of the fill factor in the Wang-Landau algorithm (key=MC trial step, value=fill factor)

property fill_factor_limit: float#

If the fill factor \(f\) falls below this value, the Wang-Landau sampling is terminated.

property flatness_check_interval: int#

Number of MC trial steps between checking the flatness condition.

property flatness_limit: float#

The histogram \(H(E)\) is deemed sufficiently flat if \(H(E) > \chi \left<H(E)\right>\,\forall E\) where flatness_limit sets the parameter \(\chi\).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run(number_of_trial_steps)[source]#

Samples the ensemble for the given number of trial steps.

Parameters:

number_of_trial_steps (int) – Maximum number of MC trial steps to run in total. The run will terminate earlier if fill_factor_limit is reached.

Raises:

TypeError – If number_of_trial_steps is not an int.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)[source]#

Updates the last state of the Wang-Landau simulation and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

mchammer.ensembles.wang_landau_ensemble.get_bins_for_parallel_simulations(n_bins, energy_spacing, minimum_energy, maximum_energy, overlap=4, bin_size_exponent=1.0)[source]#

Generates a list of energy bins (lower and upper bound) suitable for parallel Wang-Landau simulations. For the latter, the energy range is split up into a several bins (n_bins). Each bin is then sampled in a separate Wang-Landau simulation. Once the density of states in the individual bins has been converged the total density of states can be constructed by patching the segments back together. To this end, one requires some over overlap between the segments (overlap).

The function returns a list of tuples. Each tuple provides the lower (energy_limit_left) and upper (energy_limit_right) bound of one bin, which are then to be used to set energy_limit_left and energy_limit_right when initializing a WangLandauEnsemble instance.

Note

The left-most/right-most bin has no lower/upper bound (set to None).

Parameters:

  • n_bins (int) – Number of bins.

  • energy_spacing (float) – Sets the bin size of the energy grid used by the Wang-Landau simulation, see WangLandauEnsemble for details.

  • minimum_energy (float) – An estimate for the lowest energy to be encountered in this system.

  • maximum_energy (float) – An estimate for the highest energy to be encountered in this system.

  • overlap (int) – Amount of overlap between bins in units of energy_spacing. By default 4.

  • bin_size_exponent (float) –

    This parameter allows one to generate a non-uniform distribution of bin sizes. If bin_size_exponent is smaller than one bins at the lower and upper end of the energy range (specified via minimum_energy and maximum_energy) will be shrunk relative to the bins in the middle of the energy range. In principle this can be used one to achieve a more even distribution of computational load between the individual Wang-Landau simulations. By default 1.0, which gives every bin the same size.

    Note

    This is an option for advanced users. Only use this keyword if you know what you are doing.

Return type:

list[tuple[float, float]]

Thermodynamic-integration ensemble#

For analysis functions see here.

class mchammer.ensembles.ThermodynamicIntegrationEnsemble(structure, calculator, temperature, n_steps, forward, user_tag=None, boltzmann_constant=8.617330337217213e-05, random_seed=None, dc_filename=None, data_container_write_period=600, ensemble_data_write_interval=None, trajectory_write_interval=None, sublattice_probabilities=None)[source]#

Instances of this class allow one to find the free energy of the system. To this end, we use the canonical ensemble with a modified Hamiltonian,

\[H(\lambda) = (1 - \lambda) H_{A} + \lambda H_{B}\]

The Hamiltonian is then sampled continuously from \(\lambda=0\) to \(\lambda=1\). \(H_{B}\) is your cluster expansion and \(H_{A}=0\), is a completely disordered system, with free energy given by the ideal mixing entropy.

The free energy, A, of system B is then given by:

\[A_{B} = A_{A} + \int_{0}^{1} \left\langle\frac{\mathrm{d}H(\lambda)} {\mathrm{d}\lambda}\right\rangle_{H} \mathrm{d}\lambda\]

and since \(A_{A}\) is known it is easy to compute \(A_{B}\)

\(\lambda\) is parametrized as,

\[\lambda(x) = x^5(70x^4 - 315x^3 + 540x^2 - 420x + 126)\]

where \(x = (\mathrm{step} + 1) / \mathrm{n\_steps}\).

Parameters:

  • structure (Atoms) – Atomic configuration to be used in the Monte Carlo simulation. It also defines the initial occupation vector.

  • calculator (BaseCalculator) – Calculator to be used for calculating the potential changes that enter the evaluation of the Metropolis criterion.

  • temperature (float) – Temperature \(T\) in appropriate units, commonly Kelvin.

  • n_steps (int) – Number of MC trial steps over which \(\lambda\) is varied from 0 to 1.

  • forward (bool) – If this is set to True the simulation runs from \(H_A\) to \(H_B\), otherwise it runs from \(H_B\) to \(H_A\). \(H_B\) is the cluster expansion and \(H_A = 0\), is the fully disordered system.

  • boltzmann_constant (float) – Boltzmann constant \(k_B\) in appropriate units, i.e., units that are consistent with the underlying cluster expansion and the temperature units. By default eV/K.

  • user_tag (str | None) – Human-readable tag for the ensemble.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

  • dc_filename (str | None) – Name of file the data container associated with the ensemble will be written to. If the file exists it will be read, the data container will be appended, and the file will be updated/overwritten.

  • data_container_write_period (float) – Period in seconds at which the data container is written to file. Writing periodically to file provides both a way to examine the progress of the simulation and to back up the data. By default 600 s.

  • ensemble_data_write_interval (int | None) – Interval at which data is written to the data container. This includes for example the current value of the calculator (i.e., usually the energy) as well as ensembles specific fields such as temperature or the number of atoms of different species. By default the number of sites in structure.

  • trajectory_write_interval (int | None) – Interval at which the current occupation vector of the atomic configuration is written to the data container. By default the number of sites in structure.

  • sublattice_probabilities (list[float] | None) – Probability for picking a sublattice when doing a random swap. This should be as long as the number of sublattices and should sum up to 1.

Example

The following snippet illustrate how to carry out a simple thermodynamic integration. Here, the parameters of the cluster expansion are set to emulate a simple Ising model in order to obtain an example that can be run without modification. In practice, one should of course use a proper cluster expansion:

>>> from ase.build import bulk
>>> from icet import ClusterExpansion, ClusterSpace
>>> from mchammer.calculators import ClusterExpansionCalculator

>>> # prepare cluster expansion
>>> # the setup emulates a second nearest-neighbor (NN) Ising model
>>> # (zerolet and singlet ECIs are zero and only first and second neighbor
>>> # pairs are included)
>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
>>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])

>>> # prepare initial configuration
>>> structure = prim.repeat(3)
>>> for k in range(5):
>>>     structure[k].symbol = 'Ag'

>>> # set up and run MC simulation
>>> calc = ClusterExpansionCalculator(structure, ce)
>>> mc = ThermodynamicIntegrationEnsemble(structure=structure, calculator=calc,
...                                       temperature=600,
...                                       n_steps=100000,
...                                       forward=True,
...                                       dc_filename='myrun_thermodynamic_integration.dc')
>>> mc.run()
attach_observer(observer, tag=None)#

Attaches an observer to the ensemble.

If the observer has no observation interval, it is set to the number of sites in the structure.

Parameters:

  • observer (BaseObserver) – Observer to attach.

  • tag (str | None) – Tag under which the observer and, for a scalar observer, its data are stored. By default the tag of the observer.

Raises:

TypeError – If observer is not a BaseObserver.

Return type:

None

property boltzmann_constant: float#

Boltzmann constant \(k_B\) in the units of the cluster expansion.

property calculator: BaseCalculator#

Calculator attached to the ensemble.

property data_container: BaseDataContainer#

Data container associated with the ensemble.

do_canonical_swap(sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_sgc_flip(chemical_potentials, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • chemical_potentials (dict[int, float]) – Chemical potentials used to calculate the potential difference.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_thermodynamic_swap(sublattice_index, lambda_val, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • sublattice_index (int) – The sublattice the swap will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

do_vcsgc_flip(phis, kappa, sublattice_index, allowed_species=None, allowed_sites=None)#

Carries out one Monte Carlo trial step.

Parameters:

  • phis (dict[int, float]) – Average constraint parameters.

  • kappa (float) – Parameter that constrains the variance of the concentration.

  • sublattice_index (int) – The sublattice the flip will act on.

  • allowed_species (list[int] | None) – List of atomic numbers for allowed species.

  • allowed_sites (list[int] | None) – List of indices for allowed sites.

Return type:

int

Returns:

Returns 1 or 0 depending on if trial move was accepted or rejected.

property ensemble_parameters: dict#

Parameters that identify the ensemble (copy).

get_random_sublattice_index(probability_distribution)#

Returns the index of a sublattice drawn at random with the given probabilities.

Parameters:

probability_distribution (list[float]) – Probability of each sublattice, in the order of sublattices.

Raises:

ValueError – If the number of probabilities differs from the number of sublattices.

Return type:

int

property n_steps: int#

Number of trial steps over which the coupling parameter is varied.

property observer_interval: int#

Number of trial steps between observations, which is the greatest common divisor of the intervals of the observers and of the ensemble data and trajectory writing.

property observers: dict[str, BaseObserver]#

Attached observers, keyed by tag.

property random_seed: int#

Seed used to initialize the random number generator.

run()[source]#

Runs the thermodynamic integration.

Return type:

None

property step: int#

Number of trial steps carried out so far, including restarts.

property structure: Atoms#

Current configuration (copy).

property sublattices: Sublattices#

Sublattices of the configuration being sampled.

property temperature: float#

Current temperature.

update_occupations(sites, species)#

Changes the occupations of the given sites in the configuration being sampled.

The change is applied to the configuration manager and to the calculator, and either both advance or neither does.

Parameters:

  • sites (list[int]) – Indices of the sites to change.

  • species (list[int]) – New occupations of those sites by atomic number.

Raises:

ValueError – If the two lists differ in length.

Return type:

None

property user_tag: str | None#

Human-readable tag for the ensemble.

write_data_container(outfile)#

Saves the last state of the simulation in the data container and writes the data container to file.

Parameters:

outfile (str | bytes) – File to which to write.

Return type:

None

Target cluster vector annealing#

class mchammer.ensembles.TargetClusterVectorAnnealing(structure, calculators, T_start=5.0, T_stop=0.001, random_seed=None)[source]#

Instances of this class carry out simulated annealing towards a target cluster vector.

Each trial step swaps two sites in one of the supercells and accepts the swap with the Metropolis criterion applied to the objective function of the TargetVectorCalculator of that supercell. The artificial temperature decreases exponentially from T_start to T_stop over the course of the run, and the structure with the best score seen during the run is kept. Since it is impossible to know beforehand which supercell shape accommodates the best match, several supercells can be annealed at the same time.

The functions described in the section on special quasirandom structures, such as generate_sqs and generate_target_structure, set up the calculators and this class for the common cases.

Parameters:

  • structure (list[Atoms]) – Supercells to anneal. Their occupations define the initial configurations.

  • calculators (list[TargetVectorCalculator]) – One calculator per supercell in structure, in the same order.

  • T_start (float) – Artificial temperature at which the annealing starts.

  • T_stop (float) – Artificial temperature at which the annealing stops.

  • random_seed (int | None) – Seed for the random number generator used in the Monte Carlo simulation.

Raises:

ValueError – If structure is a single structure instead of a list, or if the number of structures differs from the number of calculators.

Example

The following snippet generates a special quasirandom structure for a binary alloy at equal concentrations. The cluster vector of the ideal random alloy has all elements beyond the zerolet equal to zero, which is the target:

>>> from ase.build import bulk
>>> from icet import ClusterSpace
>>> from icet.tools.structure_generation import occupy_structure_randomly
>>> from mchammer.calculators import TargetVectorCalculator
>>> from mchammer.ensembles import TargetClusterVectorAnnealing

>>> prim = bulk('Au')
>>> cs = ClusterSpace(prim, cutoffs=[6.0], chemical_symbols=['Ag', 'Au'])
>>> target_vector = [1.0] + [0.0] * (len(cs) - 1)

>>> # two supercell shapes with random starting occupations
>>> supercells = [prim.repeat((2, 2, 2)), prim.repeat((4, 2, 1))]
>>> for supercell in supercells:
...     occupy_structure_randomly(supercell, cs, {'Ag': 0.5, 'Au': 0.5})
>>> calculators = [TargetVectorCalculator(supercell, cs, target_vector)
...                for supercell in supercells]

>>> annealing = TargetClusterVectorAnnealing(supercells, calculators,
...                                          T_start=5.0, T_stop=0.001)
>>> sqs = annealing.generate_structure(number_of_trial_steps=2000)
>>> print(annealing.best_score)
property T_start: float#

Artificial temperature at which the annealing starts.

property T_stop: float#

Artificial temperature at which the annealing stops.

property accepted_trials: int#

Number of trial steps accepted so far.

property best_score: float#

Best score found so far.

property best_structure: Atoms#

Structure with the best score found so far.

property current_score: float#

Score of the current configuration.

generate_structure(number_of_trial_steps=None)[source]#

Runs the annealing and returns the best structure found.

Parameters:

number_of_trial_steps (int | None) – Total number of trial steps over all supercells. By default 3000 steps per supercell.

Returns:

The structure with the best score seen during the run.

Return type:

Atoms

property n_steps: int#

Total number of trial steps of the run.

property temperature: float#

Current artificial temperature.

property total_trials: int#

Number of trial steps carried out so far.

On this page