Suppose we want to calculate the minimum energy path of adatom diffusion on a (100) surface. We first need to choose an energy model, and in ASE, this is done by defining a "calculator". Let's choose our calculator to be Zhou's aluminum EAM potential, which we've used in previous labs.
We first import ASE's built-in EAM calculator class:
In [1]:
from ase.calculators.eam import EAM
Then set the potential file:
In [2]:
import os
pot_file = os.environ.get('LAMMPS_POTENTIALS') + '/Al_zhou.eam.alloy'
In [3]:
print(pot_file)
We just have to point ASE to the potential file. It automatically parses the file and constructs the corresponding EAM calculator for us:
In [4]:
zhou = EAM(potential=pot_file)
Next, we define our surface. Let's put the atom in the "hollow" site on the (100) surface. You can find out the adatom sites that are available for different types of surfaces in the ASE documentation: https://wiki.fysik.dtu.dk/ase/ase/build/surface.html
In [5]:
from ase.build import fcc100, add_adsorbate
slab = fcc100('Al', size=(3, 3, 3))
add_adsorbate(slab, 'Al', 2, 'hollow') # put adatom 2 A above the slab
slab.center(vacuum=5.0, axis=2) # 5 A of vacuum on either side
Let's see what the slab looks like:
In [6]:
from ase.visualize import view
view(slab, viewer='x3d')
Out[6]:
Let's set the calculator of the slab to the EAM calculator we defined above:
In [7]:
slab.set_calculator(zhou)
This lets us calculate the potential energy and forces on the atoms:
In [8]:
slab.get_potential_energy() # energy in eV
Out[8]:
In [9]:
slab.get_forces() # forces in eV/A
Out[9]:
Notice the nonzero forces, and in particular the strong force in the z-direction acting on the adatom. That's a signal that we're not relaxed.
We can use one of ASE's built-in structure optimizers to relax the structure and find the local energy minimum predicted by the EAM potential. The optimization terminates when the maximum force on an atom falls below fmax, which we set to 0.1 meV/A.
In [10]:
from ase.optimize import BFGS
dyn = BFGS(slab)
dyn.run(fmax=0.0001)
Out[10]:
To compute the activation barrier for adatom diffusion, we first need to know the endpoint of the transition path, which we can determine by looking at the atomic coordinates.
In [11]:
slab_2 = fcc100('Al', size=(3, 3, 3))
add_adsorbate(slab_2, 'Al', 2, 'hollow') # put adatom 2 A above the slab
slab_2.center(vacuum=5.0, axis=2) # 5 A of vacuum on either side
slab_2.set_calculator(EAM(potential=pot_file))
slab_2.positions
Out[11]:
In [43]:
slab_2.positions[-1][0:2] = slab_2.positions[10][0:2] # notice the adatom is directly above atom 9
view(slab_2, viewer='x3d')
Out[43]:
We again relax the structure:
In [44]:
dyn = BFGS(slab_2)
dyn.run(fmax=0.0001)
Out[44]:
We've succeeded in relaxing the endpoints. We can proceed to search for the saddle point using NEB.
In [47]:
from ase.neb import NEB
import numpy as np
In [48]:
# make band
no_images = 15
images = [slab]
images += [slab.copy() for i in range(no_images-2)]
images += [slab_2]
neb = NEB(images)
# interpolate middle images
neb.interpolate()
# set calculators of middle images
pot_dir = os.environ.get('LAMMPS_POTENTIALS')
for image in images[1:no_images-1]:
image.set_calculator(EAM(potential=pot_file))
# optimize the NEB trajectory
optimizer = BFGS(neb)
optimizer.run(fmax=0.01)
# calculate the potential energy of each image
pes = np.zeros(no_images)
pos = np.zeros((no_images, len(images[0]), 3))
for n, image in enumerate(images):
pes[n] = image.get_potential_energy()
pos[n] = image.positions
Let's plot the results.
In [51]:
import matplotlib.pyplot as plt
In [56]:
plt.plot(pes-pes[0], 'k.', markersize=10) # plot energy difference in eV w.r.t. first image
plt.plot(pes-pes[0], 'k--', markersize=10)
plt.xlabel('image #')
plt.ylabel('energy difference (eV)')
plt.show()
In this lab, we'll use the Bayesian ML code FLARE ("fast learning of atomistic rare events") that has recently been made open source. To set it up on Google cloud, vbox, or your personal machine, you'll need to pull it from github. I'll give the commands for setting everything up on an AP275 Google cloud instance, but the steps will be pretty much the same on any machine.
git clone https://github.com/mir-group/flare.git
The code is written in Python, but inner loops are accelerated with Numba, which you'll need to install with pip.
sudo apt install python3-pip
pip3 install numba
You'll see warnings from Numba that a more recent version of "colorama" needs to be installed. You can install it with pip:
pip3 install colorama
You may find it helpful to add the FLARE directory to your Python path and bash environment, which makes it easier to directly import files.
nano .profile
export FLARE=\$HOME/Software/flare
export PYTHONPATH=\$PYTHONPATH:\\$FLARE:\$FLARE/otf_engine:\\$FLARE/modules
source .profile
Let's look at a simple example to get a feel for how the code works. Let's put two atoms in a box along the x axis:
In [1]:
import numpy as np
from ase import Atoms
In [2]:
positions = np.array([[0, 0, 0], [1, 0, 0]])
cell = np.eye(3) * 10
two_atoms = Atoms(positions=positions, cell=cell)
In [3]:
from ase.visualize import view
view(two_atoms, viewer='x3d')
Out[3]:
Let's put equal and opposite forces on our atoms.
In [4]:
forces = np.array([[-1, 0, 0], [1, 0, 0]])
The FLARE code uses Gaussian process regression (GPR) to construct a covariant force field based on atomic forces, which are the training labels. GPR is a kernel based machine learning method, which means that it makes predictions by comparing test points to structures in the training set. For this simple system, we choose a two-body kernel, which compares pairwise distances in two structures.
In [5]:
from kernels import two_body, two_body_grad, two_body_force_en
from gp import GaussianProcess
The GP models are local, and require you to choose a cutoff. We'll pick 4 A. The kernel has a few hyperparameters which control uncertainty estimates and length scales. They can be optimized in a rigorous way by maximizing the likelihood of the training data, but for this lab we'll just set the hyperparameters to reasonable values.
In [6]:
hyps = np.array([1, 1, 1e-3]) # signal std, length scale, noise std
cutoffs = np.array([2])
gp_model = GaussianProcess(kernel=two_body, kernel_grad=two_body_grad, hyps=hyps,
cutoffs=cutoffs, energy_force_kernel=two_body_force_en)
The GP models take structure objects as input, which contain information about the cell and atomic coordinates (much like the Atoms class in ASE).
In [7]:
import struc
In [8]:
training_struc = struc.Structure(cell=cell, species=['A']*2, positions=positions)
We train the GP model by giving it training structures and corresponding forces:
In [9]:
gp_model.update_db(training_struc, forces)
gp_model.set_L_alpha()
As a quick check, let's make sure we get reasonable results on the training structure:
In [10]:
gp_model.predict(gp_model.training_data[0], 2) # second argument is the force component (x=1, y=2, z=3)
Out[10]:
To make it easier to get force and energy estimates from the GP model, we can wrap it in an ASE calculator:
In [11]:
from gp_calculator import GPCalculator
In [12]:
gp_calc = GPCalculator(gp_model)
Now let's test the covariance property of the model. Let's rotate the structure by 90 degrees, and see what forces we get.
In [13]:
# print positions, energy, and forces before rotation
two_atoms.set_calculator(gp_calc)
print(two_atoms.positions)
print(two_atoms.get_potential_energy())
print(two_atoms.get_forces())
In [15]:
two_atoms.rotate(90, 'z') # rotate the atoms 90 degrees about the z axis
two_atoms.set_calculator(gp_calc) # set calculator to gp model
In [16]:
# print positions, energy, and forces after rotation
print(two_atoms.positions)
print(two_atoms.get_potential_energy())
print(two_atoms.get_forces())
In the lab, we'll add a three-body term to the potential, which makes the model significantly more accurate for certain systems. Let's see how this would work for our (100) slab.
In [1]:
# initialize gp model
import kernels
import gp
import numpy as np
kernel = kernels.two_plus_three_body
kernel_grad = kernels.two_plus_three_body_grad
hyps = np.array([1, 1, 0.1, 1, 1e-3]) # sig2, ls2, sig3, ls3, noise std
cutoffs = np.array([4.96, 4.96])
energy_force_kernel = kernels.two_plus_three_force_en
gp_model = gp.GaussianProcess(kernel, kernel_grad, hyps, cutoffs,
energy_force_kernel=energy_force_kernel)
In [2]:
# make slab structure in ASE
from ase.build import fcc100, add_adsorbate
import os
from ase.calculators.eam import EAM
slab = fcc100('Al', size=(3, 3, 3))
add_adsorbate(slab, 'Al', 2, 'hollow') # put adatom 2 A above the slab
slab.center(vacuum=5.0, axis=2) # 5 A of vacuum on either side
pot_file = os.environ.get('LAMMPS_POTENTIALS') + '/Al_zhou.eam.alloy'
zhou = EAM(potential=pot_file)
slab.set_calculator(zhou)
In [3]:
# make training structure
import struc
training_struc = struc.Structure(cell=slab.cell,
species=['Al']*len(slab),
positions=slab.positions)
training_forces = slab.get_forces()
In [4]:
# add atoms to training database
gp_model.update_db(training_struc, training_forces)
gp_model.set_L_alpha()
In [5]:
# wrap in ASE calculator
from gp_calculator import GPCalculator
gp_calc = GPCalculator(gp_model)
In [6]:
# test on training structure
slab.set_calculator(gp_calc)
GP_forces = slab.get_forces()
In [8]:
# check accuracy by making a parity plot
import matplotlib.pyplot as plt
plt.plot(training_forces.reshape(-1), GP_forces.reshape(-1), '.')
plt.show()
In [ ]: