Models of orientational disorder in hybrid organic-inorganic piezo-electric materials

Hybrid organic-inorganic materials show potential as flexible and Pb-free piezoelectric materials, but ferroelectric-to-paraelectric phase transitions are responsible for a loss of properties at increased temperatures. Here, we develop a model Hamiltonian from first-principles calculations for the archetypical hybrid organic-inorganic piezoelectric TMCMCdCl 3 (TMCM = trimethylchloromethyl ammonium, Me 3 NCH 2 Cl + ) to model the order-disorder phase transition. We show that the configurational entropy from molecular orientational disorder is alone unable the quantitatively reproduce the phase transition temperature, but inclusion of vibrational entropy leads to good quantitative agreement with experiments. This highlights the importance of vibrational contributions for the phase stability of soft and flexible materials. We furthermore show that defective structures and intergrowths are readily formed at ambient temperature which is suggested to be responsible for the exceptional piezoelectric response of these materials, which cannot be explained from the conventional small-displacement limit.


Introduction
Piezoelectric materials allow for the interconversion between mechanical and electrical energy and find applications in diverse areas including as sensors, actuator, and high precision motors.][3] However, these materials are brittle, require high processing temperatures, and may contain toxic and environmentally hazardous lead, which impose challenges for flexible and biocompatible technologies.
Hybrid organic-inorganic materials are commonly based on an inorganic framework structure hosting organic dipolar cations-often in perovskite or perovskite-inspired structures-and they commonly feature several temperature induced phase transitions with increasing degrees of orientational disorder of the organic cations at increasing temperature. 13Elevated temperature may lead to either stabilisation of desired phases such as phases with high dielectric constants for increased screening of charge-carriers in photovoltaics, 14,15 or to detrimental disordered paraelectric phases, in which the ferro-and piezoelectric performance is lost. 9,13hus it is of great importance to understand the influence of molecular orientational disorder and other entropic contributions in hybrid organic-inorganic materials on their stability, phase transitions and physical properties. 16,17n this article, we focus on the hybrid ferroelectric material, TMCMCdCl 3 (TMCM = trimethylchloromethyl ammonium, Me 3 NCH 2 Cl + ) with piezoelectric performance comparable to BaTiO 3 in its pure form, 9 and to PZT for its solid-solution with the fluorine analogue of the molecular cation (TMFM = trimethylfluoromethyl ammonium). 11TMCMCdCl 3 crystallises in the non-centrosymmetric monoclinic space group Cc at low temperature with ordered molecular cations, and transitions to a disordered paraelectric phase in the centrosymmetric hexagonal space group P6 3 /mmc above 400 K (Fig. 1). 9,11esides the interest in its phase transitions, the mechanism behind the high piezoelectric performance cannot be described within the conventional small-displacement limit, and is thus still poorly understood.Originally, it was suggested that ferroelastic switching was responsible, 9 but it has been argued that this is not possible when the piezoelectric response is measured with the Berlincourt method. 18Instead, it was suggested that phase switching between the ground state and a metastable phase is involved. 18,19hrough the last decade, several studies have proposed effective model Hamiltonians to study the order-disorder transitions in hybrid organic-inorganic materials, mainly those with cubic perovskite topology.These have been based on classical dipole-dipole interactions and strain models, [20][21][22] Pottsand Ising-like models, 23,24 or pair interactions. 25However, none of these have been developed solely from first principles and then applied to model a ferroelectric-to-paraelectric phase transition.
Here, we report the construction of a model Hamiltonian based on molecular pair interactions fitted to first principles calculated energies.Applying the model Hamiltonian, we perform Monte Carlo simulations to gain insight into the orderdisorder phase transition as well as nanoscale effects and defective phases that are likely to be involved in the piezoelectric performance.Particularly, we show the importance of including contributions from vibrational entropy to model phase transition temperatures in quantitative agreement with experiment.Finally, we show that the easy development of disorder in the ab-plane results in formation of isolated and extended defects, which are likely to be responsible for the exceptional piezoelectric response.

