Advances in Modelling and Simulation of Halide Perovskites for Solar Cell Applications

Perovskite solar cells (PSCs) are attracting great attention as the most promising candidate for the next generation solar cells. This is due to their low cost and high power conversion efficiency in spite of their relatively short period of development. Key components of PSCs are a variety of halide perovskites with ABX3 stoichiometry used as a photoabsorber, which brought the factual breakthrough in the field of photovoltaic (PV) technology with their outstanding optoelectronic properties. To commercialize PSCs in the near future, however, these materials need to be further improved for a better performance, represented by high efficiency and high stability. As in other materials development, atomistic modelling and simulation can play a significant role in finding new functional halide perovskites as well as revealing the underlying mechanisms of their material processes and properties. In this sense, computational works for the halide perovskites, mostly focusing on first-principles works, are reviewed with an eye looking for an answer how to improve the performance of PSCs. Specific modelling and simulation techniques to quantify material properties of the halide perovskites are also presented. Finally, the outlook for the challenges and future research directions in this field is provided.


Introduction
As often highlighted on the public media, the issue of climate change and global warming is one of the most serious global challenges that human beings face nowadays [1]. This is an inevitable result from mass burning of fossil fuels by ourselves to run factories, drive cars and operate buildings during the past decades. In fact, this brought a release of a vast amount of greenhouse gases, such as carbon dioxide (CO 2 ) and methane (CH 4 ), and caused another challenge for energy due to the exhaustion of fossil fuel resources. To mitigate the catastrophic effects of global warming by reducing greenhouse gas emission and ensure the sustainable energy supply, more and more peoples and authorities express their growing interests in the renewable and clean (or low-carbon) energy sources, such that the majority of national and global energy policies include an exploitation of clean energy sources. Among several kinds of clean energy sources available on the globe, solar energy is the most promising source due to its abundance and eternal nature. Therefore, it is required to develop efficient techniques and devices for harnessing solar energy, which can contribute to resolving the environmental and energy challenges [2].
The most effective devices for harvesting solar energy should be undoubtedly solar cells invented in 1960s, which convert solar light directly into electricity by photovoltaic (PV) effect. So far, several kinds of solar cells have been developed, and can be found in the market; the market share is as 69.5% for single crystalline silicon (sc-Si), 23.9% for polycrystalline silicon (pc-Si) and 6.6% for thin-film polycrystalline CdTe solar cells in 2015 [3]. In fact, the most important factors determining the competence of solar cell in the market are power conversion efficiency (PCE), fabrication cost and device stability. According to the report of Green et al. [4], certified PCEs of conventional PV modules are 28.8% for sc-GaAs, 24.4% for sc-Si, 19.9% for pc-Si, 18.6% for CdTe or Cu-In-Ga-Se (CIGS) alloy, and 13% for quantum dot and dye-sensitized solar cells. For these conventional solar cells, it is difficult to get a high efficiency, a low cost and a high stability simultaneously; when the efficiency increases, the fabrication cost also increases. That is why sc-GaAs solar cells with the highest cost and high stability and dye-sensitized solar cells with the lowest efficiency and low stability could not be commercialized. For the most widely used Si solar cells, however, both the efficiency and cost are not satisfactory enough to become competitive with fossil fuels in the electricity market. Therefore, it is our ardent desire to invent an innovative new type of solar cell with high efficiency, low cost and high stability together.
In accordance to this desire, perovskite solar cells (PSCs) have emerged with an initial PCE of 3.8% in 2009 [5], and evolved rapidly with their PCEs to over 10% within a few years [6,7,8,9]. After several years of intensive and extensive research, the efficiency arrives at over 22% [10,11,12], which is an astonishing breakthrough in the PV history when compared with the development periods of the conventional solar cells mentioned above. Due to simple fabrication processes (low processing cost) and highly abundant raw materials (low materials cost) [13,14], moreover, their fabrication cost is a third time as low as c-Si solar cells, resulting in a short energy payback time and low overall CO 2 emission. Such huge success is mostly originated in their key components, halide perovskites with a chemical formula of ABX 3 , where A and B are the monovalent (organic or inorganic) and divalent (metallic) cations and X is the halide anion, respectively. These halide perovskites have been used as a solar light absorber and/or a charge carrier transmitter due to their optimal properties for these aims [15,16,17,18,19,20,21,22,23,24]. For instance, the archetype methylammonium lead iodide (CH 3 NH 3 PbI 3 , MAPI hereafter), which was firstly hired in PSCs, has appropriate band gaps, optimal photoabsorption coefficients, weak exciton binding energies, high charge carrier mobilities, and long carrier diffusion lengths. However, there are still several fundamental issues to be addressed for commercialization of PSCs: (i) high stability and long lifetime [25,26,27,28,29], (ii) scalability for large area modules [30] and (iii) low toxicity [31].
Atomistic modelling and simulation of materials become essential in materials science, as they can provide valuable insights for known material processes and properties, and further predict the shortest route for finding new functional materials that meet the requirements [2,32,33]. In particular, first-principles methods based on quantum mechanics can quantitatively determine most of material characteristics of crystalline solid, such as lattice constants, electronic structure, linear response properties, and transport properties, with a typical accuracy of 1% relative error to the experiment by using input data of only atomic informations and known (or hypothetical) lattice structure. Being an interdisciplinary subject, the power of computational materials science can be supported by two aspects: (i) in hardware, rapid and ever increase of processor speed and memory capacity, and (ii) in software, continuous progress of simulation algorithms and materials theory. As such, the state-of-the-art atomistic modelling and simulation methods such as density functional theory (DFT), many-body perturbation theory (MBPT), and pair-wise interatomic potential molecular dynamics (MD) have been successfully applied to the halide perovskites. Recently, there exist several reviews on theoretical works of halide perovskites: short reviews on the nature of chemical bonding [34] and on the electronic and ionic motions [35], perspectives on the entropy [36] and on hybrid halide perovskites [37], and reviews focusing on MD simulations [38] and on optoelectronic properties [39,40].
In this review, we inclusively discuss recent progress in first-principles atomistic modelling and simulation of halide perovskites in both hybrid organic-inorganic and purely all-inorganic forms. We focus our special attention on how the calculated material properties can be linked up with the performance of PSCs. Specific modelling and simulation techniques to elucidate material properties of halide perovskites are also discussed. We consider the following material properties: (i) crystalline structures, (ii) electronic structures and optical properties, (iii) phonon dispersions and material stability, (iv) defect physics and ionic diffusion, and (v) surface and interface. The accuracies of calculations adopted in each work is also on discussion in comparison with the available experimental data. Finally we present the challenges to modelling and simulation of halide perovskites and future research directions.

Crystalline structures and polymorphism
The first step of solid-state modelling and simulation is to define the Bravais lattice (lattice constants and angles) and structural factor (atomic coordinates) of unit cell. This can be done simply from experimental data such as X-ray diffraction (XRD) and scanning or transmission electron microscopy (SEM or TEM) with elemental analysis (e.g. energy dispersive X-ray spectroscopy: EDX or EDS) for a known material, where data mining is a crucial tool to analyze large data sets [41] and to predict the structure-property cartograms [42]. Then, the initial structure is locally refined by performing structural optimization, which can provide a preliminary accuracy estimation of adopted theoretical method by comparing the determined lattice constants with experimental data, and a supplementary information of elasticity (bulk modulus) by postprocessing the calculated data of total energy vs unit cell volume. If the structure is unknown, the crystal structure prediction should be performed in different ways: (i) simply making a set of hypothetical crystalline structures and then comparing their total energies, and (ii) global structure optimization using a variety of genetic and evolutionary algorithms [43,44]. Once the refined structure is obtained, we can perform the analysis of chemical bonding by measuring the bond length and bond angles, and proceed further calculations for material properties, getting the relation of structureproperty.

