ACSFs are a convenient way of transforming atomic coortinates and types into a computer-friendly string of numbers. Each atom gets its own set of ACSFs, computed using itself as the center, and all other atomic coordinates, which encode its chemical environment.
The ACSFs are invariant under translation and rotations of the whole system, and independent on the ordering of the atoms.
The two main type of ACSFs are two- and three-body. Each set of ACSFs becomes the input of a neural network that calculates the corresponding energy contribution. The only important quantity is the total energy of the system, given by the sum of all contributions. All neural networks are fitted to minimise the error on total energy.
Derivatives of the energy w.r.t. atomic positions (forces) can be obtained analytically with the chain rule, and this makes ACSF the method of choice for molecular dynamics simulations.
For more info see: Jörg Behler, J. Chem. Phys. 134, 074106 (2011)
In [ ]:
# --- INITIAL DEFINITIONS ---
from sklearn.neural_network import MLPRegressor
import numpy, math, random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from visualise import view
from ase import Atoms
In [ ]:
# parameters for the trimer
d1 = 1.0
d2 = 1.0
theta = 30
# convert angle to radians
theta *= math.pi/180
# generate the atomic coordinates
atoms = numpy.zeros((3,3))
atoms[1] = [math.cos(theta),-math.sin(theta),0]
atoms[2] = [-math.cos(theta),-math.sin(theta),0]
atoms = numpy.array(atoms)
atoms_sys = Atoms(positions=atoms, symbols=['H']*3)
# show the atomic position vectors - on rows
print(atoms)
# show atomic positions - drag mouse to rotate view
view(atoms_sys)
In [ ]:
# cutoff distance for all ACSFs
Rcut = 5.0
def fcut(r):
return (math.cos(math.pi * r/Rcut)+1) * 0.5
# compute the ACSF G0 (sum of cutoff)
G0 = numpy.zeros(3)
for i in range(3): # compute one ACSF for each atom i
for j in range(3): # loop over all OTHER atoms
if i == j: continue
r = numpy.linalg.norm(atoms[i]-atoms[j]) # compute distance i-j
if r > Rcut: continue # skip if too far
G0[i] += fcut(r)
print(G0)
The result is one number for each atom in the system. Larger values indicate more atoms in the vicinity. The last two atoms are a bit far from each other, thus their G0 is smaller than the central atom.
The next type of ACSF gives more information about distances to neighbouring atoms:
$$G^1_i = \sum_j e^{-\eta\left(r_{ij} - R_s \right)^2} f_c\left(r_{ij}\right)$$where $\eta$ and $R_s$ are free parameters. The ACSF is a gaussian, and we can imagine that $G^1$ quantifies how many other atoms are found $R_s$ away from the central atom.
In practice, $G^1$ is computed for several combinations of $\eta$ and $R_s$, thus each atom will have a series of $G^1$ values ($G^1$ is an array for each atom)
In [ ]:
# G1 function definition
def G1f(x, eta, Rs):
return math.exp(-eta*(x-Rs)*(x-Rs)) * ((math.cos(math.pi * x/Rcut)+1) * 0.5)
G1v = numpy.vectorize(G1f)
# this is a list of eta/Rs pairs to use in our example
p = [(0.4, 0.2),(0.4, 0.5),(0.4, 1.0),(0.5, 2.0),(0.5, 3.0),(0.5, 4.0)]
# make a nice plot of the ACSFs we are using
datar = numpy.arange(0, Rcut, Rcut/100)
for pp in p:
plt.plot(datar, G1v(datar, pp[0], pp[1]) )
plt.xlabel('Pair distance')
plt.ylabel('G1')
plt.show()
# compute the ACSF G1
G1 = numpy.zeros((3, len(p)))
for i in range(3): # compute one ACSF for each atom i
for j in range(3): # loop over all OTHER atoms
if i == j: continue
r = numpy.linalg.norm(atoms[i]-atoms[j]) # compute distance i-j
if r > Rcut: continue # skip if too far
for k in range(len(p)):
eta, rs = p[k]
G1[i,k] += G1f(r, eta, rs)
print(G1)
The result is a matrix: one row for each atom, one column for each pair or parameters $\eta$ and $R_s$.
A more complete description of the system can be achieved by adding ACSFs that depend on triplets of atoms, giving three-body terms to the energy expression. Usually they are functions of the cosine of the angle between three atoms, and the distances between them:
$$G^4_i = 2^{1-\zeta}\sum_{j,k} \left(1+\lambda\cos\left( \theta_{ijk} \right)\right)^\zeta \exp\left[-\eta\left(r_{ij}^2 +r_{ik}^2+r_{jk}^2\right)\right] f_c\left(r_{ij}\right)f_c\left(r_{ik}\right)f_c\left(r_{jk}\right)$$where $\eta$, $\zeta$ and $\lambda$ are free parameter. Usually $\zeta$ is an integer and $\lambda = \pm 1$.
In [ ]:
# G4 function definition
def G4f(cost, rij, rik, rjk, eta, zeta, lamb):
result = math.pow(2,1-zeta) * math.pow(1+lamb*cost, zeta)
result *= math.exp(-eta*(rij*rij + rik*rik + rjk*rjk))
return result * fcut(rij) * fcut(rjk) * fcut(rik)
# this is a list of eta/zeta/lambda pairs to use in our example
p = [(0.4, 1, 1),(0.4, 2, 1),(0.4, 3, 1),(0.8, 1, 1),(0.8, 2, 1)]
# compute the ACSF G4
G4 = numpy.zeros((3, len(p)))
for i in range(3): # compute one ACSF for each atom i
for j in range(3): # loop over all OTHER atoms
if i == j: continue
rij = numpy.linalg.norm(atoms[i]-atoms[j]) # compute distance i-j
if rij > Rcut: continue # skip if too far
for k in range(3):
if k == i or k == j: continue
rik = numpy.linalg.norm(atoms[i]-atoms[k]) # compute distance i-k
rjk = numpy.linalg.norm(atoms[j]-atoms[k]) # compute distance j-k
# compute cosine of angle
cost = numpy.dot(atoms[i]-atoms[k],atoms[i]-atoms[j])
cost /= rij * rik
for k in range(len(p)):
eta, zeta, lamb = p[k]
G4[i,k] += G4f(cost,rij,rik,rjk, eta, zeta, lamb)
print(G4)
These are some of the functions suggested by the authors, but more can be crafted. Two-body ACSFs must only depend on the pair distance, while three-body ones must include angular dependence.
In practice ACSFs are subdivided according to the atomic types. For example, if the system only has C and H, each atom $i$ of type C will have two sets of $G^1$:
$$G^1_{i-C} = \sum_{j\in C} e^{-\eta\left(r_{ij} - R_s \right)^2} f_c\left(r_{ij}\right)$$$$G^1_{i-H} = \sum_{j\in H} e^{-\eta\left(r_{ij} - R_s \right)^2} f_c\left(r_{ij}\right)$$and similarly for three-body ACSFs. The chemical environments w.r.t. different atomic species are separate, and go to different neural networks. This way, more information is retained in the ACSF descriptors but more NN parameters will be necessary to accurately model the system, and more training data has to be supplied.
In [ ]:
In [ ]:
In [ ]: