Random structure generation

In this notebook we will generate a large pool of randomized structures.

Random here refers to a random supercell sizes with a random occupation of the Mo-V and C-vacancy sublattices.

First some helper functions are defined to make generating randomized supercells easier.

import numpy as np
from ase.build import make_supercell


def find_sites(atoms, sym1, sym2):
    """
    Find the atomic indices corresponding to, e.g., metal sites in a structure.
    """
    mask_A = atoms.symbols == sym1
    mask_B = atoms.symbols == sym2
    mask = np.logical_or(mask_A, mask_B)
    sites = np.where(mask)[0]
    return sites


def generate_random_diagonal_P_matrix(n_atoms_in_prim, max_atoms, max_repeat_val):
    while True:
        nx = np.random.randint(1, max_repeat_val + 1)
        ny = np.random.randint(1, max_repeat_val + 1)
        nz = np.random.randint(1, max_repeat_val + 1)
        if nx * ny * nz * n_atoms_in_prim < max_atoms:
            break
    return np.diag([nx, ny, nz])


def randomly_occupy_supercell(atoms):
    """
    Randomly occupy a supercell with Mo-V on metal sites and
    C-vacancies on carbon sites.
    """

    # Find sites
    metal_sites = find_sites(atoms, *metal_symbols)
    carbon_sites = find_sites(atoms, 'C', 'Be')

    assert len(metal_sites) + len(carbon_sites) == len(atoms)
    assert len(set(metal_sites.tolist() + carbon_sites.tolist())) == len(atoms)

    # Metal occupation
    n_sites = len(metal_sites)
    sym_A, sym_B = metal_symbols
    n_A = np.random.randint(0, n_sites + 1)
    n_B = n_sites - n_A

    metal_occupations = [sym_A] * n_A + [sym_B] * n_B
    np.random.shuffle(metal_occupations)

    # Carbon occupations
    n_vacancy = np.random.randint(0, int(n_sites * x_vacancy_max) + 1)
    n_carbon = n_sites - n_vacancy

    carbon_occupations = ['C'] * n_carbon + ['Be'] * n_vacancy
    np.random.shuffle(carbon_occupations)

    # occupy supercell
    atoms.symbols[metal_sites] = metal_occupations
    atoms.symbols[carbon_sites] = carbon_occupations


def generate_random_supercell(prim, max_atoms, max_repeat_val):
    """ Generate a random (Mo, V)-(C, vacancy) supercell. """
    n_atoms_in_prim = len(prim)
    P = generate_random_diagonal_P_matrix(n_atoms_in_prim, max_atoms, max_repeat_val)
    atoms = make_supercell(prim, P)
    randomly_occupy_supercell(atoms)
    atoms.wrap()

    return atoms

We still need to define some parameters.

# number of structures to generate
n_structures = 5000

# maximum number of repetitions of the primitive unit cell per direction
max_repeat_val = 10

# maximum number of atoms
max_atoms = 50

# allowed species on the metal sublattice
metal_symbols = ['Mo', 'V']

# maximum carbon vacancy concentration
x_vacancy_max = 0.3

Now we are set up to generate configurations.

from ase.io import read

prim = read('../structures/MoC_rocksalt_prim.xyz')
n_atoms_in_prim = len(prim)

structures = []
for it in range(n_structures):
    atoms = generate_random_supercell(prim, max_atoms, max_repeat_val)
    structures.append(atoms)

We can generate histograms over the obtained concentrations.

from matplotlib import pyplot as plt

fig, axes = plt.subplots(
    figsize=(4.5, 4.5),
    dpi=120,
    nrows=2,
)

x_Mo = [atoms.symbols.count('Mo') / (atoms.symbols.count('Mo') + atoms.symbols.count('V')) for atoms in structures]
x_C = [atoms.symbols.count('C') / (atoms.symbols.count('C') + atoms.symbols.count('Be')) for atoms in structures]

ax = axes[0]
ax.hist(x_Mo, bins=50)
ax.set_xlabel('Mo concentration')
ax.set_ylabel('Number of structures')

ax = axes[1]
ax.hist(x_C, bins=50)
ax.set_xlabel('C concentration')
ax.set_ylabel('Number of structures')

fig.tight_layout()
fig.align_ylabels()