Diboron-porphyrin monolayer: A new 2D semiconductor

Computational Materials Science xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Diboron-porphyrin monolayer: A new 2D semiconductor Raphael M. Tromera, , Isaac M. Felixa, Aliliane Freitasa, Sérgio Azevedob, Luiz Felipe C. Pereiraa, ⁎

a b

⁎

Departamento de Física, Universidade Federal do Rio Grande do Norte, Natal 59078-970, Brazil Departamento de Física, Universidade Federal da Paraíba, João Pessoa 58051-970, Brazil

ARTICLE INFO

ABSTRACT

Keywords: Graphene 2D material Porphyrin Porphyrin semiconductor Boltztrap Diboron porphyrin

We investigate a new 2D semiconductor material with BC10 N2 stoichiometry based on the molecular structure of diboron porphyrin. Ab-initio molecular dynamics simulations indicate mechanical stability for temperatures up to T = 4000 K. Electronic structure calculations based on density functional theory predict a 0.6 eV direct band gap. We employ the Boltzmann transport equation within the relaxation time approximation to analyze electronic transport properties as a function of charge carrier density at room temperature. The lattice thermal conductivity at room temperature was obtained from non-equilibrium molecular dynamics simulations, and we found anisotropic conductivities given by 160 and 115 W/m·K along perpendicular in-plane directions. The thermoelectric figure of merit ZT calculated at 300 K has a maximum value of 0.008. This result indicates that, in spite of its interesting anisotropic transport properties, BC10 N2 is not a suitable candidate for thermoelectric applications.

1. Introduction In 2004 graphene layers were separated from bulk graphite through a mechanical exfoliation method [1,2]. Due to its unique physical properties, graphene has received special attention regarding the construction of transistors for flexible electronic devices [3]. However, the use of graphene in nanoelectronic devices is not possible because it is a semi-metal (sometimes called a null-gap semiconductor) which prevents the device from working in an on-off current regime. Therefore, controlling the band gap is a key factor in the production of graphenebased electronic devices. Possible ways of manipulating the electronic structure of graphene in order to create a band gap have been investigated [4,5]. A viable alternative is to synthesize semiconducting 2D materials with an appropriate energy gap in their natural form [6–8]. Recently, there has been growing interest in exploiting 2D porphyrin-based polymers [9–12]. Porphyrin molecules serve as building blocks since they can bond together to form much larger systems [13–15]. In addition, porphyrin exhibits a high coordination number and can bind to several different molecular groups. Those 2D porphyrin-based systems have attracted much interest due to their possible application in hydrogen storage [16], catalysis [17], and even spintronics [18]. Another advantage of working with porphyrin-based structures is the possibility of building large systems through a polymerization process. Abel and co-authors successfully synthesized a 2D

⁎

layer of phthalocyanine containing iron by a polymerization process [19]. Phthalocyanine is just one of the possible structures that can be formed by porphyrin molecules. A porphyrin-based system of interest, which has been previously investigated, consists of a porphyrin molecule with two boron atoms in its center forming a B2 bond. From a purely chemical point of view, this structure is not favourable since the B2 bond has a large formation energy. Nonetheless, this structure presents interesting optical properties depending on what binds to the boron atoms, such as fluorine, chlorine, oxygen, and bromide among other possibilities [20]. Arnold et al. synthesized a porphyrin structure with Li atoms, opening up a possible route to obtain porphyrin with a B2 bond [21,22]. In fact, this Li structure can interact with borane halides to produce porphyrin with B2 bonds where each boron links to two nitrogen atoms. Therefore, this chemical route produces diboron molecules, also known as dyboranylporphyrins, which were synthesyzed by Weiss et al. [23,24]. In this work we propose a new 2D semiconductor material inspired by one of the porphyrin-based structures investigated by P.J. Brothers et al. [20]. We consider a unit cell consisting of a diboron-porphyrin molecule where the C–H bonds are replaced by C–C covalent bonds, building a 2D supercell with BC10 N2 stoichometry. Through extensive computational modelling we investigate the stability and physical properties of 2D diboron-porphyrin (2DDP). First principles calculations are used to verify the mechanical stability of 2DDP at 0 K as well as finite temperatures up to 4000 K. Electronic structure calculations

