Predicting the structures and associated phase transition mechanisms in disordered crystals via a combination of experimental and theoretical methods

Abstract

Disordered materials make up a large proportion of condensed phase systems, but the difficulties in describing their structures and molecular dynamics limit their potential applications. Disordered crystalline systems, also known as plastic crystals, offer a unique perspective into these factors because the systems retain a degree of crystallinity, reducing the degrees of freedom that must be explored when interpreting the results. However, while disordered crystals do diffract X-rays, it is difficult to fully resolve meaningful crystalline structures, with the best scenario resulting in lattice parameters. In this study, we use a combination of experimental terahertz time-domain spectroscopy, and theoretical solid-state ab initio density functional theory and molecular dynamics simulations to fully elucidate the structures and associated dynamics of organic molecular solids. The results highlight that this combination provides a complete description of the energetic and mechanistic pathways involved in the formation of disordered crystals, and highlights the importance of low-frequency dynamics in their properties. Finally, with structures fully determined and validated by the experimental results, recent progress into anharmonic calculations, namely the quasi-harmonic approximation method, enables full temperature and pressure-dependent properties to be understood within the framework of the potential energy hyper-surface structure.

---

1 Introduction

Disordered (amorphous) solids are an important class of materials that have found use in a wide range of fields, from the semiconductor¹ to pharmaceutical industries.² The utility of amorphous materials lies in their metastable nature,^3,4 which ultimately leads to increased molecular dynamics compared to their crystalline counterparts.⁵ However, the actual structures and associated motions of the molecules within disordered solids remain a topic of debate, as traditional techniques, such as X-ray crystallography, provide little information that can be directly interpreted without significant approximations.⁶

Recently, the importance of low-frequency (≤330 cm⁻¹ or 10 THz) collective vibrations in bulk properties and phase transformation phenomena in molecular solids has been suggested,⁷ however the types of motions occurring in disordered materials have not yet been completely understood. Experimentally, low-frequency (terahertz) vibrational spectroscopy is a powerful tool for investigating solids, as the terahertz radiation is able to excite large-amplitude collective motions of entire molecules within the solid.⁸ This leads to both the internal conformational and intermolecular potential energy surfaces being probed simultaneously, providing a valuable measure of the bulk structure and dynamics. Terahertz time-domain spectroscopy (THz-TDS) is the far-IR extension of the more traditional FTIR technique, but unlike FTIR, which probes internal vibrations, there are no standard functional group specific terahertz vibrations, meaning each solid has a unique spectral fingerprint at terahertz frequencies.^9,10 Therefore, the interpretation of THz-TDS spectra requires significant computational modelling, and (at a minimum) the structures of the associated solids are necessary for initialising the simulations.^11,12

Disordered (or plastic) crystals are a class of materials that exhibit pseudo-order, in that the molecules occupy crystalline sites but lack true translational symmetry due to uncoupled motions within their crystallographic sites.^13,14 This leads to a scenario where the material diffracts X-rays (like a traditional crystal), but rather than yielding meaningful atomic-level structural details the data only provides a broad description of occupied sites within the solid, i.e. an average “sphere” of electron density at each lattice site.¹⁵ However, the psuedo-order limits the possible degrees of freedom that the molecules within the solid have, enabling more detailed investigation into the bulk structure and dynamics to be performed.^16,17

Because disordered crystals are ultimately disordered, diffraction methods often rely on other techniques to fully elucidate the static-ordered structures.¹⁸ While this has been successfully achieved in a number of cases, the static picture of a disordered crystal is not truly physical, since the molecules are free to move in their crystallographic sites uncoupled from their neighbours.¹⁹ In order to obtain a more physical description, molecular dynamics (MD) simulations can be used as a means of predicting the nature of the motions occurring in the solid,^20,21 as well as predicting the vibrational spectra through determination of either the molecular dipole moments (for infrared (IR) spectroscopies), or the molecular polarisabilities (for Raman spectroscopies).^22,23 However, traditional MD simulations are performed using user-defined force constants, which leads to biased vibrational results. A solution to this is ab initio molecular dynamics (AIMD), in which the forces are recalculated using density functional theory (DFT) at each individual time-step, leading to a more accurate description of the vibrational dynamics occurring within the system.^24–27

AIMD simulations also enable structural transformations to be investigated.^28–30 As many plastic crystals often undergo a transformation to an ordered crystal at low temperatures, it is possible to analyse the ordered form experimentally with diffraction methods to obtain an atomic-level structure. By using this structure as a starting point for the AIMD simulations, prediction of the transformation pathway as well as the disordered structure can be achieved.

