This page was generated from docs/src/notebooks/running-charge-balanced-gcmc.ipynb.

Running Semigrand Canonical Monte Carlo Sampling with charge-balance

import numpy as np
import json
from pymatgen.core.structure import Structure
from smol.io import load_work

0) Load the previous LNO CE with electrostatics

work = load_work('./data/basic_ce_ewald.mson')
expansion = work['ClusterExpansion']

1) Create a semigrand ensemble

The Ensemble class can also be used to run semigrand canonical MC by fixing relative chemical potentials.

In SGC chemical potential differences are set as boundary conditions. Any one of the active species can be used as reference.

from smol.moca import Ensemble

# Create the ensemble
# This specifies the size of the MC simulation domain.
sc_matrix = np.array([
    [6, 1, 1],
    [1, 2, 1],
    [1, 1, 2]
])

# relative chemical potentials are provided as a dict
chemical_potentials = {'Li+': 0, 'Vacancy': 0, 'Ni3+': 0, 'Ni4+': 0}

# this convenience method will take care of creating the appropriate
# processor for the given cluster expansion.
ensemble = Ensemble.from_cluster_expansion(
    expansion, sc_matrix, chemical_potentials=chemical_potentials
)

2) Create an MC sampler

A Sampler will take care of running MC sampling runs for a given ensemble. For a charge-balanced simulation, one can choose either a step type using: 1) the table-exchange method (TE, for historical reason, written as TableFlip in the code) 2) the square-charge bias method (SCB).

For TE, one needs to specify the step type as “table-flip”. Additional arguments to control the step type can be specified directly as keyword arguments when creating the sampler. To see what arguments are available refer to the documentation of the TableFlip class.

from smol.moca import Sampler

# For the usage of other keywords, see documentation of the TableFlip class.
sampler_table = Sampler.from_ensemble(ensemble, temperature=500, step_type="table-flip", optimize_basis=True)

print(f"Sampling information: {sampler_table.samples.metadata}")
print(f"Table exchanges in reaction formulas: {sampler_table.mckernels[0].mcusher._comp_space.flip_reactions}")

Sampling information: {'kernel': 'Metropolis', 'step': 'table-flip', 'chemical_potentials': {Species Li+: 0, Vacancy vacA0+: 0, Species Ni3+: 0, Species Ni4+: 0}, 'seeds': [97098921871036988304502408163907614831]}
Table exchanges in reaction formulas: ['1 Li+(0) + 1 Ni3+(1) -> 1 vacA0+(0) + 1 Ni4+(1)']

from pymatgen.transformations.standard_transformations import OrderDisorderedStructureTransformation

# Here we will just use the order disordered transformation from
# pymatgen to get an ordered version of a prim supercell.
# The structure will have the same composition set in the prim.
transformation = OrderDisorderedStructureTransformation()

supercell = expansion.cluster_subspace.structure.copy()
supercell.make_supercell(sc_matrix)

# this can take a bit of time....
test_struct = transformation.apply_transformation(supercell)
print(test_struct.composition)

# Obtain the initial occupancy string from the
# test structure created above.
init_occu = ensemble.processor.occupancy_from_structure(test_struct)

# The occupancy strings created by the processor
# are by default "encoded" by the indices of the species
# for each given site. You can always see the actual
# species in the occupancy string by decoding it.
print(f'The encoded occupancy is:\n{init_occu}')
print(f'The initial occupancy is:\n {ensemble.processor.decode_occupancy(init_occu)}')

Li+8 Ni3+8 Ni4+8 O2-32
The encoded occupancy is:
[0 1 1 0 0 0 1 0 0 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1 0 0 1 0 0 0 1 1 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
The initial occupancy is:
 [Species Li+, Vacancy vacA0+, Vacancy vacA0+, Species Li+, Species Li+, Species Li+, Vacancy vacA0+, Species Li+, Species Li+, Vacancy vacA0+, Species Li+, Vacancy vacA0+, Vacancy vacA0+, Vacancy vacA0+, Vacancy vacA0+, Species Li+, Species Ni3+, Species Ni3+, Species Ni4+, Species Ni4+, Species Ni4+, Species Ni4+, Species Ni3+, Species Ni4+, Species Ni3+, Species Ni3+, Species Ni4+, Species Ni3+, Species Ni3+, Species Ni3+, Species Ni4+, Species Ni4+, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-, Species O2-]

# Run a short simulation.
sampler_table.run(50000, init_occu, thin_by=10, progress=True)
energy_avg = sampler_table.samples.mean_energy(discard=sampler_table.samples.num_samples // 4)
comp_avg = sampler_table.samples.mean_composition(discard=sampler_table.samples.num_samples // 4)
print(f"Average energy of charge balanced samples (eV): {energy_avg}")
print(f"Average composition of charge balanced samples: {comp_avg}")

Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 50000/50000 [00:24<00:00, 2055.76it/s]

Average energy of charge balanced samples (eV): -599.4992302891402
Average composition of charge balanced samples: {Species Li+: 0.25, Vacancy vacA0+: 0.0, Species Ni3+: 0.25, Species Ni4+: 0.0, Species O2-: 0.5}

For SCB, one can specify a square charge bias in sampler (if other constraints are to be enforced, refer to SquareHyperplaneBias):

sampler_bias = Sampler.from_ensemble(ensemble, temperature=500,
                                      step_type="flip",
                                      bias_type="square-charge", bias_kwargs={"penalty": 1.0})

print(f"Sampling information: {sampler_bias.samples.metadata}")

Sampling information: {'kernel': 'Metropolis', 'step': 'flip', 'chemical_potentials': {Species Li+: 0, Vacancy vacA0+: 0, Species Ni3+: 0, Species Ni4+: 0}, 'seeds': [159094150126610822048360871456393092894]}

# SCB needs to filter out charge unbalanced configurations!
sampler_bias.run(100000, init_occu, thin_by=10, progress=True)

bias = sampler_bias.samples.get_trace_value("bias", discard=sampler_bias.samples.num_samples // 4)
energies = sampler_bias.samples.get_energies(discard=sampler_bias.samples.num_samples // 4)
energy_avg = np.average(energies[np.isclose(bias, 0)])
print(f"Average energy of charge balanced samples (eV): {energy_avg}")

Sampling 1 chain(s) from a cell with 64 sites...: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100000/100000 [00:11<00:00, 8584.75it/s]

Average energy of charge balanced samples (eV): -599.4992302891402