Geometric factors for 3D halide perovskites
In ABX 3 halide perovskite structure, the large monovalent cation A is bonded with the nearest twelve X halide anions, forming an AX 12 cuboctahedron, while the smaller divalent metal cation B forms corner sharing BX 6 octahedra in a three-dimensional (3D) framework. It should be noted that A cation can be organic, forming hybrid organic-inorganic halide perovskites, and inorganic, leading to all-inorganic halide perovskites. At a glance, it seems that there are a great number of combinatorial cases simply with monovalent A and divalent B cations based on the periodic table. However, there are some geometric factors restricting the formation of stable 3D halide perovskite structure. The first thing is the Goldschmidt tolerance factor [45], t, defined as, where r A , r B and r X are the ionic radii of A, B, and X ions, respectively. It was empirically known that with t = 1 for a perfect cubic lattice, 0.8 ≤ t ≤ 1.0 for most stable perovskite structures with tetragonal, orthorhombic and rhombohedral (or trigonal) lattices. Otherwise, nonperovskite structures are formed; hexagonal structure is formed when t > 1 and different structures are found when t < 0.8 [46,47].
The second constraint is known as the octahedral factor, defined as an ionic radii ratio between B cation and halide anion, µ = r B /r X [49]. This reflects the stability of BX 6 octahedron, and should be in the range of 0.44 ≤ µ ≤ 0.9 for a stable perovskite.
Together with the tolerance factor, the octahedral factor can provide a constraint to select stable halide perovskites among many possible cases. For instance, we can draw a 2D map of ionic radii of A cations (Rb + , Cs + , MA + , FA + , EA + ) and X anions (F − , Cl − , Br − , I − ), which are marked along the abscissa and ordinate in Fig. 2, for typical B-site metal cations of Pb 2+ , Sn 2+ and Ba 2+ . Their ionic radii with proper coordination numbers (CN = 12, 6, 6 for A, B, X ions, respectively) are provided in Refs. [50,51].
In Fig. 2, we can see the intersections located within the polygonal (shaded) region formed by tolerance limits of these constraint factors, which indicate a potentiality of stable perovskite formation with corresponding A, B, X elements, and those outside the region, which indicate a fail of perovskite formation. It should be noted that there is an ambiguity to define ionic radii due to the nonsphericity of organic molecular A cation and reduced electronegativity of heavy halide anions, and thus some modifications have been proposed, such as effective molecular A cation radii [48] and modified B-site metal cation radii from empirical bond lengths [47].  [40] with permission of the American Chemical Society.

Octahedral tilting or distortion and polymorphism
Once a certain compound with ABX 3 formula is proven to form a stable halide perovskite due to its proper tolerance and octahedral factors, we should also consider a possibility of modification in the shape of BX 6 octahedra. In fact, BX 6 octahedra can be tilted or distorted because of slight rotations or displacements of constituent atoms or molecules at finite temperature (and possibly pressure). It is worth noting that the octahedral tilting is known to be associated with antiferroelectricity, while the octahedral distortion caused by cation off-centering yields a ferroelectricity (i.e. creation of spontaneous electric polarization) that can enhance charge carrier separation and allow the photovoltage to exceed the band gap by thousands of times [52,53]. Such octahedral tilting or distortion may lose a symmetry of crystalline unit cell, being closely related with phase transition driven by temperature [54].
In general, the archetype cubic perovskite structure is observed at high temperature, and when decreasing temperature, lower symmetry phases are found in the order of tetragonal→orthorhombic→monoclinic and/or rhombohedral.
For the typical case of hybrid organic-inorganic halide perovskite MAPI, which has been playing the major role in advancing PSC technology, its cubic α-phase with space group Pm3m is only found at high temperature, and phase transitions to tetragonal β-phase (I4/mcm) occurs at 327.4 K, and to orthorhombic γ-phase (Pna2 1 ) at 162.2 K [55]. In this case of MAPI, the orientation of MA cation is of great importance in determining the octahedral tilting and the crystalline phase [53,56]. The organic MA + cation has two modes of reorientation: (i) methyl and/or ammonium rotation around the C−N axis and (ii) whole molecular rotation of the C−N axis itself [53]. In the cubic phase, three different orientations are identified as <100>, <110>, and <111> for the C−N axis, which lower lattice symmetry from the cubic Pm3m to the pseudo-cubic Pm and R3m phases, and the <110> model was found to be the most energetically stable configuration by ab initio MD [56] and structural optimization [57] in agreement with the experiment [53]. By performing ab initio MD simulations, the activation barriers were calculated to be 117 meV for reorientations of the axis [58] and 13.5 meV for ion rotation [59,60]. For the tetragonal phase, the C−N bonds in MA distribute in a vertical way, while for the orthorhombic phase in a parallel way, as shown in Fig. 3 [62] revealed that these inorganic halide perovskites exhibit lattice instabilities in the cubic phase, confirmed by phonon softmode (imaginary frequencies) at the Brillouin zone boundary, and thus the octahedral tilting and distortion that is by definition antiferroelectric in nature can occur spontaneously. A quantitative picture for phase diversity can be obtained by calculating free energy through phonon calculation, which will be discussed later.

Low dimensional and double lead-free halide perovskites
With respect to the environmental impact, there is a concern on the toxicity of lead in the major halide perovskites such as MAPI and CsPbI 3 , promoting an extensive search for alternatives. Although there has been extensive research to try simple substitution of Pb 2+ in ABX 3 structure with alternative divalent cations, it was turned out to be challenge. There is also research interest in another approach to search for multivalent elements. This resulted in low dimensional perovskite structures such as A I 3 B III 2 X 9 ("3-2-9") layered (2D), A I 3 B II X 5 ("3-1-5") single chained (1D) and A I 4 B II X 6 ("4-1-6") isolated octahedra (0D), and double perovskite structures A I 2 B I B III X 6 ("2-1-1-6") [66,67]. Xiao et al. [68] has reported their recent work for Pb-free halide perovskites with such different structural dimensionalities, which have a wide range of band gaps and PCEs of those-based solar cells. In particular, 3D double perovskites such as Cs 2 AgBiX 6 (X = Br, Cl), called elpasolites, have received much attention as promising absorber materials [69,70,71,72,73]. Volonakis et al. [74] has reported a series of double perovskites with B I = Cu, Ag, Au and B III = Bi, Sb, predicting their band gaps in the range from 0.5 to 2.7 eV for Bi system and from 0.0 to 2.6 eV for Sb system. So far, there is no report for double perovskite-based solar cells, which might be due to difficulties in synthesizing uniform thin films of the correct phase and chemical composition.

Solid solutions
There is an ever-increasing need towards higher performance of solar cells, represented by efficiency and stability. Although the PCE of MAPI-based solar cells reaches as high as 22%, being comparable with sc-Si solar cells, their instability in air is relatively low, which could be major barrier to commercialization. Making solid solutions or mixing at A-site and/or X-site with two or more similar elements is an efficient way to satisfy the need of enhancing the performance by increasing the stability, improving the charge carrier transport, and tuning the band gap [75,31]. Experimentally, it is much easier to realize such mixing than double perovskites. There have been several first-principles investigations for perovskite solid solutions such as hybrid perovskites, [57,76,77,78], MAPb(I 1−x Cl x ) 3 [79] and (Cs/MA/FA)PbI 3 [80], and inorganic perovskites, CsPb(X 1−x Y x ) 3 (X, Y = I, Br, Cl) [81] and (Rb x Cs 1−x )SnI 3 [82].
When mixing larger I − anion with smaller Br − or Cl − anion, the crystalline lattice can be reduced and a phase transition from the tetragonal to the cubic can occur at room temperature. Meanwhile, mixing MA + cation with larger FA + cation can tune the octahedral tilting with the proper Goldschmidt tolerance factor.
To make modelling of solid solutions, we can use supercell (SC) method at the sacrifice of high computational cost. When employing the pseudopotential plane wave (PP-PW) method, the virtual crystal approximation (VCA) approach can be used alternatively to reduce the computational cost, but the accuracy or reliability should be checked carefully before going into the main calculation. In the VCA approach, the virtual atom is introduced and its pseudopotential is constructed by averaging the relevant potentials and wave functions of constituent atoms [83,84]. There are some technical problems in these methods; random distribution of alloying atoms for SC method, while local effect from the interaction between the constituent atoms for VCA method.
It has been proven that the VCA approach can produce the reliable tendency of lattice constants in cubic phase as a linear function of mixing content, i.e. Vegard's law, for hybrid halide perovskites [78,79] (Fig. 4).