In this study, a combination of experimental THz-TDS and theoretical static solid-state-DFT (s-DFT) and AIMD simulations are used to predict the structures, dynamics, and phase-transition pathways in the order–disorder transformation of parachlorobenzamide (PCB, Fig. 1).³¹ The static structures of both the ordered and plastic phases of the solid have been previously reported experimentally using synchrotron-based X-ray diffraction,^32,33 which enables the proposed methodology to be validated by experimental results. Overall, this combination of techniques allows the static structures, as well as the associated dynamical aspects, to be fully understood, providing insight into the energetic and mechanistic forces that drive the order–disorder phase transition in plastic crystals.

Fig. 1 Skeletal and three-dimensional models of a single PCB molecule.

2 Methods

2.1 Experimental methods

Parachlorobenzamide was purchased from Sigma-Aldrich (98%, Poole, UK) and was used as-received. Samples were prepared for THz-TDS measurements by mixing each sample with polyethylene to a 10% w/w concentration and subsequently grinding with a mortar and pestle to homogenise the sample and reduce particle size. A hydraulic press was then used to press the mixture into a 13 mm diameter freestanding pellet with a thickness of 2 mm. A corresponding pellet containing pure polyethylene was also made to act as a standard blank for absorption measurements. All experimental THz-TDS spectra were obtained using a commercial TeraPulse 4000 spectrometer (TeraView Ltd, Cambridge, UK). Variable-temperature measurements were performed using a liquid nitrogen cryostat (Janis, Massachusetts, USA) equipped with an externally controlled heating element (Lakeshore 330, Ohio, USA), permitting measurements from 80–400 K. The time-domain waveforms were symmetrically zero-padded and Fourier transformed to yield terahertz power spectra. The absorption spectra presented are a result of division of the sample spectra by that of a blank spectrum.

2.2 Computational methods

2.2.1 AIMD simulations. All AIMD simulations were performed using the freely-available CP2K software package,^34,35 which incorporates periodic boundary conditions. A 1 × 3 × 2 supercell (36 molecules, 576 atoms) was generated from the experimental structure of α-PCB in order to minimise interactions with neighbouring unit cells and fully capture the physical phase transformation process. Simulations were performed within the isothermal–isobaric (NPT) ensemble, with the temperature controlled using a Nosé–Hoover chain thermostat,^36,37 and pressure maintained using the extended Nosé–Hoover barostat of Martyna et al.^38,39 All atoms were modelled using the double-ζ DZVP-GTH basis set⁴⁰ and the Goedecker–Teter–Hutter (GTH) pseudopotentials.⁴¹ The Perdew–Burke–Ernzerhof (PBE) density functional⁴² and Grimme D3-dispersion correction^43,44 were used for all AIMD simulations. A triclinic cell was chosen so that all lattice components were free to relax throughout the simulations. A 1.0 fs time-step was used for all simulations, and each trajectory was run for a total of 20 000 fs (20 ps). The convergence criteria were set to ΔE ≤ 10⁻⁹ hartree for all simulations.

IR spectra were generated by performing NVT simulations starting from the equilibrated NPT supercells. A time-step of 0.5 fs was used for all vibrational analyses, for a total trajectory length of 30 ps. Both solids underwent an initial equilibration period prior to the production simulation, which was set to 3.0 ps and 5.0 ps for the α- and γ-forms, respectively. The molecular dipoles were determined via calculation of the Wannier centres, which was performed every 2.5 fs. The resulting trajectories were analysed using the TRAVIS post-processing software package,⁴⁵ and IR spectra were calculated by performing the Fourier transform of the dipole moment autocorrelation function.^24,27

2.2.2 S-DFT simulations. All s-DFT simulations were performed using the fully-periodic CRYSTAL17 software package.⁴⁶ In order to allow comparison to the AIMD simulations, the computational parameters were kept as similar as possible between the two methods. The atomic orbitals were represented using the atom-centred def2-SVP basis set,⁴⁷ and the D3-corrected PBE density functional was used for all calculations. The experimental α-PCB structure^32,33 and the AIMD-generated unit cell of the γ-form (vide infra) were used to initialise the respective simulations. In both cases the simulations were performed in the absence of any space group symmetry in order to mimic the AIMD methodology and remove any spurious constraints. Both structures were fully optimised to a convergence of ΔE ≤ 10⁻⁸ hartree, and vibrational analyses were performed using numerical differentiation within the harmonic approximation, with an energy convergence of ΔE ≤ 10⁻¹¹ hartree.^48,49 Infrared intensities were determined using the Berry phase approach.⁵⁰ Quasi-harmonic approximation (QHA) simulations were performed by simulating the optimised structures and vibrational properties at four volumes corresponding to a −2%, 0, +2%, and +4% contraction/expansion relative to the equilibrium structure.^51–54

