Simulating icosahedral gold clusters

In this tutorial we are going to use a GAP to optimize the structure of an icosahedral gold cluster and estimate its formation energy. In this tutorial we will make use of static (single-point) calculations (turbogap predict) and geometry optimization.

Prerequisites for this tutorial

• A TurboGAP installation
• An ASE installation
• Numpy
• A utility that can perform simple least-squares regression (I will use gnuplot, which lets me plot the fit at the same time)

Computing the equilibrium lattice parameter of the reference structure

The thermodynamically stable crystal structure of gold is fcc. When we compute the formation energy of a gold nanocluster we will be referring its energy per atom to that of the fcc structure. Therefore, the first step in this tutorial is to find out the lattice parameter of fcc gold by computing the energy vs lattice parameter curve, and finding the minimum.

This script will generate a series of Au fcc structures with lattice parameter between 3.8 and 4.3 Angstrom:

1 from ase.io import write
2 from ase import Atoms
3 import numpy as np
4 
5 append = False
6 for a in np.arange(3.8,4.3,0.01):
7     atoms = Atoms("Au", positions=[[0,0,0]], cell=[[0,a/2,a/2],[a/2,0,a/2],[a/2,a/2,0]], pbc=True)
8     write("fcc.xyz", atoms, append=append)
9     append = True

We now have the fcc.xyz file containing all of these structures. Next, we retrieve the Au GAP potential from Zenodo:

1 wget https://zenodo.org/record/6302852/files/gap_files.tar.gz
2 tar -xvf gap_files.tar.gz

The following step is to set up an input file:

atoms_file = "fcc.xyz"
pot_file = "gap_files/aurum.gap"

species = Au
n_species = 1

and run Turbogap in single-point mode:

mpirun -np turbogap predict

From the resulting trajectory_out.xyz we retrive the lattice parameters and energies with the following script:

1 from ase.io import read
2 db = read("trajectory_out.xyz", index=":")
3 f = open("e.dat", "w+")
4 for atoms in db:
5     e = atoms.info["energy"]
6     a = atoms.get_cell()[0][1]*2
7     print(a, e, file=f)

Run it an save the results to a file (I save the Python script to a file called extract_energy.py here:

python3 extract_energy.py > e.dat

We can open this file in gnuplot and use gnuplot's fitting capabilities to fit a second order polynomial to the data (I restrict the fitting domain to [4.1:4.2] to get a better fit):

fit[4.1:4.2] a+b*x+c*x**2 "e.dat" via a,b,c
opt=-b/2/c
print opt # this will tell you the lattice parameter at the optimum
set xlabel "Lattice parameter (A)"; set ylabel "Energy (eV)"
plot "e.dat" w l t "Data", a+b*x+c*x**2 t "Fit"