Corresponding authors. E-mail addresses: [email protected] (R.M. Tromer), [email protected] (L.F.C. Pereira).

https://doi.org/10.1016/j.commatsci.2019.109338 Received 20 June 2019; Received in revised form 3 October 2019; Accepted 5 October 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Raphael M. Tromer, et al., Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109338

Computational Materials Science xxx (xxxx) xxxx

R.M. Tromer, et al.

predict the 2DDP is a direct gap semiconductor with an energy gap of 0.6 eV. The Boltzmann transport equation is employed to investigate electronic transport properties of 2DDP. Meanwhile, its lattice thermal conductivity is predicted by non-equilibrium molecular dynamics simulations. We find that electronic and thermal transport coefficients are anisotropic along the in-plane direction of 2DDP. Finally, the thermoelectric efficiency of 2DDP is calculated, and we obtain a maximum figure of merit of 0.008 at 300 K. Therefore, in spite of its promising physical properties, 2DDP does not seem to be viable for novel thermoelectric devices. 2. Methods

Fig. 1. Atomic structure of 2D diboron-porphyrin. The orthorrombic unit cell with 26 atoms is highlighted.

We performed ab initio calculations within two different density functional theory approaches. Localized numerical orbitals as implemented in SIESTA [25,26] were used to investigate the mechanical stability of the structures via its phonon dispersion at T = 0 K, and the dynamical stability from 1000 to 4000 K via ab initio molecular dynamics simulations. Projector-augmented plane waves [27] were employed in the calculation of electronic properties using Quantum Espresso (QE) [28]. It is worth mentioning that phonon dispersion could also be obtained from QE via density functional perturbation theory (DFPT). However, since we are dealing with a 26-atom unitcell, this calculation method becomes too computationally expensive. Therefore, we choose to work with the finite difference method implemented in SIESTA [29]. In SIESTA we use the double zeta polarized basis set (DZP), with core eletrons described by norm-conserving Troullier-Martins pseudopotentials [30]. The exchange-correlation energy is expressed in the form of Perdew-Burke-Ernzerhof (PBE) within the generalized gradient approximation (GGA) [31]. The reciprocal space is sampled on a Monkhorst-Pack scheme with a 10 × 10 × 1 k-point mesh, and an energy cutoff of 2720 eV. Similarly, in QE the exchange-correlation energy is also expressed with PBE-GGA. The energy cut-off for the electronic wave function and the charge density were chosen as 680 eV and 6800 eV, respectively. Integration over the first Brillouin zone also considered a 10 × 10 × 1 k-point mesh. In both SIESTA and QE structural relaxation was performed with a convergence criteria of 0.05 eV/ Åas the maximum force on one atom. Both codes provide very similar structures after optimization, with minor difference in atomic distances. In order to obtain the thermoelectric coefficients we use BoltzTraP, which solves the Boltzmann Transport Equation (BTE) within the relaxation time approximation (RTA) [32]. BoltzTraP takes as input the band structure calculations performed by QE, and provides the electrical conductivity, the Seebeck coefficient and the electronic contribution to the thermal conductivity. In this case, a 50 × 50 × 1 k-point mesh was used, so that all quantities extracted from BoltzTraP are converged. The methodology implemented in BoltzTraP does not provide information regarding the lattice thermal conductivity, i.e. the lattice contribution to the heat transport. Therefore, we obtain the phonon thermal conductivity via non-equilibrium molecular dynamics simulations in LAMMPS [33] with the method proposed by MüllerPlathe [34]. In this case, the atomic interaction between B, N and C was described by the Tersoff force field with a recent parameterization adapted for 2D nanostructures [35,36].

