Surface Calculations

for lithium version of Exciting


Purpose: In this tutorial you will learn some of the most important steps to simulate a surface with exciting. Simulating a surface does not mean calculating only the air-exposed part of a material, but also considering its "sub-strates", called layers. They are required for an appropriate description of the surface characteristics. Such system, combination of different layers, is called slab. Step by step, you will see here how to pass from bulk to slab calculations. Computationally, slabs are indicated as lxmxn where l,m and n indicate the number of atoms along x, y and z directions. Experiments usually indicate the size of the sample.



STEP 0. Define relevant shell variables and download scripts

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that the relevant shell variables are defined as specified in Tutorial scripts and environment variables.

As a second step, it is necessary to download the scripts and files which are used in this tutorial. Look at the instructions for downloading the script directory in Tutorial scripts and environment variables. Here is a list of the scripts and files which are relevant for this tutorial with a short description.

  • Tutorial-SURFACE.tar.gz: tar.gz archive containg files and directories which are used in this tutorial. You have already donwloaded along with the script!

From now on the symbol $ will indicate the shell prompt.

Preparation

Copy the archive Tutorial-SURFACE.tar.gz to your working directory and extract all the files and directories and enter in the tutorial directory.

$ cp $EXCITINGSCRIPTS/Tutorial-SURFACE.tar.gz ./
$ tar xzf Tutorial-SURFACE.tar.gz
$ cd Tutorial-SURFACE

From now on, Tutorial-SURFACE is identified as main tutorial directory.

You can find 4 different directories, labeled accordingly to the corresponding step in this tutorial.In general you find two kind of subfolders: To-Do where you can work, and others, labeled DONE which contain the most important output of the calculations. To visualize the data you can either make use of the experience gained in other tutorials or use any other program you may want to.


Introduction

You will be carried along the following steps:

  1. Find lattice parameter
  2. Create a slab.
  3. Set vacuum.
  4. Set a good k-space sampling.
  5. Set the number of layers.

In Details

The criteria adopted to pass from one step to the following one, are based on the numerical and the physical consistency of your calculation (e.g., computational cost, energy convergence and electronic structure analysis). The scope is to teach you some aspects you should care before starting any real calculation.

  1. You are going to find the most appropriate bulk lattice parameter (for details refer to Volume optimization for cubic systems) for Copper with two exchange-correlation functionals; choice will be done considering calculate lattice constants.
  2. Once the functional is chosen, you could either try to generate your own surface with ASE as explained in the Atomic Simulation Environment Tutorial, or use the provided input files.
  3. With a small slab, 1x1x3, you will converge the total energy with increasing amount of vacuum. The study of energy convergence, along with DOS and band structure calculations, will guide you to choose an appropriate amount of vacuum, balancing accuracy and efficiency.
  4. Fixed the slab size (still 1x1x3) and the amount of vacuum, you will converge the system with different k-points sampling, to find the most appropriate mesh for further calculations.
  5. Fixed the amount of vacuum and the k sampling on x, y and z, you will see the total energy (per atom) variations with the increasing size of the slab on z (number of layers).

STEP 1. Find the most appropriate bulk lattice parameter

This step is particularly important when you have to define the structure of your system. Nonetheless the definition of a proper lattice parameter have already been discussed in the tutorial on Volume Optimization for Cubic Systems; therefore, if you feel confident enough on such topic, you can skip this first step and pass to the discussion of the results and then to the STEP2.

Understanding the results: LDA vs GGA.

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

Question

Which xctype would you choose for create your slab?

ii) DOS and band structure for bulk (optional)

STEP 2. The surface: Cu(111)

The chosen structure is a Face Centered Cubic (fcc). In such slab you have a sequence ABC that repeats itself. Nonetheless a crystal (bulk structure) can be cut with planes in different directions.

FCCPlanes.GIF

In our case, the Cu bulk has been cut along the (111) plane, and therefore is labeled Cu(111). The numbers (Miller indices) identify the plane in coordinates of Bravais lattice.

111.png

In the (111) surface, the adsorption sites between the atoms of the surface, are not equivalent, for geometrical arrangement of the layers. Therefore 3 is the minimum number of layers to consider.

Building a slab: theoretical consideration

Go to the main tutorial directory and then in the directory STEP2, where files generated by ASE are stored.

$ cd Tutorial-SURFACE
$ cd STEP2

File name indicates the structure and the amount of vacuum considered and its unit (Angstrom). So 113-4Ang.exi indicates 1 atom on x, 1 on y, 3 on z (layers number) and 4 Angstrom of vacuum.

Exercises
  1. Try to obtain the same structure

STEP 3. The surface calculations

Go to the main tutorial directory and then in the directory STEP3

$ cd Tutorial-SURFACE 
$ cd STEP3

