Dynamical properties of organo lead-halide perovskites and their interfaces to titania: insights from Density Functional Theory.

The vibrational density of states (VDOS), electronic structure and optical properties of bulk organo lead-halide perovskites, CH3NH3PbX3 (where X = Cl, I and Br), very promising and exciting candidate materials for solar-energy applications, have been studied by means of (hybrid) Density Functional Theory (DFT), with and without spin-orbit coupling, and equilibrium Born-Oppenheimer molecular dynamics (BOMD) in the constant-volume, isothermal (NVT) ensemble at 298 K. Particular emphasis has been directed towards the detailed characterisation of optimal hybrid-DFT strategies to reproduce faithfully the band gap, band structure and optical properties vis-à-vis both experiment and more computationally demanding GW calculations (i.e., those involving the single-particle Green's function, G, and the screened Coulomb interaction, W). The VDOS was found to feature intimate coupling between the lead and halide atoms, and was dominated by acoustic phonon modes - particularly so for chlorine, suggesting this as the more effective candidate material of the considered halides. Bulk optical properties were also determined. In view of promising 'hybrid' architectures of perovskites adsorbed on titania substrates, further simulations of lead iodide in contact with titania have been performed to assess thermal stability, as well as dynamical and structural properties of these systems. It was found that lattice strain led to some atomic layers in perovskite further from the interface adopting less crystal-like structure and less pronounced phonon spectra.

Introduction

Solar cells based on methylammonium lead halide perovskite (CH3NH3PbX3 where X = Cl, I or Br) absorber materials, adsorbed usually on mesoscopic TiO2 and Al2O3 substrates, have shown very promising efficiencies of over 10% very recently [1, 2, 3, 4] - indeed, towards 15%. This suggests their potential suitability in commercial solar cells, despite the current high cost of many organic hole conductors [5]. However, optimising the entire solar cell, including the supporting mesoscale or nanoscale material, e.g., titania or alumina, and electron-conducting-layer material and open-circuit voltage, constitutes a complex multi-dimensional design challenge. Considering that this optimisation needs to be in conjunction with that of the photo-absorbing material itself (e.g., lead halide perovskite) and its nanoscale morphology, the challenge becomes formidable. A particular difficulty, or opportunity, lies in optimising hole-conducting electrolytes; solid-state candidates have shown particular promise [1,6,7].

Wang et al. synthesised lead and tin halide perovskites [8], while Kojima et al. assessed their great promise for photovoltaic cells [9]. Koh et al. have examined formamidinium-containing metal halides as an alternative for near-infrared absorption in solar cells [10]. Etgar et al. examined lead halide heterojunction cells' performance and properties [11], while Schmidt et al. have reported substantial recent progress in non-template synthesis of lead bromide nanoparticles [12]. Chen et al. have made impressive advances in photovoltaic cells’ fabrication with a vapour-assisted solution processing [13]. In all of these studies [8, 9, 10, 11, 12, 13], significant progress has been made in the reproducible synthesis of higher-quality crystals - essential for relating such experimental data to high-quality molecular simulation (a goal of the present study).

Zhao et al. have examined in some detail the crucial matter of charge transport, diffusion and recombination in perovskite solar cells [14]. Importantly, Bi et al. [15], Edri et al. [16] and Jeon et al. [17] have scrutinised various hole-conducting materials as promising electrolytes in perovskite solar cells. Arguably most important has been the focus of Edri et al. on multi-dimensional cell-design optimisation to achieve very impressive increases in open-circuit voltages [18]. All of these insightful experimental studies [14, 15, 16, 17, 18] have allowed for an increased level of process understanding of material choices; this is essential for guiding further theoretical studies in clarifying mechanistic details at the atomistic/electronic level.

Computer simulation, such as Density-Functional Theory (DFT) or Many-Body Perturbation Theory (MBPT), such as the GW approximation, has much to contribute to the structural, electronic and optical characterisation of putative organo metal halide perovskites [19]. Mosconi et al. studied lead halide (chloride, bromide and iodide) perovkites by DFT [20], assessing optimal choices of functionals. Agresti et al. have considered work functions and interfaces with DFT [21], whilst Egger has highlighted the importance of intermediate bands [22]. Gehrmann and Egger have also studied and characterised shortening of disorder potentials in perovskites [23].

Even et al. have emphasised the importance of spin-orbit coupling (SOC) [24,25], as have Lang et al. [26] – finding that SOC was needed to increase the level of (semi-quantitative) agreement with experimental quantities vis-à-vis conventional DFT. Indeed, the question of whether conventional DFT ‘disguises’ inadequacies by fortuitous cancellation of errors that SOC exposes (and does not cancel so well) is a topic of much interest - and rather lively debate [25]. Quarti et al. [27] and Park et al. [28] have computed Raman spectra and compared these to experimental data, whilst Umari et al. have gauged GW, DFT and SOC effects on replicating experimental absorption spectra of perovskites in an insightful recent study [29]. Lindblad et al. have also addressed the challenging issue of optimising the electronic structure of titania-perovskite interfaces [30], which is of great importance and relevance for practical solar-cell devices. Meanwhile, Bristow et al. have made impressive progress in parameterising a ‘classical’ (strictly, pairwise) potential model for computation of structural properties or perovskites, obtaining good agreement with available experimental measurements of various properties [31]. In terms of the important rôle of defect structures in governing perovskite behaviour, Du has used DFT to emphasise the importance of lead chemistry in lead iodide, particularly with respect to defects [32]. Buin et al. [33] and Kim et al. [34] have applied DFT to examine how intrinsic effects can have a bearing on the existence of putative, potential recombination centres - and on processing/manufacturing strategies to mitigate these. Du [32] and Giorgi et al. [35] have used DFT to consider also effective masses and how this affects transport, finding sometimes that masses can be treated to achieve better quality agreement with experimental charge-transport data. Grote and Berger have examined the effect of mechanical strain on influencing band gaps and ferroelectric poloarisation of tin- and lead-halide perovskites [36], finding that there can be effective strain-induced band-gap modulation. Giorgi et al. have consolidated earlier work on transport in ref. 36 to examine the rôle of guanidinium cations with their essentially zero dipole moment in overcoming hysteresis drawbacks in current-voltage curves - using DFT to assess thermodynamic and electronic-structure properties. They found semi-quantitative accord with experiment and confirmed that DFT is a useful tool [37]. Motta et al. have elucidated the subtle and important interplay of molecular rotations of organic cations in hybrid perovskites by van der Waals-corrected DFT, finding that electronic properties have a certain sensitivity to rotational states [38]. Leguy et al. have carried out careful calibration studies to optimise DFT performance vis-à-vis experiment for predicting optical properties [39].

Importantly, Geng et al. have used DFT-based molecular dynamics to assess the stability at finite temperature of perovskite surfaces [40], and they found a reasonable level of variability of surface-sensitive properties relative to simple geometry optimisation. In BOMD, DFT-based forces are used to propagate the system dynamically according to Newton's equations [40]. The underlying dynamical propagation itself is classical. BOMD samples systems at real, finite temperatures, rather than simple geometry optimisation to 0 K: therefore, thermal effects may be gauged on properties, and thermal motion may be gauged versus energy barriers for rotations or rearrangements in perovskites, whether in bulk or at interfaces with substrates, such as titania.

Given the importance of the nano-structured, as opposed to bulk, state for perovskites, there has understandably been less progress in the DFT and ab-initio simulation community to tackle nano-structures, given conventional DFT methods' O(N3) scaling. However, recently, Koliogiorgos et al. have made impressive progress in applying both ground-state and time-dependent DFT towards cuboid-like perovskite quantum dots (QD), studying optical spectra and band-gap properties [41]. The difficulties of DFT's O(N3) scaling in handling nano-structured systems' sizes is mirrored by challenges in their laboratory synthesis [42]. Study of the nano-structures of perovskites is both an emergent and exciting area; a good review providing an in-depth outlook on nanoscale perovskite solar cells is ref. [43]. One important aspect of realising stable nanostructured perovskite solar cells lies in handling stability with respect to moisture [44], and computational approaches to handling solubility predictions in water (e.g., Henry's Law) therefore need to be considered to afford a more expansive view of their real-world stability [45].