Model construction
A model Hamiltonian for the orientational order of TMCM cations was build using a model inspired by the pair mode concept introduced by Li et al. 25 .Using the hexagonal setting of the high temperature phase, the structure of TMCMCdCl 3 consists of one-dimensional face-sharing CdCl 6 -octahedra along the c-direction with TMCM molecules located around the 2c Wyckoff positions.Assuming that the molecular mirror plane is retained in the solid state and coincides with the mirror plane or the local 6m2 of the Wyckoff site, there are six different molecular orientations per site.
The low temperature phase is of monoclinic Cc symmetry, with its primitive cell coinciding with that of the high temperature phase.The conventional C-centered monoclinic unit cell can be constructed from the three primitive basis vectors a C = a P , b C = a P + 2b P and c C = c P , followed by a monoclinic distortion.Here, we use the primitive representation of the cell throughout.
In order to construct our model Hamiltonian, we include symmetry unique pair interactions corresponding to three different symmetry unique intermolecular vectors, (i) the interactions with the nearest neighbour at the other 2c site in the cell, i.e. with an intermolecular vector of [ 1 3 , − 1 3 , 1 2 ] (and symmetry equivalent directions), (ii) the nearest neighbour along the c-directions at an intermolecular vector of [0, 0, 1] (and symmetry equivalent directions), and (iii) the next-nearest neighbour in the ab-plane with an intermolecular vector of [1, 0, 0] (and symmetry equivalent directions).This gives rise to 13, 6 and 12 symmetry unique interactions along the first, second and third intermolecular separations.All three types of interactions are needed in order to distinguish the experimentally observed ground state from specific defective states that we will consider below.
The final energy expression for any supercell size can be written as where n i is the number of interactions of a given symmetry, and J i is their associated energy.It should be noted that in the final expression, the total number of parameters, J i , to be determined from fitting to energies calculated from DFT, is not equal to 31, but rather 23, since the parameters do not form a linearly independent basis-the simplest reduction being (considering the up-down orientation along the c-direction) that an equal number of head-to-head and tail-to-tail interactions of the molecules must be present due to the periodic boundary conditions.In practice, 8 parameters are eliminated using QR-decomposition of the original design matrix.This also means that the parameters cannot be interpreted directly as interaction strengths, as they will contain indirect information about other interaction types.Fig. 2 shows an illustration of how one of the training structures is represented as molecular vectors of the different sites.We represent the molecular cations as vectors with directions along the 1 13 -direction and its symmetry equivalents for the first 2c Wyckoff site, and the 113 -direction at its equivalents for the second one.This gives directions in agreement with experiment 9 and previous computations, 18 whereas we will use arbitrary units for the magnitude of the dipole moments and only consider their directions throughout.the 4 1x1x1 cells are chosen as the 4 symmetry unique ones possible within that cell size, and the rest are chosen randomly.
All supercells-including cell parameters-are relaxed with DFT to a convergence threshold of 0.01 eV Å−1 .The Vienna ab initio simulation (VASP) package with the projector augmented wave approximation is used for all simulations. 26,27The PBEsol functional is used, 28 the plane wave cutoff is set to 700 eV, and a 2×2×2 Γ-centered k-point grid is used for the base cell, and scaled accordingly for larger supercells.
The pair interaction parameters, J i , are fitted to the energies extracted from the DFT calculations using linear least-squares.

Phonon calculations
As we will see below, fitting our model Hamiltonian based only on internal energies from DFT leads to inaccurate predictions of the phase transition temperature in TMCMCdCl 3 .To get a better description of the energy as a function of temperature, we include the vibrational free energy calculated from harmonic phonon dispersions.
The phonon calculations are performed with PHONOPY 29 and forces are calculated with VASP.The structures are relaxed to a slightly stricter convergence criterion of 0.005 eV Å−1 , and 2×2×2 supercells (of the primitive base cell) are used for all phonon calculations.Vibrational free energies given as where ω qv are the phonon frequencies at wave vector q and band index v, are then calculated with PHONOPY.Due to the computational cost of calculating phonon dispersions for supercells, we opt to calculate them only for a subset of the structures.We then construct a surrogate model based only on the internal energy, E 0 , and equilibrium volume, V 0  relative to the ground state structure, to reproduce the relative free energy as where α and β are fitting parameters determined at each temperature of interest.The surrogate model reproduces the relative free energies quite well, especially at temperatures up to ∼600 K as seen in Fig. 4b.

