The Highly Accurate Interatomic Potential of CsPbBr3 Perovskite with Temperature Dependence on the Structure and Thermal Properties

CsPbBr3 perovskite has excellent optoelectronic properties and many important application prospects in solar cells, photodetectors, high-energy radiation detectors and other fields. For this kind of perovskite structure, to theoretically predict its macroscopic properties through molecular dynamic (MD) simulations, a highly accurate interatomic potential is first necessary. In this article, a new classical interatomic potential for CsPbBr3 was developed within the framework of the bond-valence (BV) theory. The optimized parameters of the BV model were calculated through first-principle and intelligent optimization algorithms. Calculated lattice parameters and elastic constants for the isobaric–isothermal ensemble (NPT) by our model are in accordance with the experimental data within a reasonable error and have a higher accuracy than the traditional Born–Mayer (BM) model. In our potential model, the temperature dependence of CsPbBr3 structural properties, such as radial distribution functions and interatomic bond lengths, was calculated. Moreover, the temperature-driven phase transition was found, and the phase transition temperature was close to the experimental value. The thermal conductivities of different crystal phases were further calculated, which agreed with the experimental data. All these comparative studies proved that the proposed atomic bond potential is highly accurate, and thus, by using this interatomic potential, the structural stability and mechanical and thermal properties of pure inorganic halide and mixed halide perovskites can be effectively predicted.


Introduction
Due to its unique photovoltaic characteristics, halide perovskite not only has broad prospects in the photovoltaic power generation industry, but also has potential technological applications in other aspects [1][2][3][4]. Hence, a growing stream of research involving halide perovskite has been both experimentally and theoretically conducted by technologists. It has been reported that the solar cell efficiency for halide perovskite has grown rapidly from 2.2% in 2006 to 25.2% at present [5][6][7].
As a superior photoelectric application material, premium lead perovskite (CsPbBr 3 ), which is a typical ABC 3 perovskite with a cubic phase (space group pm-3m), has been intensively studied. In the crystal structure of CsPbBr 3 , Pb 2+ occupies the lattice center, Cs + occupies the vertex and Br 2− is located in the plane center [8,9]. Since photoelectric properties are closely related to the phase transition and structural stability, CsPbBr 3 has been brought into research, and experiments have shown that the structural phase of CsPbBr 3 changes from an orthorhombic system to a tetragonal system at 373 K and to a cubic system at 403 K [10,11] due to the Jahn-Teller effect. In addition, given the high experimental cost and considerable computational time required, atomic simulations with little loss of accuracy provide a valid measure for exploring the internal mechanism of structural phase transformations. Therefore, to better characterize the long-term effect of temperature on the structural stability and phase transformation of CsPbBr 3 over a long period, we find that classical molecular dynamics (CMD) is more effective than density functional theory (DFT) [12][13][14]. Density functional theory (DFT) is one of the most reasonable physical theories to describe cell potential energy. DFT describes the lattice parameters of the dense solid system, and the binding energy and elastic modulus are consistent with the experimental measurements. Although the first-principles method has achieved success in studying the structure and the electronic and optical properties of CsPbBr 3 , the huge difficulties in studying the continuous temperature properties have promoted the development of atomic potentials and Hamiltonian potentials.
The Born-Mayer potential (referred to as the traditional BM potential in this paper) is the most common in molecular dynamics simulations of CsPbBr 3 perovskites [15][16][17]. The calculation results show that the traditional BM can accurately characterize the overall structural change in the unit cell during phase transition, but when the phase transition and internal PbBr 3 octahedral structure are distorted at the same time, the traditional BM potential cannot accurately describe the composite structural distortion, leading to the phase transition not being predicted. This is because traditional BM potentials only consider pairwise interactions. Therefore, there is an urgent need for the development of a more targeted potential function.
In this paper, the BV potential function model is employed to calculate the most intuitive structural transformation: the lattice parameter distortion and the length of the Pb-Br bond in the crystal cell. In addition to the structural and mechanical properties, the thermal properties (including thermal conductivity, thermal expansion coefficient and specific heat capacity) of CsPbBr 3 play a major role in industrial applications, which is a key feature of solar cell heat transfer. To characterize heat transfer and resistance, the equilibrium molecular dynamics (EMD) method is adopted to calculate the coefficient of thermal conductivity (κ), which agrees reasonably well with previous experimental results [18,19].

