NEB¶

Introduction¶

An algorithm to find the transition state between the initial and final images

• Reference: https://theory.cm.utexas.edu/henkelman/research/saddle/neb/
• Key papers: https://theory.cm.utexas.edu/henkelman/pubs/henkelman00_9978.pdf (https://aip.scitation.org/doi/10.1063/1.1323224)
• Climbing image NEB (CI-NEB) vs normal NEB calculations: CI-NEB can be activated using LCLIMB = True in INCAR. However, according to our experience, this algorithm is not as stable as regular NEB.

Procedure¶

1. Create structures for initial and final state -> by hand or using Pymatgen/NEB
2. Make sure that the coordinates of the initial and final images are matching with each other in sequence. It's a good practice to put the migrating specie (ion) to the top of the coordinate list.
3. Fully relax both initial and final images
4. Take the relaxed initial and final images, then interpolate to get intermediate images.
5. Before doing NEB calculations, check your trajectory and see if it is reasonable.
6. Run NEB optimization.
7. After NEB optimization, get the energy profile and plot the results.

Tips¶

1. You might need to do NEB in a supercell if your unit cell is too small (lattice parameters are less than 8-10 Å) to get rid of image interactions to the next supercell.
2. When creating initial and final images, you might need to create vacancies. Then you might have problem of charge balancing.
3. Usually initial and final images are the nearest neighboring migrating species. If the distances are large (~ 5 Å), it's better to do a bond valance calculation (SoftBV: https://www.dmse.nus.edu.sg/asn/softBV-GUI.html or ref: https://pubs.acs.org/doi/10.1021/acs.chemmater.0c03893) to estimate the diffusion path.
4. If VASP is used as a calculator, NEB can only be run on large machines such as NSCC because the resource consumption of each image is equivalent to one relaxation, so the total resource consumption will be (N_image-2) * Resource_single_relax (initial and final are excluded)
5. NEB optimizations might be hard to converge than regular relaxations, make sure your initial migration path physically make sense.

Example of NEB using MACE through ASE¶

In [ ]:

from ase.mep import NEB
from ase import Atoms
from ase.optimize import LBFGS
from mace.calculators import MACECalculator
from ase.io import Trajectory, read, write
import pandas as pd
import os

class NotConvergedError(Exception):
    pass

def geom_optimization(atoms:Atoms,
                      calculator:MACECalculator,
                      ofname:str)->bool:
    calc = calculator
    atoms.calc = calc
    # filter = UnitCellFilter(atoms)
    traj = Trajectory(ofname, 'w', atoms)
    optimization = LBFGS(atoms)
    optimization.attach(traj)
    optimization.run(fmax=0.01,steps=2000)
    return optimization.converged()

def neb(fname_initial:str,
        fname_final:str,
        nimages:int)->bool:
    initial = read(fname_initial)
    final = read(fname_final)
    # Make a band consisting of N images:
    images = [initial]
    images += [initial.copy() for i in range(nimages-2)]
    images += [final]
    neb = NEB(images)
    # Interpolate linearly the potisions of the three middle images:
    neb.interpolate()

    try:
        os.mkdir('neb')
    except FileExistsError:
        pass
    for i, image in enumerate(neb.images):
        write(filename=f'neb/image_{i}.cif',images=image,format='cif')

    # Set calculators:
    for image in images[1:nimages]:
        image.calc = MACECalculator(model_paths="mace_agnesi_medium.model", 
                                    dispersion=False, 
                                    default_dtype="float64", 
                                    device='cuda')
    # Optimize:
    print('Runing NEB...')
    optimizer = LBFGS(neb) # you can try other optimizer if LBFGS doesn't work
    optimizer.run(fmax=0.05,steps=2000) # NEB is harder to converge than geometry optimization

    try:
        os.mkdir('neb_final_images')
    except FileExistsError:
        pass
    for i, image in enumerate(neb.images):
        write(filename=f'neb_final_images/image_{i}.cif',images=image,format='cif')

    energies = neb.energies
    energies[0] = initial.get_potential_energy()
    energies[-1] = final.get_potential_energy()

    df = pd.DataFrame(energies,columns=['energies'])
    df.to_csv('neb_results.csv')
    return optimizer.converged()


def run_geometry_optimization(fname_initial_structure:str,
                              fname_final_structure:str)->bool:
    calculator = MACECalculator(model_paths="mace_agnesi_medium.model", 
                                dispersion=False, 
                                default_dtype="float64", 
                                device='cuda')
    # Initial and final structures are created manually
    initial_structure = read(fname_initial_structure)
    initial_converged = geom_optimization(initial_structure,calculator,'rlx_initial.traj')
    final_structure = read(fname_final_structure)
    final_converged = geom_optimization(final_structure,calculator,'rlx_final.traj')
    return initial_converged and final_converged