Monte Carlo simulations
Monte Carlo (MC) simulations were performed with the Metropolis algorithm 30 using the model Hamiltonians described above, i.e. both the one based only on internal energy calculated with DFT, and the one based on free energies.In the second case, the temperature of the Monte Carlo simulation was matched with the free energy from phonon calculations.
A simulation box of 8×8×8 unit cells (with two molecules per unit cell) was used for all simulations, and 2 • 10 6 MC steps were performed with the first 1 • 10 5 steps discarded as equilibration.Each simulation was repeated 10 times to improve the statistical sampling.In all cases, the simulations were started from the ground state distribution of molecular orientations, as random (and high-temperature initialised) starting configurations always led to trapping in metastable states.

Results & Discussions
3.1 Model quality Fig. 3 shows the agreement between energies from the models based on the pair mode concept and the energies calculated from DFT.For the internal energies from DFT an RMSE of 4.76 meV/f.u. is obtained.Correcting the DFT energies for vibrational entropy results in a lower spread in energies and a smaller RMSE of only 3.79 meV/f.u. at 300 K.The models coarse grain over ionic relaxations etc., and considering the simplicity and short range nature of the model, the agreement 4 a Phonon dispersion and vibrational density of states for the ground state monoclinic structure of TMCMCdCl 3 .b Vibrational free energies of a set of metastable structures relative to that of the ground state.The solid lines indicate the calculated free energies from the phonon dispersions, and diamonds show the value from the correction model based only on the internal energy and the equilibrium volume of the metastable phases given in Eq. 3.
is good.Note however, that for both athermal and thermal energies, the ground state energy (lowest energy point in Fig. 3) is predicted to be larger using the model Hamiltonian than what is obtained with DFT.We may attribute this to lack of long-range terms in the model Hamiltonian, which tend to stabilise specific local environments, but less so long range interactions.This includes strain contributions, which are likely to be underestimated by the model.

Vibrational contributions
As we shall see below, simulations based on internal energies from DFT tend to overestimate the stability of the ground state structure.However, besides configurational entropy that we will model below with our MC model, vibrational entropy will also contribute to the energetics at elevated temperatures.Combining configurational and vibrational entropic contributions has also previously been shown to be important for accurate predictions of alloy stability. 31n Fig. 4a, the phonon dispersion and density-of-states of the ground state structure is shown.The vibrational degrees of freedom give rise to vibrational entropy, or more conveniently expressed as vibrational free energy as in Eq. 2. Fig. 4b shows the relative free energy between the ground state and a set of metastable structures that are used in the model construction.It is evident that all metastable structures are stabilised compared to the ground state at increasing temperatures, with one metastable state even becoming the thermodynamic ground state above ∼900 K.][34] Fitting the model for the relative free energy given in Eq. 3 at each temperature in 100 K intervals, give rise to predicted energies shown in Fig. 4b.As it is computationally expensive to calculate the vibrational free energies for all training structures, we instead use this surrogate model to estimate the relative free energy for all other structures, which in turn is used to fit a new model Hamiltonian at each temperature, with the model quality indicated in Fig. 3b for 300 K.