Potential Model
A classical interatomic potential for CsPbBr 3 is established within the framework of the bond-valence theory and the bond-valence vector-conservation principle [20][21][22]. According to the conservation rule, the total bond valence of each atom i tends to be equivalent to its initial atomic valence V i,0 in the ground-state structure. The actual valence V i of atom i is the sum of the ionic valence v ij , as Equation (1) shows. Moreover, v ij can be modeled by using the inverse power-law correlation between bond length and ionic valence given in Equation (2), where R 0,ij and C ij are Brown's empirical parameters. Equations (1) and (2) are as follows Based on Brown's bond-valence theory, the interatomic potential of the developed bond-valence model is expressed by the following equations: in which where the total energy E tot consists of non-bond potentials (E c , E r ) and bond potentials (E BV , E BVV ), respectively. Among them, E c is Coulomb energy, Coulomb constant is 14.4 eV/Å, and q i is the ionic charge. E r is the short-range repulsive Lennard-Jones energy and A ij is a short-range repulsion parameter between different types of atoms (Pb-Cs, Pb-Br, Cs-Br and Br-Br). E BV is the bond-valence energy, B i and C i are the proportional factors of energy unit and V i is the valence electron of the ith atom computed by the bond-valence theory.
E BVV is the bond-valence vector energy, representing the measure of local symmetry breaking. U i,0 is the initial key valence vector value and U i is the absolute value of the sum of v ij vectors, calculated via the following equation [23][24][25]: As shown in Figure 1, when CsPbBr 3 exhibits a quintessential perovskites cubic structure, U Pb = 0. Instead, when no structural symmetry is present, U Pb > 0. Therefore, in the BV potential model of CsPbBr 3 , the parameters we need to fit are A ij , B i , C i (i = Cs, Pb, Br), and the specific fitting method will be introduced in the next section.
where the total energy tot E consists of non-bond potentials ( c E , r E ) and bond poten- ,0 i U is the initial key valence vector value and i U is the absolute value of the sum of ij v vectors, calculated via the following equation [23][24][25]: As shown in Figure 1, when CsPbBr3 exhibits a quintessential perovskites cubic structure, Pb U = 0. Instead, when no structural symmetry is present, Pb U > 0. Therefore, in the BV potential model of CsPbBr3, the parameters we need to fit are ij Cs, Pb, Br), and the specific fitting method will be introduced in the next section.

Method of BV Potential Parameters
The fitting of the potential parameters A ij , B i , C i (i = Cs, Pb, Br) is performed by minimizing the force and energy differences between the DFT calculation and the BV model. The complete parameters fitting flow chart is shown in Figure 2. symmetry distribution; (b) Pb U > 0 when the cell structure of CsPbBr3 is asymmetric.

