Molecular Dynamics Simulations of Perovskites : The Effect of Potential Function Representation on Equilibrium Structural Properties

The perovskites with general formula ABX3 have been widely used as for materials with their unique properties (ferroelectric, piezoelectric, dielectric, catalytic and so on). Hybrid organolead halide perovskites are a class of semiconductors with ABX3 (X = Cl, Br, and I) structures consisting of lead cations in 6-fold coordination (B site), surrounded by an octahedron of halide anions (X site, face centered) together with the organic components in 12-fold cub octahedral coordination. These hybrid perovskites have a direct band gap, a large absorption coefficient as well as high charge carrier mobility that represent a very attractive characteristic of cost-effective solar cells. Basically, these crystals are inorganic solids of CaTiO3 type held together by bonds that are either ionic or partially ionic and partially covalent. In spite of the partially covalent character of the Ti-O bond, the system is modeled by a two-body central force interatomic potential (the form of the Vashishta and Rahman interatomic potential), which has been used successfully for many materials with a perovskite structure. In the present work using molecular dynamics (MD) simulation method we investigate the dynamical and structural behavior of CaTiO3 perovskite at normal pressure and temperature conditions. The MD calculations were performed on a system of 16,000 particles (3200Ca + 3200Ti + 96,00O), initially in an orthorhombic-Pbnm structure. The orthorhombic MD box had edges Lx = 53.4 Å, Ly = 53.4 Å and Lz = 61.12 Å, which provided a density matching the experimental value of ρ = 4 g/cm3. Starting with this structure and using proposed interatomic potentials the MD system stabilizes at room temperature in its initial configuration. The aim of the present study to explore the effect of potential function representations on structural equilibrium Kh. T. Kholmurodov et al. 111 properties for the perovskite models including hybrid halide ones outlined above. Concerning the perovskite equilibrium state we elucidate the role of potential function modification on the atomic pair correlation and structural re-organization. The details of the interatomic potential representation have to be crucially important for obtaining of correct analysis data in crystallic, liquid and amorphous phases including perovskite systems.


Introduction
The cost effective solar energy production and reproducible device performance are the important subjects in today nanotechnology research.In this respect, computer design and modeling of nanostructures aimed on developing of solar cells prototypes of greater efficiencies represent a great scientific interest.The computer methods based on atomic/molecular modeling approach may provide a lot of useful information concerning the crystal chemistry of solar cell materials, the dynamical and structural data, the thermodynamic properties and phase transitions, charge transfer and diffusion processes, and so on.The computer analysis would obviously be helpful for performing experiments with fewer resources thereby indicating the rationally modifying ways and finding the best structure design for the solar cells materials.In the present article the review of structural characterization of a number solar cell systems together with a novel molecular dynamics (MD) simulation data are presented [1]- [5].
Organic-inorganic hybrid solar cells.Recently, there are new opportunities for the development of photovoltaics.In 2012-2013 in the field of perovskite semiconductor photovoltaic cells were obtained principally new results for which it was achieved an efficiency of 16%.For today, the efficiency of organic-inorganic hybrid solar elements exceeded up 20% [6].This relatively new type of photovoltaics has remarkable properties.The roots of these photovoltaic systems lie in the Grätzel' "liquid" dye-sensitized solar cells (DSSCs).In contrast, the perovskite photovoltaic cells are solid state system with all their advantages.The best results were obtained for the synthetic organic-inorganic perovskite crystals.Organic-inorganic hybrid solar cells that combine mesoporous framework are perovskite light absorber and the electron and hole transporter.Perovskite can be formed as flexible and semi-transparent solar cells.The elements of the perovskite inferior in efficiency while the silicon single-stage and other photoelectric converters, however, have a great advantage in cost and simplicity of manufacture.Some experts are already predicting new material displacement of silicon solar cells [4]- [7].
Organo-metallic perovskite semiconductor.Organo-metallic perovskite semiconductor structures are attractive because they have high charge carrier mobility and a large diffusion length, allowing the photogenerated electrons and holes effectively to travel long distances without losing energy.As a result, one can use solar cells with more thin layers, which absorb more light and hence, generate more electricity.In this regard, organo-metallic, in particular organo-lead perovskite halide solar cells has become one of the most promising candidates for next-generation solar cells [2]- [8].
Perovskite thin film solar cells.Until recently, the perovskite thin film solar cells have been used exclusively and organic polymeric materials with a relatively low mobility of the holes.In this respect the more promising as a hole conductive materials to use more stable inorganic structures, in particular copper iodide.Hole conductors based on copper previously successfully used in solar cells sensitized with dyes and quantum dots.Copper iodide is inexpensive hole conductive material that can serve as a possible alternative to spiro-OMeTAD.This may lead to the development of inexpensive, high-performance solar cell perovskite.In the future, for such solar cells can be increased voltage and efficiency, in particular by reducing the high rate of recombination.Best efficiency solar cells perovskite already competitive with modern commercial technology at a much lower cost price [1] [5]- [8].
Further improvement of photovoltaics.Further improvement of such devices is promising the most promising avenue of research in the field of photovoltaics.It is important to carry out a systematic study and optimization of organic-inorganic perovskite structure in order to improve the photoelectric conversion efficiency based on it.These problems will involve the precise calculation methods (as molecular dynamics and related techniques) for exploring the stability effect of substitution known organo-metallic perovskites of different ions to other similar ions.Parallel on the basis of quantum-mechanical calculations the spectral properties of various perovskite systems could be predicted.These results have to bring to determination and identification of the most stable crystal structures with the necessary spectral, photovoltaics, and other properties important for photoelectric conversion.
Stability of developed new perovskite solar cells.The formation of perovskite thin film photovoltaic cells with high efficiency involves their further studies from the point of view of the radiation and thermal stabilities.One believes that for the newly developed perovskite solar cells the problems of radiation damages will be on the top of the structural, spectroscopic and other physical methods research as future trends.