Electronic and optical properties
Since the halide perovskites are typically used as a light absorber in most applications, the electronic and optical properties are of primary importance. After modelling the crystalline solid using a unit cell or supercell and refining their cell by performing the structural optimization, the electronic energy band structure and the correspond-  electrons and holes). Other optical properties such as exciton binding energy and photoabsorption coefficients are determined by using the frequency-dependent dielectric constants, which can be routinely calculated within the density functional perturbation theory (DFPT) [85,86,87] without or with the effect of electron-hole interaction by using the Bethe-Salpeter approach [88].

Spin-orbit coupling and many-body effects
In general, standard DFT exchange-correlation (XC) functionals within local or semi-local approximations such as local density approximation (LDA) and generalized gradient approximation (GGA) severely underestimate the band gap of inorganic semiconductors, due to the inherent nature of DFT as ground state theory and the artificial self-interaction between electrons, although they can provide reliable structures and stabilities. Surprisingly for MAPI, the GGA functionals (PBE [89] or PBEsol [90]) can yield a band gap in excellent agreement with experimental values within ±0.1 eV [77,78,79,91]. This is due to a fortuitous error cancellation between the GGA underestimation and the overestimation by the lack of spin-orbit coupling (SOC) [92,93,94,95,96], which inversely indicates an importance of relativistic effects in halide perovskites. From the analysis of partial density of states (PDOS) for halide perovskites ABX 3 , the valence band maximum (VBM) consists of an antibonding coupling of X

p-states and B s-states, while the conduction band minimum (CBM) is dominated by B
p-states (see Fig. 5). It should be noted that for MAPI the major of occupied molecular orbitals of MA is found deep (∼5 eV) below the VBM and the minor (but not negligible) contribution is found ∼0.5 eV below the VBM, indicating the hydrogen bonding interaction between the organic MA moiety and the inorganic PbI 6 octahedra [97]. In APbX 3 perovskites the heavy Pb atom has a significant SOC effect of lowering the CBM, whereas the lighter Sn atom has a relatively weak SOC effect, resulting in again large underestimation of band gap by GGA functionals for ASnX 3 perovskites [98,99].
Moreover, the GGA functionals cannot accurately describe dispersion of the valence band even in Pb-based perovskites.
The relativistic effect can be approxmated by first-order scalar relativistic (SR) and higher order SOC contributions. Umari et al. [97] demonstrated that SR-DFT barely changes the band obtained by DFT, while SOC-DFT largely lowers the conduction bands for hybrid iodide perovskites, MAPbI 3 and MASnI 3 . When the SOC effect is in consideration, splittings of the conduction bands, known as Rashba effect, occurs as a result of interaction between the magnetic moment (spin) of the electron and the local electric field [100]. In fact, this electromagnetic force displaces electrons in the momentum space, acting on up and down spins in opposite direction. For the case of cubic halide perovskites, neglecting SOC yields direct band gaps at the band edge points of R and M in the Brillouin zone, as shown in Fig. 6. When SOC is turned on, both the valence and conduction bands around R split into symmetrical valleys for the case of MAPI (see Fig. 6(a)). More importantly, since the splitting is much more evident in the Pb 6p conduction band compared with the I 5p valence band, the CBM slightly shifts along the R→ Γ line, resulting in the vertical energy difference of 25 meV between the CBM and VBM at R point [94]. Thus the band gap changes into the indirect mode, which affect the function of light absorber as we discuss later. Motta et al. [94] revealed that the orientation of MA molecular cation has the same effect on the band gap transition. For the all-inorganic halide perovskites like CsPbI 3 , however, such relativistic spin-splitting does not occur due to a centre of inversion symmetry.
Although the SOC effect causes the separation of Pb 6p into p 1/2 and p 3/2 , no splitting of the band extrema is observed at the high symmetry points, as shown in Fig. 6(b).
Since both the standard DFT and SOC treatment based on local or semi-local XC functionals are insufficient to describe the electronic structures of halide perovskites, it is necessary to adopt more delicate approaches. One way is to use a screened hybrid  [97,103]. In addition, in-clusion of SOC in GW calculation increases band dispersion so significantly that the standard parabolic approximation can no longer be applied to the band extrema, which differs from experiment indicating parabolic band extrema in MAPI.

Tuning band gap by substitution and mixing
As the key factor of light absorber, the band gap of the halide perovskites can be tuned by different ways: (i) static volume change, (ii) temperature change, and (iii) chemical substitution. When decreasing the volume by acting a kind of small perturbation or by increasing pressure, the band gap was found to monotonically decrease in DFT calculation with an additional benefit of transition from indirect to direct band gap for MAPI [106,107]. Also, it was found that the out-of-phase band-edge states are stabilized as the lattice expands [60]. Regarding the temperature effect, there is no theoretical work yet, although the band gap was observed in temperature-dependent photoluminescence to decrease with decreasing temperature from 1.61 eV at 300 K to 1.55 eV at 150 K for MAPI [108]. We focus on the third way of chemical substitution at each site of ABX 3 , which is more genuine for materials design than the physical volume effect [34,109,110].
The A-site cation of MAPI can affect only indirectly on the band gap by changing the crystal structure due to its molecular orbitals far away from the VBM. When replacing MA cation with smaller molecular cation such as NH + 4 and H + , the band gap was found to lower by 0.3 eV with contraction of cubic lattice from 6.29 to 6.21 Å in QSGW+SOC calculation for NH 4 PbI 3 [103], and to be less than 0.3 eV with a lattice constant of 6.05 Å for HPbI 3 [60]. It is important to recognize the relationship between the ionic radii and the local structure. In fact, both NH 4 PbI 3 and HPbI 3 adopt nonperovskite chain or layer structures due to their lower tolerance factors than 0.81, which give rise to an instability of the octahedral networks with respect to tilting. Replacing MA with larger FA can also result in structural distortions [111]. Meanwhile, substituting inorganic cation Cs + with a smaller ionic radius leads to a widening of band gap as 1.73 eV in cubic CsPbI 3 [112,113,114], being suitable for top cell material in tandem cell. However, it is difficult to form its cubic phase at room temperature.
In this situation, it is desirable to adopt the mixing technique, which can promise a synergy effect of enhancing stability against moisture and keeping efficiency. Mixing can be done in various ways such as MA/FA, Cs/MA, Rb/Cs [82], Cs/MA/FA [80], and Rb/Cs/MA/FA [115].
Substitution or mixing on the B site can directly alter the conduction band dominated by B p states without the severe change of crystal structure. For example, replacing Pb with isovalent Sn in MAPI reduces the band gap by ∼0.3 eV due to a slight downward shift of CBM by weaker Sn 5p states than Pb 6p states [116]. It also induces a phase transition from the tetragonal I4cm phase for MAPI to the pseudocubic P4mm phase for MASnI 3 , being beneficial to the light absorber. However, Sn 2+ is less chemically stable in the octahedral network due to its ease of oxidation to Sn 4+ , resulting in the degradation of PSCs. Another isovalent Ge 2+ is further unstable owing to its lower binding energy of 4p electrons. When mixing Pb with Sn to form solid solutions MA(Sn 1−x Pb x )I 3 , the band gap was found to follow the quadratic function of mixing content x (see Fig. 7), occurring a band inversion associated with a systematic change in the atomic orbital composition of the conduction and valence bands [116,117,118].
Since the X-site anion dominates the valence bands, its substitution can be expected to change the VBM [119]. As going from I to Br to Cl, the valence band composition changes from 5p to 4p to 3p, resulting in a monotonic increase in electron binding energy (higher ionization potential). Accordingly, the band gap increases from 1.5 eV for MAPI to 2.10 eV on Br [78] to 2.70 eV on Cl substitution [79]. For the case of FAPbX 3 , the substitution of Br for I increases the band gap from 1.48 eV to 2.23 eV as found in experiment [120]. For the former case of MA(Sn 1−x Pb x )I 3 , mixing I with Br or Cl can give band gaps changing according to the quadratic function of mixing content (Fig. 7), where b is the bowing parameter reflecting the fluctuation degree in the crystal field and the nonlinear effect from the anisotropic binding [78,79,121]. If the bowing parameter would be zero, the band gap also could follow the Vegard's law like lattice constant.

