Local cluster expansion for kinetically resolved activation barriers

icet
1

I’m currently working on building a LCE for KRA barriers according to Van der Ven et al. (DOI: 10.1103/PhysRevB.64.184307) and Zhang and Sluiter (DOI: 10.1007/s11669-015-0427-x). My goal is to train two CEs (one for configurational (initial/final state) energies and one for KRA barriers) to calculate energy barriers for a kinetic Monte Carlo simulation in which I am modelling H diffusion on the tetrahedral interstitial sites in vacancy-free BCC Fe. I have used ICET to train a CE to predict configurational energies, but I’m currently struggling to figure out how to use a LocalOrbitListGenerator and generate_local_orbit_list to customize the cluster space to build an LCE (and I have not been able to find any tutorials besides ICET’s “Customizing cluster spaces” one, which doesn’t mention these classes/functions). I have labelled the “swapping pair” (H-occupied T site and empty adjacent T site) in the training configurations according to Zhang and Sluiter’s description (alternatively, I tried adding another sublattice whose sites are midpoints between adjacent T sites to represent transition states, but this made training the CE extremely slow), like so (partial code):

import numpy as np
from icet.core.structure import Structure
from icet.core.local_orbit_list_generator import LocalOrbitListGenerator 
from icet.tools import Constraints
from icet import ClusterSpace, StructureContainer, ClusterExpansion
structure = Structure.from_atoms(TEMPLATE.atoms)
cs = ClusterSpace(TEMPLATE.atoms, cutoffs=[3.0,3.0,3.0,3.0], chemical_symbols=chemical_symbols)
lolg = LocalOrbitListGenerator(orbit_list=cs.orbit_list,structure=structure,fractional_position_tolerance=cs.fractional_position_tolerance)
offset = (np.floor(np.round(np.dot(MID_CART_POS, np.linalg.inv(cs.primitive_structure.cell[:])),6)).astype(int).tolist())

supercell_orbit_list = lolg.generate_local_orbit_list(offset=offset, self_contained=True)

But now I am unsure of how to proceed to creating constraints based on the included/excluded orbits (this is actually my question).

sc = StructureContainer(cs)
#populate sc with training configurations
constraints = Constraints(???)
x, y = sc.get_fit_data()
x_constrained = constraints.transform(x)

opt = CrossValidationEstimator(fit_data=(x_constrained, y), fit_method='ardr', max_iter=100000)
opt.validate()
opt.train()
print(opt)

params = constraints.inverse_transform(opt.parameters)

ce = ClusterExpansion(cs, params)

While stuck on this, my alternative approach to calculating the barriers involves training a CE on initial/final state energies, and then another CE on full transition state energies (which is not an LCE) with the swapping pair labelled (in both cases, the RMSE (~0.2 eV) is significantly larger than the barrier heights (~0.043 eV) I need to calculate). For additional context: all training configurations are cubic 5 x 5 x 5 BCC Fe supercells, with a lattice parameter of 2.867 A, 0% to 5% of tetrahedral interstitials occupied by H, generated using ASE’s “crystal()” with relaxation done using LAMMPS, LBFGS (fmax=1e-5, steps=2000), Ito’s ML IAP (DOI: 10.1016/j.ijhydene.2026.155600) with ZPE correction done using Vibrations (delta=0.01) on only the H atoms (for improved efficiency, given the large mass disparity between Fe and H, as justified by Jiang and Carter, DOI: 10.1103/PhysRevB.70.064102), and saddle point searches calculated using CI-NEB (first without climbing image, fmax=0.05 and steps=1000; then with climbing image, fmax=0.02 and steps=1000) and refined using Dimer (fmax=1e-5, steps=5000); configurations that did not converge were excluded from the training set). I’m hoping that an LCE for KRA barriers would exhibit better fitting stats. I have also tried using KMCpy, but sparse documentation made it difficult to figure out. Any advice/help would be appreciated :slight_smile:

2

