Multi-scale model for the structure of hybrid perovskites: analysis of charge migration in disordered MAPbI₃ structures

We developed the PM concept based on our previous studies in MAPbI₃ structure optimization, using the all-electron numeric-atom-centered orbital DFT code FHI-aims [20–22]. In these studies we applied the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation [24] with Tkatchenko–Scheffler van der Waals (vdW) interactions [25] as the exchange-correlation functional. Based on our single-unit-cell structure optimization study [13], the MA⁺ cation has 32 discrete stable orientations in a unit cell—8 diagonal and 24 face-to-face with a small deviation angle to the lattice plane. In these low-energy states the ammonium end of the cation forms hydrogen bonds with the lattice I⁻ atoms, pulling them towards the cation. The amount of displacement is the largest in the face-to-face MA⁺ orientations, which subsequently has ca. 21 meV lower energy than the diagonal orientation. We have parametrized the structural energies of this framework distortion by describing the MA⁺ cation as a dipole vector in the C–N bond direction and studying the dipole configurations in neighboring unit cell pairs, the PMs (see figure 1(a)). The displacement of the I⁻ atom depends not only on the dipole orientation in a single unit cell, but also on the orientation in the neighboring cell. By considering the rotational and reflectional symmetries of the dipole pairs, we were able to reduce the amount of 32² = 1024 pairs to only 86 unique PMs. Furthermore, by approximating the face-to-face orientations in each unit cell to be identical (i.e. four-fold degenerate without the lattice-plane deviation angle), the amount of unique PMs in this 14-orientation framework (with 8 diagonal and 6 face-to-face orientations) was reduced to 25.

Zoom In Zoom Out Reset image size

In our recent study [14] we optimized 33 4 × 4 × 4 MAPbI₃ supercells (see figure 1(b)) with DFT and studied the effect of optimization on the PM distribution (see figure 2). We observed that PMs 1–17, in which one or both dipoles are diagonal, are extremely rare in the optimized structures; much rarer than in an 'intrinsic' random distribution (see red bars in figure 2). This is an obvious consequence of the higher energy of the diagonal structure. Comparing our optimized PM distribution to then only the face-to-face 'intrinsic' PM probabilities of a random system (green data in figure 2), we observe a notable increase in the frequency of PMs 20 and 22. These correspond to perpendicular-aligned dipole pairs. Accordingly, the frequency of the other PMs decrease. To study the dipole order in MAPbI₃ on a large scale and its effect on the charge migration properties, we then generate new structures based on this optimized PM distribution.

Zoom In Zoom Out Reset image size

We generated the structures using our own algorithm, in which the dipoles are reoriented one-by-one in a random order to reach the given PM distribution (the 'target' distribution, hereafter). In this study, the target distribution corresponds to the optimized PM distribution, shown in figure 2. The process is started from an initial state, in which the dipoles are oriented randomly in the 24 face-to-face orientations. These orientations correspond to the minimum-energy structures found in our previous study [13]. New dipole orientations are then determined by evaluating the change that any of the other orientations would induce in the PM distribution. By subtracting the present distribution from the target distribution, the algorithm calculates a difference value for each PM, which tells if the PM is too rare (positive value) or too frequent (negative value) in the present system. Each new dipole orientation is then evaluated by calculating the sum of the difference values for the 6 NNs. The most favorable orientation is the one that produces the largest sum, i.e. provides the largest increase in the deficient PMs. If a more favorable orientation is found, the dipole is rotated to the new orientation, after which the algorithm proceeds to the next random dipole. Via this process, the system evolves towards the target distribution, until a specified error threshold between the distributions has been reached. The error is calculated as the root mean square error (RMSE) of the PM frequency between the present and the target distribution. Due to a noisy variation in the present distribution, the threshold is evaluated with respect to the mean of the 10 previous RMSEs. In practice, the structure generation algorithm consists of an iteration loop, in which the following steps are executed.

• (i)  
  Select a random dipole.
• (ii)  
  Rotate the dipole through all 32 orientations to see, which one (if any) produces the most favorable change towards the target PM distribution.
• (iii)  
  If another orientation is more favorable, reorient the dipole to the most favorable orientation.
• (iv)  
  Calculate the present PM distribution.