Transportation of charge carriers
In the halide perovskites, the major charge carriers are conduction electrons and holes created by several external factors. They can couple each other by an electrostatic interaction to form a kind of quasiparticle, namely exciton. Therefore, effective masses of these charge carriers and exciton binding energy are of importance in checking the suitability for charge transport layer in solar cells.
Effective masses of conduction electrons and holes can be calculated from the lower conduction band and upper valence band near the R point for the cubic phase, which can be approximated as a parabolic function of momentum k as follows.
This parabolic approximation can be slightly broken when considering the SOC effect in the hybrid halide perovskites as discussed above. Therefore, the calculated effective masses depend somewhat on the theoretical method and the complexity of band dis-  [104]. Since the mobility is inversely proportional to the effective mass, the halide perovskites are expected to exhibit ambipolar and large diffusion lengths for electrons and holes [95].
An exciton can be viewed as hydrogen atom with differences attributed to (i) a replacement of the electron and nucleus masses by the effective masses and (ii) an introduction of a dielectric constant ε of material. Then, the exciton binding energy can be obtained using the effective masses and dielectric constant in the following equation, is the reduced effective mass. The static dielectric constant can be calculated by DFPT method without or with the electron-hole coupling as will be discussed below. If this modelling is valid as the extending radius of the lowest bound state (a * = ε m e m * r a B , a B the Bohr radius) is larger than the lattice constant, the exciton is the weak exciton, called as the Mott-Wannier exciton. Otherwise, it is a Frenkel exciton, being extremely localized. For MAPI (MAPbBr 3 ), the exciton binding energy E b and effective radius a * were calculated to be 45 (99) meV and 3.0 (1.9) nm [78], indicating that the excitons are likely to be of the Mott-Wannier type.
Charge carrier transport can be characterized by the diffusion length, defined as the average length of carrier traveling before recombination. This can be calculated using the diffusivity D and charge carrier lifetime τ by L d = √ Dτ. The diffusivity can be directly related to the mobility (µ) by the Einstein relation µ = Dq/k B T due to the fact that transport is limited by carrier scattering and carrier effective mass. The L d must be long enough so that the charge carriers can reach the contacts in solar cells, as reported to be considerably longer in MAPI than other semiconductors. This can be partly attributed to the defect-tolerance of hybrid halide perovskites [123]. Providing effective masses of <0.2m e for MAPI, the carrier mobility was calculated to be modest (<100 cm 2 /(V·s)) compared to the inorganic semiconductors Si or GaAs (>1000 cm 2 /(V·s)) [124]. Such limitation of carrier mobility can be caused by strong scattering. Zhao et al. [125] have found that by applying PBE+SOC approach the electron-acoustic phonon couplings in MAPI are weak, and charge carriers are scattered predominantly by charged defects or impurities.

Dielectric constant and photoabsorption coefficient
The frequency dependent complex dielectric functions can be calculated within DFPT without the electron-hole coupling effects. Although such excitonic effects can be calculated by solving the Bethe-Salpeter equation for the two-body Green's function without or with local field effect [126], ignoring electron-hole coupling still yields quite a reasonable result in good agreement with experiment for the small band gap semiconductors [88]. Given a self-consistent QSGW potential, the the polarizability P(q, ω) is obtained using the random phase approximation (RPA), and then the inverse Coulomb interaction [127]. The macroscopic dielectric function ε M is given as follows, Then, the frequency dependent photoabsorption coefficient is given as follows,  Fig. 8(c)) [92].
The MA molecular dipoles screen the Coulombic interaction between photo-excited electrons and holes in the Pb-I sublattice [128,129,130], leading to a reduction of exciton binding. The MA cation is situated in a cuboctahedral cage, and thus, when going from I to Br to Cl, the size of cage decreases, inducing a restriction of molecular motion and a decrease of static dielectric constant with a reduction of exciton screening.
In accordance to these tendency, the photoabsorption onset and the first peak shift to a higher photon energy, i.e. shorter wavelength light, as increasing the mixing content.

Phonon and material stability
Crystalline materials in general exhibit phase transition upon changing temperature and pressure. A certain phase with a specified crystal structure can exist in stable state in the certain range of temperature and pressure, and at the critical value, changes into another phase. Such stability at the thermodynamic condition is called phase stability.
On the other hand, crystalline materials can be chemically decomposed into several components, which may be in crystalline and/or gaseous states. This is called chemical stability. During such chemical decomposition, heat can be generated (exothermic reaction) or absorbed (endothermic reaction). For the former case, the chemical decomposition can occur spontaneously, being intrinsic chemical instability that can be considered at different thermodynamic conditions, while for the latter case the reaction can be triggered and progressed with extrinsic factors such as moisture, light and heating, being extrinsic chemical instability. Both the phase stability and intrinsic chemical instability at the certain thermodynamic condition can be predicted by performing phonon calculations to determine the entropy [36].
The Gibbs free energy of a compound as a function of temperature T and pressure P is given as where F(T, V) is the Helmholtz free energy as a function of temperature and volume V.
For the PV term, the total energies are calculated as varying the unit cell volume, fitted into the empirical equation of state for solid [135] to obtain E(V) function, and then the pressure is estimated by conducting a differentiation like P = −(∂E/∂V) T [136].
Within the adiabatic approximation, F(T, V) can be separated into ionic vibrational and electronic contributions [137,85], In the electronic Helmholtz free energy, F el (T, V) = E(T = 0 K, V) − T S el , the T S el term is ignored because the electronic temperature effect is negligible for nonmetallic systems at the room temperature vicinity [138]. Within the quasiharmonic approximation (QHA), the ionic term F vib can be calculated as follows [136,137], where ω(V) is the phonon frequency as a function of volume, M the atomic mass, g(ω) the normalized phonon DOS and ω L the maximum of the phonon frequencies. The phase stability or chemical intrinsic stability can be estimated by the Gibbs free energy difference between the two phases or the products and the reactants.