def run_neb(fname_rlx_initial:str,
            fname_rlx_final:str, 
            converged:bool, 
            nimages:int):
    if converged:
        neb_converged = False
        try:
            neb_converged = run_neb(fname_rlx_initial,fname_rlx_final, nimages)
            if not neb_converged:
                raise NotConvergedError('NEB calculation did not converge...')
        except NotImplementedError:
            print('NotImplementedError occurs')
            pass
    else:
        raise NotConvergedError('Geometry optimization did not converge, skip NEB calculation...')



rlx_converged = run_geometry_optimization('ALF_initial.cif','ALF_final.cif')
run_neb('neb_initial_rlx.traj','neb_final_rlx.traj', rlx_converged, 11)

from ase.mep import NEB from ase import Atoms from ase.optimize import LBFGS from mace.calculators import MACECalculator from ase.io import Trajectory, read, write import pandas as pd import os class NotConvergedError(Exception): pass def geom_optimization(atoms:Atoms, calculator:MACECalculator, ofname:str)->bool: calc = calculator atoms.calc = calc # filter = UnitCellFilter(atoms) traj = Trajectory(ofname, 'w', atoms) optimization = LBFGS(atoms) optimization.attach(traj) optimization.run(fmax=0.01,steps=2000) return optimization.converged() def neb(fname_initial:str, fname_final:str, nimages:int)->bool: initial = read(fname_initial) final = read(fname_final) # Make a band consisting of N images: images = [initial] images += [initial.copy() for i in range(nimages-2)] images += [final] neb = NEB(images) # Interpolate linearly the potisions of the three middle images: neb.interpolate() try: os.mkdir('neb') except FileExistsError: pass for i, image in enumerate(neb.images): write(filename=f'neb/image_{i}.cif',images=image,format='cif') # Set calculators: for image in images[1:nimages]: image.calc = MACECalculator(model_paths="mace_agnesi_medium.model", dispersion=False, default_dtype="float64", device='cuda') # Optimize: print('Runing NEB...') optimizer = LBFGS(neb) # you can try other optimizer if LBFGS doesn't work optimizer.run(fmax=0.05,steps=2000) # NEB is harder to converge than geometry optimization try: os.mkdir('neb_final_images') except FileExistsError: pass for i, image in enumerate(neb.images): write(filename=f'neb_final_images/image_{i}.cif',images=image,format='cif') energies = neb.energies energies[0] = initial.get_potential_energy() energies[-1] = final.get_potential_energy() df = pd.DataFrame(energies,columns=['energies']) df.to_csv('neb_results.csv') return optimizer.converged() def run_geometry_optimization(fname_initial_structure:str, fname_final_structure:str)->bool: calculator = MACECalculator(model_paths="mace_agnesi_medium.model", dispersion=False, default_dtype="float64", device='cuda') # Initial and final structures are created manually initial_structure = read(fname_initial_structure) initial_converged = geom_optimization(initial_structure,calculator,'rlx_initial.traj') final_structure = read(fname_final_structure) final_converged = geom_optimization(final_structure,calculator,'rlx_final.traj') return initial_converged and final_converged def run_neb(fname_rlx_initial:str, fname_rlx_final:str, converged:bool, nimages:int): if converged: neb_converged = False try: neb_converged = run_neb(fname_rlx_initial,fname_rlx_final, nimages) if not neb_converged: raise NotConvergedError('NEB calculation did not converge...') except NotImplementedError: print('NotImplementedError occurs') pass else: raise NotConvergedError('Geometry optimization did not converge, skip NEB calculation...') rlx_converged = run_geometry_optimization('ALF_initial.cif','ALF_final.cif') run_neb('neb_initial_rlx.traj','neb_final_rlx.traj', rlx_converged, 11)

Plot results.

In [ ]:

import matplotlib.pyplot as plt
import pandas as pd
df = pd.read_csv('neb_results.csv')
energies = df['energies'].to_list()
e_ref = (energies[0] + energies[-1])/2
df['energies'] = (df['energies'] - e_ref) * 1000 # to meV and adjust to reference
fig, ax = plt.subplots(1,1, figsize=(4,4))
ax.set_xlabel('Images')
ax.set_ylabel('Energies (meV)')
plt.plot(df['energies'],'o--')
fig.tight_layout()
fig.savefig('neb.pdf')

import matplotlib.pyplot as plt import pandas as pd df = pd.read_csv('neb_results.csv') energies = df['energies'].to_list() e_ref = (energies[0] + energies[-1])/2 df['energies'] = (df['energies'] - e_ref) * 1000 # to meV and adjust to reference fig, ax = plt.subplots(1,1, figsize=(4,4)) ax.set_xlabel('Images') ax.set_ylabel('Energies (meV)') plt.plot(df['energies'],'o--') fig.tight_layout() fig.savefig('neb.pdf')