This page was generated from
docs/src/notebooks/running-semigrand-mc.ipynb.
Running Semigrand Canonical Monte Carlo Sampling#
[1]:
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#
[2]:
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.
[3]:
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
)
print(f'The supercell size for the processor is {ensemble.processor.size} prims.\n')
print(f'The ensemble has a total of {ensemble.num_sites} sites.\n')
print(f'The active sublattices are:')
for sublattice in ensemble.sublattices:
print(sublattice)
The supercell size for the processor is 16 prims.
The ensemble has a total of 64 sites.
The active sublattices are:
Sublattice(site_space=Li+0.5 vacA0+0.5, sites=array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]), active_sites=array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]), encoding=array([0, 1]))
Sublattice(site_space=Ni3+0.5 Ni4+0.5 , sites=array([16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]), active_sites=array([16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]), encoding=array([0, 1]))
Sublattice(site_space=O2-1 , sites=array([32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]), active_sites=array([], dtype=float64), encoding=array([0]))
2) Create an MC sampler#
A Sampler will take care of running MC sampling runs for a given ensemble. The sampler allows many different options for MC sampling most importantly setting the MCMC algorithm and the type of MC steps taken. However the defaults are usually enough for almost all use cases.
[4]:
from smol.moca import Sampler
# This will take care of setting the defaults
# for the supplied canonical ensemble
sampler = Sampler.from_ensemble(ensemble, temperature=500)
print(f"Sampling information: {sampler.samples.metadata}")
Sampling information: Metadata(chemical_potentials={Species Li+: 0, Vacancy vacA0+: 0, Species Ni3+: 0, Species Ni4+: 0}, cls_name='SampleContainer', kernels=[Metadata(seed=253072746018418677829111918777968458889, step=Metadata(sublattices=[(Species Li+, Vacancy vacA0+), (Species Ni3+, Species Ni4+), (Species O2-,)], sublattice_probabilities=array([0.5, 0.5]), cls_name='Flip'), cls_name='Metropolis')])
3) Create an initial structure and get occupancies#
[5]:
from smol.capp.generate import generate_random_ordered_occupancy
init_occu = generate_random_ordered_occupancy(ensemble.processor)
print(f"The disordered structure has composition: {ensemble.processor.structure.composition}")
print(f"The initial occupancy has composition: {ensemble.processor.structure_from_occupancy(init_occu).composition}")
The disordered structure has composition: Li+8 Ni3+8 Ni4+8 O2-32
The initial occupancy has composition: Li+11 Ni4+10 Ni3+6 O2-32
[6]:
# 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)}')
The encoded occupancy is:
[0 0 0 0 0 0 0 0 0 1 0 1 1 1 1 0 1 1 0 1 0 1 1 0 1 1 0 1 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+, Species Li+, Species Li+, Species Li+, Species Li+, Species Li+, Species Li+, Species Li+, Species Li+, Vacancy vacA0+, Species Li+, Vacancy vacA0+, Vacancy vacA0+, Vacancy vacA0+, Vacancy vacA0+, Species Li+, Species Ni4+, Species Ni4+, Species Ni3+, Species Ni4+, Species Ni3+, Species Ni4+, Species Ni4+, Species Ni3+, Species Ni4+, Species Ni4+, Species Ni3+, Species Ni4+, 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-]
4) Sweeping chemical potentials#
Using a semigrand ensemble, we can sweep through relative chemical potentials and get average energies at various compositions.
Here we will sweep only through the relative chemical potential of vacancies.
[7]:
from smol.cofe.space import Vacancy
dmus = np.arange(4, 7, .25)
discard, thin_by = 50, 10
concentration_avgs = []
energy_avg = []
sampler.clear_samples()
# initialize sampler
chemical_potentials['Vacancy'] = dmus[0]
ensemble.chemical_potentials = chemical_potentials
sampler.run(10000, init_occu, thin_by=thin_by, progress=True)
energy_avg.append(sampler.samples.mean_energy(discard=discard))
# this gets the composition in all active sites
# excludes inactive sites (i.e. O2- in this case)
# to get the composition per sublattice only
# mean_sublattice_composition
concentration_avgs.append(sampler.samples.mean_composition(discard=discard))
# sweep through rest
for dmu in dmus[1:]:
# update to get only averages of last samples
# you can also clear samples and start fresh at each dmu
discard_tot = discard + len(sampler.samples)
chemical_potentials["Vacancy"] = dmu
ensemble.chemical_potentials = chemical_potentials
sampler.run(100000, thin_by=thin_by, progress=True)
energy_avg.append(sampler.samples.mean_energy(discard=discard_tot))
concentration_avgs.append(sampler.samples.mean_composition(discard=discard_tot))
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 10000/10000 [00:04<00:00, 2375.34it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:42<00:00, 2374.55it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:41<00:00, 2385.60it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:43<00:00, 2305.62it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:43<00:00, 2324.53it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:44<00:00, 2232.78it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:42<00:00, 2335.12it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:43<00:00, 2297.89it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:42<00:00, 2359.55it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:41<00:00, 2382.00it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:54<00:00, 1843.77it/s]
Sampling 1 chain(s) from a cell with 64 sites...: 100%|██████████| 100000/100000 [00:50<00:00, 1988.32it/s]
[8]:
# Plot some results
import matplotlib.pyplot as plt
# We need to use this since species keys are pymatgen Specie objects
from smol.cofe.space import get_species
plt.plot([c[get_species("Li+")] for c in concentration_avgs],
energy_avg, 'o')
plt.xlabel('Li+ concentration')
plt.ylabel('Average total energy (eV)')
plt.show()