Phonon dispersion and phase stability
It is crucial to get a precise insight of phonon dispersions in the halide perovskites for understanding material processes such as ionic transport, and recombination and scattering of charge carriers, as well as material stabilities. The phonon dispersion is computationally accessible via lattice dynamics calculations, in which the dynamic matrix is constructed via the Hessian matrix as the second derivatives of the total energy with respect to the atomic displacements that can be obtained directly from DFT calculations. The Hessian matrix (or interatomic force constants: IFCs) can be calculated either in real space by performing force calculations on a series of symmetry-inequivalent displaced structures with an use of supercell [139,140] or in reciprocal space through perturbation theory (e.g. DFPT [85]). Diagonalizing the dynamic matrix yields a set of eigenvectors (phonon modes) and eigenvalues (phonon frequencies) [36]. The lattice dynamics in DFT is one of the highest expensive calculations, and thus, classical or semi-classical MD based on interatomic potential function [141] can be used to overcome the limitations of DFT.
For the hybrid perovskite MAPI in the cubic, tetragonal and orthorhombic phases, the phonon bands were calculated in the harmonic approximation [142,143]. Confirming that the calculated phonon spectra were in quite well agreement with the experimental data, these calculations indicate that the lowest bands below 100 cm  [53,144,145,146]. It was found that MAPI exhibits double-well instabilities at the Brillouin zone boundary and remarkably short phonon quasiparticle lifetime, being associated with low thermal conductivity, and the optical phonon scattering is stronger than the acoustic one at room temperature in these materials [52,145,147].
Vibrational entropy was proved to play a crucial role in determining the stable phase, e.g. between the pseudocubic and tetragonal phase for MAPI, such that the materials favour structures that maximize the number of soft intermolecular interactions as the temperature rises [148].  It is important to understand the anharmonic effects in the phase transition of halide perovskites. When the potential energy of a crystalline solid is expanded as a Taylor series of ionic displacements, only the second term (d 2 U/dr 2 ) is considered in the harmonic approximation with the temperature-independent frequencies and infinite lifetime, while the effects of temperature and first-order anharmonicity can be considered in QHA [150]. It is associated with the imaginary (or negative) frequency, i.e. soft mode. While all the phonon modes have positive frequencies in the tetragonal and orthorhombic phases ( Fig. 9(a)), there are two imaginary frequency acoustic modes in the cubic phase, of which centres are around the R[ 2π a ( 1 2 , 1 2 , 1 2 )] and M[ 2π a ( 1 2 , 1 2 , 0)] points, corresponding to the <111> and <110> directions ( Fig. 9(b)). As a common feature of the perovskites, these soft modes indicate a double-well characteristic in the potential energy surface, and thus a dynamic instability of the cubic structure as observed in inelastic X-ray scattering measurements [52]. Such instabilities represent antiferroelectric distortions associated with collective rotation and tilting of the corner-sharing octahedral framework that can be observed directly in MD [59,151]. In the double-well potential, the cubic phase is a saddle point between two equivalent broken-symmetry phases ( Fig. 9(d)), and the barriers are 37 and 19 meV for the R and M modes, being comparable to k B T at room temperature [147]. The phase transitions from cubic structure can be understood as a condensation of these soft modes, as proved by solving the time-dependent Kohn-Sham equation describing the nuclear motion in the double-well potential [147,152,153].
The soft phonon modes are also found in the all-inorganic halide perovskites. For the case of CsSnI 3 perovskite, the temperature dependent lattice dynamics were carried out for the four different phases within QHA, revealing the strong anharmonic effects such as soft modes [65,154]. However, the calculated phase transition temperatures fail in quantitative agreement with experimental values, requiring further studies to eliminate these discrepancies. Yang et al. [62] systematically investigated the phonon dispersions of CsSnX 3 and CsPbX 3 (X = F, Cl, Br, I) perovskites, verifying the phase instabilities of cubic structures related to the spontaneous octahedral tilting. In all the eight cases, double-well potentials as a function of distortion amplitude (Q), given by with fitting parameters a and b, were found for soft phonon modes with barrier heights ranging from 108 to 512 meV, assessing the chemical and thermodynamic driving forces for these instabilities. When compared with the hybrid halide perovskites, the soft mode is found surprisingly even at the zone centre point Γ (Fig. 9(c)), which is related to the ferroelectric distortion. Marronnier et al. [155,156] deeply investigated this anomaly at Γ for the cubic and tetragonal phases of CsPbI 3 . As a polar mode, the soft phonon mode found at Γ for cubic CsPbI 3 is linked to the displacements of Cs + cation in one direction and of I − anion in the opposite direction. In spite of such displacements, ferroelectricity has not been observed at a macroscopic scale due to the oscillations along this polar soft mode. Volume relaxation with tight convergence thresholds (10 −4 Ry/bohr for the force and 10 −14 for the phonon self-consistency) and frozen phonon calculations remove the soft modes at Γ. Based on those results, it was concluded that the strongly anharmonic mode at Γ will not condensate at lower temperatures, whereas the remaining phonon instabilities at the edge points of M and R are responsible for soft modes that condensate at lower temperatures to induce the phase transition [155].

Thermodynamic miscibility in solid solutions
Forming solid solutions or alloys by mixing equivalent two or three elements at the same site can provide a direct and easy route to improve the performance of materials.
As such, perovskite solid solutions were found to have a controlled band gap and improved stability; in lead-based hybrid halide perovskites APbX 3 , the best efficiency can be obtained by mixing of organic cations (MA/FA) on the A-site and halides on the X site. In such mixing cases, it is important to understand whether the solid solution is stable against phase separation in the entire range of composition, where configurational entropy can play an important role in describing the disorder of site occupancy [36,77].
Thermodynamic miscibility for mixing two components, e.g. (1 − x)ABX 3 + xA ′ BX ′ 3 , can be calculated by the Helmholtz free energy difference given as The first term, internal energy difference, is written as follows [157], where E k (x) is the DFT total energy of the alloy with the mixing content x and the configuration number k, and g k is the degeneracy representing the number of configurations with the same energy E k . At room temperature, the internal energy of mixing can be well fitted into the subregular solution expression [76,157], with the fitting parameters, α and β. Then, the entropy of alloy is given in the ideal solution limit as follows, which is expected for a random alloy at high temperatures.
Brivio et al. [76] have provided the thermodynamic insight for photoinstability in the hybrid halide solid solutions MAPb(I 1−x Br x ) 3 . By generating different configurations (using the SOD code [158]) and applying the generalized quasichemical approximation method [159], they calculated the internal energy, configurational entropy, Helmholtz free energy of solid solutions and phase diagram as functions of the alloy composition ( Fig. 10). At low temperatures, the free energy curve is asymmetric and positive, indicating the existence of a miscibility gap. When increasing the temperature, the curve becomes symmetric and negative since the probability of sampling   3 . In part (a), the solid lines with blue, green, and red colours refer to 200, 300 K and high T limit, and dashed line represents the convex hull. In part (c), the purple and pink regions are the binodal and spinodal regions. Stable solid solution can be formed in the white region. Adapted from [76] with permission of the American Chemical Society. (d) Helmholtz free energy of Cs 1−x Rb x PbI 3 in cubic phase. Adapted from [149] with permission of The Royal Society of Chemisty.

