Introduction

We will investigate some properties of crystalline silver: the relative surface energies for the 100 and 111 cleaved faces.  In this post, we will use DFT calculations using CASTEP [1] software to calculate the energy of bulk silver as well as relaxing a cleaved surface. These properties can reveal the preferred behavior in surface Ag formation.  Note that silver is an FCC crystal with a lattice constant of 4.09 Å [2] and that this experimental result will be used in the calculations.

Surface Energy Calculation

First, in order to calculate the surface energy for a particular surface, we need to refer to the following equation found in “Density Functional Theory : A Practical Introduction”:

So for both 100 and 111 surfaces, we need to calculate the energy of bulk Ag.

Using CASTEP with the GGA PBE functional and OTFG ultrasoft pseudopotentials, we run through different energy cutoffs to find an appropriate energy cutoff for energy convergence.  Ag has the electron configuration of 1s2 2s2 2p6 3s2 3p6 3d10 4s2 4p6 4d10 5s1, and the pseudopotential treats 4s2 4p6 4d10 5s1 as the valence electrons.  The convergence tolerances were chosen somewhat arbitrarily to be small:  energy at 2.0e-5 eV/atom, force at 0.05 eV/Å, stress at 0.1 GPa, and displacement at 0.002 Å.  For all the rest of the options, the defaults were used.

Energy convergence with respect to ENCUT

Change in energy from the previous ENCUT energy calculation

Seeing this, it is sufficient for our purposes to use a cutoff of 900 eV.

Next, we check for k-point energy convergence at our cutoff energy of 900 eV.

k-point convergence with energy as a function of k-points

Here we see that 88 k-points is sufficient.

From this we find that -4006.239 eV/atom.

Next, we need to relax the surface of 111 Ag.  We choose a 7 atom thick slab, and a 10 Å thick vacuum.  The ENCUT was set to be 900 eV once again and an nxnX1 grid of k-points was selected to keep the k-point density equivalent in the plane of the surface, where n is varied as energy is converged with respect to k-points.   It should be noted that fewer k-points are needed in the direction of the vacuum which is why only 1 k-point is needed in the direction of the vacuum.  Both sides were relaxed while the central layers were fixed to simulate bulk.  3 atomic layers were fixed and the “top” 2 on each side were allowed to relax.

The exact same procedure is repeated for the 100 surface with a 900 eV basis cutoff energy, and an nxnx1 grid of k-points.

The (111) unit cell. Here the yellow highlighted atoms are held fixed while the outer layers are allowed to relax.

Energy convergence on the (111) surface with respect to k-points.

The change in energy from the previous k-point Energy calculation.

The (100) unit cell. Here the yellow highlighted atoms are held fixed while the outer layers are allowed to relax.

Energy convergence on the (100) surface with respect to k-points.

The change in energy from the previous k-point energy calculation.

The following Energies were found:  -28043.265 and -28043.196.  This gave the following results for the surface energy densities:

0.0244 eV/Å^2

0.0286 eV/Å^2

Conclusion

The results are in agreement with the known result [4] that says that the (111) surface is energetically favorable to the (100) surface of Ag. The energy of the (111) surface is found to be 15% lower than the (100) surface.  In previous comparisons against multiple functionals [4], the percentage difference was lower.  This result shows that the (111) surface is the preferred face for crystal growth.

With better computational ability, these results should be checked against varying layers for the surface calculations as well as a higher ENCUT for both the bulk calculation and the surface calculations.

References

1. “First principles methods using CASTEP” Zeitschrift fuer Kristallographie 220(5-6) pp. 567-570 (2005) S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson, M. C. Payne
2. https://periodictable.com/Elements/047/data.pr.html
3. Density Functional Theory: A Practical Introduction. (2009)  David S. Sholl, Janice A. Steckel
4. Patra, A., Bates, J. E., Sun, J., & Perdew, J. P. (2017). Properties of real metallic surfaces: Effects of density functional semilocality and van der Waals nonlocality. Proceedings of the National Academy of Sciences of the United States of America, 114(44), E9188–E9196. https://doi.org/10.1073/pnas.1713320114

Lingjie Zhou

Introduction

Silicon(Si) is a widely used substrate for molecule beam epitaxy growth due to its high quality and low cost. However, the real surface of Si doesn’t resemble the plane directly cut from bulk material. Due to the recombination of dangling bonds, usually there will be reconstruction patterns formed at the surface. Such phenomenon have been broadly investigated and confirmed through Scanning Tunneling Microscope. To explain and predict such phenomena, we used plane-wave basis sets with ultrasoft pseudopotentials to perform DFT methods to calculate the energy of the surface with and without reconstruction.

Method

The CASTEP [1] package is used to carry out the DFT calculation. The exchange and correlation function were calculated using the Perdew, Burke and Ernzerhof(PBE) functional described within the generalized gradient approximation(GGA) [3]. The ‘on the fly’ generated ultrasoft pseudopotential for Si has a core radius of 1.8 Bohr(0.95 Angstroms) and was generated with 4 electrons in the valence panel with (3s2 3p2) as the electronic configuration.

First, the optimization of bulk material is done to make sure there is no artificial stress in the model [2]. The calculation starts from the experimentally reported result with a=5.381 A. 700eV cutoff energy and 7×7×7 kpoint set is used and the optimized value is a=5.468 Å. The number of kpoint and cutoff is initially hypothesized as sufficient before checking convergence and the convergence will be further checked.

Energy cutoff convergence

The convergence of energy cutoff is first checked by carrying out the optimization of structures with a=5.381 Å and 7×7×7 k points but different cutoff energy. The optimized result (a=5.481) is well converged if we use cutoff energy higher than 600 eV. For the rest of the calculation we would use 700 eV as cutoff energy.

Fig 1 Convergence of energy cutoff

Kpoint convergence

Here is the result if the optimization all start from experimental result(a=5.381 Å) and 700eV cutoff but use the different number of kpoints(6×6×6, 8×8×8, 10×10×10, 15×15×15). The result is well-converged so that we will pick 6×6×6 for our later kpoint mesh.

Fig 2 Convergence of number of k points

Si(001) surface

First, we calculate the surface of Si(001) without reconstruction. Slabs of 3 layers, 4 layers and 5 layers of Si atoms and 10 Å vacuum is used for the calculation. During the calculation, only the atoms of the top layer are allowed to move while the rest are fixed. The optimized structure is very similar to the bulk. To calculate the reconstruction pattern, top two atoms are shifted towards each other. The optimized structure is similar to the experimentally verified reconstruction pattern.

Fig 3 Surface without reconstruction

Fig 4 Surface with reconstruction

The layer dependence of surface energy is plotted below. Due to the limitation of time and computational resources, it’s hard to reach a well-converged result. But the surface energy of reconstructed structure is always smaller than the one without reconstruction. Thus, the reconstructed pattern is preferred. With limited time and computation resources, only one reconstruction pattern is checked and other reconstruction patterns remain to be checked.

