Methodological Issues in First-Principle Calculations of CH₃NH₃PbI₃ Perovskite Surfaces: Quantum Confinement and Thermal Motion

Methodological Issues in First-Principle Calculations of CH3NH3PbI3 Perovskite Surfaces: Quantum Confinement and Thermal Motion Lucas Lodeiro,*Felipe Barría-Caceres, Karla Jimenez, Renato Contreras, Ana L. Montero-Alejo,* and Eduardo Menendez-Proupin* Cite This: ACS Omega 2020, 5, 29477−29491 Read Online ACCESS Metrics & More Article Recommendations * sıSupporting Information ABSTRACT: Characterization and control of surfaces and interfaces are critical for photovoltaic and photocatalytic applications. In this work, we propose CH3NH3PbI3(MAPI) perovskite slab models whose energy levels, free of quantum confinement, explicitly consider the spin−orbit coupling and thermal motion. We detail methodological tools based on the density functional theory that allow achieving these models at an affordable computational cost, and analytical corrections are proposed to correct these effects in other systems. The electronic state energies with respect to the vacuum of the static MAPI surface models, terminated in PbI2and MAI atomic layers, are in agreement with the experimental data. The PbI2-terminated slab has in-gap surface states, which are independent of the thickness of the slab and also of the orientation of the cation on the surface. The surface states are not useful for alignments in photovoltaic devices, while they could be useful for photocatalytic reactions. The energy levels calculated for the MAI- terminated surface coincide with the widely used values to estimate the MAPI alignment with the charge transport materials, i.e., −5.4 and −3.9 eV for valence band maximum and conduction band minimum, respectively. Our study offers these slab models to provide guidelines for optimal interface engineering. ■INTRODUCTION Hybrid halide perovskites (HHPs) are a family of original semiconductor materials, and HHP solar cells under research achieve efficiencies that compete with silicon photovoltaic cells. 1 The semiconductor-unique features include long carrier diffusion lengths, low recombination losses, low material cost, and band gap-tuning capacity. Optimizing interfaces is one of the critical tasks that has a direct effect on the performance of solar cells. 2 The selection of optimal interfaces should consider commensurable vectors to avoid stress, favorable atomic junctions that prevent dangling bonds, and finally, barrierfree charge transport, i.e., electronic states aligned between both materials. Therefore, a first step to the correct design of an interface requires knowing details of the structure and energy of the states of each part. Experimental advances have been critical to produce HHP thin films with the required degree of crystallinity to provide high cell performance. However, in the practice of building these solar cells, very little information is available on the characteristics of the surfaces obtained, neither their composition nor their energy levels. The same is true for their interfaces. By reviewing the electronic energy levels of perovskite films, determined essentially by combining ultraviolet photoelectron spectroscopy (UPS) and inverse photoemission spectroscopy (IPES), it is possible to realize that the reports differ greatly. 3 The differences are associated not only with the preparation methods and the ratio of precursors but also with the substrates and the exposure of the prepared films, air, light, and the thermal conditions. 2,3 Particularly in the Received: September 9, 2020 Accepted: October 21, 2020 Published: November 9, 2020 Article http://pubs.acs.org/journal/acsodf © 2020 American Chemical Society 29477 https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes. Downloaded via UNIV DE SEVILLA on January 24, 2024 at 13:00:51 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. methylammonium lead iodide perovskite (MAPI, MA+= CH3NH3 +), the HHP prototype, variations of up to 0.7 eV in energy levels are reported, changing the ratio of their precursors. However, most of the studies looking for alignments with MAPI levels still consider the reported VBM energy for bulk (−5.4 eV) and project the CBM energy at −3.7 or −3.9 eV, depending on the band gap considered. 2−6 By scanning tunneling microscopy (STM) and spectroscopy, it was possible to elucidate the atomic rearrangement of the surface and the density of state of the valence band of an ultrathin film of MAPI. 7 The measurements corresponded to the orthorhombic phase, and it was concluded that the surface ends in a layer rich in methylammonium iodide (MAI) in the plane (001), with two rearrangements of iodine coexisting. With this experiment, a real image of defects is shown for the first time, although large areas of the surface without defects are mostly observed. The computational model explains the influence of organic cations on the rearrangements of the observed atoms. The STM technique also provided the surface description of MAPbBr3and other HHP-mixed compositions. 8−10 However, in all cases, the information is limited to the low-temperature phase that is far from the experimental conditions in which the solar cells are manufactured. The elemental composition of the MAPI surface obtained by X-ray photoelectron spectroscopy (XPS), under controlled humidity conditions, has recently been reported for photocatalytic applications. 11 The analysis indicates that below 35% humidity, the MAI-type composition of the surface predominates, and at higher humidity, the surface is a PbI2type. However, the authors found no photoconversion activity on the MAPI surface in humidity conditions above 35%. The report supports that the MAI-type surface exhibits photocatalytic behavior because moisture generates surface defects exposing surface lead atoms as active sites. Prior to the STM 7 and XPS 11 results, theoretical models of the surfaces of each of the three MAPI phases had already been attained, including some defect proposals. For the orthorhombic phase, the (100) and (001) slab models appear with the lowest surface energies, which correlates with the largest coordination number of the surface Pb atoms. 12 In agreement to the subsequent STM experiment, the surfaces show MAI- type terminations, and the interactions between organic cations and iodine atoms are a stabilizing factor. Surface models corresponding to a pseudo-cubic unit cell have also been anticipated, this time, to study the adsorption of a molecule that simulate hole transport layers in solar cells. 13 The slabs representative of the (001) surfaces included the MAI- and PbI2-type terminations, both stoichiometric. Again, the orientation of the surface MA+influences the stabilization mechanism of the slabs. Binding energy calculations predict the preferred conformation of the adsorbed molecule on the surface, although energy level alignment was not considered. The surface of the tetragonal phase has been more studied since it is the phase most stable at room temperature, and several stoichiometric and nonstoichiometric (with defects) models have been proposed. 6,14−17 These studies suggest that the MAI-terminated (110) and (001) surfaces are energetically favored, although PbI2-type terminations are still possible. The greater stability of the MAI-terminated slab is again associated with the NH3 +group of MA+interacting with the iodine atoms, configuration where the CH3group is oriented outward, namely, top C. 17 The relative instability of the PbI2-terminated surface at vacuum could favor its energy of adhesion with any transporting layers in the cell. 15 In addition to the differences in the stability of the surface termination, reports show that electronic properties of each surface are highly sensitive to the size and the relative orientation of MA+within each slab. An initial report estimated the convergence in surface energy by modeling slabs up to five layers thick, where each layer extent is a half-unit cell. 15 The surface terminated by MAI did not reach convergence with respect to the band gap for these slab thickness, whereas the PbI2-terminated surface revealed a constant band gap. Other models show that both the band gaps and the energies of the frontier states with respect to the vacuum level change if the CH3or NH3 +ends of the MA+are regularly oriented toward the vacuum or if an apolar configuration of these cations is imposed within the entire slab. 17 Despite finding these variations, the agreement between the calculated frontier electronic states and the experimental data was reasonably achieved including spin−orbit coupling (SOC) effects as a posteriori correction. How to get accurate slab models? When surfaces and interfaces are studied by ab initio methods, computational cost restricts the slab thickness to a few nanometers, compared to the tens or hundreds of nanometers in experimental cases. In this quantum well regime, the electronic states are affected by quantum confinement. The correct description of the geometry and the electronic structure of the surface and the inside, requires a study of the convergence of certain parameters such as electronic state energies, surface energies and band gap, with respect to the thickness of the slab. Also, vacuum thickness is an important parameter when periodic boundary conditions are used for the three directions due to the fact that the model has replicas in the cleavage direction, and the vacuum needs to be long enough to avoid interaction with periodic replicas. In this work, we approach different computational methodologies in search of accurate representations of MAPI surfaces from slab models. The surface energy levels are modeled by means of a method that refers the band energies of the bulk to the energy levels of the slab, considering hybrid functional and SOC. Furthermore, the results were validated with the self-consistent wave function of sufficiently large slabs under the same level of theory. Correspondingly, the convergence of these energy levels with respect to the size of the slab and its corresponding vacuum layer is also studied. Finally, the surface structures have been characterized including thermal motion, so we can include features that static models cannot address. The displacement of the band edges due to thermal motion is quantified to finally discuss its implication in the alignment of MAPI bands with other materials. ■MODELS AND COMPUTATIONAL METHODS We modeled MAPI surfaces as slab models constructed by replication of the MAPI unit cell. We start from a model of the tetragonal crystal structure of MAPbI3, which is a 48-atom unit cell, named Unpol in ref 18. The name Unpol refers to the null polarization associated to the orientation of the C−N bonds in all (001) MAI layers. This structure represents a polymorphic configuration of the phase since the distortions are similar to those obtained within the configurational ensemble of this phase including thermal motion. 18 The Unpol unit cell was optimized by means of variable-cell relaxation. Total energies, forces, and stress tensors were ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29478 obtained from density functional theory (DFT) calculations, as implemented in the Quantum ESPRESSO code. 19 Two exchange and correlation functionals were used, the improved generalized gradient approximation (GGA) PBEsol 20 and the nonlocal vdW correlation optB88-vdW 21,22 (vdW). These functionals have different strengths. Dispersion-corrected vdW functionals are essential to describe binding of layered materials and molecular crystals. Although their use is not mandatory for HHP, they allow to obtain lattice parameters in very excellent agreement with the experiment. 23,24 However, vdW functionals cannot be used with the SOC, which modifies the band structure of the HHP. Hybrid functionals, including a fraction of Hartree−Fock exchange, together with the SOC, allow to accurately describe the band structure of CH3NH3PbI3. 23 However, the computational cost associated precludes its use in more computationally demanding simulation techniques like molecular dynamics. SOC is crucial to describe the conduction band minimum, which is composed mainly of Pb 6p atomic orbitals, and presents interesting effects like the Rashba splitting. An intermediate solution is in the PBEsol functional, which allows to obtain good lattice parameters, and can be combined with SOC to obtain suitable conduction bands. The main drawback of PBEsol is that the band gaps are still too low. Part of this problem can be overcome using the scissors operator or by expensive calculations using hybrid functionals. The Brillouin zone was sampled using a Γ-centered 3 ×3× 2k-point grid. The volume deviation with respect to the experimental value 25 reached −3.6 and −0.2% for PBEsol and vdW functionals, respectively. The lattice vectors were orthogonalized, neglecting small deviations from right angles. The spin−orbit coupling was not included here because its effect on the geometries is insignificant. Details on cutoffs, convergence threshold, and pseudopotentials employed are given in the Supporting Information. Slab models of MAPI(001) surfaces were obtained, cleaving bulk supercells at selected atomic planes with orientation (001). As shown in Figure 1, the simulation supercell contains a crystal volume (the slab) and a vacuum volume. The crystal volume is continuous in two directions, namely, xand y, and the surface is normal to the third direction, z. Periodic boundary conditions are used for the three directions. Therefore, these represent two infinite surfaces, which must be sufficiently separated not to interact. The separation of the surfaces is determined by both the slab and vacuum thicknesses, the effect of which will be studied in subsections Slab Thickness Convergence and Vacuum Width Convergence. In order to avoid using dipole correction, the slabs and their surfaces were constructed and relaxed, canceling the dipoles presented in surfaces, 26,27 mainly, by MA+rotations (see in subsection Surface Structure and Thermal Motion). We did not enforce the use of symmetry elements (mirror planes, inversion, for example) on the calculation, albeit the crystal geometry remained very close to its symmetric structure, which presents these elements of symmetry. Let us define the slab stoichiometric unit as a couple of adjacent MAI and PbI2(001) planes that span half of a unit cell with thickness c/2. Mind that c/2 roughly equals the lattice parameter of the cubic phase of MAPI. Henceforth, the number of stoichiometric units defines the slab thickness, as well as the number of MAI and PbI2planes. An additional PbI2 or MAI plane was added to each slab in order to have two symmetric surfaces (breaking the system stoichiometry) and to avoid internal electric field. Hence, a MAPI(001)-PbI2slab of n stoichiometric units has a total of nMAI layers and n+1 PbI2 layers, while for MAPI(001)-MAI, the numbers are permuted. The lattice vectors parallel to the surface can be linear combinations of the primitive vectors, allowing for surface reconstruction. Hence, we employ the slab nomenclature (m× m)×nfor the slab models, where mindicates the unit cell multiplication in each orthogonal direction along the surface. The index nindicates the number of stoichiometric units along the direction [001]. All the atomic coordinates within the slabs were relaxed with the aand bcell vectors fixed at their bulk values. The cell vector perpendicular to the surface was also kept fixed, and its optimal value was studied in subsection Slab thickness convergence. The electronic state energies were computed for the (1 ×1) ×nslab without and with the SOC, using a Γ- centered 3 ×3×1k-point grid. With this setup, a convergence study was performed (subsection Vacuum Width Convergence) in order to set the optimal values of slab thickness and vacuum width. Once the convergence criteria were established, the surface physical properties were studied using two, MAI- and PbI2- terminated, (2 ×2) ×6 slabs with a vacuum width of 35 Å, and the Γk-point. For these production slabs, the atomic coordinates were relaxed using the vdW functional. The GGA and vdW functionals seriously underestimate the band gaps and give erroneous values of the band energies, although the band gap error is fortuitously compensated when the SOC is not included. 23 To obtain accurate band energies in the bulk crystal, the HSE hybrid functional 28,29 was employed. The band edges were aligned with respect to the out-of-surface vacuum level using a two-step procedure explained in section Alignment with Respect to the Vacuum Level. In addition, selfconsistent slab calculations with vdW and the HSE hybrid functional were performed with the VASP code. 30,31 To assess the effect of thermal motion on the surface properties, ab initio molecular dynamics simulations were performed with the CP2K code, 32,33 for the (2 ×2) ×6 slabs. We took advantage of the efficiency of this code by using the Figure 1. Planar averaged (V(z)) and macroscopic (V(z)) electrostatic potentials, VBM and CBM energies with respect to the vacuum level for the (1 ×1) ×8 MAPI(001)-PbI2slab computed with PBEsol. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29479 hybrid Gaussian and plane wave method (GPW) to represent the wave functions and electronic density and the approach to extrapolate the density matrix implemented in the DFT QuickStep module. 34 Force evaluations were performed using the PBE functional with the Grimme correction scheme (DFT- D3) to account for vdW interactions. 35,36 The Born− Oppenheimer scheme, microcanonical ensemble, and periodic boundary conditions were used. Further details are given in the Supporting Information. The simulations started from the relaxed slab structures, with an initial temperature of 600 K. A time step of 0.5 fs was used to integrate the Newton equations of motion. The thermalization time was 5 ps. Afterward, the production runs were extended up to 40 ps. ■ALIGNMENT WITH RESPECT TO THE VACUUM LEVEL The energies of the valence band maximum (VBM) and the conduction band minimum (CBM), with respect to the vacuum level, are the negatives of the ionization potential (IP) and the electron affinity (EA), respectively. With DFT calculations, the VBM and CBM are usually determined by means of a two-step procedure. 37−39 First, the average electrostatic potential is obtained for the slab model, as shown in Figure 1 and Figure S2 for PbI2and MAI surfaces, respectively. From the slab calculation, the electrostatic potential V(x,y,z) is extracted and analyzed by means of two kinds of averages. First, the planar-averaged electrostatic potential is given as ∬ =+Vz AVxyz xy E() 1(, , )dd , xy 0 (1) where Axy is the transverse area of the slab. The constant E0is an energy shift applied to make V(z) = 0 in the vacuum region. In this way, the zero of energy is set to the electrostatic potential energy in vacuum. The vacuum reference determination is a current research topic. 40,41 Second, the convoluted macroscopic electrostatic potential is defined as ∫ =+′′ − Vz cVz z z() 1()d, c c /2 /2 (2) where cis the lattice constant in the direction perpendicular to the slab surface. As V(z) is periodic inside the crystal, V(z) should be constant at the central part of the slab. Also, V(z)is strictly constant and independent of the functional when it is computed for a unit cell of the crystal, i.e., without vacuum, as Figure 2 shows. The PBE0 42,43 hybrid functional was used in this case to compare. V(z) is precisely the reference used to transfer the VBM and CBM energies, from a three-dimensional (3D) crystal to a slab with vacuum. From a crystal calculation are obtained the average electrostatic potential V3D, and the energies EVBM 3D and ECBM 3D . Finally, the VBM and CBM energies with respect to the vacuum level are obtained as =+−  αα EVz E V() ( ), in 3D 3D (3) where α= VBM, CBM, and zin is any point in the region of constant V(z) near the slab center, as shown in Figure 1. This two-step process, with the difference Eα 3D −V3D transferred from a bulk calculation with a primitive cell, allows to correct the errors in VBM and CBM energies (leading to gap underestimation in GGA) by means of higher level calculations, such as hybrid functionals or the Green function method GW, over a primitive bulk cell. The slab calculation is prohibitively expensive with such high-level methods, but it can be done with GGA functionals, assuming that the electrostatic potential is well described by GGA. This assumption is true in most cases, also in MAPI, as demonstrated Figure 2 for non-SOC and SOC calculations with three functionals, over a bulk unit cell. Particularly, the different levels of theory and functionals do not change V(z) and V3D. There is just a little change in the V(z) profile (<60 meV, as shown in Figure 2) due to little changes in electron density. SOC shifts electron density from MAI layers to the PbI2layer, while for both hybrid functionals, the shift is from the interlayer space to MAI and PbI2layers. ■RESULTS AND DISCUSSION Slab Thickness Convergence. In this subsection, we study the convergence of the electronic state energies, band gaps, surface energies, and structural parameters with respect to the slab. We have classified the different slab states, through inspection of spatial localization of squared wave functions, into surface states and bulk-confined states (VBM and CBM). The electronic states were studied at the Γpoint, obeying the position of the bulk VBM and CBM (see the band-structure diagrams in Figure S1). The first checked item was the convergence of the energy levels with respect to the slab thickness. Figure 3 shows, for MAPI(001)-PbI2and MAPI(001)-MAI slabs, the dependence of the frontier energy levels (relative to the vacuum level). The blue and sky-blue symbols correspond to the quantum confined states of CBM (top plot) and the VBM (bottom plot). The CBM is the lowest unoccupied crystal orbital (LUCO) for both surface types. In MAPI(001)-PbI2, the higher occupied crystal orbital (HOCO) and the next (HOCO −1) are surface states, which are doubly degenerate because this slab model has two equivalent surfaces. Thus, the VBM is the HOCO −2. This is consistent with previous reports. 14 The surface states are caused by the cleavage of iodine octahedra around the surface Pb atom, which generates naked Pb atoms with dangled bonds. For MAPI(001)-MAI, the VBM-confined state is the HOCO, and no surface state has energy within the gap. The red and black symbols correspond Figure 2. Planar-averaged and macroscopic electrostatic potential for the MAPI unit cell, along z, for different levels of theory and functionals. (A) PBEsol. (B) SOC and other functionals, as difference with respect to (A). ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29480 to the CBM and VBM energies, respectively, obtained from the two-step procedure explained in section Alignment with Respect to the Vacuum Level. Its dependence with the slab thickness follows the behavior of the average electrostatic potential at slab center V(zin), from which they differ in constant values (eq 3). Geng et al. 15 reported that the MAPI(001)-PbI2band gap depends weakly on slab thickness, also being smaller than the band gap of MAPI(001)-MAI. However, they did not identify the HOCO and HOCO −1 in MAPI(001)-PbI2as surface states. This smaller band gap is consistent with our results (see Figure 4). Moreover, the surface state (HOCO) and CBM (LUCO) present a similar energy variation with respect to the slab thickness, leading to a nearly independent HOCO− LUCO gap (in the thinner slab thickness regime). Figures S3 and S4 show several squared wave functions (HOCO, LUCO, etc.) and their energies for (1 ×1) ×8, MAPI(001)-PbI2and MAPI(001)-MAI, respectively, averaged like eq 1. Here, the surface characters of the HOCO and HOCO −1 are clearly appreciated in the MAPI(001)-PbI2 case. These states are localized mainly on the first PbI2layer. For the other states, the bulk character and the quantum confinement effect are clearly appreciated in Figures S3 and S4. In the envelope function approximation,44 the band edge wave functions are products of Bloch functions (oscillating at the atomic scale) and slowly varying envelope function, which are of sinusoidal type in the quantum wells. 44 The envelope function tends to zero at the surface and has lnodes inside the slab for VBM −land CBM + lstates. There are also pseudonodes with smaller periodicity length along the slab due to the Bloch function oscillation. We see that, regardless of surface termination, the variation of slab state energies due to quantum confinement with respect to the slab thickness (see Figure 3) is almost the same for non- Figure 3. Convergence of energy slab levels, at the Γpoint, with respect to slab thickness. (A, C) PBEsol. (B, D) PBEsol+SOC. Energies are with respect to the vacuum level. Figure 4. Convergence of MAPI(001)-PbI2and MAPI(001)-MAI band gaps (for inner slab bands and surface slab bands) with respect to slab thickness. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29481 SOC (top) and SOC (bottom) calculations. The main difference between non-SOC and SOC calculations (beyond states energies changes and decreased band gap) is in CBM (and higher states) squared wave function shapes. The SOC wave functions have smoother profiles. The VBM- and CBM-confined states energies in Figure 3 exhibit quantum confinement behavior, as expected, in all cases. The CBM and VBM exhibit down (positive effective mass) and up (negative effective mass) energy shifts, respectively, for thicker slabs. While the energy convergence is faster in the VBM case, CBM presents quantum confinement even for the thickest calculated slab. This explains the VBM− CBM band gap overestimation (with respect to the bulk band gap) in all cases (see Figure 4). When n< 4, some VBM or CBM energies do not present smooth behavior, and the slabs are too thin to host bulk-like states. This shows that the bulkinner slab description is poor at this regime. Also, the two-step procedure VBM and CBM are not well described in this ultrathin slab regime. For n= 4, the energies are roughly converged, and for n≥6, they are tightly converged, in all cases. The two-step procedure is, in principle, free of quantum confinement. This conception was checked, and for n≥6, the VBM and CBM present asymptotic behavior. Finally, the surface states (MAPI(001)-PbI2slab models) exhibit energy splitting in the thinner regime. When the slab is sufficiently thick (n= 6 and n= 7 for non-SOC and SOC calculations, respectively), it is degenerated. This behavior shows that splitting is due to the surface state interaction (repulsion, like for the He2dimer) through the inner slab space and not through vacuum space. The last parameter studied is the surface energy. Figure 5 displays the surface energy as a function of the slab thickness. For a β= PbI2- or MAI-surface type, this energy can be computed as =−− ββ β S AEnEE 1 22 , n xy n slab, MAPI i k j j jy { z z z (4) where Eslab, β nis the total energy of an n-stoichiometric unit slab with two surfaces of type β,EMAPI is the total energy of tetragonal bulk MAPI, and Eβis the total energy of solid bulk PbI2or MAI. The latter is the sum of chemical potentials of the out-of-stoichiometry species in a β-rich condition. The double surface counting is corrected with the 1/2 factor and normalized with respect to surface area Axy. The surface energy, for both surface types and for non-SOC and SOC calculations, is converged for all thicknesses, as shown in Figure 5. In the MAPI(001)-MAI case, the surface energy presents periodic oscillations for even and odd number of stoichiometric units. This is due to the existence of little dipole asymmetries between MA+in both surfaces, but these differences are within the error bars. About 80% of error is associated with the k-point grid, and the rest is associated with higher cutoff. Also, the negative value in surface energy is explained by MA+inversion and independent of the MAI bulk reference. 45−47 Before geometrical relaxation, the surface energy is positive as expected (∼110 mJ m−2). In the MAPI(001)-PbI2case, the surface energy behaves smoothly because both surfaces are highly symmetric for even and odd number of stoichiometric units. The higher surface energy in this case is explained by the contribution of surface states. We have used the PbI2bulk structure of ref 48. The inclusion of SOC in the calculations does not significantly change the surface energies for both surface types. SOC has a little effect in total energies, reinforcing not having included it in the atomic relaxation calculations. In the MAPI(001)-MAI case, the SOC effect is negligible. Meanwhile, in the MAPI(001)-PbI2case, surface energy has rigid downshifts. As this case has surface states located in Pb and I atoms (heavier atoms), it is understandable that SOC shifts are not totally negligible. For all the thickness range considered, the geometrical parameters (see section Surface Structure and Thermal Motion) undergo minimal variations, matching the constancy of the surface energies. Square Well Approximation. As has been discussed before, inner slab states present quantum confinement. They are 2-D delocalized states in the (001) plane but confined along the [001] direction (zaxis). Moreover, the trends in the electronic state energies can be fitted with an infinite potential well model. Other considerations follow. The squared wave function envelopes, which modulate the Bloch functions, have sinusoidal behavior. The state energies’behavior seems consistent with a squared hyperbola curve, with asymptotic convergence at a high stoichiometric unit number. Also, the macroscopic electrostatic potential presents a pseudo-square well shape. Blancon et al. 49 published a work on Ruddlesden−Popper halide perovskites, where the effect of the PbI3inorganic framework thickness was studied both experimentally and theoretically. These perovskites are dielectric well structures formed by two-dimensional layers of halide perovskite separated by butylammonium spacer layers. The optical absorption spectra of these dielectric wells (which are analogues to our thinnest vacuum-slab well models) exhibit lower exciton energies for thicker perovskite thickness. Also, these dielectric wells show the quantum confinement effect on band gaps. Therefore, we propose the following analytical function, based on infinite potential well energies, to obtain the correct behavior relative to curvature and thickness dependence γ σ = + +∞ ≥ αα () En n En() 2 () , ( 4) c 2 2 (5) Figure 5. Convergence of MAPI(001)-PbI2and MAPI(001)-MAI surface energy with respect to slab thickness. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29482 where Eα(n) is the state energy as a function of the number of stoichiometric units n, with αdefined in Table 1. c 2 is the stoichiometric unit thickness (6.745 Å). γ,σ, and Eα(∞), are the adjustable parameters. The first parameter, γ, is proportional to the curvature of the function. The second, σ,isa correction factor to account for the effect of non-infinite and a pseudo-square well. Finally, Eα(∞) is the asymptotic value of state energy at infinite thickness. The parameter values are shown in Table 1. We used the confined state energy data from n=4−12 for all fits. Figure S5 shows that the functions Eα(n)fit the data extremely accurately, where the root mean square of residues and χ2values are negligible, with orders of 10−3and 10−5, respectively. The Eα(∞) values show that the CBM (in all cases) exhibit quantum confinement even for the thickest slab models (n= 12), while the VBM are less affected, as previously discussed. Furthermore, the two-step process for n≥4 provides CBM and VBM energies close to Eα(∞), showing faster convergence than the slab state energies. The band gaps calculated from the asymptotic energy values without SOC are 1.40 and 1.43 eV for PbI2and MAI surfaces, respectively, which match the value 1.43 eV computed in the bulk system. With SOC, the asymptotic energy values are 0.67 and 0.71 eV for PbI2and MAI surfaces, respectively, matching the bulk value 0.71 eV. The SOC correction (ΔSOC)isdefined as the difference between the SOC and non-SOC state energy function, Δ =− αα nEnEn() () () , SOC SOC (6) while the quantum confinement correction (ΔQC)isdefined as the difference between the asymptotic value of the state energy and the state energy function, γ σ Δ =∞− = − + αα () nE En n () ( ) () 2 . c QC 2 2 (7) Figures S6 and S7 show eqs 6 and 7, respectively, using data from Table 1.ΔSOC presents near-constant behavior for all cases in the fitting range. Also, the ΔSOC is practically equal to the value computed in the bulk system. ΔQC is the negative of the quantum confinement energy. This correction scheme allows to obtain the quantum confinement free state energies from a thin slab calculation performed with or without SOC. Vacuum Width Convergence. The convergence of the state energies with respect to the width of the vacuum region was studied using n= 5 slabs, both for PbI2and MAI surfaces. We have considered three vacuum thicknesses: 6.4, 17.9, and 42.3 Å. To begin with, we have verified that the relaxed atomic positions do not change with the vacuum width. Figure 6 shows the behavior of the slab state energies and the macroscopic electrostatic potential (V(zin)). The variations of all magnitudes are smaller than 20 meV, which is negligible. As shown in Figure S8, the surface energy variations are negligible (<0.2 mJ m−2), being much smaller than the error bars. Convergence Model Criteria. The previous subsections show the convergence studies for slab and vacuum thickness using a set of energetic and state variables. Here, we comment as a summary the most important results and their implications. In all slabs studied as MAPI(001) models, the VBM and CBM state energies present quantum confinement (greatest in CBM states), artificially increasing the band gap. We propose a correction scheme (ΔSOC and ΔQC), based on the square well approximation, to obtain energies of quantum confinement free states. The MAPI(001)-PbI2models exhibit occupied surface states caused by cleavage of a surface iodine octahedron, which triggers a reduction in band gap. The state energies obtained through the two-step process depend on the macroscopic electrostatic potential, and this parameter rapidly converges with the slab thickness. According to this method, the slab requires a minimum thickness of n= 4, but n=6is considered the optimal thickness to obtain state energies numerically close to the converged eigenvalues. The two-step process does not show the existence of in-gap surface states. The surface energies depend weakly on the slab thickness. Also, all parameters present weak dependence with the vacuum width. Surface Structure and Thermal Motion. It is well known that distortion effects suffered by MAPI perovskites, unlike the low-temperature orthorhombic phase, should be considered to adequately estimate their electronic properties. Many times, the calculations show that a single conformation, usually those with greater symmetry, unequivocally represents the minimum of internal energy. However, thermal motion actually breaks this symmetry, and the system must be represented by a set of distorted conformations representative of the phase. Table 1. Fitted Eq 5 Parameters for VBM and CBM (Non- SOC and SOC Calculations) αγ(eV Å2)σ(Å) Eα(∞) (eV) MAPI(001)-PbI2 CBM +808.65 12.111 −4.690 VBM −225.05 12.420 −6.088 CBMSOC +859.81 15.020 −5.219 VBMSOC −392.29 17.935 −5.891 MAPI(001)-MAI CBM +439.55 5.682 −3.096 VBM −97.64 10.965 −4.521 CBMSOC +338.31 7.218 −3.639 VBMSOC −181.56 16.131 −4.349 Figure 6. Convergence of MAPI(001)-PbI2and MAPI(001)-MAI state energies (in meV) with respect to vacuum thickness and relative to standard vacuum thickness. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29483 We simulate dynamic trajectories of the (2 ×2) ×6 slabs, and with them we assess the distortions by thermal motion. Starting at 600 K in the microcanonical ensemble, our surface models reach equilibrium at 291 and 294 K for MAPI(001)- PbI2and MAPI(001)-MAI slabs, respectively. With these trajectories, we were able to evaluate the structure of the surfaces, validate our supercell models obtained at 0 K, and calculate the band edge shifts from representative structures of the tetragonal phase surfaces. Finally, we also verify if a small deformation in the ab plane of the slab models could affect the band edge states. We begin the analysis, verifying that both static (0 K) and dynamic slab models preserve the tetragonal phase in the slab interior. As explained in ref 24, the tetragonal phase is geometrically described by a set of parameters. These include (1) dihedral, rotation, apical, and equatorial angles of the inorganic framework (see Table S2) and (2) angles defining the orientation of the MA+cations. First, we measured these parameters on the surface and internal layers of the slab models. As a reference, we consider reported mean values and their standard deviations for bulk MAPI. 24 The dynamic slab models have the inorganic framework parameters within the reported range for the bulk phase. Thermal motion mainly cancels out the differences in these parameters between the surface layers and bulk (inner) layers. Furthermore, there are no significant differences in the parameters between the two types of slabs. Similar to the dynamic models, the static slab models represent the tetragonal phase, although the apical angle (Pb−I−Pb along the zaxis) shows deviations. The orientation of the MA+differs significantly at surface and bulk layers. Allow us to consider the angle between the vector representing the C−N bond and the lattice ab plane, called polar angle. Figure 7 shows the distributions of polar angles for all the layers of MAPI(001)-MAI and MAPI(001)- PbI2slabs. The distributions indicate that the MA+tends to have a preferential dipole orientation in the surface layers that disappears in the bulk. For the MAPI(001)-MAI surface, the surface MA+has all the CH3groups pointing outward (to the vacuum) all the time. On the other hand, for MAPI(001)-PbI2 surfaces, the subsurface MA+presents both orientations, but the higher peak of the distribution corresponds to NH3 +groups pointing outward. In this case, it is observed that MA+changes its orientations during the trajectory. For inner (bulk) layers under both types of surface, the polar angles are distributed with two almost symmetrical peaks. There is also a surface effect on the inclination of the cation that is independent of the surface termination; the polar angle is 22°in the bulk, while on the surface it is greater than 50°degrees. The static model of MAPI(001)-MAI reproduces the preferred MA+orientation in the surface and agrees with previous reports with static 17 and dynamic 16 calculations. These reports associated those changes with weak interactions accounted through vdW functionals and Grimme dispersion correction, respectively. In contrast, we found this MA+ orientation even without considering the vdW corrections while using the PBEsol functional in convergence studies. We realize that this orientation is maintained by electrostatic interactions and mainly by hydrogen bonds of each surface iodine atom with two nearby cations simultaneously. Hence, each surface iodine atom increases its coordination number, which decreases surface energy and increases surface stability compared to other configurations. Under these interactions, the I−I interatomic distances are significantly reduced compared to the other situations. The dynamic model suggests, as will be discussed later, that it is possible to observe a pattern of rearrangement of iodine atoms at the surface due to these interactions. The static model of MAPI(001)-PbI2does not reflect the trends of the MA+orientations observed during the dynamics. In this case, the model in which no MA+layer has a net dipole moment (unpolarized slab) corresponds to the conformation of the minimum (internal energy). The formation of polarized slabs, in which all the cations of the first layer have the NH3 + group facing outward, is energetically feasible (ΔE= 0.1 eV/ MA+with respect to the unpolarized slab). In both structures, it is observed that NH3 +groups form hydrogen bonds, creating a network of similar interactions with the nearest I atoms. Unlike the MAPI(001)-MAI surface, the MA+orientation preference observed in this case as the polarized model could be justified by the electrostatic interaction between NH3 + groups and high electronic density in surface (due to the surface states). Figure 8 shows that both types of MAPI surfaces undergo relaxation along the caxis but not in the ab plane. There is a reduction in the Pb−I bond distance distribution of the outermost surface layers as compared to the bulk (inner) layers. In the figure, the labels Surf out and Surf in mean the Pb−I distance distributions in the outermost layer (out) and the next closest to it (in), while the Bulk label means the average distribution of two layers in the center of the slab. On average, the difference between Surf out and Surf in is around 0.25 on both types of surfaces. However, the distance Pb−Iin MAPI(001)-PbI2is reduced by 0.16 with respect to bulk, while this difference is 0.08 for the MAPI(001)-MAI slab. There is no significant difference between the Pb−I average distance in the bulk of both models, as expected. Notice that the corresponding distance of the static models, represented by vertical lines, reproduces the dynamic behavior in the case of the MAPI(001)-PbI2slab. The static model of the MAPI(001)-MAI slab hides this surface relaxation phenomenon observed while considering the thermal motion. This means that the static model is one of the possible Figure 7. Distributions of C−N polar angles for every layer of MAPI(001)-MAI and MAPI(001)-PbI2slabs. Vertical lines indicate the average angles for static slab models. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29484 configurations that represent the surface properties. However, distorted configurations should be also used to model the surface electronic properties. As noted above, our MAPI(001)-MAI model also shows regular rearrangements of iodine atoms at the surface compared to those in bulk, which are related to stabilizing interactions with MA+. According to the I−I distance distributions, the surface layers show iodine arrangements in a certain order depending on the crystallographic orientation on the surface (Figure 9, bottom). In contrast, the iodine atoms in the bulk show a unique (6.2 Å as average) although wide distribution due to thermal motion. During dynamic equilibrium, organic cations experience different orientations within the lattice ab plane, and the I−I distance distributions have a maximum at 5.4/7.1 Å in the upper surface and at 5.9/ 6.7 Å in the lower surface. Perpendicular to these regular distributions, both surfaces show that the iodines are distributed on average at 6.2 Å in wide distributions, as in the bulk. Notice that two different surfaces are obtained for this slab. The surface of the static model (Figure 9, top-left), as ideally ordered reference, allows to understand that, if two iodine atoms share hydrogen bonds with the same cations, then they approach to a distance of 5.4 Å (or move away to 7.1 Å), in average. The rest of the distributions correspond to situations where hydrogen bonds with organic cations are also presented but without a regular order as in a snapshot of the dynamic (Figure 9, top right). Notice that in the static model, the hydrogen bonds are maximized either in the surface or in the bulk. The iodine arrangement patterns found on the surface are consistent with the results obtained from the in situ STM experiment of this type of surface but for the lowtemperature orthorhombic phase. 7 Despite the limited size of our models, we believe that our evaluation can contribute to the analysis of surface phenomena in similar experiments in the tetragonal phase. It should be noted that in a convergence study, all thickness models show similar surface structural features as described here for static models. The fact that the geometry of the static models does not show appreciable changes explains the small variations found in the calculated surface energy based on the thickness of the slabs. Slab Energy Levels and Band Edge Shifts. Hereinafter, we summarize our results on the electronic properties of modeled MAPI(001)-PbI2and MAPI(001)-MAI surfaces. First, the electronic structure of static models is described, for which it is feasible to explore different theoretical approaches. Then, given the dynamic models, we can assess the influence of the structural distortions on the slab energy levels and the band energy shifts that may occur. Figure 10 shows the DOS and the plane-projected DOS for the MAPI(001)-PbI2surface comparing the optB88-vdW and HSE+SOC methods. Based on dynamic analysis, we consider adding a static model of this surface that we call polarized. This polarized model is a configuration in which the NH3 +groups of the surface layers are oriented outward. Figure 10 shows the DOS of both models, unpolarized (solid lines) and polarized (dashed lines), to compare. For both methods the HOCO is a surface state (one spin degenerate for each surface). The VB-confined states are below in energy, close to the bulk VBM, as expected from Figure 3. The method optB88-vdW shows more surface states with energies close to VBM, while with HSE+SOC, these surface states are a bit deeper in energy but still close to the VBM. The DOS rises more abruptly near the CBM with optB88-vdW because the neglect of SOC causes the CB to be flat in some directions. The quantum confinement is apparent Figure 8. Distributions of the Pb−I bond distances along the [001] direction and parallel to the ab plane for MAPI(001)-MAI and MAPI(001)-PbI2slabs. Surf out and Surf in refer to the outermost bonds and the next below, respectively. Bulk refers to bonds at slab centers. Vertical lines correspond to the average bond distances of the static slab. Figure 9. Top left: surface layer of relaxed static slab. Top right: MD snapshot of the surface layer. Bottom: distribution of I−I distances for I atoms in several (001) planes: upper surface, central planes (bulk), and lower surface. Vertical lines are the average I−I distances in static slabs. ACS Omega http://pubs.acs.org/journal/acsodf Article https://dx.doi.org/10.1021/acsomega.0c04420 ACS Omega 2020, 5, 29477−29491 29485