Part 1: A nickel surface¶
By the end of this tutorial, you will be able to:
- use the selective dynamic mode to relax specific layers with the residual minimization method with direct inversion in the iterative subspace (RMM-DIIS)
- set the $\mathbf{k}$ points and unit cell to model surfaces
- compute the surface energy of a monoatomic system
1.1 Task¶
Perform a spin-polarized calculation to relax two of five layers representing a nickel (100) surface using RMM-DIIS. Determine the relaxation energy by comparing the surface energy and the bulk calculation.
To state the obvious, surfaces are the frontier of any realistic bulk system. It is possible to model surfaces in VASP by extending the structure's unit cell along one direction, such that one part is filled with the bulk material and the other part is empty to represent the vacuum. Such a geometry is usually called a surface slab or just a slab. Only one $\mathbf{k}$ point is considered along the extended direction because any additional $\mathbf{k}$ point describes the interaction with neighboring unit cells. This also implies that both parts, the material, and the vacuum region, must be sufficiently thick to simulate bulk properties adequately in the middle of the slab. In other words, along the elongated direction, you cannot rely on periodic boundary conditions; Instead, you must annul their effects to model a surface.
1.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e01_Ni100-relaxation. Check them out!
fcc Ni (100) surface
3.53
.50000 .50000 .00000
-.50000 .50000 .00000
.00000 .00000 5.00000
5
Selective Dynamics
Cartesian
.00000 .00000 .00000 F F F
.00000 .50000 .50000 F F F
.00000 .00000 1.00000 F F F
.00000 .50000 1.50000 T T T
.00000 .00000 2.00000 T T T
fcc Ni
3.53
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
1
Cartesian
0 0 0
SYSTEM = clean Ni (100) surface
ISTART = 0
ICHARG = 2
ENCUT = 270
ISMEAR = 2
SIGMA = 0.2
ALGO = Fast
EDIFF = 1E-6
ISPIN = 2
MAGMOM = 5*1
NSW = 100
POTIM = 0.8
IBRION = 1
ISIF = 2
SYSTEM = fcc Ni
ISTART = 0
ICHARG = 2
ENCUT = 270
ISMEAR = 2
SIGMA = 0.2
ALGO = Fast
EDIFF = 1E-6
ISPIN = 2
MAGMOM = 1
LORBIT = 11
K-Points
0
Monkhorst Pack
9 9 1
0 0 0
K-Points
0
Monkhorst Pack
9 9 9
0 0 0
Pseudopotential of Ni.
The POSCAR file defines a structure with a lattice constant of 3.53 Å. In the bulk case, only one atom exists in the unit cell. In the surface calculation, how many Ni layers are defined? How much vacuum is between the Ni surfaces? Which layers are allowed to relax their ionic positions?
Click to see the answer!
There are five layers, with one atom per layer. Then, there are 3 $\times$ 3.53 Å = 10.59 Å vacuum. The selective dynamics mode allows two layers to relax and holds three layers fixed.
Make comments after each tag in the INCAR file! This calculation takes spin-polarization into account on the level of collinear magnetism within the framework of spin-density-functional theory (SDFT), see ISPIN = 2 and MAGMOM = 1. In other words, the interaction between spin components is effectively considered by separately computing the Kohn–Sham orbitals for both spin components for an effective Hamiltonian taking both spin components' charge densities into account. The initial charge density is obtained from overlapping atoms, see ICHARG and the occupation function is approximated using the Methfessel-Paxton method [Phys. Rev. B, 1989, 40, 3616]. The goal is to perform a geometry relaxation, see IBRION, NSW, POTIM and ISIF. That is, to find the ionic positions with zero forces and stress acting on them. On the level of pseudocode, the performed task reads
- given the ionic position, the electronic ground state is computed within SDFT,
- forces and stress acting on the ions are computed as the expectation value of the gradient of the electronic Hamiltonian (Hellmann-Feyman theorem),
- ions are relaxed to their instantaneous ground state using the RMM-DIIS,
- repeat 1.-3. until the convergence criterion is met.
In the KPOINTS file, an equally spaced $\mathbf{k}$ mesh is defined using the Monkhorst-Pack method. The odd number of $\mathbf{k}$ points results in a Γ-centered $\mathbf{k}$ mesh. For the surface calculation, we define only one $\mathbf{k}$ point normal to the surface direction. Explain why!
1.3 Calculation¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e01_*
mpirun -np 2 vasp_std
cd $TUTORIALS/surface/e01_*/bulk
mpirun -np 2 vasp_std
Visualize the bulk and the surface structure using py4vasp! Passing 3 to the plot function allows plotting a $3 \times 3 \times 3$ grid of the unit cell. You can measure distances in the figure by right-clicking the starting point and right-clicking the endpoint twice.
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_Ni100-relaxation/bulk")
mycalc.structure.plot(3)
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_Ni100-relaxation")
mycalc.structure.plot(3)
What is the relative change of the relaxed surface layers compared to their bulk position in the z direction?
Click to see the answer!
fcc Ni (100) surface
3.53000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
-0.5000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.0000000000000000
Ni
5
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.5000000000000000 0.5000000000000000 0.1000000000000014 F F F
0.0000000000000000 0.0000000000000000 0.2000000000000028 F F F
0.5000000000000000 0.5000000000000000 0.3002101608664644 T T T
0.0000000000000000 0.0000000000000000 0.3955075210077033 T T T
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
Inward relaxation of surface layers: $$ \Delta d_{{12}} = ((0.3955075210077033-0.3002101608664644)-0.1)/0.1*100=-4.70\%, $$ $$ \Delta d_{{23}} = ((0.3002101608664644-0.2000000000000028)-0.1)/0.1*100=+0.21\%. $$
Use py4vasp to check the convergence of the surface calculation!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e01_Ni100-relaxation")
mycalc.energy[:].plot()
Click to see the answer!
First ionic step:
FORCES acting on ions
electron-ion (+dipol) ewald-force non-local-force convergence-correction
-----------------------------------------------------------------------------------------------
-.551E-14 0.577E-13 0.285E+04 0.909E-14 0.234E-14 -.285E+04 -.856E-23 -.137E-23 -.248E+01 -.260E-13 -.131E-12 0.280E-02
-.484E-14 -.331E-13 0.145E+04 -.112E-13 -.202E-13 -.145E+04 -.854E-23 -.139E-23 -.125E+00 0.265E-13 -.720E-13 0.688E-02
0.129E-14 0.824E-14 0.196E-02 0.116E-13 0.872E-14 -.719E-10 -.139E-23 0.139E-23 -.548E-02 0.122E-13 -.227E-13 -.685E-03
0.133E-13 -.135E-12 -.145E+04 -.764E-15 -.531E-14 0.145E+04 0.229E-24 0.308E-23 0.105E+00 0.253E-13 0.432E-13 -.706E-02
-.259E-13 -.610E-13 -.285E+04 -.879E-14 -.131E-13 0.285E+04 -.620E-23 -.869E-23 0.249E+01 -.516E-13 0.578E-13 -.231E-02
-----------------------------------------------------------------------------------------------
-.217E-13 -.163E-12 0.579E-02 -.934E-16 -.275E-13 0.000E+00 -.245E-22 -.697E-23 -.135E-01 -.136E-13 -.125E-12 -.370E-03
POSITION TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
0.00000 0.00000 0.00000 0.000000 -0.000000 0.419195
0.00000 1.76500 1.76500 -0.000000 -0.000000 -0.401602
0.00000 0.00000 3.53000 0.000000 0.000000 -0.002587
0.00000 1.76500 5.29500 -0.000000 -0.000000 0.392846
0.00000 0.00000 7.06000 0.000000 0.000000 -0.407853
-----------------------------------------------------------------------------------
total drift: -0.000000 -0.000000 -0.008063
In the first step the forces on layer 1 and 2 are comparable to layer 4 and 5. The selective dynamics mode only allows layer 4 and 5 to relax.
In the 5th and final step, we obtain:
FORCES acting on ions
electron-ion (+dipol) ewald-force non-local-force convergence-correction
-----------------------------------------------------------------------------------------------
-.526E-14 0.109E-13 0.286E+04 0.907E-14 0.232E-14 -.286E+04 -.468E-25 -.437E-25 -.248E+01 0.640E-14 0.161E-14 0.483E-02
0.254E-13 -.337E-13 0.147E+04 -.120E-13 -.206E-13 -.147E+04 0.178E-25 0.806E-26 -.146E+00 0.213E-14 0.152E-13 0.325E-02
-.139E-13 -.154E-13 0.131E+02 0.478E-15 0.126E-13 -.131E+02 -.598E-27 0.598E-27 0.447E-01 -.441E-13 -.285E-13 -.125E-01
-.217E-12 0.775E-13 -.145E+04 0.173E-14 -.121E-14 0.145E+04 0.384E-25 0.650E-25 0.819E-01 0.542E-14 0.233E-13 -.941E-02
-.380E-14 0.600E-13 -.289E+04 0.620E-15 -.207E-13 0.289E+04 0.302E-26 0.345E-26 0.235E+01 -.321E-13 -.475E-13 0.842E-02
-----------------------------------------------------------------------------------------------
-.214E-12 0.993E-13 0.140E+00 -.934E-16 -.275E-13 -.455E-12 0.118E-25 0.334E-25 -.154E+00 -.622E-13 -.359E-13 -.543E-02
POSITION TOTAL-FORCE (eV/Angst)
-----------------------------------------------------------------------------------
0.00000 0.00000 0.00000 -0.000000 -0.000000 0.425897
0.00000 1.76500 1.76500 -0.000000 -0.000000 -0.402956
0.00000 0.00000 3.53000 0.000000 -0.000000 -0.022978
0.00000 1.76500 5.29871 0.000000 0.000000 -0.001056
0.00000 0.00000 6.98071 -0.000000 -0.000000 0.001092
-----------------------------------------------------------------------------------
total drift: -0.000000 0.000000 -0.020002
Here, the force on layer 4 and 5 are very small!
EDIFFG = 10 $\times$ EDIFF = 1e-5
Finally, compute the surface energy! Follow the link to learn how.
Click to see the answer!
The unrelaxed surface energy yields: $$ \sigma^\mathrm{unrel} = \frac{1}{2A} \left( E_\mathrm{slab}^\mathrm{unrel} - N E_\mathrm{bulk} \right) = \frac{1}{2\cdot3.53^2/2} \left[ -25.5567 - 5 \cdot \left(-5.4681\right) \right] = 143.0 \,\mathrm{meV/A^2} $$ where we used the total energy of the unrelaxed slab $E_\mathrm{slab}^\mathrm{unrel}=-25.5567\, \mathrm{eV}$.
The average surface energy is based on the calculation with relaxed top two layers on one side and yields: $$ \sigma^\mathrm{avg} = \frac{1}{2A} \left( E_\mathrm{slab} - N E_\mathrm{bulk} \right) = 141.5 \,\mathrm{meV/A^2} $$
The approximated relaxed surface energy finally yields $$ \sigma^\mathrm{rel} = 2\sigma^\mathrm{avg} - \sigma^\mathrm{unrel} = 140.1 \,\mathrm{meV/A^2} $$
The relaxation energy is quite small, at $~3 \,\mathrm{meV/A^2}$, but some materials relax a lot more than a close-packed metal like fcc Nickel, and the methodology, in that case, remains the same. Literature gives a PBE surface energy of $138 \,\mathrm{meV/A^2}$ for 100 Ni [Sci. Data, 2016, 3, 160080], so the partial relaxation method provided a good result with minimal computational cost.
Note that the expression for the surface energy
$$\sigma = \frac{E_\mathrm{slab} - N E_\mathrm{bulk}}{2A}$$
is only valid for symmetric and stoichiometric slabs. The factor of 1/2 in the formula is there because we always create 2 surfaces when we cleave a bulk. The formula only works if these are equivalent, i.e., the slab is symmetric. The multiplicator $N$ signifies that the slab contains $N$ times the number of formula units of the bulk. This is only true if the slab is stoichiometric, i.e., the ratio of its constituents is equal to the ratio in the bulk.
What does it mean that the surface energy is positive?
Click to see the answer!
The atoms energetically prefer the bulk crystal environment; thus, it costs energy to form a surface. Semi-local exchange-correlation functionals within DFT tend to predict surfaces to be more stable than they are in the experiment. That is, the surface energy is often smaller in the simulation than in the experiment.
And what is the interpretation if the (111) surface has a smaller energy than the (211) surface?
Click to see the answer!
This means that creating the (211) surface costs more energy. Therefore, it is more difficult to make the cut along that facet.
By the end of this tutorial, you will be able to:
- compute local quantities, such as the local charge and local spin magnetization
- compute and plot the local density of states (DOS)
2.1 Task¶
Compute the charge on each nickel atom and the local density of states for a nickel (100) surface.
To compute the charge on each nickel atom, we must ask what part of the charge density belongs to which atom. Or, more generally, to compute local quantities, we need to access localized Kohn-Sham (KS) orbitals.
One way is to specify a volume that is considered to belong to a specific atom. That is we define the Wigner–Seitz radius that is set with the RWIGS tag and gets a default value in the POTCAR file. It is used for the DOS when LORBIT<10. The problem with this approach is that the volume is ambiguous. It is preferable, to define a localized basis indexed by local quantum numbers, such as the angular momentum $l$, magnetic quantum number $m$ and spin $\sigma$.
A different approach is to perform a unitary transformation, usually of a subset of KS orbitals, to construct a set of localized states. This is the basis of Wannier function based methods, such as the maximally-localized Wannier functions. In practice, such a unitary transformation is difficult to find, when KS bands of different characters are strongly entangled.
Finally, there is a way to use projector operators that project out only KS states with a desired character. This is the basis of projected localized orbitals [Phys. Rev. B, 2008, 77, 205112], and is used for the DOS when LORBIT>10. The projectors and phase factors are written into a file called PROCAR.
2.2 Input¶
The input files are in $TUTORIALS/surface/e02_Ni100-DOS. Check them out!
POSCAR from Example 1
fcc Ni (100) surface
3.53000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
-0.5000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.0000000000000000
Ni
5
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.5000000000000000 0.5000000000000000 0.1000000000000014 F F F
0.0000000000000000 0.0000000000000000 0.2000000000000028 F F F
0.5000000000000000 0.5000000000000000 0.3002101608664644 T T T
0.0000000000000000 0.0000000000000000 0.3955075210077033 T T T
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
SYSTEM = clean (100) Ni surface
ENCUT = 270
ISMEAR = -5
ALGO = Normal
ISPIN = 2
MAGMOM = 5*1
LORBIT = 11
K-Points
0
Monkhorst-Pack
13 13 1
0 0 0
Pseudopotential of Ni.
Click to see the answer!
ISMEAR = -5 calls the tetrahedron method with Blöchl corrections. Mind that this method is not variational with respect to partial occupancies and yields erroneous forces and stress for systems with partially filled bands. Therefore, it should not be used when performing a geometry optimization in metals! LORBIT = 11 allows the DOSCAR and PROCAR files to be written. The PROCAR file contains the partial DOS for KS states of a specific character, i.e., for a specific site and quantum numbers $l$, $m$. The partial DOS is projected out using the projectors which are written to the PROCAR file.
2.3 Calculation¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e02_*
mpirun -np 2 vasp_std
Plot the DOS of Ni for each layer using py4vasp! Executing the Python code will show the documentation of the py4vasp DOS plot function. Read it and change the code to plot a specific layer! The projected DOS near the surface should show a larger exchange splitting compared to the most bulk-like layer, i.e., layer 3.
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e02_Ni100-DOS")
#Plot the DOS for layer 3 (most bulklike) and 5 (relaxed surface)
mycalc.dos.plot(selection="3, 5")
Open the OUTCAR file and find information about the local charge and magnetization! Is the local charge and local magnetization reduced or enhanced towards the surface layer?
Click to see the answer!
At the end of the OUTCAR file information on the local charge and magnetization is given.
total charge
# of ion s p d tot
------------------------------------------
1 0.466 0.328 8.315 9.109
2 0.488 0.483 8.311 9.283
3 0.492 0.486 8.314 9.291
4 0.499 0.505 8.318 9.322
5 0.477 0.349 8.325 9.152
--------------------------------------------------
tot 2.422 2.152 41.582 46.156
magnetization (x)
# of ion s p d tot
------------------------------------------
1 -0.003 -0.020 0.751 0.727
2 -0.007 -0.025 0.668 0.635
3 -0.007 -0.026 0.675 0.642
4 -0.008 -0.026 0.669 0.636
5 -0.004 -0.021 0.744 0.720
--------------------------------------------------
tot -0.029 -0.117 3.506 3.360
The surface layers 1 and 5 have a reduced charge density, but enhanced magnetization.
The on-site spin magnetization can also be visualized using py4vasp. Plot the structure and compare the size of the arrows on each site!
import py4vasp
mycalc = py4vasp.Calculation.from_path("./e02_Ni100-DOS")
mycalc.local_moment.plot()
Note there are difficulties visualizing the arrows in html, so we have included an image here:
Change the INCAR file to read LORBIT = 1 and RWIGS = 1.286. Then run a fresh VASP calculation! Does the computed DOS, local charge, and local magnetization change?
Click to see the answer!
volume of typ 1: 40.5 %
total charge
# of ion s p d tot
------------------------------------------
1 0.410 0.283 8.254 8.947
2 0.430 0.417 8.249 9.095
3 0.433 0.419 8.250 9.102
4 0.439 0.435 8.255 9.129
5 0.421 0.301 8.262 8.984
--------------------------------------------------
tot 2.13 1.85 41.27 45.26
magnetization (x)
# of ion s p d tot
------------------------------------------
1 -0.002 -0.017 0.757 0.738
2 -0.005 -0.021 0.675 0.648
3 -0.005 -0.022 0.682 0.655
4 -0.005 -0.022 0.675 0.648
5 -0.002 -0.018 0.750 0.731
--------------------------------------------------
tot -0.02 -0.10 3.54 3.42
2.4 Questions¶
- What strategies exist to compute local quantities within the projector-augmented wave method?
- What does the tag RWIGS control?
- What information is written into the PROCAR file?
By the end of this tutorial, you will be able to:
- compare the Wigner–Seitz cell of a bulk and surface calculation
- plot the band character
3.1 Task¶
Compute the band structure along Γ-X-M-Γ of a nickel (100) surface.
Recall that for a surface calculation of some compound, the crystal structure is elongated along one direction in real space and some padding with vacuum is added. This has dramatic effects on the Brillouin zone defined in reciprocal space. In the bottom row of the figure, you can see the first Brillouin zone, the so-called Wigner-Seitz cell, of an fcc structure. This corresponds to the bulk nickel computed in Example 1. You can check this using the tool SeeK-path by uploading the POSCAR file of the bulk calculation used in Example 1!
You can use the same tool, to obtain the Wigner-Seitz cell of the structure of the surface calculation. You can see that the Wigner-Seitz cell is squashed along one direction and the labels of the high-symmetry points have changed. Since the lattice vectors defining the Bravais lattice have changed, you must be extra careful when identifying the high-symmetry points of the bulk and surface calculations with each other. In the figure above, this is done for you. The top row shows the high-symmetry points in the Wigner-Seitz cell of the surface calculation.
3.2 Input¶
Check out the input at $TUTORIALS/surface/e03_Ni100-band!
fcc Ni (100) surface
3.53000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
-0.5000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.0000000000000000
Ni
5
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.5000000000000000 0.5000000000000000 0.1000000000000014 F F F
0.0000000000000000 0.0000000000000000 0.2000000000000028 F F F
0.5000000000000000 0.5000000000000000 0.3002101608664644 T T T
0.0000000000000000 0.0000000000000000 0.3955075210077033 T T T
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
SYSTEM = clean (100) nickel surface
ICHARG = 11
ENCUT = 270
ISMEAR = 2
SIGMA = 0.2
ALGO = Normal
ISPIN = 2
MAGMOM = 5*1
LORBIT = 11
kpoints for band structure Γ-X-M-Γ
15
line-mode
reciprocal
0.0000000000 0.0000000000 0.0000000000 Γ
0.0000000000 0.5000000000 0.0000000000 X
0.0000000000 0.5000000000 0.0000000000 X
0.5000000000 0.5000000000 0.0000000000 M
0.5000000000 0.5000000000 0.0000000000 M
0.0000000000 0.0000000000 0.0000000000 Γ
Pseudopotential of Ni.
CHGCAR from Example 1
Check the meaning of all tags used in the INCAR file and make comments in the file! What does ICHARG = 11 mean?
3.3 Calculation¶
Open a terminal, navigate to this example's directory, copy the charge density from Example 1, and run VASP by entering the following:
cd $TUTORIALS/surface/e03_*
cp ../e01_*/CHGCAR .
mpirun -np 2 vasp_std
Check that, in the OUTCAR file and the standard output to the terminal, it states that the charge density is not adjusted during the calculation!
Click to see the message you should find!
...
Static calculation
charge density remains constant during run
spin polarized calculation
...
Plot the band structure using py4vasp! Try to specify a quantum number $l$ to see which bands have which character.
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e03_Ni100-band")
# Plotting the bandstructure and visualizing band character
mycalc.band.plot(selection="Ni(s, p, d)")
3.4 Questions¶
- How does the Wigner–Seitz cell of a surface calculation change compared to a bulk calculation?
- What is the character of a band?
By the end of this tutorial, you will be able to:
- explain the purpose and application of the Methfessel-Paxton method
- use py4vasp to get the angle or distance
4.1 Task¶
Perform a non-spin-polarized calculation to relax two of five layers representing a nickel (111) surface using RMM-DIIS and determine the surface energy.
Many quantities, including the charge density computed in each electronic iteration, require numerical integration over the first Brillouin zone. In a metallic system, we must integrate a discontinuous function, because the occupation at $T=0$ is a step function. If the same band is occupied at some $\mathbf{k}$ point and empty at another $\mathbf{k}$ point, these are partially filled bands and are called conduction bands. For systems featuring conduction bands, numerical integration over the first Brillouin zone requires a very fine $\mathbf{k}$ mesh, which takes huge computational effort. The different strategies to face this problem are controlled by the ISMEAR tag. One method to approximate the occupation function is the Methfessel-Paxton method. The scheme is based on smooth approximants to the step functions constructed to give the exact result when integrating polynomials of a prescribed order. This order is set by ISMEAR $> N$.
4.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e04_Ni111-relaxation. Check them out!
Ni(111)
3.53000000000000
0.7071067811865476 0.0000000000000000 0.0000000000000000
-0.3535533905932738 0.6123724356957946 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.7735026918962570
Ni
5
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.3333333333333333 0.6666666666666666 0.1000000000000000 F F F
0.6666666666666666 0.3333333333333333 0.2000000000000000 F F F
0.0000000000000000 0.0000000000000000 0.3000000000000000 T T T
0.3333333333333333 0.6666666666666666 0.4000000000000000 T T T
bulk/POSCAR
fcc Ni
3.53000000000
0.0000000000 0.5000000000 0.5000000000
0.5000000000 0.0000000000 0.5000000000
0.5000000000 0.5000000000 0.0000000000
1
Cartesian
0.000000000 0.000000000 0.000000000
SYSTEM = clean (111) surface
ISTART = 0
ICHARG = 2
ENCUT = 270
ISMEAR = 2
SIGMA = 0.15
ALGO = Fast
EDIFF = 1E-6
NBANDS = 34
LMAXMIX = 6
NSW = 20
POTIM = 0.8
IBRION = 1
bulk/INCAR
SYSTEM = fcc Ni
ISTART = 0
ICHARG = 2
ENCUT = 270
ISMEAR = 2
SIGMA = 0.15
ALGO = Fast
EDIFF = 1E-6
LMAXMIX= 6
K-Points
0
Monkhorst-Pack
11 11 1
0 0 0
bulk/KPOINTS
K-Points
0
Monkhorst Pack
11 11 11
0 0 0
Pseudopotential of Ni.
Click to see the answer!
There are again five layers with one atom per layer. Then, there are $$ (1-0.4)\cdot 5.774\cdot 3.53Å \approx 12.23Å $$ of vacuum. The selective dynamics mode allows two layers to relax and holds three layers fixed.
4.3 Calculation¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e04_*
mpirun -np 2 vasp_std
cd $TUTORIALS/surface/e04_*/bulk
mpirun -np 2 vasp_std
Then, plot the structure using py4vasp and determine the area of the surface! To take the distance of two atoms right-click one of them once, and the other twice. To measure the angle between atoms, you need to right-click three atoms and the last one twice.
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e04_Ni111-relaxation" )
mycalc.structure.plot()
What is the surface area for the (111) surface?
Click to see the answer!
Either you know the area of a parallelepiped with side length $a = 2.5Å$ to be $$ (2.5 Å)^2 \cdot \sin(120/360*\pi) = 5.4126587737 Å^2, $$ or take twice the area of an equilateral triangle $A = (a^2\sqrt{3})/4$. The Ni 111 surface is more compact than the Ni 100 surface.
Which surface area remains closer to the bulk structure? The 100 or the 111 surface? How does this relate to the surface energy? Compute the surface energy of both the unrelaxed and the relaxed structure to find out!
Click to see the answer!
As the forces are already rather small at the beginning for the (111) surface, it remains closer to its bulk structure. We can follow the same procedure to compute the surface energies as in example 1 for the unrelated surface energy $\sigma^\mathrm{unrel}$: $$ \sigma^\mathrm{unrel} = \frac{1}{2A} \left( E_\mathrm{slab}^\mathrm{unrel} - N E_\mathrm{bulk} \right) = 122.7 \,\mathrm{meV/A^2} \quad , $$
the average surface energy $\langle\sigma\rangle$ after relaxation of the top two layers,
$$ \langle\sigma\rangle = \frac{1}{2A} \left( E_\mathrm{slab} - N E_\mathrm{bulk} \right) = 122.1 \,\mathrm{meV/A^2} \quad , $$
and the fully relaxed surface energy $\sigma^\mathrm{rel}$: $$ \sigma^\mathrm{rel} = 2\langle\sigma\rangle - \sigma^\mathrm{unrel} = 121.6 \,\mathrm{meV/A^2} \quad . $$
The gain of the energy per area due to relaxation is even smaller for the (111) surface than for the (001) surface due to the lower interlayer distance and the densely packed surface of the (111) slab. (111) has the lowest surface energy of all fcc Ni surfaces. Again the error we make compared to the literature result of $\sim 120 \,\mathrm{meV/A^2}$ is minimal, even though we disregard magnetism here [Sci. Data, 2016, 3, 160080].
Disclaimer: The calculations here are examplatory to show specific steps in the workflow. Therefore, some comments are in order:
It is always necessary to perform convergence studies with respect to k-mesh density (KPOINTS), cutoff energy (ENCUT), and for surface calculations also the vacuum thickness and the slab thickness.
In the OUTCAR file you can find a rather large value for the
external pressure. This is Pulay stress that can be avoided by setting a higher cutoff energy (ENCUT). Mind that it is absolutely necessary to use the same value for ENCUT in order to compare the total energy of two calculations. Thus, all calculations in this notebook should be done at a higher cutoff energy, to compute the surface energy without Pulay stress.Furthermore, we took the bulk material as a given. As the lattice constants vary depending on the exchange-correlation functional, it is expedient to perform a structure relaxation for the bulk material as a preparatory step.
4.4 Questions¶
- What problem does the Methfessel-Paxton method address?
- How can you get the distance between two atoms from a structure plot in py4vasp?
- How can you verify that your surface slab is thick enough for an accurate surface energy calculation?
Good job! You have finished Part 1!¶
Get started with Part 2.
In the file browser, navigate to $TUTORIALS/surface and open part2.ipynb!