In any event, building on experience of computational design and treatment of solar-cell materials [46, 47, 48, 49] and on selection of appropriate functionals in perovskites (albeit not organo-lead-halide-based) [50], here, we carry out DFT studies of promising candidate perovskite materials. In so doing, this study has a number of goals:

Characterising dynamical behaviour, and assessing the importance of SOC thereon, noting recent advances in sophistication of DFT simulation of perovkites [51, 51, 52, 53].

Treating with optimal accuracy the simulation of perovskite-titania interfaces, and establishing an accurate and robust means of high-fidelity DFT simulation. Here, we apply DFT-based molecular dynamics (BOMD) to study perovskite phonon modes, which has been largely neglected to date; vibrational coupling is important for electron/hole transport.

Given interest in perovskites mounted on light-absorbing metal-oxide substrates, such as titania [1, 2, 3, 4,21], another focus is the study of the structural and dynamical properties of such interfaces.

Simulation methodology

A cubic unit cell (length ~6.49 Å) was created for both proton- and methyl-terminated forms, and the three different halide species were prepared therein. Some representative bulk-state coordinates may be found in the Supporting-Information file “bulk_perovskite_coords.txt”. For BOMD (vide infra), in the proton-terminated case, 2 × 2 × 2 and 3 × 3 × 3 supercells were used. These unit cells and supercells were used for bulk calculations (cf. Figure 1, with space group P1). We also prepared a 48-atom 2 × 2 × 1 methyl-terminated iodide ‘slab’ form of x-y box length 13.78 Å, in contact with 3 × 3 × 3 bulk rutile [48,54] supercell (with identical x-y box length 13.78 Å, given the a = b ~ 4.58 Å lattice parameters for rutile [48,54]), adsorbing additional 12 lead and iodine atoms on the interface layer between the rutile as a ‘template’ for connection to the perovskite system, leading to a system composed of 222 atoms; this constituted a ‘single-layer’ adsorbed system (cf. Figure 2a). To this singly-layered system, we added a further perovskite layer, creating a ‘double-layered’ adsorbed composite system featuring 282 atoms (cf. Figure 2b). Some representative titania-perovskite coordinates may be found in the Supporting-Information file “tio2_interface_with_lead_halide_coords.txt”.

