Author Archives: Vishal Jindal

1-D band structure of polythiophene using different functionals

1 Reply

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.

Unit PT Cell

Fig 1. All-trans polythiophene chain with two rings in each unit cell having dimensions 8 x 15 x 15 (in angstroms), here two periodic unit cells are shown.

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

12

21

Fig 2. Band structures for polythiophene using LDA, GGA-PBE, GGA-PW91 and GGA-BLYP functionals.

Density of State

Upd1

upd2

Fig 3. The density of state plot for different XC functionals

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:

eff mass e

The denominator of the above equation is calculated by fitting a parabola (y = kx2 + 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-PZ1.471.231.731.461.39
GGA PBE1.491.251.721.471.39
GGA PW911.461.251.721.491.38
GGA BLYP1.461.251.721.491.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.

DFT Frequency Calculations of Thiophene using VASP

1 Reply

by Vishal Jindal

1. Introduction

Thiophene is a heterocyclic aromatic compound with the formula C4H4S, 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.

Model TP

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

Thio MS Crystal

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.

Thiophene Optimized

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 (C4H4S) 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.

Gr1

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

gr2

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 Eo  = – 54.237 eV:

Equation 5.12

 

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

Displacement value used for Hessian Matrix (δb)Inverse of displacement value
(1/δb)		Zero-point Energy (ZPE) Correction (eV)ZPE Corrected Energy (eV) (Eo = -54.237 eV)Absolute sum of Rotational/ Translational modes (cm-1)
0.0052001.783-52.454518.82
0.0061671.775-52.462467.79
0.011001.771-52.466358.96
0.025401.769-52.468223.58
0.05201.770-52.467135.83
0.07141.772-52.465131.04
0.1101.778-52.459173.87

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.

Determination of preferred crystal structure of Platinum and comparison with experimental observation.

Leave a reply

by Vishal Jindal

1. Introduction

Platinum (symbol Pt and atomic number 78) is a noble metal. It is used in catalytic converters, electrodes, and many other laboratory equipments, dentistry and jewelry [6]. Experimentally, platinum (Pt) is found to exist in cubic close-packed (ccp) structure with a lattice constant of a = 3.9242 Å [1]. In this post, we are using Density Functional Theory [DFT] [2] calculations to determine the preferred crystal structure for platinum metal amongst simple cubic (sc), cubic close-packed (ccp), also known as face-centered cubic (fcc) and hexagonal close packing (hcp) crystal structures.

2. Methodology

We used a plane-wave basis set with ultrasoft pseudopotentials as implemented in the CASTEP code [3] to perform DFT calculations to get total energy for different crystal structures and lattice parameters. Materials Studio was used as a builder, visualizer, and user interface for the CASTEP calculations. The exchange and correlation energies were calculated using the Perdew, Burke, and Ernzerhof (PBE) [4] functional described within the generalized gradient approximation (GGA) [5]. SCF tolerance is taken to be 2.0e-6 eV/atom for all the calculations. The “on the fly” generated ultrasoft pseudopotential for Pt has a core radius of 2.4 Bohr (1.27 Angstroms) and was generated with 32 electrons in the valence panel with (4f14 5s2 5p6 5d9 6s1) as the electronic configuration.

3. K-Points and Cut-off Energy Optimization

Before varying the lattice parameter to find the most stable crystal structure, we need to make sure that the plane-wave basis sets give convergent results for energies with respect to mesh size/ k-points and cut-off energy (ENCUT). To perform k-points and ENCUT optimization we used the experimental value of lattice parameter, i.e. a = 3.9242 Å.

3.1 K – Point Optimization

3.1.1 Simple Cubic (sc)

To optimize the k-points, we calculated total energy of simple cubic (sc) lattice for different mesh sizes (keeping ENCUT fixed at material studio default i.e. 272.1 eV). Fig. 1 shows the variation of total energy when we change the mesh size. We choose 8 x 8 x 8 lattice having 20 k-points as total energy changed less than 0.001 eV as we further increased the mesh size.

SC Vishal

Figure 1. k-point optimization for simple cubic crystal lattice

3.1.2 Face-centered cubic (fcc)

FCC Vishal

Figure 2. k-point optimization for face-centered cubic crystal lattice

3.1.3 Hexagonal Close Packed (hcp)

HCP k vishal

Figure 3. k-point optimization for hexagonal close-packed crystal lattice

3.2 Energy cut-off

Now, using the optimized number of k-points from simple cubic, we varied ENCUT from 250-450 eV. Again, using the same convergence criterion of |ΔE| less than 0.001 eV on increasing the cut-off energy. In Fig. 4, we found that the energy cut-off of 425 eV sufficiently converged the total energy of the system. Similarly for FCC and HCP, 425 eV was taken to be the cut-off energy after checking the convergence of total energy at 400 eV, 425 eV, and 450 eV.

Energy Cut off

Figure 4. Total energy vs energy cut-off [ENCUT]. Total energy tends to converge at 425 eV.

4. Results

To optimize the lattice parameters, we minimize the total energy per atom (eV/atom) with respect to lattice constant (a) using respective k-points and energy cut-off for SC and FCC crystal lattice.

4.1 Simple Cubic (SC)

SC a

Figure 5. Total energy vs lattice parameter plot for SC

4.2 Face Centered Cubic (FCC)

VJ FInal

Figure 6. Total energy vs lattice parameter plot for FCC

4.3 Hexagonal Closed Pack (HCP)

VJ HCP Final

Figure 7. Total energy vs lattice parameter (a) plot for different values of c/a ratio along with trendlines connecting the same c/a ratio (c/a = 1.5, 1.67, 1.85)

5. Conclusion

The table below summarizes the optimized lattice parameter for all three types of crystal lattices studied in our report.

Lattice System Lattice Parameter/s
(Angstroms)		Total Energy per atom
(eV/atom)		k-points
(M x M x N)		Energy Cut-off
(eV)
Simple Cubic (SC)a = 2.66-13049.8920 (8 x 8 x 8)425
Face Centered Cubic (FCC)a = 3.97-13050.96110 (10 x 10 x 10)425
Hexagonal Close Packed (HCP)a = 2.8
c = 4.2		-13050.9076 (12 x 12 x 8)425

We can observe that among the three types of lattice structure listed above, the minimum total energy per atom is found to be for the Face Centered Cubic (FCC) crystal structure. Therefore, our DFT calculations show that the Platinum(Pt) crystal is most stable in the FCC lattice structure with a lattice parameter (a) of 3.97 Å. This is in line with the experimentally observed Face Centered Cubic (FCC) crystal structure of Platinum with a = 3.9272 Å.

 

6. References

[1] https://www.webelements.com/platinum/crystal_structure.html

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

[3] 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)

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

[6] https://en.wikipedia.org/wiki/Platinum