Slab surface tutorial for SeqQuest: CO/Ru(0001)

Introduction

The purpose of this tutorial is to introduce the user to the effective use of the code for two-dimensional slab calculations for the study of surface properties of materials. This tutorial uses a series of examples used in setting up a calculation culminating in the adsorption CO on the Ru(0001) surface.

It will be helpful (but not necessary) in the following that the user already be somewhat familiar with use of the code for 3D bulk calculations (see the Bulk Tutorial). The appropriate lattice parameters to use in the surface plane are dictated by the corresponding bulk lattice parameters theoretically computed in the bulk material using the same level of approximation, so that 3D bulk calculations are a natural prerequisite to any 2D slab calculation to study chemistry at the surface of a material.

Two-dimensional slab calculations have special issues to consider in comparison to calculations for bulk materials, because the slab is periodic in two dimensions and is finite in one. The code must treat the boundary conditions differently in the periodic vs. the non-periodic directions. In addition, the code runs a formally a supercell calculation, treating a slab calculation as an infinite series of slabs separated by vacuum, and needs to be told where the “vacuum boundary” is in order to set up the boundary conditions correctly. These needs impose special requirements on a slab calculation that are reflected in the input conventions described below. First, note that the periodic directions of the slab unit cell must be in the x-y plane (the first two primitive lattice vectors), and that the third lattice vector (which must be specified) must be normal to the slab (i.e. along the z-axis). Second, while the code is formally a supercell code, and therefore subject to supercell effects such as the interaction between the planar dipoles of different slab images, Quest automatically and exactly eliminates the dipole interaction using the Local Moment CounterCharge (LMCC) method:
“Local electrostatic moments and periodic boundary conditions”,
P.A. Schultz, Phys. Rev. B 60, 1551 (1999)
resulting in slabs that are Coulombically isolated from one another. With a localized basis set, that does not cross supercell boundaries, this means that the calculation does indeed represent an isolated slab.

Units

The default units for energy are Rydberg (1 Rydberg = 0.5 Hartree = 13.605 eV), and for distances are bohr (1 bohr = 0.52918 Ang = 0.052918 nm). All input coordinates will be in these units, though one can use the scaling functions in the code to aid in the conversions.

Input files

The code uses one main input file that the user must construct. Atom (potential+basis) files, which will need to be listed in the input file, can be obtained from the atom libraries. The input file is a text-rich keyword-driven document that, once constructed, is self-documenting. The input is a sequence of keyword lines with the indicated data on following lines. Only the first few (6 or 8) characters of the keyword line are significant. The remainder of a keyword line can be used for additional documentation. The data itself, in general, is free-format. The input in the setup data section is highly structured, and requires its input in a very specific order. Other sections take their data in any order. If you miss putting something in the right place in the input, not to worry. The code will very happily process the input file, and stop when it cannot find required input in the right place. It echoes the input file as it reads it, and will politely tell you what it was looking for when it runs into a problem. Interspersed among the required inputs, are a variety of usually invisible optional input sections. The code uses a variety of defaults that can be overridden in the input file. The use of many optional inputs will be illustrated in the tutorial.

Example: A basic input file

Illustrated in this first example is a very basic file for doing a 7-layer Ru(0001) slab calculation (i.e. the hexagonal surface of obtained by cleaving along the basal plane of bulk hcp Ru). In preparation to setting up this calculation, we had done a 3D bulk calculation for hcp Ru to determine the bulk lattice parameters (using the cell optimization capability of the code) using the same flavor of density functional theory we will be using for the slab calculations, the Perdew/Zunger parameterization of the Ceperley/Alder electron gas results (CAPZ) within LDA. Using a double-zeta-d basis for the Ru atom (ru2.atm in the current library), and a k-point sample of 12x12x8 in the bulk Brillouin Zone, the theoretical lattice parameters are a=5.025 bohr and c=7.94 bohr, and we use these values to construct our initial slab as seven unrelaxed layers cleaved from the bulk material.

Important: While a calculation with this input file will run to completion, this input file is flawed and would not be a viable production calculation because of an incorrect use of the Brillouin Zone sampling, which will be discussed in the next example.

Command options

The input file begins with a sequence of instructions to the code about the sort of calculation it will be doing. In this example, we will “do” a “setup” (the iteration independent integrals), and “do iters” (self-consistency), but “no force” and no relaxation. The “setup data” is the instruction to the code to read the data that describes the system. Please do not try to “do iters” with “no setup”!

Setup phase data

This section begins with a title line(s). The keyword “notes” signals that the next line is a text line that describes the nature of the problem. The keyword “end_notes” signals the end of the title section. The notes keyword and notes section are optional, but highly encouraged. Since the entire input file will be echoed into the output listing file, this provides a useful record of the nature of the calculation.

The “dimension” statement is the first required data in the input. We are doing a 2-dimensional slab problem. If we were doing a molecular system, we would enter a “0”, or when we did the bulk hcp Ru problem before this slab calculation, we would have entered a “3” here.

The next required input is “primitive lattice” – supercell vectors. The periodic directions must be in the x-y plane, and the first two lattice vectors define the (1×1) cell for the (0001) surface of Ru. This surface is hexagonal, and we use the bulk theoretical lattice parameter a=5.025, and we start with the bulk interlayer spacing of 7.94 bohr.