Figure 1

(a) Schematic of bulk crystal-state bulk chloride-form 3 × 3 × 3 supercell; (b) enlargement of a unit cell.

Figure 2

Representative schematics of the (a) single- and (b) double-layer iodide/perovskite interfaces, prior to substantial relaxation and thermalisation via BOMD. Dark grey represents Ti, red is O, large light grey is Pb, yellow is N, blue is H, and medium-sized light grey is I.

As a representative example, to gauge system-energy/volume relationship in the ‘proton-terminated’ bulk case, we used CRYSTAL09 [55] with Gatti [56] basis sets and B3LYP [57,58] to estimate lowest-energy structures dues to constant-volume, cubic-box geometry optimisation; a series of such calculations for different box lengths was adopted systematically, and this is reported in the Discussion section. We also performed a series of constant-volume, ‘slab-form’ two-dimensional geometry-optimisation steps of the composite 222- and 282- atom rutile/perovskite interface in CRYSTAL09 to reduce interfacial strain, for subsequent finite-temperature study under vacuum conditions using BOMD (vide infra).

Three-dimensional periodic DFT was used for each step of BOMD in the proton-terminated cases for each bulk halide system, in addition to the composite interface system, as implemented in the Vienna Ab initio Simulation Package (VASP 5.2) [59,60]. The projected augmented wave (PAW) method was chosen to represent electron-ion interactions [60]. The exchange-correlation interactions were handled by the Perdew-Burke-Ernzerhof generalised gradient approximation (GGA) approach [61], so as to investigate other functionals in addition to B3LYP. Γ-point sampling was used for MD, given the reasonably large size of the periodic supercells (vide infra). The plane wave cut-off energy was set to 400 eV. Prior to BOMD and after molecular mechanics optimisation, the atomic positions were optimised using the GGA method until the atomic forces were smaller than 0.02 eV/Å. A Nosé-Hoover NVT ensemble using Verlet integration with a 0.5 fs timestep was used, with a relatively mild thermostat coupling period of 0.5 ps [62]. Around 35 ps of NVT BOMD was run at 298 K for 3 × 3 × 3 bulk perovskite supercells, with about 10 ps for 2 × 2 × 2 case (at 298 K); no discernible size effects were detected on the nuclear density of states (vide infra). For the singly- and doubly-perovskite-layered composite systems, relaxed via two-dimensional slab-geometry DFT in CRYSTAL09, we wished to characterise thermal (structural) stability in contact with a vacuum layer via BOMD. Unfortunately, the not-unimportant point of thermal stability under experimental, finite-temperature conditions is often neglected when using DFT or other molecular-simulation techniques to study composite interfaces, relying instead on zero-temperature geometry optimisation (neglecting the crucial rôle of entropic interactions and thermal vibrations). We placed the interfacial system in supercells with z-dimensions of ~50 and 60 Å for the respective 222- and 282- atom single- and double-perovskite-layer systems (i.e., roughly three times the z-axis ‘height’ of the composite system, allowing vacuum spacings such that periodic interactions between z-periodic replicae are minimal and not allowing ‘artificial’ structural stabilisation afforded by two-dimensional slab geometry). However, we also wished to assess the impact of a lack of vacuum spacing on such single- and double-layer composite systems, and so we also created systems with no vacuum layer whatsoever, effectively creating ‘multiple-decker sandwich’-like structures replicated periodically in the z-direction ad infinitum. We carried out BOMD for ~17 ps at 10, 50 and 298 K, but found the singly-layered systems less structurally stable at 298 K (vide infra); those at 10 and 50 K remained stable, and we used sampled configurations thereof for further electronic- and optical-property analysis (vide infra). In the case of doubly-layered composite structures, it was found that the immediately-adsorbed perovskite layer retained a higher degree of crystalline order, with the second one more prone to less order at higher temperatures (e.g., 298 K). In any event, the phonon modes and structural properties of the composite systems are explored in the Results & Discussion section, along with optical properties for the singly-layered case.

