Abstract

DFT was used to study the surface diffusion of a silver atom between two hollow sites on the Ag 100 surface. Diffusion by hopping was inspected using the GGA PBE functional and VASP to perform a transition state search. The diffusion barrier was found to be 0.576 eV which agrees well with similar studies.

Introduction

Materials not at absolute zero will have some diffusion of surface atoms. The energy barrier to diffusion is an important property that predicts rates of diffusion.  In particular, these rates are useful for kinetic Monte Carlo (KMC) calculations.  We will inspect a method for finding the minimum energy path of a silver atom diffusing by hopping on a silver (100) surface.  In particular, we will look at the results using the GGA (generalized gradient approximations) PBE (Perdew–Burke-Ernzerhof) functional.  These calculations are made using VASP (Vienna Ab initio Simulation Package).

Method

A 3×3 supercell of 100 is chosen to ensure minimal influence from periodicity.  A slab thickness of 3 atomic layers is chosen with the bottom 2 layers fixed.  These are chosen for computational efficiency.  A vacuum slab of 12 Angstroms is chosen.

We use PAW (Projector augmented-wave) potentials [1]. Ag has the electron configuration of 1s2 2s2 2p6 3s2 3p6 3d10 4s2 4p6 4d10 5s1, and the pseudopotential treats 4d10 5s1 as the valence electrons.  The convergence tolerances were chosen to be:  energy at 1.0e-5 eV/atom, force at 0.05 eV/Å, stress at 0.1 GPa, and displacement at 0.002 Å.  A basis cutoff energy of 900 eV is chosen.

The initial and final states are geometry optimized first, then a 5 frame trajectory is made.  An odd number of frames is chosen to avoid missing the peak of the energy barrier by symmetry.

Calculations were run for various k-point values against the barrier height to test for energy convergence.

We see here that by 242 k-points we have reached k-point convergence within 1.0e-2 eV Which is a reasonable level of accuracy for our purposes.

The calculation results at 242 k-points are shown below.

We find the energy peak to be 0.576 eV, which is in reasonably good agreement with similar studies [4] which found a value of 0.53 eV: a difference of about 9%.

Conclusion

While the results using the parameters chosen gave good results, for additional accuracy, a thicker slab should be used.  The slab thickness used was chosen to cut down on computation time.  Additionally, convergence with respect to the cutoff energy should be performed beyond 900 eV to confirm energy convergence.

References

1. P.E. Blöchl, “Projector augmented-wave method”, Phys. Rev. B 50, 17953 (1994).
   G. Kresse, and J. Joubert, “From ultrasoft pseudopotentials to the projector augmented wave method”,
   Phys. Rev. B 59, 1758 (1999).
2. Density Functional Theory: A Practical Introduction. (2009)  David S. Sholl, Janice A. Steckel
3. https://periodictable.com/Elements/047/data.pr.html
4. “Anisotropy of Growth of the Close-Packed Surfaces of Silver” Yu, Byung Deok, Scheffler, Matthias, PhysRevLett.77.1095, p 1095-1098

Author: Fan Zhang

Introduction

Bulk \(FeSe\) is a conventional superconductor with PbO structure while \(FeTe_{1-x}Se_{x}\) is a candidate for topological superconductor.[1][2] It is useful to perform first principle density functional theory (DFT) calculation on the band structure and density of states of \(FeSe\) and \(FeTe_{1-x}Se_{x}\) with various x values since they provide information such as the role of Te on the whole structure, carrier density of the whole material, whether there would be critical feature near Fermi level. In order to ensure our DFT calculation could provide any useful insight on the real material, it is essential to choose the correct exchange correlation (XC) functional. This work aims at testing the band structure and density of states calculation by different functionals to fix the best XC functional among chosen ones. Also provide a simple step on the effect of Te substitution on the band structure and density of state.

Method

The code used in this work is from CASTEP implementing plane-wave density functional theory. The XC functional used in the following calculations are Generalized Gradient Approximation (Perdew-Burke-Ernzerholf), Generalized Gradient Approximation (Perdew-Wang 91), and Local Density Approximation (Ceperley–Alder–Perdew–Zunger). The pseudopotential used is OTFG ultrasoft which treats 3d6 4s2 for Fe, 3d10 4s2 4p4 for Se, 5s2 5p4 for Te as the valence electrons. Spin polarization is taken into account in all cases. The self-consistent-field (SCF) tolerance of the calculated energy is \(2^{-6}\) eV with 400 maximum cycles. The smearing of the orbital occupancy is fixed at 1eV with 20% empty bands. K points sampling used in the calculation is 5x5x4 (12 irreducible k points for the most symmetric structure). The cut off energy is 440eV. The convergence of both k points sampling and cut off energy were achieved with  \(\Delta\) E \(\leq\) 0.01eV.

Result

1.Convergence of cut off energy and k points sampling

The convergence of cut off energy and k points sampling is done by using the optimized FeSe crystal with PbO structure. The k points configuration used is 5x5x4 mesh as default when performing the convergence of cut off energy. Then the ideal 440eV of cut off energy is used in the convergence of k points sampling. It is found that a 440eV cut off energy with 5x5x4 mesh k points sampling (corresponding to 12 irreducible k points in this structure) could give \(\leq\) 0.01eV. The corresponding graph of these data is also shown in Figure 1 and Figure 2.

2. Structures used in this work

This work used three different structures. Each structure is  geometrically optimized with convergence energy tolerance to bw \(2^{-5}\) eV/atom, the maximum force to be 0.05eV/A and the maximum displacement to be 0.002A.

The top view, side view, and overall view of the structure of optimized FeSe with PbO are shown in Figure3. In this structure, a=b=3.947A, c=4.988A,  \(\alpha\) = \(\beta\) = \(\gamma\) = \(90.000^{\circ}\). Fe atoms are put on sites (0, 0, 1) and (0.5, 0.5, 1). Se atoms are put on sites (0, 0.5, 1.2385) and (0.5, 1, -1.2385).

By replacing one \(Se\) atom by \(Te\), \(FeTe_{0.5}Se_{0.5}\) is achieved. The top view, side view, and overall view of the structure of optimized FeTe0.5Se0.5 are shown in Figure 4. In this structure, a=b=3.947A, c=4.988A,  \(\alpha\) = \(\beta\) = \(\gamma\) = \(90.000^{\circ}\). Fe atoms are put on sites (0, 0, 0) and (0.5, 0.5, 0). Se atom is put on site (0.5, 0, 0.7615). Te atom is put on site (0, 0.5, 0.2385).

By combining 4 conventional cell of \(FeSe\) into a supercell and replacing one of the eight \(Se\) atoms by \(Te\), \(FeTe_{0.125}Se_{0.875}\) is achieved. The top view, side view, and overall view of the structure of optimized FeTe0.125Se0.875 are shown in Figure 5. In this structure, a=b=7.894A, c=4.988A,  \(\alpha\) = \(\beta\) = \(\gamma\) = \(90.000^{\circ}\). Fe atoms are put on sites (0, 0, 0), (0.25, 0.25, 0), (0.5, 0, 0), (0.75, 0.25, 0), (0, 0.5, 0), (0.25, 0.75, 0), (0.5, 0.5, 0) and (0.75, 0.75, 0). Se atoms are put on site (0, 0.25, 0.2385), (0.25, 0, -0.2385), (0.75, 0, -0.2385), (0, 0.75, 0.2385), (0.25, 0.5, -0.2385), (0.5, 0.75, 0.2385), and (0.75, 0.5, -0.2385). Te atom is put on site (0.5, 0.25, 0.2385).