Method of BV Potential Parameters
The fitting of the potential parameters ij A , i B , i C ( i = Cs, Pb, Br) is performed by minimizing the force and energy differences between the DFT calculation and the BV model. The complete parameters fitting flow chart is shown in Figure 2. In order to obtain as much structural data as possible to fully fit the BV potential parameters, the vc-md calculation method (molecular dynamics method under the firstprinciples framework) is used. The initial parameters of CsPbBr3 crystal lattices originate from experimental values given in Refs. [26,27]. According to scale, the initial lattice parameters in equal proportion range from −10% to 10%, the multiple equilibrium structures calculated are selected as the initial database on the condition that cell size change and atomic movement are allowed. In addition, Quantum-Espresso (QE) is a first-principles computing software package based on DFT. With this software package, 260 CsPbBr3 structures were obtained, including orthogonal structures, tetragonal structures and basic cubic structures. The pseudopotentials in this work were generated using the Perdew-Burke-Ernzerhof (PBE) [28] version of Generalized Gradient Approximation (GGA) to estimate the exchange correlation. The convergence tests for cutoff energy and k-point mesh samplings were fulfilled one by one. A plane wave cutoff energy of 46 Ry and Brillouin zone sampled by a 12 × 12 × 12 Monkhorst-Pack k-point mesh was employed for energy and force calculations. Then, the annealing algorithm (SA) was optimized for the linear combination of the BV potential function model and DFT calculation of structural energy E and atomic force F: The following formula can obtain the atomic force from the potential energy of the BV model: Then, the iterative conditions were analyzed. To ensure the rationality of the energy and atomic force differences calculated by our BV model, we selected 30 comparative samples with different temperatures and structures for calculations. The next step is to bring the parameters into molecular dynamics calculations. If the energy and force differences between MD simulation and DFT results are large, the new structure obtained by MD In order to obtain as much structural data as possible to fully fit the BV potential parameters, the vc-md calculation method (molecular dynamics method under the firstprinciples framework) is used. The initial parameters of CsPbBr 3 crystal lattices originate from experimental values given in Refs. [26,27]. According to scale, the initial lattice parameters in equal proportion range from −10% to 10%, the multiple equilibrium structures calculated are selected as the initial database on the condition that cell size change and atomic movement are allowed. In addition, Quantum-Espresso (QE) is a first-principles computing software package based on DFT. With this software package, 260 CsPbBr 3 structures were obtained, including orthogonal structures, tetragonal structures and basic cubic structures. The pseudopotentials in this work were generated using the Perdew-Burke-Ernzerhof (PBE) [28] version of Generalized Gradient Approximation (GGA) to estimate the exchange correlation. The convergence tests for cutoff energy and k-point mesh samplings were fulfilled one by one. A plane wave cutoff energy of 46 Ry and Brillouin zone sampled by a 12 × 12 × 12 Monkhorst-Pack k-point mesh was employed for energy and force calculations. Then, the annealing algorithm (SA) was optimized for the linear combination of the BV potential function model and DFT calculation of structural energy E and atomic force F: ∑(|Ede f − E(qi, Aij, Bi, Ci)| + |Fde f − F(qi, Aij, Bi, Ci)|). The following formula can obtain the atomic force from the potential energy of the BV model: Then, the iterative conditions were analyzed. To ensure the rationality of the energy and atomic force differences calculated by our BV model, we selected 30 comparative samples with different temperatures and structures for calculations. The next step is to bring the parameters into molecular dynamics calculations. If the energy and force differences between MD simulation and DFT results are large, the new structure obtained by MD relaxation will continue to be input into the initial database, and finally, include the balanced multi-state and high-energy unbalanced structure. After several cycles, the energy and force between MD and DFT are approached. The optimized BV potential parameters are shown in Table 1.

