MATSE 598 HW #6: VASP, Part II November 22, 1999 Due: 7 December, 1999.

1 Introduction
In this last of our VASP exercises we will explore some of the more advanced features in the code. We will do this via a simple example in which we calculate the surface properties of the (111) plane of Al. Some of the topics we will touch on include: K-point convergence issues Finite-temperature DFT Geometry optimizations Analyzing charge densities and electrostatic potentials Densities of States You will nd the les you need to begin these exercises in ~dsiegel/vasp2.tar.gz. Copy this le to your directory and uncompress it and untar it as you did in the last exercise. Finally, we will be making extensive use of the VASP manual in these exercises. It can be found on the web at: http://tph.tuwien.ac.at/~vasp/.

2 Planewave Convergence
Using the supplied POTCAR, POSCAR bulk, and KPOINTS les, construct an appropriate INCAR le and conduct a planewave convergence test. Use the same convergence criteria as before: 1-2 meV per atom. Determine whether you need to set the ISPIN tag for a material susch as Al, and explain your reasoning. Generally speaking, one never needs to run a spin polarized calculation when conducting a planewave convergence test. This is because planewave convergence depends only on the pseudopotential, which is not altered by the inclusion of spin polarization.

3 K-point Convergence Testing

Once you have determined the appropriate setting for ENCUT, you can now begin the second phase of convergence testing. This involves calculating the total energy per atom as a function of the number of k-points in the irreducible wedge of the Brillouin zone. As we saw in class, this procedure involves setting up a regular grid in reciprocal space which is then reduced by symmetry operations to a few \special" k-points. As a user our task in testing for k-points convergence is relativly easy. Since VASP takes care of all the symemtry operations, all we must do is incrementally increase the size of the grid until convergence is attained. (However, it may also be necessary in some circumstances to investigate the e ect of altering the origin of the grid.) 1

3.1 Bulk
ISMEAR = 1 SIGMA = 0.1

We will begin our testing on the bulk Al system. First o , you must replace the ISMEAR INCAR with the following:

= -5

line in your

See section 7.33 in the manual for more details on what these settings do; we will have more to say about them later. Secondly, you need to edit the KPOINTS le. Before you do this read section 5.5.2 in the manual. We will be conducting two independent k-points tests|one for the origin of the grid taken to be at the -point, and one for the origin shifted o the -point. These are the only two important shifts. To switch to the -centered grid change line 3 in the KPOINTS le to read Gamma-centered. (Obviously, the other grid is obtained by saying Monkhorst-Pack). Using one of the two choices for grid origin, and beginning with the smallest possible grid (1 1 1), incrementally increase the values on line 4 until you get convergence to within 1-2 meV. That is, increase the dimensions of the mesh along each of the reciprocal lattice vectors in unit increments: 111 ! 222 ! 333 (1) By looking in the OUTCAR le one can identify the explicit coordinates of the k-points generated by the given mesh. Look for a series of lines with the following format:
Subroutine IBZKPT returns following result: =========================================== Found 110 irreducible k-points:

Run a separate k-point convergence test for each choice of grid origin, and report your results in tabular form. Each table must include the size of the grid, the number of irreducible k-points (generated by the grid), and the total energy per atom. Indicate which origin and grid density gives a converged result. Do you notice any di erence in convergence speed between the odd-numbered grids vs. even-numbered grids? Is the choice of grid origin signi cant? Especially for the case of metals, calculating an accurate ground state energy when using a relatively small number of k-points can be numerically di cult. This is due to the fact that metals have partially lled bands near the Fermi Energy, which require a large number of k-points to accurately determine the Fermi surface. One common way around this problem is discussed in sections 7.33 and 9.4 of the manual, and involves broadening the 0 K Fermi function into a smooth function. The trick is to determine how much broadening is necessary, since too little broadening requires a large number of k-points, and too much broadening leads to incorrect results. In this section we investigate this e ect. Using your converged set of k-points, calculate the size of the dE term reported in the OSZICAR as you change the value of SIGMA in the INCAR. I suggest starting from a value of 0.1 and increasing in 0.1 increments until you get to 0.5. Identify the largest SIGMA value that gives a dE per atom that is less than 1 meV. Use this SIGMA for the remainder of these exercises. Report the size of dE per atom as a function of SIGMA in tabular form. Now that you have a converged k-point set for the bulk system, we can begin another convergence test for the Al (111) surface. (Generally, we must perform a new k-point convergence test for each new supercell we investigate.) Since our slab is bulk-like in the directions parallel to the surface, we can simply use the same 2

3.2 Finite-temperature DFT

3.3 Surface

origin and grid density along these directions as was determined for the bulk system. For example, if you found that a -centered grid with dimensions 8 8 8 was su cient for the bulk, you can re-use the rst two grid densities and must only check for convergence in grid points along the long axis (which in our case is z ). Run an additional k-point test, starting with a grid density given by b1 b1 1, where b1 represents the optimal bulk grid density. Tabulate your results like you did in the previous exercise. For most (all?) surfaces it turns out that only 1 grid division is necessary along the long axis. Explain why is this is so.