Fig 5 layer dependence of surface energy

Reference

Introduction

Adsorption of small molecules on transition metal surfaces is of great interest among chemists, since adsorption behavior could be essential to the catalytic process happening on transition metal surfaces. For adsorption on (111) surface of fcc metals, there are 4 different kinds of sites, as shown in figure 1. It is worth noticing that there are two different possible sites for threefold adsorption situations, which are called fcc and hcp site respectively. Since the bulk structure of Rh has been studied in previous work, in this work, the preference of adsorption of CO of low surface coverage on Rh (111) is assumed to be on-top site according former computational study[1]. Intuitively, CO adsorption on-top sites at low coverage is similar to the situation that CO bonds to central metal atom as a ligand, which should have a lower energy comparing to other adsorption sites.

Figure 1 (a) fcc site, (b) on-top site, (c) hcp site, (d) bridge site

Method

In this work, plane-wave based DFT code CASTEP is used. The exchange-correlation functional is described by a GGA-PBE functional, and pseudopotentials are kept default setting. As for pseudopotential, the “on the fly” generated ultrasoft pseudopotential for Rh is used with a core radius of 1.6 Bohr(0.847Å) and was generated with a 4d8 5s1 valence electronic configuration . The SCF tolerance is set to be 2e-6 eV per atom for all the optimization and single point energy calculations.

Optimized bulk structure of Rh is taken from previous work[2]. According to previous computational work[3], a 2×2 Rh(111) surface slab model of 4 layers of atoms is used to simulate the behavior of Rh surface. As for the optimization of surface slab model, energy cutoff is kept consistent with that for bulk calculation, while k points gird setting is set to 7x7x1. In order to get comparable energies for molecule and surfaces, the same setting is used for optimization of CO in vacuum. CO adsorption is modeled with single CO adsorbed on different sites of optimized surface model. The corresponding adsorption energy is calculated with the equation[4]

\begin{equation} E_{ads}=E_{Rh-CO}-E_{surface}-E_{CO} \end{equation}

Results

1. k points grid optimization

With an energy cutoff of 625 eV, single point energies of 1×1 slab models are calculated. In figure 2, the energies relative to 9x9x1 k points gird is shown. As a trade off of accuracy and efficiency, 7x7x1 k points grid is used to get valid results for surface calculations.

Figure 2 k points convergence test, energy relative to maximum k points grid is convergent for grid larger than 7x7x1 gird.

2. Surface optimization

Although 2×2 slab model will be used for adsorption of CO calculations in order to make sure a 1/4 ML coverage, here surface structure optimization is performed for 1×1 slab model. Due to the limitation of computation resources, 2×2 slab model optimization won’t converge in reasonable time, a 1×1 slab model is optimized instead. Since Rh surface slab model is periodic system, a 1×1 slab won’t introduce significant error to final results. Meanwhile, the energy of the supercell should be scalable, which makes it reasonable to get energy of 2×2 slab out of 1×1 slab. Only top two layers are relaxed while two layers at bottom are kept fixed to simulate bulk effect on surface atoms.

BFGS algorithm is used to perform the optimization with energy convergence of 2e-5 ev/atom and force convergence of 0.05 eV/A. The optimized structure has an energy of -12123.0410 eV.

3.Optimization of CO in vacuum

In order to obtain consistent value for molecular energy, CO configuration is optimized in the same lattice as the surface slab model. The optimized bond length is 1.138 A, with an energy of -596.123 eV.

4. Optimization of CO adsorption

A 2×2 slab model for Rh(111) surface is used to get a optimized structure of one single CO adsorption on- top site, which corresponds to a 1/4 ML coverage. Such systems would significantly increase the cost of computation if full optimization is performed. Thus, only the distance between CO molecule and Rh atom is allowed to relax to find the most stable configuration of CO adsorption on-top site. Such modelling of CO adsorption is based on the assumption that adsorption of CO with low coverage won’t induce significant surface structure reconstruction.

The optimization algorithm and convergence setting is consistent with those for surface structure optimization. The optimized distance between C and Rh atom is 1.843 A with a total energy of -49089.5603 eV. The CO adsorption energy for Rh on-top site is -1.273 eV. Meanwhile, C-O bond is stretched to 1.160 A, which suggests the weakening of bond strength between C and O.

Conclusion

The adsorption energy of CO on on-top site of Rh (111) surface is -1.273 eV with an optimized bond distance of 1.843 A, which is consistent with former computational results at 1/4 ML coverage, 1.83A [1] and 1.842 A[3], respectively. Since different levels of theory or different exchange-correlation functionals have been used, the bind energies are not comparable, it is reasonable to compare the bonding length between metal atoms and adsorbate to see if the optimized geometry is consistent with former study. Further work may focus on the influence of surface coverage on the adsorption energy and bond weakening of C-O bonds. On the other hand in order to get adsorption values comparable to experimental ones, a test for different exchange-correlation functionals and thermodynamic parameters are required.

Reference

[1] Mavrikakis, M., Rempel, J., Greeley, J., Hansen, L. B. & Nørskov, J. K. Atomic and molecular adsorption on Rh(111). J. Chem. Phys. 117, 6737–6744 (2002).

[2]Prediction of preferred structure and lattice parameters of Rh https://sites.psu.edu/dftap/2020/02/03/prediction-of-preferred-structure-and-lattice-parameters-of-rh/

[3] Krenn, G., Bako, I. & Schennach, R. CO adsorption and CO and O coadsorption on Rh(111) studied by reflection absorption infrared spectroscopy and density functional theory. J. Chem. Phys. 124, 144703 (2006).

[4] Sholl, D. & Steckel, J. A. Density functional theory: a practical introduction. (John Wiley & Sons, 2011).

Da Zhou

Department of Physics, Penn State University

Introduction

Transition metal dichalcogenides (TMD) are materials with the composition of MX2, where M is a transition metal like tungsten, molybdenum, etc, and X is a chalcogen atom such as sulfur, selenium, or tellurium. TMDs have layered crystal structures, where the interlayer bonding is formed by the Van der Waals forces. Therefore, monolayer TMDs can be obtained from either mechanical exfoliation or chemical vapor deposition (CVD). A recent experimental breakthrough is doping vanadium into WS2 and WSe2 and bring magnetism into the picture. Experimentally, the doping concentration of vanadium can be up to 10%. Here, I build a toy model, where one vanadium atom is doped into a 3 by 3 WSe2 supercell (and gives a doping concentration of 11%), and performed energy calculations using DFT to solve the Hessian matrix and find the vibration frequencies.

Methods