Lattice Parameters and Elastic Constants
According to the above fitting process, the BV potential function parameters of CsPbBr 3 were obtained and their reliability was verified by comparing them with the average atomic energy difference of DFT. However, the accuracy of the description of specific structural characteristics still needs to be further verified. The Lattice parameters, as the important structural parameters of crystals, are often employed to characterize molecular system changes and cell phase transitions. Moreover, the elastic constants that express the stress-strain relationship within the stress range are employed to evaluate the mechanical stability of CsPbBr 3 . In this section, the validity of the optimized BV potential parameters was tested via molecular dynamic (MD) simulations with the lattice parameters and elastic constants fixed to theoretical and experimental values, as implemented in the Linear Atomic/Molecular Massive Parallel Simulator (LAMMPS) [29] open-source software package.
CsPbBr 3 undergoes the orthorhombic-tetragonal-cubic phase transitions at room temperature, 373 K and 403 K [27], respectively. Therefore, the lattice parameters and elastic constants of the three different structures were calculated by the MD simulation. Then, we carried out NPT ensembles using a 7 × 7 × 6 orthorhombic supercell with 5880 atoms at 300 K, a 7 × 6 × 6 tetragonal supercell with 5040 atoms at 380 K and a 6 × 6 × 6 cubic phase supercell with 4320 atoms at 600 K. The temperature was controlled by a Nose-Hoover thermostat and the pressure was maintained at 1 atm by a Parrinello-Rahman thermostat. The force field cutoff distance and time-step were set as 8 angstroms and 2 fs, respectively. There was a total of 120,000 running steps.
To compare the correctness of the lattice parameters and elastic constants calculated by MD simulation, Birch-Murnaghan's equation of state [30] and Born stability criteria [31] were employed to obtain the parameters of all three compounds by DFT method. The detailed calculation method is as follows.
Birch-Murnaghan's equation of state [30], as shown in Equation (10), was employed to evaluate the lattice parameters of CsPbBr 3 by fitting the unit-cell energy E versus the unit-cell volume V of each compound. In this equation, we have where E 0 and V 0 are the equilibrium unit-cell energy and volume, respectively. B 0 is the bulk modulus and B 0 is its pressure derivative, which ranges from 3.3 to 5.5. These calculations were performed by the DFT method that was implemented in the QE package. The convergence threshold of total energy is 1.0 × 10 −5 Ry and that on forces is 1.0 × 10 −4 Ry/Bohr for ionic minimization. To achieve an accurate convergence, the Monkhorst-Pack schemes with 12 × 12 × 8, 16 × 16 × 16 and 16 × 16 × 16 k-mesh were employed for the calculation of orthorhombic, tetragonal and cubic structures of bulk CsPbBr 3 , respectively. In the optimization procedure, the calculated lattice parameters a converted from volumes V were plotted against corresponding energies E. The fitting curves for all three compounds are plotted in Figure 3. The optimum lattice parameters correspond to the minimum energy of the unit cell, and then, lattice parameters a in the ground state were evaluated.
calculation of orthorhombic, tetragonal and cubic structures of bulk CsPbBr3, respectively. In the optimization procedure, the calculated lattice parameters a converted from volumes V were plotted against corresponding energies E. The fitting curves for all three compounds are plotted in Figure 3. The optimum lattice parameters correspond to the minimum energy of the unit cell, and then, lattice parameters a in the ground state were evaluated. On the other hand, elastic constants were obtained using ElaStic software [32], which provided a tool to figure out elastic constants from total energy and stress calculations with QE codes. In this study, we preferred total energy for elastic constant calculations. The values of elastic coefficients (C11, C12, C13, C22, C23, C33, C44, C55 and C66) need to be determined using the Projector Augmented Wave (PAW) with the Perdew-Burke-Ernzerhof (PBE) [28] approach. The essential Born stability criteria [31] are listed below for the three different structures of CsPbBr3: To evaluate the mechanical behavior of CsPbBr3, it is necessary to determine the elastic modulus under applied stress. According to the Hill approximations, the bulk modulus (B) and Shear modulus (G) can be obtained from the values of elastic coefficients, as shown in the following formulas [33,34]  On the other hand, elastic constants were obtained using ElaStic software [32], which provided a tool to figure out elastic constants from total energy and stress calculations with QE codes. In this study, we preferred total energy for elastic constant calculations. The values of elastic coefficients (C 11 , C 12 , C 13 , C 22 , C 23 , C 33 , C 44 , C 55 and C 66 ) need to be determined using the Projector Augmented Wave (PAW) with the Perdew-Burke-Ernzerhof (PBE) [28] approach. The essential Born stability criteria [31] are listed below for the three different structures of CsPbBr 3 : (i). Combined with Equations (10) and (11), the values of Young's modulus (E) are calculated as follows The comparison of the lattice parameters and elastic constants of CsPbBr 3 given by the MD simulation with experimental values, LDA, DFT calculation and the BV model potential, is shown in Table 2. It can be seen that the lattice parameters calculated by the DFT method are close to the experimental values (error less than 3%), indicating good reliability of data sources in an initial database. When CsPbBr 3 presents the cubic phase under high temperature, the lattice parameters calculated by the BV potential are in good agreement (error less than 2%) with the experimental values. When CsPbBr3 is orthogonal at room temperature, the accuracy of the BV model is higher than that of other models and the error is within 1%. On the other hand, the bulk modulus is calculated to implement the characterization of the material's ability to resist compression in Table 2, which is in good agreement with the reported value. This suggests that CsPbBr 3 is strong enough to return to its original volume after applied stress is removed. Moreover, the strength of CsPbBr 3 against plastic deformation and cracking or breakage is, respectively, demonstrated by the shear modulus and Young's modulus, which is consistent with theoretical (DFT method) and experimental results (error less than 20%). The Pugh ratio (B/G), calculated from the ratio of the bulk modulus to the shear modulus, can be employed to evaluate the ductility/brittleness, and a Pugh ratio of more than 1.75 indicates that CsPbBr 3 is ductile. Through the above analysis, the mechanical and structural stability of CsPbBr 3 at different temperatures can be proven to be excellent. Based on the analysis of the lattice parameters and elastic constants of CsPbBr 3 , we conclude that Jahn-Teller distortion is the most complex when CsPbBr 3 is in the orthogonal phase, and there is a large deviation from the characteristic of a typical Pnma space group. At this time, LDA cannot accurately reflect the distortion caused by the electronic structure level. The BV potential function model constructed for the perovskite structure is more accurate to describe the deformation and distortion of a regular octahedron, which is a major factor contributing to the phase transformation.
In order to comprehensively observe the structural distortion trend of CsPbBr 3 during phase transformation, an 11 × 11 × 11 cubic CsPbBr 3 supercell of 6655 atoms was constructed and relaxed under the action of BV potential function. The simulated environment was the NPT ensemble with a pressure of 1 atm and a temperature range of 200 K~600 K. The time step was set to 2fs and a total of 100,000 steps were run. The results after relaxation are shown in Figure 4. The equilibrium structure under multiple temperature fields was simulated, and the values of a, b and c axes were significantly different. From Figure 4, CsPbBr 3 cells present an orthonormal phase at 200 K. With the increase in temperature, the a and c axes tend to be equal, and CsPbBr 3 cells transform into a tetragonal phase at 310 K. At 380 K, the triaxial values are almost equal and the cells transform into cubic phase.
Based on the analysis of the lattice parameters and elastic constants of CsPbBr3, we conclude that Jahn-Teller distortion is the most complex when CsPbBr3 is in the orthogonal phase, and there is a large deviation from the characteristic of a typical Pnma space group. At this time, LDA cannot accurately reflect the distortion caused by the electronic structure level. The BV potential function model constructed for the perovskite structure is more accurate to describe the deformation and distortion of a regular octahedron, which is a major factor contributing to the phase transformation.
In order to comprehensively observe the structural distortion trend of CsPbBr3 during phase transformation, an 11 × 11 × 11 cubic CsPbBr3 supercell of 6655 atoms was constructed and relaxed under the action of BV potential function. The simulated environment was the NPT ensemble with a pressure of 1 atm and a temperature range of 200 K~600 K. The time step was set to 2fs and a total of 100,000 steps were run. The results after relaxation are shown in Figure 4. The equilibrium structure under multiple temperature fields was simulated, and the values of a, b and c axes were significantly different. From Figure 4, CsPbBr3 cells present an orthonormal phase at 200 K. With the increase in temperature, the a and c axes tend to be equal, and CsPbBr3 cells transform into a tetragonal phase at 310 K. At 380 K, the triaxial values are almost equal and the cells transform into cubic phase. Compared with BM potential [17], our BV potential function model is more accurate because of the failure to show the structural transition with temperature for the BM potential function. In the MD simulation, it can be seen that the phase transition temperature is underestimated. This underestimation has been observed in other previous DFT fitting simulations [37] due to systematic errors in the exchange-correlation function used for field optimization. Compared with BM potential [17], our BV potential function model is more accurate because of the failure to show the structural transition with temperature for the BM potential function. In the MD simulation, it can be seen that the phase transition temperature is underestimated. This underestimation has been observed in other previous DFT fitting simulations [37] due to systematic errors in the exchange-correlation function used for field optimization.