4 Vacuum Thickness Convergence

As we did for the Jellium calculations|and as must be done for any surface calculation|we must check to see that the amount of vacuum in our slab supercell is su cient to prohibit interactions between surfaces as a result of using periodic boundary conditions. Starting from a vacuum thickness of 6A, calculate the total slab energy as a function of vacuum thinkness in increments of 2Auntil you reach 14A. Convergence is attained when the total energy is stable to within 10 meV. Use the converged vacuum thinkness for the remainder of your calculations.

5 Surface Relaxation
When any surface is created by cleaving from the bulk material there will be a rearrangement of atoms near the surface. The extent of this rearrangement depends upon the geometry of the surface and type of bonding. One of the goals of this exercise is to illustrate how to calculate an accurate surface energy. In order to do this we must take account of these surface rearrangements|commonly called reconstructions or relaxations |by allowing the surface atoms to relax to their minimum energy positions. This can be accomplished by instructing VASP to perform a geometry optimization on our Al(111) surface. Performing a geometry optimization calculation requires the addition of several tags to the INCAR. For our purposes the following additions are su cient: POTIM = 0.5 Sets the scaling constant for the forces. See 7.20. NSW = 25 Sets the max number of atomic moves. See 7.17. IBRION = 2 Use the conjugate-gradients method for the minimization. See 7.19. EDIFFG = -0.03 Convergence criteria for the forces in eV/A. See 7.16. MAXMIX = 40 Instructs the charge density mixer to re-use information from previous iterations. See 7.40. Once you have set these tags in the INCAR you are set to begin your relaxation calculation. VASP will start and begin the rst series of SCF cycles to calculate the energy of the starting structure. Once it has this, it calculates the forces on each atom using the Hellman-Feynmann theorem, and updates the positions of the atoms in accordance with these forces. VASP then begins a new series of SCF cycles to calculate the electronic groundstate for the new atomic positions. Forces are again calculated, the atoms are moved, and the whole procedure repeats until the magnutude of the forces on all atoms is less than the setting of EDIFFG, namely 0.03 eV/A. During the course of the run you can monitor the forces on the atoms by looking for the lines:
POSITION TOTAL-FORCE (eV/Angst) -----------------------------------------------------------------------------------

At the end of the run you will have the relaxed structure and its corresponding energy. The last force statement in the OUTCAR contains the nal positions of the atoms. You can also get this information from the CONTCAR, although in this le the atomic positions are reported in direct coordinates. Make a table of the percentage change in interlayer spacing (relative to the bulk spacing) between the atoms as a function of distance into the slab. 3

Using the nal relaxed energies from the surface calculation Esurf , and the optimal bulk energy Ebulk , calculate the surface energy according to the following formula: = 21 (Esurf Ebulk ); (2) A where A is the surface area of the slab. (Note the factor of 2 since we have two surfaces.) Compare your

6 Surface Energy

answer with the experimental value.

7 Charge Density Pro le

In the next two sections we will calculate quantities similar to those that we analyzed for the Jellium model, namely the planar averaged charge density and electrostatic potential. By default, VASP writes out the charge on a regular grid for use in plotting. There are two les which contain this information, the CHGCAR and the CHG les. The shape of the grid follows the symmetry of the supercell, with cubic supercells having cubic grids, and hexagonal supercells having hexagonal grids, etc. The size of the grid along each lattice vector of the supercell is equal to the dimensions of the ne FFT grid, which is identi ed by the labels NGXF, NGYF, NGZF in the OUTCAR. These dimensions are also listed in the CHGCAR and CHG les after the atomic coordinates. Write a brief program to read in either the CHGCAR or CHG les and perform a planar average of the charge density as a function of position perpendicular to the surface. Note that you will need to divide the raw charge from the CHGCAR or CHG les by the supercell volume to get a charge density. It does not matter which le you choose to read. Read sections 5.10 and/or 5.11 in the manual for more information. Plot your planar-averaged charge density for both the bulk and surface systems as a function of position perpendicular to the surface. Indicate the positions of the atoms in your plots.

8 Workfunction Calculation
As you did for the charge, make a planar averaged plot of the electrostatic potential as a function of position perpencicular to the surface. You will need to set the tag LVTOT = .TRUE. in the INCAR, and read the potential from the LOCPOT le. Make plots for both the bulk and surface systems, identifying the positions of the atoms and the Fermi energies. (The Fermi energy is given in the OUTCAR). Make one additional plot in which the minimum of the potential corrugations for the surface system are aligned with those from the bulk. The workfunction is then given by the di erence between the bulk Fermi energy and the top (or vacuum level) of the electrostatic potential for the slab. (See Figure 1 for an example plot.) What value do you calculate? What is the experimental value?

9 Extra Credit: Layer-Projected Density of States

Read Sections 5.15 and 7.30 of the manual. Choose an appropriate RWIGS setting for Al and plot the density of states as a function of ion position perpendicular to the surface. Compare with what you see for the bulk system.