To probe dynamical properties, autocorrelation functions (ACFs) were defined for velocities for each atom therein. In general, a normalised ACF of a quantity q(t) is [62].where brackets denote ensemble-averaging. This gives an indication of the probability that a quantity is correlated with itself as a function time; once sufficient time has elapsed, this probability will have decayed from 1 initially to essentially zero [62]. The normalised velocity-ACF (VACF) of an atom (or particular groups of atoms, α), cα(t), is

The ‘power spectrum’ of a (V)ACF is the real (cosine) part of the numerical Fourier transform; the underlying frequency modes characteristic of time-variations may be gleaned therefrom [62]. In (crystalline) solids, this is often termed the ‘vibrational density of states’ (VDOS), and it describes phonon modes, which may also be probed by lattice-dynamics methods (although not adjusted for finite temperatures, often making over-simplifications in potential energy). It is also possible to investigate coupling between various quantities, e.g., velocities of lead and halide atoms by computing velocity cross-correlation functions (VCCFs), where one of the sets of atoms in Eq. (2) is replaced with another, e.g., β.

It ought to be noted that for VACF-spectra, the level of expected variability of vibrational frequencies with respect to experiment is estimated to be significant (of the order of 50 cm−1, perhaps more). This is reflected in part by statistical limitations and also the appropriateness of the functional (PBE); the lack of a more refined implementation of dispersion interactions is a drawback. An objective of this study is to probe how the halide type alters dynamical properties and phonon modes.