3.Band structure of \(FeTe_{0.5}Se_{0.5}\) calculated by different XC functionals

The reciprocal lattice and various directions used in this work is shown in Figure 6. In order to compare with previous work, we focus on GZ direction (usually called \(\Gamma\)Z direction).[2]

The band structure calculated by various XC functionals are shown in Figure 7 and the calculation done by Peng Zhang et al. are shown in Figure 8. [2]

GGA PBE result is the only one that shows a band cross feature along \(\Gamma\)Z direction ,the only one shows spin polarized structure, and is the only one does not have a band gap. It is thus clear that GGA PBE is the most suitable XC functionals among the three tested. The difference between this calculation and the previous work may caused by spin-orbit coupling (SOC). In this work, we did not count SOC while SOC is claimed to be the reason for the opened gap in the paper by Peng Zhang et al.[2]

For the following analysis, we will choose GGA PBE as our XC functional.

4.Effect of Te substitution on band structure and density of states

The band structure of \(FeSe\), \(FeTe_{0.125}Se_{0.875}\), and \(FeTe_{0.5}Se_{0.5}\) are shown in Figure 9.

The density of states of \(Fe\), \(Se\) and \(Te\) in \(FeSe\), \(FeTe_{0.125}Se_{0.875}\), and \(FeTe_{0.5}Se_{0.5}\) are shown in Figure 10, Figure 11, and Figure 12.

It could be seen clearly from the calculation that \(Te\) substitution greatly changes the band structure of original \(FeSe\) instead of just moving chemical potential. By substituting a small amount of \(Te\), density of states around fermi surface could be doubled, which indicate an increase in the carrier density. However, increasing \(Te\) will decrease the density of states around fermi surface again. This change in density of states is contributed mostly by p and d orbitals of \(Fe\).

5.Comparison of calculation time of various XC functionals

The calculation time of band structures and density of states by different XC functionals are shown in Table 1. Results are normalized with respect to the shortest time.

GGA PBE is again the best choice for doing this calculation as it requires the least time per iteration.

Conclusion

By comparing the band structure of \(FeTe_{0.5}Se_{0.5}\) calculated using various XC functionals with the result in previous work, GGA PBE is the best choice among all three XC functionals chosen here. Comparison of computational time also suggests GGA PBE to be the best choice. By the calculation done with GGA PBE, adding \(Te\) to \(FeSe\) can greatly change the band structure. As the weight of \(Te\) increasing, the carrier density is expected to increase following by a decrease. The difference between the band structure of this work and previous work indicate the importance of SOC in this material.

For the next step of this work, I would add SOC to my material and increase the number of k points when calculating band structures to allow a more detailed comparison. Also, I would proceed to more possible x values in \(FeTe_{1-x}Se_{x}\) to see how \(Te\) would affect the material. Then, transport properties could be calculated to allow more possible comparisons with experiments.

References

by Vishal Jindal

1. Abstract

This post calculates the band structure and density of state for an isolated polythiophene chain using density functional theory (DFT) and compares the various electronic properties such as bandgap (Eg), bandwidth and effective charge mass of the valence and conduction bands for different exchange-correlation functionals.

2. Introduction

Polythiophene (PT) is an organic polymer which has the simplest structure among a general class of polymerized thiophenes known as polythiophenes (PTs). In the presence of electron-acceptor materials PTs become conductive and demonstrate interesting optical properties resulting from their conjugated backbone which have made them standard materials for organic semiconductor devices [1,2]. The development of polythiophenes and related conductive organic polymers was recognized by the awarding of the 2000 Nobel Prize in Chemistry to Alan J. Heeger, Alan MacDiarmid, and Hideki Shirakawa “for the discovery and development of conductive polymers” [3]. Electronic transport properties, and optical properties such as light absorption and emission, are dictated by the valence and conduction bands, which are the closest bands to the Fermi level, and the gap between them. The undistorted PT chain has a 1-D semiconductor band structure [4] and the reported experimental value for bandgap is 2.0eV [8].

3. Methodology

We used a plane-wave basis set with ultrasoft pseudopotentials as implemented in the CASTEP code [5] to perform DFT calculations to get band structure and density of states for different exchange-correlation functionals such as LDA CA-PZ, GGA-PBE, GGA-PW91, and GGA-BLYP. Materials Studio was used as a builder, visualizer, and user interface for the CASTEP calculations.

First, the structure of the PT chain was optimized using the Perdew, Burke, and Ernzerhof (PBE) [6] functional described within the generalized gradient approximation (GGA) [7] and the same optimized structure was used to get band structure using different functionals. Based on the work by Joel et al [4], a basis set cutoff energy of 450 eV and k-point sampling of 16 x 1 x 1 (instead of 7 x 1 x 1 used in the paper to increase the k-point sampling) were used. SCF tolerance is taken to be 2.0e-6 eV/atom for all the calculations. The “on the fly” generated ultrasoft pseudopotential for H has a core radius of 0.6 Bohr (~0.32 Å), 1.4 Bohr (~0.74 Å) for carbon (C)  and 1.8 Bohr (~0.95 Å) for sulfur (S) and was generated with 4 valence electrons (2s2 2p2) for C and 6 valence electrons (3s2 3p4) for S.

While calculating the band structures, an approximate separation of 0.005 1/Å between k-points was used for each of the different exchange-correlation functionals.

4. Results

Band Structure

Density of State

The band structures and density of states for polythiophene using 4 different exchange-correlation functionals are given in Fig 1. and Fig 2., respectively. The energy band gaps, bandwidth and effective charge mass for valence and conduction bands are calculated for each exchange-correlation functional (XC).

The energy between the valence and conduction band of a polymer is related to the lowest energy of its monomer units and to the bandwidth resulting from the overlap between the
monomer orbitals [9]. A bandgap (𝐸g) is defined as the difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) energy levels in the polymer [10]:  𝐸g = 𝐸(LUMO) − 𝐸(HOMO) (eV).

In accordance with the band-like theory, the bandwidth (BW) and electron effective
mass (m*) are good parameters for predicting the hole and electron-transporting ability of polymers [11]. Valence (conduction) bandwidth is defined as the energy difference between the first band located right below (above) the Fermi level at Γ-point (k=0) to the M-point (k=1) [12].

The effective mass (m*) of the carrier at the band edge representing mobility was
obtained as the square of ħ (= h/2π) multiplied by the reciprocal of the curvature from E(k) with k (wave vector), and the formulation is defined as:

The denominator of the above equation is calculated by fitting a parabola (y = kx² + c) using the first 5 data points (which lie between k = 0 to 0.33) of energy vs k plot (band structure) for both conduction and valence for each of the different functional used.

Table 1. (below) shows the bandgap, bandwidths, and effective charge mass calculated using different exchange-correlation (XC) functionals.

Correlation functional (XC)	Band gap (eV)	Valence band effective mass (-m*)	Valence bandwidth (eV)	Conduction band effective mass (m*)	Conduction bandwidth (eV)
LDA CA-PZ	1.47	1.23	1.73	1.46	1.39
GGA PBE	1.49	1.25	1.72	1.47	1.39
GGA PW91	1.46	1.25	1.72	1.49	1.38
GGA BLYP	1.46	1.25	1.72	1.49	1.38

