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 construct a supercell and then create a calculator by combining our cluster expansion model with the supercell.

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

# step 1: set up the structure to simulate
ce = ClusterExpansion.read('mixing_energy.ce')
chemical_symbols = ce.cluster_space.chemical_symbols
atoms = ce.cluster_space.primitive_structure.repeat(3)
calculator = ClusterExpansionCalculator(atoms, ce)

In this example the sampling will be carried out in the semi-grand canonical (SGC) ensemble, which requires a temperature and a set of chemical potentials as input. Accordingly we set up a loop over different temperature and chemical potential values. In the body of the two nested loops we instantiate a SGC ensemble object using the calculator configured before and then run a MC simulation for a number of trial steps.

ntrials = 500
equil = 200
for temperature in [600, 900, 300]:
    for dmu in arange(-2.0, 2.0, 0.1):
        mc = SemiGrandCanonicalEnsemble(
            calculator=calculator, atoms=atoms,
            data_container='sgc.dc',
            random_seed=42, temperature=temperature,
            chemical_potentials={chemical_symbols[0]: 0,
                                 chemical_symbols[1]: dmu},
            ensemble_data_write_interval=10)
        mc.run(ntrials)

Here, the results of the simulation are stored in a data container, which is written to a file named sgc.dc. In the next step these data will be analyzed to generate e.g., a map of the chemical potential difference vs composition.

Source code

The complete source code is available in tutorial/basic/5_run_monte_carlo.py
from icet import ClusterExpansion
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import SemiGrandCanonicalEnsemble
from numpy import arange

# step 1: set up the structure to simulate
ce = ClusterExpansion.read('mixing_energy.ce')
chemical_symbols = ce.cluster_space.chemical_symbols
atoms = ce.cluster_space.primitive_structure.repeat(3)
calculator = ClusterExpansionCalculator(atoms, ce)

# step 2: set up calculator and ensemble
ntrials = 500
equil = 200
for temperature in [600, 900, 300]:
    for dmu in arange(-2.0, 2.0, 0.1):
        mc = SemiGrandCanonicalEnsemble(
            calculator=calculator, atoms=atoms,
            data_container='sgc.dc',
            random_seed=42, temperature=temperature,
            chemical_potentials={chemical_symbols[0]: 0,
                                 chemical_symbols[1]: dmu},
            ensemble_data_write_interval=10)
        mc.run(ntrials)