We performed band-structure calculations for both the proton- and methyl-terminated forms using the PBE functional in VASP, in conjunction with a 12 × 12 × 12 k-point Monkhorst-Pack grid [63]. We also investigated SOC using PBE with VASP (due to previously-remarked importance of SOC [24, 25, 26,29]) to estimate band structure and band gaps.

Results and discussion

Bulk perovskite

In the proton-terminated case, the minimum-energy box lengths for B3LYP were 5.758, 6.186, and 6.66 Å for chloride, bromide, and iodide forms, respectively. In the case of the methyl-terminated systems, the box lengths were typically around 2–3 % larger, as would be expected for their explicit inclusion. These results are in reasonable accord with respective experimental box lengths of 5.68, 5.90 and 6.33 Å [64]; it is to be expected that ‘tuning’ the level of hybrid DFT would lead to better agreement, and this approach is adopted below in the guise of optical properties (vide infra).

The vibrational density of states (from maximum-peak-normalised VACF power spectra, cf. eqn. (1)) of each set of atoms are shown in Figure 3 for the proton-terminated case. It can be seen that the lowest-lying acoustic-translational phonon peak at circa 100 cm−1 for nitrogen atoms (cf. Figure 3b) coupling with (secondary) acoustic peaks of halide atoms at ~50-150 cm−1, especially in the case of the lead chloride (cf. Figure 3d), with a further peak for bromine atoms in the lead bromide at ~100-110 cm−1. Examination of the power spectra of the corresponding VCCFs for the nitrogen and halide atoms demonstrated coupling with peaks in similar ‘rattling’ frequency ranges, underlying the importance of acoustic phonons in the crystal lattice. The observation of more comparatively enhanced coupling of nitrogen atoms' translational motion with chlorine, relative to bromine and iodine, may be rationalised in terms of chlorine being lighter in mass than the other two halides, and therefore closer in mass to nitrogen. Of particular note, however, is the significant overlap of lead atoms' translational-acoustic modes at circa 20-40 cm−1 (cf. Figure 3c) coupling with halide atoms at ~50-150 cm−1, especially in the case of the lead chloride (cf. Figure 3d). This enhanced level of acoustic coupling of both nitrogen and lead with Cl is more prominent than for Br and I - resulting in greater inter-mode energy transfer upon vibrational excitation, meaning that Cl is the more effective candidate of the halides considered. This domination of VDOS by acoustic phonon modes, highlighting intimate coupling between the halide and lead atoms, is nuanced by the presence of optic modes dominated by H–N–H bending and stretch interactions of protons in the circa 1,500 and 2,800–3,500 cm−1 range, respectively (cf. Figure 3a). Although there is an appreciable extent of translational coupling with acoustic phonons in the ~150-400 cm−1 range, due to covalent bonding to the neighbouring nitrogen atoms, the proton DOS (Figure 3a) depicts the substantial rôle of optic modes; such an ‘echo’ of H–N–H bending interactions is seen readily in the case of nitrogen atoms (cf. Figure 3c). However, it has been established by thermal-conductivity investigations that optic modes in crystalline solids are typically of lesser importance [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69], and so the translational modes tend to dominate the overall VDOS.

Figure 3

Vibrational density of states (from peak-normalised power spectra of the respective atomic-species’ VACFs) of all three proton-terminated bulk lead-halide perovskites at 298 K: (a) hydrogen atoms, (b) nitrogen atoms, (c) lead atoms, and (d) halide (I, Br or Cl) atoms. Born-Oppenheimer molecular-dynamics (BOMD) simulations were carried out under NVT conditions, and the velocity-autocorrelation functions (VACFs) were extracted therefrom (cf. eqn. (2)); Fourier transformation and peak-normalisation was then applied to extract these Density of States (DOS).