As evident from the results in the table, values of the energy bandgap, valence and conduction bandwidths, and effective masses for carrier charge in valence and conduction bands are approximately the same whether we use the LDA or GGA exchange-correlation functionals for the DFT calculation of band structure.

5. Conclusion

Local density (LDA) and generalized gradient approximation (GGA) density functionals generally underestimate band gaps for semiconductors by about 40% of the actual value [13, 14], but can be valuable to study the shape of the valence and conduction bands. Hybrid functionals such as B3PW91 and B3LYP are known to be better predict the experimental band gaps for semiconductors and insulators [14]. The experimental value of the bandgap for polythiophene was found to be 2.0 eV, which indeed agrees well with the calculated bandgap in the framework of the hybrid (B3LYP) as shown by Kaloni et al. [12].

CASTEP code as implemented using Material Studio, for the calculation of properties such as band structure, does not allow for use of hybrid functionals along with the ultrasoft pseudopotential. Comparing the above results for bandwidths and effective mass using LDA, GGA functionals with the results from using hybrid functionals for our calculations can strengthen the argument that less computationally intensive calculations (involving LDA, GGA functionals) can be used reliably for studying the band properties other than the energy bandgap.

6. References

[1] G. Li, R. Zhu, and Y. Yang, Nat. Photonics, 2012, 6, 153–161.

[2] S. P. Rittmeyer and A. Groß, Beilstein J. Nanotechnol., 2012, 3, 909–919.

[3] https://en.wikipedia.org/wiki/Polythiophene

[4] Joel H. Bombile, Michael J. Janik and Scott T. Milner, Tight binding model of conformational disorder effects on the optical absorption spectrum of polythiophenes, Phys. Chem. Chem. Phys., 2016, 18, 12521.

[5] J. Clark Stewart, D. Segall Matthew, J. Pickard Chris, J. Hasnip Phil, I.J. Probert Matt, K. Refson, C. Payne Mike, First-principles methods using CASTEP,  Zeitschrift für Kristallographie – Crystalline Materials, 220(5-6) pp. 567-570 (2005).

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

[7] 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.

[8] Kobayashi M, Chen J, Chung T-C, Moraes F, Heeger A J, Wudl F. Synth Met. 1984;9:77–86. doi: 10.1016/0379-6779(84)90044-4.

[9] U. Salzner, J. B. Lagowski, P. G. Pickup, and R. A. Poirier, Design of low band gap polymers employing density functional theory—hybrid functionals ameliorate band gap problem, Journal of Computational Chemistry, vol. 18, no. 15, pp. 1943–1953, 1997.

[10] E. Bundgaard and F. C. Krebs, Low band gap polymers for organic photovoltaics, Solar Energy Materials and Solar Cells, vol. 91, no. 11, pp. 954–985, 2007.

[11] X.-H. Xie, W. Shen, R.-X. He, and M. Li, A Density Functional Study of Furofuran Polymers as Potential Materials for Polymer Solar Cells, Bulletin of the Korean Chemical Society, vol. 34, no. 10, pp. 2995–3004, Oct. 2013.

[12] Thaneshwor P. Kaloni, Georg Schreckenbach, and Michael S. Freund, Band gap modulation in polythiophene and polypyrrole-based systems, Sci. Rep. 6, 36554; doi: 10.1038/srep36554 (2016).

[13] J. P. Perdew, Density functional theory and the band gap problem, Int. J. Quantum Chem., 1996, 28, 497–523.

[14] X. Hai, J. Tahir-Kheli and W. A. Goddard, Accurate Band Gaps for Semiconductors from Density Functional Theory, J. Phys. Chem. Lett., 2011, 2, 212–217.

by Da Zhou

Abstract

In this post, ab initio density functional theory (DFT) calculations with different functionals and spin-orbit coupling (SOC) parameters were performed to study the band structure of bulk and monolayer WSe2.

Introduction

Recently, more and more research interest has been drawn into the monolayer transition metal dichalcogenides (TMDs) for their novel physical properties and potential device applications. For example, the bulk MX2 (M=Mo, W; X=S, Se) has indirect bandgaps whereas the monolayer MX2 has direct bandgaps located at the K points. Therefore it is highly desirable to use DFT to numerically calculate the band structure of the TMDs and compare the calculation results with the experimental ones. Here in this post, the band structure of bulk and monolayer WSe2 were numerically calculated by the plane-wave based CASTEP, and local density approximation (LDA) and generalized gradient approximation (GGA) exchange-correlation functionals were used.

Method details

For the bulk WSe2, the lattice constants were 3.327 angstroms in a/b axis and 15.069 in c axis (the data was directly imported from materialsproject.org and no relaxation was performed). The monolayer slab was cleaved in the 001 direction and a vacuum layer of 20 angstroms was added to eliminate interlayer interactions.

Based on the work by  Julia Gusakova et al [1], 50 Ry which is about 680eV was enough for the energy cutoff. The energy cutoff for all of the calculations in this post was 700eV. For both bulk and monolayer WSe2 band structure calculations, a sampling separation of 0.015 1/angstrom was used. SCF tolerance was 1.0e-5 eV/atom and electronic minimizer was all bands/EDFT. Usually, the density mixing option is more recommended for the choice of electronic minimizer. But to perform calculations with the SOC effect included, all bands/EDFT were required. For the purpose of consistency, the calculations without SOC also used all bands/EDFT as the electronic minimizer to make sure the comparison of results being meaningful. The Monkhorst-Pack k-point set was 15 by 15 by 3 for the bulk and 16 by 16 by 2 for the monolayer. And norm-conserving pseudopotentials were used. The pseudo atomic calculation performed for Se: 4s2 4p4. The pseudo atomic calculation performed for W: 5d4 6s2.

Results

A series of calculations were performed on the bulk and monolayer WSe2, with LDA CAPZ or GGA PBE functionals and also with or without SOC. Only the two bands closet to the Fermi energy were plotted in the figures. And the results were summarized at the table below.

It is quite exciting to see that all the calculations, other than the one with GGA PBE functional and without SOC, successfully captured the indirect(direct) bandgap feature for the bulk(monolayer). It is also noticeable that after including SOC, the bandgap tends to be about 0.3eV smaller than the bandgap without SOC for both bulk and monolayer, which makes the calculated band gap values for the bulk very close to the experimental reference[1]. For the monolayer, the results will not be discussed here in this post since there is the experimental controversy over its bandgap being direct or indirect[5], and also the value of its bandgap varies from different measurements[1][6].

Conclusion

Most GGA and LDA calculated band structures matched well with the traditional experimental reference in terms of the bulk WSe2 having indirect bandgap and the monolayer WSe2 having direct bandgap located at the K points. And within the numerical accuracy, the calculated band gap values for the bulk were very close to the reference value when SOC was included. The work can be more meaningful if more exchange-correlation functionals can be used, and if the experimental reference results for the monolayer WSe2 can converge. Nonetheless, the calculated band structures of the monolayer WSe2 in this post are all direct bandgaps, which certainly supports one kind of opinion.

References

[1] Julia Gusakova et al, Electronic Properties of Bulk and Monolayer TMDs: Theoretical Study Within DFT Framework (GVJ‐2e Method), Phys. Status Solidi A 2017, 214, 1700218

