Quantifying Thermal Disorder in Metal-Organic Frameworks: Lattice Dynamics and Molecular Dynamics Simulations of Hybrid Formate Perovskites.

Hybrid organic-inorganic materials are mechanically soft, leading to large thermoelastic effects which can affect properties such as electronic structure and ferroelectric ordering. Here we use a combination of ab initio lattice dynamics and molecular dynamics to study the finite temperature behavior of the hydrazinium and guanidinium formate perovskites, [NH2NH3][Zn(CHO2)3] and [C(NH2)3][Zn(CHO2)3]. Thermal displacement parameters and ellipsoids computed from the phonons and from molecular dynamics trajectories are found to be in good agreement. The hydrazinium compound is ferroelectric at low temperatures, with a calculated spontaneous polarization of 2.6 μC cm-2, but the thermal movement of the cation leads to variations in the instantaneous polarization and eventually breakdown of the ferroelectric order. Contrary to this the guanidinium cation is found to be stationary at all temperatures; however, the movement of the cage atoms leads to variations in the electronic structure and a renormalization in the bandgap from 6.29 eV at 0 K to an average of 5.96 eV at 300 K. We conclude that accounting for temperature is necessary for quantitative modeling of the physical properties of metal-organic frameworks.

Introduction

Hybrid organic–inorganic materials constitute a diverse class of materials, due to the many possible combinations of metal and organic ligand. Current applications include gas adsorption[1] and catalysis,[2] typically possible in porous systems, while more dense crystals are investigated due to properties such as magnetism[3] and ferroelectricity.[4,5]

A class of hybrid materials that have received significant attention as potential multiferroic materials are the amine-templated metal formate networks,[6] in particular those with the perovskite composition ABX3. Here A is a small amine cation, B is a metal atom with charge +2, and X is the formate anion HCO2–. The metal ions are linked by the formate groups, allowing the mediation of magnetic interactions; however, the coupling is weak with typical ordering temperatures around 10 K.[7−12] The ferroelectricity on the other hand is associated with an ordering of the cations, residing in the cavities of the metal formate cage, which can persist at and above room temperature.[8,13,14] The exact combination of metal and molecular cation can lead to the presence of one or both of these properties, and though they arise from different parts of the framework, they can be coupled.[15−19] In addition some of the formate perovskites have interesting physical properties such as negative thermal expansion[8,20] and large dielectric constants.[21]

The formate perovskites are more flexible than typical inorganic ferroelectric and magnetic materials, and understanding how this flexibility affects their properties is important. Therefore, we study the formate perovskites based on the hydrazinium cation, [NH2NH3][Zn(CHO2)3][8] (FP-Hy), and the guanidinium cation, [C(NH2)3][Zn(CHO2)3][9] (FP-Gua), the structures of which are shown in Figure . The Gua+ ion does not have a dipole moment in itself and is placed in the middle of the cage, meaning that FP-Gua is not ferroelectric; however, it has a relatively high thermal stability with the first decomposition step occurring at 503 K.[9] FP-Hy on the other hand is ferroelectric below 352 K but decomposes at slightly higher temperatures (380 K). We use a combination of ab initio lattice dynamics (LD) and molecular dynamics (MD) to investigate the finite-temperature behavior of the two structures and quantify the displacements of the atoms from the equilibrium positions. Furthermore, we investigate how the dynamics affects interesting properties such as the bandgap and the polarization of the materials. While both materials can be made with different metals, we choose the Zn variant because it is nonmagnetic; however, many of the results are expected to carry over to the magnetic variants of these structures.

Figure 1

Crystal structures of [C(NH2)3][Zn(CHO2)3] (FP-Gua) in the (a) ac-plane and (b) ab-plane and [NH2NH3][Zn(CHO2)3] (FP-Hy) in the (c) ab-plane and (d) ac-plane. Oxygen is red, carbon gray, hydrogen white, nitrogen blue, and zinc purple.