supercell in Fig. 1 contains rings with five atoms (4 carbon and 1 nitrogen), 6 atoms (all carbon), 6 atoms (3 carbon, 2 nitrogen and 1 boron), 7 atoms (3 carbon, 2 nitrogen and 2 boron) and 8 atoms (all carbon). The formation energy of 2DDP is calculated at 8.7 eV/atom, which is comparable to that of graphene at 8.8 eV/atom. Besides graphene, the formation energy of 2DDP is also comparable to other carbon allotropes and carbon-based structures [37,38]. It is important to point out that 2DDP was not built on the basis of a particular structure in which we attempted different atomic configurations to investigate stability. Instead, we rely on a single previously synthesized molecule and replace the C–H bonds by C–C covalent bonds. This well established experimental technique is very common in porphyrin systems, where large molecular blocks, such as polymers, are built on the basis of the individual molecule. Even so, we performed some stability tests maintaining the shape of porphyrin and replacing B–B covalent bonds by B–C, B–N, C–C, N–C and N–N. Table 1 presents the per atom formation energy of alternative structures to 2DDP, in which bonds other than the B–B bond have been considered. According to the calculated values, the formation energy of the alternative structures is at most equal to the one for 2DDP (B–C and C–C cases). We take this observation as an indication that 2DDP is among the most stable BNC structures obtained from porphyrin. Therefore, we believe this to be an indication that synthesizing 2DDP could indeed be energetically favourable. In order to verify if the predicted structure is mechanically stable at T = 0 K we calculate the phonon dispersion with SIESTA. Displacements of 0.04 Åwere used to obtain the inter-atomic force constants. Fig. 2 shows the phonon dispersion along the reciprocal = (0, 0, 0), X = (0.5, 0, 0), U = (0.5, 0.5, 0.0), Y = space path: (0.0, 0.5, 0.0) and back to . The absence of imaginary frequencies (which would appear as negative values in Fig. 2) is an indication that the predicted structure should be stable at T = 0 K. Around there are three acoustic modes, two of them representing transversal vibrational modes and a third one for the longitudinal mode. The longitudinal mode and one of the transverse modes show the usual linear dispersion relation, while the out-of-plane transverse mode displays a non-linear behavior, characteristic of 2D materials. The lowest optical modes appear close to 200 cm 1. The absence of imaginary frequencies is not enough to guarantee the dynamical stability of 2DDP at finite temperatures. Therefore, we performed ab initio molecular dynamics simulations with SIESTA in the NPT ensemble at T = 1000, 2000, 3000 , and 4000 K. The simulations were performed with a 2 × 2 × 1 supercell, obtained by replication of

3. Results 3.1. Mechanical stability

Table 1 Formation energy per atom for 2DDP and alternative BNC structures in which bonds other than the B–B bond have been considered.

The unit cell of 2DDP is composed of 26 atoms, where there are 20 carbon atoms, 4 nitrogen atoms, and 2 boron atoms. Fig. 1 shows the optimized structure predicted by QE. The lattice vectors are given by a1 = 8.16 Å x and a2 = 8.64 Å y . For comparison, the lattice vectors predicted by SIESTA are a1 = 8.19 Å x and a2 = 8.68 Å y . The 2DDP 2

Structure

2DDP

B–C

B–N

C–C

N–C

N–N

Eformation (eV/atom)

−8.7

−8.7

−8.6

−8.7

−8.6

−8.4

Computational Materials Science xxx (xxxx) xxxx

R.M. Tromer, et al.

Fig. 2. Phonon dispersion for the 26-atom orthorrombic unit cell. The absence of negative frequencies indicate structural stability at T = 0 K.

the 26-atom unit cell. A time step of 0.1 fs was used to integrate the equations of motion for a total simulation time of 2 ps. Fig. 3 shows snapshots of the final configuration at each specified temperature. The crystalline structure of 2DDP is preserved at temperatures as high as 4000 K, which would make the material suitable for high-temperature applications. So far we verified the mechanical stability of 2DDP not only at T = 0 K, but also for high temperatures. Next, we analyze its electronic structure.

Fig. 4. (a) Electronic band structure, and (b) projected density of states. According to PBE calculations, 2DDP presents a 0.6 eV direct band gap.