The Perovskites Structure Properties and Design Aspects
Understanding the basic properties and functioning of organo-metallic perovskites solar cells represent a great scientific and technological interest due to a number of their unique properties (ferroelectric, piezoelectric, dielectric, catalytic and so on).The radiation and thermal stability of perovskite developed solar cells in this respect are the problems of a great research interest as well.Synthesis of new optimal perovskite semiconductor systems inevitably demands the studies of their spectral, structural and photovoltaic properties.The formation of new organic-inorganic perovskite solar cells of high efficiency characteristics correlate with the theoretical research investigation.So far, the development of theoretical and computational methods, as like as molecular dynamics and quantum chemical calculations, are important issues for the search of optimizing perovskite structures, studying the melting and phase transitions, diffusion and conductivity phenomena, etc., that are aimed on the innovation of new photovoltaic systems unknown for today [1] [4] [9].In the present study we have performed molecular dynamics (MD) simulation on CaTiO 3 system as a basic perovskite structure.The above mentioned system represents the most studied perovskite model.We have been aimed to compare the structural behavior of the CaTiO 3 perovskite under different intermolecular potential representations.Such analysis and comparison between different simulation approaches allow one to extend the studies on more complicated perovskite structures (as like as hybrid organolead ones) on stronger motivated basis.
The perovskites with general formula ABO 3 have been widely used for crystal materials.A typical representative of the class of modern crystal materials known as perovskites is calcium-titanate (CaTiO 3 ) The perovskite CaTiO 3 is feroelastic with an orthorhombic symmetry at room temperature (space group Pbnm) and undergoes two phases transitions at respectively T 1 ~1398 K (space group Mmcm) and T 2 ~1523 K (space group Pm3m).CaTiO 3 can be prepared by the combination of CaO and TiO 2 at temperatures >1300˚C.The sol-gel processes have been used to make a more pure substance, as well as lowering the synthesis temperature.These compounds synthesized are more compressible due to the powders from the sol-gel process as well and bring it closer to its calculated density (~4.04 g•mol −1 ).Below in Figure 1 several typical perovskite structures are shown.
In the basic structure, CaTiO 3 (Ca 2+ , Ti 4+ , O 2− ), or general representation ABX 3 , alkali atoms occupy A sites, A (Cs, Rb, CH 3 NH 3 ), Pb atoms occupy B sites, B (Pb), and halogen atoms occupy X sites, X (Cl, Br, I).In this regard, the perovskite system of Cs + Pb 2+ F − type represents the most interesting object.We emphasize the modern trend in developing of more complicate perovskites as 3 demonstrates schematic aspects of the perovskite structure models.
Basically, hybrid organolead perovskites are inorganic solids of CaTiO 3 type structure held together by the ionic or partially ionic and partially covalent bonds.In spite of the partially covalent character of the Ti-O bond, the system is modeled by a two-body central force interatomic potential, which has been used successfully for many ceramics materials of hybrid halide perovskites type.In the present work a series of the MD simulations were performed to investigate the effect of potential function representations on structural equilibrium properties for the CaTiO 3 model structure possessing similar behavior as hybrid halides and other complicated perovskites [3]- [10].

