Monte Carlo simulations

We are now in a position to carry out a series of Monte Carlo (MC) simulations to sample the cluster expansion model that was constructed and validated in the previous steps. To set up the simulation we first construct a supercell and initialize an associated calculator by combining our cluster expansion model with the supercell.

ce = ClusterExpansion.read('mixing_energy.ce')
chemical_symbols = ce.cluster_space.chemical_symbols[0]
atoms = make_supercell(ce.cluster_space.primitive_structure,
                       3 * array([[-1, 1, 1],
                                  [1, -1, 1],
                                  [1, 1, -1]]))
# TODO: Remove this line once atoms is not longer decorated with H atoms
atoms.numbers = [47]*len(atoms)

calculator = ClusterExpansionCalculator(atoms, ce)

In this example the sampling will be carried out in the semi-grand canonical (SGC) ensemble. To this end, we set up a SGC ensemble object object and loop over both temperatures and chemical potential differences.

We carry out a rather long MC run, anticipating that the analysis will only include the latter part of the simulation after equilibration. After the run the results are written to file (in the form of a DataContainer object). The latter will be used in the next step to analyze the runs. Note that the ensemble object is only initialized once for each temperature. Thereby the configuration evolves gradually and the period needed for equilibration is shortened.

for temperature in [900, 600, 300]:
    # Evolve configuration through the entire composition range
    for dmu in arange(-0.6, 0.51, 0.05):
        # Initialize MC ensemble
        mc = SemiGrandCanonicalEnsemble(
            atoms=atoms,
            calculator=calculator,
            temperature=temperature,
            data_container='sgc-T{}-dmu{:.3f}.dc'.format(temperature, dmu),
            chemical_potentials={chemical_symbols[0]: 0,
                                 chemical_symbols[1]: dmu})

        mc.run(number_of_trial_steps=len(atoms) * 30)

On an Intel i7-950 CPU the set up of the calculator takes about 10 seconds, whereas the Monte Carlo simulation takes about 0.2 seconds per MC trial step.

Source code

The complete source code is available in tutorial/basic/5_run_monte_carlo.py
from ase.build import make_supercell
from numpy import arange, array

from icet import ClusterExpansion
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import SemiGrandCanonicalEnsemble

# step 1: set up the structure to simulate and the calculator
ce = ClusterExpansion.read('mixing_energy.ce')
chemical_symbols = ce.cluster_space.chemical_symbols[0]
atoms = make_supercell(ce.cluster_space.primitive_structure,
                       3 * array([[-1, 1, 1],
                                  [1, -1, 1],
                                  [1, 1, -1]]))
# TODO: Remove this line once atoms is not longer decorated with H atoms
atoms.numbers = [47]*len(atoms)

calculator = ClusterExpansionCalculator(atoms, ce)

# step 2: carry out Monte Carlo simulations
for temperature in [900, 600, 300]:
    # Evolve configuration through the entire composition range
    for dmu in arange(-0.6, 0.51, 0.05):
        # Initialize MC ensemble
        mc = SemiGrandCanonicalEnsemble(
            atoms=atoms,
            calculator=calculator,
            temperature=temperature,
            data_container='sgc-T{}-dmu{:.3f}.dc'.format(temperature, dmu),
            chemical_potentials={chemical_symbols[0]: 0,
                                 chemical_symbols[1]: dmu})

        mc.run(number_of_trial_steps=len(atoms) * 30)