Chemical stability
It is well known that the hybrid halide perovskites are weak on external actions such as humidity, ultraviolet light and heat. Also, whether they are stable compounds with respect to the chemical decomposition or not is an important issue. If they are intrinsically stable, the degradation of PSCs can be suppressed by thoroughly protecting the device from the external actions. To check the chemical stability, the following reaction is generally suggested for chemical decomposition of ABX 3 perovskites, Then, the Gibbs free energy difference between the products and the reactants, ∆G = G ABX 3 − (G BX 2 + G AX ), is an estimation of the chemical stability.
At zero temperature and zero pressure condition, simply the DFT total energy difference, i.e. formation energy, can be calculated for this aim [78,79,160]. It was found that for MAPI all the possible phases are unstable for chemical decomposition with PBE functional but including the vdW correction lets the orthorhombic phase stable. As shown in Fig. 11(a) [79] (Fig. 11(b)), which agreed well with the experimental findings. At finite temperature, the vibrational contributions were considered by performing the lattice dynamics calculations (mostly via DFPT [85]), revealing that interestingly for MAPI and MASI there is no significant change of Gibbs free energy between the reactant and products while for CsSnI 3 its free energy decreases faster than that of CsI and SnI 2 . In other work for MAPI using experimental data for chemical potentials, the free energy difference was calculated to be 0.16 eV/f.u. at 300 K, indicating that the finite temperature effects tend to stabilize the perovskite structure against the spontaneous chemical decomposition under standard thermodynamical condition [161]. For mixing Cs with Rb in CsPbI 3 , the free energy differences tend to increase as the Rb content x increases with a turning point. This turning point decreases as the temperature increases, as shown in Fig. 11(c). The phase diagram for chemical decomposition of CsPbI 3 into PbI 2 and CsI can be obtained as shown in Fig. 11(d), indicating that CsPbI 3 is stable in the temperature range from 0 to 600 K and the pressure range from 0 to 4 GPa [149], being well fitted with experiment. impurity and interstitial impurity, or higher-dimensional defects including dislocations, grain boundaries [162] and precipitates. Such defects can mediate the ion diffusion or can be mobile themselves inside the solid, which are suggested to be responsible for a current-voltage hysteresis of device and degradation related to moisture and light exposure [163,164].

Defect formation and ion diffusion
Simulations of defect-related phenomena are in highly cost from the viewpoint of memory capacity and computational time. In these simulations, sufficiently large supercells must be adopted to ensure a reliable accuracy of (charged) defect formation energy, because of the errors caused by the finite sizes of supercells [165,166]. To determine the charge transition levels formed by charged defects, the electronic structures of perfect and defect-containing crystals (modelled by supercells) are calculated with a high accuracy, e.g. using hybrid functionals such as HSE06 [102] and/or considering SOC effects coupled with many-body theory (GW). Formation energies of defects depend on the chemical potential of ingredient elements, which need to consider thermodynamics. The formation enthalpy of a point defect with a charge state q is calculated using the grand canonical expression [167,168,169], where each term can be calculated by first-principles method coupled with thermodynamics. On the other hand, modelling and simulation of ion diffusion, which are mostly vacancy mediated, need to suggest preliminary migration paths based on the empirical bond valence sum (BVS) approach [170], and perform structural optimizations to fix the local configurations, and finally the minimum energy path calculations using the nudged-elastic-band (NEB) approach [171] to determine the energy barrier.

Defect formation energy and transition levels
The most important defects in the halide perovskites ABX 3  reflecting the growth conditions. In general, equilibrium conditions for the product with its pure constituents (e.g. MAPI is in equilibrium with MA, Pb and I 2 ) give the constraints for the relation between their chemical potentials, resulting in two different scenarios; halide-rich and halide-poor growth conditions. Therefore, two diagrams for defect formation energy can be obtained under these conditions. From these defect formation energy diagrams, it is possible to identify the dominant defect that has the lowest formation energy and thermodynamic charge transition levels.
In the case of tetragonal MAPI, the dominant point defects were found to be the acceptor-type lead vacancy (V −2 Pb ) under I-rich condition [169,172,173], while donortype MA +1 i under I-poor condition [172]. Buin et al. [174] have got a little different findings, such that the major acceptor defects are V −2 Pb , V −1 MA and I −1 i , whereas donor defects are V +1 I and Pb +2 i . Among these, defects V +1 I , MA +1 i , and V −1 MA were found to possess the lowest formation energies over the entire band gap. The low formation energy of V Pb is related to the energetically unfavourable s-p antibonding coupling [172].
In the case of MAPbBr 3 and MAPbCl 3 with the cubic phases, Pb −3 Br(Cl) has somewhat comparable formation energy with V Pb [173,175]. Interestingly, the dominant defects that have the low formation energies exhibit shallow transition levels in the band gap, while the defects that have higher formation energy (and thus difficult to form) would have deep-trap levels. In addition, all vacancies yield shallow traps and resonances within the band, implying that charge carriers can still move easily to VBM and CBM [174]. It has been reported that only interstitial iodine defect I i is a deep trap and non-radiative recombination centre for hole [176,177], pointing out the strong dependence of transition levels on the selection of XC functional (HSE+SOC should be used) [178]. Therefore, it can be suggested that point defects in the hybrid halide perovskites should not contribute to the density of deep traps (recombination centres) that directly controls the diffusion length of charge carriers. However, this can be changed according to the growth condition. For example, the formation energy of deep-level antisite Pb I may be low enough to significantly contribute to the density of recombination centres under I-rich condition [173,174]. In this sense, the optimal growth conditions should be halide-poor for MAPI, and halide-rich for MAPbBr 3 and MAPbCl 3 [173].
The experimentally observed unintentional doping for both n-type and p-type MAPI without any added dopant can be explained by creation of neutral vacancy pair defects, i.e. Schottky defects and Frenkel defects. It was found that Schottky defects such as V PbI 2 and V MAI , which may dominate the defect formation in the halide perovskites under stoichiometric growth conditions, do not make a trap state [179]. The elemental defects derived from Frenkel defects play the role of unintentional doping sources, implying that n-type and p-type doping can be controlled by carefully choosing the proper atomic composition in the growth procedure [180]. Such vacancy pair defects were also considered in the water-intercalated and mono-hydrated MAPI in order to reveal the mechanism of PSC degradation upon exposure to moisture [169]. From the calculated binding energy of the vacancy pair defects, it was concluded that the formation of V PbI 2 from the vacancy point defects V −1 I and V +2 Pb is spontaneous in these three compounds, while the formation of V MAI is less favourable than the formation of individual vacancies V −1 I and V +1 MA in the hydrous compounds. Moreover, all these Schottky defects exhibit deep trap levels in the hydrous compounds (Fig. 12), giving an evidence of degradation of PSCs under humid condition [169].
Defect calculations have been performed for other halide perovskites. By using GW+SOC method, Li et al. [181] investigated the defect physics of cubic CsPbI 3 .
They found that under the Pb-rich conditions the vacancy point defects V Pb and V I are the dominant acceptor and donor defects and can pin the Fermi energy in the middle of the band gap due to their comparable formation energies. Under the Pb-poor conditions, the acceptor V Pb has a high density and thus can generate a high density of hole carriers. While the antisite Pb Cs acts as the recombination centre under the Pb-rich condition, there are no such centres with a low-energy formation energy under the Pb-poor condition [181]. There exist recent theoretical works of defect physics and chemistry for the hybrid layered perovskite (CH 3 NH 3 ) 2 Pb(SCN) 2 I 2 [182], bismuth-based lead free double perovskites [183], and defect passivation in the hybrid perovskites using quaternary ammonium halide anions and cations [184].