The MD Simulation and Interaction Potential
The MD calculations were performed on a system of 16,000 particles (3200Ca + 3200Ti + 9600O), initially in   Table 1 and Table 2 demonstrate the atomic position details of the perovskite CaTiO 3 structure model.The orthorhombic MD box had edges L x = 53.4Å, L y = 53.4Å and L z = 61.12Å, which provided a density matching experimental value of ρ = 4.0 g•cm −3 as in work [1].(Comparing with [1] where the MD parameters were used for a system of 10,240 particles = 2048Ca + 2048Ti + 6144O of an orthorhombic-Pbnm structure within the box L x = 43.022Å, L y = 43.494Åand L z = 61.107Å, providing the density ρ = 4 g/cm 3 ).Starting with the initial structure as described above and using the proposed interatomic potential, the MD system stabilizes at room temperature in the relaxed configuration.The MD simulation have been performed on the basis of the DL_POLY general-purpose code [11] [12].Various combinations of the integration algorithm were used   (NPT and NVT ensembles) controlling the pressure/temperature of the system with the standard termostat and barostat relaxation procedures (see, Appendices 1 and 2 for the FIELD and CONTROL files).
It should be noted that perovskite crystals are inorganic solids held together by bonds that are either ionic or partially ionic and partially covalent [1].In spite of the partially covalent character of the Ti-O bond, the system is modeled by a two-body central force interatomic potential, based on the form of the Vashishta and Rahman (VR) interatomic potential, which has been used successfully for many different systems.The total potential energy of the system is written as ( [1], see also Figure 4 [1]): The first term in the above formula is the Coulomb interaction potential between the ions Z i , Z j (in units of the electron charge |e|), r ij = r i − r j is he interatomic distance between ions i and j, and λ is the screening length for the Coulombic interactions.The second term represents the steric effects of the ions sizes, where H ij and η ij are the strength and exponent of steric repulsion, respectively (in [1] the following values were used η ij = 11 (for Ca-Ca and Ti-Ti pairs), 9 (for Ca-Ti, Ca-O and Ti-O), and 7 (for O-O)).The third term represents the charge-induced dipole interaction, to include the electronic polarizabilities of the atoms, where D ij is the strength of the charge-dipole attraction (O 2 − is a highly polarizable ion).The last is the induced dipole-dipole potential Figure 4. Two-body interaction potential as a function of distance.The total interaction potential is the sum of all two-body potentials [1].
based on the van der Waals interaction, where W ij is its strength.Parameter ξ is the screening length for chargedipole interactions, respectively.
We neglect both screening terms for the Coulombic (i.e.exp(−r/λ) = 1) and charge-dipole (ξ → 0) interactions.Also for the ions steric repulsion we consider a fixed value η ij = 12 for all interacting atomic pairs.In such approximation the above mentioned VR interatomic potential approximates the well-know Lennard-Jones (LJ) or (12 -6) potential types that are being widely used for the MD simulations of the condensed molecular systems similar to perovskites [12]- [16]: Table 3 and Table 4 show the atomic and potential parameters of the perovskite CaTiO 3 model.