For the benefit of time, the calculation quality was rather coarse. The DMol3 package and an LDA-PWC functional were used for the calculations. The spin was unrestricted. SCF energy tolerance was 0.0001Ha. Dipole slab correction was applied though whether applying this correction or not would not change the binding energy. And a 3 by 3 by 1 kpoints grid was used. The global cutoff was 3.7Å. The DFT semi-core pseudopotential was used for the core treatment. And a DN 3.5 basis set was used.

Results and Discussions

A bulk WSe2 structure was first built and then cleaved along the 001 direction to get the monolayer WSe2. A vacuum layer with a thickness of 20Å is added. A 3 by 3 supercell is built and the central tungsten atom is replaced by a vanadium atom.

Then, a geometry optimization was performed to relax the supercell. After a fine optimization which took about 2 hours, the system was still not at its energy minimum. This is found out when the vanadium atom is displaced for 0.03Å to solve the Hessian matrix. This reflects certain limitations to the effectiveness of the geometry optimization. To better relax the supercell and find the energy minimum, the vanadium atom was first displaced for 0.03Å and then the geometry optimization was performed. This time, the energy was minimized and a reasonable Hessian matrix was calculated.

A top view of the supercell. The grey atom is the doped vanadium atom.

To calculate the Hessian matrix, the vanadium atom needs to be displaced in x, y, and z directions. The vibration frequencies will increase when the displacements are increased, so a roughly 1% displacement (0.03Å, since the a/b lattice constant for WSe2 is 3.327Å) is used.

The 3 by 3 Hessian matrix is listed below in units of Ha/(Å^2):

0.614333, 0.023111, -0.01

0.023111, 0.668555, -0.006222

-0.01 -0.006222, 0.794555.

The eigenvalues of the Hessian matrix were solved by the Wolfram Mathematica: 0.79556, 0.676342, 0.605541. The mass of one vanadium atom is 8.4590343 × 10^-23 g. Therefore, the vibrational frequencies can be converted in units of cm^-1, and are 339.95cm^-1, 313.44cm^-1, 296.58cm^-1. For comparison, an independent calculation which directly uses the frequency analysis from the geometry optimization was also performed. All the other atoms but vanadium were fixed since only the vanadium atom was displaced during the calculation of the Hessian matrix. The vibration frequencies calculated from this method are 372.96cm^-1, 342.12cm^-1, and 330.61cm^-1. The 2 sets of results are quite different, but the ones calculated from the Hessian matrix depends on how far the vanadium atom is displaced.

Conclusion and Potential Improvements

The vibration frequencies of the vanadium atom doped in WSe2 were numerically solved in a 3 by 3 supercell and the results were compared with results calculated from a different method.  For future improvements, a better model or a larger supercell can be used to simulate different doping levels. One can also change the magnitude of the displacement and repeat the calculations to solve the new Hessian matrix and compare the frequencies calculated from different displacements if time permits.

Author: Andrew Wong

1. Introduction

The goal of this post is to determine the preferred adsorption site of atomic oxygen (O) on the iron (Fe) (110) surface using DFT techniques. Due to their catalytic and magnetoelectric properties, iron and iron oxides are vital in various applications [1]. On the BCC Fe (110) surface, the most stable termination of Fe [2], there are four potential adsorption sites that O can adsorb onto; on top (OT), long bridge (LB), short bridge (SB), and the pseudo-three fold bridge (TB). To better understand adsorption of atomic oxygen on the Fe (110) surface, plane wave basis Density Functional Theory (DFT) was implemented with the Vienna Ab initio Simulation Package (VASP) [3] to calculate the corresponding adsorption energy of each of the four sites. By comparing these energies, the most preferred adsorption site of atomic oxygen on the Fe (110) surface can be determined by which site had the most negative adsorption energy.

2. Methodology

2.1 Calculation Parameters

The DFT calculations performed in this post utilized the plane wave basis set with pseudo potentials method in Vienna Ab initio Simulation Package (VASP). The Perdew, Burke, and Ernzerhof functional (PBE) was implemented to model the electron-electron exchange and correlation energies. To represent the ion-core electron interactions, the projector augmented-wave (PAW) method was implemented [4]. The Monkhorst-Pack was utilized to model the K-point grid. An Fe core radius of 2.4 Bohrs (1.27 Å) with a panel of 8 valence electrons ( 3d6 and 4s2) were implemented with an electornic convergence tolerance of 2E-06 eV in all calculations. For all calculation, a cutoff energy of 400 eV was implemented as this was seen to be optimal energy cutoff value during the geometry optimization calculations.

2.2 Optimization of Triplet Oxygen in Gas Phase

Before adsorption energies of atomic oxygen on Fe (111) surface can be calculated, the energy of diatomic oxygen must be calculated first. A geometry optimization of diatomic oxygen in the gas phase was performed in a 15 x 15 x 15 Å vacuum cube at an energy cutoff value of 400 eV and a 1 x 1 x 1 K Point set.

2.3 Construction and Optimization of Fe(110) Vacuum Slab

Based on literature that also utilized VASP for electronic structure calculations ,the 3 x 3 Fe (110) slab model was determined to be the most stable termination of Iron [2]. Specifically, a vacuum region of approximately 10 Å with a k-point mesh of 5 x 5 x 1 was utilized to constructed the slab model. A cutoff energy of 400 eV was also implemented as well. Dipole corrections and selective dynamics were implemented as the bottom two layers were fixed to imitate the bulk iron and the two top layers were relaxed. Figure 1 below presents the Fe (110) surface slab model constructed.

Figure 1a: Side View of 3 x 3 Fe (110) Slab Model

Figure 1b: Top View of 3 x 3 Fe (110) Slab Model

2.4 Adsorption of Atomic Oxygen on Fe(110) Surface

For the BCC Fe (110) surface, there are four potential sites of adsorption: On top (OT), Short bridge (SB), Long bridge (LB), and Pseudo-Three Fold Bridge (TB) [5]. A top view figure of these four sites are presented in figure 2 below.

Figure 2: Top View of Potential Adsorption Sites of BCC Fe (110)

In order to determine the most preferred site of atomic O on Fe (110) surface, adsorption energies of O onto the Fe (110) surface must be calculated from geometry optimization calculations. The adsorption site corresponding with the lowest adsorption energy is the most thermodynamically favorable site for O on Fe (110).  In order to calculate the adsorption energies for each site, the following equation was used. [6]

Three terms are needed to calculate the adsorption energy. The energy of the Fe surface with O adsorb onto it will need to be determined first. To calculate the energy of a O adsorbate, half of the energy of diatomic oxygen will be utilized as it is easier to perform a geometry optimization on diatomic oxygen rather than atomic oxygen and there are more readily available experimental values for diatomic oxygen. Lastly, the energy of the bare Fe (110) surface without the adsorbate will be needed. As a result, the adsorption energies can then be compared to determine which adsorption site is the preferred site O to adsorb on.

3. Results and Discussion

