Geometry Optimization

In order to run growing string calculations, first you need the optimized geometry of the reactant and product. For surface reactions, we need to define the slab and adsorbate species. The atomic simulation environment (ASE) is the tool we are going to use. ASE is implemented in Python and automates the process to some extent (e.g, if you are using VASP as the calculator, it creates all the input files). The following script is the general form of a geometry optimization script.

There are comments added to the code explaining each part. If VASP is used as the calculator (which happens almost always), make sure to run each job in a separate directory. To run an optimization job, copy this script and a queue submission script (example given in general information section) in a directory. There are also a couple of settings needed in your bashrc file (check the general information section).

For simple diatomic molecules like CO, the position and bond length of adsorbate (CO) can be defined using the script (like the one below), but there are cases where you have a more complicated adsorbate (like propanoic acid) and you want to adjust the position of the molecule on surface. In this case, use ASE to create the slab and optimize it, then using Avogadro (check the general information section) add molecule to the surface. After this step, surface and adsorbate should be optimized together, so copy (see general information section on how to copy from/to the cluster) the structure created by Avogadro as input.xyz file to the same directory and add the molecule to the surface in the script (lines 18 & 23). The order of atoms in input.xyz file and on the script (line 18) should be conserved. To set up the surface with new coordinates read from input.xyz, the lines 30-32 should be in the script.

More examples of setting up metal and metal oxide surfaces are given in surface set up section.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/python
#Set up a simple hcp0001 slab of Ru

from ase import Atoms
from ase.calculators.vasp import Vasp
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.lattice.surface import hcp0001
from ase.io import write,read

#surface type, name of the atom, size of the slab which
#has 3 atoms in length and width and 2 layers of atoms
slab = hcp0001('Ru', size=(3,3,4), vacuum=10)

#The molecule adsorbed on surface. The first set of
#coordinates correspond to C atom. 1.16 is the bond
#length of CO molecule.
molecule = Atoms('CO', [(0., 0., 0.), (0., 0., 1.16)])

#add molecule on the slab. The distance netween C atom
#and surface is 1.7 A. ontop is the position the molecule
#will be added.
add_adsorbate(slab, molecule, 1.7, 'ontop', offset=1)

#If you have a structure prepared using a visulization
#software like Avogadro, it can be given as an input
#file to re-define the structure. Make sure the input
#structure and slab and molecule created in the first
#code block match (number and order of atoms)
fName = 'input.xyz'
newSlab = read(fName)
slab.set_positions(newSlab.get_positions())

#set the calculator. For information on each tag, see the VASP manual.
calc = Vasp(xc='PBE',lreal='Auto',kpts=[2,2,1],ismear=1,sigma=0.2,algo='fast',istart=0, npar=8)
slab.set_calculator(calc)

#set constraints like frozen layers. There are two ways to do so; 1. freezing the whole layer or 2. freezing specific atoms
#method 1. 2 refers to the top 2 ayers. this means layers 3 and below are frozen
mask = [atom.tag > 2 for atom in slab]
slab.set_constraint(FixAtoms(mask=mask))
#method 2. Only the atoms specified here are frozen. note that it starts from 0
c = FixAtoms(indices=[0, 1, 2, 3, 4, 12, 13, 14, 15, 16])
slab.set_constraint(c)

#running optimization and creating a trajectory file
dyn = QuasiNewton(slab, trajectory='CO-Ru.traj')
dyn.run(fmax=0.05) #optimization condition (eV/A)

#writing the output structure in xyz format
write('outPut.xyz', slab)