[2] Hohenberg, P. and Kohn, W., Inhomogeneous electron gas, Phys. Rev. 1964

[3] Kohn, W. and Sham, L. J., Self-consistent equations including exchange and correlation effects, Phys. Rev. 1965

[4] Clark, S. J., and Segall, M. D., and Pickard, C. J. and Hasnip, P. J. and Probert, M. J. and Refson, K. and Payne, M.C., First-principles methods using {CASTEP}, Z. Kristall. 2005

[5]Hsu, W., Lu, L., Wang, D. et al. Evidence of indirect gap in monolayer WSe2. Nat Commun 8, 929 (2017)

[6]Zhao, W., Huang, Y., Shen, C. et al. Electronic structure of exfoliated millimeter-sized monolayer WSe2 on silicon wafer. Nano Res. 12, 3095–3100 (2019)

Author: Mihir Parekh

Abstract: In this post, Density Functional Theory (DFT) was used to study the surface diffusion of Li adatom between two adjacent hollow sites on Li 001 surface. Two diffusion mechanisms, namely hopping and substitution, were studied by using Perdew-Burke-Ernzerhof (PBE) exchange correlational functional. Complete Linear/Quadratic Synchronized Transient method (Complete LST/QST) was used to perform the transition state search. Vibrational analysis of the transition states gave one imaginary frequency along the reaction co-ordinate, which effectively confirmed that the obtained states were indeed transition states. The diffusion barriers were found to be 0.271 and 0.667 eV for hopping and substitution, respectively, suggesting that hopping is more favorable compared to substitution.

Introduction: In this post, the self diffusion of lithium has been studied on 001 surface by using Density Functional Theory (DFT).  Lithium was chosen as a subject of study because of its high importance in high energy density lithium metal batteries [1]. This is because lithium metal is not only the lightest metal but it also has a highly negative reduction voltage with respect to the standard hydrogen electrode. Dendrites tend to grow from the lithium metal electrode surface during charging and self diffusion of lithium adatoms on lithium metal surfaces is an important phenomenon that can help in controlling dendrite growth. The 001 surface (which is equivalent to 100 surface for BCC lattice) was chosen as it is the most stable lithium surface and thus on electrodeposition Li would tend to deposit as 001 surface.

Methods: Lithium crystallizes as BCC with a lattice constant of 3.51 Å [2] at room temperature. For the simulations, a BCC lattice with a lattice constant of 3.46 Å was used as an input and geometry optimization was performed by using Generalized Gradient Approximation- Perdew Burke Ernzerhof (GGA-PBE) exchange correlation functional in the DMol package. A Double Numerical+ D functions (DND) basis set with the 3.5  version of the basis set file was used for the geometry optimization. The optimization was performed using 4x4x4 and 12x12x12 k point configuration and both yielded a BCC lattice with a lattice constant of 3.51 Å. Hence, 3.51 Å was chosen as the lattice constant.

Various k point configurations were tried out for the optimized geometry till convergence was achieved for energy values. A convergence criteria of 0.0001 Ha was used. The SCF convergence criteria was set to 0.00001 eV. Fig. 1 shows that 10x10x10 configuration is sufficient enough to satisfy the convergence criteria. The y axis in Fig. 1 plots the energy relative to the energy for 11x11x11 k point configuration.

Fig. 1: Optimization of k points

A 001 surface was cleaved from  the optimized geometry. The number of layers in the slab were varied till convergence was achieved and a vacuum slab of 10 Å was used for all simulations. 0.0017 J/m² was used as the convergence criteria for surface energy. This convergence criteria corresponds with the 0.0001 Ha convergence criteria that was used earlier, when the surface area and other conversion factors are taken into account. As shown in Fig. 2, convergence was achieved for 9 layers and the converged surface energy was about 0.54 J/m² which matches well with the surface energy of 0.47 J/m² calculated from DFT calculations by Gaissmaier et al. [1]. The slight difference in the values could be explained by the fact that Gaissmaier et al. [1] had frozen the bottom two layers during their surface energy calculation, whereas in the calculations performed for this post all the layers were left unfrozen. Also, freezing the bottom 2 layers introduces an external asymmetry in the system which can cause the surface energy to change slightly.

Fig. 2: Surface energy versus number of layers for 001 surface of Li

A 2x2x1 supercell was created and a lithium adatom was then added to the hollow site above the top layer. All the 9 layers were frozen and geometry optimization was performed to find the optimum position of the lithium adatom in the hollow site. The optimized structure was then used to study the self diffusion of Li adatom from the hollow site to its neighboring hollow site on the Li 001 surface. The reactant and product have been shown below in figures 3 and 4 respectively.

Fig. 3 Reactant

Fig. 4: Product

The self diffusion process has two possible mechanisms, namely, hopping and substitution. Complete Linear Synchronized Transient/ Quadratic Synchronized Transient (LST/QST) search was used to search for the transition states for both the mechanisms. Vibrational analysis was performed on the obtained transition state structure and an imaginary frequency was obtained for vibration along the direction of reaction co-ordinate. The structures of transition states for the hopping and substitution process have been shown in figures 5 and 6, respectively.

Fig. 5: Transition state for the hopping process

Fig. 6: Transition state for the substitution process

Diffusion barriers of 0.271 eV and 0.667 eV were obtained for the hopping and substitution processes respectively. These values differ from the barrier of 0.04 eV and 0.14 eV obtained via Climbing Image -Nudged Elastic Band (CI-NEB) method by Gaissmaier et al. [1].

Discussion: Although, the diffusion barriers are different when compared to the values obtained by Gaissmaier et al. [1], hopping has been found to be favorable compared to substitution in both calculations. Different supercell sizes could possibly be the reasons for difference in answers.

Another thing to note is that the vibrational analysis gave more than 1 imaginary frequency. However, there was a single imaginary frequency in the direction of reaction co-ordinate. The extra imaginary frequencies could possibly have come from the fact that even the Li atoms in the slab were allowed to vibrate. Moreover, harmonic approximation may not be a good approximation to calculate vibrations and that could also give rise to the extra imaginary frequencies.

Conclusion: Between hopping and substitution, hopping is a more favorable mechanism for self diffusion of Li adatom on Li 001 surface.

References:

1. Gaissmaier, D., Fantauzzi, D., & Jacob, T. (2019). First principles studies of self-diffusion processes on metallic lithium surfaces. The Journal of chemical physics, 150(4), 041723.

2. https://www.webelements.com/lithium/crystal_structure.html

Author: Jeremy Hu

Abstract

Density functional theory (DFT) calculations of TiO₂ anatase, a commonly used catalyst support, typically requires a DFT+U correction term to account for the electron self-interaction error in Ti. DFT as implemented in the Vienna Ab Initio Simulation Package (VASP) was used to determine the minimum value of the Hubbard’s U parameter to accurately represent the electronic density of states of TiO₂ and partially-reduced TiO₂. This study determined that a minimum Hubbard’s parameter of U = 3 eV is sufficient to accurately determine the density of states of TiO₂ and partially-reduced TiO₂.

Introduction