3 Results and discussion

3.1 Terahertz time-domain spectroscopy

In order to determine the vibrational response of PCB over a wide temperature range, THz-TDS spectra were acquired from 100–350 K and the results are shown in Fig. 2. At low temperatures, there is at least one distinct feature at 1.23 THz, with a number of potential features present at higher frequencies, which is characteristic of a well-ordered crystalline solid. As the temperature is increased, the main absorption feature significantly broadens and shifts to lower frequencies, likely due to thermal expansion leading to a softening of the vibrational potential energy and an increased number of thermally-activated vibrational relaxation pathways, the latter of which leads to spectral broadening. It is important to note that while many materials experience these phenomena upon heating, the magnitude of the broadening observed here is rather large in comparison to most molecular crystals.^27,55 The origins of the spectral broadening might be related to the ordered solid deviating from an ideal crystalline geometry, as minor displacements from equilibrium positions would lead to slightly different vibrational potentials for each unit molecule. This would manifest in the experimental spectra as broadened features, as observed in the continuous phase transition of crystalline N-methyl-4-carboxypyridinium chloride.⁷ Additionally, a decrease in the intensity of the spectra can also be attributed to the increased vibrational motion with increasing temperature, as the IR intensity is proportional to the change in dipole moment with respect to the atomic displacement, .²⁷ Overall, it is likely that the experimental spectral broadening is due to a combination of all of these phenomena, which results in the large effect observed.

Fig. 2 Experimental THz-TDS spectra of PCB at various temperatures ranging from 100 K (blue) to 350 K (red). There is a marked transition between 310 K and 320 K, corresponding to the order–disorder transition, with all spectra corresponding to the α-form coloured blue, while all spectra corresponding to the γ-form are coloured red. The inset shows the absorption at 1.5 THz as a function of temperature, with the experimental transition temperature of 317 K noted by a dotted vertical line.

Upon increasing the temperature above 310 K, there is a clear change in the vibrational response of the material, with the spectra exhibiting a noticeable increase in intensity across the entire spectral region that continues to increase with rising temperature, a common trait of disordered solids and in contrast to crystalline solids, which typically undergo a decrease in absorption as temperature is increased. While there is apparently little difference in the shape of the individual spectra (i.e. both spectra appear to be featureless with an increasing baseline), there is an obvious difference in the temperature-dependence of the spectra. This is highlighted in the Fig. 2 inset, which shows the clear change in the absorption coefficient at 1.5 THz that is observed between 310 K and 320 K. This discontinuity in the spectral trends corresponds to the order–disorder phase transition in PCB, which is in good agreement with the published transition temperature of 317 K.^32,33,56 However, unlike fully-disordered materials, which do not exhibit any resolvable low-frequency vibrational features, the spectra of the disordered γ-PCB solid seemingly exhibit at least one absorption at ∼1.72 THz. It is important to note that experimentally the discontinuity in the THz-TDS measurements is reproducible upon heating and cooling, indicating that the kinetics of the transformation are faster than the experimental method can probe (minutes, due to the equilibration time required for temperature changes).

3.2 Structural elucidation and phase transformation pathway

In order to determine the structure, and associated dynamics, of the disordered γ-PCB solid, AIMD simulations were performed. Initially, a 1 × 3 × 2 supercell of the α-PCB solid was generated and used as the starting geometry for the AIMD simulations. Fully periodic NPT simulations were performed using CP2K, with the pressure set to 1 bar and the temperature varied from 100–400 K, with a time-step of 1.0 fs. In the high temperature simulations (above 300 K), the system quickly undergoes a large-scale transformation, with all of the PCB molecules experiencing large amplitude hindered-rotational type motions about their long axis (see Fig. 3). It is important to note that the AIMD simulations underestimate the transition temperature by ∼20 K, however, given the limitations of the exact implementation of thermal energy in the model, this represents a relatively minor error (20 K ≈ 0.04 kcal mol⁻¹).

