View source for 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 ([[Documentation#turbogap predict|turbogap predict]]) and [[geometry optimization]].

== Prerequisites for this tutorial ==

* A '''TurboGAP''' installation
* An [https://wiki.fysik.dtu.dk/ase/ 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:

<syntaxhighlight lang="python" line="line">
from ase.io import write
from ase import Atoms
import numpy as np

append = False
for a in np.arange(3.8,4.3,0.01):
    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)
    write("fcc.xyz", atoms, append=append)
    append = True
</syntaxhighlight>

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

<syntaxhighlight lang="bash" line="line">
wget https://zenodo.org/record/6302852/files/gap_files.tar.gz
tar -xvf gap_files.tar.gz
</syntaxhighlight>

The following step is to set up an <code>input</code> 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 <code>trajectory_out.xyz</code> we retrive the lattice parameters and energies with the following script:

<syntaxhighlight lang="python" line="line">
from ase.io import read
db = read("trajectory_out.xyz", index=":")
f = open("e.dat", "w+")
for atoms in db:
    e = atoms.info["energy"]
    a = atoms.get_cell()[0][1]*2
    print(a, e, file=f)
</syntaxhighlight>

Run it an save the results to a file (I save the Python script to a file called <code>extract_energy.py</code> 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; E_opt = a+b*opt+c*opt**2
 print opt, E_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"

I get <math>a_\text{opt} = 4.1536 \text{Å}</math> and <math>E_\text{opt} = -3.0448 \text{eV}</math>, and the plot looks like this:

[[File:tutorial4_01.png|center|480px|]]

== Setting up the gold cluster ==

Now that we know what is the reference fcc lattice parameter and its energy, let's construct an icosahedral gold cluster and relax it. We use ASE for constructing the cluster. ASE has a lot of useful functions for constructing clusters, surfaces and other common atomic structures:

<syntaxhighlight lang="python" line="line">
from ase.cluster import Icosahedron
from ase.io import write

atoms = Icosahedron("Au", 8, 4.1536)
atoms.center(vacuum=10)

write("ico.xyz", atoms)
</syntaxhighlight>

In this case we built an icosahedron with 8 atomic layers and surrounded it by 10 Angstrom of vacuum in all directions.

After running the Python script, we will have an <code>ico.xyz</code> file with the unrelaxed gold cluster in it:

[[File:tutorial4_02.png|center|480px|]]

We will now proceed to relax this cluster using <code>turbogap md</code>, but in this case not to perform molecular dynamics but to carry out a [[gradient descent]] geometry optimization. This is the needed <code>input</code> file:

 atoms_file = "ico.xyz"
 pot_file = "gap_files/aurum.gap"
 
 species = Au
 n_species = 1
 
 [[optimize]] = "gd"
 write_xyz = 1
 md_nsteps = 1000

Run '''TurboGAP''' to relax the structure (this will take 1-2 minutes):

 mpirun -np 4 turbogap md

In the <code>trajectory_out.xyz</code> file we can find the energy for this structure. Check the energy of the relaxed cluster, i.e., the last frame in the trajectory file. I get <math>E = -4085.948634 \,\text{eV}</math> for <math>N = 1415</math> atoms. Therefore, the cohesive energy per atom is <math>e = -4085.948634/1415 \,\text{eV} + 3.0448 \,\text{eV} = 0.157 \,\text{eV}</math> versus fcc.