Question: Why the slab provided come with vacuum?
Answer To understand the choice some considerations must be made.

In the directory STEP3 you have 3 subdirectories;
To-Do, where you will run your calculations
DONE-221, where your calculations have been already run
DONE-441, where you find more precise calculations.
In the directory To-Do, you can find all input files (.exi) as generated by ASE, and the following files:
script_113-vacuum.bash (executable bash script that edits the files for you)
HEAD (new heading for the input files)
TAIL (new bottom for the input files)
BAND_DOS.txt (part for the input files for band structure and DOS calculations)

Let's see what the script does

str="1 5 8"
for a in $str
do
mkdir 113-$a

The most important variable is $a: the amount of vacuum (in Angstrom). Note: $a refers to the file name. The vacuum has been set with ASE, and therefore is expressed in Angstrom and not in Bohr (as inside the file). Do all these values make sense?

cp HEAD input.xml
tail -12 113-"$a"Ang.exi | head -9 >> input.xml
cat TAIL >> input.xml
mv input.xml 113-$a
cd 113-$a

The script has changed the .exi file, edited, moved to the appropriate directory, previously created and, once in, run exciting.

excitingser

It creates a backup of your input and edits it for DOS and band structure calculations

cp input.xml input.xml.backup
head -16  input.xml > new.input
cat ../BAND_DOS.txt >> new.input
mv new.input input.xml

It starts exciting

excitingser

The density of states and the band structures are treated with the visualization tools

xsltproc $EXCITINGROOT/xml/visualizationtemplates/xmlband2agr.xsl bandstructure.xml 
xsltproc --stringparam ID "t///" $EXCITINGROOT/xml/visualizationtemplates/xmldos2grace.xsl dos.xml >tdos.agr
cd ..
done
done

Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT (see Volume optimization for cubic systems).
Now you can execute the bash script.

$ ./script_113-vacuum.bash &

While you are running the calculations, you can analyze the results in the directories DONE-221 and DONE-441. To do so type

$ cd DONE-221

or
$ cd DONE-441
Exercises
  1. Study the total energy variation as function of vacuum. Such information is contained in the file TOTENERGY.OUT in each directory. Divide it by 3 and compare with the bulk total energy (-1655.03363786 Ha).