Using the optimal K-point grid of 5 x 5 x 1 and cutoff energy of 400 ev, Figure 3 was constructed below to represent the top view of the four adsorption sites of Fe (110) with 1/9 ML surface coverage of O after geometric optimization done in VASP.

Figure 3. Optimization of atomic O on Fe (110) top view for the following four sites a. On Top b. Long Bridge c. Short Bridge d. Pseudo-Three Fold Bridge

The table below lists the results of the adsorption energy and distance of the O atom from the Fe (110) surface for the four potential adsorption sites. By comparing the adsorption energies between each site, it can be seen that the adsorption O favors the long bridge site as this site significantly had the most negative adsorption energy compared to the other three sites. It is important to note that, although not as favorable by 0.42 eV, the short bridge site is also a  thermodynamically favorable site due to its negative adsorption energy value. This may be due to the multiple bond interactions between multiple oxygen atoms seen in the short and long bridge sites. Although the Pseudo-Three Fold bridge has multiple oxygen atom interactions, it required multiple attempts to converge exactly to this position, as it requires it to be between three oxygen atoms. As a result, the TB site is thermodynamically unfavorable for the adsorption of O onto the Fe (110) surface due to it having the highest adsorption energies. The OT Adsorption site was also favorable but it tend to converge towards the long bridge and short bridge adsorption site and did not have nearly as strong of an adsorption energy by approximately 1 – 1.5 eV.

4. Conclusion

The thermodynamically favorable adsorption sites for the 1/9 ML surface coverage of atomic O on the Fe (110) were the long bridge, short bridge sites, and the on top sites, with the long bridge adsorption site being the most preferred adsorption site. Additionally, the pseudo-three fold bridge sites was not a thermodynamically favorable site due to their positive adsorption energies. As a result, they converged to the local minimums of the long bridge and short bridge sites during the geometric optimization of the adsorbed Fe (110) surface, rarely converging towards the on top configuration. Although the DFT calculations creates a notion that the preferred adsorption site of O on the Fe (110) is the long bridge site, additional considerations, such as surface coverage variations, non-vacuum conditions, and impurities within the bulk surface model, should be noted when studying this particular adsorption chemistry in a more practical and experimental setting. An experimental study examining the adsorption of atomic oxygen on the Fe (110) surface with the use of standard surface characterization methods, such as Auger electron spectroscopy and conversion electron Mossbauer spectroscopy, has confirming agreement that the long bridge site is the preferred adsorption site [5].

5. Citations

[1] Ossowski, Tomasz, and Adam Kiejna. “Oxygen Adsorption on Fe(110) Surface Revisited.” Surface Science, North-Holland, 13 Mar. 2015, www.sciencedirect.com/science/article/pii/S0039602815000618.

[2] Maheshwari, Sharad, et al. “Elementary Kinetics of Nitrogen Electroreduction on Fe Surfaces.” AIP Publishing, AIP Publishing LLC, 28 Jan. 2019, aip.scitation.org/doi/full/10.1063/1.5048036.

3] Init.at. “Vienna Ab Initio Simulation Package.” VASP, www.vasp.at/.

[4] Rostgaard, and Carsten. “The Projector Augmented-Wave Method.” ArXiv.org, 12 Oct. 2009, arxiv.org/abs/0910.1921.

[5] Freindl, K., et al. “Oxygen on an Fe Monolayer on W(110): From Chemisorption to Oxidation.” Surface Science, North-Holland, 13 July 2013, www.sciencedirect.com/science/article/pii/S0039602813002021.

[6] Nguyen, Angela, and Angela Nguyen. Density Functional Theory and Practice Course, 2 Apr. 2019, sites.psu.edu/dftap/2019/04/02/exploring-the-effects-of-surface-coverage-on-the-binding-site-and-adsorption-energy-of-atomic-o-on-pt111/.

Introduction

The interaction of hydrogen with transition-metal surfaces is of great importance in heterogeneous catalysis, metallurgy, energy storage, and fuel cell technology [1]. One can find differing studies in the literature regarding the preferred adsorption site of atomic hydrogen. Periodic DFT-GGA calculations by Eder et al. [1] suggested that H-atom favors the 2-fold bridge site on Fe(100), and the 4-fold hollow site is slightly higher in energy. Whereas, the study by Jiang and Carter [2] concludes that the fourfold-hollow site is the most stable while both the 2-fold bridge and the 4-fold hollow sites are true minima. In this post, spin-polarized density functional theory calculations are performed to characterize the atomic hydrogen (H) adsorption on the Fe (100) surface using the Vienna Ab initio Simulation Package (VASP).

Methodology 

We use the plane wave DFT calculations to evaluate the adsorption energy of atomic hydrogen on the Fe (100) surface on the bridge, hollow and on the top site as shown in Fig.1. All the calculations are performed using the Vienna Ab initio Simulation Package (VASP). The exchange and correlation energies were calculated using the PBE functional described within the projector augmented wave method (PAW). Electronic convergence tolerance of 1E-05 was used for all the calculations. The core electrons were treated using pseudo-potential with a core radius of 2.3 Bohrs (1.22 Angstrom) generated with a panel of 8 valence electrons ( 3d6 4s2). All the calculations are performed with the spin-polarized calculations. 

Fig. 1 H adsorption sites shown on the Fe (1 0 0) 3×3 cell

Cut-off Energy and K-Point Optimization 

The essential step before performing a plane-wave basis set calculation is to optimize the k-points and cut-off energy. This process is performed as follows:

K-point  

The numerically determined optimal lattice parameter for Fe bcc crystal structure is 2.835 Angstrom which agrees reasonably with the experimental value of 2.856 Angstrom [3]. The K-point optimization is performed at this lattice parameter with a 2×2 Fe (100) surface with 5 atomic layers containing 20 Fe atoms in the slab.

Table 1: K-point grid vs the number of K-points used

Fig.2  shows the absolute change in the consecutive total energy (|ΔE|) with respect to the number of irreducible k-points. The corresponding K-point grid is shown in Table 1. As a convergence criterion of |ΔE| less than 0.002 eV is used for the K-point convergence. This corresponds to the optimal grid of 7x7x1. 

Cut-off energy 

In Fig. 3, the absolute change in the consecutive total energy (|ΔE|) is plotted as a function of the cut-off energy. As a convergence criterion, it is assumed that the total energy is converged if |ΔE| is less than 0.005 eV. From Figs. 3 and 4, it is clear that the cut-off energy of 600 eV sufficient considering both the functionals. A lattice parameter of 2.835 Angstrom and K-point grid of 7x7x1 is used for this convergence. 

H adsorption energy 

As mentioned earlier a 5 atomic layer 2×2 Fe (100) slab with a 10 Angstrom vacuum layer containing  20 Fe-atoms is constructed (see Fig.1). Later, the H-atom is placed on several probable surface sites and a geometry optimization is performed for 0.25 ML coverage. The adsorption energy of the H-atoms is calculated as per the following expression