Computational Details

Density Functional Theory Calculations

DFT calculations were performed using the Vienna ab initio simulation package (VASP)[22] with PAW frozen-core pseudopotentials. The PBEsol[23] functional was used for the exchange-correlation energy with the D3[24] correction and Becke–Johnson damping[25] to account for dispersive interactions. For the optimization of the equilibrium unit cells and the calculation of phonons, a plane wave cutoff of 700 eV and a 2 × 2 × 2 k-point mesh was used, and the structures were relaxed until all forces were below 0.01 eV Å–1. For the polarization calculations a plane wave cutoff of 500 eV was found to be sufficient.

Optimization of the high-temperature structure of FP-Hy (Pnma spacegroup) is complicated as it appears as a result of dynamics. In ref (8) it was shown that the change to the Pnma structure is related to libration of the cation, which would result in it having a larger effective volume. Such libration would not be accounted for in an athermal DFT calculation, and it thus seems likely that the unit cell size would be underestimated. Furthermore, for the experimental structure to fit in the Pnma space group the molecule must occupy each of the two equilibrium positions with a probability of 0.5. In the DFT structure optimization we must choose one of these positions for each molecule. We choose the ferroelectric ordering of the molecules which has Pna21 symmetry. To prevent the structure from returning to the low-temperature structure during optimization, the ratio between the unit cell vectors was fixed and the symmetry was reduced. Optimization with these settings results in a reorientation of the molecules such that the final structure is ferrielectric, belonging to the P21 space group, but with all angles equal to 90°.

Formate perovskites have many soft degrees of freedom, and therefore phonon calculations are performed for the optimized structures to ensure that they are true local minima. For the FP-Gua structure negative frequencies were found. To remove these the structure is displaced along the corresponding phonon eigenvectors and reoptimized, repeating until no negative frequencies remain. This procedure leads to small deviations from the Pnna space group (the structure is optimized as P1), mainly for the carbons and hydrogens in the formate groups along the b-direction. Thus, all other atoms are still within the Pnna space group for a tolerance of 0.1 Å, while a tolerance of 0.8 Å is required for the entire structure to be Pnna. We will bear this in mind when the structure is analyzed within this space group in the following.

Lattice Dynamics

The phonon spectra of the equilibrium structures were calculated with Phonopy,[26] using the finite displacement method as implemented in VASP. As phonon calculations require a tight convergence to minimize noise in the force constants, we converge the self-consistent energy to 10–8 eV. The vibrations are calculated in the 1 × 1 × 1 unit cell, and thus only the frequencies at the Γ-point are considered.

The phonon spectrum can be used to calculate the thermal ellipsoids as a function of temperature in the harmonic approximation. This is done by first calculating the mean square displacement matrix Uwhere κ is the index of the atom considered which has mass mκ and N is the number of unit cells. q is the wavevector, and eκ is the eigenvector of the vibrational mode j with frequency ω. n is the phonon population at a given temperature, T:The relative lengths of the ellipsoid axes are then obtained as the square root of the eigenvalues of U(κ), and their directions are defined by the eigenvectors. Assuming a normal distribution, the axes are scaled by a factor (from a χ2 distribution with 3 degrees of freedom) so that the ellipsoids contain 50% of the data.

By performing phonon calculations for the structure at a number of volumes close to the DFT optimized volume, it is possible to calculate the thermal expansion in the quasi-harmonic approximation.[26,27] The volume at a given temperature (T) is found as the volume that minimizes the Helmholtz free energy at this temperature to give the Gibbs free energy G(T, P):Here U(V) is the electronic energy calculated by DFT and Fphonon(T, V) is the phonon energy at the given temperature and volume. The phonons are still considered to be harmonic in this approach, but thermal expansion arises from the volume dependency of the phonon frequencies. We have calculated the thermal expansion of FP-Hy by optimizing the structure at eight different volumes and fitting the Helmholtz free energy curves at different temperatures to the Vinet equation of state[28] to obtain G(T) and V(T).