In the file DONE-441/total-energy-difference_respect-to-bulk.agr you can see the convergence for 4 {{k}-points (mesh 4 4 1), as shown in the following picture. Your calculations are done with the faster mesh 2 2 1.

  1. Reproduce a similar graphic for your mesh and compare the two.
total-energy-difference_respect-to-bulk.png
Understanding the output

Question: Why 1 Angstrom does not make sense at all?
Answer Try to visualize what you did

$ xcrysden --exciting 113-1/input.xml

The slabs are too near! So near that the distance between the slabs (1 Angstrom) is even inferior to the distance between layers (2.09 Angstrom)! Next time, think before acting and save your computational time.

Message

You don't have to be careful only of what is the meaning of the results coming out from the "box", but also to what you put in. Lobsters do not come out from potatoes!

Criteria to choose the appropriate vacuum

1. Exclusively on total energy consideration, which vacuum size would you consider? Assume 20 Angstrom as an optimal reference. The energy difference per atom between 5 and 20 Angstrom is .001365947 Ha, while between 8 and 20 Angstrom is .00094837 Ha. Could 5 Angstrom be enough?

2. Computational time as function of vacuum. To evaluate it, use the following command in your konsole. Between " you have a string that will be searched in the file INFO.OUT in all directories which name is, at least 113-

$ grep " total                                 :" 113-*/INFO.OUT

you will obtain values in this order (your values may differ)

113-1/INFO.OUT: total                                 :        14.61
113-15/INFO.OUT: total                                 :      1236.13 
113-20/INFO.OUT: total                                 :      2507.53
113-3/INFO.OUT: total                                 :        55.84 
113-5/INFO.OUT: total                                 :       100.30 
113-8/INFO.OUT: total                                 :       256.19

Therefore do no always redirect your grep into a file, be careful! We can easily see that 15 Angstrom needs already more than twice the time required by 8 Angstrom of vacuum. Up to now 5 layers of vacuum seems a good choice.

3. Let's analyze the DOS. In the file all.agr you have the 6 DOS together.

cd DONE-441

If you open the file with xmgrace, you can "navigate" in the various energy ranges. At low energies the differences are negligible, but they increase approaching the Fermi energy. Try to zoom in x between -0.31 and 0.31 and between -0.3 and 350 in y. You will notice the differences. As you can see from the following figure, 5 Angstrom are not enough (green line) to reproduce a 20 Angstrom (violet) like behaviour, while 8 Angstrom is already an acceptable value.
DOS_113_increasing-vacuum.png
Exercises

Compare these results with the one you have obtained with the less precise k mesh.

  1. What do you understand?

You can compare the DOS of the bulk structure with the one we have now obtained.

  1. Which differences do you see? What are they related to?
  2. You can do the same with the band structures
Message

When you start your calculations and set up your system, starting with a poor k mesh is never a good idea!


STEP 4. K space sampling

We have chosen the lattice parameter and 8 Angstrom of vacuum, now we have to optimize the total energy as function of k}-points. Go to the main tutorial directory and then in the directory {{STEP4

$ cd Tutorial-SURFACE
$ cd STEP4/To-Do

In the directory To-Do, copy the following file as a bash script.

str="4 12 16"
for a in $str
do
mkdir $a
cd $a
cat>input.xml<<EOF
<input >
 
  <title>Surface</title>
 
  <structure speciespath="$EXCITINGROOT/species/">
    <crystal>
      <basevect>4.83661030535009 0.00000000000000 0.00000000000000</basevect>
      <basevect>2.41830515267505 4.18862739263879 0.00000000000000</basevect>
      <basevect>0.00000000000000 0.00000000000000 38.13376688639916</basevect>
    </crystal>
    <species chemicalSymbol="Cu" speciesfile="Cu.xml">
      <atom coord="0.33333333333333 0.33333333333333 0.39644149791542"/>
      <atom coord="0.66666666666667 0.66666666666667 0.50000000000000"/>
      <atom coord="0.00000000000000 0.00000000000000 0.60355850208458"/>
    </species>
  </structure>
  <groundstate  ngridk="$a $a 1"
               xctype="GGAPerdew-Burke-Ernzerhof"
               rgkmax="8.0d0">
  </groundstate>
    <properties>
    </properties>
</input>
EOF
excitingser
tail -1 TOTENERGY.OUT > ../0$a.dat
cd ..
done

Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT (see Volume optimization for cubic systems).

Change properties to this file and execute it.

chmod +x nameofyourfile
./nameofyourfile &
Understanding the output

While the job runs, let's analyze the results provided. Move to the directory DONE_113-8_k.

$ cd STEP4/DONE_113-8_k

The subdirectories contain the output of calculations with different k-meshes. The total-energy difference respect to the calculation obtained with 32x32x1 mesh (102 k-points) is shown in the following image. Which is the most appropriate choice?
convergence_with_k_points.png

Question
  • Why there is only 1 k-point on the direction perpendicular to the surface?
  • The k-space sampling in x and y does not change increasing the number of layers, therefore we can pass to the last step of this tutorial.

STEP 5. Increasing the number of layers

Up to now, we have found the most appropriate lattice parameter, we have defined the vacuum size and studied the convergence respect to the k-ponts mesh. Last step of this tutorial is to understand the criteria to determine the number of layers.

Go to the main tutorial directory and then in the directory STEP5

$ cd Tutorial-SURFACE
$ cd STEP5
Question Are 3 layers enough?

Answer No, let' s see why.

The first step is to consider that we have an fcc structure, therefore we want to repeat ABC as many times possible. Nonetheless an approximation is required to save computational time. Sometimes ABCA or ABCAB could be acceptable approximations or not. This last part of the tutorial focuses on this choice. In the directory STEP5 you find 2 directories: To-Do and DONE. You can choose an input from the second directory, copy in the first one and run it.

Understanding the output

In DONE you find calculations already completed for slab with 3, 4, 5, 6 and 9 layers (directories 113 114 115 116 and 119 respectively) with 8 Angstrom of vacuum with 1 (directory 111) and 4 k-points (directory 441). You can visualize the total energy per atom and the computational time for these runs. In the following images the total energy convergence with number of layers is plotted for 1 and 4 k-points.

totalenergy_numberoflayers_GAMMA.png
totalenergy_numberoflayers_4K.png

What would you understand from the first graph? ABCA (4 layers) is not an acceptable approximation, while the total energy suggests that the ABCAB sequence (5 layers) is already acceptable. Nonetheless the second graph shows that increasing the number of k-points the previous considerations are not valid anymore. Therefore the number of k-points is an important parameter to set during our calculations when you want to establish the number of layers for your calculations.

Nonetheless this is not the only criteria, otherwise we could think to use a lot of k-points and layers from the beginning. Hereby I stress once more on the computational cost is important to balance efficiency and reliability of your results. As shown in the following image, the time does not scale linearly with the system size, therefore your choice must be compatible with your computational resources.

times_441.png

Resume

In this tutorial we have seen the following necessary steps to define a slab calculations:
1) find the appropriate lattice parameter with bulk calculations
2) build your surface with ASE
3) with a minimum number of layers define the vacuum size
4) increase the number of
k-points
5) increase the number of layers

For each step, remember to look at the same time at:
a) Density of States and Band Structures
b) Total Energy Convergence
c) Computational Time

… and good luck !


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