Anharmonic energy surfaces

In this tutorial we study the potential energy surface (PES) for a (110) Cu slab. The atoms at the surface experience a strong anharmonic PES due to the asymmetry in bonds normal to the surface. We test how well 2nd and 3rd order models can describe this PES by examining their capacity to predict the energy landscape when shifting the entire top surface layer up and down.

Energy surface obtained from the full model (EMT) as well as a 2nd and 2nd+3rd order force constant potential.

By carrying out ensemble fitting (bagging) we obtain a distribution of force constant models. Using these models we can estimate how sensitive our model is with respect to a complex property such as the PES.

Several 2nd+3rd models compared with the ensemble average (bold black line).

This also illustrates the advantage of ensemble fitting and using averaged parameters. The average model is much more stable and less prone to overfitting.

It is important to remember though to feed the model with good data, i.e. if we want to sample the PES with displacements up to 0.5 A we should not train the model with configurations rattled with 0.01 A.

Source code

from hiphive.structure_generation import generate_rattled_structures
from ase.build import fcc110
from ase.io import write
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from hiphive.utilities import prepare_structures

# parameters
number_of_structures = 5
rattle_std = 0.05
surface_size = (3, 4, 8)
structures_fname = 'rattled_structures.extxyz'

# setup atoms and calculator
atoms_ideal = fcc110('Cu', size=surface_size)
calc = EMT()

atoms_ideal.center(vacuum=20, axis=2)
atoms_ideal.pbc = True

# relax structure
atoms_ideal.set_calculator(calc)
dyn = BFGS(atoms_ideal)
converged = dyn.run(fmax=0.0001, steps=1000)


# generate rattled structures
structures = generate_rattled_structures(atoms_ideal, number_of_structures, rattle_std)
structures = prepare_structures(structures, atoms_ideal, calc)

# save structures
write(structures_fname, structures)

from ase.io import read
from hiphive import StructureContainer
from hiphive import ClusterSpace


# setup
rattled_structures = read('rattled_structures.extxyz@:')

# Second order model
cs = ClusterSpace(rattled_structures[0], [5.1])
sc = StructureContainer(cs)
for structure in rattled_structures:
    sc.add_structure(structure)
sc.write('structure_container2')


# Third order model
cs = ClusterSpace(rattled_structures[0], [5.1, 3.0])
sc = StructureContainer(cs)
for structure in rattled_structures:
    sc.add_structure(structure)
sc.write('structure_container3')

import numpy as np
import matplotlib.pyplot as plt
from ase.io import read
from ase.calculators.emt import EMT
from hiphive import StructureContainer, ForceConstantPotential
from trainstation import Optimizer
from hiphive.calculators import ForceConstantCalculator
from tools import compute_energy_landscape


# parameters
dz_vals = np.linspace(-0.4, 0.4, 21)
atoms_ideal = read('rattled_structures.extxyz')

# read data
sc2 = StructureContainer.read('structure_container2')
sc3 = StructureContainer.read('structure_container3')
cs2 = sc2.cluster_space
cs3 = sc3.cluster_space

# fit models
opt = Optimizer(sc2.get_fit_data(), train_size=1.0)
opt.train()
fcp2 = ForceConstantPotential(cs2, opt.parameters)

opt = Optimizer(sc3.get_fit_data(), train_size=1.0)
opt.train()
fcp3 = ForceConstantPotential(cs3, opt.parameters)

# test models
emt_calc = EMT()
calc_fc2 = ForceConstantCalculator(fcp2.get_force_constants(atoms_ideal))
calc_fc3 = ForceConstantCalculator(fcp3.get_force_constants(atoms_ideal))

pes_emt = compute_energy_landscape(atoms_ideal, emt_calc, dz_vals)
pes_fc2 = compute_energy_landscape(atoms_ideal, calc_fc2, dz_vals)
pes_fc3 = compute_energy_landscape(atoms_ideal, calc_fc3, dz_vals)


# plot pes
fs = 14
lw = 1.5
ms = 8

plt.plot(pes_emt[:, 0], pes_emt[:, 1], '-o', lw=lw, ms=ms, label='EMT')
plt.plot(pes_fc2[:, 0], pes_fc2[:, 1], '-o', lw=lw, ms=ms, label='2nd')
plt.plot(pes_fc3[:, 0], pes_fc3[:, 1], '-o', lw=lw, ms=ms, label='2nd+3rd')
plt.xlim([np.min(dz_vals), np.max(dz_vals)])
plt.ylim(bottom=0.0)

plt.xlabel('Surface layer shift ($\\mathrm{\\AA}$)', fontsize=fs)
plt.ylabel('Potential Energy (eV)', fontsize=fs)

plt.gca().tick_params(labelsize=fs)

plt.legend(loc=9, fontsize=fs)
plt.tight_layout()
plt.savefig('energy_surface.svg')

import numpy as np
import matplotlib.pyplot as plt
from ase.io import read
from hiphive import StructureContainer
from hiphive.force_constant_model import ForceConstantModel
from trainstation import EnsembleOptimizer
from hiphive.calculators import ForceConstantCalculator
from tools import compute_energy_landscape


# parameters
ensemble_size = 20
dz_vals = np.linspace(-0.4, 0.4, 21)
atoms_ideal = read('rattled_structures.extxyz')


# read data
sc = StructureContainer.read('structure_container3')
cs = sc.cluster_space

# training
eopt = EnsembleOptimizer(sc.get_fit_data(), train_size=0.5, ensemble_size=ensemble_size)
eopt.train()

# average model
fcm = ForceConstantModel(atoms_ideal, cs)
fcm.parameters = eopt.parameters
calc_ave = ForceConstantCalculator(fcm.get_force_constants())
pes_ave = compute_energy_landscape(atoms_ideal, calc_ave, dz_vals)

# ensemble models
ensemble_pes = []
for parameters in eopt.parameters_splits:
    fcm.parameters = parameters
    calc = ForceConstantCalculator(fcm.get_force_constants())
    pes = compute_energy_landscape(atoms_ideal, calc, dz_vals)
    ensemble_pes.append(pes)

# plot pes
fs = 14
lw = 2.2
alpha = 0.7

for pes in ensemble_pes:
    plt.plot(pes[:, 0], pes[:, 1], '-', lw=lw, alpha=alpha)
plt.plot(pes_ave[:, 0], pes_ave[:, 1], 'k', lw=lw+1, label='Average model')

plt.xlim([np.min(dz_vals), np.max(dz_vals)])
plt.ylim(bottom=0.0)
plt.xlabel('Surface layer shift ($\\mathrm{\\AA}$)', fontsize=fs)
plt.ylabel('Potential Energy (eV)', fontsize=fs)
plt.gca().tick_params(labelsize=fs)

plt.legend(loc=9, fontsize=fs)
plt.tight_layout()
plt.savefig('ensemble_energy_surface.svg')

import numpy as np


def compute_energy_landscape(atoms, calc, dz_vals, tol=1e-3):
    """ Compute energy landscape for shifting the top surface layer """

    # compute reference energy
    atoms.set_calculator(calc)
    E0 = atoms.get_potential_energy()

    # find surface atom
    z = atoms.positions[:, 2]
    atom_indices = np.where(z + tol > np.max(z))

    # run displacement path
    data = []
    for dz in dz_vals:
        atoms_tmp = atoms.copy()
        atoms_tmp.set_calculator(calc)
        atoms_tmp.positions[atom_indices, 2] += dz
        data.append([dz, atoms_tmp.get_potential_energy() - E0])
    return np.array(data)