Calculations of the spontaneous polarization density, Ps, were performed using the Berry-phase formalism.[29] An important factor in such calculations is the choice of a reference state, as only differences in polarization are physically meaningful. This is because the polarization for a given structure is not unique but can take on a set of values differing by the so-called quantum of polarization, corresponding to the displacement of one electron (or ion) by one unit cell vector. To ensure that we calculate the correct polarization difference between two structures, we need to verify that no such displacement occurs, i.e., that the change in polarization is continuous when one structure is transformed into the other.[30] Here we optimize the structure with polarization +Ps and invert this structure to get the opposite polarization, −Ps. The polarization difference between those two structures, equal to 2Ps, is calculated, and we verify that the change in polarization is continuous by also calculating the polarization for a number of configurations connecting the two structures, i.e., with their coordinates r obtained aswhere the parameter λ is a number between 0 and 1.

For bandgap calculations the HSE06[31,32] hybrid functional was used with a cutoff energy of 500 eV. Γ-point calculations were found to be sufficient for convergence of the gap to within 0.02 eV, as tested for both an equilibrium and a nonequilibrium structure.

Molecular Dynamics

Ab initio molecular dynamics simulations (MD) were performed in the NVT ensemble with a Nosé thermostat,[33] using a cutoff energy of 500 eV. The simulations were performed in the 1 × 1 × 1 unit cell which contains a total of 80 and 92 atoms for FP-Hy and FP-Gua, respectively, using 2 × 2 × 2 k-points. While this unit cell is minimal, containing only four molecules, an investigation of the cross-correlation between the Hy+ molecules within the unit cell shows that there is little correlation even between nearest neighbors. We thus conclude that this unit cell size is sufficient for our investigations, which are mainly focused on the local molecular motion. For studies of long-range correlated effects, a larger unit cell would however be desirable.

The simulations were run at 150 and 300 K for both structures and also at 450 K for FP-Gua. The FP-Hy structure is not stable at 450 K, but a simulation is run at 375 K, which is between the experimentally determined ferroelectric transition temperature of 352 K and the decomposition temperature of 380 K.[8] At each temperature 3 ps of equilibration simulation were run before trajectories of 10 ps were recorded for the statistical analysis. A time step of 0.5 fs was used, and the configuration was recorded every 5 fs, resulting in a total of 2000 configurations for analysis of the trajectory.

The mean square displacement matrix corresponding to the matrix in eq is readily calculated from the trajectory, with matrix elements U obtained aswhere i and j refer to different coordinate axes and is the average position along the ith axis. The length and direction of the thermal ellipsoid axes are then found following the procedure described for the phonon derived matrix.

To get better data for the statistics, the atoms that are related by symmetry are grouped together and used to calculate the matrix elements. The 2000 configurations thus give rise to at least 8000 data points used for each thermal ellipsoid.

Results and discussion

Table shows the experimental and the DFT optimized lattice constants for the FP-Gua and the FP-Hy structures. The calculated values are generally smaller than the experimental values; however, thermal expansion is neglected in these calculations, and therefore perfect agreement with the experimental structures, which are recorded at 293 K (FP-Gua) and 110 and 375 K (FP-Hy), is not expected. Indeed, if the thermal expansion as calculated in the quasi-harmonic approximation is included the agreement with experiment becomes almost perfect (c.f. section 3.3).

Table 1

Experimental Unit Cell Parameters Recorded at 293 K (FP-Gua), 110 K (FP-Hy, Pna21), and 375 K (FP-Hy, Pnma) Together with the Lattice Constants Obtained from DFT (PBEsol+D3)

	exp. FP-Gua[9]	DFT FP-Gua	exp. FP-Hy[8] <352 K	DFT FP-Hy	exp. FP-Hy[8] >352 K	DFT FP-Hy
