From Bulk to Surface: how to

Purpose: In this tutorial you will learn some of the most important steps to simulate a surface (from now onw slab) with exciting. Step by step you will see how to pass from bulk to a surface.
You will be carried along the following steps:

  1. find the most appropriate lattice parameter (not discussed in details) for the next step;
  2. with the chosen lattice parameter, you will create your surface;
  3. with a minimal number of layers you will converge the minimal vacuum required
  4. you will converge increasing the k space sampling
  5. converge the number of layers
  6. learn how to relax

How
The physical consistency or your results and the reasons of the adopted choices (number or layers, k space sampling etc) will be given studying total energy converge and variation in the density of states (DOS) and in the band structure.

First you are going to find the lattice parameters with two exchange correlation functionals, and look for such values at their band structure and density of states.
Adopting GGA functional you can try to generate your own surface with ASE as from this tutorial. Nonetheless the structures will be provided in the file downloadable here.
Nonetheless the 1x1x3 surface will be provided for step 3. and 4.
Studying the convergence of total energy with increasing amount of vacuum, along with DOS and band structure calculations, you will find the minimum vacuum amount required for the next step: total energy convergence as function of k points.
Fixed the amount of vacuum and the k sampling on x, y and z, increasing the slab size on z (number of layers) you will see the total energy variations with the increasing number of layers.


Table of Contents

STEP 1. Set up the calculation
Consider an fcc transition metal, you can choose between Au, Cu, Ag, Pt etc, it's your choice. My pick is Cu and therefore from now on everything system specific is related to copper. Nonetheless the changes required for the others are minimal.
At this stage you should know what to edit in the input.xml.
Nonetheless, here you have a simple bash script that with two variables, $a and $c changes the lattice parameter

str2="3.3 3.38 3.43 3.5"
for a in $str2;
do
...

and the exchange correlation functionals…

str="GGAPerdew-Burke-Ernzerhof LDA-X-alpha"
for c in $str;
do

creates the appropriate directories and entering into

mkdir $a-$c
cd $a-$c

creates the input file

cat > input.xml<<EOF
<input>
  <title>Copper: Bulk Calculations</title>

Here the fcc base vectors are set

  <structure speciespath="/users/sol/giulio/EXCITING/exciting/species">
    <crystal scale="$a">
      <basevect> 1.0     1.0     0.0 </basevect>
      <basevect> 0.0     1.0     1.0 </basevect>
      <basevect> 1.0     0.0     1.0</basevect>
    </crystal>
    <species speciesfile="Cu.xml" >

along with the atom coordinate

      <atom coord="   0.0  0.0  0.0 "/>
    </species>
  </structure>