• (v)  
  Calculate the RMSE between the present and target distributions. Stop if the specified accuracy is reached.
• (vi)  
  Return to (i).

In this process, we also analyzed the electrostatic dipole–dipole interaction energy of the structure. The energy is

$$\begin{eqnarray}&&U=-\displaystyle \sum _{i}\displaystyle \sum _{j\ne i}{{\boldsymbol{p}}}_{i}\cdot {{\boldsymbol{E}}}_{{ij}}=-\displaystyle \frac{1}{4\pi {\epsilon }_{0}{\epsilon }_{r}}\displaystyle \sum _{i}\displaystyle \sum _{j\ne i}{{\boldsymbol{p}}}_{i}\cdot \displaystyle \frac{3({{\boldsymbol{p}}}_{j}\cdot {\hat{{\boldsymbol{r}}}}_{{ij}}){\hat{{\boldsymbol{r}}}}_{{ij}}-{{\boldsymbol{p}}}_{j}}{| {{\boldsymbol{r}}}_{{ij}}{| }^{3}},\end{eqnarray} \tag{ 1 }$$

in which ${{\boldsymbol{p}}}_{i}$ and ${{\boldsymbol{p}}}_{j}$ are the electric dipole moments of dipoles i and $j,{{\boldsymbol{E}}}_{{ij}}$ is the electric field at dipole i by dipole $j,{{\boldsymbol{r}}}_{{ij}}$ is the distance vector from dipole j to dipole i, ₀ is the vacuum permittivity and _r is the dielectric constant. The results of the structure generation are presented in section 3.1.

We analyzed the effect of dipole order on the charge migration using a self-developed algorithm (see figures 3 and 4), in which a fully localized charge-carrier executes a classical random walk in a simulation box, based on PM-specific NN hopping probabilities. The QM behavior of the charge-carrier is evaluated at each hop based on an analytical solution of the hopping probabilities and the time step (see section 2.3.1), without the need to compute the time evolution of the QM wave-function explicitly. The analytical solution applies the PM-specific electronic coupling energies between the NNs, which are pre-computed transfer integrals (TIs) from the DFT band structure of MAPbI₃ (see section 2.3.2). Importantly, the analytical solution results in an extremely low computational simulation cost. This facilitates the calculation of a large number of migration paths in large systems. For example, in this study, each sample set consists of 10⁴ paths, which were calculated in systems of 10⁶ unit cells and above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Each hopping path was started from a random initial location at the bottom plane (z = 0) of the system (see figure 4(a)). To increase the number of the paths that migrate through the system and not exit from the sides, the initial sites were selected from the middle half of the bottom plane (in range $[d/4,3d/4]$, in which d is the dimension of the system in x or y direction). The hopping was then executed (see figure 4(b)) until the path exited the system, from any of the 6 sides of the system (see figure 4(c)). To simulate the force of an external electric field on the charge-carrier, we applied a constant bias multiplier to the hopping probabilities in the direction of the electric field, which was set to +z. The multiplier increased the number of the paths that migrated through the system, which was essential in obtaining a sufficiently large set of migration times for the statistical charge velocity. For this model, we first derived the analytical solution for the NN hopping probabilities (see section 2.3.1). The solution applies the electronic coupling energies V_n between the NNs, which we calculated from the DFT band structure of MAPbI₃ (see section 2.3.2).

2.3.1. Analytical solution for hopping probabilities

The NN hopping probabilities and the time step have an analytical QM solution (see appendix A for details). The Hamiltonian consists of the NN coupling energies V_n between the initial site 0 and the neighboring sites (i.e. the 6 NNs that share a common face of the cubic unit cell with the initial site) n ∈ [1,6], such that

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{n=1}^{6}(| 0\rangle {V}_{n}\langle n| +| n\rangle {V}_{n}\langle 0| ).\end{eqnarray} \tag{ 2 }$$

We defined the squared sum of the coupling values as

$$\begin{eqnarray}&&V={\left(\displaystyle \sum _{n=1}^{6}{V}_{n}^{2}\right)}^{1/2}.\end{eqnarray} \tag{ 3 }$$

From these, we derived the time evolution of the wave function and the probabilities for the charge to be at the initial site or at a neighboring site. The probability for the charge to be at the initial site at time t, starting from a fully localized state of the wave function, is