space group	Pnna	Pnnaa	Pna21	Pna21	Pnma	P21b
a (Å)	8.349	8.17	8.664	8.57	8.596	8.50
b (Å)	11.728	11.51	7.716	7.63	11.644	11.52
c (Å)	8.909	8.87	11.482	11.34	7.847	7.76
V (Å3)	872.34	833.43	767.58	742.20	785.40	759.67

The relaxed structure deviates slightly from the Pnna space group.

The Pnma space group requires partial occupancy of the molecules. We choose the ferroelectric combination of molecular orientations (Pna21) but reduce the symmetry to P21 and fix the ratio between the unit cell vectors (see section 2.1).

Lattice Vibrations

The phonon densities of states for the optimized structures are shown in Figure a for FP-Gua and Figure b for FP-Hy. The partial density of states (PDOS) is similar for the two structures, with a low-frequency band below 500 cm–1 containing modes involving the Zn atoms, a band from 500 cm–1 to 1700 –1 mainly consisting of vibrations within the formate units and the cations, and a high-frequency band around 3000 –1 consisting of the C–H and N–H stretch frequencies. Animations of selected modes from each of the three bands can be found in the Supplementary Information.

Figure 2

Phonon partial density of states for (a) FP-Gua and (b) FP-Hy. Oxygen is red, carbon gray, hydrogen cyan, nitrogen blue, and zinc purple. The atomic contributions are assigned according to the sign contribution to each phonon eigenvector.

Thermal Displacements

The phonon calculations are used to extract the harmonic thermal displacements as described in the Computational Details (eqs –2). The dimensions of the thermal ellipsoid along the longest (Llong), medium (Lmed), and shortest (Lshort) axis are given for the symmetry inequivalent atoms in Tables and 3 for the FP-Gua and FP-Hy structures, respectively. We furthermore give the asphericity parameter calculated asThe corresponding values calculated from the MD trajectory are also given. All values are for a temperature of 300 K, while displacements at 150 and 375/450 K are given in the Supplementary Information.

Table 2

Dimensions of the Thermal Ellipsoids Calculated from Harmonic Lattice Dynamics (LD) and Molecular Dynamics (MD) Simulation for the FP-Gua Structure at 300 Ka

	Llong LD	Llong MD	Lmed LD	Lmed MD	Lshort LD	Lshort MD	b LD	b MD
Zn	14	13	13	12	12	11	0.4	0.1
O1	26	25	15	16	13	11	–4.6	–2.0
O2	27	26	17	18	13	11	–2.5	–0.7
O3	24	22	18	18	13	12	–0.2	1.1
C1	29	29	14	15	13	11	–7.6	–5.4
C2	25	22	17	17	14	13	–2.7	–0.8
H1	74	74	25	20	17	19	–21	–26
H2	67	56	26	22	20	18	–17	–15
molecule
C11	23	18	16	16	15	15	–2.7	–0.5
H11	36	32	26	22	20	18	–1.9	–2.6
H12	35	30	26	24	20	16	–1.4	0.8
H13	34	31	25	22	20	18	–2.2	–2.3
N1	37	29	17	18	17	15	–9.7	–3.6
N2	28	26	18	18	15	17	–3.0	–3.0

The longest dimension (Llong) is listed first, followed by the medium (Lmed) and the shortest dimension (Lshort). Also given is the asphericity, b, as calculated from eq . All values are in units of 10–2 Å.

Table 3

Dimensions of the Thermal Ellipsoids Calculated from Harmonic Lattice Dynamics (LD) and Molecular Dynamics Simulation for the FP-Hy Structure at 300 Ka

	Llong LD	Llong MD	Lmed LD	Lmed MD	Lshort LD	Lshort MD	b LD	b MD