I think I might have answered my own question, but any critique/corrections are welcome; I am here to learn :blush:

Much helpful information is available on ICET’s “Training with constraints and weights” advanced tutorial. I also found this article super helpful: https://doi.org/10.1103/PRXEnergy.3.042001

In case helpful to others, here is my solution:

import numpy as np
from ase import Atoms
from icet import ClusterSpace, StructureContainer, ClusterExpansion
from icet.tools import Constraints
from trainstation import CrossValidationEstimator


def train_kra_ce(configurations: list[tuple[Atoms,float,dict[str,float]]]):
    print("Building cluster space and structure container...")

    chemical_symbols: list[list[str]] = [[""]] * len(TEMPLATE.atoms)

    for idx in TEMPLATE.host_indices:
        chemical_symbols[idx] = [HOST_SYMBOL]

    for idx in TEMPLATE.tetrahedral_interstitial_indices:
        chemical_symbols[idx] = [VACANCY_SYMBOL, INTERSTITIAL_SYMBOL, SWAPPING_PAIR_SYMBOL]

    cs = ClusterSpace(TEMPLATE.atoms, cutoffs=TRANSITION_CUTOFFS, chemical_symbols=chemical_symbols)
    cs.write(KRA_CS_NAME)

    from icet.core.structure import Structure
    from icet.core.local_orbit_list_generator import LocalOrbitListGenerator
    structure = Structure.from_atoms(TEMPLATE.atoms)
    lolg = LocalOrbitListGenerator(orbit_list=cs.orbit_list, structure=structure,
                                   fractional_position_tolerance=cs.fractional_position_tolerance)
    orbit_list = lolg.generate_full_orbit_list()

    active_site_indices = {SRC_IDX, DST_IDX}
    exclude_orbit_indices = set()
    orbits = []
    for orbit_idx, orbit in enumerate(orbit_list.orbits):
        orbits.append([])
        for cluster in orbit.clusters:
            cluster_lattice_site_indices = {ls.index for ls in cluster.lattice_sites}
            if active_site_indices.issubset(cluster_lattice_site_indices):
                orbits[orbit_idx].append(cluster)
        if not orbits[orbit_idx]:
            exclude_orbit_indices.add(orbit_idx)

    print(f"Excluding: {len(exclude_orbit_indices)} / {len(cs.orbit_list.orbits)} orbits")

    exclude_correlation_func_indices = set()
    for correlation_func in cs.as_list:
        if correlation_func['orbit_index'] in exclude_orbit_indices:
            exclude_correlation_func_indices.add(correlation_func['index'])

    print(f"Excluding: {len(exclude_correlation_func_indices)} / {len(cs.as_list)} correlation functions")

    sc = StructureContainer(cs)

    print(f"Adding structures to structure container...")
    for atoms, energy, data in configurations:
        try:
            init_egy = data['init zpe corr']
            fin_egy = data['fin zpe corr']
            ts_egy = data['ts zpe corr']
            kra_egy = ts_egy - (init_egy + fin_egy) * 0.5
            sc.add_structure(atoms, properties={"energy": kra_egy})
            print(f"{atoms}; KRA barrier: {kra_egy}; data: {data}")
        except Exception as e:
            print(f"ERROR, SKIPPED: {e}")

    sc.write(KRA_SC_NAME)

    n_params = len(cs)
    n_constraints = len(exclude_correlation_func_indices)
    matrix = np.zeros((n_constraints, n_params))
    for i, param_idx in enumerate(exclude_correlation_func_indices):
        matrix[i, param_idx] = 1.0

    constraints = Constraints(n_params=n_params)
    constraints.add_constraint(matrix)

    x, y = sc.get_fit_data()
    x_constrained = constraints.transform(x)

    print("Training optimizer...")
    opt = CrossValidationEstimator(fit_data=(x_constrained,y), fit_method='ardr', max_iter=100000)
    opt.validate()
    opt.train()
    print(opt)

    params = constraints.inverse_transform(opt.parameters)

    print("Constructing cluster expansion for transition states...")
    ce = ClusterExpansion(cs, params)
    ce.write(KRA_CE_NAME)
    print(ce)