E_{adsorption, H-atom} = E_Fe-H – E_Fe – 0.5E_H2

Based on the K-point and cut off energy convergence study, cut-off energy of 600 eV and a K-point grid of 7x7x1 is employed in these calculations.  

Results

Table 2 The adsorption energies of H-atom on a Fe (100) surface in eV.

Table. 2 shows the adsorption energy for the different sites on a Fe (100) surface. It also contains data from the references cited earlier. The current calculations match quite well with the literature data. The current calculations predict that the hollow site is energetically favored for H-atom adsorption. Additionally, calculations with more number of layers of Fe are also performed to assess its effect on the adsorption energy. Table. 3 lists the adsorption energies for different Fe layers for all the sites. The typical differences converge to within 0.01 eV or 2-5% of the adsorption energy value for 9 layers. 

Table 3 The adsorption energies of H-atom on a Fe (100) surface with different atomic layers.

Conclusions

This post investigates the preferred adsorption site of H-atom on the Fe (100) surface using DFT calculations. The calculations show consistent results with the previous theoretical studies. The results of current calculations predict that the hollow site to be slightly favorable compared to the bridge site for H-atom adsorption. The calculations also demonstrate the effect of using a bigger metal slab for the calculations. The adsorption energy change becomes less than 0.01 eV from 9 layers onwards.

References

1. Eder, M.; Terakura, K.; Hafner, J. Phys. ReV. B 2002, 64, 115426
2. Jiang, D. E., and Emily A. Carter. “Diffusion of interstitial hydrogen into and through bcc Fe from first principles.” Physical Review B 70.6 (2004): 064102.
3. Greenwood, Norman Neill, and Alan Earnshaw. Chemistry of the Elements. Elsevier, 2012.

Author: Jeremy Hu

Introduction

Titanium dioxide (TiO2) is a frequently-used metal oxide for depositing metal catalysts onto, such as those for hydrodeoxygenation (HDO) and hydrogenation reactions [1].  Since catalytically-active metals are generally expensive, depositing them onto a metal oxide support such as TiO2 helps to stabilize and maximize the usage of the active metal sites. Although several polymorphs and surfaces of TiO2 exist, the anatase form has been shown to be favorable under reduction conditions, while the (001) surface has demonstrated increased catalytic activity [1] [2]. Additionally, silver (Ag) single atom catalysts (SACs) on TiO2 anatase (001), which have highly dispersed Ag in the single-atom limit on the surface of TiO2, have been of interest due to their increased selectivity and activity in HDO reactions [1].  In this study, the surface energy of TiO2 anatase (001) with 1, 2, and 3 layered slabs was calculated. Additionally, the preferred binding sites of Ag single atoms on TiO2 anatase (001) were identified.

Methods

Electronic Methods

Density functional theory (DFT) analysis of TiO2 anatase (001) and Ag on TiO2 anatase (001) was calculated using the plane-wave basis set in the Vienna Ab Initio Simulation Package (VASP) [3]. The Perdew–Burke-Ernzerhof (PBE) exchange correlation functional was used as well as the Generalized Gradient Approximation (GGA) for the exchange correlation functional type [4]. PBE+D3 was used for dispersion corrections, along with the projector augmented-wave method (PAW) to correct for core-valence interactions [5] [6]. The Ti d-electrons were accounted for using the DFT+U approach, which corrects for the self-interaction error and underestimation of the band gap [5]. The Hubbard’s parameter (U) was set for Ti to be U = 4 eV. The self-consistent field tolerance of the calculations used a convergence criteria of < 0.05 eV/Å.

The Monkhorst-Pack k-point mesh for the reciprocal space was 3 x 3 x 1 and the cutoff energy was 450 eV, both of which are above the minimum for convergence of similar TiO2 DFT calculations in the literature [7] [8]. The valence electrons considered for each atom type were O (2s² 2p⁴), Ti (3s² 3p⁶ 4s² 3d²), and Ag (5s¹ 4d¹⁰). To minimize dipole interactions between periodic slabs in the z-direction (i.e., normal to the surface), a vacuum space of 10 Å was used between slabs. For the surface energy calculations, a 1 x 1 slab was used. One “layer” was considered to be the smallest thickness of TiO2 with the stoichiometric number of atoms (i.e., Ti_xO_2x). The adsorption energy of Ag was calculated using a 2 x 2 supercell of TiO2 anatase (001).

Calculation Methods

The surface energy for the different number of layers was calculated as follows, where A = the surface area of the slab (i.e., a² where a is the lattice parameter) and n is the number of TiO2 units in the slab [9]:

$$E_{surf} = \frac{1}{2A} (E_{slab} – n E_{bulk})$$ (Eqn. 1)

The 2A in the denominator accounts for the fact that there are two surface areas of interest from the slab (i.e., the top and bottom of the cell). All the atoms in the slab were allowed to fully relax during the geometric optimization. The bulk energy of TiO2 anatase was taken as -25.539 eV/unit of TiO2 [11]. The surface area of each slab (A) was 14.26 Ų.

For the calculation of adsorption energy of Ag on various sites on TiO2 anatase (001), the adsorption energy was calculated as follows, with the isolated energy of the Ag atom taken under vacuum:

$$E_{ads} = E_{Ag,TiO2} – E_{TiO2} – E_{Ag}$$ (Eqn. 2)

Results and Discussion

1. Slab models and surface energy calculations

The geometry optimization of a 1 x 1 single-layered slab of TiO2 anatase (001) was performed, with one layer defined as the smallest thickness for stoichiometric TiO2 (Fig. 1).

 Figure 1. Single-layered slab of TiO2 anatase (001).

Next, a geometry optimization (i.e., energy minimization) of a 1 x 1 two-layered cell of TiO2 anatase (001) was performed (Fig. 2).

Figure 2. Two-layered slab of TiO2 anatase (001).

Finally, the surface energy of a 1 x 1 three-layered slab was calculated (Fig. 3).

Figure 3. Three-layered slab of TiO2 anatase (001).

The surface energy for various layered slabs is displayed in Table 1.

Table 1. Surface energy (eV/Ų) with respect to number of layers in a 1 x 1 slab of TiO2 anatase (001)

It would appear that the single, two, and three-layered slab models of TiO2 have reasonably consistent surface energies. Additionally, the three calculations relatively approximate the literature values of the surface energy of TiO2 anatase (001) of 0.0562 eV/Ų [10]. The similarity of surface energy values suggest that slab models beyond the single-layered slab are sufficient for the calculation of the surface energy of TiO2 anatase (001).

2. Ag adsorption sites

The following sites of Ag single atoms on TiO2 anatase (001) were identified as possible energy minima (Fig. 4).

Figure 4. Four possible adsorption sites of Ag single atoms on TiO2 anatase (001) (top view).

The O2C and O3C refer to the coordination number of the oxygen atom of interest. The following adsorption energies of Ag were calculated using Eqn. 2 (Table 2).

