# 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
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))