Hybrid ensembles

In systems with multiple sublattices, it can sometimes be advantageous to use different thermodynamical ensembles for different sublattices in a single Monte Carlo simulation. This can for example be relevant in systems in which one sublattice does not exchange atoms with the environment (closed system) while the other does (open system). Metallic alloys exposed to hydrogen are examples where this is often the case (for an example, see [RahLofFra21]). This tutorial demonstrates how sublattice-specific ensembles can be used in mchammer using the HybridEnsemble class.

As in any Monte Carlo simulation with mchammer, the first steps are to define a simulation cell and to construct a ClusterExpansionCalculator object. To this end, we first construct a toy ClusterExpansion for a system with two sublattices, one occupied by Pd/Au and one occupied by H and vacancies.

prim = bulk('Pd', a=4.0)
prim.append(Atom('H', position=(2, 2, 2)))
cs = ClusterSpace(prim, cutoffs=[3], chemical_symbols=[('Au', 'Pd'), ('H', 'X')])
ce = ClusterExpansion(cluster_space=cs, parameters=[-0.15, 0, 0, 0, 0.1, 0.05])
structure = prim.repeat(5)
calculator = ClusterExpansionCalculator(structure, ce)

We then define the parameters that will enter our Monte Carlo simulation. Here, we will run a simulation in which the Pd-Au sublattice is sampled in the canonical ensemble and the H-vacancy sublattice in the semi-grand canonical ensemble. The ensembles are specified via a list of dictionaries, which define the parameters specific to each ensemble. Moreover, since the concentrations on the Pd-Au sublattice are fixed once a starting configuration is defined, we must also create a supercell with the desired concentration (note that the concentrations on the H-vacancy sublattice will change during the simulation since concentrations are not conserved in the semi-grand canonical ensemble, hence the choice of starting concentrations on the H-vacancy sublattice is unimportant).

# Make sure output directory exists
output_directory = 'monte_carlo_data'
try:
    mkdir(output_directory)
except FileExistsError:
    pass

muH = -0.1
temp = 300
cAu = 0.2
cH_start = 0.2
ensemble_specs = [{'ensemble': 'canonical', 'sublattice_index': 0},
                  {'ensemble': 'semi-grand', 'sublattice_index': 1,
                   'chemical_potentials': {'H': muH, 'X': 0}}]
occupy_structure_randomly(structure=structure,
                          cluster_space=cs,
                          target_concentrations={'A': {'Pd': 1 - cAu, 'Au': cAu},
                                                 'B': {'H': cH_start, 'X': 1 - cH_start}})

Finally, we define our HybridEnsemble and commence the simulation.

mc = HybridEnsemble(
    structure=structure,
    calculator=calculator,
    temperature=300,
    ensemble_specs=ensemble_specs,
    dc_filename=f'{output_directory}/hybrid-T{temp}-muH{muH:+.3f}-cAu{cAu:.3f}.dc')

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

The HybridEnsemble can also be used together with the variance-constrained semi-grand canonical (VCGSC) ensemble. The following is a valid specification for a simulation with the canonical ensemble on the Pd/Au sublattice and VCSGC on the H/vacancy sublattice:

ensemble_specs = [{'ensemble': 'canonical', 'sublattice_index': 0},
                  {'ensemble': 'vcsgc', 'sublattice_index': 1,
                   'phis': {'H': phiH}, 'kappa': 200}]

Source code

The complete source code is available in examples/advanced_topics/hybrid_ensembles.py

from ase import Atom
from ase.build import bulk
from icet import ClusterSpace, ClusterExpansion
from icet.tools.structure_generation import occupy_structure_randomly
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import HybridEnsemble
from os import mkdir

# step 1: Set up cluster expansion, structure and calculator
prim = bulk('Pd', a=4.0)
prim.append(Atom('H', position=(2, 2, 2)))
cs = ClusterSpace(prim, cutoffs=[3], chemical_symbols=[('Au', 'Pd'), ('H', 'X')])
ce = ClusterExpansion(cluster_space=cs, parameters=[-0.15, 0, 0, 0, 0.1, 0.05])
structure = prim.repeat(5)
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

muH = -0.1
temp = 300
cAu = 0.2
cH_start = 0.2
ensemble_specs = [{'ensemble': 'canonical', 'sublattice_index': 0},
                  {'ensemble': 'semi-grand', 'sublattice_index': 1,
                   'chemical_potentials': {'H': muH, 'X': 0}}]
occupy_structure_randomly(structure=structure,
                          cluster_space=cs,
                          target_concentrations={'A': {'Pd': 1 - cAu, 'Au': cAu},
                                                 'B': {'H': cH_start, 'X': 1 - cH_start}})

# step 3: construct ensemble and run
mc = HybridEnsemble(
    structure=structure,
    calculator=calculator,
    temperature=300,
    ensemble_specs=ensemble_specs,
    dc_filename=f'{output_directory}/hybrid-T{temp}-muH{muH:+.3f}-cAu{cAu:.3f}.dc')

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