Titanium dioxide (TiO₂) is a common metal oxide for depositing catalytically-active metals onto, such as those used in biochemical production [1].  The metal oxide support provides structure for the active metal sites to be deposited onto. TiO₂ supported-catalysts for biochemical syntheses generally undergo reduction conditions, with the formation of oxygen vacancy sites in the presence of H₂ [2]. The anatase polymorph of TiO₂ most readily forms oxygen vacancies, with the (001) facet widely considered the most catalytically active [1]. Understanding the chemistry of TiO₂ using DFT typically requires a DFT+U correction term to account for the self-interaction error between electrons. This correction term is typically necessary for Ti due to its numerous d-orbital valence electrons. Thus, the DFT+U correction penalizes the delocalization of d-orbital electrons in Ti [3]. Previous literature suggests that the U parameter can range anywhere from 2-10 eV, with the U parameter “sufficient” when the density of states of TiO₂ behaves as an insulator with KS orbitals appearing in the band gap in the presence of oxygen vacancies [3]. Thus, understanding the effects of the U parameter on the density of states of TiO₂ could bring insight into the minimum value of U necessary to confirm the localization of electrons in bare TiO₂ and TiO₂ under oxidation conditions.

Methods

Electronic Methods

Density functional theory (DFT) analysis of TiO₂ anatase (001) was calculated using the plane-wave basis set in the Vienna Ab Initio Simulation Package (VASP) [4]. The Perdew–Burke-Ernzerhof (PBE) exchange correlation functional was used [5]. PBE+D3 was used for dispersion corrections and the projector augmented-wave method (PAW) corrected for core-valence interactions [6] [7]. The Hubbard’s parameter (U) for the DFT+U correction was iterated for Ti from U = 0-3 eV [8]. Each structure was reoptimized for each value of U before calculating the density of states. The forces on the atoms from the geometric optimizations used a convergence criteria of < 0.05 eV/Å. The self-consistent field tolerance for all calculations was 10^-5 eV.

Similar work on TiO₂ suggest that a Monkhorst-Pack k-point mesh of 3 x 3 x 1 and a cutoff energy of 450 eV are above the minimum for convergence [9] [10]. The valence electrons considered for each atom type were O (2s² 2p⁴) and Ti (3s² 3p⁶ 4s² 3d²).

A vacuum space of 15 Å between slabs was used to minimize dipole interactions in the z-direction (i.e., normal to the surface). A 2 x 2 supercell of anatase TiO₂ (001) was used for the DFT calculations, with the bottom three atomic layers fixed while the rest of the atoms were allowed to freely relax (Fig. 1). The slab’s thickness was two layers thick, with a single layer defined as the minimum thickness of atoms to have a stoichiometric ratio of Ti_xO_2x. Additionally, the termination of the surface was chosen to be the same as the termination on the bottom of the slab (i.e., oxygen atoms on the top and bottom) to minimize any large dipole moments that would otherwise occur through asymmetry.Figure 1. A 2 x 2 supercell of TiO₂ anatase (001) used for the density of states calculation.

The second-coordinated surface oxygen atom on the surface, which required the lowest amount of energy to remove, was removed for the density of states calculations of anatase TiO₂ with an oxygen vacancy (Fig. 2)

Figure 2. A 2 x 2 supercell of TiO₂ anatase (001) with a second-coordinated bridging oxygen (O_2C ) removed used for the density of states calculation. The symmetrically-equivalent O_2C behind the removed atom was removed in this figure for clarity.

Results and Discussion

The density of states for the bare TiO₂ anatase (001) was plotted for values of U from 0 to 3 (Fig. 3a-d).

Figure 3a. The total density of states vs. the Fermi level subtracted from the energy (eV) for U=0. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 3b. The total density of states vs. the Fermi level subtracted from the energy (eV) for U=1. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 3c. The total density of states vs. the Fermi level subtracted from the energy (eV) for U=2. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 3d. The total density of states vs. the Fermi level subtracted from the energy (eV) for U=3. The blue line represents the spin up states and the black line denotes the spin down states.

There was slight noise and variation in the density of states when U was iterated from 0 to 3. However, the band gap of approximately ~2 eV stayed relatively constant in all cases, which relatively agrees with the experimentally determined band gap of ~3 eV in TiO₂ [11]. The nuances in the density of states at U=3 (i.e., the clearer differentiation of peaks) suggest that a higher U value may be more appropriate for the calculation of TiO₂ density of states.

The density of states for the partially-reduced TiO₂ (i.e., TiO₂ with a surface oxygen vacancy) was plotted as well from U = 0 to 3 (Fig. 4a-d).

Figure 4a. The total density of states vs. the Fermi level subtracted from the energy (eV) for TiO₂ with a surface oxygen vacancy (O_2C) at U=0. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 4b. The total density of states vs. the Fermi level subtracted from the energy (eV) for TiO₂ with a surface oxygen vacancy (O_2C) at U=1. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 4c. The total density of states vs. the Fermi level subtracted from the energy (eV) for TiO₂ with a surface oxygen vacancy (O_2C) at U=2. The blue line represents the spin up states and the black line denotes the spin down states.

Figure 4d. The total density of states vs. the Fermi level subtracted from the energy (eV) for TiO₂ with a surface oxygen vacancy (O_2C) at U=3. The blue line represents the spin up states and the black line denotes the spin down states.

As expected, there was noise and variation in the density of states ranging from U = 0 to 3. The peak in the band gap only clearly appeared in the cases where U=0 and U=3. As mentioned previously, the density appearing in the band gap for TiO₂ with an oxygen vacancy likely represents occupied KS orbitals following reduction. Recent experimental and DFT studies of TiO₂ with an oxygen vacancy suggest a band gap of approximately 3.0 eV with a gap state appearing around 0.7 eV below the conduction band [12].  In the case where U was set to 0, the gap state was ~0.50 eV below the conduction band, whereas when U=3 the gap state was ~0.70 eV below the conduction band. Since the peak in the band gap more closely matches experimental results in the case where U=3, it is likely that a U parameter of 3 and above is necessary to correctly determine the density of states of TiO₂ with an oxygen vacancy.

 Conclusion

The density of states was calculated for anatase TiO₂ (001) and anatase TiO₂ (001) with a surface oxygen vacancy at varying U correction values from 0 to 3. The data suggest that for bare TiO₂, a Hubbard’s U correction of 3 eV and above may prove appropriate for clearer peak density. However, the band gap of around 2 eV stayed relatively constant in all four cases. In the case of TiO₂ (001) with an oxygen vacancy, the appearance of a gap state near the Fermi level correlates to occupied KS orbitals following the reduction of TiO₂. Although gap states appeared in both cases where U=0 and U=3, the gap state peak which occurred when U=3 more closely agreed with experimental data. Thus, the results suggest that a Hubbard’s U correction of 3.0 eV and above is appropriate for calculating the density of states of TiO₂ (001) and TiO₂ (001) with an oxygen vacancy.

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] Z. Hu, H. Metiu, J. Phys. Chem. C 2011, 115, 13, 5841-5845

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

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

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

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

[8] B. Himmetoglu et al., International Journal of Quantum Chemistry (2014), 114, 14–49

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

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

[11] G. Rocker, J.A. Schaefer, W. Göpel, Phys. Rev. B, 30 (1984), p. 3704

[12] B.J. Morgan, G.W. Watson, Surface Science, 601 (2007) 5034–5041

Abstract

In this post, the transition state search for a elementary reaction in carbon dioxide reduction on Cu(111) surface is performed using CASTEP, and the activation barrier for this reaction is calculated.

