Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on-the-fly with Bayesian inference

TL;DR

This paper introduces an on-the-fly machine learning approach to generate force fields for molecular dynamics, enabling accurate, large-scale simulations of hybrid perovskites' phase transitions with minimal human intervention.

Contribution

The authors develop a versatile, automated machine learning scheme for force field generation during simulations, specifically applied to complex hybrid perovskites to accurately model phase transitions.

Findings

Successfully simulated entropy-driven phase transitions in hybrid perovskites.

Linked phase transition temperatures to ionic radii of involved species.

Determined the order of phase transitions using Landau theory.

Abstract

Realistic finite temperature simulations of matter are a formidable challenge for first principles methods. Long simulation times and large length scales are required, demanding years of compute time. Here we present an on-the-fly machine learning scheme that generates force fields automatically during molecular dynamics simulations. This opens up the required time and length scales, while retaining the distinctive chemical precision of first principles methods and minimizing the need for human intervention. The method is widely applicable to multi-element complex systems. We demonstrate its predictive power on the entropy driven phase transitions of hybrid perovskites, which have never been accurately described in simulations. Using machine learned potentials, isothermal-isobaric simulations give direct insight into the underlying microscopic mechanisms. Finally, we relate the phase transition temperatures of different perovskites to the radii of the involved species, and we determine the order of the transitions in Landau theory.

Equations

$$U=\sum _{i=1}^{N_{\mathrm{a}}}{U}_i=\sum _{i=1}^{N_{\mathrm{a}}}\sum _{i_{\mathrm{B}}=1}^{N_{\mathrm{B}}}{w}_{i_{\mathrm{B}}}K\left({\mathbf{X}}_i,{\mathbf{X}}_{i_{\mathrm{B}}}\right).$$

$$\mathbf{w}$$

$$\mathbf{\sigma }={\sigma }_{\mathrm{v}}^2\mathbf{I}+{\mathbf{\varphi }}^{\mathrm{T}}{\left[\mathbf{I}/{\sigma }_{\mathrm{w}}^2+{\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }/{\sigma }_{\mathrm{v}}^2\right]}^{-1}\mathbf{\varphi }.$$

Full text

Predicting the finite temperature properties of materials from first principles (FP) has always been a dream of materials scientists but it has hardly been achieved except for the simplest of solids. The main obstacle is that the required system sizes and simulation times are simply not attainable using standard FP techniques. Training force fields using machine learning (ML) techniques is an obvious solution to the problem. However, what has prevented ML from being widely applied is the construction of suitable reference structures. In conventional approachesBehler and Parrinello (2007); Bartók et al. (2010), training structures are selected using chemical intuition, FP calculations are performed for them, and machine-learned force fields (MLFF) are fitted. Later, when the user realizes that structures outside the present training set need to be included, additional structures are added, and the force field is retrained. This is a time-consuming trial and error process often taking months for a single material, and it is practically untraceable for multi-elemental complex materials. On-the-fly machine learning has been suggested as an alternative possibly reducing human interventionLi et al. (2015). The prime progress in the present work is that the predicted (Bayesian) error is used to decide whether FP calculations are required or can be bypassed. We put our generally applicable algorithm to the test by applying it to a puzzling material exhibiting very fast hydrogen dynamics, as well as very slow rotational dynamics. The on-the-fly ML allows us to predict phase diagrams with FP quality so far unprecedented efficiency.

We have chosen hybrid perovskites as a first application of our scheme, because the slow rotational dynamics of the molecules makes straightforward FP molecular dynamics exceedingly time-consuming. Furthermore, hybrid perovskites possess a huge scientific and technological potential. Methylammonium (MA)PbI3 is a promising solar cell materialHirasawa et al. (1994); Kojima et al. (2009); sol with a high charge-carrier mobilityStranks et al. (2013). Many experimental and theoretical studies have been performed on its atomic structure and dynamical propertiesWeber (1978); Onoda-Yamamuro et al. (1990); Baikie et al. (2013); Kawamura et al. (2002); Stoumpos et al. (2013); Whitfield et al. (2016); Poglitsch and Weber (1987); Wasylishen et al. (1985); Chen et al. (2015); Filippetti and Mattoni (2014); Mattoni et al. (2015); Lahnsteiner et al. (2016); Bokdam et al. (2017); Lahnsteiner et al. (2018) and they have revealed that experimentally this material exhibits two entropy-driven phase transitions from an orthorhombic to a tetragonal phase at 160 K, and from a tetragonal to a cubic phase at 330 K. Estimates of the transition temperatures from FP have not been reported to date and are elusive to be obtained using standard FP techniques alone, owing to the fact that the transitions are entropy driven. Although in the orthorhombic phase the molecules and the cage are essentially frozen, in the cubic phase the MA molecules and the cage reorient rapidly exploring a large phase space. Moreover, recent theoretical studiesLahnsteiner et al. (2018) indicate that the available semiempirical force fields are not accurate and that only few fairly expensive density functionals describe the instabilities of the cage and the interaction between the molecules and the cage with sufficient accuracyBokdam et al. (2017).