Polarisation & local correlations
From the MC simulations, we can extract expectation values of relevant structural parameters as a function of temperature.Of main interest for probing the ferroelectric-to-paraelectric phase transition in TMCMCdCl 3 is the total polarisation of the simulation cell.In Fig. 5, the polarisation as a function of temperature is shown for model Hamiltonians based on DFT internal energies and free energies including vibrational contributions.The polarisation is calculated from the molecular orientations using the (normalised) 1 13 and 113 family of directions as stated above for the orientation of the dipole of two different sites, and averaged over the number of molecules.
First, we note that the transition to a paraelectric phase, characterised by zero polarisation, is observed to occur in the 600-900 K range when neglecting effects of vibrational entropy, whereas it is predicted to occur in the range of 400-600 K based on free energies.The results including vibrational contributions are in good agreement with experiment, where the transition is observed close to 400 K. 9 Thus, the interplay between vibrational and configurational entropy is of great importance for the relative phase stability in hybrid organic-inorganic materials.The gradual transition observed here compared to the sharper transition in experiment is most likely a result of finite size effects of the simulation cell.
The ferroelectric-to-paraelectric phase transition is observed when considering the z-component of the polarisation, whereas the polarisation in the xy-plane gradually disappears below room temperature.This is in disagreement with experiment, where only one phase transition is observed in which all components of the polarisation are lost simultaneously. 9,11e suggest that this is a result of neglecting long-range interactions in the xy-plane.Thus, in our model, it is possible to form local regions where the favourable relative local orientations are retained, but at relatively little energetic cost, defectcenters and grain boundaries can be formed.Put in another way, the strain energy arising from lattice mismatch of different domains is only included indirectly as an energetic penalty of a local relative orientation, but not from its overall effect on the strain field.
Interestingly, the observed behaviour of the model around room temperature is reminiscent of the experimentally observed room temperature phases of the related TMCMCdBr 3 , 35 and the TMCM 1−x TMFM x CdCl 3 solidsolution at concentrations, x, between 0.3 and 0.5, 11  solid solution. 36In all three cases, a structure with P6 3 mc symmetry is observed, which corresponds to disorder between the three sites with the same z-polarisation.This suggests a strong competition between a fully ordered Cc-phase stabilised by long-range strain contributions, and an entropically stabilised-but still ferroelectric-P6 3 mc-phase.For TMCMCdBrCl 2 the Cc-phase is observed at lower temperature.Whether the two first mentioned systems also condense to the fully ordered Cc-phase at lower temperature is unknown.
The local correlations shown in Fig. 5 are defined as the average value of the dot-product between nearest neighbours in the xy-plane (the first interaction described in Section 2.1), and along the z-direction (the second interaction in Section 2.1).We note that local correlations are retained well above the paraelectric transition.Interestingly, the correlation between neighbours in the xy-plane remain strong after the overall polarisation in the plane is lost.This supports our claim that, locally, the structure remains close to the ground state structure, but the long-range order disappears due to neglecting strain contributions.
Furthermore, it is evident that strong correlations exist along the z-direction, i.e. even when the overall polarisation is zero, each molecule in a chain along the z-direction will prefer to have similar z-components of their polarisation.It would be of interest to probe the local order in the high temperature phase using diffuse X-ray (or neutron) scattering to further understand the merits of the model Hamiltonian beyond prediction of the phase transition temperature.

Intergrowths, defects and enhanced piezoelectric response
The large piezoelectric response of TMCMCdCl 3 cannot be explained by intrinsic effects as calculated from DFT. 9,18 Furthermore, it has been shown experimentally that the piezoelectric response is improved at the morphotrobic phase boundary between the P6 3 mc and Cc phases for its solid solutions. 11s mentioned above, there has been a discussion in the literature about the origin of the large piezoelectric response, which is beyond the standard displacive regime.Both effects of ferroelasticity 9 and metastable phases have been suggested. 18Based on the easy accessibility of xy-disorder observed in our MC model, we instead suggest that the origin may be related to defects and intergrowths of differently ordered domains in an ordered Cc-symmetry matrix.
In Fig. 6, we show two different models, one of an isolated defect corresponding to a different molecular orientation, and one of an intergrowth of a "plane" of differently ordered molecules.In both cases, the local orientations remain favourable, i.e. nearest neighbour interactions (in the xy-plane for the isolated defect, and in all directions for the intergrowth) remain the same as in the ground state structures, and should thus be expected to be of low formation energy.
For the isolated orientational defect, a defect formation energy of 41 (50) meV is found when unit cell relaxation is included (excluded).This corresponds to a thermal population at room temperature of 20 ( 14) % when considering Boltzmann statistics.
For the intergrowth model, the formation energy is 26 (39) meV with (without) unit cell relaxation.As this is an extended defect, it is not directly comparable to the thermal energy, but by normalising it by the area of the intergrowth plane, we can compare it to other calculated stacking fault energies.For our case we get 43 (62) meV/nm 2 , which, for comparison, is slightly lower than the formation energy of hexagonal stacking faults in cubic metal halide perovskites, 38 which are experimentally observed to be present. 39n the real solids, a complex mixture of defects and intergrowths will be present, and thus we shall not make an attempt to calculate the piezoelectric coefficient arising from these.However, it is very likely that applied stress can lead to reorientations of defects, movement of intergrowths etc., which will lead to a change in polarisation and therefore a piezoelectric response beyond the small-displacement limit.Furthermore, these defects are of significantly lower energy than the metastable phase suggested by Ghosh et al. 18 , and therefore more likely to be involved in the piezoelectric mechanism.

