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
atoms = make_supercell(ce.cluster_space.primitive_structure,
                       3*array([[-1, 1, 1],
                                [1, -1, 1],
                                [1, 1, -1]]))
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]:

    # Initialize MC ensemble
    # TODO: remove chemical_potentials once possible
    mc = SemiGrandCanonicalEnsemble(
        calculator=calculator,
        atoms=atoms,
        ensemble_data_write_interval=10,
        temperature=temperature,
        chemical_potentials={chemical_symbols[0]: 0,
                             chemical_symbols[1]: 0})

    # Evolve configuration through the entire composition range
    for dmu in arange(-0.6, 0.51, 0.05):
        mc.chemical_potentials = {chemical_symbols[0]: 0,
                                  chemical_symbols[1]: dmu}

        mc.reset_data_container()
        mc.run(number_of_trial_steps=len(atoms)*30)
        # TODO: change the next line (and the tutorial) once mc.data_container
        # is writable
        mc.data_container.write('sgc-T{}-dmu{:.3f}.dc'
                                .format(temperature, dmu))

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
atoms = make_supercell(ce.cluster_space.primitive_structure,
                       3*array([[-1, 1, 1],
                                [1, -1, 1],
                                [1, 1, -1]]))
calculator = ClusterExpansionCalculator(atoms, ce)

# step 2: carry out Monte Carlo simulations
for temperature in [900, 600, 300]:

    # Initialize MC ensemble
    # TODO: remove chemical_potentials once possible
    mc = SemiGrandCanonicalEnsemble(
        calculator=calculator,
        atoms=atoms,
        ensemble_data_write_interval=10,
        temperature=temperature,
        chemical_potentials={chemical_symbols[0]: 0,
                             chemical_symbols[1]: 0})

    # Evolve configuration through the entire composition range
    for dmu in arange(-0.6, 0.51, 0.05):
        mc.chemical_potentials = {chemical_symbols[0]: 0,
                                  chemical_symbols[1]: dmu}

        mc.reset_data_container()
        mc.run(number_of_trial_steps=len(atoms)*30)
        # TODO: change the next line (and the tutorial) once mc.data_container
        # is writable
        mc.data_container.write('sgc-T{}-dmu{:.3f}.dc'
                                .format(temperature, dmu))