For the description of the machine learned potential energy surface, we use a variant of the Gaussian approximation potential (GAP) pioneered by Bart́ok and coworkersBartók et al. (2010). In GAP, the potential energy $U$ of a system with $N_{\mathrm{a}}$ atoms is described as a summation of local atomic potential energies $U_{i}$,

[TABLE]

Each $U_{i}$ is expressed as a linear combination of the kernel function $K\left(\mathbf{X}_{i},\mathbf{X}_{i_{\mathrm{B}}}\right)$ and weight factors $w_{i_{\mathrm{B}}}$. The kernel measures the similarity between the local configuration around atom $i$ and the reference local configuration $i_{\mathrm{B}}$. For the descriptor $\mathbf{X}_{i}$ and the kernel $K$, we adopted a variant of the Smooth Overlap Atomic Positions (SOAP)Bartók et al. (2013) (see Supplemental Materials (SM)SM ). The Eq. (1) allows to describe the energy, forces and stress tensor (EFS) for a given structure as $\bm{\phi}\mathbf{w}$. Here, $\mathbf{w}=\{w_{i_{\mathrm{B}}}\}$ and $\bm{\phi}$ is a matrix containing $K\left(\mathbf{X}_{i},\mathbf{X}_{i_{\mathrm{B}}}\right)$ and its derivatives with respect to the coordinates and lattice vectors. Similarly, the EFSs on all training structures can be summarized as $\bm{\Phi}\mathbf{w}$, where $\bm{\Phi}$ collects $\bm{\phi}$ for all training structures. The Bayesian theorem Bishop (2006) allows us to determine $\mathbf{w}$ and the uncertainty $\bm{\sigma}$ in the predicted EFS as

[TABLE]

[TABLE]

The FP data of the training structures enter in the vector $\mathbf{T}$. $\mathbf{I}$ is the identity matrix. The parameters $\sigma_{\mathrm{v}}^{2}$ and $\sigma_{\mathrm{w}}^{2}$ are determined to balance the accuracy and robustness of the MLFF using the evidence approximation Bishop (2006). The on-the-fly scheme has been integrated within the VASP code Kresse and Furthmüller (1996a, b) (section A in SMSM ).

The actual training was performed using a state-of-the-art meta-gradient corrected functionalSun et al. (2015) and running extensive FP simulations in all three experimentally known phases using $2\times 2\times 2$ unit cellsSM . To determine whether FP calculations are required, the Bayesian error Eq. (3) is used (see Fig. 1). If the error is above a certain threshold, adjusted also on the fly, FP calculations are performed as indicated by the white circles. The FP data are used to refine the force field at the MD steps shown as black dots. In this manner, if the system stays in a local minimum, most FP calculations are bypassed, but after reordering of the PbI cage or the MA molecules FP calculations are performed as illustrated in Fig. 1. The efficiency of the on-the-fly learning is demonstrated by the fact that during the training 99 % of the FP calculations are skipped reducing the computational time by almost a factor of 100 even during learning. This enables us to extensively explore the phase space involving slow molecular reorientations in the PbI3 framework occurring on a ps time scale. The generated regression model predicts energies, forces and stress tensors with near-FP quality of 2.6 meV/atom, 0.07 eV/Å and 0.82 kbar, respectively (Section C in SMSM ).

