Part 2: A CO molecule adsorbed on a nickel surface¶
By the end of this tutorial, you will be able to:
- relax a molecule on a surface with 100% coverage.
- apply and understand dipole corrections.
- evaluate geometry changes in a slab-adsorbate system.
5.1 Task¶
Perform a relaxation of a CO molecule and two of five layers representing a nickel (111) surface using RMM-DIIS. Determine the geometry changes compared to a molecule in a box.
Adsorption of molecules on surfaces is a common and complex problem in surface science. The configuration space is large, even for simple molecules, and the resulting high-dimensional potential-energy surface can be flat. Generally, the molecule can absorb in different orientations, and equilibrium bond angles can change due to interaction and charge transfer to/from the surface. CO adsorption is a common test case for catalyst reactivity because hazardous CO often needs to be oxidized to CO2 [New J. Chem., 2023, 47, 20222-20247].
The adsorbate may significantly change the charge distribution of the slab, and as a result, the system develops a finite dipole moment. For periodic boundary conditions, this finite dipole moment results in a gradient in the potential across the vacuum region. However, the aim of the calculation is to simulate a single isolated slab-adsorbate system without interaction with periodic images. To annul the effects of the undesired electric field, we must introduce so-called dipole corrections. These correct the potential, energies, and forces and are recommended for any system with broken inversion symmetry, i.e., allowed dipole moment.
5.2 Input¶
5.2.1 Input for CO on Ni[111]¶
The input files to run this example are prepared at $TUTORIALS/surface/e05_CO_on_Ni111_relaxation. Check them out!
from ase import Atoms
from ase.build import fcc111, add_adsorbate, bulk
from ase.io import read, write
adsorbate = Atoms('CO')
adsorbate[1].z = 1.1421714441457986
adsorbate.center(vacuum=4.0)
write("./e05_CO_on_Ni111_relaxation/CO/POSCAR",adsorbate, direct=True)
slab = fcc111('Ni', (1, 1, 5), a=3.53)
add_adsorbate(slab, adsorbate, 1.8, 'ontop')
slab.center(vacuum=4.0, axis=2)
write("./e05_CO_on_Ni111_relaxation/POSCAR",slab, direct=True)
There are different adsorption sites on the surface. Here, we could conveniently add the adsorbant to the so-called "ontop" configuration, marked with a red dot in the figure. It is the most stable adsorption site for high coverage [Surf. Sci., 2003, 526, 332-340]. High coverage refers to the fact that the surface is completely covered in CO molecules. To simulate smaller coverage, we would need a supercell. Refer to the ASE documentation on surfaces for an overview of other high-symmetry adsorption sites.
Open the POSCAR files and add the selective dynamics such that the CO molecule and two surface layers can relax!
Ni C O
1.0000000000000000
2.4960869375885126 0.0000000000000000 0.0000000000000000
1.2480434687942563 2.1616746980061543 0.0000000000000000
0.0000000000000000 0.0000000000000000 19.0943572451033177
Ni C O
5 1 1
Selective dynamics
Direct
-0.3333333333333333 0.6666666666666666 0.2094859726700561 F F F
0.0000000000000000 0.0000000000000000 0.3162215084138439 F F F
0.3333333333333334 0.3333333333333333 0.4229570441576317 F F F
-0.3333333333333333 0.6666666666666666 0.5296925799014195 T T T
0.0000000000000000 0.0000000000000000 0.6364281156452074 T T T
0.0000000000000000 0.0000000000000000 0.7306968033467327 T T T
0.0000000000000000 0.0000000000000000 0.7905140273299439 T T T
The DIPOL should be set to the center of charge. As this is not known a priory, a good estimate is to use the center of mass instead. Compute the center of mass!
slab = read("./e05_CO_on_Ni111_relaxation/POSCAR")
print("DIPOL = 0.5 0.5 ", slab.get_center_of_mass(scaled=True)[2])
SYSTEM = CO adsorption on Ni(111)
ISTART = 0
ENCUT = 400
ISMEAR = 0
SIGMA = 0.15
ALGO = Fast
EDIFF = 1E-4
AMIX = 0.4
LMAXMIX = 6
SYSTEM = CO adsorption on Ni(111)
ISTART = 1
# electronic minimization
ENCUT = 400
ISMEAR = 0
SIGMA = 0.15
ALGO = Fast
EDIFF = 1E-6
AMIX = 0.4
LMAXMIX = 6
# structure relaxation w RMM-DIIS
NSW = 20
POTIM = 0.8
IBRION = 1
EDIFFG = -0.05
# dipole correction
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = 0.5 0.5 0.452747064661102
K-Points
0
Monkhorst Pack
11 11 1
0 0 0
Pseudopotentials of Ni, C, and O.
There are two INCAR files for the slab-absorbate system: INCAR.1 to obtain a charge distribution for the system without dipole corrections, and INCAR.relax that restarts the calculation by reading the WAVECAR (ISTART=1) and performs a structure relaxation (IBRION, etc.) with dipole corrections (LDIPOL, IDIPOL, and DIPOL). Why is it better to apply the dipole corrections after a preparatory electronic minimization?
Click to see the answer!
The charge distribution of the periodic system will yield a certain dipole moment. To suppress the long-range interaction of the dipole moments, a discontinuity is introduced in the potential to create a saw touth. The discontinuity depends on the size of the dipole moment and thus introducing it early on can slow down the convergence of the Kohn-Sham eigenvalues. Without the correction, however, one yields a $1/L^2$ dependence of the total energy on the vacuum thickness $L$, thus, inhibiting the convergence with respect to vacuum thickness and effectively simulating a "bulk of layered slabs".
Click to see the answer!
The Oxygen requires a rather large cutoff as you can see from the value of ENMAX in the POTCAR file. Here, we even choose a cheap pseudopotential. For small molecules, i.e., CO, we recommend testing if hard potentials, i.e., C_h and O_h, lead to different results. Generally, short bond lengths are better described by hard potentials; however, these are computationally more expensive. Refer to the VASP Wiki for recommendations on choosing pseudopotentials!
The Ni (111) surface is magnetic, so it would be better to perform a magnetic calculation (ISPIN=2 or LNONCOLLINEAR=True). As CO does not have a magnetic moment, the nonmagnetic calculation may already be a good approximation. So, to demonstrate the procedure, we perform a nonmagnetic calculation here.
Long unit cells with broken inversion symmetry are prone to charge sloshing. This issue is addressed by the charge-density mixing. So, in case of convergence issues, please check the various parameters and options for the charge-density mixing!
If there are tags you are unfamiliar with, refer to the VASP Wiki and add comments after the tags in the INCAR file!
The KPOINTS files define 1 k point in the directions vacuum padding is added and a fine k mesh within the slab.
5.2.2 Input for CO molecule¶
Let us also have a look at the inputs of the CO molecule
co = read("./e05_CO_on_Ni111_relaxation/CO/POSCAR")
print("DIPOL = 0.5 0.5 ", co.get_center_of_mass(scaled=True)[2])
CO/POSCAR
C O
1.0000000000000000
8.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 8.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 9.1421714441457986
C O
1 1
Direct
0.5000000000000000 0.5000000000000000 0.4375328142157524 F F T
0.5000000000000000 0.5000000000000000 0.5624671857842477 F F T
CO/INCAR
SYSTEM = CO dimer in a box
ISMEAR = 0 ! Gaussian smearing
EDIFF = 1E-6
ALGO = Fast
ENCUT = 400
SIGMA = 0.15
AMIX = 0.4
# dipole correction
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = 0.5 0.5 0.509
# structure relaxation w conjugate gradient
IBRION = 2
POTIM = 0.1
EDIFFG = -0.05
NSW = 10
CO/KPOINTS
Gamma-point only
0
Monkhorst Pack
1 1 1
0 0 0
CO/POTCAR
Pseudopotentials of C and O.
5.3 Calculation¶
Make sure you added the selective dynamics in the POSCAR! Then, open a Terminal, navigate to this example's directory, and run VASP for the CO on Ni(111). This involves first converging to a charge distribution and, then, a structure relaxation with an appropriate dipole correction.
cd $TUTORIALS/surface/e05_*
cp INCAR.1 INCAR
mpirun -np 2 vasp_std
cp INCAR.relax INCAR
mpirun -np 2 vasp_std
You can find the final dipole moment and the energy corrections in the OUTCAR. Search for dipolmoment and dipol+quadrupol energy correction! Are these large?
Click to see the answer!
No, it is small, as expected for a metallic slab with a rather unpolar adsorbent. The dipole moment is much larger in a system with a charged ion, e.g., Potassium or radical Oxygen, as an adsorbent.
Vizualize the final structure using py4vasp! Click the top most Ni atom and then twice the Ni underneath. Did the distance change compared to the clean slab in Example 4?
import py4vasp
CO_on_Ni = py4vasp.Calculation.from_path( "./e05_CO_on_Ni111_relaxation")
CO_on_Ni.structure.plot((2, 2, 1))
The structure relaxation in Example 4 was done with a lower cutoff energy ENCUT. Is this a problem when comparing the interlayer distances?
Click to see the answer!
The cutoff energy ENCUT must be the same when comparing total energies of different calculations/structures. However, in case of geometric parameters, e.g., interlayer distances, it is not strictly necessary. It is sufficient if the quantity of interest is converged w.r.t. the cutoff energy for each calculation individually.
Compute the relative change of the interlayer distance w.r.t. the clean surface. In other words, compute the inward or outward relaxation the Ni surface experiences due to the CO molecule.
import py4vasp
def layer_relaxation(s1,s2):
d_clean = s1.cartesian_positions()[4] - s1.cartesian_positions()[3]
d_final = s2.cartesian_positions()[4] - s2.cartesian_positions()[3]
return 100*(d_final[2]-d_clean[2])/d_clean[2]
CO_on_Ni = py4vasp.Calculation.from_path("./e05_CO_on_Ni111_relaxation")
clean_Ni = py4vasp.Calculation.from_path("./e04_Ni111-relaxation")
d_CO_Ni = CO_on_Ni.structure.cartesian_positions()[5] - CO_on_Ni.structure.cartesian_positions()[4]
print("Distance CO - Ni: ",d_CO_Ni[2], " Angstrom.")
print("The two top most Ni layers, relax outwards by: ")
print(layer_relaxation(clean_Ni.structure,CO_on_Ni.structure), "%")
Let us analyze how the CO molecule changes when adsorbed to the surface compared to the isolated molecule. First, run the ionic minimization for the isolated molecule by entering the following in the Terminal:
cd $TUTORIALS/surface/e05_*/CO
mpirun -np 2 vasp_std
In the case of the isolated CO molecule, we immediately introduce dipole corrections. In fact, many calculations converge even if the dipole corrections are introduced from the beginning. This is particularly common if the molecule is not strongly polar, i.e., has a small dipole moment in the charge distribution.
Compute the change of the bond length!
import py4vasp
CO_on_Ni = py4vasp.Calculation.from_path( "./e05_CO_on_Ni111_relaxation")
CO = py4vasp.Calculation.from_path( "./e05_CO_on_Ni111_relaxation/CO")
# For the adsorbed molecule, the indices are 5 for C and 6 for O,
# while they are 0 and 1, respectively, for the molecule in the box.
d_CO_adsorbed = CO_on_Ni.structure.cartesian_positions()[6] - CO_on_Ni.structure.cartesian_positions()[5]
d_CO_free = CO.structure.cartesian_positions()[1] - CO.structure.cartesian_positions()[0]
d_CO_rel = 100*(d_CO_adsorbed[2] - d_CO_free[2])/d_CO_free[2]
print("Adsorbed CO: ", d_CO_adsorbed[2], " Angstrom.")
print("Isolated CO: ", d_CO_free[2], " Angstrom.")
print("Changes by ", d_CO_rel, "%")
Click to see the answer!
The bond length of the adsorbed CO is increased by about 1%.
5.4 Questions¶
- How many symmetrically inequivalent, high-symmetry adsorption sites are possible on an fcc (111) surface?
- What is the effect of a dipole moment of the charge density on the potential in the vacuum?
- Name four quantities with respect to which you need to perform a convergence study for a molecule adsorbed to a slab!
- Does the cutoff energy need to be the same in order to compare the structure, e.g., interlayer distances, of the clean slab and the adsorbate-slab system?
By the end of this tutorial, you will be able to:
- compute adsorption energies of molecules on surfaces
- compute the work function of surfaces
- visualize the average local potential
6.1 Task¶
Compute the adsorption energy of a CO molecule adsorbed on a nickel (111) surface in the on-top position.
The adsorption energy $E_\mathrm{ads}$ of a molecule on a surface is calculated as
$$ E_\mathrm{ads} = E_\mathrm{total} - E_\mathrm{slab} - E_\mathrm{molecule}, $$
where $E_\mathrm{total}$, $E_\mathrm{slab}$, and $E_\mathrm{molecule}$ correspond to the total energy of the slab with the adsorbed molecule, the clean slab, and the isolated molecule, respectively. The structure of these systems should be relaxed. If $E_\mathrm{ads}>0$ the molecule tends to be adsorbed to the surface. The value of $E_\mathrm{ads}$ depend on various factors, e.g., the adsoption site, the coverage, the surface facet, the chemical decomposition, etc.
Read about computing the work function on the VASP Wiki!
6.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e06_CO_on_Ni111_adsorption. Check them out!
We computed the clean slab in Example 4 and the total system and isolated CO in Example 5, but these calculations use different values for the plane-wave cutoff (ENCUT). While many properties of interest are well-converged for the individual calculations, for the comparison of total energies, it is paramount to manually set the ENCUT tag to the same value for all calculations! Execute the following to copy the CONTCAR of Example 4 to POSCAR for this example! Then, we can use the same value for ENCUT as in Example 5 to calculate the total energy of the clean slab und reuse the total energies computed in Example 5.
! cp ./e04_Ni111-relaxation/CONTCAR ./e06_CO_on_Ni111_adsorption/POSCAR
POSCAR CONTCAR from Example 4 Part 1
Ni(111)
3.5299999999999998
0.7071067811865476 0.0000000000000000 0.0000000000000000
-0.3535533905932738 0.6123724356957946 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.7735026918962573
Ni
5
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.3333333333333357 0.6666666666666643 0.1000000000000014 F F F
0.6666666666666643 0.3333333333333357 0.2000000000000028 F F F
-0.0000000000000000 0.0000000000000000 0.2987324140345219 T T T
0.3333333333333357 0.6666666666666643 0.3969599428265009 T T T
SYSTEM = clean nickel (111) surface
ENCUT = 400
ISMEAR = 0
SIGMA = 0.15
ALGO = Fast
AMIX = 0.4
LMAXMIX = 6
EDIFF = 1E-6
NBANDS = 34
# work function
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = 0.5 0.5 0.199
WRT_POTENTIAL= hartree ionic
LVHAR = .TRUE.
K-Points
0
Monkhorst-Pack
11 11 1
0 0 0
Pseudopotential of Ni.
Check the meaning of LVHAR, WRT_POTENTIAL!
Why do we apply a dipole correction in this calculation?
Click to see the answer!
It is very unlikely that charge will accumulate on a pure metal surface unless some external perturbation, e.g., an adsorbent or an applied field, causes this accumulation. Nevertheless, the computation of the work function relies on a well-defined potential in the vacuum region. It performs a line average in the chosen direction and tries to identify the vacuum region by determining where the potential is constant. A (small but finite) dipole moment can cause issues for the search algorithm.
Use the code below to find the center of mass and set DIPOL appropriately!
from ase.io import read
ni111 = read("./e06_CO_on_Ni111_adsorption/POSCAR")
ni111.get_center_of_mass(scaled=True)[2]
Click to see the answer!
DIPOL = 0.5 0.5 0.1990837
6.3 Calculation¶
Open the Terminal, navigate to this example's directory, and run the calculation by entering the following:
cd $TUTORIALS/surface/e06_*
mpirun -np 2 vasp_std
Calculate the adsorption energy using the calculations for the total system and the isolated molecule from the previous exercise, as well as the clean surface calculation that just ran. How large is the adsorption energy?
import py4vasp
total_calc = py4vasp.Calculation.from_path( "./e05_CO_on_Ni111_relaxation/")
molecule_calc = py4vasp.Calculation.from_path( "./e05_CO_on_Ni111_relaxation/CO/")
clean_slab_calc = py4vasp.Calculation.from_path( "./e06_CO_on_Ni111_adsorption/")
#Get the total energies of the last steps
E_slab = clean_slab_calc.energy.to_numpy('without_entropy')
E_mol = molecule_calc.energy.to_numpy('without_entropy')
E_tot = total_calc.energy.to_numpy('without_entropy')
E_ads = E_tot - E_mol - E_slab
print("Adsorption energy of CO on Ni (111): ", E_ads*1000, " meV")
LDA strongly overestimates the adsorption energy. Semilocal xc functionals, i.e., GGA, slightly improve the prediction but still overestimate the adsorption energy [Phys. Rev. B 59, 7413–7421 (1999)]. It is curious that semilocal xc functionals underestimate the cost to form a surface, but overestimate the adsorption energy.
Let us look at the work function $\boldsymbol{\phi}$. Detailed step-by-step instructions and best practices on computing the work function are available on the VASP Wiki. Plot the work function using py4vasp!
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e06_CO_on_Ni111_adsorption/")
workfunction_value = mycalc.workfunction.to_dict()['vacuum_potential'][1] - mycalc.workfunction.to_dict()['fermi_energy']
print(mycalc.workfunction)
print("The work function is: ", workfunction_value, " eV")
mycalc.workfunction.plot()
Why do we set LVHAR = T and not LVTOT = T. Or in other words, why is the exchange-correlation (xc) potential not included in the potential to evaluate the vacuum level for the work function? Why do we get two slightly different vacuum levels on either side of the slab?
Click to see the answer!
The evaluation of the exchange-correlation potential in the vacuum region is problematic for numerical reasons. But from a physical point of view, it is clear that the contribution from the xc potential to the vacuum potential should be vanishing. Therefore, it is more accurate to compute the vacuum potential without the exchange-correlation potential. Nevertheless, all the exchange-correlation effects are accounted for since the self-consistent field and the Fermi energy fully include the exchange-correlation potential.
A production calculation using thicker slabs and more accurate settings reports $\boldsymbol{\phi}$ = 5.10 eV [Surf. Sci., 2019, 687, 332-340].
The values of the vacuum potential are slightly different on the two sides of the slab because the slab is not symmetric. The lower side is unrelaxed, while the upper is relaxed.
6.4 Questions¶
- When should dipole-related tags be used in surface calculations?
- Should the exchange-correlation potential be included when determining the vacuum level of the potential? Why?
By the end of this tutorial, you will be able to:
- compute the partial DOS of a molecule on a surface.
- plot the planar average of the potential for the molecule on the surface and analyze the changes to the work function.
7.1 Task¶
Compute the partial DOS and work function of the Ni(111) CO adsorbate system to identify the bonding and anti-bonding states.
The isolated CO atom form a $\sigma$ bond between the $sp_z$-hybrid orbital of C and the $p_z$ orbital of O. In addition, the $p_x$ and $p_y$ orbitals of C and O combine to two $\pi$ bonds. Finally, there is a pair of electrons in the O $2s$ orbital and a pair of electrons in a C $sp_z$-hybrid orbital. The so-called CO $5\sigma$ is the highest occupied molecular orbital (HOMO) and the CO $2\pi^*$ is the lowest unoccupied molecular orbital (LUMO). The HOMO and LUMO interact with the metal surface to form bonding and antibonding CO–metal hybrid orbitals. As the antibonding CO $5\sigma$–metal hybrid orbital is above the Fermi level, it corresponds to a charge donation to the surface. Similarly, the bonding $2\pi^*$–metal hybrid orbital become populated which is refered to as back-donation.
Above explanation corresponds to the Blyholder model [J. Phys. Chem. 68, 2772–2777 (1964).] of CO bonding to metal surfaces. Our aim is to observe this charge transfer in the partial DOS.
7.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e07_CO_on_Ni111_pDOS. Check them out!
In Example 5, we performed the relaxation for the adsorbed CO on Ni(111). We will restart the calculation from that structure and charge density, but increase the k mesh and add some INCAR tags to obtain a the partial DOS.
! cp ./e05_CO_on_Ni111_relaxation/CONTCAR ./e07_CO_on_Ni111_pDOS/POSCAR
CO on Ni(111)
3.53000000000000
0.7071067811865475 0.0000000000000000 0.0000000000000000
-0.3535533905932737 0.6123724356957946 0.0000000000000000
0.0000000000000000 0.0000000000000000 11.5470053837925146
Ni C O
5 1 1
Selective dynamics
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 F F F
0.3333333333333357 0.6666666666666643 0.0499999999999972 F F F
0.6666666666666643 0.3333333333333357 0.1000000000000014 F F F
0.0000000000000000 0.0000000000000000 0.1495895086893959 T T T
0.3333333333333357 0.6666666666666643 0.1998803094841342 T T T
0.3333333333333357 0.6666666666666643 0.2430926282202450 T T T
0.3333333333333357 0.6666666666666643 0.2714058495025823 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
0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00
SYSTEM = CO adsorption on Ni(111)
ICHARG = 1
# electronic minimization
ENCUT = 400
ISMEAR = -5
ALGO = Fast
EDIFF = 1E-6
AMIX = 0.4
LMAXMIX = 6
# partial DOS
NEDOS = 5001
EMIN = -10
EMAX = 10
LORBIT = 11
NBANDS = 48
# dipole correction
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = 0.5 0.5 0.452747064661102
LVHAR = .TRUE.
K-Points
0
Gamma
13 13 1
0 0 0
Pseudopotential of Ni, C, O
The structure defined in the POSCAR file is the result written to the CONTCAR file of Example 5.
The mesh in the KPOINTS file was increased to generate a smoother DOS.
We have previously used all tags in the INCAR file, please make sure you remember them all from previous examples. Nevertheless, you should review the LORBIT and ISMEAR tags and understand why we use the specific values in this example.
Click to see the answer!
LORBIT = 11 uses the PAW projectors to decompose the Kohn-Sham orbitals into local quantum numbers. It can be used to quickly find good estimates for the character of the site-projected DOS. ISMEAR = -5 selects the tetrahedron method with Blöchl corrections for setting the partial occupancies of the Kohn-Sham orbitals and is the recommended setting for highly accurate total-energy calculations and calculations of the DOS.
7.3 Calculation¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e07_*
cp $TUTORIALS/surface/e05_CO_on_Ni111_relaxation/CHGCAR $TUTORIALS/surface/e07_CO_on_Ni111_pDOS/.
mpirun -np 4 vasp_std
It is not possible to plot the molecular orbitals directly via the simple projection scheme employed via LORBIT=11. However, using a linear combination of atomic orbitals (LCAO) within molecular orbital (MO) theory to create an MO diagram yields an understanding of which atomic orbitals contribute to which bond. Construct the MO diagram for the CO molecule and plot the contributions to the CO $5\sigma$ state!
import py4vasp
from scipy import signal
win = signal.windows.hann(30)
def smooth(dos):
return signal.convolve(dos, win, mode='same') / sum(win)
mycalc = py4vasp.Calculation.from_path( "./e07_CO_on_Ni111_pDOS/")
dos = mycalc.dos.read("C(pz) O(pz) C(s) 5(d)")
en = mycalc.dos.read()['energies']
c_spz = dos["C_pz"] + dos["C_s"]
o_pz = dos["O_pz"]
surface = dos["Ni_5_d"]
# optionally smooth the surface DOS
ssurface = smooth(surface)
(py4vasp.plot(en,c_spz, "Cspz", xlabel="Energy (eV)", ylabel="DOS 1/eV")
+py4vasp.plot(en,o_pz, "Opz")
+py4vasp.plot(en,ssurface, "surface Ni d"))
You can see that the CO $5\sigma$ state gets split and forms a bonding and an anti-bonding state with the surface. There is a tiny peak above the fermi surface which means $5\sigma$ is not fully occupied when adsorbed to the surface. It corresponds to a charge transfer to the surface.
Above we used a smoothing function on the surface DOS. The k-mesh density would need to be increased to an extremely fine density in order to get a smooth DOS. Here, we use a smoothing based on a simple convolution. This helps to recognize the features in the DOS, but mind that this may also introduce artefacts or reduces peaks that are actually present in the true system. Hence, applying a convolution should not be done as a standard postprocessing, but decided on a case-by-case basis.
Have a look at the effect of the smoothing!
(py4vasp.plot(en,surface, "surface Ni d", xlabel="Energy (eV)", ylabel="DOS 1/eV")
+py4vasp.plot(en,ssurface, "surface Ni d"))
Next, plot the contributions to the CO $\pi$ states! The $1\pi$ state is dominated by oxygen, while the $2\pi^*$ state is dominated by carbon.
dos = mycalc.dos.read("C(px) C(py) O(px) O(py) 3(d)")
en = mycalc.dos.read()['energies']
o_pxy = dos["O_px"] + dos["O_py"]
c_pxy = dos["C_px"] + dos["C_py"]
bulk = dos["Ni_3_d"]
o_pxy = smooth(o_pxy)
c_pxy = smooth(c_pxy)
bulk = smooth(bulk)
(py4vasp.plot(en,o_pxy, "Opx+Opy", xlabel="Energy (eV)", ylabel="DOS 1/eV")
+py4vasp.plot(en,c_pxy, "Cpx+Cpy")
+py4vasp.plot(en,surface, "surface Ni d")
+py4vasp.plot(en,bulk, "bulk Ni d"))
By comparing the surface and bulk Ni $d$ states, we see that the surface states form additional weight around -8.5eV and -6.2eV, as well as unoccupied states around 1eV. This causes a depletion of the Ni $d$ surface states near the Fermi level. The CO states that are dominated by C above the Fermi level correspond to the CO-$2\pi^*$ state that causes back-donation. The CO states that are dominated by O below the Fermi surface are the CO-$1\sigma$ states. These form a fully occupied bonding and anti-bonding manifold which passify the surface by moving the Ni $d$ states away from the Fermi level.
Finally, consider how the occupation of the individual Ni $d$ orbitals changes due to the formation of a surface and the adsorbtion? Is there more weight in the $d_{3z^2-r^2}$ states or Ni $d_{xy}$? Which states hybridize with the CO states?
dos = mycalc.dos.read("1(dz2), 1(dxy), 3(dz2), 3(dxy), 5(dz2), 5(dxy)")
clean_dz2 = dos["Ni_1_dz2"]
clean_dxy = dos["Ni_1_dxy"]
bulk_dz2 = dos["Ni_3_dz2"]
bulk_dxy = dos["Ni_3_dxy"]
adsu_dz2 = dos["Ni_5_dz2"]
adsu_dxy = dos["Ni_5_dxy"]
(py4vasp.plot(en,smooth(clean_dz2), "clean_dz2", xlabel="Energy (eV)", ylabel="DOS 1/eV")
+py4vasp.plot(en,smooth(bulk_dz2), "bulk_dz2")
+py4vasp.plot(en,smooth(adsu_dz2), "adsu_dz2"))
The Ni $d_{3z^2-r^2}$ states hybridize with the CO molecule. The Ni $d_{xy}$ states lie within the plane of the slab and are barely influenced.
Compute and plot also the work function. What changes compared to the work function of pure Ni(111) in Example 6? In which way does the CO adsorbate modify the local potential?
import py4vasp
mycalc = py4vasp.Calculation.from_path( "./e07_CO_on_Ni111_pDOS")
workfunction_value = mycalc.workfunction.to_dict()['vacuum_potential'][1] - mycalc.workfunction.to_dict()['fermi_energy']
print(mycalc.workfunction)
print(f"The workfunction of CO on Ni (111) is {workfunction_value:.2f} eV")
mycalc.workfunction.plot()
Click to see the answer!
The work function changes only minimally since the Fermi energy changes by nearly the same amount as the vacuum level on the CO side of the slab.
We now see seven minima for the potential, corresponding to the five Ni positions as well as C and O. The dipole corrections ensure that we have two separate and well-defined vacuum potentials for our different surfaces. Setting DIPOL correctly to the center of mass, ensures that the step in the potential is located in the middle of the vacuum region. If we set LDIPOL = .FALSE., the potential will instead just smoothly and linearly connect the opposite sides of the cell, prohibiting a correct definition of the vacuum potential. A warning about a missing free-field region in the cell is printed out.
7.4 Questions¶
- Do bonding and anti-bonding orbitals form when a molecule adsorbs to a surface?
- Is it possible to plot molecular orbitals using LORBIT=11 and py4vasp?
By the end of this tutorial, you will be able to:
- calculate vibrational frequencies at the $\Gamma$ point of a molecule on a surface
8.1 Task¶
Perform a phonon calculation using the finite differences method and adjust the degrees of freedom for the molecule.
Lattice vibrations are the most important contribution of thermal excitations at surfaces. Thus vibrational spectra are critical in understanding thermodynamic properties of surfaces. Together with experimental techniques like high-resolution electron-energy-loss spectroscopy or Raman spectroscopy, surface phonon calculations provide detailed insights into the bonding situations of adsorbates. More information can be found e.g. in the review by Wöll, Appl. Phys. A, 1991, 53, 377–387.
More general examples of phonon calculations are given in a separate phonon tutorial. For an overview of theory, methods, and related tags, check out the phonons page on the VASP Wiki. Please read the how-to article on the finite differences method to understand the tags used in this example and learn about best practices.
8.2 Input¶
The input files to run this example are prepared at $TUTORIALS/surface/e08_CO_on_Ni111_Freqs. Check them out!
Consider how the POSCAR file uses selective dynamics to change the bond length between the C and O atom only along the z direction.
Ni C O
1.0000000000000000
2.4960869375885126 0.0000000000000000 0.0000000000000000
1.2480434687942563 2.1616746980061543 0.0000000000000000
0.0000000000000000 0.0000000000000000 19.0943572451033177
Ni C O
5 1 1
Selective dynamics
Direct
0.6666666666666643 0.6666666666666643 0.2108879566630485 F F F
0.0000000000000000 0.0000000000000000 0.3160636237749635 F F F
0.3333333333333357 0.3333333333333357 0.4227454972643674 F F F
0.6666666666666643 0.6666666666666643 0.5292842768315449 F F F
0.0000000000000000 0.0000000000000000 0.6373438393530619 F F F
0.0000000000000000 0.0000000000000000 0.7296118454904720 F F T
-0.0000000000000000 -0.0000000000000000 0.7900590120873849 F F T
SYSTEM = CO adsorption on Ni(111)
ICHARG = 1
# electronic minimization
ENCUT = 400
PREC = Acc
ISMEAR = 0
SIGMA = 0.05
ALGO = Fast
EDIFF = 1E-7
AMIX = 0.4
LMAXMIX = 6
# finite differences
IBRION = 5
POTIM = 0.015
NFREE = 2
# dipole correction
LDIPOL = .TRUE.
IDIPOL = 3
DIPOL = 0.5 0.5 0.452747064661102
K-Points
0
Monkhorst-Pack
11 11 1
0 0 0
Pseudopotential of Ni, C, O
The important INCAR tags are IBRION, NFREE, and POTIM. Add comments to the file once you have understood the meaning of these tags. EDIFF needs to be reduced compared to the default of $10^{-4}$ eV to ensure accurate forces. It is also recommended to use PREC = Accurate.
Given the settings of selective dynamics in the POSCAR file and NFREE in the INCAR file, how many different frequencies will we calculate and how many static calculations will VASP perform in total?
Click to see the answer!
We allow the two atoms in the molecule, C and O, to move along the z-direction only. They can either move in the same direction as each other, i.e., up and down with respect to the surface (CO-metal vibration), or move in opposite directions to each other, the so-called C-O stretch mode. We will calculate 2 frequencies. Since we set NFREE = 2, we have two displacements per atom per allowed direction, the Cartesian z-direction in this case. Including the calculation in the relaxed position, we end up with five calculations in total.
Why is POTIM set to such a small value?
Click to see the answer!
POTIM is the atomic displacement in Å. The calculation of the vibrational frequencies assumes that the harmonic approximation holds. If the displacement is too large the restoring force may not be directly proportional to the displacement.
8.3 Calculation¶
Open a terminal, navigate to this example's directory, and run VASP by entering the following:
cd $TUTORIALS/surface/e08_*
cp $TUTORIALS/surface/e05_*/CHGCAR .
mpirun -np 4 vasp_std
To view the results we look at the OUTCAR and search for the string Eigenvectors and eigenvalues of the dynamical matrix.
Click to reveal the frequency table.
Eigenvectors and eigenvalues of the dynamical matrix
----------------------------------------------------
1 f = 64.012322 THz 402.201282 2PiTHz 2135.221233 cm-1 264.733574 meV
X Y Z dx dy dz
2.496087 1.441116 4.026770 0 0 0
0.000000 0.000000 6.035032 0 0 0
1.248043 0.720558 8.072054 0 0 0
2.496087 1.441116 10.106343 0 0 0
0.000000 0.000000 12.169671 0 0 0
0.000000 0.000000 13.931469 0 0 -0.780749
0.000000 0.000000 15.085669 0 0 0.624845
2 f = 12.127787 THz 76.201133 2PiTHz 404.539432 cm-1 50.156475 meV
X Y Z dx dy dz
2.496087 1.441116 4.026770 0 0 0
0.000000 0.000000 6.035032 0 0 0
1.248043 0.720558 8.072054 0 0 0
2.496087 1.441116 10.106343 0 0 0
0.000000 0.000000 12.169671 0 0 0
0.000000 0.000000 13.931469 0 0 -0.624845
0.000000 0.000000 15.085669 0 0 -0.780749
Now modify the POSCAR so that all degrees of freedom are active for both the C and the O atoms by changing selective dynamics from F F T to T T T. Restart the calculation and figure out what additional modes were added.
Click to see the answer!
We get four additional modes. In two of those the C and O atoms move in opposite directions parallel to the surface, effectively tilting the CO molecule back and forth in the Cartesian x and y direction. The frequency of these modes is slightly higher than the collective movement normal to the surface we have seen before. The other two modes describe a collective movement in x and y directions respectively with less than 10% of the frequency.
8.4 Questions¶
- Is selective dynamics considered for IBRION=5?
- What issue may occur if the step size in the finite difference method for phonon calculations is set too large?