$$\begin{eqnarray}&&P(0;t)={\cos }^{2}\left(\displaystyle \frac{{Vt}}{{\hslash }}\right).\end{eqnarray} \tag{ 4 }$$

Importantly, the probability P(0; t) is a simple function of V and does not require solving the QM time evolution of the wave function explicitly in each hopping step, which is computationally expensive. As can be seen from the cosine function, the probability for the charge to be at the initial site starts to decrease as time progresses. Subsequently, the charge has a non-zero probability to exist at 7 different sites, i.e. at the initial site or at any of the 6 neighboring sites. Thus, we approximated that the charge will hop, on average, once its initial-site probability has decreased below 1/7. The hopping time is therefore

$$\begin{eqnarray}&&{t}_{\mathrm{hop}}=\displaystyle \frac{{\hslash }}{V}{\cos }^{-1}\left(\displaystyle \frac{1}{\sqrt{7}}\right).\end{eqnarray} \tag{ 5 }$$

At the time t_hop, the probability of the charge to be at site n is

$$\begin{eqnarray}&&P(n;{t}_{\mathrm{hop}})={\left(\displaystyle \frac{{V}_{n}}{V}\right)}^{2}{\sin }^{2}\left(\displaystyle \frac{{{Vt}}_{\mathrm{hop}}}{{\hslash }}\right).\end{eqnarray} \tag{ 6 }$$

We then applied the bias multiplier to the probability of the neighbor n in the direction of the electric field and determined the site to which the charge hopped according to the thus-obtained probability distribution. Next, we will introduce the calculation of the electronic coupling energies V_n from the DFT band structure of MAPbI₃.

2.3.2. Electronic coupling between neighboring unit cells

The DFT band structure of MAPbI₃ (see figures 6(a) and C1) exhibits a cosine-type dispersion relation, which is known from the tight-binding (TB) theory [26]. Based on this, we adopted the TB model in the calculation of the electronic coupling between the neighboring unit cells. The coupling energies correspond to the TIs, from which we calculated the NN hopping probabilities in the charge-migration model via the analytical solution.

In the TB approximation, the wave function Ψ_k of a particle in a periodic lattice, in which k is the momentum quantum number, is written as a linear combination of the atomic wave functions ϕ_n, in which n denotes the site index. In the one-dimensional single-band TB model:

• (i)  
  ϕ_n are considered to form a complete orthonormal set, i.e. $\langle {\phi }_{n}| {\phi }_{n^{\prime} }\rangle ={\delta }_{{nn}^{\prime} }$.
• (ii)  
  All sites have the same on-site energy, i.e. $\langle {\phi }_{n}| \hat{H}| {\phi }_{n}\rangle =\mathrm{constant}$. This constant can be set to 0 without loss of generality.
• (iii)  
  Only the couplings between the NNs are non-zero, which we will denote by −Λ, i.e. $\langle {\phi }_{n}| \hat{H}| {\phi }_{n^{\prime} }\rangle =-{\rm{\Lambda }}{\delta }_{{nn}^{\prime} \pm 1}$.

The one-dimensional TB Hamiltonian can thus be written as

$$\begin{eqnarray}&&\hat{H}=-{\rm{\Lambda }}\displaystyle \sum _{n=1}^{\infty }(| {\phi }_{n}\rangle \langle {\phi }_{n+1}| +| {\phi }_{n+1}\rangle \langle {\phi }_{n}| ).\end{eqnarray} \tag{ 7 }$$

In practice, we will work on finite systems, in which the total number of sites is N. For this, we will apply periodic boundary conditions by coupling the first and the last (Nth) site with −Λ, such that

$$\begin{eqnarray}&&\hat{H}=-{\rm{\Lambda }}\left[\displaystyle \sum _{n=1}^{N-1}(| {\phi }_{n}\rangle \langle {\phi }_{n+1}| +| {\phi }_{n+1}\rangle \langle {\phi }_{n}| +| {\phi }_{N}\rangle \langle {\phi }_{1}| +| {\phi }_{1}\rangle \langle {\phi }_{N}| )\right].\end{eqnarray} \tag{ 8 }$$