Conclusions
To conclude, we have build a model Hamiltonian for the orientational disorder in the hybrid organic-inorganic ferroelectric material TMCMCdCl 3 fitted to energies calculated from first principles both with and without inclusion of vibrational contributions.We show that the ferroelectric-to-paraelectric phase transition temperature can be quantitatively reproduced when vibrational entropic contributions are included, highlighting the importance of the interplay between vibrational and configurational degrees of freedom for soft and flexible crystalline materials.
While the temperature dependence of the z-axis polarisation is in great agreement with experiment, the xy-plane polarisation is not reproduced as well with the current model.Instead, it resembles phases observed in the related TMCMCdBr 3 , and the TMCM 1−x TMFM x CdCl 3 and TMCMCdBrCl 2 solid solutions.This shows that a close competition between strain stabilised long-range order and entropic stabilisation is present in this structural family.
Based on the observation of this competition, we calculate the formation energy of isolated and extended defects corresponding to different molecular orientations, and we show that these are readily formed.This suggests that reorientation of molecular defects and movement of intergrowth planes are involved in the exceptional piezoelectric response.
To experimentally validate the presence of large concentrations of both isolated and extended defects, accurate Xray diffraction experiments-possibly modelled with the maximum entropy method 40 -should be able reveal significant disorder on the molecular sites even in the low temperature phase, 41,42 and diffuse scattering experiment may reveal local correlations among these orientations. 37Furthermore, high resolution microscopy techniques such as transmission electron microscopy could lead to direct visualisation of (extended) defects, possibly through mapping of the associated strain fields. 43,44

Fig. 1
Fig. 1 Crystal structures of TMCMCdCl 3 .a and b show the high temperature disordered phase with P6 3 /mmc symmetry, and c and d show the low temperature monoclinic Cc phase, which is responsible for the ferroelectric and piezoelectric properties.Hydrogen atoms are omitted from b for clarity.

2. 2 Fig. 2
Fig. 2 Illustration of the structural representation in the pair mode concept.a Example of an input structure with randomly generated TMCM orientations, and b the corresponding representation of the structure with orientation vectors.Full (dotted) arrows note a vector pointing out of (into) the plane, and arrows in purple (turquoise) tiles are located at a z-coordinate of 1 4 ( 3 4 ).

Fig. 3
Fig. 3 Comparison of model energies based on Eq. 1 and DFT energies.a Model energies and internal energies from DFT, b model energies and temperature dependent energies corrected using Eq. 3 at 300 K.The orange lines are of unity slope and indicate perfect agreement between model and calculated energies.

Fig. 5
Fig. 5 Polarisation and local correlations as a function of temperature from MC simulations.a and b show the total polarisation, and c and d show the local correlations along the z-direction and in the xy-plane.a and c show the results for the model based on internal energies from DFT, and b and d show the results based on the energies corrected for vibrational free energies.

Fig. 6
Fig. 6 Defect and intergrowth structures in TMCMCdCl 3 .a An isolated orientational defect and c its schematic representation in terms of pair modes, b an extended defects corresponding to a two-dimensional intergrowth of differently oriented molecules and d its schematic representation.Note how the local nearest neighbour pair modes in the xy-plane remain the same in ground state and defect structures.