The VASP PBE-determined band gaps are provided in Table 1 for the methyl- and proton-terminated cases, respectively, whilst, using lead iodide as a representative example, band structures for the both methyl- and proton-terminated forms are provided in Figure 4 (without SOC), and in Figure 5 with SOC (for the methyl-terminated case). It was found that the band gaps were direct. The band structures for the methyl- and proton-terminated cases are broadly similar (cf. Figure 4 a and b) and in reasonable agreement with recent other DFT results [29,32]; proton termination appears to affect electronic properties relatively little. The use of SOC, however, has a rather dramatic effect of reducing the predicted band gaps substantially, in line with previously reported results [24,25,29,32,70]; indeed, it is thought that SOC is important for treating effectively systems containing heavy 6p cations, e.g., TlBr [[71], [72]]. The band gaps without SOC are in better agreement with experimental values for iodide (1.51 eV [71] and 1.52 eV [73]), bromide of about 2.25 eV [10], and 2.86 for CsPbCl3 (the nearest available to the methylammonium form) [74], albeit fortuitously so, given that GGA tends to under-predict band gaps. There is little substantial difference between the methyl- and proton-terminated forms (cf. Table 1). SOC results in lower band gaps than experiment (cf. Table 1), consistent with other studies [27,33], but refs. 28 & 34 show that the application of the HSE functional, in conjunction with SOC, has resulted in better agreement with experiment – showing that suitably-applied hybrid DFT can be very helpful. Naturally, the absolute alignment of bands is vital for photovoltaic activity, as pointed out by Butler et al. with acuity [75].

Table 1

PBE band gaps for methyl- and proton-terminated systems (eV), with results without SOC in brackets.

Termination	Cl	Br	I
Methyl	1.63 (2.41)	1.27 (2.18)	0.70 (1.80)
Proton	1.59 (2.36)	1.21 (2.10)	0.67 (1.74)

Figure 4

Band structure for bulk lead iodide with (a) CH3 form and (b) proton-terminated, both along directions: Γ (0,0,0) → M (0.5,0.5,0); Γ → Z (0,0,0.5); Γ → X (0,0.5,0); Γ → A (0.5,0.5,0.5); Γ → R (0,0.5,0.5). In both cases, the energy is set to nought for the highest occupied state, and this zero-energy denoted with a solid horizontal line. These were obtained by sampling some 298 K configurations from PBE-based BOMD, so they reflect conditions at 298 K.

Figure 5

Band structure for bulk lead iodide with CH3 form, with spin-orbit coupling. Notation identical to Figure 4. These were obtained by sampling some 298 K configurations from PBE-based BOMD, so they reflect conditions at 298 K.

Interfaces of single-layered perovskite with titania

As mentioned previously, the single-layered composite system at 298 K became unstable (whether with vacuum spacing or in ‘sandwich’ arrangement, but more so in thevacuum-spacing case). This occurred after ~8 ps when the lead-iodine intermediate layer began to break up, appearing to adopt a motif close to the PbI2 form (cf. Figure 6) – more dramatically so for the case featuring the vacuum spacing; unsurprisingly, the system did not revert to its prior configuration (similar to Figure 2) thereafter. At 10 and 50 K, the structures remained stable, and we report phonon analysis from VACFs in Figures 7 and 8. Power spectra for Ti and O atoms in the ‘slab’ are depicted in Figure 7. It can be seen that the acoustic phonon peak at circa 110-120 cm−1 for Ti atoms is comparatively damped in the ‘upper layer’ in contact with the adsorbed perovskite (and intermediate lead/iodine layer) at higher frequencies than this. Considering O, the lowest-lying acoustic peak lies at around 220-230 cm−1, and this is again damped in the upper layer in contact with perovskite. In both cases for the upper layer, the peaks are generally shifted by up to 10 cm−1 with respect to more bulk-like rutile [76]; this arises primarily due to lack of crystal coordination in the upper layer and also partly to potential coupling of certain possible vibrational modes with those of the physically adsorbed lead and iodine atoms. This is confirmed by ‘echoes’ of sub-50 cm−1 longitudinal acoustic (LA) modes present in lead and iodine in bulk perovskite (cf. Figure. 3c and d) in the spectra of Ti and O at 50 K in the interfacial system (cf. Figure 7), that are absent in more bulk-like rutile surfaces: such modes are largely absent deeper in rutile-110 surfaces, for example, in ref. 77 sampled from BOMD with the PBE functional. VCCF spectra also confirmed this sub-50 cm-1 LA mode-coupling between rutile Ti and O atoms and the adsorbed lead and iodine atoms (in the interfacial layer, especially, between the rutile and perovskite).