Bond Lengths
For CsPbBr 3 perovskite, with the increase in system temperature, the crystal structure will undergo orthogonal-tetragonal-cubic phase transformation. This process is accompanied by Jahn-Teller structural deformation and torsion of the PbBr 3 octahedron, which is the main reason for the complex phase transformation of CsPbBr 3 . In the PbBr 3 octahedron, the deformation is mainly reflected in the stretching of the Pb-Br ionic bond length along the a, b and c axes. By processing the position information of the structure atoms obtained from the relaxation, the changes in the ionic bond length in the structural distortion can be described.
The radial distribution function (RDF) is the ratio of the local density around an ion to the bulk density when determining ion coordinates [38,39]. It is employed to describe the density distribution of particles around the central particle in the multi-particle system. RDF curves are generally spectral where a peak denotes the presence of a particle and the first peak denotes the location of the nearest particle. We calculate the bond length of Pb-Br and perform RDF analysis at a cutoff distance of 8 angstroms. The BV model is used to calculate the average Pb-Br bond length in the ac plane and at the b-axis of the orthogonal Materials 2023, 16, 2043 9 of 13 structure at 300 K, in the ac plane and at the b-axis of the tetragonal structure at 380 K, and in the cubic structure at 600 K, respectively.
As can be seen from the RDF curve in Figure 5 and the Pb-Br bond length in Table 3, when CsPbBr 3 presents an orthogonal structure, the change in the Pb-Br bond length along the ac plane and b axis suggests the occurrence of phase transformation and octahedral distortion, which is in good agreement with the experimental observation results [14,39] within a 2% error. When CsPbBr 3 presents a cubic symmetric phase, the three axial Pb-Br bonds are consistent, and the error between the calculated value with BV potential function and the experimental value is less than 0.2%, indicating that the phase distortion in the expansion of the Pb-Br ionic bond is simpler during the phase transformation of CsPbBr 3 .
The radial distribution function (RDF) is the ratio of the local density around an ion to the bulk density when determining ion coordinates [38,39]. It is employed to describe the density distribution of particles around the central particle in the multi-particle system. RDF curves are generally spectral where a peak denotes the presence of a particle and the first peak denotes the location of the nearest particle. We calculate the bond length of Pb-Br and perform RDF analysis at a cutoff distance of 8 angstroms. The BV model is used to calculate the average Pb-Br bond length in the ac plane and at the b-axis of the orthogonal structure at 300 K, in the ac plane and at the b-axis of the tetragonal structure at 380 K, and in the cubic structure at 600 K, respectively.
As can be seen from the RDF curve in Figure 5 and the Pb-Br bond length in Table 3, when CsPbBr3 presents an orthogonal structure, the change in the Pb-Br bond length along the ac plane and b axis suggests the occurrence of phase transformation and octahedral distortion, which is in good agreement with the experimental observation results [14,39] within a 2% error. When CsPbBr3 presents a cubic symmetric phase, the three axial Pb-Br bonds are consistent, and the error between the calculated value with BV potential function and the experimental value is less than 0.2%, indicating that the phase distortion in the expansion of the Pb-Br ionic bond is simpler during the phase transformation of CsP-bBr3. .