Introduction

Electrochemical reduction of CO_2 to valuable chemicals and fuels such as hydrocarbon and alcohols is of great interest. However, experimental results show that most metal electrodes except for copper could only reduce CO_2 mainly to CO. Previous computational work[1][2] showed that the transition state energy of reduction of CO is a good descriptor that could help us to understand the trends of catalytic activity of transition metal surfaces. To understand the process that adsorbed CO reduced to COH on Cu(111) surface, the transition state geometry and energy of *CO reduction to *COH on Cu(111) surface will be calculated and compared to reported results, here * signs are used to show the species are absorbed on the surface.

Method

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

A 2×2 four-layer slab model is used to simulate Cu(111) surface structure with 10Å of vacuum. The optimized bulk with lattice parameter of 3.63Å is used to build the surface model. When the optimization and transition state search are performed, bottom two layers are kept fixed and remaining atoms are relaxed to simulate bulk behaviors. A optimized 5x5x1 k points grid and energy cutoff of 500 eV are kept consistent for all calculations.

Results

k points grid and energy cutoff optimization

With an energy cutoff of 500 eV, single point energies of 2×2 slab models are calculated. In figure 1, the energies of different k points gird relative to 8x8x1 k points grid is shown.

The energy cutoff convergence optimization is carried out with a k point grid of 8x8x1. In figure 2 energies of different energy cutoff relative to 550 eV are shown.

According to those two convergence test results, a 5x5x1 k points grid with an energy cutoff of 500 eV is reasonable for calculation of a 2×2 slab model of Cu(111) surface. Thus these settings are kept consistent for subsequent optimization and search for transition state.

Cu(111) surface optimization

In order to build a reasonable model of reaction happens on Cu(111) surface, the optimization of surface structure is carried out first. A 2×2 Cu(111) surface slab model is used with bottom 2 layers fixed to represent bulk effects on metal surface. BFGS algorithm is used to perform the optimization with energy convergence of 2e-5 ev/atom and force convergence of 0.05 eV/A.

Transition state search

According to previous work[2],  co-adsorption of CO and H on Cu(111) surface will be neighboring hcp sites. In figure 3, co-adsorbed H and CO on neighboring Cu(111) surface is set as reactant geometry. In figure 4, adsorbed COH on Cu(111) hcp sites is set as product geometry. The structures of reactant and product are optimized using the same convergence setting for surface optimization as shown in figure 3 and 4.

The transition state search is performed using a complete LST/QST method implemented in the CASTEP module of Material Studio. The RMS(Root Mean Square) convergence is set to 0.24 eV/A. Followed by a transition state confirmation of higher precision with energy convergence of 2.0e-5 eV/atom, maximum force of 0.1 eV/Å and maximum displacement of 0.005Å. In figure 5, the transition state geometry is shown. Although in both reactant and product structure, adsorbed species are perpendicular to Cu surface, the adsorbed species in transition state structure is rotated towards the initial position of adsorbed H atom.

Figure 6 is the relative energy diagram. Here energy values are reported relative to that of initial reactant state. As shown in the figure, activation barrier  Consequently the activation energy from reactant is found to be 2.48 eV. Comparing to reported activation energy in [2-4], which is less than 1.0 eV, this activation energy is significantly larger. Since in their models, they include explicit water molecules to model what happens in solution. The structure of reactant and product used in this work actually modeled the reaction happens in gas phase. Meanwhile, they also argued that possible proton and electron transfer between adsorbed CO molecules and solvent molecules, which won’t induce surface atoms rearrangement as much as in gas phase.

Conclusion

The transition state geometry of *CO to *COH is found with a reaction barrier from reactant of 2.08 eV. Comparing to the reported electrode model, gas phase reduction of *CO to *COH is significantly more difficult. If given more computational time and resources, similar procedure could be used to address more complex surface reaction network starting from *CO and include explicit solvent molecules to model more realistic reaction conditions.

Reference

[1] Liu, X.; Xiao, J.; Peng, H.; Hong, X.; Chan, K.; Nørskov, J. K. Understanding Trends in Electrochemical Carbon Dioxide Reduction Rates. Nat. Commun. 2017, 8 (1), 15438. https://doi.org/10.1038/ncomms15438.

[2] Hussain, J.; Jónsson, H.; Skúlason, E. Calculations of Product Selectivity in Electrochemical CO2 Reduction. ACS Catal. 2018, 8 (6), 5240–5249. https://doi.org/10.1021/acscatal.7b03308.

[3] Nie, X.; Esopi, M. R.; Janik, M. J.; Asthagiri, A. Selectivity of CO2 Reduction on Copper Electrodes: The Role of the Kinetics of Elementary Steps. Angew. Chemie – Int. Ed. 2013, 52 (9), 2459–2462. https://doi.org/10.1002/anie.201208320.

[4] Nie, X.; Luo, W.; Janik, M. J.; Asthagiri, A. Reaction Mechanisms of CO2 Electrochemical Reduction on Cu(1 1 1) Determined with Density Functional Theory. J. Catal. 2014, 312, 108–122. https://doi.org/10.1016/j.jcat.2014.01.013.

Author: Andrew Wong

1. Abstract

The purpose of this post is to study the dissociation of NO2 to NO and O on Fe (110) surface by performing a transition state search using the nudged elastic band method. A reaction energy diagram was constructed to evaluate the activation barrier for this particular reaction. The properties and geometry of the transition state were then studied by comparing it to the proposed initial and final states.

2. Introduction

NO2 reduction is an important area of research as its reaction path leads to the production of several industrially vital products, such as NH3 and N2O [1]. The reduction of NO2 into NH3 is important for environmental applications, such as water remediation projects. Additionally, research on Fe as a catalyst in NO2 reduction is of particular interest, especially in the applications of electrocatalysis in environmental systems. [2]. Using the nudged elastic band method, the purpose of this post is to utilize density functional theory calculations (DFT) to determine the transition state to reduce NO2, nitrogen dioxide, to NO, nitric oxide, and O, atomic oxygen, on an Fe (110) surface. Fe as a catalyst was selected for its catalytical desirable properties and is the most optimal pure candidate for the Haber-Bosch process of reducing N2 to NH3 [3].

3. Methodology

3.1 Calculation Parameters

A plane wave basis set with pseudo potentials in the Vienna Ab initio Simulation Package (VASP) was used to perform all DFT calculations. The Perdew, Burke, and Ernzerhof functional (PBE) was selected to model the electron-electron exchange and correlation energies. The projector augmented-wave (PAW) method was implemented to represent the ion-core electron interactions [5]. The Monkhorst-Pack method 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) was chosen with an electronic convergence tolerance of 2E-06 eV in all calculations. For all calculations, a cutoff energy of 450 eV was selected as this was seen to be optimal energy cutoff value during the geometry optimization calculations of the molecules and the Fe surface.

3.2 Model Parameter

From previous post, a 4 layer 3 x 3 Fe (110) surface was used and optimized as it is the most stable termination facet of Fe [2]. The slab was modeled with a vacuum region of approximately 10 Å with a k-point mesh of 5 x 5 x 1 and cutoff energy of 450 eV. To imitate the bulk Fe, the bottom two layers of the slab were fixed. Lastly, dipole corrections and selective dynamics were used.

3.3 Transition State Search Algorithm