direct gap semiconductor. It is important to keep in mind that PBE underestimates the energy gap, and a more precise calculation would probably result in a larger band gap [39]. In fact, we have estimated the gap with the HSE06 functional and obtained a value of 0.9 eV [40,41]. In Fig. 4(b) we present the projected density of states, where it can be seen that the largest contribution to band formation comes from the carbon p orbital (red line). In fact, the top of the valence band and the bottom of the conduction band are mostly due to this orbital. We have also considered the effect of mechanical strain on the electronic properties of 2DDP. To this end, uniaxial tensile strains were applied along the x and y directions, ranging from 2% up to 10%. The top panels of Fig. 5 show the electronic band structure of 2DDP under a 6%

3.2. Electronic structure Our electronic structure analysis was performed with QE since its plane-wave description is more conveniently integrated with BoltzTraP. Fig. 4(a) shows the DFT band structure calculation along the k-point X U Y path . There is a direct band gap at the U point with an energy of 0.6 eV. Therefore, this novel structure would be a

Fig. 3. Final configuration of ab initio molecular dynamics simulations at high temperatures. The crystalline structure is preserved at temperatures as high as 4000 K.

Fig. 5. Electronic band structure under a 6% strain along (a) x - and (b) y -directions respectively. (c) Band gap energy as a function of applied strain along each in-plane direction. 3

Computational Materials Science xxx (xxxx) xxxx

R.M. Tromer, et al.

Fig. 6. Localized density of states for the top of the valence band and the bottom of the conduction band.

uniaxial strain along x and y directions. Meanwhile Fig. 5(c) shows the estimated DFT band gap as a function of tensile strain. We have verified that the band gap is always direct within the strain range considered. Interestingly, the magnitude of the energy band gap may increase or decrease, depending on the direction in which the strain is applied. In fact, this dependence is due to the geometrical arrangement of the atomic orbitals in 2DDP. The orbitals with greatest contribution to the band structure of 2DDP become closer when tensile strain is applied along the x direction, approximating the valence and conduction bands. Meanwhile, the opposite happens in the case of strain along y. From Fig. 5(c) we notice that the band gap increase saturates close to 1 eV in the case of strain along y, while it seems to decrease linearly in the case of strain along x. This behavior offers the possibility of mechanically controlling the band gap of 2DDP, similarly to a recent proposal for graphyne heterostructures [42]. Localized density of states for the top of the valence band and the bottom of the conduction band are shown in Fig. 6. We observe the formation of delocalized charges on the carbon atoms due to bonds formed by pz atomic orbitals. This is a typical behavior of 2D system containing carbon atoms in hybridization states [43]. The data in Figs. 4(b) and 6 are in good agreement, since both indicate that carbon’s p atomic orbitals have a greater contribution to valence and conduction states. Our electronic structure calculations predict that, in the absence of tensile strain, 2DDP is a direct band gap semiconductor with an energy gap of at least 0.6 eV. Therefore, it is worth investigating its thermoelectric properties. Next, we feed information about the electronic structure to BoltzTraP in order to predict the thermoelectric properties of 2DDP.

Fig. 7. (a) and (b) Electrical conductivity as a function of carrier concentration for electrons and holes. (c) and (d) Power factor as a function of electron and hole concentration. (e) and (f) Electronic contribution to the thermal conductivity.

pronounced in the high doping region. The so-called figure of merit ZT is a good predictor of the thermoelectric efficiency of a material. This quantity is defined as

ZT =

(

ele

S2 +

ph )

T

(1)

where ele and ph are the electronic and the lattice (phonon) contributions to the thermal conductivity, respectively. The product S 2 defines the Power Factor (PF) of the material. In order to maximize ZT, one can either increase PF or decrease the total thermal conductivity. Fig. 7(c) and (d) shows the dependence of the PF on the carrier concentration for both types of doping. The highest PF value at 300 K, 2.8 mW/m·K2, is obtained for n-type doping with a carrier concentration of 28 × 1013 electrons/cm2. Meanwhile, the maximum PF for p-type doping does not exceed 1 mW/m·K2. According to BoltzTraP, ele is proportional to , and takes values lower than 10 W/m·K at 300 K, as shown in Fig. 7(e) and (f). Similar 2D materials present ph on the order of 10 2 W/m·K, so we can already speculate that the thermal conductivity is dominated by the lattice contribution [47,7,48]. Non-equilibrium molecular dynamics simulations predict an anisotropic lattice thermal conductivity for 2DDP samples from 20 to 300 nm, as shown in Fig. 8. The size-independent lattice thermal conductivity of 2DDP along its in-plane directions can be obtained by fitting the data according to [49]