Table 2. Adsorption energy (eV) of Ag single atoms on various sites on the TiO2 anatase (001) surface

Regardless of the binding site, the binding of Ag single atoms on the surface of TiO2 anatase (001) is energetically favorable, since all the adsorption energies are negative. Furthermore, the DFT calculations suggest that position 1 is the most energetically favorable position for Ag single atoms to bind, since its adsorption energy is the lowest. Thus, we propose that deposition of Ag on TiO2 anatase (001) in the single-atom limit likely results in a majority of positions where the Ag atom is bridged between the two second-coordinated oxygen atoms on the surface (Fig 5).

Figure 5. Most energetically favorable position of Ag SAC on TiO2 anatase (001) (side view).

Conclusion

The surface energy calculations of single, two, and three-layered slabs of TiO2 anatase (001) suggest that models beyond a single-layered slab are sufficient to approximate the surface energy. The adsorption of Ag single atoms on a 2 x 2 supercell of TiO2 anatase (001) demonstrated possible energy minima of adsorption sites. The adsorption study suggests that the deposition of Ag in the single-atom limit results in a state where Ag is bridged between two O2C oxygen atoms. These results are promising for understanding catalyst behavior and structure of Ag single atoms on TiO2 anatase (001) supports.

Citations

[1] Liu et al., Journal of Catalysis 2019, 369, 396-404

[2] Wan et al., J. Phys. Chem. C 2018, 122, 17895-17916

[3] Kresse, G.; Furthmueller, J., Physical Review B: Condensed Matter (1996), 54 (16), 11169-11186

[4] J. P. Perdew, K. Burke, M. Enzerhof, Phys, Rev. Lett., 1996, 77, 3865.

[5] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys., 132 (2010), Article 154104

[6] P.E. Blöchl, Phys. Rev. B, 50 (1994), pp. 17953-17979.

[7] H. J. Monkhorst, J.D. Pack, Physical Review B, 1976, 13, 5188-5192.

[8] H.V. Thang et al., Journal of Catalysis, 367 (2018), 104-114

[9] Sholl, D.S. Steckel, J.A. (2009) “Density Functional Theory: A Practical Introduction.” Wiley, 96-98.

[10] Lazzeri, M.; Vittadini, A.; Selloni, Phys. Rev. B 2001, 63, 155409

[11] A. Kiejna, T. Pabisiak, S.W. Gao, J. Phys.: Condens. Matter 18 (2006) 4207-4217

by Yingdong Guan

Introduction

Surface energy is an important property for crystals. In most cases, the lowest surface energy determines the most stable surface and will thus influence the crystal growth process. [1] Also, there are many interesting phenomena happening on the surface of thin films and on the interfaces of multiple different thin films. Some of those phenomena are related to surface energy. In this post, we used the Density Functional Theory (DFT) calculation to calculate surface energy for different surfaces of Au crystal, namely (100) and (111). Slab models [2] are used to simulate different surfaces of Au crystal. The surface energy of different surfaces are calculated based on convergence tests and are compared to determine which surface is energetically favored.

Bulk Au

Au lattice with the cubic close-packed structure was constructed. The crystal structure of Au is shown in Figure 1.

 Figure 1 Crystal structure of FCC Au.

Parameters used to build Au lattice is listed below: [3]

space group: Fm-3m (No.225)

Cell parameters:

a: 4.0782 Å

b: 4.0782 Å

c: 44.0782 Å

α: 90.000°

β: 90.000°

γ: 90.000°

 Figure 2.1 Energy cutoff convergence test for bulk Au.

 Figure 2.2 K-point convergence test for Bulk Au. Blue and red dots are data points acquired from DFT calculation and the curves are B-spline fitting of corresponding data points.

Our calculation used plane-wave bases with on the fly generated ultrasoft (OTFG-ultrasoft) pseudopotentials in CASTEP. PBE-GGA was used as the functional. Convergence tests for both k-point mesh and energy cutoff are done. The results of convergence tests is shown in Figures 2.1 and 2.2. Blue and red dots are data points acquired from DFT calculation and the curves are B-spline fitting of corresponding data points. Total energy for the conventional cell is calculated using a k-point mesh of 10*10*10 and an energy cutoff of 500 eV. The core radius for ultrasoft pseudopotential for Au is 2.4 Bohr (~1.27 Å). The ultrasoft pseudopotential was generated with 25 valence electrons (4f14 5s2 5p6 5d10 6s1). The calculated total energy is -56515.8988 eV and energy for each gold atom is -14128.9747 eV. We used equation Ecoh = ( Ebulk – N*Eatom ) / N to calculate cohesive energy. Ebulk is the final energy for bulk Au, Eatom is the final energy for Au atom and N is the number of atoms in the unit cell. The cohesive energy here is -5.6961 eV.

(100) and (111) surfaces

FCC unit cell was built and was cleaved along the Miller indices of (100) and (111). The thickness of atom layers was chosen to be 2.5 for (100) surface and 5 for (111) surface so that five layers of Au atoms were in the slab. This difference in thickness is due to the different stacking style of Au atoms in different directions. The vacuum distance was set to 10 Å to prevent interactions between the top layer and the bottom layer of the slab. The bottom two layers of Au atoms were chosen to be frozen layers to simulate bulk Au crystal. The slabs are shown in Figure 3.

Figure 3 Structure of (111) slab (upper) and (100) slab (lower).

Convergences tests for k-point mesh and energy cutoff were done for both (100) slab and (111) slab. Results were shown in Figures 4.1, 4.2, 5.1, 5.2. Three to four calculations were done for each convergence tests here due to the limited computational power of my laptop. However, it is clear the both (111) slab and (100) slab converges with a k-point mesh of 9*9*1 and an energy cutoff of 500 eV.

Figure 4.1 Energy cutoff convergence test for (100) slab

Figure 4.2 K-point convergence for Au (100) slab.

Figure 5.1 Energy cutoff convergence for Au (111) slab.

Figure 5.2 K point convergence for Au (111) surface.

The energy was calculated for (100) surface slabs with different vacuum thickness (10 Å, 15 Å, 20 Å). The result of surface energy dependence on vacuum thickness is shown in Figure 6. As can be seen, the convergence problem does exist here. However, the difference in surface energy decrease to below 0.0005 eV after vacuum thickness reached 15 Å. Thus, we chose to use a vacuum thickness of 15 Å to make sure both accuracy and calculation time are in an acceptable range.

Figure 6 Surface energy dependence on vacuum thickness for Au (100) slab.

Attempts on calculating layer dependence of energy were also made. However, after we increase atom layers in the slab, the calculation became extremely slow. We waited for a long period of time but the second optimization step in geometry optimization still didn’t show up. Thus, we have to stop the calculation.

For (111) surface, we simply change vacuum thickness to 15 Å to determine its surface energy due to the extremely long calculation time.