Zn	16	14	13	14	12	13	–1.4	0.2
O1	27	27	18	19	14	13	–2.4	–0.5
O2	24	28	18	18	13	13	–1.1	–2.0
O3	26	29	16	17	14	13	–3.9	–4.0
O4	23	25	16	17	14	13	–2.7	–1.5
O5	24	25	16	18	12	13	–2.2	–1.9
O6	29	30	17	20	12	14	–3.3	–1.9
C1	27	27	17	17	14	13	–3.9	–3.0
C2	31	28	15	15	14	13	–8.3	–5.5
C3	24	24	16	18	13	13	–2.7	–1.0
H1	55	58	25	22	19	19	–13	–17
H2	64	63	25	21	19	20	–16	–21
H3	76	64	26	21	18	18	–21	–20
molecule
N1	22	25	18	19	17	18	–1.5	–2.4
N2	29	34	24	28	20	23	–1.1	0.2
H11	42	62	35	54	22	30	2.9	7.5
H12	40	58	33	50	23	32	0.8	5.0
H13	29	30	27	27	22	22	1.7	0.6
H14	31	31	27	31	21	21	0.6	4.6
H15	31	36	29	29	21	19	2.6	1.4

The longest dimension (Llong) is listed first, followed by the medium (Lmed) and the shortest dimension (Lshort). Also given is the asphericity, b, as calculated from eq . All values are in units of 10–2 Å.

There is generally good agreement between the values obtained by the two different methods, which is reassuring given the underlying differences in the two approaches. The LD method assumes a harmonic potential, while the MD simulations should also account for anharmonic effects. In addition the mode occupations differ, as the modes are populated according to Bose–Einstein statistics in the LD method, while the MD simulation follows a Boltzmann distribution. Comparison with the thermal displacement parameters for the C, N, O, and Zn atoms given in the single crystal X-ray diffraction (SC-XRD) data of refs (8) and (9) show that the experimental ellipsoids obtained at room temperature are generally larger than the calculated ones, with calculated volumes typically 60%(Hy)/70%(Gua) of the experimentally determined volumes. The agreement is reasonable, and differences can be attributed to finite size effects (limited supercell size and simulation time) and/or limitations in the description of electron exchange and correlation at the PBEsol level. No thermal displacement data is currently available for the hydrogen atoms.

The thermal ellipsoids of the FP-Gua cage atoms are plotted at different temperatures in Figure (the cage of FP-Hy looks very similar and is shown in the Supplementary Information). For clarity the molecular positions are plotted separately in Figure for Gua+ and Figure for Hy+.

Figure 3

Thermal positions of the FP-Gua cage atoms at (a) 150 K, (b) 300 K, and (c) 450 K as obtained from MD simulations. The zinc ellipsoids are colored purple, oxygen ellipsoids red, carbon ellipsoids gray, and hydrogen ellipsoids white.

Figure 4

Thermal ellipsoids of the Gua+ position at (a) 150 K, (b) 300 K, and (c) 450 K as obtained from MD simulations.

Figure 5

Thermal ellipsoids of the Hy+ positions at (a) 150 K and (b) 300 K as obtained from MD simulations.

The most notable feature of the cage positions is the large variation in the position of the formate unit, corresponding to the partial rotation of the group around an axis going through the two oxygen atoms. While this motion becomes larger with increasing temperature, full rotations are not observed, even at 450 K. This is the case for both of the frameworks despite their differences in unit cell volume and cation shape and size.

In the FP-Gua structure the long axis of the formate hydrogen thermal ellipsoid (MD), which is a descriptor for the partial rotation of the formate unit, is longer for the formate unit binding the cage along the (long) b-direction (H1, 1.5 Å) than it is for those binding in the ac-plane (H2 1.1 Å). This is not surprising since the deviation from the Pnna space group is particularly large for the carbons and hydrogens of the formate groups along the b-axis; however, we note that the Zn–Zn distance along this axis is also shorter (5.76 Å) than in the ac-plane (6.03 Å). For the FP-Hy cage the thermal ellipsoids are very similar in size (1.2, 1.3, and 1.3 Å for H1, H2, and H3. respectively), and the Zn–Zn distances are almost equal (5.71, 5.73, and 5.73 Å).