3.3. Thermoelectric properties BoltzTraP relies on the RTA to calculate the electrical conductivity , the Seebeck coefficient S, and the electronic contribution to the thermal conductivity ele . In its formalism, and ele are given in terms of a relaxation time , such that their units are not ( ·m) 1 and W/m·K but actually ( ·m·s)−1 and W/m·K·s respectively. A precise calculation of is an arduous task involving several calculation steps [44]. In what follows we consider an ansatz where the relaxation time will take the value recently obtained for a similar 2D carbon-based semiconductor [45,44,46]. This system presents a C2 N stoichiometry and its formation energy is close to that of 2DDP. We expect the ansatz to be mostly valid when the charge carrier density is high, and take the relaxation time to be 23 fs at 300 K. Fig. 7 (a) and (b) present the electrical conductivity as a function of carrier density for electrons and holes, respectively. In the case of ntype doping, the electrical conductivity is almost isotropic along x and y directions for any concentration. Meanwhile, for p-type doping, the electrical conductivity is anisotropic, and the difference is more

1 ph (L)

=

1 ph

1+

eff

L

,

(2)

where eff is an effective phonon mean free path and ph is the intrinsic x · lattice thermal conductivity. We find ph = 160 W/m K and y y x · = 42.8 = 115 = 59.1 nm, and W/m K and nm. Interestingly, eff eff ph the lattice thermal conductivity is higher along the x direction, while the electrical conductivity is higher along y. This indicates that the electronic and phononic transport processes are decoupled in 2DDP. Once we have obtained the PF and the total thermal conductivity of 2DDP, we can calculate its ZT via Eq. (1). Fig. 9 shows the dependence of ZT on the concentration of charge carriers. Due to the anisotropy in 4

Computational Materials Science xxx (xxxx) xxxx

R.M. Tromer, et al.

simulations that the proposed material in mechanically stable for temperatures up to T = 4000 K. Electronic structure calculations based on the PBE functional predict a 0.6 eV direct band gap, while calculations based on the HSE06 functional predict a larger gap of 0.9 eV. The LDOS at the top of the valence band and the bottom of the conduction band show the formation of bonds between carbon pz orbitals. This sp2 hybridization provides a robust mechanical stability for 2DDP comparable to that of graphene. The Boltzmann transport equation within the relaxation time approximation has been used to analyze electronic transport properties as a function of charge carrier density at room temperature. We found anisotropic electronic transport coefficients along in-plane directions of 2DDP. Non-equilibrium molecular dynamics simulations yield anisotropic lattice thermal conductivities along perpendicular in-plane directions. The thermoelectric figure of merit ZT calculated at 300 K has a maximum value of 0.008 for n-type doping and transport along the y direction. The result indicates that, in spite of its interesting anisotropic transport properties, BC10 N2 is not a suitable candidate for thermoelectric applications.

Fig. 8. Lattice thermal conductivity of 2DDP as a function of sample length. Notice the anisotropy in ph along the in-plane directions.

CRediT authorship contribution statement Raphael M. Tromer: Conceptualization, Formal analysis, Investigation, Writing - original draft. Isaac M. Felix: Formal analysis, Investigation. Aliliane Freitas: Formal analysis, Investigation. Sérgio Azevedo: Conceptualization, Supervision, Writing - review & editing. Luiz Felipe C. Pereira: Conceptualization, Supervision, Writing - original draft, Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors thank Leonardo D. Machado for a critical reading of the manuscript. We acknowledge financial support from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for the project “Thermal and electronic transport in 2D materials” (Grant 309961/2017). We are grateful for the provision of computational resources by the High Performance Computing Center (NPAD) and the CLIMA cluster at UFRN.