In Figure 2(a), we show the simulated lattice constants of MAPbI3 as a function of the temperature and compare them with experimentWhitfield et al. (2016). For this simulation $4\times 4\times 4$ unit cells were used (for determination of lattice constants see section D in SMSM ). The constructed force field accurately reproduces the structure of all three phases. An analysis on the tetragonal distortion, $(c-a)/c$, shown in Fig. 2(b) allows to pin-point the tetragonal to cubic phase transition temperature at 353 K. The critical exponent of 0.24 agrees well with the reported experimental results of 0.22-0.285Whitfield et al. (2016) and a theoretically expected value of $1/4$ for a tri-critical point on the basis of Landau theory. A careful free energy “umbrella sampling”Kästner (2011) analysis (details in the SMSM ) was used to determine the transition temperature of 215$\pm$10 K between the orthorhombic and tetragonal phase. Furthermore, for the orthorhombic to tetragonal transition, the change of the entropy at this phase transition is 1.3$\pm$0.6 $k_{\rm B}$ per MA molecule agreeing reasonably with the experimental value of 2.3 $k_{\rm B}$Onoda-Yamamuro et al. (1990). The theoretical results compare well with the available experimental data and are in essence only limited by the accuracy of the density functional.

Contrary to experiment, our simulations readily provide atomic-scale insight into the entropy-driven phase transitions. As a first step, we have analyzed the orientation of the molecular C-N axis. Figure 2(c) shows 3-dimensional polar plots of the probability distribution of the molecular orientation in the PbI3 framework at 200-400 K. In the orthorhombic phase at 200 K, the polar distribution exhibits two specific orientations predominantly along the $x$- or $y$-axis. In this phase the molecules are frozen, and their orientation alternates only spatially. In the tetragonal phase at 250 K, the molecules are also canting in $+z$ and $-z$ direction out of the $xy$-plane, so that eight lobes are visible. It should be noted that our previous FP MD simulations did not describe the molecular order in the tetragonal phase accurately at lower temperatures either because of the short simulation time or the fixed volumeBokdam et al. (2017). In the present simulations, at 250 K the short range molecular order is consistent with the order shown in the snapshot of Fig. 2(d) for the tetragonal phase. Between 300 and 350 K a spherical probability distribution gradually develops, indicating that the molecules continuously obtain additional rotational freedom close to the tetragonal to cubic phase transition. In the cubic phase at 400 K, the molecules realize a nearly free rotation, whereas in the tetragonal phase, the molecules exhibit hindered reorientations. Specifically, at 350 K we predict a reorientation rate of 6.4 ps (Section F in the SM SM ) agreeing well with the experimental results of 1.0-5.4 ps at 300-350 KOnoda-Yamamuro et al. (1990); Poglitsch and Weber (1987); Chen et al. (2015).

To obtain more insight on the atomic scale mechanism of the phase transitions, we performed simulations on $4\times 4\times 4$ unit cells using a slow heating rate of 0.5 K/ps. Even with state-of-the-art massively-parallel computers such a simulation would take several years using FP techniques alone. To unravel the microscopic mechanism of the phase transition, we introduce order parameters denoted as $\mathbf{O}$ (octahedron) and $\mathbf{M}$ (molecular) as sketched in Fig. 3(b). The vector $\mathbf{O}$ measures the angular correlation between adjacent PbI6 octahedra along the $x$-, $y$- and $z$-axis and resolves the ordering of the frame. $O^{x}$, for example, approaches unity when adjacent octahedra along the $x$-axis are tilted in the same direction, while it approaches zero when adjacent octahedra are tilted oppositely. Similarly, the vector $\mathbf{M}$ measures the angular correlation between adjacent MA molecules and resolves the molecular orderSM . Figure 3 (b) shows that below 220 K adjacent octahedra are tilted in the same direction along the $z$-axis, while they are oppositely rotated in the $xy$-plane. A corresponding trend appears in the molecular order parameter. This is the typical pattern of the orthorhombic phase as illustrated in Fig. 2(d). Between 220 K and 270 K, thermal fluctuations allow the molecules to reorient, and the three elements of the molecular order parameter gradually merge to 0.5, which corresponds to the 90-degree angle between neighboring molecules in the tetragonal phase. The octahedral order parameter changes rather abruptly to the tetragonal pattern after the transition in the molecular order parameter has finished (drop of $O^{y}$ in Fig. 3(b) at 270 K). These results clearly indicate that the transition of the PbI6 octahedra occurs after the reorganization of the molecules has completed, in other words, the molecules seem to inhibit the transition. The exact same trend is observed in a constant temperature MD at 220 K shown in Section D in the SMSM .