The equation we use to calculate surface energy is from “Density Functional Theory, A Practical Introduction” book:

σ(surface)= 1/2A*[E(slab)-n*E(bulk)]

Here, σ(surface) is the surface energy per area. A is the area of used slab and n is the number of atoms in the slab. E(slab) is the total energy of the slab and E(bulk) is the energy per atom in bulk Au.

The final result is listed below:

The surface energy for (100) surface is 0.0544 eV/Å^2.

The surface energy for (111) surface is 0.0316 eV/Å^2.

Here we compare to the surface energy calculated by Tran, R., Xu, Z., Radhakrishnan, B. et al. Their surface energy for (100) surface is 0.054 eV/Å^2 and their surface energy for (111) surface is 0.044 eV/Å^2. The difference in surface energy for (111) surface may due to difference in vaccum thickness and atom layer.

The formation of a new surface by cleaving form a bulk crystal requires the breaking of atom bonds on two side of the cleavage plane. Surface with lower energy is usually considered more stable since lower surface energy means smaller amount of work to be done when cleaving along this surface. Thus, the Au (111) surface is energetically more favored than the (100) surface. This is reasonable since the more closely stacked surface is usually more favored.

References:

[1] Tran, R., Xu, Z., Radhakrishnan, B., Winston, D., Sun, W., Persson, K. A., & Ong, S. P. (2016). Surface energies of elemental crystals. Scientific data, 3(1), 1-13.

[2] Sholl, David, and Janice A. Steckel. Density functional theory: a practical introduction. John Wiley & Sons, 2011.

[3] Maeland, A., & Flanagan, T. B. (1964). Lattice spacings of gold–palladium alloys. Canadian journal of physics, 42(11), 2364-2366.

[4] Tran, R., Xu, Z., Radhakrishnan, B. et al. Surface energies of elemental crystals. Sci Data 3, 160080 (2016).

Author: Mihir Parekh

Abstract: In order to obtain the surface energies of 001 and 100 surfaces of zinc crystal, plane wave basis set calculations have been performed by using CASTEP in Materials Studio. Generalized Gradient Approximation (GGA) was used along with Perdew-Burke-Ernzerhof (PBE) exchange correlation functional and On The Fly Generated (OTFG) ultrasoft pseudopotential. The surface energies of 001 and 100 surfaces are 0.667 J⁄m² and 0.962 J⁄m². Thus, the 001 surface is more stable compared to the 100 surface.

Introduction:

Metal anode batteries have been gaining traction due to their high energy densities. Along with alkali metal anodes, zinc anodes have also been investigated thoroughly in the past [1]. While charging, just like alkali metals, charging current densities above the limiting current densities leads to unstable plating which is the cause for dendrite formation [2]. One of the factors governing the morphology of electrodeposited metal is surface energy. Hence, in this project, surface energies for 001, and 100 surfaces of zinc were obtained.

Methods:

CASTEP [3] was used to perform calculations using plane wave basis set in Materials Studio. Generalized Gradient Approximation (GGA) was used with Perdew-Burke-Ernzerhof (PBE) exchange correlation functionals and On The Fly Generated (OTFG) ultrasoft pseudopotential. A core radius of 2.2 Bohr and an outer shell electronic configuration of 3d10 4s2 was used for calculations. The crystal structure of Zn was imported from the Materials Studio library. Geometry optimization was performed for the imported zinc structure. 1000 eV of cut-off energy and 18×18×16 k points were used for geometry optimization. The lattice parameters obtained from geometry optimization were as follows: a= 2.6649 Å, b=2.6649 Å, c=4.9468 Å, α= 90 , β=90, γ= 120.  The fractional co-ordinates of zinc atoms were (0.333,0.333, 0.25) and (-0.333, -0.667, 0.75). The mentioned lattice parameters were used for optimizing the number of k points and cut-off energy. 0.005 eV was chosen as the convergence criteria for k points and cut-off energy optimization.

Fig. 1 shows the magnitude of change in energy with respect to the previous energy versus the number of k points. It shows that k points more than 520 k points (22×22×20) satisfy the convergence criteria. So, the configuration of 22×22×20 was used in simulations that were done to obtain the converged cut-off energy.

Fig. 1: Optimizing number of k points

Fig. 2 shows the energy relative to the lowest energy obtained during k points and cut-off energy optimization. Cutoff energies more than 500 eV satisfy the chosen convergence criteria. So, 500 eV was chosen as the cut-off energy for further simulations.

Fig. 2 Optimization of cut-off energy

Simulations were performed to derive the surface energies for 001 and 100 surface. Since, the lattice parameters and atom positions in x and y directions are symmetric, the 010 surface gives exactly the same results as 100 surface. So, in this post, only the results of 001 and 100 surface have been mentioned.

In order to calculate the surface energy of 001 and 100 surfaces multiple simulations with varying number of layers were performed. In all calculations, all the layers except the top and bottom layers were frozen. A vacuum slab of 10 Å in ‘c’ direction was used for all simulations. The vacuum slab thickness was varied and the changes were not significant. So, 10 Å of vacuum slab thickness was used for all simulations. 22×22×1 and 22×20×1 were the k point configurations that were used for simulations related to 001 and 100 surfaces respectively.   These configurations were chosen based on the converged k point configuration. In the direction of vacuum slab, only 1 k point was chosen to minimize the computational expense. The surface energy was calculated by using the formula,

σ= (E-n×E_bulk)⁄(2×A),

where E is the energy of the structure with the cleaved surface, E_bulk is the energy of the bulk material per atom,  n is the number of  atoms, and A is the area of the surface. 0.0005 eV⁄Ų was chosen as the tolerance criteria. In order to account for the order of magnitude change due to division by surface area, the chosen tolerance criteria is 1 order of magnitude lower than the tolerance criteria that was chosen while optimizing k points and cut-off energy. The surface energies (converted to J⁄m²) are shown in Fig. 3.

Fig. 3: Surface energies of 001 and 100 surface

The converged surface energies for 001 and 100 surfaces are 0.667 J⁄m² and 0.962 J⁄m² respectively. The converged surface energies were compared with the surface energies obtained  from calculations in which all the layers were unfrozen and no significant change was observed. The surface energy obtained in this calculation for 001 surface is roughly 2/3rd of the value (0.99 J⁄m²) obtained from DFT calculations [4] and liquid surface tension measurements [5].  However, Tyson et al. [5] do  not clearly mention the surface for which the calculations were done.  Vitos et al. [4] use Full Charge Density (FCD) in GGA method to calculate the surface energies. However, the calculations performed for this post use the ‘Self Consistent Field’ (SCF) technique and that might be a reason for the difference. Moreover, the lattice constant used by Vitos et al. [4] (2. 684 Å, with c/a=1.86) is slightly different than the lattice constant used in our calculations (2.665 with c/a=1.856). This also may lead to some minor difference. Moreover, Zheng et al. [1] show that zinc electrodeposition leads to the creation of zinc platelets which are perpendicular to the [001] direction i.e. zinc prefers to deposit as [001] surface. This indicates that [001] surface should have a lowest surface energy amongst all zinc surfaces. This agrees qualitatively with our calculations as we obtain that 001 surface has a lower surface energy than 100 surface. The number of layers required to achieve converged surface energies were 15 and 17 for 001 and 100 surfaces respectively.