Transition State searches in this post employ the Climbing Image Nudged Elastic Band Method (NEB) to determine the transition state [4]. Geometry optimizations of the initial and final states were performed at a k-point mesh of 5 x 5 x 1 and cutoff energy of 450 eV. The reactant and product states were defined as NO2 to NO and O. Figure 1 below displays the optimized initial, transition, and final states.

Based on the optimized structures of the reactant and product states, eight linearly interpolated images were created to determine the transition state. Once the tangential forces of the highest energy image is less than the absolute value of 0.05 eV/Å, the transition state search is complete. An animation of this particular transition state search is shown below.

Movie 1: Transition state search of reducing NO2 to NO and O using eight linearly interpolated images from NEB Method. Red: O, Blue: N, Purple: Fe

4. Results and Discussion

Once the transition state search was complete, a reaction energy diagram was constructed to determine the transition state and the activation barrier of this reaction. Using the eight linearly interpolated images, energies were plotted below in figure 2 relative to the reactant state.

From the transition state search, This transition state was confirmed as it had a single imaginary zero-point vibrational frequency value. From figure 2, the activation barrier was determined to be a value of 0.72 eV. Additionally, this reaction is energetically downhill favorable as the difference between the product and reactant state is approximately 2 eV.

During the trial transition state search, It is important to note that when using four linearly interpolated images, the NEB method missed a saddle point between the initial image and the first interpolated image. This was confirmed as the first image did not have an imaginary vibration frequency value during zero-point vibrational energy calculations. As a result, another transition state search was performed by creating an additional four images by linearly interpolating the original initial state and the first interpolated image. A reaction diagram including the 8 linearly interpolated images is shown below in figure 2.

As seen in figure 2, the transition state was concluded to be more similar the initial state than final state. Figure 3 was constructed to compare the bond lengths between the reactant and transition state.

As seen in figure 3, the NO bond length between oxygen atom 1 in the reactant and transition state is approximately the same. However, the NO bond of oxygen atom 1 rotated out of the plane approximately 90 degrees in the transition state, primarily contributing to the increase of 0.72 eV in energy. Additionally, the N atom began to adsorb between the adjacent Fe atom. Both of these observations most likely occurred to allow the oxygen atom 1 to rotate above the nitrogen atom like in the final proposed state. For the NO bond interaction between oxygen atom 2, elongation of the NO bond and stronger adsorption towards the Fe surface occurred. Specifically, the bond length increased by 0.2 Å to allow this atom to dissociate from the NO2 and adsorb more strongly to the Fe surface atoms.

5. Conclusion

Using the NEB method, a transition state with an activation barrier of 0.72 eV was found for the reduction of NO2 to NO and O. Although a particular transition state for the NO2 reduction was found, it is important to note that the NEB method are local minimization calculations. As local minimization calculations, it cannot confirm whether other transition states that are between the initial and final states also exist [7]. Particular caution should be considered when using the NEB method around the saddle point as using not a sufficient number of images can neglect a transition state during the transition state search. Based on a particular DFT adsorption study, the reaction barrier differed by approximately 0.1 eV [8]. Improvements on this transition search can be made by allowing the NO to adsorb onto the hollow site in the product state, the most stable site, rather than the short bridge site. Overall, the dissociation of NO2 onto the Fe (110) is an energetically favorable reaction with an activation barrier of 0.72 eV.

 

6. Citations

[1] Liu, Jin-Xun, et al. “Activity and Selectivity Trends in Electrocatalytic Nitrate Reduction on Transition Metals.” ACS Catalysis, vol. 9, no. 8, 2019, pp. 7052–7064., doi:10.1021/acscatal.9b02179.