All in all, our simulations suggest that the orthorhombic to tetragonal transition is first order. The transition region in Figs. 3(a) and (b) is a result of the still fairly fast heating compared to experiment and slow dynamics of the MA molecules. The second transition from the tetragonal to the cubic phase, however, is continuous, since the lattice parameters and order parameters evolve smoothly. This also agrees with the reversibility, i.e. upon cooling the cubic to tetragonal transition is readily observed (see Fig. D1 in the SMSM ) and the hysteresis between heating and cooling is small.

The phase diagram of MAPbI3 is well known experimentally. In order to test the methodology for a less well studied material, we move to CsPbI3, where the experimentally observed phases strongly depend on the process of heating and annealing as well as sample preparation, and thus, their interpretations are still controversialMarronnier et al. (2018); Sutton et al. (2018). Figure 3(c) shows the simulated lattice constants as a function of the temperature. The force-field generation and MD simulations are described in the SMSM . Our MD simulations show two now clearly continuous phase transitions between orthorhombic and tetragonal at 403$\pm$13 K, and between tetragonal and cubic at 550$\pm$5 K. These transitions are also fully reversible and are observed both under heating (shown as colored lines) and cooling (shown as gray lines). The simulated lattice constants agree well with the experimental results reported by Marronier et al.Marronnier et al. (2018) The difference between MAPbI3 and CsPbI3 is related to the slow, extra rotational degrees of freedom of the MA molecule, which become accessible only in the tetragonal and cubic phases, whereas for Cs only lateral “rattling” movements become progressively more pronounced upon heating.

Although both perovskites exhibit similar crystallographic structures, the transition temperatures of MAPbI3 are significantly lower than those of CsPbI3. In order to better understand the origin of this difference, we further extend our on-the-fly scheme to other inorganic perovskites $ABX_{3}$ ($A$=Cs or Fr, $B$=Pb, and $X$=I, Br or Cl). Comparison of the transition temperatures with the Goldschmidt tolerance factorGoldschmidt (1927), $t=(r_{A}+r_{X})/(\sqrt{2}(r_{\rm Pb}+r_{X}))$, indicates that the different transition temperatures relate reasonably well with the size of the ions. Here, $r_{A}$, $r_{\rm Pb}$ and $r_{X}$ denote the ionic radii of cation $A$, Pb and halogen $X$, respectively. The determination of these values is described in Section G in the SMSM . Figure 3(d) illustrates that the transition temperatures of the hybrid perovskites decrease with increasing tolerance factor. This trend can be explained by the mechanism described below. The calculated tolerance factors of the inorganic perovskites indicate that the PbX6 octahedra need to tilt in order to allow the formation of multiple bonds between the cation $A$ and the halogen $X$. Without tilting, the Pb-halogen distance ($r_{\rm Pb}+r_{X}$) is too large compared to the optimal cation-halide distance ($r_{A}+r_{X}$). The prediction is consistent with the fact that all six perovskites exhibit the orthorhombic phase at low temperature. The situation, however, changes with rising temperature. The effective radius of the cation $A$ increases by thermal fluctuations (as schematically shown in Fig. 3(e)) and the tilting of the octahedra becomes unnecessary. Thus, the phase transition occurs. The necessary amount of the thermal fluctuation depends on the radii. The larger the cation $A$ or the smaller the halide $X$ radius, the smaller is the necessary fluctuation. Therefore, the transition temperatures decrease with increasing tolerance factor.

In summary, we have shown here that the combination of FP calculations with on-the-fly ML has the potential to lead to a paradigm shift in the modelling of complex materials at finite temperature. On-the-fly ML enables an exceedingly efficient sampling of structures over a large phase space with very little human intervention. Our scheme is straightforwardly applicable to complex multi-elemental materials and can be easily applied to other materials science problems such as ionic diffusion or catalytic reactions. For hybrid perovskites, we obtain excellent qualitative agreement with experiment as well as useful new insights. We observe the temperature driven transition from the orthorhombic, over the tetragonal to the cubic phase to be omnipresent in Pb-chalcogenides. However, MAPbI3 is in fact unique with a first-order orthorhombic to tetragonal phase transition, related to the unfreezing of the molecular rotational degrees of freedom.