Fig. 3 Bulk static structures of the two forms of PCB, the low-temperature ordered α-form (left) and the high-temperature plastic γ-crystal (right). The main structural difference between the two crystals is highlighted in orange, along with arrows denoting the hindered-rotational motion that drives the transformation between the two forms and corresponds to the vibrations discussed later in the text.

The phase transformation can be clearly observed through analysis of the angle between the planes formed by the two distinct molecular environments in the α-form. In low-temperature AIMD simulations (below 300 K) the molecules do not undergo any large-scale rearrangement, rather they simply oscillate within their crystallographic sites about their equilibrium positions, typical of an ordered molecular crystal (Fig. 4). However, at higher temperatures, there still exist oscillations, but with a much larger magnitude compared to the low-temperature case, as well as an overall trend towards a vanishing angle between the planes of the molecules (i.e. all molecules trending towards forming planar sheets). Additionally, the increased molecular motion in the high-temperature simulation results in an initial expansion of the simulation cell, which is a direct example of vibrational dynamics leading to thermal expansion. After ∼10 ps, the high-temperature simulation relaxes to a new bulk structure that generally exhibits an overall planar arrangement of all of the molecules within the system, however the individual molecules continue to undergo rather large-amplitude hindered-rotational type motions that are seemingly uncoupled from the bulk. Moreover, in the high-temperature simulation, the molecules are generally not as well-ordered as in the low-temperature simulation, with many deviations from the static crystal structure readily apparent. However, there do exist a number of local ordering events that make extracting a molecular basis possible, and this process was used to generate a crystallographic unit cell for the γ-form of PCB.

Fig. 4 AIMD simulated angles between adjacent molecular planes at 100 K (blue) and 300 K (red). A horizontal dashed line is drawn through the origin as a guide.

Six molecules (two stacked rows of three molecules each, with alternating dipole moments) were extracted from the equilibrated AIMD simulation and a unit cell was fit around them. The system was loaded into the fully-periodic CRYSTAL17 software package and optimised within the constraints of the input lattice vectors, in order to minimise any possibility of relaxation back into the α-form. Upon successful optimisation, the system was then allowed to fully relax with no constraints (lattice vectors and atomic positions). The resulting structure (Fig. 3) was found to be triclinic with P1̄ space group symmetry, lattice parameters of a = 4.934 Å, b = 5.383 Å, c = 13.794 Å, α = 93.931°, β = 116.421°, and γ = 93.498°, and a Z = 2 (Z′ = 1). This structure is in good general agreement with the previously suggested static-structure of γ-PCB (bond, angle, and dihedral angle RMSDs of 0.035 Å, 0.656°, and 2.552°, respectively), although with slightly contracted lattice vectors due to the s-DFT simulations effectively simulating a 0 K crystal.

In both simulations, the individual molecules undergo predominately rocking motions with respect to their neighbours, which is evident in the periodic fluctuations present in Fig. 3. Even in the plastic crystal, in which the molecules generally move freely within their sites, there is a clear correlated rocking-type motion present. Interestingly, these oscillations have periods of about 1–2 ps, which is on the order of the period of terahertz vibrations (1 THz = 10¹² Hz = 1 ps). With the structures of the two forms known, a more rigorous analysis of the vibrational dynamics occurring within the solids can be performed.

3.3 Vibrational dynamics

The vibrational dynamics of both PCB solids were predicted using the AIMD simulations, and the results are shown in Fig. 5. In both cases, the AIMD simulations generally reproduce the experimental spectra, yielding a number of features that correspond to those experimentally observed. However, the AIMD spectra provide little specific atomic-level detail surrounding the nature of each individual vibrational mode, leaving further analyses necessary. In this regard, s-DFT simulations are able to determine the precise atomic-motions that produce each vibrational transition, and this method was employed once the static structure of the γ-form was determined. The static crystal structures of both forms of PCB were fully optimised with no constraints using CRYSTAL17, and subsequently the phonon properties were calculated using a finite difference method within the harmonic approximation. The results of the vibrational simulations are shown in Fig. 5 for both crystals.