The solutions to this eigenvalue problem are Bloch waves Ψ_k, and the corresponding eigenvalues ε_k are cosine functions:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{k}=\displaystyle \frac{1}{\sqrt{N}}\displaystyle \sum _{n=1}^{N}\exp ({\rm{i}}{kna}){\phi }_{n}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\varepsilon }_{k}=-2{\rm{\Lambda }}\cos ({ka}),\end{eqnarray} \tag{ 10 }$$

in which a is the distance between the NN sites. With (10), we can calculate the TIs as the coupling energies −Λ from the amplitude of the conduction band in the DFT band structure of the material.

The effective mass ${m}^{* }$ of the charge-carrier can be calculated from the curvature of the conduction band (see for example [27]) as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{m}^{* }}=\displaystyle \frac{1}{{{\hslash }}^{2}}\displaystyle \frac{{{\rm{d}}}^{2}{\epsilon }_{k}}{{\rm{d}}{k}^{2}}.\end{eqnarray} \tag{ 11 }$$

Using the cosine function of the energy bands (10) we get

$$\begin{eqnarray}&&{m}^{* }=\displaystyle \frac{{{\hslash }}^{2}}{2{a}^{2}{\rm{\Lambda }}},\end{eqnarray} \tag{ 12 }$$

in which we have applied ka = 0 in the conduction-band minimum. With the effective mass, we can then calculate the mobility of the charge-carrier as

$$\begin{eqnarray}&&\mu =\displaystyle \frac{q\tau }{{m}^{* }}=\displaystyle \frac{2{a}^{2}{\rm{\Lambda }}q\tau }{{{\hslash }}^{2}},\end{eqnarray} \tag{ 13 }$$

in which q is the charge of the charge-carrier and τ is the scattering time. The processes that affect the scattering time in MAPbI₃, for example the electron–phonon coupling, are outside the scope of our charge-migration model. Thus, we calculated the charge-carrier mobility based on an experimental reference value for the scattering time in the normal-temperature tetragonal phase. We obtained the mobilities for the other analyzed structures by scaling the computed velocities in our charge-migration model linearly in relation to the tetragonal-phase mobility. The results of this analysis are presented in section 3.3.

We performed the band-structure DFT calculations with the all-electron numeric-atom-centered orbital code FHI-aims [20–23]. For all calculations we used tier 2 basis sets and the PBE exchange-correlation functional [24] augmented with long-range vdW corrections based on the Tkatchenko–Scheffler method of Hirshfeld electron-density partitioning [25]. In addition, scalar relativistic effects were included by means of the zero-order regular approximation [28]. A Γ-centered 4 × 4 × 4 k-point mesh was used for the supercell models a and b (see figure C1), and a 8 × 8 × 8 k-point mesh for the primitive-cell model c (see figure 6(a)). The results of the band-structure calculations are available from the Novel Materials Discovery (NoMaD) repository [29].

First, we studied the performance of the structure-generation algorithm by generating 50 structures of 20 × 20 × 20 unit cells (see figure 5(b)). With the PM distribution of the 33 energy-optimized MAPbI₃ 4 × 4 × 4 supercells (see figure 2) as the target distribution, we consider the generated structures to describe a low-energy MAPbI₃ structure on a large scale. In the process, we followed the electrostatic dipole–dipole interaction energy (1) within a cutoff radius of 3 unit cells (see figure 5(a)), using the experimental dielectric constant of 25.7 [30] and the dipole moment of 2.3 D [9]. Since the electrostatic dipole–dipole energy decreases with respect to the distance between the dipoles as r⁻³, we consider the dipole–dipole interactions to be negligible in this study beyond the radius of 3 unit cells (ca.19 Å). We evaluated the convergence of the solution by following the RMSE between the present and the optimized PM distributions. The moderately small system size allowed us to update the PM distribution after each dipole reorientation. The RMSE was updated once in every 100 reorientations, and the system generation was terminated once the mean value of the 10 previous RMSEs decreased below 10⁻⁷.

Zoom In Zoom Out Reset image size

The results show that the optimized PM distribution is reached well within 10⁴ generation steps (i.e. dipole reorientations). The electrostatic energy, which is initially ca. 0 in a system of randomly oriented dipoles, decreases in the generation to ca. −0.08 meV per unit cell.

