\maketitle \tableofcontents
This document contains supporting data for the work, “Neural Network and ReaxFF comparison for Au properties.” Much of the code utilized to produce this work is recorded here for reproducibility and transparency. Much of the molecular simulation code is produced using tools from the Atomic Simulation Environment (ASE) python modules: https://wiki.fysik.dtu.dk/ase/. This includes all of the databasing techniques, which are conducted through ASEs database module: https://wiki.fysik.dtu.dk/ase/ase/db/db.html#module-ase.db.
Due to the large number of calculations, data utilized in this work is stored in SQLite format (.db files) for ease of use. These data files are embedded into the PDF and can be accessed through the following links: Calculations specific to the 6 atom MD simulations performed in the paper: \attachfile{md-6atom.db}{(double-click to open)}. All other results: \attachfile{data.db}{(double-click to open)}. Examples of using these data files are provided in the following sections.
Large portions of the supporting information file function as demonstrations of how the data can be accessed and manipulated from these database files. For ease of search-ability the supporting information file is organized into sections which reflect the organization of the manuscript. That way, references from the manuscript can be directly referenced in the corresponding section of the supporting information file. The code used to generate the figures in each section of the manuscript are also included.
The supporting information document was prepared in org-mode (http://orgmode.org) syntax, which was subsequently exported to LaTeX and converted to a PDF. Briefly, org-mode is a plain editing format that enables intermingling of text, code and figures, with markup for typical document elements such as headings, links, tables, etc…, and arbitrary inclusion of LaTeX for equations. With the Emacs editor, the code in an org-mode document can be executed in place, and the output captured in the document. For example, tables can be generated by code, or the code for generating a figure can be embedded in the document. The data in tables can be used as input to other code blocks in the document as well. Org-mode enables selective export of the content to various formats including html, LaTeX, and PDF. These features, and many others, make org-mode a convenient platform for reproducible research, where all of the steps leading to conclusions drawn in the work can be documented in one place, but where it may be desirable not to show all the details in every view. For example, in an exported manuscript where code should usually not be visible, or supporting information document such as this one where it is desirable to see the code. Nevertheless, it may still be valuable to go back to the source, for example, to figure out how some analysis was done, especially if all the code is not exported.
The org-mode source for this document can be found here: \attachfile{supporting-information.org}.
All DFT calculations were performed using VASP (https://www.vasp.at). Due to the large amount of data produced from each individual VASP calculation, the results from each calculation are summarized in the databases discussed in the introduction. The following section demonstrates how these database can be searched for easy access to each of the 9972 calculation referred to in the manuscript. The DFT directories themselves can also be found on Figshare: http://dx.doi.org/10.6084/m9.figshare.1603022.
The VASP calculation directories are too large to be directly incorporated into the supporting information file, but critical information from these files can be organized into an ASE database and embedded into this PDF (See introduction for how to access these files). The ASE database is an extremely valuable tool and others may find it useful to know how we implement it into this work. For that reason the Python code used to generate the database is shown below for demonstrative purposes. Of course, this code is dependent upon files local to the Kitchin group, and can not be directly utilized unlike the other samples of code provided in this document.
import os, sys
from ase.io.trajectory import Trajectory
# Functions produced by Jacob Boes for work in computational catalysis
# These are freely available at: https://github.com/jboes/jbtools.git
import jbtools.gilgamesh as jb
# JASP is a wrapper used to streamline VASP calculation in VASP.
# More information can be found: https://github.com/jkitchin/jasp.git
from jasp import *
JASPRC['restart_unconverged'] = False
# We want to collect enegies and postions for all ionic steps in the
# geometry minimization. So first we compile a list of out.traj files.
for r, d, f in os.walk('DFT/'): # The location of the DFT calculations
if 'OUTCAR' in f and 'XDATCAR' in f:
try:
with jasp(r) as calc:
# Create a trajectory file
jb.compile_trajectory(calc)
# Access the created file
traj = Trajectory('out.traj')
n = len(traj)
# Add the contents to the database
for i, atoms in enumerate(traj):
jb.makedb(calc,
atoms=atoms,
dbname='~/research/neural/data.db',
keys={'traj': int(n-i-1)})
except:
print(r)
print(sys.exc_info()[0])
print('')Access to the VASP directories used to generate the database incorporated in this work is also available through a publicly hosted repository on the FigShare website (http://figshare.com/). The repository can be found here: http://dx.doi.org/10.6084/m9.figshare.1603022.
We organize information from the VASP calculations using descriptive keywords, which define the purpose of each calculation. A full list of the keywords used to categorize each calculation can be found in the results produced from the code block below.
from ase.db import connect
# Connect to the ASE database
db = connect('data.db')
# Selects all calculations in the database
data = db.select(['traj=0'])
# Create a dictionary of all the keys
keys, cnt = {}, 0
for entry in data:
cnt += 1
for k, v in entry.key_value_pairs.iteritems():
# Add all possible values to each dictionary
if k in keys:
keys[k] += [v]
else:
keys[k] = [v]
# Iterate through each key and print the values
print('{0:15s} {1:15s} {2} calculations total'.format('keyword', 'value', cnt))
print('---------------------------------------------------------')
for k, v in keys.iteritems():
vals = list(set(v))
# Only print the first 5 values
if len(vals) <= 5:
val = ", ".join(str(e) for e in vals)
print('{0:15s}: {1}'.format(k, val))
else:
val = ", ".join(str(e)[:5] for e in vals[:5])
print('{0:15s}: {1}, etc...'.format(k, val))
NOTE: NEB calculators do not store INCAR, KPOINT, or POTCAR parameters, as these files are kept in the parent directory. Keys generated from these files are place holders, and thus are not correct for NEB calculations. i.e. there are no PW91 calculations, the value was simply not recorded.
The database is broken into three main categories: bulk, surface, and cluster calculations. These categories are organized by the correspondingly named key. For example, to find all bulk calculations, one would search for the ‘bulk’ key:
from ase.db import connect
db = connect('data.db')
# Now we select only bulk calculations
data = db.select(['bulk'])
# 'bulk' = all bulk calculations
# 'surf' = all surface calculations
# 'cluster' = all cluster calculations
# And again, we print all the possbile keys
keys, cnt = {}, 0
for entry in data:
cnt += 1
for k, v in entry.key_value_pairs.iteritems():
if k in keys:
keys[k] += [v]
else:
keys[k] = [v]
print('{0:15s} {1:15s} {2} calculations total'.format('keyword', 'value', cnt))
print('---------------------------------------------------------')
for k, v in keys.iteritems():
vals = list(set(v))
if len(vals) <= 5:
val = ", ".join(str(e) for e in vals)
print('{0:15s}: {1}'.format(k, val))
else:
val = ", ".join(str(e)[:5] for e in vals[:5])
print('{0:15s}: {1}, etc...'.format(k, val))
This can also be done for surface and cluster calculations as noted in the comments above. Information about specific data sets can be found in the following section.
The reactive force field produced in this work utilizes 3-body interaction as shown in Table ref:tbl-3body. The ReaxFF itself can be opened using the following attachment: \attachfile{ffield.reax.mitch}{(double-click to open)}.
| Valence Angle | θ$o$ (degrees) | k$a$ (kcal/mol) | k$b$ (radians$2$) | p$v,1$ | p$v,2$ |
|---|---|---|---|---|---|
| Au-Au-Au | 12.2362 | 6.788 | 0.1388 | 0 | 0.9168 |
***Comparing Force Fields with and without 3-Body Terms We parameterized force fields with and without 3-body terms using identical training sets to test the effect of incorporating 3-body terms in Au force fields. Our preliminary results in ref:si-2B3BEOS show that the curvature, optimal volume, and cohesive energy of the Au fcc, sc, and diamond EOS all significantly improve by adding 3-body terms, which appear to shift the onset of the convex regions to larger volumes.
We note that this treatment also brings significantly higher costs to the calculations ref:si-Reax3Bcost. While the increase in computational cost due to 3-body terms is minimal in small systems, large systems prove to be much more expensive.
MCFFopt is a stochastic force field optimization process that relies on minimizing an objective function (Equation ref:eq-1).
\begin{equation} Total Error = ∑i=1n\left[\frac{EFField - EQM}{Weight}\right] \label{eq-1} \end{equation}
The optimizer program does this by randomly changing force field parameters within a range defined by the user and recalculating the objective function. Any parameter change that decreases the total error is instantly accepted, while moves that increase the total error have a probability of being accepted determined by Equation ref:eq-2.
\begin{equation} Probability = min[1,exp(-β Δ Error)] \label{eq-2} \end{equation}
Where
The MCFFopt tool does not have a built in force field parameter convergence criteria. However, over the course of an optimization calculation, the degree that parameters can change decreases. This means that after a certain number of steps parameters will begin to change by insignificant amounts, and the total error will stagnate. Thus, after a single run with the MCFFopt tool parameters will appear converged, but we can achieve still lower total errors, as shown in ref:tbl-mcffopt, by simply restarting the MCFFopt process until differences in total errors between runs becomes small.
Restarting the MCFF optimization procedure resets these parameters, and the process will go through another annealing phase, and then the total error will continue to decrease. We stopped optimizations once the final error was less than one percent of the initial error.
The training process requires user-defined parameters to control the MCFF optimization (contained in table ref:tbl-mcffopt), and a range of maximum and minimum values for each optimizable ReaxFF parameter. More information on these specific files is available in the ADF documentation.
| Parameter | Value |
|---|---|
| ($βinitial$) (mcbeta) | 0.0100 |
| Increase |
1.0001 |
| Fraction of variables active per step (mcacof) | 0.2 |
| Target acceptance rate (mctart) | 0.30 |
| Max. acceptance rate (mcmart) | 0.70 |
| Divide parameter range into y steps (mcrxxd) | 100 |
ReaxFF training requires files that contain the geometries and energies that the force field will be trained to, weights of the relative importance of those geometries and energies, and an initial set of force field parameters taken from Ref. citenum:keith-2010-react. That ReaxFF did not contain 3-body interaction parameters, so we used an average of 3-body interaction terms available for other atom types as our initial guess.
We initially used a weighting scheme nearly equivalent to that used in Ref. citenum:keith-2010-react. We increased and/or decreased the weights of different geometry types (bulk vs. cluster vs. surface) until the force field adequately reproduced all of the data in the training set. Briefly, EOS structures were weighted on scale sliding from 0.2 to 10. Geometries near the minimum energy bulk structures were weighted highly (0.2), while structures further away from the equilibrium geometries were weighted lower (10.0). All surface and cluster calculations had weights of 1.0. In principle, changing these weights can result in force fields that are better suited for different geometry types.
ReaxFF is a module included in the LAMMPS simulation suit (http://lammps.sandia.gov), and must be run inside of this framework. However, a wrapper has been previously developed which allows LAMMPS calculations to be started and managed from a python interface. A python function was developed to rapidly calculate the energy of an ASE atoms object. The code for this function is demonstrated below.
# Personal function for interaction with LAMMPS in python
# Can be found at: https://github.com/jboes/jbtools
import jbtools.utils as jb
from ase.db import connect
db = connect('data.db')
# Select all entries
for d in db.select():
atoms = db.get_atoms(d.id)
Rnrg = jb.reax_potential_energy(atoms)
db.update(d.id, reax_energy=Rnrg)With this function, the energies of all of the structures in the database could then be rapidly calculated and incorporated into the database for ease of access in the analysis portion of this work. This way, ReaxFF predicted energies for all structures included in the database can be accessed through the keyword “reax_energy”.
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from ase.db import connect
import numpy as np
from scipy.stats import norm
from matplotlib import gridspec
import matplotlib.patches as mpatches
db = connect('data.db')
S, Re, Qe = [], [], []
for d in db.select(['train_set=True']):
S += [d.structure]
Qe += [d.energy / d.natoms]
Re += [d.reax_energy / d.natoms]
S = np.array(S)
Qe = np.array(Qe)
Re = np.array(Re)
cmap, hdl = {}, []
for s, c in zip(set(S), ['b', 'r', 'g']):
cmap[s] = c
hdl += [mpatches.Patch(color=c, label=s)]
c = []
for s in S:
c += [cmap[s]]
RMSE = np.sqrt(sum((Re - Qe) ** 2) / len(Re - Qe))
(mu, sigma) = norm.fit(Re - Qe)
fig = plt.figure(figsize=(6, 4))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1])
ax0 = plt.subplot(gs[0])
ax0.plot([min(Qe), 0], [0, 0], 'k--', lw=2)
ax0.scatter(Qe, Re - Qe, marker='o', color=c)
ax0.text(min(Qe) + 0.15, 1.08,
'RMSE: {0:1.3f}'.format(RMSE),
fontsize=12, va='top', ha='left')
ax0.set_xlim(min(Qe), 0)
ax0.set_ylim(-1.2, 1.2)
ax0.set_xlabel('DFT potential energy (eV/atom)')
ax0.set_ylabel('Residual error (eV/atom)')
ax0.legend(loc=3, handles=hdl, fontsize=12, frameon=False)
ax1 = plt.subplot(gs[1])
n, bins, patches = ax1.hist(Re - Qe, 50,
range=(-1.2, 1.2),
weights=np.ones_like(Re - Qe)/len(Re),
facecolor='k',
alpha=0.5,
orientation='horizontal')
y = mlab.normpdf(bins, mu, sigma)
ax1.text(0.05, 1.12, '$\mu$: {0:1.3f}'.format(mu), fontsize=12,
va='top', ha='left')
ax1.text(0.05, 0.96, '$\sigma$: {0:1.3f}'.format(sigma), fontsize=12,
va='top', ha='left')
ax1.plot(y / sum(y), bins, 'k--', lw=2)
ax1.plot([0, 10], [0, 0], 'k--', lw=2)
ax1.set_xlabel('Probability')
ax1.set_ylim(-1.2, 1.2)
ax1.set_xlim(0, 0.6)
ax1.set_yticklabels([])
ax1.set_xticks(ax1.get_xticks()[1::2])
plt.tight_layout(w_pad=-0.5)
for ext in ['png', 'eps']:
plt.savefig('./images/fig-reax-train.{0}'.format(ext), dpi=300)import matplotlib.pyplot as plt
from ase.db import connect
import numpy as np
from matplotlib import gridspec
import matplotlib.patches as mpatches
db = connect('data.db')
S, Re, Qe = [], [], []
for d in db.select(['train_set=False']):
S += [d.structure]
Qe += [d.energy / d.natoms]
Re += [d.reax_energy / d.natoms]
S = np.array(S)
Qe = np.array(Qe)
Re = np.array(Re)
cmap, hdl = {}, []
for s, c in zip(set(S), ['b', 'r', 'g']):
cmap[s] = c
hdl += [mpatches.Patch(color=c, label=s)]
c = []
for s in S:
c += [cmap[s]]
RMSE = np.sqrt(sum((Re - Qe) ** 2) / len(Re - Qe))
fig = plt.figure(figsize=(6, 4))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1])
ax0 = plt.subplot(gs[0])
ax0.plot([min(Qe), 0], [0, 0], 'k--', lw=2)
ax0.scatter(Qe, Re - Qe, marker='o', color=c)
ax0.text(min(Qe) + 0.15, 0.7,
'RMSE: {0:1.3f}'.format(RMSE),
fontsize=12, va='top', ha='left')
ax0.set_xlim(min(Qe), 0)
ax0.set_ylim(-0.8, 0.8)
ax0.set_xlabel('DFT potential energy (eV/atom)')
ax0.set_ylabel('Residual error (eV/atom)')
ax0.legend(loc=1, handles=hdl, fontsize=12, frameon=False)
ax1 = plt.subplot(gs[1])
n, bins, patches = ax1.hist(Re - Qe, 50,
range=(-1.2, 1.2),
weights=np.ones_like(Re - Qe)/len(Re),
facecolor='k',
alpha=0.5,
orientation='horizontal')
ax1.plot([0, 10], [0, 0], 'k--', lw=2)
ax1.set_xlabel('Probability')
ax1.set_ylim(-1.2, 1.2)
ax1.set_xlim(0, 0.6)
ax1.set_yticklabels([])
ax1.set_xticks(ax1.get_xticks()[1::2])
plt.tight_layout(w_pad=-0.5)
for ext in ['png', 'eps']:
plt.savefig('./images/fig-reax-valid.{0}'.format(ext), dpi=300)BPNNs were produced using the Neural code developed by the Peterson group at Brown University (https://bitbucket.org/andrewpeterson/neural). This code is no longer supported since the time the work was completed. However, the Peterson group is currently developing a sister code called AMP which is capable of all the same functionality as Neural (https://bitbucket.org/andrewpeterson/amp). The parameter files used by AMP are not backwards compatible with Neural, so the parameter files included here have been manually updated to function with AMP.
The parameters file needed to run the BPNN produced in this work is attached here: \attachfile{neural-parameters.json}{(double-click to open)}.
The code developed by the Peterson group is already integrated into ASE, making calculation of energy using the BPNN developed in this work trivial. The code used to calculate these energies and add them to the database is included below.
from ase.db import connect
from amp import Amp
db = connect('data.db')
# Establish the BPNN as the calculator for our energies
calc = Amp(load='neural-parameters.json')
# Select all entries
for d in db.select():
atoms = db.get_atoms(d.id)
atoms.set_calculator(calc)
Nnrg = atoms.get_potential_energy()
db.update(d.id, neural_energy=Nnrg)This allows the BPNN predicted energy of a structure to be accessed easily by first defining the structure of interest and then utilizing the keyword “neural_energy”.
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from ase.db import connect
import numpy as np
from scipy.stats import norm
from matplotlib import gridspec
import matplotlib.patches as mpatches
db = connect('data.db')
S, Qe, Ne = [], [], []
for d in db.select(['train_set=True']):
S += [d.structure]
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
S = np.array(S)
Qe = np.array(Qe)
Ne = np.array(Ne)
cmap, hdl = {}, []
for s, c in zip(set(S), ['b', 'r', 'g']):
cmap[s] = c
hdl += [mpatches.Patch(color=c, label=s)]
c = []
for s in S:
c += [cmap[s]]
RMSE = np.sqrt(sum((Ne - Qe) ** 2) / len(Ne - Qe))
(mu, sigma) = norm.fit(Ne - Qe)
fig = plt.figure(figsize=(6, 4))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1])
ax0 = plt.subplot(gs[0])
ax0.plot([min(Qe), 0], [0, 0], 'k--', lw=2)
ax0.scatter(Qe, Ne - Qe, marker='o', color=c)
ax0.text(min(Qe) + 0.1, 0.14,
'RMSE: {0:1.3f}'.format(RMSE),
fontsize=12, va='top', ha='left')
ax0.set_xlim(min(Qe), 0)
ax0.set_ylim(-0.15, 0.15)
ax0.set_xlabel('DFT potential energy (eV/atom)')
ax0.set_ylabel('Residual error (eV/atom)')
ax0.legend(loc='best', handles=hdl, fontsize=12, frameon=False)
ax1 = plt.subplot(gs[1])
n, bins, patches = ax1.hist(Ne - Qe, 50,
range=(-0.15, 0.15),
weights=np.ones_like(Ne - Qe)/len(Ne),
facecolor='k',
alpha=0.5,
orientation='horizontal')
y = mlab.normpdf(bins, mu, sigma)
ax1.text(0.05, 0.142, '$\mu$: {0:1.3f}'.format(mu), fontsize=12, va='top', ha='left')
ax1.text(0.05, 0.122, '$\sigma$: {0:1.3f}'.format(sigma), fontsize=12, va='top', ha='left')
ax1.plot(y / sum(y), bins, 'k--', lw=2)
ax1.plot([0, 50], [0, 0], 'k--', lw=2)
ax1.set_xlabel('Probability')
ax1.set_ylim(-0.15, 0.15)
ax1.set_xlim(0, 0.6)
ax1.set_yticklabels([])
ax1.set_xticks(ax1.get_xticks()[1::2])
plt.tight_layout(w_pad=-0.5)
for ext in ['png', 'eps']:
plt.savefig('./images/fig-neural-train.{0}'.format(ext), dpi=300)import matplotlib.pyplot as plt
from ase.db import connect
import numpy as np
from matplotlib import gridspec
import matplotlib.patches as mpatches
db = connect('data.db')
S, Qe, Ne = [], [], []
for d in db.select(['train_set=False']):
S += [d.structure]
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
S = np.array(S)
Qe = np.array(Qe)
Ne = np.array(Ne)
cmap, hdl = {}, []
for s, c in zip(set(S), ['b', 'r', 'g']):
cmap[s] = c
hdl += [mpatches.Patch(color=c, label=s)]
c = []
for s in S:
c += [cmap[s]]
RMSE = np.sqrt(sum((Ne - Qe) ** 2) / len(Ne - Qe))
fig = plt.figure(figsize=(6, 4))
gs = gridspec.GridSpec(1, 2, width_ratios=[2, 1])
ax0 = plt.subplot(gs[0])
ax0.plot([min(Qe), 0], [0, 0], 'k--', lw=2)
ax0.scatter(Qe, Ne - Qe, marker='o', color=c)
ax0.text(min(Qe) + 0.1, 0.14,
'RMSE: {0:1.3f}'.format(RMSE),
fontsize=12, va='top', ha='left')
ax0.set_xlim(min(Qe), 0)
ax0.set_ylim(-0.15, 0.15)
ax0.set_xlabel('DFT potential energy (eV/atom)')
ax0.set_ylabel('Residual error (eV/atom)')
ax0.legend(loc='best', handles=hdl, fontsize=12, frameon=False)
ax1 = plt.subplot(gs[1])
n, bins, patches = ax1.hist(Ne - Qe, 50,
range=(-0.15, 0.15),
weights=np.ones_like(Ne - Qe)/len(Ne),
facecolor='k',
alpha=0.5,
orientation='horizontal')
ax1.plot([0, 50], [0, 0], 'k--', lw=2)
ax1.set_xlabel('Probability')
ax1.set_ylim(-0.15, 0.15)
ax1.set_xlim(0, 0.6)
ax1.set_yticklabels([])
ax1.set_xticks(ax1.get_xticks()[1::2])
plt.tight_layout(w_pad=-0.5)
for ext in ['png', 'eps']:
plt.savefig('./images/fig-neural-valid.{0}'.format(ext), dpi=300)import matplotlib.pyplot as plt
from ase.utils.eos import EquationOfState
from ase.db import connect
import numpy as np
from ase.units import kJ
db = connect('data.db')
print('#+caption: Comparison of EOS metrics for DFT, ReaxFF, and NPNN fits as shown in Figure ref:fig-bulk-eos.')
print('#+attr_latex: :placement [H]')
print('#+tblname: tbl-eos')
print('|Structure|Minimum volume (\AA^{3})|Minimum energy (eV)|Bulk Mod. (GPa)|')
print('|-')
f, ax = plt.subplots(1, 3, figsize=(6, 5))
tag = ['Face Centered\nCubic', 'Simple Cubic', 'Diamond']
for i, key in enumerate(['fcc', 'sc', 'diam']):
V, Qe, Re, Ne = [], [], [], []
for d in db.select(['bulk={0}'.format(key), 'factor']):
V += [d.volume / d.natoms]
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
sel = V[Qe.index(min(Qe))]
ind = (np.array(V) > sel - 15) & (np.array(V) < sel + 15)
x = np.linspace(min(V), max(V), 250)
V = np.array(V)[ind]
for nrg, name, col in zip([Qe, Ne, Re],
['DFT', 'BPNN', 'ReaxFF'],
['k-', 'r--', 'b:']):
nrg = np.array(nrg)[ind]
eos = EquationOfState(V, nrg)
v0, e0, B = eos.fit()
fit = np.poly1d(np.polyfit(V**-(1.0 / 3), nrg, 3))
ax[i].plot(x, fit(x**-(1.0 / 3)), col, lw=2, label='{0}'.format(name))
ax[i].set_xlim(min(V), max(V))
ax[i].set_ylim(-3.5, -1.0)
ax[i].set_title('{0}'.format(tag[i]))
if i > 0:
ax[i].set_yticklabels([])
print('|{0}-{1}|{2:1.2f}|{3:1.2f}|{4:1.0f}'.format(name, key, v0, e0,
B / kJ * 1.0e24))
print('|-')
ax[0].set_xticks([14, 19, 24, 29])
ax[1].set_xticks([17, 22, 27, 32])
ax[2].set_xticks([21, 26, 31, 36, 41])
ax[0].set_ylabel('Potential energy (eV/atom)')
ax[1].set_xlabel('Volume ($\AA$/atom)')
ax[2].legend(loc='best', fontsize=12)
plt.tight_layout(w_pad=-0.3)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-bulk-eos.{0}'.format(ext), dpi=300)| Structure | Minimum volume (Å3) | Minimum energy (eV) | Bulk Mod. (GPa) |
|---|---|---|---|
| DFT-fcc | 17.97 | -3.23 | 147 |
| BPNN-fcc | 17.99 | -3.23 | 145 |
| ReaxFF-fcc | 17.60 | -3.22 | 122 |
| DFT-sc | 20.73 | -3.02 | 110 |
| BPNN-sc | 20.66 | -3.02 | 110 |
| ReaxFF-sc | 21.29 | -2.96 | 84 |
| DFT-diam | 29.04 | -2.51 | 56 |
| BPNN-diam | 28.98 | -2.51 | 57 |
| ReaxFF-diam | 31.92 | -2.54 | 37 |
To conserve space in the manuscript, the EOS for hcp and bcc structures are shown in Figure ref:si-bulk-eos. The corresponding data resulting from the fits to the EOS can be found in Table ref:tbl-eos2
import matplotlib.pyplot as plt
from ase.utils.eos import EquationOfState
from ase.db import connect
import numpy as np
from ase.units import kJ
db = connect('data.db')
print('#+caption: Comparison of EOS metrics for DFT, ReaxFF, and NPNN fits as shown in Figure ref:si-bulk-eos.')
print('#+attr_latex: :placement [H]')
print('#+tblname: tbl-eos2')
print('|Structure|Minimum volume (\AA^{3})|Minimum energy (eV)|Bulk Mod. (GPa)|')
print('|-')
f, ax = plt.subplots(1, 2, figsize=(6, 5))
tag = ['Body Centered\nCubic', 'Hexagonal Close\nPacking']
for i, key in enumerate(['bcc', 'hcp']):
V, Qe, Re, Ne = [], [], [], []
for d in db.select(['bulk={0}'.format(key), 'factor']):
V += [d.volume / d.natoms]
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
sel = V[Qe.index(min(Qe))]
ind = (np.array(V) > sel - 15) & (np.array(V) < sel + 15)
x = np.linspace(min(V), max(V), 250)
V = np.array(V)[ind]
for nrg, name, col in zip([Qe, Ne, Re],
['DFT', 'BPNN', 'ReaxFF'],
['k-', 'r--', 'b:']):
nrg = np.array(nrg)[ind]
eos = EquationOfState(V, nrg)
v0, e0, B = eos.fit()
fit = np.poly1d(np.polyfit(V**-(1.0 / 3), nrg, 3))
ax[i].plot(x, fit(x**-(1.0 / 3)), col, lw=2, label='{0}'.format(name))
ax[i].set_xlim(min(V), max(V))
ax[i].set_ylim(-3.5, -1.0)
ax[i].set_title('{0}'.format(tag[i]))
if i > 0:
ax[i].set_yticklabels([])
print('|{0}-{1}|{2:1.2f}|{3:1.2f}|{4:1.0f}'.format(name, key, v0, e0,
B / kJ * 1.0e24))
print('|-')
ax[0].set_xticks([14, 19, 24, 29])
ax[1].set_xticks([14, 19, 24, 29])
ax[0].set_ylabel('Potential energy (eV/atom)')
ax[0].set_xlabel('Volume ($\AA$/atom)')
ax[1].set_xlabel('Volume ($\AA$/atom)')
ax[1].legend(loc='best', fontsize=12)
plt.tight_layout(w_pad=-0.3)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-bulk-eos.{0}'.format(ext), dpi=300)| Structure | Minimum volume (Å3) | Minimum energy (eV) | Bulk Mod. (GPa) |
|---|---|---|---|
| DFT-bcc | 18.02 | -3.21 | 145 |
| BPNN-bcc | 18.00 | -3.21 | 146 |
| ReaxFF-bcc | 18.36 | -3.11 | 107 |
| DFT-hcp | 17.98 | -3.23 | 147 |
| BPNN-hcp | 17.93 | -3.23 | 148 |
| ReaxFF-hcp | 17.60 | -3.22 | 122 |
import matplotlib.pyplot as plt
from ase.utils.eos import EquationOfState
from ase.db import connect
import numpy as np
from ase.units import kJ
db = connect('data.db')
f, ax = plt.subplots(1, 5, figsize=(12, 5))
tag = ['Face Centered\nCubic', 'Body Centered\nCubic',
'Hexagonal Close\nPacking', 'Simple Cubic', 'Diamond',]
for i, key in enumerate(['fcc', 'bcc', 'hcp', 'sc', 'diam',]):
V, Qe, Re, Ne = [], [], [], []
for d in db.select(['bulk={0}'.format(key), 'factor']):
V += [d.volume / d.natoms]
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
srt = [j[0] for j in sorted(enumerate(V), key=lambda x:x[1])]
V = np.array(V)[srt]
Qe = np.array(Qe)[srt]
Ne = np.array(Ne)[srt]
Re = np.array(Re)[srt]
ax[i].plot(V, Qe, 'k-', lw=2, label='DFT')
ax[i].plot(V, Ne, 'r--', lw=2, label='BPNN')
ax[i].plot(V, Re, 'b:', lw=2, label='ReaxFF')
if i > 0:
ax[i].set_yticklabels([])
ax[i].set_ylim(-3.5, 0.2)
ax[i].set_xlim(0, 200)
ax[i].set_title('{0}'.format(tag[i]))
ax[0].set_ylabel('Potential energy (eV/atom)')
ax[2].set_xlabel('Volume ($\AA$/atom)')
ax[4].legend(loc='best', fontsize=12)
plt.tight_layout(w_pad=-1.3)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-full-eos.{0}'.format(ext), dpi=300)from ase.db import connect
import numpy as np
import matplotlib.pyplot as plt
db = connect('data.db')
# Energy of single atom fcc reference
ref = db.get(['bulk=fcc', 'traj=0', 'factor=1'])
C, Qe, Re, Ne = [], [], [], []
for d in db.select(['bulk=fcc', 'type=vacancy', 'traj=0', 'image=0']):
C += [d.concentration]
Qe += [d.energy - (d.natoms * ref.energy)]
Ne += [d.neural_energy - (d.natoms * ref.neural_energy)]
Re += [d.reax_energy - (d.natoms * ref.reax_energy)]
plt.figure(figsize=(6, 4))
plt.plot([0, 0.13], [0.42, 0.42], 'k--')
plt.text(0.065, 0.42, 'Literature DFT', va='bottom', fontsize=14)
plt.plot([0, 0.13], [0.93, 0.93], 'k--')
plt.text(0.065, 0.93, 'Experimental', va='bottom', fontsize=14)
plt.text(0.125, Qe[0]+0.03, 'DFT', color='k', fontsize=14, ha='right')
plt.scatter(C, Qe, marker='o', color='k')
plt.text(0.125, Ne[0]+0.03, 'BPNN', color='r', fontsize=14, ha='right')
plt.scatter(C, Ne, marker='s', color='r')
plt.text(0.125, Re[0]+0.03, 'ReaxFF', color='b', fontsize=14, ha='right')
plt.scatter(C, Re, marker='^', color='b')
plt.xlabel('Vacancy concentration (vacancies/atom)')
plt.ylabel('Vacancy formation energy (eV/vacancy)')
plt.xlim(0, 0.13)
plt.ylim(0, 1.0)
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-vacancy-formation.{0}'.format(ext), dpi=300)As mentioned in the manuscript, the unit cell used to calculate the vacancy concentration at ≈ 0.015 vacancies/atom reconfigured to a different, less favorable, local minima. This is demonstrated in Figure ref:si-vacancy-reconfig which depicts the energies of the relaxation pathways for the DFT calculation. The reconfiguration is shown as the last structure in the trajectory along side the minimum energy structure.
from ase.db import connect
import numpy as np
import matplotlib.pyplot as plt
from ase.io import write
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
import os
db = connect('data.db')
trajectory = db.select(['bulk=fcc', 'type=vacancy', 'image=0',
'concentration=0.015625'])
E, t = [], []
for traj in trajectory:
E.append(traj.energy)
t.append(traj.traj)
E = np.array(E) - min(E)
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
for i, a in enumerate([3, 0]):
atoms = db.get_atoms(['bulk=fcc', 'type=vacancy',
'traj={0}'.format(a), 'image=0',
'concentration=0.015625'])
write('./images/temp.png', atoms, show_unit_cell=2)
image = mpimg.imread('./images/temp.png')
imagebox = OffsetImage(image)
ax.add_artist(AnnotationBbox(imagebox,
xy=(a, E[t[a]]),
xybox=(a + 4, E[t[a]] + 0.04*(i + 1)),
pad=-0.2,
frameon=False,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=-90,armA=0')
)
)
os.unlink('./images/temp.png')
ax.plot(t, E)
ax.invert_xaxis()
ax.set_xticklabels([])
plt.ylabel('Relative potential energy (eV)')
plt.xlabel('Reaction coordinate (a.u.)')
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-vacancy-reconfig.{0}'.format(ext), dpi=300)from ase.db import connect
from ase.visualize import view
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import fmin
db = connect('data.db')
ref = db.get(['bulk', 'NEB=True', 'image=0', 'type=vacancy', 'traj=0'])
I, Qe, Re, Ne = [], [], [], []
for d in db.select(['bulk', 'NEB=True', 'type=vacancy', 'traj=0']):
I += [d.image]
Qe += [d.energy - ref.energy]
Ne += [d.neural_energy - ref.neural_energy]
Re += [d.reax_energy - ref.reax_energy]
sort = [i[0] for i in sorted(enumerate(I), key=lambda x:x[1])]
I = np.array([I[i] for i in sort]) + 1
x = np.linspace(I.min(), I.max())
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(1,1,1)
for nrg, name, s in zip([Qe, Ne, Re],
['DFT', 'BPNN', 'ReaxFF'],
['ko-', 'rs--', 'b^:']):
nrg = np.array([nrg[i] for i in sort])
f = interp1d(I, -nrg, kind='cubic', bounds_error=False)
xmax = fmin(f, 3.0, disp=False)
ax.plot(I, nrg, s[:2], label=name)
ax.plot(x, -f(x), color=s[0], ls=s[2:])
ax.annotate('{0} $E^\ddag$= {1:1.2f}'.format(name, float(-f(xmax))),
xy=(3.0, -f(xmax)), xytext=(3.7, -f(xmax)),
ha='left', va='center', color=s[0],
arrowprops=dict(arrowstyle="->",
shrinkB=10,
color=s[0]))
ax.set_xticklabels([])
plt.xlabel('Reaction coordinate (a.u.)')
plt.ylabel('Potential energy (eV)')
plt.xlim(1, 5)
plt.ylim(0, 0.62)
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-vacancy-diffusion.{0}'.format(ext), dpi=300)from ase.db import connect
import matplotlib.pyplot as plt
import numpy as np
db = connect('data.db')
f, ax = plt.subplots(2, 2, figsize=(6, 5))
for i, k in enumerate([['single1', 'terrace', 12],
['dimer1', 'dimer', 6]]):
Qe, Re, Ne = [], [], []
for j, d in enumerate(db.select(['config={0}'.format(k[0]),
'group=timo'])):
Qe += [d.energy]
Ne += [d.neural_energy]
Re += [d.reax_energy]
Qe = np.array(Qe[:k[2]])
Ne = np.array(Ne[:k[2]])
Re = np.array(Re[:k[2]])
if i == 0:
Qe = np.hstack([Qe[::-1], Qe])
Ne = np.hstack([Ne[::-1], Ne])
Re = np.hstack([Re[::-1], Re])
m = list(Qe).index(Qe.min())
n = range(len(Qe))
Nerr = (Ne - Ne[m]) - (Qe - Qe[m])
Rerr = (Re - Re[m]) - (Qe - Qe[m])
# Plotting energy points
ax[0, i].plot(n, Ne - Ne[m], 'rs', label='BPNN')
ax[0, i].plot(n, Re - Re[m], 'b^', label='ReaxFF')
ax[0, i].plot(n, Qe - Qe[m], 'ko', fillstyle='none', label='DFT')
# Plotting the residuals
ax[1, i].plot([min(n), max(n)], [0, 0], 'k:')
ax[1, i].plot(n, Nerr, 'r.')
ax[1, i].plot(n, Rerr, 'b.')
ax[0, i].set_xlim(min(n), max(n))
ax[1, i].set_xlim(min(n), max(n))
ax[1, i].set_ylim(-0.3, 0.3)
ax[0, i].set_ylim(0, 1.0)
ax[0, i].set_title('{0} diffusion'.format(k[1]))
ax[1, i].set_xlabel('Reaction pathway (a.u.)')
ax[0, i].set_xticklabels([])
ax[1, i].set_xticklabels([])
if i != 0:
ax[0, i].set_yticklabels([])
ax[1, i].set_yticklabels([])
ax[1, 0].set_yticks(ax[1, 0].get_yticks()[:-2])
ax[1, 0].set_ylabel('Residual error (eV)')
ax[0, 0].set_ylabel('Total energy (eV)')
ax[0, 0].legend(loc='best', numpoints=1, fontsize=12, frameon=False)
plt.tight_layout(w_pad=0.2, h_pad=-0.3)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-full-diffusion.{0}'.format(ext), dpi=300)from ase.db import connect
import xlrd
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
import numpy as np
db = connect('data.db')
C = {}
for d in db.select(['group=timo', 'config!=sl']):
if d.config not in C.keys():
C[d.config] = {}
if d.identity not in C[d.config].keys():
C[d.config][d.identity] = []
C[d.config][d.identity] += [[d.energy, d.neural_energy, d.reax_energy]]
rQe, rNe, rRe = [], [], []
for cfg, data0 in C.iteritems():
for lbl, data1 in data0.iteritems():
Qe = np.array(data1).T[0]
m = list(Qe).index(Qe.min())
Qe = Qe - Qe[m]
Ne = np.array(data1).T[1] - np.array(data1).T[1][m]
Re = np.array(data1).T[2] - np.array(data1).T[2][m]
Qe = np.delete(Qe, m)
Ne = np.delete(Ne, m)
Re = np.delete(Re, m)
rQe += list(Qe)
rNe += list(Ne - Qe)
rRe += list(Re - Qe)
pct = 0.1
verts = [(0., pct), (1.2, pct), (1.2, -pct), (0., -pct), (0., 0.)]
codes = [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY]
path = Path(verts, codes)
Nc = path.contains_points(zip(rQe, rNe))
Ni = float(sum(Nc)) / len(Nc)
Rc = path.contains_points(zip(rQe, rRe))
Ri = float(sum(Rc)) / len(Rc)
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
patch = patches.PathPatch(path, facecolor='y', edgecolor='y', alpha=0.3)
ax.add_patch(patch)
ax.plot([0, 1.2], [0, 0], 'k--', zorder=1)
ax.scatter(rQe[:8] + rQe[30:], rNe[:8] + rNe[30:], color='r', marker='s', zorder=10, s=30)
ax.scatter(rQe[:8] + rQe[30:], rRe[:8] + rRe[30:], color='b', marker='^', zorder=20, s=30)
ax.scatter(rQe[8:30], rNe[8:30], color='r', marker='s', facecolor='none', zorder=30, s=20)
ax.scatter(rQe[8:30], rRe[8:30], color='b', marker='^', facecolor='none', zorder=40, s=30)
ax.text(0.02, 0.38, 'Within $\pm$ 0.1 eV error:', va='top', ha='left')
ax.text(0.02, 0.32,
'ReaxFF: {0:1.1f}%'.format(Ri*100, pct*100),
va='top', ha='left', color='b', zorder=100)
ax.text(0.02, 0.26,
'BPNN: {0:1.1f}%'.format(Ni*100, pct*100),
va='top', ha='left', color='r', zorder=100)
plt.xlabel('DFT Potential Energy (eV)')
plt.ylabel('Residual Error (eV)')
plt.xlim(0, max(rQe))
plt.ylim(-0.4, 0.4)
plt.tight_layout()
for ext in ['png', 'eps']:
plt.savefig('./images/fig-barrier-residuals.{0}'.format(ext), dpi=300)Not all of the fcc(100) surface diffusion pathways could be directly shown in the manuscript. Instead, we show them here to demonstrate how each of the potentials performs in each case. The residuals of each method are also shown. These residuals are the same values incorporated in Figure fig-barrier-residuals included in the manuscript.
from ase.db import connect
import matplotlib.pyplot as plt
import numpy as np
db = connect('data.db')
C = {}
for d in db.select(['group=timo', 'config!=sl']):
if d.config not in C.keys():
C[d.config] = {}
if d.identity not in C[d.config].keys():
C[d.config][d.identity] = []
C[d.config][d.identity] += [[d.energy, d.neural_energy, d.reax_energy]]
f, ax = plt.subplots(2, 14, figsize=(20, 6))
i = 0
for cfg, data0 in C.iteritems():
for lbl, data1 in data0.iteritems():
Qe = np.array(data1).T[0]
m = list(Qe).index(Qe.min())
n = range(len(Qe))
Qe = Qe - Qe[m]
Ne = np.array(data1).T[1] - np.array(data1).T[1][m]
Re = np.array(data1).T[2] - np.array(data1).T[2][m]
Nerr = Ne - Qe
Rerr = Re - Qe
# Plotting energy points
ax[0, i].plot(n, Ne - Ne[m], 'rs', label='BPNN')
ax[0, i].plot(n, Re - Re[m], 'b^', label='ReaxFF')
ax[0, i].plot(n, Qe - Qe[m], 'ko', fillstyle='none', label='DFT')
# Plotting the residuals
ax[1, i].plot([min(n), max(n)], [0, 0], 'k:')
ax[1, i].plot(n, Nerr, 'r.')
ax[1, i].plot(n, Rerr, 'b.')
ax[0, i].set_xlim(min(n), max(n))
ax[1, i].set_xlim(min(n), max(n))
ax[1, i].set_ylim(-0.4, 0.4)
ax[0, i].set_ylim(-0.25, 1.6)
ax[0, i].set_title('{0}\n{1}'.format(cfg, lbl))
ax[0, i].set_xticklabels([])
ax[1, i].set_xticklabels([])
if i != 0:
ax[0, i].set_yticklabels([])
ax[1, i].set_yticklabels([])
i += 1
ax[1, 0].set_yticks(ax[1, 0].get_yticks()[:-2])
ax[1, 0].set_ylabel('Residual error (eV)')
ax[0, 0].set_ylabel('Total energy (eV)')
ax[0, 0].legend(loc='best', numpoints=1, fontsize=12, frameon=False)
plt.tight_layout(w_pad=0.1, h_pad=-0.3)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-full-diffusion.{0}'.format(ext), dpi=300)from ase.db import connect
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import fmin
from ase.io import write
import matplotlib.image as mpimg
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import os
db = connect('data.db')
Qe, Re, Ne = [], [], []
for d in db.select(['miller=111', 'diffusion=slipping', 'config!=double', 'traj=0']):
Qe += [d.energy]
Ne += [d.neural_energy]
Re += [d.reax_energy]
Qe = np.array([Qe[-1]] + Qe) - Qe[-1]
Ne = np.array([Ne[-1]] + Ne) - Ne[-1]
Re = np.array([Re[-1]] + Re) - Re[-1]
x = np.linspace(0, len(Qe) - 1)
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
for nrg, name, c, o in [[Re, 'ReaxFF', 'b^', 0.0],
[Ne, 'BPNN', 'rs', 0.0],
[Qe, 'DFT', 'ko', 0.0]]:
f = interp1d(range(len(nrg)), -nrg, 'cubic')
xmax = fmin(f, len(nrg) / 2., disp=False)
ax.plot(x, -f(x), c[0] + '--')
ax.plot(nrg, c, label=name)
ax.annotate('{0} $E^\ddag$= {1:1.2f}'.format(name, float(-f(xmax))),
xy=(xmax, -f(xmax)), xytext=(4.0, -f(xmax) + o),
ha='right', va='center', color=c[0],
arrowprops=dict(arrowstyle="->",
shrinkB=10, color=c[0]))
for i, a in enumerate([[0, 1.9], [2, 5.1], [7, 8.3]]):
atoms = db.get_atoms(['miller=111',
'diffusion=slipping',
'config!=double',
'traj=0', 'image={0}'.format(a[0])])
del atoms[0]
atoms *= (3, 3, 1)
write('./images/temp.png',
atoms,
colors=['#333333', '#999999', '#CCCCCC',] * 27,
scale=20,
show_unit_cell=2,
radii=0.75)
image = mpimg.imread('./images/temp.png')
imagebox = OffsetImage(image, zoom=1.5)
ax.add_artist(AnnotationBbox(imagebox,
xy=(a[0], nrg[a[0]]),
xybox=(a[1], 0.25),
pad=-0.2,
frameon=False,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=-90,armA=0')
)
)
os.unlink('./images/temp.png')
ax.set_xticklabels([])
plt.xlabel('Reaction coordinate (a.u.)')
plt.ylabel('Potential energy (eV)')
plt.ylim(0, 0.3)
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-111-slipping.{0}'.format(ext), dpi=300)To conserve space in the manuscript, only the fcc(111) single-layer surface slipping barrier is shown, while the fcc(100) single- and double-layer barriers are depicted in this section.
from ase.db import connect
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import fmin
from ase.io import write
import matplotlib.image as mpimg
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import os
db = connect('data.db')
Qe, Re, Ne = [], [], []
for d in db.select(['miller=100', 'diffusion=slipping', 'config!=double', 'traj=0']):
Qe += [d.energy]
Ne += [d.neural_energy]
Re += [d.reax_energy]
Qe = np.array([Qe[-1]] + Qe) - Qe[-1]
Ne = np.array([Ne[-1]] + Ne) - Ne[-1]
Re = np.array([Re[-1]] + Re) - Re[-1]
x = np.linspace(0, len(Qe) - 1)
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
for nrg, name, c, o in [[Qe, 'DFT', 'ko', 0.0275],
[Ne, 'BPNN', 'rs', 0.04],
[Re, 'ReaxFF', 'b^', -0.010]]:
f = interp1d(range(len(nrg)), -nrg, 'cubic')
xmax = fmin(f, len(nrg) / 2., disp=False)
ax.plot(x, -f(x), c[0] + '--')
ax.plot(nrg, c, label=name)
ax.annotate('{0} $E^\ddag$= {1:1.3f}'.format(name, float(-f(xmax))),
xy=(xmax, -f(xmax)), xytext=(6.7, -f(xmax) + o),
ha='left', va='center', color=c[0],
arrowprops=dict(arrowstyle="->",
shrinkB=10, color=c[0]))
for i, a in enumerate([[0, 1.9], [5, 5.1]]):
atoms = db.get_atoms(['miller=100',
'diffusion=slipping',
'config!=double',
'traj=0', 'image={0}'.format(a[0])])
del atoms[0]
atoms *= (3, 3, 1)
write('./images/temp.png',
atoms,
colors=['#333333', '#999999', '#CCCCCC',] * 27,
scale=20,
show_unit_cell=2,
radii=0.75)
image = mpimg.imread('./images/temp.png')
imagebox = OffsetImage(image, zoom=1.5)
ax.add_artist(AnnotationBbox(imagebox,
xy=(a[0], Ne[a[0]]),
xybox=(a[1], 0.53),
pad=-0.2,
frameon=False,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=-90,armA=0')
)
)
os.unlink('./images/temp.png')
ax.set_xticklabels([])
plt.xlabel('Reaction coordinate (a.u.)')
plt.ylabel('Potential energy (eV)')
plt.ylim(0, 0.64)
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-100-slipping.{0}'.format(ext), dpi=300)
from ase.db import connect
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
from scipy.optimize import fmin
from ase.io import write
import matplotlib.image as mpimg
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import os
db = connect('data.db')
Qe, Re, Ne = [], [], []
for d in db.select(['miller=100', 'diffusion=slipping', 'config!=single', 'traj=0']):
Qe += [d.energy]
Ne += [d.neural_energy]
Re += [d.reax_energy]
Qe = np.array([Qe[-1]] + Qe) - Qe[-1]
Ne = np.array([Ne[-1]] + Ne) - Ne[-1]
Re = np.array([Re[-1]] + Re) - Re[-1]
x = np.linspace(0, len(Qe) - 1)
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
for nrg, name, c, o in [[Qe, 'DFT', 'ko', 0.0375],
[Ne, 'BPNN', 'rs', 0.06],
[Re, 'ReaxFF', 'b^', 0.06]]:
f = interp1d(range(len(nrg)), -nrg, 'cubic')
xmax = fmin(f, len(nrg) / 2., disp=False)
ax.plot(x, -f(x), c[0] + '--')
ax.plot(nrg, c, label=name)
ax.annotate('{0} $E^\ddag$= {1:1.3f}'.format(name, float(-f(xmax))),
xy=(xmax, -f(xmax)), xytext=(6.7, -f(xmax) + o),
ha='left', va='center', color=c[0],
arrowprops=dict(arrowstyle="->",
shrinkB=10, color=c[0]))
for i, a in enumerate([[0, 1.9], [5, 5.1]]):
atoms = db.get_atoms(['miller=100',
'diffusion=slipping',
'config!=single',
'traj=0', 'image={0}'.format(a[0])])
del atoms[0]
atoms *= (3, 3, 1)
write('./images/temp.png',
atoms,
colors=['#333333', '#999999', '#CCCCCC',] * 27,
scale=20,
show_unit_cell=2,
radii=0.75)
image = mpimg.imread('./images/temp.png')
imagebox = OffsetImage(image, zoom=1.5)
ax.add_artist(AnnotationBbox(imagebox,
xy=(a[0], Ne[a[0]]),
xybox=(a[1], 0.53),
pad=-0.2,
frameon=False,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=-90,armA=0')
)
)
os.unlink('./images/temp.png')
ax.set_xticklabels([])
plt.xlabel('Reaction coordinate (a.u.)')
plt.ylabel('Potential energy (eV)')
plt.ylim(0, 0.64)
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-100-slipping-2.{0}'.format(ext), dpi=300)from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
from ase.io import write
import matplotlib.pyplot as plt
import numpy as np
from ase.db import connect
import os
db = connect('md-6atom.db')
Qe = dict.fromkeys(range(99, 2000, 100), 0)
I, Re, Ne = [], [], []
for d in db.select():
I += [d.image]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
# Only every 100th structure is DFT validated
if d.image in Qe.keys():
Qe[d.image] = d.energy / d.natoms
sort = [i[0] for i in sorted(enumerate(I), key=lambda x:x[1])]
I = np.array([I[i] for i in sort])
Ne = np.array([Ne[i] for i in sort])
Re = np.array([Re[i] for i in sort])
m = list(Ne).index(min(Ne))
# Create an image of the DFT predicted minimum
db0 = connect('data.db')
gm_atoms = db0.get_atoms(['cluster=plane-fcc', 'traj=0', 'config=2', 'natoms=6'])
gm = db0.get(['cluster=plane-fcc', 'traj=0', 'config=2', 'natoms=6'])
write('./images/tmp-gm.png', gm_atoms)
# Create an image of the BPNN predicted minimum
write('./images/tmp-pm.png', db.get_atoms(m), rotation='45y')
# Establish baseline for global minimum
Q_min = gm.energy / gm.natoms
N_min = gm.neural_energy / gm.natoms
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
fig0 = OffsetImage(mpimg.imread('./images/tmp-pm.png'))
ax.add_artist(AnnotationBbox(fig0,
xy=(m, Ne[m]),
xybox=(1600, -1.8),
pad=0,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=0')
))
fig = OffsetImage(mpimg.imread('./images/tmp-gm.png'))
ax.add_artist(AnnotationBbox(fig,
xy=(m, Q_min),
xybox=(1400, -2.22),
pad=0,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=0')
))
# Remove the temporary images
os.unlink('./images/tmp-gm.png')
os.unlink('./images/tmp-pm.png')
ax.text(1700, -2.16, 'DFT', color='k')
ax.text(200, -2.27, 'ReaxFF', color='b')
ax.text(200, -1.9, 'BPNN', color='r')
ax.plot(range(len(Ne)), Re, 'b-')
ax.plot([0, len(Ne)], [N_min, N_min], 'r--')
ax.plot([0, len(Ne)], [Q_min, Q_min], 'k--')
for k, v in Qe.iteritems():
plt.plot([k, k], [v, Ne[k]], 'k:')
ax.plot(range(len(Ne)), Ne, 'r-')
ax.scatter(Qe.keys(), Qe.values(), color='k', marker='o', zorder=10)
ax.set_xlim(0, 2000)
ax.set_ylim(-2.3, -1.7)
plt.xlabel('Time step')
plt.ylabel('Potential energy (eV/atom)')
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-6atom-md.{0}'.format(ext), dpi=300)To demonstrate how the BPNN “learns” the potential energy surface as the number of training points increases, we have included a 6 atom MD simulation from one of the earlier BPNNs created. This BPNN was trained to all of the cluster calculations included in the full database up to 13 atoms large. However, no 6 atom clusters were used in the training set in order to observe how well these structures could be extrapolated. This is a training set of approximately 2000 calculations. The MD trajectory can be found here: \attachfile{md-6atom-old.db}{(double-click to open)}.
The resulting 6 atom MD simulation was performed identically to the one described in the paper. The results of the 4000 step MD are shown in Figure ref:si-6atom-md. For easier comparison, the energies are plotted on the same scale as the 6 atom MD simulation reported in the paper which used the full database as a training set. The result is a substantial improvement in the performance of the full BPNN. This demonstrates how a BPNN can be made to obtain arbitrary levels of accuracy even with a large variety of different structure types being used for training.
Although the errors of the early 6 atom MD simulation shown here are significantly larger, the BPNN still accurately predicts the planer structure to be the lowest energy configuration. Perhaps even more impressive is that it manages to do so without any information about 6 atom structures.
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
from ase.io import write
import matplotlib.pyplot as plt
import numpy as np
from ase.db import connect
import os
db = connect('md-6atom-old.db')
Qe = {}
I, Ne, = [], []
for d in db.select():
I += [d.image]
Ne += [d.neural_energy/ d.natoms]
if d.DFT:
Qe[d.image] = d.energy / d.natoms
sort = [i[0] for i in sorted(enumerate(I), key=lambda x:x[1])]
I = np.array([I[i] for i in sort])
Ne = np.array([Ne[i] for i in sort])
m = list(Ne).index(min(Ne))
# Create an image of the DFT predicted minimum
db0 = connect('data.db')
gm_atoms = db0.get_atoms(['cluster=plane-fcc', 'traj=0', 'config=2', 'natoms=6'])
gm = db0.get(['cluster=plane-fcc', 'traj=0', 'config=2', 'natoms=6'])
write('./images/tmp-gm.png', gm_atoms)
# Create an image of the BPNN predicted minimum
write('./images/tmp-pm.png', db.get_atoms(m), rotation='45y')
# Establish baseline for global minimum
Q_min = gm.energy / gm.natoms
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
fig0 = OffsetImage(mpimg.imread('./images/tmp-pm.png'))
ax.add_artist(AnnotationBbox(fig0,
xy=(m, Ne[m]),
xybox=(1200, -2.0),
pad=0,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=0')
))
fig = OffsetImage(mpimg.imread('./images/tmp-gm.png'))
ax.add_artist(AnnotationBbox(fig,
xy=(m, Q_min),
xybox=(3400, -2.22),
pad=0,
arrowprops=dict(arrowstyle='->',
color='0.5',
zorder=5,
connectionstyle='arc,angleA=0')
))
# Remove the temporary images
os.unlink('./images/tmp-gm.png')
os.unlink('./images/tmp-pm.png')
ax.text(2500, -2.16, 'DFT', color='k')
ax.text(200, -1.9, 'BPNN', color='r')
ax.plot([0, len(Ne)], [Q_min, Q_min], 'k--')
for k, v in Qe.iteritems():
plt.plot([k, k], [v, Ne[k]], 'k:')
ax.plot(range(len(Ne)), Ne, 'r-')
ax.scatter(Qe.keys(), Qe.values(), color='k', marker='o', zorder=10)
ax.set_xlim(0, 4000)
ax.set_ylim(-2.3, -1.7)
plt.xlabel('Time step')
plt.ylabel('Potential energy (eV/atom)')
plt.tight_layout()
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/si-6atom-md.{0}'.format(ext), dpi=300)import matplotlib.pyplot as plt
import numpy as np
from ase.db import connect
db = connect('data.db')
Qe, Re, Ne = [], [], []
for d in db.select(['post=minima']):
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
Re = np.array(Re)
Nre = np.array(Ne) - np.array(Qe)
Rre = Re - np.array(Qe)
C = np.array(range(len(Qe))) + 1
f, ax = plt.subplots(2, sharex=True, figsize=(6, 4))
ax[0].plot(C, Ne, mec='none', mfc='r', marker='s', lw=0, label='BPNN')
ax[0].plot(C, Re + 0.11, mec='none', mfc='b', marker='^', lw=0, label='ReaxFF')
ax[0].plot(C, Qe, mec='k', mfc='none', marker='o', lw=0, label='DFT')
ax[0].text(12, -2.6, 'ReaxFF offset by +0.11 eV/atom',
va='bottom', color='b', fontsize=12)
ax[1].plot([C[0], C[-1]], [0, 0], 'k--')
ax[1].scatter(C, Nre, c='r', marker='s', edgecolor='none')
ax[1].scatter(C, Rre, c='b', marker='^', edgecolor='none')
ax[0].set_yticks([-2.65, -2.63, -2.61, -2.59])
ax[1].set_xlabel('MD minima')
ax[0].set_ylabel('Potential energy\n(eV/atom)')
ax[1].set_ylabel('Residual error\n(eV/atom)')
ax[0].set_xlim(C[0], C[-1])
ax[1].set_xlim(C[0], C[-1])
ax[0].legend(loc='best', numpoints=1, fontsize=12, frameon=False)
plt.tight_layout(h_pad=-0.0)
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/fig-38atom-minima.{0}'.format(ext), dpi=300)This is the Table of Contents graphic.
from ase.db import connect
import numpy as np
import matplotlib.pyplot as plt
from ase.io import write
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
import os
db = connect('data.db')
Qe, Re, Ne = [], [], []
for d in db.select(['post=minima']):
Qe += [d.energy / d.natoms]
Ne += [d.neural_energy / d.natoms]
Re += [d.reax_energy / d.natoms]
DFT_time = db.get(['post=minima', 'config=121']).calc_time
pos = [(0.5, 0.8), (0.2, 0.2), (0.8, 0.2)]
pos2 = [(0.5, 0.6), (0.35, 0.4), (0.65, 0.4)]
rot = ["-15x, -130z", "", ""]
col = ['k', 'r', 'b']
pad = [-0.1, -0.7, -0.7]
fig = plt.figure(figsize=(5.5, 5))
ax = fig.add_subplot(111)
for i, nrgs in enumerate([Qe, Re, Ne]):
ind = nrgs.index(min(nrgs))
matoms = db.get_atoms(['post=minima',
'config={0}'.format(ind)])
write('./images/temp.png'.format(ind), matoms, rotation=rot[i])
image = mpimg.imread('./images/temp.png')
imagebox = OffsetImage(image, zoom=1.75)
ax.add_artist(AnnotationBbox(imagebox,
xy=pos2[i],
xybox=pos[i],
pad=pad[i],
frameon=False,
arrowprops=dict(arrowstyle='<-',
color=col[i],
zorder=5)
)
)
os.unlink('./images/temp.png')
plt.text(pos2[0][0], pos2[0][1], 'DFT',
ha='center', size=20, va='top')
ax.annotate('Neural\nNetwork', xy=(0.45, 0.58), xytext=pos2[1],
va='bottom', ha='center', color='r', size=20,
arrowprops=dict(arrowstyle="<|-", fc='r',
connectionstyle="arc3,rad=-0.5",))
ax.annotate('ReaxFF', xy=(0.55, 0.58), xytext=pos2[2],
va='bottom', ha='center', color='b', size=20,
arrowprops=dict(arrowstyle="<|-", fc='b',
connectionstyle="arc3,rad=0.6",))
plt.text(0.5, 0.98, '{0:1.1f} hrs'.format(DFT_time/3600.),
va='center', ha='center', size=17)
plt.text(0.2, 0.01, '0.14 s', color='r',
va='center', ha='center', size=17)
plt.text(0.8, 0.01, '0.01 s', color='b',
va='center', ha='center', size=17)
fig.patch.set_visible(False)
ax.axis('off')
plt.tight_layout(rect=[-0.15, -0.03, 1.1, 1.02])
for ext in ['png', 'eps', 'pdf']:
plt.savefig('./images/toc.{0}'.format(ext), dpi=300)bibliographystyle:unsrt bibliography:./manuscript.bib