Fig. 5 Experimental THz-TDS spectra (blue) and s-DFT simulated IR spectra (black) for both PCB solids. The α-PCB spectrum (80 K) is predicted to have a number of IR-active modes which primarily correspond to hindered-rotational motions of the individual molecules within the solid. The γ-form contains fewer IR-active modes (due to the smaller static unit cell compared to the α-form), however there are additional IR-inactive modes (dotted lines) predicted that might be activated in the experimental disordered structure. In both panels, the simulated AIMD spectra are also shown as dotted red lines, and are in good general agreement with the experimental spectra. Note that the experimental spectra have been offset by +1.5 cm⁻¹ for clarity.

In the α-form, there are a number of low-frequency vibrations predicted that result in good agreement with the experimental spectrum. Overall, the predicted spectrum is slightly over-estimated in terms of frequency, which is likely due to the simulations being performed at 0 K while the experimental spectrum was recorded at 80 K, as in some cases it has been shown that spectra continue to shift and sharpen upon cooling to temperatures as low at 4 K.⁵⁷ Indeed, the s-DFT simulations are also blue-shifted compared to the 100 K AIMD simulations, further indicating that the spectrum is strongly influenced by temperature. Additionally, AIMD vibrational analysis performed at higher temperatures confirms that the spectra become much more diffuse and red-shifted as the transition temperature is approached (Fig. 6). The temperature-dependent AIMD-simulated vibrational spectra generally reproduce the experimental results, predicting a red-shift of most features, as well as an increase in absorption intensity between ∼1.25 THz and 1.85 THz, which originates from the combination of spectral shifting and broadening.

Fig. 6 AIMD-simulated IR spectra of α-PCB at 100, 200, and 250 K (blue, purple, and red, respectively). The spectra exhibit a significant red-shift and spectral broadening as the temperature is increased.

The vibrational mode-types are mainly hindered rotations of the individual molecules with respect to one another, along both the long molecular axis (i.e. the Cl-ring bond, see Fig. 3) as well as librations orthogonal to the plane of the molecule, with the exception of the weak predicted feature at 1.645 THz which is a translational-type motion. The well-resolved experimental transition at 1.23 THz, and predicted at 1.337 THz, involves a coupled hindered rotation of all of the molecules within the solid along their long-axes (see Fig. 3), which at its maximum displacement yields a totally planar arrangement of all of the molecules within the solid. Thus, this mode corresponds well to the phase-transformation pathway predicted from the AIMD simulations. Interestingly, this mode is strongly influenced by temperature both experimentally (see Fig. 2) and theoretically (Fig. 6), and this assignment helps to explain the significant broadening and shifting that occurs with increasing temperature, as the vibrational potential likely becomes much softer and more anharmonic (leading to more relaxation pathways) with increasing thermal expansion (and thermal energy).

In the γ-form, there are fewer IR-active modes predicted due to fewer molecules within the crystallographic unit cell. The resolvable features at 1.72 THz can be assigned to an intense mode predicted at 1.755 THz. This is a surprising result given most disordered materials do not show discernible low-frequency absorption features, and it is even more surprising that the static-model is able to reproduce the disordered experimental results. Investigation of the vibrational mode type highlights that the vibrational motion is an asymmetric rocking of two molecules about their long axes with respect to one another. This directly corresponds to the coupled oscillations present in the AIMD simulations between pairs of neighbouring molecules (see Fig. 4), and also corresponds well to the oscillation period, indicating that this particular feature is indeed the local hindered-rotational motion suggested. Furthermore, the observation of this mode experimentally confirms that this vibration is not a long-range collective vibration, which would lead to a very broad experimental absorption, but rather is predominantly a local effect, being primarily dependent on a single molecular neighbour.

However, the s-DFT predicted spectrum fails to reproduce the correct IR density of states in the terahertz region as shown by the experimental and AIMD results. This is not surprising, as the static model is not a true representation of the disordered solid. However, there are a number of IR-inactive modes predicted which could become activated upon the removal of crystalline symmetry, and interestingly these transitions are predicted at similar frequencies to those predicted by the AIMD simulations. These modes, shown in Fig. 5 as dotted sticks, also correspond to hindered-rotational motions, with the lowest-frequency mode (1.296 THz) corresponding well with the mode type of the 1.337 THz vibration in the α-form. Indeed, the occurrence of a transition in both the experimental and AIMD-predicted spectra at ∼1.17 THz confirms that in the absence of true crystalline symmetry this vibration can become activated.