The generated structures suggest that, instead of the polarized domains (see e.g. Frost et al [16]), the dipole order in MAPbI₃ on a large scale may predominantly consist of non-aligned dipoles, which can be in perpendicular formations, corresponding to PMs 20 and 22, and the out-of-plane PM 25. This resembles the tetragonal-phase MAPbI₃ structure (see for example [12]), which however is not static at normal temperatures due to dipole reorientations [8, 10]. The dynamic structure may exhibit the tetragonal-phase dipole order in local domains on a small scale, while remaining macroscopically disordered.

To obtain the NN hopping probabilities in the charge-migration model, we calculated the TIs from the DFT band structure of MAPbI₃. We analyzed the charge migration in structures, which consist of the face-to-face PMs 18–25, thus excluding the rare diagonal dipoles. We obtained the TIs for these PMs from the band structures of one primitive-cell model and two 2 × 2 × 2 supercell models. We calculated the band structure of MAPbI₃ in the primitive-cell model (see figure 6(a)) using lattice parameters $a=b=c=6.313\,\mathring{\rm A}$ [31]. In the relaxed structure, the dipole is oriented approximately in the +y direction. Here we use the usual notation ${R}\equiv \left(\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2}\right)$ for the reciprocal lattice of a simple cubic crystal system, in which the conduction-band minimum is located. In addition, we denote ${{M}}_{x}\equiv \left(0,\tfrac{1}{2},\tfrac{1}{2}\right)$, ${{M}}_{y}\equiv \left(\tfrac{1}{2},0,\tfrac{1}{2}\right)$ and ${{M}}_{z}\equiv \left(\tfrac{1}{2},\tfrac{1}{2},0\right)$. The band dispersion along R–M_x/y/z reflects the effects of the inter-site coupling along x/y/z. As shown in figure 6(a), they all exhibit the typical cosine character and can be described with the one-dimensional TB theory. We can thus calculate the TI from the amplitude of the conduction band using(10). For example, the TI of the extending PM 18 can be calculated from the primitive-cell model along R–M_y as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{y}=\displaystyle \frac{| {E}_{R}-{E}_{{M}_{y}}| }{4},\end{eqnarray} \tag{ 14 }$$

in which E_R and ${E}_{{M}_{y}}$ are the conduction band energies at points R and M_y, respectively. Also, in these points, ka = 0 and ka = π in (10), respectively. To obtain the TI for the parallel-aligned PM 23, we calculated the mean value of the TIs in the primitive-cell model along R–M_x and R–M_z.

Zoom In Zoom Out Reset image size

For the remaining PMs, we calculated the TIs from the DFT band structures of two cubic 2 × 2 × 2 supercell models, which mimic the dipole orientations in the low-temperature orthorhombic phase and in the mid-temperature tetragonal phase (see figure C1). In these supercell models, we used the lattice constant of $2a=12.626\,\mathring{\rm A}$. In both models, the conduction band minimum is located at Γ ≡ (0, 0, 0), which is different from the primitive-cell model due to band folding. Instead of the common notation for the high-symmetry points in the Brillouin zone, we denote here ${X}\equiv \left(\tfrac{1}{2},0,0\right),Y\equiv \left(0,\tfrac{1}{2},0\right)$ and $Z\equiv \left(0,0,\tfrac{1}{2}\right)$.

From the conduction band dispersion along Γ–X and Γ–Y, we calculated the TIs associated with the perpendicular PMs 20 and 22. Since both the head-perpendicular and the tail-perpendicular PMs exist in both models and cannot be singled out in the band structure, we took the mean value of the calculated TIs and assigned the same value to both PMs. Finally, the antiparallel PM 24 and the out-of-plane PM 25 correspond to the conduction band dispersion along Γ–Z in models a and b, respectively. The opposite-aligned PMs 19 and 21 are extremely rare in the energy-optimized structures (see figure 2) and difficult to obtain via structural relaxation in DFT. Therefore, they were excluded from this TI analysis and assigned a coupling value of 0 eV in the charge-migration model. The TIs for all the calculated PMs are presented in figure 6(b).