The thermal ellipsoids of the molecules are shown in Figures and 5. If one looks first at the Gua+ molecules (Figure ), the positions are clearly fixed in space, and though the movement becomes larger at 450 K, the molecule is far from making a full rotation. This molecular behavior is very different from that of the methylammonium (MA) cation in the hybrid halide perovskite MAPbI3 which undergoes full rotations at room temperature[34] and the dimethylammonium cation in the formate perovskite series [(CH3)2NH2][M(HCO2)3] (M = Mg, Mn, Fe, Co, Ni, Cu, Zn) which has been found to rotate around an axis through the two carbon atoms at temperatures above 156 K for M = Zn.[35] Neutron diffraction data of the M = Fe compound showed large thermal displacements in the low-temperature phase, indicating that even below the transition temperature the molecule is only weakly fixed within the cavity.[36] The reason for the immobility of the guanidinium cation is the six strong N–H–O hydrogen bonds between the amine groups of the cation and the oxygens of the formate units (average H–O distance of 1.83 Å in the 0 K structure), resulting in a significant barrier for the rotation. The topology of the molecule furthermore means that all hydrogen bonds have to be broken simultaneously for the molecule to rotate. This is possibly also the reason for the relatively high mechanical strength of this structure as suggested in ref (37). Interestingly, the addition of Gua+ has been found to suppress the formation of defects in MAPbI3 solar cells.[38]

Hydrogen bonding is also preventing free rotation of the Hy+ ion; however, this molecule shows larger movements than the Gua+ molecule, as seen from the thermal ellipsoids at 150 and 300 K shown in Figure (at 375 K the cation is no longer well represented by a single position, and therefore the thermal ellipsoids at this temperature are not shown). The NH3+ group forms three short hydrogen bonds (H–O distances of 1.80, 1.77, and 1.73 Å in good agreement with previous calculations[39]) which are shorter than the two bonds formed by the NH2 group (1.97 and 2.08 Å). Consequently the NH2 group can move more than the NH3+ group, and at 300 K it is possible to break both hydrogen bonds of the NH2 group simultaneously, allowing for a full rotation around the N–N axis (an occurrence is seen approximately halfway through the MD trajectory given as Supplementary Information). This is probably the main reason for the discrepancy between the LD and the MD values for the thermal ellipsoids of the atoms in the NH2 group (N2, H11, and H12 in Table ), as such movements would not be described in the LD calculations. Rotation of the Hy+ cation has previously been proposed based on the line broadening in 1H magic angle spinning NMR spectra at ∼320 K,[39] though this technique cannot distinguish between rotation of the NH2 group and rotations of the whole molecule. Our simulations show that, on the time scale of the simulation, the NH2 group can rotate, while the NH3+ group remains fixed in a position displaced from the center of the cage. This displacement leads to a spontaneous polarization which we will discuss in the following.

Ferroelectric Properties

The spontaneous polarization, Ps, is calculated by performing Berry phase calculations of the polarizations along a path connecting the minimum energy structures with polarization Ps and −Ps, through a nonpolarized state. Figure a shows the result with the DFT optimized structure as end points and the coordinates of the intermediate configurations defined by different values of λ in eq . The spontaneous polarization is calculated to be 2.6 μC cm–2, which is somewhat less than the 3.48 μC cm–2 estimated in ref (8) based on the static position of the NH3+ group relative to the cage alone. This suggests that the polarization arising from the molecule is partially compensated by the cage as previously found for [NH4][Cd(HCO2)3][13] and [CH3CH2NH3][Mn(HCO2)3].[40]

Figure 6