Result and Discussion
Configuration snapshots.In Figure 5 and Figure 6 we present the computer generated CaTiO 3 perovskite structures at initial and relaxed (equilibrium) states.During the equilibrization the sample structure undergoes the recrystallization modification due to the cooling and melting processes.It is well know that in the single crystal the recrystallization occurs in the orthorhombic structure though the amorphous regions to be formed during the relaxation procedure.The present results obtained by present MD calculations agree with the simulation and experimental reported for the polycrystalline material similar to CaTiO 3 perovskite structure [1].Pair distribution functions.The MD simulation results for the RDF (radial distribution function) g αβ (r) are summarized in the Figures 7-10.For the comparison we also present the results reported in [1] and obtained by the VR potential.In Figure 7. the graphs of RDF (radial distribution function) g αβ (r) vs. r are presented for the ionic pair Ca-Ca: (a) 12 -6 potential and (b) VR potential [1].The comparison indicate a similar RDF behave for both 12 -6 and VR potentials.We see existing of the three largest g αβ (r) peaks for both (a) and (b) models.However, for the 12 -6 potential the position of the g αβ (r) peaks are located closer to the origin of r axis.This means that the Ca-Ca ionic pair forms in 12 -6 model relatively stronger bond than in the VR model.In Figure 8. the RDF graphs are presented for the ionic pair Ca-O.The comparison with the MD results of [1] indicate that the location of the g αβ (r) peak for both 12 -6 and VR potentials are very close to each other.However, for the 12 -6 potential we observe a very large g αβ (r) peak (the value of the g αβ (r) peak in model (a) is four times larger than in model (b)).This implies for the Ca-O a stronger ordering and ionic pair correlation for 12 -6

Conclusion
The paper is aimed on molecular-dynamics (MD) simulation of the CaTiO 3 perovskite model concerning the effect of interatomic potential function modification on the relaxed equilibrium states.The proposed earlier in [8] the Vashishta and Rahman (VR) interatomic potential has shown to successfully work for many different systems including perovskites.In [1] the VR interaction potential was used for the MD simulation of CaTiO 3 system and it was proved that the VR potential to be very effective for the describing of structural phase transitions under the temperature change; the MD results obtained in [1] were the same as in experimental observations.In