Figure 6

Snapshot of ‘break-up’ event of composite, singly-layered system after ~8 ps of PBE-based BOMD at 298 K (with the vacuum later in place), with lead-iodine intermediate layer began to break up, appearing to adopt a motif close to the PbI2 form. Atomic representation is the same as Figure 2.

Figure 7

Normalised power spectra of Ti and O atoms' VACFs in the rutile-titania layers of the interfacial system at 50 K (with the vacuum later in place) – sampled from PBE-based BOMD.

Figure 8

Power spectra at 50 K of (a) nitrogen and (b) lead and iodine atoms' VACFs adsorbed onto the rutile system (with the vacuum later in place) – obtained from PBE-based BOMD.

In Figure 8, the spectra for nitrogen (part a) and lead and iodine (part b) atoms are shown at 50 K in the adsorbed perovskite part of the interfacial system. They are broadly similar in basic features to those of the bulk perovskite in Figure 3, albeit with more minor peaks in the 100-200 cm−1 region redolent of Ti and O peaks (cf. Figure 7), respectively, and some higher-frequency acoustic modes up to ~400 cm−1. This coupling with Ti and O acoustic modes was confirmed by respective VCCF spectra, underscoring the contribution of host-substrate phonons towards modulating dynamical properties (and, by extension, electronic, light-absorbing properties) of perovskite.