Ion diffusion
Once the intrinsic ionic defects formed inside solids, they can migrate according to the minimum energy pathways, as the hybrid perovskites exhibit ionic charge transport as well as electronic conduction in experiments. Such ionic diffusion in the halide perovskites is closely related to the long-term stability and power conversion efficiency of PSCs. In particular, the anomalous hysteresis of the current-voltage (J-V) curves, which is a severe disadvantage of PSCs, and giant switchable photovoltaic effects, are importantly invoked by ion migrations [185]. The chemical decomposition of MAPI 0 0.5 1 1.5 upon moisture exposure can also be explained partly by ion diffusion. In the simulations, the main tasks are to identify the migrating species and their migration pathways, and to quantitatively determine the underlying energetics, on condition that the formation of such ionic defects have been already clarified.
Eames et al. [186] investigated the diffusions of vacancy points (ions), not interstitials, in cubic MAPI, based on the fact that in ABX 3 halide perovskites vacancymediated ion diffusion is the most acceptable diffusion process while interstitial migra-tion has not been observed in experiments. It should be noted that experiments probing ion migration, however, typically cannot distinguish between different migration pathways. They suggested three vacancy-mediated ion migration pathways including I − migration along the octahedron edge, Pb 2+ migration along the diagonal direction and MA + hopping into a neighbouring vacant site ( Fig. 13(a)). The activation energies for these ionic migrations were calculated to be 0.58, 2.31, and 0.84 eV, respectively [186].
When compared with the measured values for hysteresis, the calculated value for I − is in good agreement with the measured value, 0.60−0.68 eV, whereas for Pb 2+ and MA + vacancies they are higher than the measured values, confirming that these ion migrations are difficult to occur [186]. It should be noted that there is a wide spread of activation energies for ion migrations in both theory and experiments [187]. Azpiroz In these calculations, the common conclusion is that halide anions and MA or FA cations are the major ionic carriers due to their relatively low activation energies, which is consistent with the experimental findings and can explain the J-V hysteresis [190,191]. Based on the dilute diffusion theory, it was suggested that replacement of MA with larger cation FA can suppress the hysteresis and prevent the aging of PSC performance [189].
Recently, we investigated the vacancy-mediated ion migrations in cubic MAPbX 3 (X = I, Br, Cl), and further water molecular diffusion in their water-intercalated and mono-hydrated phases, aiming at uncovering the role of water molecule in the chemical decomposition and thus the degradation of PSC performance [192]. It was found that the insertion of water into the MAPI reduces the activation energies of ion migrations ( Fig. 13(c)), indicating the easy decomposition of the hybrid halide perovskite under humid condition. Also, the activation barriers for ion and water migrations become higher from X = I to Br to Cl. These results indicate that the decomposition of halide perovskite occur through a multi-step process such as from the water intercalation to hydration and to decomposition, identifying the crucial role of water molecule in this process. Egger et al. [193] investigated the hydrogen migration in the tetragonal MAPI to explain the hysteresis and material stability. Using the supercells containing a hydrogen impurity with different charge states, they found that the crystal structure may relaxed significantly for charged hydrogen impurities, collective iodide displacements can enhance proton diffusion, and the migration barriers for proton transfer are relatively low. Based on these findings, it was suggested that in the hybrid halide perovskites the hydrogen-like defects, introduced either extrinsically or intrinsically, may be mobile and play an important role in explaining the hysteresis effects and stability issues [193].

Surface and interface
The device performance and long-term stability of PSCs can be greatly affected by material processes occurring at surfaces, grain boundaries and interfaces as well.
In reality, all of the bulk materials have the surface due to their finite size, and as such, the halide perovskite bulk crystals have several types of surfaces classified hklindices and terminations. Moreover, the perovskite films are fabricated mostly through solution processing methods, which make them polycrystalline and thus the formation of grain boundaries unavoidable. With respect to the device structure, it is necessary to adopt electron and/or hole transporting layers that contact to the photoabsorbing layer, halide perovskite, creating the interfaces such as TiO 2 /MAPI for interface with electron transporting layer. Possibly important material processes at surfaces and interfaces are halide diffusion, ion accumulation, charge carrier transport and recombination at the defect states. For these phenomena, first-principles modelling and simulations can also provide useful microscopic insights as well.

Surface phase diagram and electronic states
In the standard simulation techniques using the three dimensional periodic boundary condition, the surfaces can be modelled by slab with a supercell, which consists of atomic layers and vacuum layer. The number of atomic layers and thickness of vacuum layer must be checked to ensure that the DFT total energy difference, e.g. surface formation energy, is not influenced by these modelling parameters. By applying ab initio atomistic thermodynamics, the formation energies of surfaces with various indices and terminations are calculated as functions of chemical potentials of constituent elements at finite temperature and pressure, producing a surface phase diagram.
The structural and electronic properties of the tetragonal MAPI surfaces with various indices and terminations have been investigated by using DFT calculations with the rev-vdW-DF functional and the inclusion of SOC effect [194,195]. Among low-index surfaces considered in that work, tetragonal (110) and (001)  Identifying the constraint relations between the chemical potentials, the surface phase diagrams were drawn for the four types of these surfaces by calculating the Gibbs free energy differences as functions of chemical potentials of iodine and lead. Through these diagrams, stable terminations in each index surface can be determined at different growth conditions, i.e. Pb-rich (poor) and I-rich (poor) conditions. Also, the relaxed structures and electronic properties were analyzed in detail. Based on those calculations, it was concluded that a vacant termination is more stable than the PbI 2 -rich flat termination on all the surfaces under thermodynamic equilibrium growth conditions of bulk MAPI. The flat terminations on the (001) and (110) surfaces were found to have surface states above the valence band of bulk (possibly attracting the photogenerated holes), which can be efficient intermediates of hole transfer to the adjacent hole transport materials with smaller energy loss [194].
As discussed above, the trapping of charge carriers at defects on the surfaces and grain boundaries is one of the detrimental factors in PSC performance, because it causes nonradiative recombination loss, deteriorates the carrier lifetime, and originates the J-V hysteresis. Uratani and Yamashita [196] investigated the types of surface defects responsible for carrier trapping in tetragonal MAPI (001) surface by performing the HSE+SOC calculations. They constructed the (2×2) surface models with a vacuum region of 15 Å and three types of terminations, i.e. MAI, flat and vacant ( Fig. 14(a)).
Several kinds of intrinsic point defects such as vacancies (V I , V MA , V Pb ), interstitials (I i , Pb i ) and antisites (Pb I , Pb MA ) were identified to be formed on each terminatedsurface, and their formation energies were calculated at I-rich, moderate and Pb-rich conditions. The energy levels of defect states were also drawn comparing with the VBM and CBM of each system. Using the calculation results, they concluded that under the I-rich condition iodine interstitials I i on flat and vacant surfaces are responsible for the carrier trapping, while under the Pb-rich condition V I on vacant surfaces and Pb i defects are carrier traps. Therefore, the Pb-rich condition is better than the I-rich condition in terms of PV performance, which is consistent with the experiments.  [196]. Copyright 2017 American Chemical Society. (b) Lowest-energy structures of water under the first layer of the PbI 2 -terminated P + surface (left) and of water initially placed at the hollow site of the PbI 2 -terminated P − surface, and the corresponding reaction energy profiles. Reprinted with permission from [197]. Copyright 2015 American Chemical Society.
Surface modelling and simulations are also important for understanding the interaction of halide perovskites with water, which can be related with the chemical decomposition of especially MAPI under humid condition. Such water-assisted chemical decomposition of MAPI can be explained by hydrolysis [60] or hydration [192]. In these processes, the initial step is the adsorption of water on the MAPI surface and its diffusion (or penetration) into the bulk. Liu and co-workers explored the water effects on the detailed structure and properties of MAPI surfaces based on ab initio molecular dynamics using the Gaussian-type double-ζ polarized basis sets (DZVP) and the PBE+vdW XC functional [198]. Tetragonal MAPI (001) surfaces with a MA + termination were considered as the base surfaces. It was found that the water adsorption energy on the MAPI (001) surface is about 0.30 eV and the diffusion barrier of water molecule from the surface to inside region is only about 0.04 eV, indicating that the water can easily penetrate into the bulk. Koocher et al. [197] found that the water adsorption is greatly affected by the orientation of the MA + cations close to the surface.
The activation barriers of water penetrations were calculated to be 0.023 eV in the case of P + (CH 3 -end of the MA molecule towards top surface) and 0.271 eV in the case of P − (NH 3 -end) on the MAI-terminated (001) surface ( Fig. 14(b)). Therefore, it could be suggested that controlling orientation via poling or interfacial engineering can enhance the moisture stability. Choosing the MAI-terminated (110) surface in MAPI and the MABr-terminated (110) surface in MAPbBr 3 , Zhang et al. [199,200] calculated the adsorption energies of water molecules to be 1.74 and 1.51 eV with PBE while 1.87 and 1.67 eV with PBE+vdW on MAPI and MAPbBr 3 surfaces. They suggested that the deprotonation of the MA cations followed by the desorption of the CH 3 NH 2 molecules is a primary step in the degradation mechanism.