.. the K-space grid, and the rgmax. The exchange correlation funcionational is chosen by xctype here set equal to $c (variable set above.

  <groundstate ngridk="10 10 10"
               xctype="$c"
               rgkmax="8.0d0">
  </groundstate>
    <properties>
    </properties>
</input>
EOF

here exciting is started, the last total energy value copied in an appropriate file

/users/sol/giulio/EXCITING/exciting/bin/excitingser 
tail -1 TOTENERGY.OUT >> ../$a-$c.dat
cd ..

done
done
grep "  " *alpha.dat >energy.latticeLDA-X-alpha
grep "  " *rhof.dat >energy.latticeGGAPerdew-Burke-Ernzerhof

Now you can visualize the two files energy.latticeLDA-X-alpha and energy.latticeGGAPerdew-Burke-Ernzerhof with gnuplot or xmgrace.
You can even take a look at the files with

more energy.lattice*

You find that the minimum for LDA is at 3.38 while for GGA is 3.43.


optional

If you want you can repeat the previous exercise (but in anther directory), choosing one exchange-correlation functional and increasing the kpoint mesh

  <groundstate ngridk="10 10 10"

with something higher, for example

  <groundstate ngridk="20 20 20"

You can use this script:

str2="3.0 3.3 3.4 3.43 3.46 3.5 3.6"
for a in $str2;
do
str="10 16  24"
for c in $str;
do
mkdir $a-$c
cd $a-$c
cat > input.xml<<EOF
<input>
  <title>Copper</title>
<structure speciespath="/users/sol/giulio/EXCITING/exciting/species">
    <crystal scale="$a">
      <basevect> 1.0     1.0     0.0 </basevect>
      <basevect> 0.0     1.0     1.0 </basevect>
      <basevect> 1.0     0.0     1.0</basevect>
    </crystal>
    <species speciesfile="Cu.xml" >
      <atom coord="   0.0  0.0  0.0 "/>
    </species>
  </structure>
  <groundstate ngridk="$c $c $c"
               xctype="GGAPerdew-Burke-Ernzerhof"
               rgkmax="8.0d0">
  </groundstate>
    <properties>
    </properties>
</input>
EOF

/users/sol/giulio/git/exciting2/bin/excitingser
tail -1 TOTENERGY.OUT >> ../$a-$c.dat
cd ..

done
done
grep "  " *-10.dat >energy.lattice_10
grep "  " *-16.dat >energy.lattice_16
grep "  " *-24.dat >energy.lattice_24

Now you can compare the runtime (file INFO.OUT) with the previous cases for the minimum 3.43. How does the computational time scale? How much is the total energy difference in the various cases at the minimum?

Computationally: the smallest grid is the cheapest. But cheap is not always good. Quality and quantity must find an average. Experience and in particular the specific problem you are considering will guide you.


STEP 1a DOS AND BAND STRUCTURE

QUESTION: Which xctype would you choose for further calculation?
Answer In our case we choose the lattice parameter calculated with xctype="GGAPerdew-Burke-Ernzerhof ", 6.86 Bohr (3.63 Angstrom).
Why GGA and not LDA? Experimental value is 3.615 Angstrom, and LDA says 6.76 (3.58 Angstrom). This underestimation is relevant.
If you are considering a more complex systems, as the adsorption of large molecules, the matching between substrate and adsorbate is extremely important: what, experimentally, is alligned, it may not be so in the theoretical calculation.

Ok, no lets's go now into the subfolder 3.42-GGAPerdew-Burke-Ernzerhof /
We must edit the file input.xml with the following parts

  <groundstate do="skip" ngridk="10 10 10"
               xctype="GGAPerdew-Burke-Ernzerhof "
               rgkmax="8.0d0">

.. we are telling exciting to look at the restart files do="skip" and to skip the ground state calculation that already exists.

Then we tell the kpoint path for the band structure…

    <properties>
        <bandstructure>
            <plot1d>
                <path steps="100">
                    <point coord="0.00000   0.00000   0.00000" label="GAMMA"></point>
                    <point coord="0.50000   0.50000   0.00000" label="X"></point>
                    <point coord="0.50000   0.75000   0.25000" label="W"></point>
                    <point coord="0.00000   0.50000   0.00000" label="L"></point>
                    <point coord="0.00000   0.00000   0.00000" label="GAMMA"></point>
                </path>
            </plot1d>

.. and the bandstructure (in Bohr)

        </bandstructure>
        <dos nsmdos="0" winddos="-5.67d0 5.67d0"/>
    </properties>

we run again
$EXCITINGDIR/bin/excitingser
and here we have the output files.
Now you can generate the DOS and the band structure, typing in the console

xsltproc /users/sol/giulio/EXCITING/exciting/*/*/xmlband2agr.xsl bandstructure.xml >xmlgracebandstructure.agr
xsltproc --stringparam ID "t///" /users/sol/giulio/EXCITING/exciting/xml/visualizationtemplates/xmldos2grace.xsl dos.xml >tdos.agr

Store both DOS and bands.


STEP 2. Why the surface: (111)

I have chosen an fcc (111) because in such surface you have a sequence ABC that repeats itself. A fcc adsorption site is different from a hcp. Therefore 3 is the minimum amount of layers to consider/
When we create a slab, in our case a simple 1x1xn, we are reproducing the sequence ABC n times. If we consider a 1x1x6 slab, we have 2 surfaces, the top and the bottom layers. Therefore in such case, we have 2 slabs of 3 layers each. Then, for examples, all computed adsorption energies must be divided by factor 2. The internal layers are considered bulk like.
This is necessary because in reality, going under the surface, the layers become more or less slowly similar to a bulk structure. The thickness (or the number of layers) differs by element, and for each chemical specie by the crystallographic direction considered.
Logically a 2 layers slab will never really resemble the properties of a surface + bulk layer, thefore, ideally, the more layers we can consider, the better! Practically it has a cost and therefore the number of layers must be minimized.

The surface calculations
In the file provided here, you can find all

STEP 5
You have now to consider that meanwhile you have your slab, you need to add some vacuum. But vacuum has some computational cost that must be reduced to minimum.
There is not only a mere time problem, but also a physical one. On one side the supercell apporoach is useful and allow us to consider infinite materials (we have done that with the bulk calculations), but in a surface, while useful for the in place interactions, on the direction vertical to the surface generate another surface. A top layer faces a bottom one of the upper unit cell. This may bring a lort of undesired effects.
Then before relaxing anything on a surface or doing whatever we must be sure to have
1) appropriate sample of kspace (now on z, with vacuum we set 1)
2) a minimum amount of vacuum

STEP6
We well see the impotance of the point 2 with the help of the DOS and the bandstructure.
Here you have n different slab wth different amount of vacuum.

run the calculation and evaluate the total energy per atom as function of the vacuum thickness.
Now look at the BANDstuctures and the DOS and how they changes according to the vacuum
Remember when you set the vacuum, do it in units of layers distances.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License