A typical configuration for a thermalised and relaxed double-layered composite system is shown in Figure 9, for the vacuum-layer case at 298 K, following ~15 ps of BOMD. More order was evident in the immediately-adsorbed layer at 298 K, with significant adhesion of Pb and I atoms to the TiO2. However, the outer perovskite layer (further from the titania surface) exhibits less crystalline order – bordering on quasi-amorphous, and this was seen in the case of ‘vacuum’ and ‘sandwich’ structures. However, at 50 K, there was somewhat less departure from crystallinity in the outer layer, albeit with less evident structuring vis-à-vis the adsorptive layer. Given that the level of mechanical strain (for both singly- and doubly-layered composite structures) introduced by the direct placement of the perovskite onto the rutile substrate, in terms of length ‘mismatch’ (around 0.5 Å) being of the order of ~3–4%, one might expect some degree of structural distortion, especially for the single-layered structure, possibly with some subtle temperature effects moderating this (as has indeed been found in the present study). The gross instability of the single-layered structure at 298 K sits in stark contrast to the relative mechanical stability of the double-layered composite, albeit with the outer layer displaying, especially at 298 K, a comparative loss of crystalline structure. To probe these structural effects more closely at 298 K, radial distribution functions (RDFs) [62] were computed for the double-layered composite system, and a typical example at 298 K is depicted in Figure 10 for N–N and Pb–Pb, in comparison with respective bulk-structure RDFs at 298 K (not dissimilar to those using classical pairwise potentials from MD in ref. [77]. The greater crystal-like character of the adsorbed layer is evident, with more quasi-amorphous structuring evident in the layer further from the titania template. This shows that, despite the mechanical strain at the interface, the additional perovskite layer lends structural stability to the adsorbed one, allowing it to retain a more crystal-like nature.

Figure 9

Typical configuration for a thermalised and relaxed double-layered composite system (here, at 298 K and with the vacuum layer), from PBE-based BOMD; more order is evident in the immediately-adsorbed layer, with significant adhesion of Pb and I atoms to the TiO2. However, the outer perovskite layer exhibits somewhat less crystalline order. Atomic representation is the same as Figure 2.

Figure 10

Radial distribution functions in the double-layered composite system (at 298 K from PBE-based BOMD and with the vacuum layer) for both layers (first being adsorbed, second the one further from the titania substrate) for (a) N–N, (b) Pb–Pb, in comparison with respective bulk-structure RDFs at 298 K. The greater crystal-like character of the adsorbed layer is evident, with more quasi-amorphous structuring evident in the layer further from the titania template.

Conclusions

The vibrational density of states (VDOS) and electronic-structure properties of organo lead-halide perovskites, CH3NH3PbX3 (where X = Cl, I and Br), very promising and exciting candidate materials for solar-energy applications, have been studied by means of (hybrid) Density Functional Theory (DFT), with and without spin-orbit coupling, and equilibrium Born-Oppenheimer molecular dynamics (BOMD) in the NVT ensemble at 298 K. Particular emphasis has been directed towards the detailed characterisation of optimal hybrid-DFT strategies to reproduce faithfully the band gap, and band-structure properties vis-à-vis both experiment and more computationally demanding GW calculations. The VDOS was found to feature intimate coupling between the lead and halide atoms, and was dominated by acoustic phonon modes.

In terms of titania-perovskite composite systems, the additional structural stability conferred by the outer perovskite layer to the immediately-adsorbed one is a key finding as regards interfacial structural properties. Developing our understanding of these interfaces is crucial to allow for molecular simulation to become a genuinely predictive design tool for perovskite-based PV devices: certainly, the use of DFT-based MD is important to study subtle temperature and size effects in such composite systems, where ‘classical’ (pairwise) potential models may be less valid (notwithstanding the impressive efforts of refs. 32 and 78). The future deployment and tuning of hybrid functionals as an alternative to more expensive spin-orbit coupling and (especially) GW methods is sine qua non in order to allow for optimal computational tractability. Given simulation's progress in tackling bulk perovskites, this is a key new challenge – as is the accurate simulation of perovskite junctions with substrate materials, such as titania.

It must be acknowledged that the titania may appear, at first glance, to be a less suitable substrate at room temperature from the results of the present study; indeed, the lower-temperature results appear more encouraging vis-à-vis structural stability, although lower-temperature solar cells may be less energetically viable in an “energy round-trip” sense. However, lattice-mismatch effects at the titania-perovskite interface evident in the present study for the small systems examined here would be less of a problem in macroscopic, experimental interfaces – although such interfaces on titania substrates have been less reported. With careful surface design, it may be the case that molecular simulation can serve as a convenient prototyping tool for optimal substrate design for perovskites in terms of their potential solar-cell screening.

Declarations

Author contribution statement

Niall English: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

Funding statement

This work was supported by the Science Foundation Ireland – National Science Foundation of China (SFI-NSFC) bilateral funding scheme (grant number SFI/17/NSFC/5229).

Competing interest statement

The authors declare no conflict of interest.

Additional information

Samples of VASP input-coordinate files for bulk- and interfacial- systems are provided in Supporting Information (with similar coordinates used in the former case for CRYSTAL also); these were sampled from molecular dynamics in VASP. Further VASP/CRYSTAL material cannot be shared openly due to licence and copyright restrictions, but the author is open to reasonable requests to share further material with licence-holders upon contact with he.

The names of the SI files are as follows:

“bulk_perovskite_coords.txt”: sample bulk perovskite coordinates

“tio2_interface_with_lead_halide_coords.txt”: sample interfacial titania/perovskite coordinate