Thermal Conductivity
It is important for well-designed hybrid perovskite-based energy-conversion devices [40] to understand the thermal transport in the CsPbBr 3 . The knowledge of thermodynamic properties of CsPbBr 3 has been limited by the difficulty of attaining an experimental measurement of thermal conductivity below room temperature [41,42]. Therefore, the theoretical prediction of thermal conductivity for CsPbBr 3 is quite necessary.
Molecular dynamics is often used as an effective theory to predict the thermal conductivity of materials. The commonly used methods are equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD). However, both methods require a high degree of precision, and slight deviations can make a difference in orders of magnitude. In this section, the thermal conductivity of CsPbBr 3 is calculated in LAMMPS using the EMD method to verify the applicability and accuracy of the BV potential function.
The EMD method is based on the Green-Kubo formula for thermal conductivity κ. This equation relates the transport index of the non-equilibrium structure to the fluctuations of the physical quantities of the equilibrium structure. The transport index is defined as the integral of the correlation time of the autocorrelation function, shown in Equation (13) as where κ uv and C uv denote the thermal conductivity tensor and heat current autocorrelation function (HCACF), respectively. It is assumed that CsPbBr 3 has isotropic characteristics in different phase states, and the 4 × 3 × 3 orthogonal supercell of 720 atoms, 4 × 2 × 3 tetragonal supercell of 480 atoms and 3 × 3 × 3 cubic supercell of 540 atoms are established as the calculation model. The initial equilibrium temperature was set at 300 K and the ambient pressure was set at 1 atm (NVT ensemble). In the potential field, the potential field cutoff energy is 8 Å, the time step is set as 0.04 fs and the balance steps are 150,000, with a total time of 6 ps. According to the orthogonal CsPbBr 3 model at 300 K, heat current autocorrelation function and thermal conductivity were calculated, as shown in Figure 6. Here, the red line represents the average values of ten simulated calculations. Heat current autocorrelation function flow becomes close to zero as time passes, which is consistent with the theoretical trend. The thermal conductivity curve changes with simulation time, and the value of thermal conductivity tends to be stable at last, which is 0.497 W·m −1 K −1 . The thermal conductivity of CsPbBr 3 can be experimentally measured as 0.42 ± 0.04 W·m −1 K −1 at room temperature [18,19]. By contrast, the error between simulation values obtained by BV potential and experiment values is within 30%. The BV potential function calculation results are better than the DFT calculation results [43].  As shown in Figure 7, the thermal conductivity of tetragonal phase at 380 K and cubic phase at 600 K are 0.412 W‧m −1 K −1 and 0.395 W‧m −1 K −1 , respectively. The red line represents the average values of ten simulated calculations. Unfortunately, we did not find experimental data on thermal conductivity at high temperatures for comparison. However, we compared the thermal conductivity of materials with a similar structure to CsPbBr3, shown in Table 4 (the unit is W‧m −1 K −1 ). It can be seen that the thermal conductivity of CsPbBr3 gradually decreases with increasing temperature, which is in accordance with the trends of CsPbI3 and CsSnI3. Moreover, the BV model results are closer to the experimental values than the DFT results for CsPbBr3 [43], and the BV model results are quite close to the experimental values of CsPbI3 and CsSnI3 [44] within a 10% error range. Similar to As shown in Figure 7, the thermal conductivity of tetragonal phase at 380 K and cubic phase at 600 K are 0.412 W·m −1 K −1 and 0.395 W·m −1 K −1 , respectively. The red line represents the average values of ten simulated calculations. Unfortunately, we did not find experimental data on thermal conductivity at high temperatures for comparison. However, we compared the thermal conductivity of materials with a similar structure to CsPbBr 3 , shown in Table 4 (the unit is W·m −1 K −1 ). It can be seen that the thermal conductivity of CsPbBr 3 gradually decreases with increasing temperature, which is in accordance with the trends of CsPbI 3 and CsSnI 3 . Moreover, the BV model results are closer to the experimental values than the DFT results for CsPbBr 3 [43], and the BV model results are quite close to the experimental values of CsPbI 3 and CsSnI 3 [44] within a 10% error range. Similar to DFT simulation methods, the accuracy of thermal conductivity κ in the MD simulation largely depends on the accuracy of the potential function. To obtain more accurate BV potential, there are more factors to be considered in the process of parameter fitting. For example, thermal conductivity calculations are continuously added to the database when fitting parameters. Another method is to add a thermal conductivity influence term to the BV potential model, which can be expressed as a function of temperature. The BV potential function provides a new idea for the thermal conductivity simulation of CsPbBr 3 .
phase at 600 K are 0.412 W‧m −1 K −1 and 0.395 W‧m −1 K −1 , respectively. The red line represents the average values of ten simulated calculations. Unfortunately, we did not find experimental data on thermal conductivity at high temperatures for comparison. However, we compared the thermal conductivity of materials with a similar structure to CsPbBr3, shown in Table 4 (the unit is W‧m −1 K −1 ). It can be seen that the thermal conductivity of CsPbBr3 gradually decreases with increasing temperature, which is in accordance with the trends of CsPbI3 and CsSnI3. Moreover, the BV model results are closer to the experimental values than the DFT results for CsPbBr3 [43], and the BV model results are quite close to the experimental values of CsPbI3 and CsSnI3 [44] within a 10% error range. Similar to DFT simulation methods, the accuracy of thermal conductivity κ in the MD simulation largely depends on the accuracy of the potential function. To obtain more accurate BV potential, there are more factors to be considered in the process of parameter fitting. For example, thermal conductivity calculations are continuously added to the database when fitting parameters. Another method is to add a thermal conductivity influence term to the BV potential model, which can be expressed as a function of temperature. The BV potential function provides a new idea for the thermal conductivity simulation of CsPbBr3.

Conclusions
In this article, for CsPbBr3 perovskites, the new BV potential function was developed based on I. D. Brown's bond-valence theory, in which the potential parameters were obtained by minimizing the energy and force differences between the DFT method and the BV model. The calculated lattice parameters, elastic constants, bond lengths and thermodynamic properties of CsPbBr 3 by the MD method verified the reliability and accuracy of potential parameters for the BV model through a comparison with the experimental data or DFT values. The final results and conclusions were as follows. We used the BV potential energy model to study the orthorhombic, tetragonal and cubic structures of CsPbBr 3 . Compared with LDA, the BV potential had higher accuracy in describing the lattice parameters, elastic constants and Pb-Br bond lengths, as well as their changing trends in different phase states, and the thermal conductivity calculated by the BV potential was very close to the