Fig. 9. Figure of merit ZT as a function of carrier density at 300 K considering ntype (a) and p-type (b) doping.

, ele and ph , ZT also depends on the transport direction. In the case of n-type doping, we predict a maximum ZT of 0.004 for transport along the x direction and 0.008 when the transport is along the y direction. Both maxima happen at a charge carrier concentration of 2.8 × 1014 electrons/cm2. The charge carrier concentration for maximum ZT in 2DDP is in accordance with what was recently predicted for nitrogenated holey graphene [46]. Meanwhile, for p-type doping ZT does not exceed 0.003 at any carrier concentration, as shown in Fig. 9(b). It is clear that ZT values for 2DDP are low when compared to standard thermoelectrics such as Bi2 Te3 [50], as is the case for many carbon-based materials. Furthermore, our results indicate that despite its interesting anisotropic transport properties, 2DDP is not a suitable candidate for thermoelectric applications. Nonetheless, 2DDP still presents one possible advantage since its structure presents a high coordination number and is able to bind with many molecular groups [20].

References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666. [2] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.V. Morozov, A.K. Geim, Proc. Natl. Acad. Sci. 102 (2005) 10451. [3] H. Gwon, H.-S. Kim, K.U. Lee, D.-H. Seo, Y.C. Park, Y.-S. Lee, B.T. Ahn, K. Kang, Energy Environ. Sci. 4 (2011) 1277. [4] F. Xia, D.B. Farmer, Y.-M. Lin, P. Avouris, Nano Lett. 10 (2010) 715. [5] E. Cadelano, P.L. Palla, S. Giordano, L. Colombo, Phys. Rev. B 82 (2010) 235414 . [6] J. Mahmood, E.K. Lee, M. Jung, D. Shin, I.-Y. Jeon, S.-M. Jung, H.-J. Choi, J.M. Seo, S.-Y. Bae, S.-D. Sohn, et al., Nat. Commun. 6 (2015) 6486. [7] B. Mortazavi, O. Rahaman, T. Rabczuk, L.F.C. Pereira, Carbon 106 (2016) 1. [8] B. Mortazavi, M. Shahrokhi, M. Raeisi, X. Zhuang, L.F.C. Pereira, T. Rabczuk, Carbon 149 (2019) 733. [9] C. Zou, C.-D. Wu, Dalton Trans. 41 (2012) 3879. [10] S. Hamad, N.C. Hernandez, A. Aziz, A.R. Ruiz-Salvador, S. Calero, R. Grau-Crespo, J. Mater. Chem. A 3 (2015) 23458. [11] J. Tan, W. Li, X. He, M. Zhao, RSC Adv. 3 (2013) 7016. [12] H.K. Singh, P. Kumar, U.V. Waghmare, J. Phys. Chem. C 119 (2015) 25657, https:// doi.org/10.1021/acs.jpcc.5b09763. [13] J. Tan, W. Li, X. He, M. Zhao, RSCAdv. 3 (2013) 7016. [14] L.E. Jiménez-Ramírez, D.C. Camacho-Mojica, E. Muñoz-Sandoval, F. López-Urías, Carbon 116 (2017) 381. [15] Y. He, M. Garnica, F. Bischoff, J. Ducke, M.-L. Bocquet, M. Batzill, W. Auwärter, J.V. Barth, Nat. Chem. 9 (2017) 33.

4. Conclusion In summary, we proposed a new 2D semiconductor material with BC10 N2 stoichiometry based on the molecular structure of diboron porphyrin. We have shown via ab initio molecular dynamics 5

Computational Materials Science xxx (xxxx) xxxx