Conclusion:

Since, the surface energy of 001 surface is lower than the surface energy of 100 surface, 001 surface is more stable compared to 100 surface.

References:

[1] Zheng, J., Zhao, Q., Tang, T., Yin, J., Quilty, C. D., Renderos, G. D., … & Jaye, C. (2019). Reversible epitaxial electrodeposition of metals in battery anodes. Science, 366(6465), 645-648.

[2] Desai, D., Wei, X., Steingart, D. A., & Banerjee, S. (2014). Electrodeposition of preferentially oriented zinc for flow-assisted alkaline batteries. Journal of Power Sources, 256, 145-152.

[3] Clark SJ, Segall MD, Pickard CJ, Hasnip PJ, Probert MI, Refson K, Payne MC. First principles methods using CASTEP. Zeitschrift für Kristallographie-Crystalline Materials. 2005 May 1;220(5/6):567-70.

[4] Vitos, L., Ruban, A. V., Skriver, H. L., & Kollar, J. (1998). The surface energy of metals. Surface science, 411(1-2), 186-202.

[5] Tyson, W. R., & Miller, W. A. (1977). Surface free energies of solid metals: Estimation from liquid surface tension measurements. Surface Science, 62(1), 267-276.

by Vishal Jindal

1. Introduction

Thiophene is a heterocyclic aromatic compound with the formula C₄H₄S, consisting of a planar five-membered ring. The thiophene ring is the most widely used building block for the construction of conjugated polymers mainly because of its high environmental stability [1]. Thiophene-based π-conjugated systems have attracted much attention in developing high performance organic solar cells [2]. Here, in this post, we are using Density Functional Theory (DFT) [3] calculations to determine the vibrational frequencies of an isolated Thiophene molecule.

2. Methodology

Before finding the vibrational frequencies of thiophene, we first optimized the structure of the thiophene ring. A model of the thiophene ring was constructed using Material Studio.

Fig 1. Model of Thiophene on Material Studio

In order to use a plane-wave basis set, we need to make the system periodic, so we used a box with side 15 Å and K-point mesh of 1x1x1 for the isolated thiophene molecule (as shown below).

Fig 2. Isolated Thiophene molecule in a periodic box of side 15 Å.

Calculation Details

Structure optimization calculations were performed using the Vienna Ab initio Simulation Package (VASP), a plane wave basis set pseudo-potential code. The exchange and correlation energies were calculated using the Perdew, Burke, and Ernzerhof (PBE) [4] functional described within the generalized gradient approximation (GGA) [5]. A plane-wave basis set cutoff energy of 450 eV was used. The calculations were considered optimized when the force on every relaxed atom ( ionic convergence limit) was less than 0.02 eV Å^-1. The electron densities were self consistently solved for energy with the convergence limit (electronic convergence) set to 10^-5 eV.

Fig 3. Structurally optimized thiophene ring using VASP (Molden is the visualization tool used)

The energy of the optimized structure is -54.237 eV (zero-point energy not corrected). After getting the final structure of the thiophene ring, we calculated the vibrational frequencies by varying the displacements (δb) used in finite-difference calculations for vibrational frequencies [3]. The displacement values (δb) used were 0.005, 0.006, 0.01, 0.025, 0.05, 0.07 and 0.1 Å.

3.  Results

Fig 4. The 5th vibrational mode (f5) for the thiophene molecule (~1500 cm^-1)

Thiophene (C₄H₄S) has 9 atoms. So, it has 3N – 6, i.e. 21 vibrational modes. The rest of the 3 rotational modes and 3 translation modes should be theoretically zero, but they have some small or imaginary value because of the small numerical inaccuracies that inevitably occur in calculating a Hessian via finite difference expressions and a calculation method with finite numerical accuracy [3]. Below are the plots of vibrational frequencies as a function of the inverse of displacement value (δb) used for calculating the Hessian matrix.

Fig 5. Plot for vibrational frequency (f) vs inverse displacement value (1/δb) used (f1, f2, f3, f4 modes)

Fig 6. Plot for vibrational frequency (f) vs inverse of displacement value (1/δb) (f5 – f11 modes)

Fig 7. Plot for vibrational frequency (f) vs inverse of displacement value (1/δb) (f12 – f21 modes)

Zero-point energy (ZPE) corrected energy for isolated thiophene was also calculated for various displacement values (δb) using the equation given below taking into account 21 vibrational modes and E_o  = – 54.237 eV:

Table 1. Corrected zero-point energy and the absolute sum of rotational/ translational modes for different displacement values (δb).

4. Conclusion

From the graphs above, we can notice that as the displacement magnitude is lowered, the frequencies of the modes converge to a stable value before they again diverge a little when the displacement value becomes too small. This is expected as harmonic approximation is valid at smaller displacements and can have huge deviations from reality for larger separations. But as the displacement value becomes too small, the energy differences also become too small relative to the SCF convergence which results in slight divergence from the “converged” value.

Also, we know that theoretically, the frequencies of translational/rotational modes should be 0. For DFT vibrational calculations, we can compute the sum of the absolute values of these frequencies (taking sum of last 6 normal modes), and closeness of that sum to 0 can be taken as a measure for better results. We did this for different values of displacement (δb), as shown in the last column of the table in the previous section.

Considering the factors mentioned in the previous two paragraphs, displacement values (δb) of 0.025, 0.01 and 0.006 Å (i.e. 1/δb = 40, 100 and 167 respectively), all give convergent and better results than other values of displacement used to calculate Hessian matrix. Therefore,  0.01 Å can usually be used for getting normal vibrational modes, once we have an optimized structure of the molecule of interest.

5. References

[1] Stergios Logothetidis, Handbook of Flexible Organic Electronics (2015) Materials, Manufacturing, and Applications, Section 4.2.1 (https://doi.org/10.1016/C2013-0-16442-2).

[2] Fan Zhang, Dongqing Wu, Youyong Xu, and Xinliang Feng, Thiophene-based conjugated oligomers for organic solar cells, J. Mater. Chem., 2011,21, 17590-17600.

[3] D. Sholl, J.A. Steckel, Density Functional Theory: A Practical Introduction, Wiley 2009

[4] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 77 (1996) 3865-3868.

[5] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Atoms, Molecules, Solids, And Surfaces – Applications of the Generalized Gradient Approximation for Exchange and Correlation, Phys. Rev. B, 46 (1992) 6671-6687.