(a) Polarization along a path connecting the equilibrium structures with polarization Ps and −Ps. The polarization at λ = 0.5 is set to 0. (b) Calculated thermal expansion (blue line) and experimental data points (magenta) taken from ref (8). (c) Equilibrium polarization as a function of temperature. Error bars indicate one standard deviation as estimated from the sum of displacements of the NH3+ groups from the center of their respective cages. (d) Distribution of displacements of the NH3+ group of one molecule from the center of the cage along the polarization axis at 150 K (top), 300 K (middle), and 375 K (bottom).

We investigate two possible effects of temperature on the polarization density: the changes arising from thermal expansion and the time-dependent variations in polariztion arising from the thermal movement of the atoms. The thermal expansion as calculated in the quasi-harmonic approximation is shown in Figure b. Note that the 0 K volume calculated here (761.91 Å3) is not the same as the 0 K volume found by DFT optimization of the unit cell size given in Table (742.20 Å3). This is because the DFT optimized volume is a minimization of the electronic energy only, while the quasi-harmonic optimized volume is found by minimizing the free energy which includes the zero point energy (ZPE). The quasi-harmonic volume is larger since a larger unit cell leads to softer phonons and thereby a lower ZPE, compensating a slightly higher electronic energy. Comparison of the quasi-harmonic volume with experimentally determined volumes at three different temperatures (magenta points in Figure b) shows very good agreement (less than 0.8% deviation). Furthermore, the slope of the curve, i.e., the volume expansion coefficient, follows the experimental trend closely.

The effect of volume on the equilibrium polarization, arising from thermal expansion of the lattice, was investigated. Using the calculated thermal expansion, the equilibrium polarization as a function of temperature can be extracted. The resulting curve, shown in Figure c, shows that the equilibrium polarization changes very little in the temperature interval between 0 and 352 K, falling by less than 0.1 μC cm–2.

The above calculation does not consider the variation in the instantaneous polarization arising from vibrations of the molecule and cage. Berry phase calculations of configurations from the MD simulation are made difficult by the movement of the atoms relative to the reference point. Instead we quantify the variation in polarization by the displacement of the NH3+ unit along the polarization axis (c for the Pna21 structure and b for the Pnma structure) relative to the midpoint between the Zn planes along this direction. (Similar results are obtained by defining the midpoint from the positions of the Zn atoms in the 0 K structure or from their instantaneous positions in the MD simulation.)

In the 0 K structure, which has polarization Ps, each NH3+ group is displaced 0.23 Å along the c-direction relative to the center of the cage. The average displacements at 150 and 300 K are almost identical (0.23 and 0.21 Å, respectively), indicating that the molecule oscillates in a nearly harmonic potential around the equilibrium position. Assuming a linear relationship between the displacement and the polarization, we can estimate the average polarization at those temperatures (blue dots, Figure c) again finding only a small change compared to the 0 K structure. The distribution of displacements during the 10 ps simulation is shown for one of the four molecules in the unit cell in Figure d for temperatures of 150 K (top), 300 K (middle), and 375 K (bottom). Each distribution is fitted with a Gaussian, from which a standard deviation can be extracted, to indicate how much the displacement deviates from the average value over time as a result of vibrations. A similar distribution can be made for the sum of the displacements for the four molecules, and the standard deviation can be used to estimate the standard deviation of the polarization at the given temperature. (The sum of variances should equal the variance of the sum in the case of completely uncorrelated molecules (no covariance). This is not exactly true in this case, indicating a small correlation between the four molecules.) We have plotted one standard deviation as error bars in Figure c for the two simulations at 150 and 300 K.

