SemiGrandCanonicalTransitionMatrixHook

class SemiGrandCanonicalTransitionMatrixHook(elements, chemical_potential=None, temperature=None, swap_ratio=None, random_seed=None, call_interval=None, pnpd_update_interval=None, measurement_interval=None)

A hook function that implements semi-grand canonical Monte Carlo simulation

Parameters:
  • elements (Sequence of PeriodicTableElement) – The elements being swapped between.

  • chemical_potential (PhysicalQuantity of type energy) – The chemical potential of swapping from the first to the second atom.
    Default: 0 eV

  • temperature (PhysicalQuantity of type temperature) – The Monte Carlo selection temperature
    Default: 300 Kelvin

  • swap_ratio (float) – The fraction of available atoms that are attempted to be swapped at each call.
    Default: 0.01

  • random_seed (int) – The seed for the random number generator used to select the Monte Carlo moves.

  • call_interval (int) – The number of dynamics steps between each Monte Carlo step.
    Default: 1

  • pnpd_update_interval (int) – The number of Monte Carlo steps between each update of the particle number probability distribution. Note that this is based on Monte Carlo steps, not dynamics steps.
    Default: 1

  • measurement_interval (int) – The number of dynamics steps between recording the particle number probability distribution.
    Default: 1

callInterval()
Returns:

The call interval of this hook function.

Return type:

int

chemicalPotential()
Returns:

The chemical potential difference used in calculating swap probabilities.

Return type:

PhysicalQuantity of type energy

collectionMatrix()
Returns:

The tri-banded collection in which sampled transition probabilities are collected.

Return type:

ndarray

elements()
Returns:

List of the elements being swapped in the Monte Carlo step.

Return type:

list of PeriodicTableElement

measurementInterval()
Returns:

The number of dynamics steps between writing data to the trajectory.

Return type:

int

monteCarloQuantity()

The value being averaged in the Monte Carlo simulation. In this simulation this is the fraction of the second atom in the configuration.

Returns:

The ratio of second atoms to the first in the configuration.

Return type:

float

particleNumberProbabilityDistribution()
Returns:

The relative probability of being in each composition state.

Return type:

ndarray

pnpdUpdateInterval()
Returns:

The number of dynamics steps between measurements of the PNPD.

Return type:

int

randomSeed()
Returns:

The random seed used to create the Monte Carlo random number generator.

Return type:

int

sampleCount()
Returns:

The number of transition samples from each state.

Return type:

ndarray

swapElementCount(configuration)

Count the number of atoms that have been swapped to the second element.

Parameters:

configuration (BulkConfiguration) – The current configuration.

Returns:

The number of elements

Return type:

int

swapRatio()
Returns:

The fraction of atoms to attempt to swap in each Monte Carlo cycle.

Return type:

float

temperature()
Returns:

The temperature of the simulation.

Return type:

PhysicalQuantity of type temperature

uniqueString()

Return a unique string representing the state of the object.

In this example the semi-grand canonical transition matrix Monte Carlo method is used to study the composition of a gold-silver alloy. The result of the calculation is a ConfigurationDataContainer that contains the particle number probability distribution (PNPD) for the system. Since these metals are miscible, the PNPD should have a single maximum at the average composition value of around 25% gold and 75% silver.

A workflow for the calculation can be downloaded here with a corresponding python script Au_Ag_Transition_Matrix.py. This simulation may take around several hours to run, given the number of samples required.

Notes

The SemiGrandCanonicalTransitionMatrixHook defines a transition matrix Monte Carlo simulation in the semi-grand canonical ensemble. This is a complimentary technique to the SemiGrandCanonicalMonteCarloHook method. Both simulations are carried out in the semi-grand canonical ensemble and can be used to determine the composition of materials at different chemical potentials. While in normal Monte Carlo methods the property of interest is directly sampled, in transition matrix Monte Carlo the probability of transitioning from one state to another is sampled[1]. Here this is the number of each element in the material. As only one atom is swapped at a time, the probability of being in any state can be calculated. This is known as the particle number probability density (PNPD). In order to effectively calculate the PNPD it is necessary to sample all states. This is done by biasing the selection by the negative change in the PNPD. Since the PNPD is needed to effectively sample the states, this is calculated and updated during the simulation. Once the converges the simulation samples evenly across all states.

The PNPD relates straightforwardly to the relative free energy of each state such that:

\[W(N; \mu, V , T) = −k_BT \ln \Pi(N; \mu, V , T)\]

Here \(\Pi(N; \mu, P, T)\) is the PNPD at given particle number \(N\) and chemical potential, temperature and pressure and \(W\) is similarly the free energy. To get the average system property from the PNPD, that property is averaged over the available states. For instance, the average number of elements in the material \(\left<N\right>\) can be given as:

\[\left<N\right> = \frac{\sum N \Pi(N; \mu, P, T)}{\sum \Pi(N; \mu, P, T)}\]

The PNPD can also be re-weighted to different chemical potential values. This allows exploring the whole range of chemical potentials with one simulation. The PNPD at a different chemical potential can be given as:

\[\ln \Pi(N; \mu, P, T) = \ln \Pi(N; \mu_0, P, T) + \frac{N (\mu - \mu_0)}{k_BT} + C\]

Here \(\mu_0\) is the chemical potential that the simulations was originally carried out at, \(\mu\) is the new chemical potential and \(C\) is a normalization constant.

One of the main advantages of transition matrix simulations is that it samples the entire space of compositions, which works well with systems that have multiple stable states. In the case of immiscible metals, there can be a minimum in the free energy close to each pure state, and then a large free energy barrier as the elements mix. As the chemical potential changes the relative energies of the minima change, leading to one being more populated than the other. This situation is very difficult to model with traditional Monte Carlo techniques, but can be explored with transition Matrix techniques. The trade-off here is that transition Matrix Monte Carlo methods require at least an order of magnitude more steps than a traditional Monte Carlo simulation to reach convergence. This is due to the transition Matrix simulation scanning the entire space of configurations rather than just a small amount around a minimum in the free energy.

To practically carry out a transition Matrix Monte Carlo simulation, it is recommended that the chemical potential should be near the difference in free energy between the two materials. The PNPD stored in the trajectory is the natural log of the PNPD, as this is numerically much easier to handle. In addition to the arguments to the arguments from the SemiGrandCanonicalMonteCarloHook, the SemiGrandCanonicalTransitionMatrixHook has two additional arguments that control the frequency of updating and storing the PNPD. Frequent updating of the PNPD can slow the simulation in some cases. This is controlled by the argument pnpd_update_interval, which gives the number of Monte Carlo cycles between PNPD updates. Similarly measurement_interval gives the frequency with which the PNPD is written to the trajectory. This should be done often enough so that convergence in the PNPD can be observed, but not so much as to take large amounts of space in the trajectory. This frequency is also based on the number of dynamics steps between recording each PNPD. As well as the PNPD the trajectory also records the number of samples from each state. This is stored as the measurement sample_count. For a converged PNPD each state should be sampled enough times to have an accurate probability of transitioning from that state.