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 Rahm et al. (2021)). This tutorial demonstrates how sublattice specific ensembles can be used in mchammer using the HybridEnsemble class. It is also shown how this class can be used to control the involved species and sites in a Monte Carlo simulation.

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.

from ase import Atom
from ase.build import bulk
from icet import ClusterSpace, ClusterExpansion
from mchammer.calculators import ClusterExpansionCalculator

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

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

from os import mkdir
from icet.tools.structure_generation import occupy_structure_randomly

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

Finally, we define our HybridEnsemble and commence the simulation.

from mchammer.ensembles import HybridEnsemble

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

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': 0.5}, 'kappa': 200}]

Controlling the species and sites involved

Within the HybridEnsemble, it is possible to also specify the allowed_species and allowed_sites, i.e., which chemical species and sites of the supercell that are allowed to be involved in a trial step. Here, we show how these functionalities can be used for a ternary surface system where only the surface region is allowed to change and only two of the chemical species are allowed in the surface region.

First, we setup a toy cluster expansion for a 6 layer AuCuPd surface slab.

from icet import ClusterSpace, ClusterExpansion
from ase.build import fcc111

prim = fcc111('Au', size=(1, 1, 6), a=4.0, vacuum=10, periodic=True)

# Set up toy Cluster Expansion
cs = ClusterSpace(structure=prim, cutoffs=[0.0], chemical_symbols=['Au', 'Cu', 'Pd'])
ce = ClusterExpansion(cluster_space=cs,
                      parameters=[0, -0.25, 0.5, 0.05, 0.02, 0.03, -0.02])

The primitive structure is a 6 layer surface slab, where the sites of lower and upper surface have indicies 0 and 5, respectively. We repeat the primitive structure to obtain a \(3\times3\times6\) supercell for MC. The sites of the lower surface of the supercell will then have indices 0, 6, 12, …, 48 and sites of the upper surface will have indices 5, 11, 17, …, 53.

from icet.tools.structure_generation import occupy_structure_randomly

# Supercell for MC with random occupation
structure = prim.repeat((3, 3, 1))
occupy_structure_randomly(structure=structure,
                          cluster_space=cs,
                          target_concentrations={'Au': 10/54, 'Cu': 10/54, 'Pd': 34/54})

# Set up a list of the allowed sites corresponding to the surface sites
allowed_sites = [i + j*6 for j in range(9) for i in [0, 5]]
print('Allowed sites: ', allowed_sites)

# Change any surface Cu to Pd
for i in allowed_sites:
    if structure[i].symbol == 'Cu':
        structure[i].symbol = 'Pd'

Allowed sites:  [0, 5, 6, 11, 12, 17, 18, 23, 24, 29, 30, 35, 36, 41, 42, 47, 48, 53]

Next, we set up a HybridEnsemble where we specify the allowed_sites from above as well as the allowed_symbols: ['Au', 'Pd'] to avoid Cu in the surface region. Then we run the simulation.

from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import HybridEnsemble

calculator = ClusterExpansionCalculator(structure, ce)

ensemble_specs = [{'ensemble': 'canonical', 'sublattice_index': 0,
                   'allowed_sites': allowed_sites, 'allowed_symbols': ['Au', 'Pd']}]

mc = HybridEnsemble(structure=structure, calculator=calculator,
                    temperature=300, ensemble_specs=ensemble_specs)

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

Lastly, we check all sites that were changed during the simulation and verify that only surface sites were affected and that no Cu were involved.

changed_sites = []
for i in range(len(structure)):
    if (structure[i].symbol != mc.structure[i].symbol):
        print(f'site {i}: {structure[i].symbol} -> {mc.structure[i].symbol}')
        changed_sites.append(i)

site 18: Au -> Pd
site 23: Pd -> Au