Interface
In PSCs, the halide perovskite photoabsorbers are sandwiched between the electron transport layer (ETL) and hole transport layer (HTL), so that the photoexcited electrons and holes can be efficiently transported to the collecting electrodes through these layers. In the case of MAPI, for example, it is in contact with mesoporous TiO 2 scaffold as ETL in one side and with organic Spiro-OMeTAD material as HTL in other side, thus creating heterojunction interfaces at these contacts [201]. Such interfaces can be considered as the large defects, where nonradiative recombination of electrons and holes may occur, resulting in the decrease of PV performance. In fact, charge carrier extaction from the absorber to the transporter can be obstructed at the interfaces by several factors such as interfacial energy barriers originated from imperfect band align-ment and charge carrier recombination driven by interfacial defect traps [202]. On the other hand, the interfaces can be weak against the infiltration of external particles such as water molecule and hydrogen atom or proton, which can trigger or facilitate the degradation of PSCs. Therefore, understanding the interface phenomena and further engineering the interface structure and composition are very important for improving the PSC performance.
Mosconi and co-workers investigated the interfaces of TiO 2 to MAPI and also MAPbI 3−x Cl x with a small Cl content (∼4%) using SOC-DFT calculations [188,203,204]. The most stable tetragonal (110) and pseudocubic (001)  surface. From the simulation results, they conclude that MAPI and MAPbI 3−x Cl x tend to grow (110) surfaces on TiO 2 due to their higher binding energies to TiO 2 than (001) surfaces, and interfacial Cl atoms increase the binding energy of the MAPbI 3−x Cl x (110) surface compared to MAPI. Through the electronic structure calculations, it was revealed that the interaction of the perovskite with TiO 2 changes the interface electronic structure to a stronger interfacial coupling between the Ti 3d and Pb 6p conduction band states and to a slight upshift of TiO 2 conduction band energy [203]. In order to investigate the effect of interface defects, vacancies V I and V MA were created at the perovskite surface exposed to the vacuum (HTL contact) and at the oxide-contacting surface, based on the already established fact that under working conditions V I should diffuse towards the HTL while V MA towards the ETL side [188]. When compared with the non-defective interface, the outermost valence band states of the perovskite intrude into the band gap of the TiO 2 substrate, while the conduction band states are found inside the manifold of the oxide's conduction states. The presence of defects at in-terfaces was found to modify the band alignment, cause the bending of the perovskite bands close to the TiO 2 surface, and create the trap states [188].
For the case of all-inorganic perovskites, CsPbBr 3 /TiO 2 heterostructure has been investigated using the PBEsol+HSE06 functional [205]. From the experimental observations for CsPbBr 3 /TiO 2 interfaces, the interface modified CsBr and PbBr 2 layers were confirmed to be made during the synthesis, leading to the corresponding CsBr/TiO 2 and PbBr 2 /TiO 2 interfaces. From the analysis of the charge population, it was found that large charge transfer occurs between Cs and O atoms, while much smaller charge transfer between Br and Ti atoms, being caused by the interface effects.
Accordingly, the charge difference before and after the interface formation was drawn; the charge accumulation occurs on Ti and O atoms while the depletion on Cs and Br atoms, implying the electron transfer form perovskite to the TiO 2 . From the integrated local density of states projected along the direction perpendicular to the interface, the staggered gap offset junction was observed for both interfaces, revealing the driving force of the charge transfer through the interfaces. Moreover, the band gap in the regions far away from the interface is almost consistent with that of individual CsPbBr 3 or TiO 2 , whereas in the interface vicinity it decreases. The interface electronic states were found to locate on the conduction band edge of CsPbBr 3 , indicating the shallow gap states. It was suggested that these shallow gap states may increase the absorbed photons and promote the optical absorption [205].
Regarding the degradation issue, the interaction between the MAPI surface and liquid water environment was investigated by making the slab interface model consisted of MAPI surface and water molecules and performing ab initio MD simulations [206].
The perovskite surfaces were made by cutting (2 × 2) slabs from the bulk tetragonal MAPI crystal, leading to (001) surfaces with three kinds of terminations: MAI, PbI 2 and PbI 2 -defective terminations. To generate the defective surface, six PbI 2 units were removed out of a total eight PbI 2 units from the PbI 2 -terminated (001) surface [194].
The vacuum region was filled with water molecules, the numbers of which were determined from the experimental density of liquid water as 284, 226 and 235 for MAIterminated, PbI 2 -terminated and PbI 2 -defective surfaces, respectively. By analyzing the MD trajectories, the MAI-terminated slabs were found to undergo rapid solvation according to a nucleophilic substitution of I by H 2 O, desorption of MA molecules, and thus a net dissolution of a MAI molecule. On the contrary, the PbI 2 -terminated surfaces were stable against such solvation process, although the percolation of a water molecule can occur to form a hydrated bulk phase. It was also concluded that PbI 2 defects may facilitate the solvation of the surface and initiate the degradation of the entire perovskite.

Summary and outlook
Being inspired by astonishing experimental findings for PSCs, a lot of computational works for the halide perovskites have been carried out within the latest five years in order to reveal the underlying mechanisms of perovskites' physics and chemistry and indicate a way for finding new perovskites with an improved performance. In this review, we have presented a critical discussion about the theoretical findings and conclusions for crystalline structures, electronic and optical properties, lattice dynamics and material stabilities, defect physics and ionic diffusions, and surfaces and interfaces in the hybrid organic-inorganic and purely all-inorganic halide perovskites, which have been made based on first-principles calculations. Specific theories, methods, and modelling approaches for these calculations have been provided, trying to point out the discrepancy and limit of the modelling and calculations when compared with the available experimental data. Special attention has been paid to making it clear the structure-property relationship in these materials, which was in close association with PV performance (efficiency and stability). We have demonstrated the predicting power of first-principles materials design with some predictions for new halide perovskites unexplored by experiment or optimal mixing ratios in some solid solutions.
In order to commercialize PSCs with a realistic impact on global energy market share, the efficiency and stability of PSCs need to be further improved, as the efficiency is still under the maximum theoretical efficiency of 25∼27% estimated by first-principles modelling [207]. In particular, inherent chemical instability of halide perovskites under working conditions such as humid, light and thermal environment should be of the primary urgency in commercialization. On one hand, the performance enhancement can be realized by identifying the most ideal device architecture, as some materials-architecture combinations have shown high efficiency and stability. On the other hand, due to a rich versatility of halide perovskites with respect to the chemical composition and structural phase, it is desirable to find the most optimal choice of composition and structure that can give the best performance of PSC. Relying on computational or coupled experimental-computational method will be the best way to achieve this aim with respect to the time and cost.
Band gap engineering is the first thing to do. Once the chemical composition and crystalline structure are identified, the band gap can be calculated in almost autonomous way with tolerable wasting time and demand. However, much more timeconsuming methods beyond standard DFT such as hybrid functional and GW should be needed to get a reliable band gap and especially band alignment. Then, lattice dynamics should be performed to estimate the material stability at finite temperature and pressure. Although quasi-harmonic approximation can be used to get phonon dispersion, anharmonicity that the perovskites exhibit inherently should be considered; this is not yet well-established. Anharmonicity in lattice dynamics could be related to the ferroelectric nature of perovskite, which is still in debate. To save the time and cost, inexpensive interatomic potential functions can be used in lattice dynamics but with a careful consideration of accuracy. Finally, defect and interface engineering are also of great importance in the meaning of direct relation to synthesis process and device architecture. Despite the progress in these respects, many important questions remain untouched, e.g. proton percolation into the interface and influence on performance.