In a 2-D slab calculation, the third lattice vector is the non-periodic direction and must be along the z-axis, normal to the slab plane. This lattice vector must be large enough so that the densities from the different periodic images of the slab do not interact. The effective range of the density about any atom is typically between 9-12 bohr (smaller for “hard” atoms like O, larger for metals), so that our lattice vector should be at least our slab thickness (3×7.94=~24 bohr) plus the range of the density on either side (2 x 11.0 used here = ~22) for a total of 46.0 bohr. I’ll choose 48.0 to have a nice “round” number for the grid interval …

The next required input is “grid dimensions”. The code evaluates many of its integrals on a regular space (fft) grid. This regular grid is defined by taking the primitive lattice vectors, and dividing them into grid intervals. Here we use 15×15 in the plane, for a grid spacing of 0.335/point, and 150 grid intervals along the third lattice vector, a spacing of 0.320/point. The more grid intervals, the more accurately the integrals are evaluated, but, also, the more expensive the calculation. A general rule of thumb: you want spacing between points of 0.30-0.40 bohr for most systems, and 0.20-0.30 for finer accuracy in systems containing “hard” atoms. Hard atoms include late first row atoms like N, O, F, or late 3d atoms. Here we have a Ru slab, relatively soft, but later we will be adsorbing CO on the surface, which will require a denser mesh for the same accuracy. For a hexagonal slab, it is a good idea to use a grid dimension in the plane that is divisible by three, so that the various high-symmetry points in the surface (i.e., atop, hcp and fcc hollows) are at the same position with respect to the grid.

The next required input is the number of “atom types”. For the clean slab there is only one type, the Ru. Later in this tutorial, we will add carbon and oxygen for a total of three atom types. The code then requires input to specify the atom file to be used to describe each atom type, and here we are reading the atom potential/basis information from the file “ru.atm”.

Next, the input requires the number of atoms in the unit cell (7), and then a list of each “atom, type, and position”. The arguments for an atom input are the atom index (the order of the atom in the list, #1-#7 here), which type (in the order in which the types were listed above), and the positions of the atoms in Cartesian coordinates (in bohr). Note that the slab is centered about Z=0. The code internally places the Z=0 position at the center of the supercell, so that the outer edges of the supercell is defined as vacuum. The code detects if the atoms and their density impinge on the cell boundary – giving a warning, or even stopping with an error message if the atoms are much to close to a vacuum boundary.

The final required input in the setup section (for a periodic calculation) is definition of the Brillouin Zone sampling. For molecular calculations this is unnecessary – one has only real integrals, or the “gamma-point”, to worry about. In a (periodic) slab calculation, the Schroedinger equation generally needs to be solved for a sampling of points in reciprocal space. It is possible to input the k-point sampling as an explicit list of k-vectors and their weights (using a different keyword/data sequence than shown here), but this example takes advantage of an optional feature to automatically generates a simple k-point sample. The “kgrid ” keyword tells the code to generate a regular mesh in reciprocal space, and here we specify an 6×6 grid in the x-y plane and must use a “0” for the non-periodic direction. A sample with just the gamma-point would be defined by 0x0x0. The kgrid command generates this grid, and offsets this grid from the reciprocal space origin (i.e., the gamma point) by (1/2,1/2,0) in the reciprocal space mesh, i.e., the body-centered position.

The final, required statement in the setup phase section is the “end setup phase” line.

Run phase data

The next section begins with the line “run phase” and ends with the line “end run phase”, and allows you to modify parameters that affect the run (SCF) phase of the calculation. In this example, we modify the convergence criterion to be 0.0005. This criterion does NOT refer to the convergence of the energy in the SCF cycle, but rather refers to the maximum change in an Hamiltonian matrix element (integrals between basis functions over potentials). The value specified here will converge the total energy to roughly 1.d-6 Ry for this problem, but this will vary between different systems (e.g., a covalent system with a big energy gap will converge differently than a metallic system). Note that the convergence criterion is an optional parameter, and that the code will provide a default if you do not override it with this input.

do setup
do iters
no force
no relax
setup data
notes
 Ru(0001) surface: 7-layer slab, unrelaxed cleaved bulk positions
end_notes
dimension of system (0=cluster ... 3=bulk)
 2
primitive lattice vectors
  5.0250     0.000000000    0.00
  2.5125     4.351777654    0.00
  0.0000     0.000000000   48.00
grid dimensions
  15 15  150
atom types
   1
atom file
 ru.atm
number of atoms in unit cell
   7
atom, type, position vector (cleaved bulk positions)
   1    1    2.5125     1.450592551   -11.91
   2    1    0.0000     0.000000000    -7.94
   3    1    2.5125     1.450592551    -3.97
   4    1    0.0000     0.000000000     0.00
   5    1    2.5125     1.450592551     3.97
   6    1    0.0000     0.000000000     7.94
   7    1    2.5125     1.450592551    11.91
kgrid
 12 12 0
end setup phase data
run phase input data
convergence criterion
 0.000500
end run phase data

Slab surface tutorial for SeqQuest: CO/Ru(0001)

Introduction

Units

Input files

Example: A basic input file

Command options

Setup phase data

Run phase data

Hexagonal symmetries and BZ sampling

Picking a density functional (and spin)

Symmetry

Scaling functions I: periodic scaling

Scaling functions II: direction-specific

SCF parameters

Atomic relaxation

Relaxation parameters

Lattice coordinates

Cell optimization

Clean relaxed surface

Adsorbed atoms: (1×1) CO/Ru(0001)

Aliasing atom types

Relaxing atoms

External electric field