In order to further predict the role that thermal expansion plays in the structures and vibrational dynamics of the PCB crystals, QHA simulations were performed. The QHA method performs volume-constrained geometry optimisations (all lattice vectors and atomic positions were allowed to relax while maintaining a constant volume) and frequency analyses at four distinct volumes corresponding to the equilibrium (no volume constraint) structure, and −2, +2, and +4% contracted/expanded volumes. Based on these results, the temperature dependence of the vibrational features can be estimated within the context of thermal expansion effects. The results indicate that all of the low-frequency vibrations in both crystals are dependent on the volume to a large extent, with most modes experiencing shifts on the order of 10 cm⁻¹ between the two extreme volumes. This is also reflected in the temperature-dependent AIMD simulations (see Fig. 6), which also predict significant red-shifting upon increasing the temperature. Remarkably, the vibrational modes that involve motions corresponding to the phase-transformation mode (i.e. hindered rotations about the molecular plane) involve even larger shifts, for example the mode at 1.337 THz in the α-form is predicted to have a 26 cm⁻¹ red-shift. These data highlight the importance of unit cell dimensions in the low-frequency vibrational potentials in molecular solids.

3.4 Thermodynamic analysis

The accurate modelling of the vibrational spectrum by the s-DFT simulations lends confidence in the theoretical methodology employed, and enables the prediction of thermodynamic parameters to be performed. It is important to note that while the s-DFT predicted IR (and Raman) activities might not be completely accurate for certain modes in the disordered-crystal, thermodynamic functions depend more simply on the correct identification of the vibrational states, regardless of their symmetry. The harmonic canonical vibrational partition function was determined, and related quantities such as the entropy and free energy were calculated. All the energies discussed have been corrected for basis set superposition error using the counterpoise correction method.⁵⁸ When strictly considering the electronic energy, which has no intrinsic temperature or pressure dependence, the α-crystal is the more stable system by 0.56 kJ mol⁻¹. However, when the vibrational results are used to determine the temperature-dependent entropic and free energy parameters, the results suggest that the γ-solid becomes the more thermodynamically stable crystal above 301 K, close to the experimental transition temperature of 317 K (Fig. 7).

Fig. 7 Harmonic (top) and quasi-harmonic (bottom) free energy curves for the α- and γ-PCB crystals (red and blue, respectively). The vertical dotted line is drawn at the intersection of the two curves, occurring at 301 K and 326 K for the harmonic and quasi-harmonic simulations, respectively. In both charts the data have been offset relative to the value of the α-form free energy at 50 K.

In an attempt to increase the accuracy of the calculated thermodynamic constants, the QHA vibrational results were used to determine the temperature-dependent free energies. Because these simulations explicitly take thermal expansion into account, they are able to provide a more physical description of the thermodynamic forces involved in these materials. The simulated QHA free energy curves, also shown in Fig. 7, not only show a larger deviation in energy between the two solids, but also provide a more accurate description of the phase-transition temperature, predicted at 326 K, only 9 K higher than the experimentally observed temperature. The greater accuracy achieved when explicitly considering the thermal expansion properties of the solids, while not surprising, highlights the importance of taking into account these forces when attempting to describe the thermodynamics of weakly bound molecular crystals.

3.5 Interpretation and further developments of proposed method

The combination of experimental and theoretical results enables a great deal of information surrounding the phase transition behaviour in molecular solids to be understood. Thus, such a methodology is generally applicable to a large variety of materials, particularly those that undergo diffusionless transformations. The sensitivity of THz-TDS measurements to the bulk packing arrangement, and any minor deviations from ideal crystalline ordering, leads to it being a powerful probe of the processes and dynamics occurring within solid materials. While X-ray diffraction is often used for interpreting this type of behaviour, it is in reality a complement to the THz-TDS experiment, as the former provides information related to the average structure, while the latter probes molecular dynamics. The addition of computational results provides a means of interpreting the experimental spectra, indicating what types of motions are occurring within the solids. Additionally, the accurate reproduction of the experimental THz-TDS spectra indicates that the simulations are correctly reproducing the forces, including the weak intermolecular interactions, providing a high level of confidence in them. However, not all systems are well suited for study using this methodology, particularly in predicting large-scale structural reorganisations, which would require much longer AIMD simulations in order to capture the effect. Nonetheless, the proposed methodology is well suited for order–disorder transitions, particularly in molecular solids where such transformations occur near ambient conditions.

4 Conclusions