At 375 K the average displacement of each NH3+ group is calculated to be 0.05 Å. This is close to the average displacement in the optimized Pnma structure of 0.04 Å per molecule, arising from two molecules with a negative displacement and two molecules with a slightly larger positive displacement. Indeed, a calculation of the spontaneous polarization of the DFT optimized structures gives a value of 2.2 μC cm–2. The reason that this structure can be fitted into the centrosymmetric (zero polarization) space group Pnma can be understood by studying the distribution of displacements for one of the molecules as extracted from the MD simulation (Figure d). This molecule has an average displacement of 0.00 Å, arising from a change in the direction of the displacement during the simulation. Similar plots for the other three molecules show a distribution centered above (below) zero with a tail that extends into the region below (above) zero (cf. Supplementary Information for image). Thus, at this temperature the molecules can change the direction of the displacement on the time scale of the simulation; however, such events are not frequent enough that the total displacement averages to 0 over the simulation.

Electronic Structure

The dynamics of the crystal structure does not only affect the ferroelectric properties but also the electronic structure. To investigate the effect of electron–phonon coupling at room temperature, we calculate the bandgap of FP-Gua for 100 different configurations on the MD trajectory with a time separation of 0.1 ps, using the hybrid functional HSE06. The result at 300 K is shown in Figure . The bandgap varies in time with an average value of 5.96 eV, considerably lower than the value calculated for the 0 K equilibrium structure of 6.29 eV. The energy difference between the maximum and the minimum gap is found to be 0.65 eV with a root-mean-square (rms) deviation from the average of 0.14 eV. Looking only at the highest occupied crystalline orbital (HOCO) the difference between the maximum and the minimum value is 0.44 eV, while the rms deviation is 8.5 × 10–2 eV. As expected the variation and rms deviation of the bandgap are larger than those of the HOCO, since the instantaneous bandgap includes variations in both the HOCO and the lowest unoccupied crystalline orbital (LUCO) energy. The HOCO is located mostly on the Zn and O atoms while the LUCO is located on the formate groups with the largest weight on the carbon atoms (cf. Supplementary Information for images), and thus the fluctuations in the HOCO and LUCO are not expected to be strongly correlated.

Figure 7

Variation in the bandgap over time for the FP-Gua structure at 300 K. The histogram on the right shows the distribution around the average value.

While the bandgaps of the formate perovskites are too large to be of interest for most applications, variations are likewise expected to be found for the hybrid halide perovskites which are currently investigated for applications in solar cell materials. Indeed, the variations in the HOCO energy of MAPbI3 were investigated in ref (41), and a difference between the maximum and the minimum value of 0.47 eV at 319 K was found, with a rms deviation of 7.8 × 10–2 eV. These values are very similar to the ones we find for FP-Gua at 300 K, suggesting that the variation in the electronic structure is similar for the two structures, despite the different connection between the metal nodes. We note that the periodic boundary conditions may affect the variations in the bandgap, as the distorted structures occur periodically. The results do however stress the importance of correcting for the effect of temperature when matching experimental and theoretically calculated bandgaps for mechanically soft materials.

Conclusion

We have studied the dynamic behavior of the formate perovskites with the guanidinium and the hydrazonium cation. We find a similar behavior in the movements of the framework for the two structures, with small differences related to the size of the cations. We compare the thermal ellipsoids obtained from the MD simulations to those obtained from a harmonic phonon calculation and find that they are generally in good agreement, except for those of the hydrazonium cation, which is moving too much to be well described by the harmonic approximation.

We further considered the movement of the molecules within the cage. The Gua+ ion is static, with no rotations observed even at 450 K. The position of the Hy+ ion within the cage leads to ferroelectricity, and we calculate the spontaneous polarization density to be 2.6 μC cm–2. The equilibrium polarization shows only a weak temperature dependence; however, the variation in polarization due to molecular vibrations increases with temperature. Above the ferroelectric transition temperature the thermal energy allows the cation to change orientation on the time scale of our simulation. For sufficiently large simulation times this is expected to lead to the experimentally observed lack of polarization.

Finally, we look at the variation in bandgap during the simulation for the FP-Gua structure. We find that the average value of the bandgap is signifcantly lower than the value calculated for the 0 K structure, showing the importance of including temperature effects when comparing experimental and theoretical values.