[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.

DOI: 10.1021/acscentsci.6b00272

[4] G. Henkelman, H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, The Journal of Chemical Physics, 113 (2000) 9978-9985.

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

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

[7] Scholl, David S. and Steckel, Janice A. “Density Functional Theory: A Practical Introduction” Wiley (2009).

[8] Xu, Lang, et al. “Atomic and Molecular Adsorption on Fe(110).” Surface Science, North-Holland, 12 Sept. 2017, www.sciencedirect.com/science/article/pii/S0039602817305988?via=ihub.

by Yingdong Guan

Abstract

In this post, we studied band structures and DOS for Cu and CuO using DFT calculation. Two exchange-correlation functionals (GGA-PBE and LDA-CA-PZ) were used for comparison. Calculation results show that Cu is metal as expected. A gap of 1 eV was opened in the CuO band structure, which is experimentally observed values. The error in bandgap is a result of our naïve simulation of the complex magnetic order in CuO.

Introduction

Multiferroics are single-phase crystals that exhibit both ferromagnetism (or antiferromagnetism) and ferroelectricity (or antiferroelectricity) under the same external conditions like temperature and pressure.[1] Multiferroics are currently gathering more and more interest due to its role in studying coupling between spin and lattice and their great potential in electronic applications [1]. Copper oxide (CuO) was first purposed to be a multiferroic material in 2008 and more studies have been carried out recently [1]. In this post, we report the density functional theory (DFT) calculation of band structure and density of states (DOS) for both CuO and Cu. We also compared our results with the experimental and computational results and discuss the reason for the difference.

Calculation details

The structures for copper and copper oxide crystals are shown in Figure 1 and Figure 2. Copper crystal is in a face-centered cubic structure and the space group is Fm-3m. The lattice parameter is 2.561 Å [2].  For CuO, the experimentally observed crystal structure is a monoclinic structure crystal with C2/c symmetry [3]. Lattice parameters are a=4.653 Å, b=3.410 Å and c=5.108 Å [4]. Those crystal constants are experimentally observed values and are directed used in the calculation without reoptimization.

Our calculation used plane-wave bases with on the fly generated ultrasoft (OTFG-ultrasoft) pseudopotentials in CASTEP. Exchange-correlation functionals of both CA-PZ local density approximation (LDA) and PBE generalized gradient approximation (GGA) are used for calculation of DOS and band structure. The core radius for ultrasoft pseudopotential is 2.2 Bohr (~1.16 Å) for Cu and 1.1 Bohr (~0.58 Å) for O. The ultrasoft pseudopotential was generated with 3d10 4s1 valence electrons for Cu and 2s2 2p4 valence electrons for O.

Convergence tests were done for both functionals. From the convergence test, 550 eV is enough for Cu calculation and 750 eV is used for CuO calculation. DOS and band structure calculations using different k-point mesh were performed for Cu and CuO. For copper, the density of states was calculated using a k-point mesh of 15*15*15 and its band structure was calculated using 0.0184 Å^-1 spacing in k-space. For CuO, DOS was calculated using a k-point mesh of 8*12*8. A separation in k-space of 0.0245 Å^-1 was used to calculate the band structure. SCF tolerance was chosen to be 2.0*10^-6 and the maximum SCF cycle is chosen to be 100. We choose path W-L-G-X-W-K for Cu and L-M-A-G-Z-V for CuO.

CuO was reported to be an antiferromagnetic material and the Neel temperature is 231K. Since there are two layers of Cu atoms in the conventional cell, we used the initial spin of 0 to simulate the antiferromagnetic order in CuO. DFT + U method was used to add Coulomb interaction in calculations for CuO.

Results and discussion

Band structures and DOS diagrams for Cu calculated by GGA-PBE functional and LDA-CA-PZ functional are shown in Figures 3 and 4. Both band structure and DOS calculated by different functionals are the same. Our calculated result is also similar to other reported results. [5] There is no gap opening at the Fermi level, indicating Cu is a metal.

Diagrams for CuO are shown in Figures 5 and 6. There is a small gap opening above the Fermi level. The calculated band gap is 1.0 eV, which is pretty close to the experimentally observed value of 1.0 eV [3]. Some other experimental reported results vary from 1.0 eV to 1.9 eV [3]. This difference can result from the complex magnetic structure of CuO. The magnetic structure of CuO is shown in Figure 7 [1]. Our error in calculation is a result of the failure of simulating cycloidic spin arrangement in CuO. 2*2 supercell along a and c axis and two 1*2 supercells along a or c axis were also built for the calculation they all return the same result.

Conclusion

In this post, we have calculated band structures and DOS for Cu and CuO. Cu was calculated to be a metal, which is consistent with the experimental result. However, the calculation for band structure and DOS for CuO turns into a failure due to our naïve simulation of the magnetic order in CuO. More efforts are necessary to be paid to the band structure calculation for CuO.

 

References

[1] Scott, J. F. (2013). Room-temperature multiferroic magnetoelectrics. NPG Asia Materials, 5(11), e72-e72.

[2] Jain, A., Ong, S. P., Hautier, G., Chen, W., Richards, W. D., Dacek, S., … & Persson, K. A. (2013). Commentary: The Materials Project: A materials genome approach to accelerating materials innovation. Apl Materials, 1(1), 011002.

[3] Ekuma C E, Anisimov V I, Moreno J, et al. Electronic structure and spectra of CuO[J]. The European Physical Journal B, 2014, 87(1): 23.

[4] Pierre Villars (Chief Editor), PAULING FILE in: Inorganic Solid Phases, SpringerMaterials (online database), Springer, Heidelberg (ed.) SpringerMaterials

[5] Jedidi, Abdesslem & Rasul, Shahid & Masih, Dilshad & Cavallo, Luigi & Takanabe, Kazuhiro. (2015). Materials Chemistry A. Journal of Materials Chemistry A. 10.1039/C5TA05669A.

 

Lingjie Zhou

 

Abstract

In this post, the band structure is calculated for two phases of WTe2. The effect of spin-orbital coupling(SOC) on the band structure is checked and when the SOC is turned on, there is a gap opening at the band touching point.

Introduction

Tungsten ditelluride(WTe2) has two different phases. Shown in figure 1, the 1T phase has the hexagonal structure while the 1T’ phase comes from the lattice distortion from the 1T phase. Although the structures are differed only by distortion, the 1T WTe2 is a normal semimetal but the 1T’ WTe2 is predicted to have nontrivial topological properties such as quantum spin Hall effect. This material has sparked intense research as it is a two-dimensional material and can be easily exfoliated[4]. In this post, we used plane-wave basis sets with norm-conserving pseudopotentials to perform DFT methods to calculate the band structure and density of states(DOS) for the 1T and 1T’ WTe2.

Method

 

The CASTEP[1] package is used to carry out the DFT calculations. The exchange and correlation functional we used is the Perdew, Burke and Ernzerhof(PBE) functional described within the generalized gradient approximation(GGA) [2]. The ‘on the fly’ generated norm conserving pseudopotential for Te was generated with 6 electrons in the valence panel with (5s2 5p4)  and for W was generated with 14 electrons in the valence panel with (5s2 5p6 5d4 6s2) as the electronic configuration. The self-consistent-field tolerance of the calculated energy is 2^-6 eV and the maximum self-consistent-field tolerance cycle is 100. The monolayer structure is created with 15 A vacuum slab.

 

k point Convergence 

 

The convergence of k point is checked on the geometrically optimized 1T and 1T’  WTe2 monolayer structure with 15A vacuum slab. The Ecut 900eV is hypothesized as sufficient and will be further checked. For 1T phase, after the configuration goes beyond 13×13×1, the energy difference is below 0.01eV, so we would use 13×13×1 for later calculation. For 1T’ phase, after the configuration goes beyond 15×15×1, the energy difference is below 0.01eV, so we would use 15×15×1 for later calculation.

Ecut	Configuration	number of irreducible k points	energy(eV)	energy difference(eV)
900	10×10×1	14	-647.6695202	--
900	12×12×1	19	-647.6771917	-0.0076715
900	13×13×1	21	-647.6894468	-0.0122551
900	15×15×1	27	-647.6928398	-0.003393
900	20×20×1	44	-647.69014	0.0026998

Table 1.a Convergence check of k points for 1T WTe2

 

 

Ecut	configuration	number of irreducible k points	energy(eV)	energy difference(eV)
900	10×10×1	25	-1309.782212	--
900	12×12×1	36	-1309.753284	0.028928
900	15×15×1	64	-1309.769129	-0.015845
900	20×20×1	100	-1309.766136	0.002993

Table 1.b Convergence check of k points for 1T' WTe2

 

 

plane-wave basis set cutoff energy(Ecut) Convergence 

 

The convergence of Ecut is checked using k point configuration determined from the previous section(13×13×1 for 1T WTe2 and 15×15×1 for 1T’ WTe2 ). For 1T and 1T’ phase, the energy difference drops below 0.01 eV when the Ecut exceeds 900eV. So we will use 900eV as Ecut for the band calculation.

Ecut	configuration	number of irreducible k points	energy(eV)	energy difference(eV)
600	15×15×1	21	-647.6775926	--
900	15×15×1	21	-647.6894468	-0.0118542
1200	15×15×1	21	-647.6942889	-0.0048421

Table 2.a Convergence of Ecut for 1T WTe2

 

 

Ecut	configuration	number of irreducible k points	energy(eV)	energy diffference(eV)
600	15×15×1	64	-1309.741619	--
900	15×15×1	64	-1309.769129	-0.02751
1200	15×15×1	64	-1309.769712	-0.000583

Table 2.b convergence of Ecut for 1T' WTe2

 

 

Band Structure

Band structure is calculated with 0.015A^-1 separation in the k space. For 1T WTe2, we choose the path G-M-K-G while for 1T’ WTe2, we choose Y-G-Y. Only the bands close to the Fermi surface is shown. For the 1T’ phase, without SOC, the bands will closely touch each other. However, if we take SOC into account, there will be a gap opening at the touching point.

 

Fig 5 is the band structure of 1T’ WTe2 from the published papers[3]. The band structure near the Fermi surface is qualitatively similar to what we’ve calculated. The gap opening due to the SOC is also confirmed in our post.

 

Conclusion

In this post, we calculate the band structure for monolayer 1T and 1T’ WTe2. We saw that SOC will open a gap when bands touching each other. However, to further determine if there are topological properties such as band inversion, the further calculation to track how the band evolves from 1T to 1T’ is needed.

 

Reference

1   Burke, K. The ABC of DFT.  (2007).

2   Clark, S. First principles methods using CASTEP. Z. Kristallogr 220, 567-570 (2005).

3   Wu, S. et al. Observation of the quantum spin Hall effect up to 100 kelvin in a monolayer crystal. Science 359, 76-79, doi:10.1126/science.aan6003 (2018).

4   Tang, S. et al. Quantum spin Hall state in monolayer 1T’-WTe2. Nature Physics 13, 683-687, doi:10.1038/nphys4174 (2017).