Overall, the proposed combination of experimental terahertz time-domain spectroscopy and theoretical ab initio molecular dynamics and static-solid-state-density functional theory methods provides a complete means to elucidate the structures, phase transition mechanisms, and bulk properties (e.g. thermodynamics and mechanochemical effects such as thermal expansion, and other related quantities) of molecular solids. The results shed light on the importance of low-frequency vibrational modes in many aspects of the properties of disordered crystals, and the atomic-level details provided by the s-DFT simulations enable experimental observations regarding temperature-dependent spectral trends to be fully rationalised. And while traditional harmonic approximation vibrational analyses were shown to be a suitable means of determining the thermodynamic properties of molecular crystals, the results can be significantly improved by taking into account thermal expansion effects via the use of quasi-harmonic approximation simulations, which provide a more physical (and accurate) description of both the temperature-dependent vibrational and thermodynamical forces present within such materials.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the UK Engineering and Physical Sciences Research Council (EPSRC) for funding (EP/N022769/1). Simulations were possible via our membership of the UK’s HEC Materials Chemistry Consortium, which is funded by the EPSRC (EP/L000202), and this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). JAZ and JK thank the EPSRC Cambridge Centre for Doctoral Training in Sensor Technologies and Applications (EP/L015889/1). JK thanks the German Academic Scholarship Foundation for funding.

Notes and references