R.M. Tromer, et al. [16] D.H. Lee, W.J. Lee, W.J. Lee, S.O. Kim, Y.-H. Kim, Phys. Rev. Lett. 106 (2011) 175502 . [17] M.E. Ali, B. Sanyal, P.M. Oppeneer, J. Phys. Chem. C 113 (2009) 14381. [18] W.J. Cho, Y. Cho, S.K. Min, W.Y. Kim, K.S. Kim, J. Am. Chem. Soc. 133 (2011) 9364. [19] M. Abel, S. Clair, O. Ourdjini, M. Mossoyan, L. Porte, J. Am. Chem. Soc. 133 (2010) 1203. [20] P.J. Brothers, Inorgan. Chem. 50 (2011) 12374. [21] J. Arnold, J. Chem. Soc. Chem. Commun. 976 (1990). [22] H. Brand, J.A. Capriotti, J. Arnold, Inorg. Chem. 33 (1994) 4334. [23] A. Weiss, M. Hodgson, P. Boyd, W. Siebert, P. Brothers, Chem. Eur. J. 13 (2007) 5982. [24] W.J. Belcher, M.C. Hodgson, K. Sumida, A. Torvisco, K. Ruhlandt-Senge, D.C. Ware, P.D.W. Boyd, P.J. Brothers, Dalton Trans. 1602 (2008). [25] D. Sánchez-Portal, P. Ordejón, E. Artacho, J.M. Soler, Int. J. Quant. Chem. 65 (1997) 453. [26] J.M. Soler, E. Artacho,J.D. Gale, A. García, J. Junquera, P. Ordejón, D. SánchezPortal, 14 (2002) 2745. [27] P.E. Blöchl, Phys. Rev. B 50 (1994) 17953. [28] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, et al., J. Phys.: Condensed Matter 21 (2009) 395502 . [29] Y. Jaluria, S.N. Atluri, Comput. Mech. 14 (1994) 385. [30] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. [31] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [32] G.K. Madsen, D.J. Singh, Comput. Phys. Commun. 175 (2006) 67. [33] S. Plimpton, J. Comput. Phys. 117 (1995) 1.

[34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50]

6

F. Müller-Plathe, J. Chem. Phys. 106 (1997) 6082. J. Tersoff, Phys. Rev. B 37 (1988) 6991. A. Kínací, J.B. Haskins, C. Sevik, T. Çağín, Phys. Rev. B 86 (2012) 115410 . Z. Wang, X.-F. Zhou, X. Zhang, Q. Zhu, H. Dong, M. Zhao, A.R. Oganov, Nano Lett. 15 (2015) 6182. A. Freitas, L.D. Machado, C.G. Bezerra, R.M. Tromer, L.F.C. Pereira, S. Azevedo, RSC Adv. 8 (2018) 24847. K.A. Johnson, N.W. Ashcroft, Phys. Rev. B 58 (1998) 15548. J. Heyd, G.E. Scuseria, J. Chem. Phys. 118 (2003) 8. A.F.I. Aliaksandr, V. Krukau, Oleg A. Vydrov, G.E. Scuseria, J. Chem. Phys. 125 (2006) 224106. J. Yun, Y. Zhang, Y. Ren, M. Xu, J. Yan, W. Zhao, Z. Zhang, Phys. Chem. Chem. Phys. 20 (2018) 26934. D.R. Cooper, B. D’Anjou, N. Ghattamaneni, B. Harack, M. Hilke, A. Horth, N. Majlis, M. Massicotte, L. Vandsburger, E. Whiteway, et al., ISRN Condensed Matter Phys. 2012 (2012), https://doi.org/10.5402/2012/501686. L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V.P. Dravid, M.G. Kanatzidis, Nature 508 (2014) 373. R.M. Tromer, M.G.E. da Luz, M.S. Ferreira, L.F.C. Pereira, J. Phys. Chem. C 121 (2017) 3055. Y. Zhao, Z. Dai, C. Lian, S. Meng, RSC Adv. 7 (2017) 25803. L.F.C. Pereira, B. Mortazavi, M. Makaremi, T. Rabczuk, RSC Adv. 6 (2016) 57773. I.M. Felix, L.F.C. Pereira, Scientific Rep. 8 (2018) 2737. P.K. Schelling, S.R. Phillpot, P. Keblinski, Phys. Rev. B 65 (2002) 144306 . J. Zhang, H.J. Liu, L. Cheng, J. Wei, J. Shi, X.F. Tang, C. Uher, J. Appl. Phys. 116 (2014) 023706 .