However, fitting using 1500 configurations and their KRA barriers yields a few high correlation number warnings (~1e+18) and unfortunately the cross-validation stats are still quite poor (RMSE is comparable to barriers and R2 is <<1). I’ve also noticed these stats aren’t much different when I include all orbits/correlation functions without zeroing anything.

The code below is my initial simpler approach of fitting to transition state configurational energies instead of local KRA barriers. Perhaps I should rather stick to this?

def train_transition_ce(configurations: list[tuple[Atoms,float,dict[str,float]]]):
    print("Building cluster space and structure container...")

    chemical_symbols: list[list[str]] = [[""]] * len(TEMPLATE.atoms)

    for idx in TEMPLATE.host_indices:
        chemical_symbols[idx] = [HOST_SYMBOL]

    for idx in TEMPLATE.tetrahedral_interstitial_indices:
        chemical_symbols[idx] = [VACANCY_SYMBOL, INTERSTITIAL_SYMBOL, SWAPPING_PAIR_SYMBOL] 

    cs = ClusterSpace(TEMPLATE.atoms, cutoffs=TRANSITION_CUTOFFS, chemical_symbols=chemical_symbols)
    cs.write(TRANSITION_CS_NAME)

    sc = StructureContainer(cs)

    print(f"Adding structures to structure container...")
    for atoms, energy, data in configurations:
        try:
            sc.add_structure(atoms, properties={"energy": energy})
            print(f"{atoms}; {energy}; data: {data}")
        except Exception as e:
            print(f"ERROR, SKIPPED: {e}")

    sc.write(TRANSITION_SC_NAME)

    fit_data = sc.get_fit_data()
    print(fit_data)

    print("Training optimizer...")
    opt = CrossValidationEstimator(fit_data=fit_data, fit_method='ardr', max_iter=100000)
    opt.validate()
    opt.train()
    print(opt)

    print("Constructing cluster expansion for transition states...")
    ce = ClusterExpansion(cs, opt.parameters)
    ce.write(TRANSITION_CE_NAME)
    print(ce)
3

Hi @Hayley_Britz

I am happy to see that you took such a deep dive into the icet core. Incidentally we are currently working on documenting and cleaning up the C++ core (see this MR on gitlab), since as you saw there documentation gaps, and we actually are working on extending icet for kMC simulations. One thing you might have to look out for when using the LOLG is whether your cluster vectors are additive or not. In the MR linked above is now a function for that purpose are_local_cluster_vectors_additive.

4

Thank you for your reply, Prof Erhart! I appreciate your insight.

In my code, replacing:

from icet.core.structure import Structure
from icet.core.local_orbit_list_generator import LocalOrbitListGenerator
structure = Structure.from_atoms(TEMPLATE.atoms)
lolg = LocalOrbitListGenerator(orbit_list=cs.orbit_list, structure=structure,
                               fractional_position_tolerance=cs.fractional_position_tolerance)
orbit_list = lolg.generate_full_orbit_list()

with:

orbit_list = cs.orbit_list.get_supercell_orbit_list(structure=TEMPLATE.atoms,
    fractional_position_tolerance=cs.fractional_position_tolerance)

yields the exact same orbit list, and testing the code snippet from the MR you linked:

for orbit in orbit_list.orbits:
    for cluster in orbit.clusters:
        indices = [site.index for site in cluster.lattice_sites]
        if len(set(indices)) != len(indices):
            return False
return True

returns True, which I assume is good?

5

Sorry the late reply. I had missed the post. Yes. In principle that is the correct approach. Please note that in icet 4.0 (which was released yesterday) the core has been rewritten and there are some breaking changes, which affect the code above. A key motivation for the rewrite was in fact to make it easier for users to write additional calculators, including support for local CEs and kMC.