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')
structure = make_supercell(ce.get_cluster_space_copy().primitive_structure,
                           3 * np.array([[-1, 1, 1],
                                         [1, -1, 1],
                                         [1, 1, -1]]))
calculator = ClusterExpansionCalculator(structure, 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 stored on disk in the form of a DataContainer object. The latter will be used in the next step to analyze the runs. At the end of each iteration we save the last state of the system and provide it as input to the ensemble object in the next iteration. Thereby the configuration evolves gradually and the period needed for equilibration is shortened.

# Make sure output directory exists
output_directory = 'monte_carlo_data'
try:
    mkdir(output_directory)
except FileExistsError:
    pass
for temperature in [900, 300]:
    # Evolve configuration through the entire composition range
    for dmu in np.arange(-0.7, 0.51, 0.05):
        # Initialize MC ensemble
        mc = SemiGrandCanonicalEnsemble(
            structure=structure,
            calculator=calculator,
            temperature=temperature,
            dc_filename='{}/sgc-T{}-dmu{:+.3f}.dc'.format(output_directory, temperature, dmu),
            chemical_potentials={'Ag': 0, 'Pd': dmu})

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

On an Intel i5-6400 CPU the set up of the calculator takes about 1 second, whereas the Monte Carlo simulation takes about 1 to 2 milliseconds per MC trial step, but note that for a production simulation, the supercell should be larger and the sampling both longer and more dense.

The VCSGC ensemble¶

A simulation in the VCSGC ensemble can be carried on an analogous fashion.

# Make sure output directory exists
output_directory = 'monte_carlo_data'
try:
    mkdir(output_directory)
except FileExistsError:
    pass
for temperature in [900, 300]:
    # Evolve configuration through the entire composition range
    for phi in np.arange(-2.1, 0.11, 0.08):
        # Initialize MC ensemble
        mc = VCSGCEnsemble(
            structure=structure,
            calculator=calculator,
            temperature=temperature,
            dc_filename='{}/vcsgc-T{}-phi{:+.3f}.dc'.format(output_directory, temperature, phi),
            phis={'Pd': phi},
            kappa=200)

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

Here, we chose \(\kappa=200\), which usually provides a good compromise between acceptance probability and fluctuations in variance.

Source code¶

from ase.build import make_supercell
from icet import ClusterExpansion
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import SemiGrandCanonicalEnsemble
import numpy as np
from os import mkdir

# step 1: Set up structure to simulate as well as calculator
ce = ClusterExpansion.read('mixing_energy.ce')
structure = make_supercell(ce.get_cluster_space_copy().primitive_structure,
                           3 * np.array([[-1, 1, 1],
                                         [1, -1, 1],
                                         [1, 1, -1]]))
calculator = ClusterExpansionCalculator(structure, ce)

# step 2: Carry out Monte Carlo simulations
# Make sure output directory exists
output_directory = 'monte_carlo_data'
try:
    mkdir(output_directory)
except FileExistsError:
    pass
for temperature in [900, 300]:
    # Evolve configuration through the entire composition range
    for dmu in np.arange(-0.7, 0.51, 0.05):
        # Initialize MC ensemble
        mc = SemiGrandCanonicalEnsemble(
            structure=structure,
            calculator=calculator,
            temperature=temperature,
            dc_filename='{}/sgc-T{}-dmu{:+.3f}.dc'.format(output_directory, temperature, dmu),
            chemical_potentials={'Ag': 0, 'Pd': dmu})

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

from ase.build import make_supercell
from icet import ClusterExpansion
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import VCSGCEnsemble
import numpy as np
from os import mkdir

# step 1: Set up structure to simulate as well as calculator
ce = ClusterExpansion.read('mixing_energy.ce')
structure = make_supercell(ce.get_cluster_space_copy().primitive_structure,
                           3 * np.array([[-1, 1, 1],
                                         [1, -1, 1],
                                         [1, 1, -1]]))
calculator = ClusterExpansionCalculator(structure, ce)

# step 2: Carry out Monte Carlo simulations
# Make sure output directory exists
output_directory = 'monte_carlo_data'
try:
    mkdir(output_directory)
except FileExistsError:
    pass
for temperature in [900, 300]:
    # Evolve configuration through the entire composition range
    for phi in np.arange(-2.1, 0.11, 0.08):
        # Initialize MC ensemble
        mc = VCSGCEnsemble(
            structure=structure,
            calculator=calculator,
            temperature=temperature,
            dc_filename='{}/vcsgc-T{}-phi{:+.3f}.dc'.format(output_directory, temperature, phi),
            phis={'Pd': phi},
            kappa=200)

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