With our hopping model, we simulated the charge migration in various MAPbI₃ structures and analyzed the effect of dipole configurations on the charge mobility. We estimated the charge mobility based on the computed charge velocities, the effective mass of the charge-carrier calculated from the DFT band structure (see (12) and figure 6(a)), and an experimental reference value for the scattering time. Since we computed the velocities from the average charge-migration time through a finite-sized systems, they were dependent on the system size. Therefore, we first analyzed the velocities in systems of different sizes to determine the required size for a good approximation of a bulk material.

We generated dipole configurations of 10–130 unit cells with cubic system dimensions, which are analogous to (dis)ordered MAPbI₃ structures and known crystal phases in the coarse-grained model. The structures (see figures 7 and B1) included the following.

• (i)  
  The low-temperature orthorhombic phase in two different orientations with respect to the charge migration direction (i.e. the direction of the bias multiplier in the NN hopping probabilities).
• (ii)  
  The mid-temperature tetragonal phase, in two different orientations, as above.
• (iii)  
  Parallel and antiparallel dipoles in z and x directions, which correspond to the parallel and perpendicular orientations with respect to the charge migration direction, respectively.
• (iv)  
  Randomly oriented dipoles.
• (v)  
  Optimized dipole configurations, generated with the structure-generation method based on the optimized PM distribution (see section 2.2).

In all the charge-migration calculations, we set the bias multiplier to 5, which resulted in ca. 80% of the migration paths going through the system. Since the paths were started within the center half of the bottom plane, they rarely exited at the side-faces of the cubic system, particularly in large systems (see figure 8). Typically, a path was terminated due to exiting through the bottom plane during the first few hops in the migration. We computed 10⁴ paths in each of the structures, which based on test calculations on larger sample sets, provided a sufficient accuracy for the statistical velocity calculation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results (see figure 9(a)) show that in all the analyzed structures, a good approximation of a bulk system was obtained with a system dimension of ca. 100 unit cells or larger. We observed the highest velocities in the mid-temperature tetragonal phase in the xy direction. With respect to the charge migration in the +z direction, the tetragonal xy structure consists of the out-of-plane PMs 25, which have the strongest electronic coupling, i.e. the highest TI value (see figure 6(b)). Conversely, the tetragonal xz structures had the lowest computed velocities. This can be attributed to the low electronic coupling of the perpendicular PMs 20 and 22 in the direction of the charge migration. In this case, the highly coupled out-of-plane PMs in the x and y directions effectively deviate the charge sideways, which results in a prolonged migration time through the system and thus a low velocity. The same deviation effect can be seen in the parallel and antiparallel z structures, in which the higher coupling of the parallel PM 23 results in more deviation in the sideways direction and thus a lower velocity than in the antiparallel structure. With respect to the disordered structures, the velocities were slightly higher in the generated optimized structures than in the fully random structures. This likely results from the reduced number of the non-coupling PMs 19 and 21 in the optimized structures.

Zoom In Zoom Out Reset image size

We analyzed the effect of structural defects on the charge-migration velocity by introducing volumes of randomly oriented dipoles in the tetragonal xy structures of 100 × 100 × 100 unit cells. We studied the effect of two different types of defects: planar defects and precipitates (see figure 10(a)). In the planar defects, we introduced a system-wide plane in the middle of the system, such that the orientation of the plane (xy) is perpendicular to the charge migration direction (+z). We then varied the thickness of the plane and analyzed the charge velocity with respect to the defect volume in the system, i.e. ${\widetilde{V}}_{\mathrm{defect}}/{\widetilde{V}}_{\mathrm{system}}$. In the precipitates, we introduced 25 randomly located non-overlapping cubic volumes of randomly oriented dipoles in the system. To reduce the noise that originates from the randomly located precipitates, we generated 20 systems of each defect concentration in 0%–20% volume range, and computed the velocities as their mean value.

Zoom In Zoom Out Reset image size

The results (see figure 10(b)) show that the velocity decreases fairly linearly with defect concentration. By extrapolating the linear velocity decrease, which starts from the tetragonal xy phase velocity of ca. 7.4 Å fs⁻¹ at 0% defect concentration, the velocity approaches the fully random system with a velocity of ca. 5.8 Å fs⁻¹. The linear velocity decrease is qualitatively similar in both defect types, with the precipitates having slightly lower velocities.