Figure 1 . 3 O
Figure 1.Several typical perovskite structures with a chemical formula ABX 3 are shown.The red spheres are X atoms (usually oxygens 2 3 O − ), the blue spheres are B-atoms (a smaller metal cation, such as Ti 4+ ), and the green spheres are the A-atoms (a larger metal cation, such as Ca 2+ ).(Refs: A public domain of the World Wide Web and https://en.wikipedia.org/wiki/Perovskite_(structure)).

Figure 3 .
Figure 3.The structure representation of more complicate perovskite system.an orthorhombic Pbnm structure.The structure was built using the data base of Institute of Experimental Mineralogy of Russian Academy of Sciences, http://database.iem.ac.ru/mincryst (Card No: 3594, PEROVSKITE Ca-TiO(3), Orthorhombic Pbnm, Z = 4; Ref.: Wyckoff R.W.G. (1963), Crystal Structures, 2, 410).The lattice parameters were as: a = 5.37 Å, alpha = 90o; b = 5.44 Å, beta = 90o; c = 7.64 Å, gamma = 90˚; unit cell volume = 223.19Å 3 ; molar volume = 33.61cm 3 /mol.Table1and Table2demonstrate the atomic position details of the perovskite CaTiO 3 structure model.The orthorhombic MD box had edges L x = 53.4Å, L y = 53.4Å and L z = 61.12Å, which provided a density matching experimental value of ρ = 4.0 g•cm −3 as in work[1].(Comparing with[1] where the MD parameters were used for a system of 10,240 particles = 2048Ca + 2048Ti + 6144O of an orthorhombic-Pbnm structure within the box L x = 43.022Å, L y = 43.494Åand L z = 61.107Å, providing the density ρ = 4 g/cm 3 ).Starting with the initial structure as described above and using the proposed interatomic potential, the MD system stabilizes at room temperature in the relaxed configuration.The MD simulation have been performed on the basis of the DL_POLY general-purpose code[11] [12].Various combinations of the integration algorithm were used

Figure 7 .
Figure 7.The RDF (radial distribution function) g αβ (r) vs. r for the ionic pair Ca-Ca from the MD simulations at T = 300 K obtained by (a) 12 -6 potential and (b) VR potential [1].

Figure 8 .
Figure 8.The RDF (radial distribution function) g αβ (r) vs. r for the ionic pair Ca-O from the MD simulations at T = 300 K obtained by (a) 12 -6 potential and (b) VR potential [1].

Figure 9 .
Figure 9.The RDF (radial distribution function) g αβ (r) vs. r for the ionic pair O-O from the MD simulations at T = 300 K obtained by (a) 12 -6 potential and (b) VR potential [1].potential comparing to the VR one.In Figure 9 the RDF graphs are shown for the ionic pair O-O.The comparison of the left and right hand-side graphs indicate that the amplitudes of the g αβ (r) for both 12 -6 and VR potentials are not so much differ from each other.Though even in (a) 12 -6 model a secondary g αβ (r) peak seem to appear which does not exist in the (b) VR poteitial model.Thus, for the O-O pair interactions negleting the Coulombic and charge-dipole screening potentail terms do not effect O-O ordering even so the O 2− is a highly polarizable ion.In Figure 10 the RDF graphs are shown for the ionic pair Ti-Ti.The results show that the Ti-Ti ordering and interaction for the (a) 12 -6 potential visibly are weakened in comparison with the (b) VR potential.

Figure 10 .
Figure 10.The RDF (radial distribution function) g αβ (r) vs. r for the ionic pair Ti-Ti from the MD simulations at T = 300 K obtained by (a) 12 -6 potential and (b) VR potential [1].
the present study work we have made modification of the VR potential, neglecting both screening terms for the Coulombic (i.e.exp(−r/λ) = 1) and charge-dipole (ξ → 0) interactions, thereby approaching the VR potential to the well-know Lennard-Jones (LJ) or(12 -6)  ones.With such modifications the MD generated configuration snapshot remain similar to existing experimental and simulation observations.Nevertheless, the behavior of the RDF (radial distribution function) g αβ (r) undergoes essential changes.For the ionic pair Ca-Ca the comparison indicate a similar RDF behave for both 12 -6 and VR potentials, showing existing of three largest g αβ (r) peaks, but for the 12 -6 potential the position of the g αβ (r) peaks are located closer to the origin of r axis.For the ionic pair Ca-O the g αβ (r) peak for both 12 -6 and VR potentials are very close to each other, but for the 12 -6 potential we observe a very large g αβ (r) peak where the value of the g αβ (r) in (12 -6) model has to be four times larger than in VR model.That just implies the Ca-O stronger ordering and ionic pair correlation for 12-6 potential comparing to the VR one.On the contrary, for the ionic pair O-O the g αβ (r) for both 12 -6 and VR potentials are not so much differ from each other, implying that for the O-O pair interactions negleting the Coulombic and charge-dipole screening potentail terms do not effect O-O ordering even so the O 2− is a highly polarizable ion.Also for the ionic pair Ti-Ti our MD results indicate that the Ti-Ti ordering and interaction for the (12 -6) potential have to weaken in comparison with the VR potential.In conclusion, the details of the atomic pair re-organization due to interatomic potential modifications elucidate on the role of potential function representation for cystallic, liquid and amorphous phases including the perovskite systems.

Table 2 .
Co-ordinates for all atomic positions.

Table 3 .
The atomic masses and charges of the perovskite system CaTiO 3 .

Table 4 .
Parameters of the two-body interaction potential used in the MD simulations.