1. V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Chem. Rev., 2007, 107, 926–952.
2. S. B. Murdande, M. J. Pikal, R. M. Shanker and R. H. Bogner, J. Pharm. Sci., 2010, 99, 1254–1264.
3. A. Cavagna, Phys. Rep., 2009, 476, 51–124.
4. M. Oguni and N. Okamotot, Solid State Commun., 1996, 99, 53–56.
5. M. T. Ruggiero, M. Krynski, E. O. Kissi, J. Sibik, D. Markl, N. Y. Tan, D. Arslanov, W. van der Zande, B. Redlich, T. M. Korter, H. Grohganz, K. Löbmann, T. Rades, S. R. Elliott and J. A. Zeitler, Phys. Chem. Chem. Phys., 2017, 19, 30039–30047.
6. E. Shalaev and A. K. Soper, J. Phys. Chem. B, 2016, 120, 7289–7296.
7. M. T. Ruggiero, J. Axel Zeitler and T. M. Korter, Phys. Chem. Chem. Phys., 2017, 19, 28502–28506.
8. E. P. J. Parrott and J. A. Zeitler, Appl. Spectrosc., 2015, 69, 1–25.
9. J. A. Zeitler, P. F. Taday, D. A. Newnham, M. Pepper, K. C. Gordon and T. Rades, J. Pharm. Pharmacol., 2007, 59, 209–223.
10. G. J. O. Beran, Chem. Rev., 2016, 116, 5567–5613.
11. M. D. King, T. N. Blanton, S. T. Misture and T. M. Korter, Cryst. Growth Des., 2011, 11, 5733–5740.
12. M. T. Ruggiero, J. Sibik, R. Orlando, J. A. Zeitler and T. M. Korter, Angew. Chem., Int. Ed., 2016, 55, 6877–6881.
13. J. Timmermans, J. Phys. Chem. Solids, 1961, 18, 1–8.
14. J. M. Pringle, P. C. Howlett, D. R. MacFarlane and M. Forsyth, J. Mater. Chem., 2010, 20, 2056–2062.
15. M. Descamps, N. T. Correia, P. Derollez, F. Danede and F. Capet, J. Phys. Chem. B, 2005, 109, 16092–16098.
16. F. Affouard, J. F. Willart and M. Descamps, J. Non-Cryst. Solids, 2002, 9, 307–310.
17. A. Hédoux, Y. Guinet, P. Derollez, J.-F. Willart, F. Capet and M. Descamps, J. Phys.: Condens. Matter, 2003, 15, 8647–8661.
18. L. Jin, K. M. Nairn, C. M. Forsyth, A. J. Seeber, D. R. MacFarlane, P. C. Howlett, M. Forsyth and J. M. Pringle, J. Am. Chem. Soc., 2012, 134, 9688–9697.
19. R. Brand, P. Lunkenheimer and A. Loidl, J. Chem. Phys., 2002, 116, 10386–10401.
20. F. Affouard and M. Descamps, Phys. Rev. Lett., 2001, 87, 035501.
21. F. Affouard, E. Cochin, R. Decressain and M. Descamps, Europhys. Lett., 2001, 53, 611–617.
22. O. B. M. H. Duparc and M. Meyer, J. Chem. Phys., 1998, 93, 1313–1324.
23. J. Adebahr, F. C. Grozema, S. W. deLeeuw, D. R. MacFarlane and M. Forsyth, Solid State Ionics, 2006, 177, 2845–2850.
24. M. Thomas, M. Brehm, R. Fligg, P. Vöhringer and B. Kirchner, Phys. Chem. Chem. Phys., 2013, 15, 6608–6622.
25. M. Yamaguchi and A. Ohira, Comput. Theor. Chem., 2016, 1089, 54–58.
26. M. Thomas, M. Brehm and B. Kirchner, Phys. Chem. Chem. Phys., 2015, 17, 3207–3213.
27. M. T. Ruggiero and J. A. Zeitler, J. Phys. Chem. B, 2016, 120, 11733–11739.
28. G. Kresse and J. Hafner, Phys. Rev. B, 1994, 49, 14251–14269.
29. A. S. Larsen, M. T. Ruggiero, K. E. Johansson, J. A. Zeitler and J. Rantanen, Cryst. Growth Des., 2017, 17, 5017–5022.
30. S. Tsuneyuki, H. Aoki, M. Tsukada and Y. Matsui, Phys. Rev. Lett., 1990, 64, 776–779.
31. A. Kubo, R. Ikeda and D. Nakamura, Z. Naturforsch., A: Phys. Sci., 1986, 41, 270–274.
32. A. Schönleber, P. Pattison and G. Chapuis, Z. Kristallogr. – Cryst. Mater., 2003, 218, 47.
33. Y. Takaki, K. Nakata, T. Taniguchi and K. Sakurai, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1978, 34, 2579–2586.
34. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 4, 15–25.
35. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
36. S. Nosé, J. Chem. Phys., 1998, 81, 511–519.
37. S. Nosé, Mol. Phys., 2002, 100, 191–198.
38. G. J. Martyna, D. J. Tobias and M. L. Klein, J. Chem. Phys., 1994, 101, 4177–4189.
39. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Mol. Phys., 1996, 87, 1117–1157.
40. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
41. S. Goedecker and M. Teter, Phys. Rev. B, 1996, 54, 1703–1710.
42. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
43. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
44. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
45. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
46. R. Dovesi, R. Orlando, A. Erba, C. M. Zicovich-Wilson, B. Civalleri, S. Casassa, L. Maschio, M. Ferrabone, M. De La Pierre, P. D’Arco, Y. Noël, M. Causà, M. Rérat and B. Kirtman, Int. J. Quantum Chem., 2014, 114, 1287.
47. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
48. F. Pascale, C. M. Zicovich-Wilson, F. L. Gejo, B. Civalleri, R. Orlando and R. Dovesi, J. Comput. Chem., 2004, 25, 888–897.
49. C. M. Zicovich-Wilson, F. Pascale, C. Roetti, V. R. Saunders, R. Orlando and R. Dovesi, J. Comput. Chem., 2004, 25, 1873–1881.
50. Y. Noël, C. Zicovich-Wilson, C. M. Zicovich-Wilson, R. Dovesi, B. Civalleri, P. D’Arco and R. Dovesi, Phys. Rev. B, 2001, 65, 014111–014119.
51. A. Erba, J. Maul and B. Civalleri, Chem. Commun., 2016, 52, 1820–1823.
52. A. Erba, M. Shahrokhi, R. Moradian and R. Dovesi, J. Chem. Phys., 2015, 142, 044114.
53. M. T. Ruggiero, J. A. Zeitler and A. Erba, Chem. Commun., 2017, 53, 3781–3784.
54. J. G. Brandenburg, J. Potticary, H. A. Sparkes, S. L. Price and S. R. Hall, J. Phys. Chem. Lett., 2017, 8, 4319–4324.
55. M. T. Ruggiero, T. Bardon, M. Strlič, P. F. Taday and T. M. Korter, Phys. Chem. Chem. Phys., 2015, 17, 9326–9334.
56. P. G. Karamertzanis and C. C. Pantelides, Mol. Phys., 2007, 105, 273–291.
57. Y. C. Shen, P. C. Upadhya, E. H. Linfield and A. G. Davies, Appl. Phys. Lett., 2003, 82, 2350–2352.
58. M. Gutowski and G. Chalasiński, J. Chem. Phys., 1998, 98, 5540–5554.

---

Footnote

† Present address: Department of Chemistry, University of Vermont, 82 University Place, Burlington, VT 05403, USA.

---