Finally, we calculated the charge mobility in the analyzed structures, based on (i) the computed velocities in the largest computed structures of size 130 × 130 × 130 unit cells, (ii) the effective mass of the charge-carrier, calculated from the tetragonal-phase DFT band-structure (12), and iii) an experimental reference value of the scattering time in the tetragonal-phase MAPbI₃. Karakus et al [7] have studied the electron–phonon scattering in MAPbI₃ with time-resolved terahertz spectroscopy, measuring a scattering time of ca. 4 fs in the tetragonal phase at 300 K. We associated this value with the best-performing tetragonal xy structure in our charge-migration model, which consists of the out-of-plane PMs 25 with a TI of 0.895 eV. With these values and the lattice constant of a = 6.313Å, we estimated the best charge mobility in our model from (13) as

$$\begin{eqnarray}&&{\mu }_{\mathrm{tet}}=\displaystyle \frac{2{a}^{2}{\rm{\Lambda }}| {q}_{{\rm{e}}}| \tau }{{{\hslash }}^{2}}=65.865\,{{\rm{c}}{\rm{m}}}^{2}\,{{\rm{V}}}^{-1}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 15 }$$

in which q_e is the charge of an electron. The resulting mobility is more than twice the experimental reference mobility of ca. 27 cm² V⁻¹ s⁻¹, measured by Karakus et al [7]. The difference may arise e.g. from the calculation of the effective mass, which we have derived from the DFT band structure of the tetragonal-phase 2 × 2 × 2 supercell model with the TB approximation. In reality, the atomic structure of MAPbI₃ is likely disordered on a macroscopic scale at room temperature. Despite the simplified model in this study, the result is of the correct order of magnitude.

We obtained the mobilities for the other structures by relating the computed velocities to the tetragonal-phase mobility in (15). First, we associated the tetragonal-phase mobility μ_tet with the computed velocity of the corresponding structure v_tet. Then, we assigned mobilities for the other structures in relation

$$\begin{eqnarray}&&{\mu }_{i}={\mu }_{\mathrm{tet}}\displaystyle \frac{{v}_{i}}{{v}_{\mathrm{tet}}},\end{eqnarray} \tag{ 16 }$$

in which μ_i is the charge mobility in structure i, and v_i is its computed velocity in the charge-migration model. The results (see figure 9(b)) show that the mobilities vary in a range of ca. 50–66 cm² V⁻¹ s⁻¹. The relative difference between the best- and worst-performing structures is ca. 30%. Also, the variation of 20 cm² V⁻¹ s⁻¹ is of the order of the experimentally measured charge mobility in the reference study by Karakus et al. The result suggests that the dipole order has a significant effect on the charge mobility in MAPbI₃. Interestingly, the highest and the lowest charge mobilities arise from the same tetragonal-phase structure, the difference being the direction of the charge migration. In MAPbI₃ preparation, correct orientation of the structure is therefore important to achieve optimal PV performance.

Overall, the computation times in this study are extremely low, measured in seconds or minutes. In comparison, calculations with the conventional DFT methods can typically take hours or days. The structure-generation algorithm (see section 2.2) scales approximately as ${ \mathcal O }({n}^{2})$ with respect to the number of unit cells n. The scaling factor arises non-trivially from the way that the system size relates to reaching the target PM distribution within the set RMSE threshold. Since the PM distribution of the whole system is updated at regular intervals during the generation process, the computation time depends heavily on the system size. As an example, a system of a cubic dimension of 40 unit cells takes approximately 1 min to generate, whereas a system of a dimension of 50 unit cells requires ca. 3.5 min.

The charge-migration algorithm (see section 2.3) relies on the analytical solution for the NN hopping probabilities (see section 2.3.1). Due to this, the process is extremely rapid, and the system size affects directly only on the path length in the migration direction. In large systems, in which the overhead of initiating the process is negligible, the algorithm scales approximately linearly as ${ \mathcal O }(n)$ with respect to the number of unit cells. For a given system, the charge-migration analysis can be executed in a fraction of the time that it takes to generate the system. For example, computing 10⁵ migration paths in a system of a cubic dimension of 50 unit cells (with the optimized PM distribution) takes ca. 2 s.

The computation times presented here were calculated on a single core of a modern workstation processor at the time of this research.