Pressure-driven phase transition mechanisms revealed by quantum chemistry : L-serine polymorphs †

The present study delivers a computational approach for the understanding of the mechanism of phase transitions between polymorphs of small organic molecules. By using state of the art periodic DFT calculations augmented with dispersion corrections and an external stress tensor together with gasphase cluster calculations, we thoroughly explained the reversible phase transitions of three polymorphs of the model system, namely crystalline L-serine in the pressure range up to 8 GPa. This study has shown that at the macroscopic level the main driving force of the phase transitions is the decrease in the volume of the crystal unit cell, which contributes to the enthalpy difference between the two forms, but not to the difference in their internal crystal energies. At the microscopic level we suggest that hydrogen bond overstrain leads to a martensitic-like, cooperative, displacive phase transition with substantial experimental hysteresis, while no such overstrain was found for the ‘‘normal type’’, atom per atom, reconstructive phase transition. The predicted pressures for the phase transitions deducted by the minimum enthalpy criterion are in reasonable agreement with the observed ones. By delivering unambiguous explanations not provided by previous studies and probably not accessible to experiment, this work demonstrates the predictive and explanatory power of quantum chemistry, confirming its indispensable role in structural studies.

The control of polymorphism is becoming increasingly important to understand the basics of crystal structure formation, but also to ensure control over the properties of various molecular materials, including non-linear optics, explosives, hybrid materials, etc. This is a crucial point especially in the pharmaceutical industry. [1][2][3][4] The progress in lattice energy calculations achieved in recent years makes it possible to explain, and in some cases -even predict the existence of polymorphs and their crystal structures. [5][6][7][8][9][10][11] However, it remains a challenge to find the experimental conditions to obtain a predicted polymorph in the laboratory. 12 The main routes for obtaining a new polymorph include crystallization (under very different conditions) as well as solid-state phase transitions or solid-state chemical reactions. The mechanisms of these processes and the role of thermodynamic and kinetic factors in the formation of a given polymorph often remain unclear. In particular, though many pressure-induced phase transitions in molecular crystals have been documented, it remains difficult to rationalize a priori if and when a solid-state transition may occur for any selected compound. [13][14][15][16][17][18][19] In the present contribution we demonstrate the capability of density functional theory (DFT) calculations to elucidate the mechanisms of pressure-induced phase transitions and greatly improve our understanding of this compelling topic. In particular, the present work focuses on polymorphs of molecular compounds with conformationally flexible molecules that are linked in a hydrogen bond (H-bond) network.
L-Serine polymorphs have been chosen as a case study. The behavior of this crystalline amino acid with increasing pressure has been extensively studied experimentally by optical microscopy, Raman spectroscopy, powder and single-crystal X-ray diffraction, and neutron powder diffraction. [20][21][22][23][24][25][26] Two reversible isosymmetric (space group P2 1 2 1 2 1 ) phase transitions have been reported, giving polymorph II at 5.3 GPa, and polymorph III at 7.8 GPa. ‡ The transition to phase II is of a martensitic type with a pronounced hysteresis. 20,24 In a single crystal, as soon as the phase transition starts, it propagates rapidly through the crystal until the transformation is complete, and the single crystal is preserved. At the same time, in a polycrystalline sample two phases -polymorphs I and II -can co-exist in a wide pressure range. 24,26 Apparently, nucleation is the rate-limiting step of this transformation. The phase transition to phase III also preserves a single crystal intact; no visible changes occur, whereas Raman spectra and X-ray diffraction patterns change. Structural data (including atomic coordinates) at multiple pressure points were obtained from X-ray diffraction and later from neutron powder diffraction enabling us to follow also the changes in molecular conformations and H-bonding. 22,24,25 With increasing pressure, hydrogen bonds were supposed to act as ''springs'' and -CH 2 OH side chains -as ''joints'': the two pressure-induced phase transitions are accompanied by distortion, switching over of existing hydrogen bonds (from the OHÁ Á ÁOH type in L-serine I to more preferable OHÁ Á ÁOOC bonds in L-serine II), as well as by formation of new hydrogen bonds (additional NHÁ Á ÁO hydrogen bonds in L-serine III). These changes in the hydrogen bond network were supposed to be the driving force of the phase transitions. However, there was no proof that optimizing the H-bonds (all or some of them), rather than molecular packing and non-directional van der Waals interactions, is the main driving force of the phase transitions.
Periodic DFT calculations enhanced by dispersion corrections and an external stress tensor were performed to model L-serine solid state behavior at different pressures. The models for all calculations were built on the basis of the experimental X-ray crystal and neutron powder diffraction data collected at different pressures. All solid state calculations were performed by the program package VASP 27-30 using the functional of Perdew, Bruke and Ernzerhof (PBE), 31 a plane-wave basis set with a kinetic energy cutoff of 500 eV and the Projector Augmented Wave atomic pseudopotentials. 32,33 The integrals in the reciprocal space were calculated on a Monkhorst-Pack mesh of 8 Â 8 Â 4 k-points. 34 The effects of external pressure were enforced by the stress tensor (PSTRESS keyword) corresponding to a selected pressure value in the range between 0.0001 GPa and 8.1 GPa.
Gas phase cluster models were used to study intermolecular interactions between pairs of molecules. All molecular pairs were extracted from crystal structures optimized at corresponding pressure and the distance between these molecules was changed while the geometry remained fixed. The Gaussian 09 program 35 was used to study the intermolecular interactions between pairs of molecules. The calculations were performed using the M06-2X/6-311++G(2d,2p) 36 level of theory. Technical details of the performed calculations, accuracy validation and additional data are described in the ESI. †

Crystal energies and enthalpies of the three polymorphs
The values of crystal energy and enthalpy calculated for L-serine crystal structures (optimized experimental structures documented to dominate at these pressures) are plotted in Fig. 1.
The I -II phase transition near 5 GPa manifests itself in the changes in the calculated internal crystal energy, which increases jump-wise at the transition point. The II -III phase transition is not evident from the same plot ( Fig. 1). § This type of calculation is usually performed in studies devoted to polymorphism or high-pressure phase transitions. However, this technique cannot help to explain why the phase transitions in fact occur, since we cannot compare the enthalpies and crystal energies of different phases at the same pressure. Therefore, as the next step we have calculated the energy characteristics for all polymorphs in the same pressure range, beyond the range in which they have been observed experimentally. Even if a polymorph could never be obtained at a selected pressure in reality, its structure and energy can be still simulated by DFT calculations by optimizing it to the selected pressure conditions. To predict a hypothetical structure out of the region of structural stability, it was possible to take the structure at experimental pressure and provide optimization with programmed pressure (e.g., taking initial experimental atom coordinates and cell parameters for ambient pressure and running optimization at 5.0 GPa). This technique definitely does not consider phase transitions. Instead, it allows for prediction of which hypothetical structure the same polymorph would adapt at any pressure if no phase transition took place. Such an approach is beyond experimental assessment but can reveal very important information about phase transitions. Unit cell parameters of the optimized structures can be found in Table S1 (ESI †).
Laborious verification of correctness of the employed computational approach was done (see Accuracy check, ESI †). After the calculations proved to reproduce correctly the structures in the range where they have been observed experimentally, the same methodology was used to predict hypothetical structures of all three polymorphs out of the range of their existence.
A comparison of the changes in lattice energies and enthalpies of the three L-serine polymorphs in a wide pressure range 0-9 GPa (neglecting the entropic term) shows that the decrease in volume gives a major contribution to the total decrease in free energy according to these calculations. For I -II phase transition the volume decrease can be considered as the main driving force. As for II -III transition, the same conclusion cannot be made unambiguously because of the lower accuracy of the computational technique with respect to the prediction of volume changes over pressure for polymorphs II and III. These results could be adjusted with better accuracy if the structures of polymorph II and especially polymorph III could be refined at more pressure points experimentally. Up to now the contributions of energy and PV terms to the enthalpy change during II to III transition seem to be at least comparable and in agreement with previous experimental and theoretical work. 22,24,25 § Phase transitions are not evident from the enthalpy plot because the PV term has more impact (hundreds of kJ mol À1 ) on the enthalpy in comparison to the crystal energy (tens of kJ mol À1 ), so the PV term masks crystal energy changes; the enthalpy change versus pressure looks therefore nearly linear (R 2 = 0.9863) in the whole range. This result is similar to the one reported previously. 19 The phase transitions from I to II and from II to III can be expected when the free energies (in the approximation of the calculations used here -the enthalpies) of the corresponding phases become equal. These calculations predict this to happen at 3.7 GPa and 5.3 GPa, respectively (Fig. 2).
For both phase transitions the calculated transition pressures are smaller than the experimental values. This can be a consequence of the limitations of the model, but can also reflect the kinetic control over the phase transitions, so that the transition is observed at pressures higher than the equilibrium value. It is known from experiments that L-serine I, II, III can co-exist in a wide pressure range, and the transitions could be triggered at different pressure values, depending on the experimental conditions. 24,26 Energy-volume profiles were calculated for all three polymorphs of L-serine, to estimate the zero-pressure bulk modulus (Fig. S3, ESI †). L-Serine I showed the highest bulk modulus value of 23.4 GPa, whilst those for L-serine II and III appear to be significantly less and very similar to each other À 14.7 GPa and 13.9 GPa respectively (Table S2, ESI †). This reflects similarity in the H-bond network in polymorphs II and III and a significant difference in comparison to L-serine I.

Modeling effect of pressure on individual H-bonds
Changes in the phase transitions of L-serine are accompanied by significant changes in the H-bond networks: some H-bonds are broken, whereas new ones are formed (Fig. 3).
Such phenomena have been documented for many pressureinduced phase transitions in molecular crystals. However, it always remains an open question whether the observed change in an H-bond is directly caused by increasing pressure to optimize its energy, or if it is an indirect consequence of some other structural changes, for example of increasing the packing density. 37 In the present work we have modeled the effect of pressure on individual H-bonds in the approximation of pair-wise interactions (see Techniques in the ESI †). The changes in the interaction energy vs. the distance between the proton donor and the proton acceptor for all H-bonds were compared for dimers extracted from the three polymorphs (if the selected type of bond was present for the given polymorph) (Fig. 4). The donoracceptor distances in experimental structures for all H-bonds are summarized in Table S2 in the ESI. † The bottom of the energy well for a selected hydrogen bond corresponds to the optimum distance between the donor and the acceptor, while any change results in the weakening of this interaction. The bond can still be compressed with gain in energy if the donor-acceptor distance is longer than the minimum value (the right side of the well), whereas the distances shorter than the optimum distance (the left side of the well) correspond to over-straining of the bond (Table S3 in    vector components. The highest overstrain for all three phases was found along the x vector (see Table S3 in the ESI †).
The profiles of some hydrogen bonds (#4 and #6) were almost identical for all three polymorphs (see Fig. 4(a) for H-bond #4). One can suppose that these hydrogen bonds can hardly have any significant effect on the pressure-induced transitions. In contrast some hydrogen bonds (#1 and #2, #3, #8 and #5) displayed significantly different energy wells, with notable differences in the minimum for different polymorphs, suggesting that it is the overstrain of these bonds that drives the microscopic mechanism of the structural transition resulting from the increasing pressure ( Fig. 4(b) for H-bond #5). All H-bond profiles are plotted in Fig. S4 in the ESI. † The two phase transitions in L-serine are significantly different in terms of hydrogen bond overstrain.
If L-serine I is compressed, the OHÁ Á ÁO hydrogen bond (#5) is overstrained along the a vector up to 11 kJ mol À1 . This is a large value, especially remembering that the energy difference between the polymorphs is several kJ mol À1 (ca.10 kJ mol À1 as calculated in this work). The chains of zwitter-ions linked by OHÁ Á ÁO H-bonds are directed along the a axis -the same direction as the direction in which the overstrain of this bond is observed. The overstrain along this hydrogen bond upon compression seems to be the reason for the I -II phase transition at the microscopic level. This bond acts as the ''spring'' which cannot be compressed further and therefore stretches back to the initial length at a sudden, triggering the phase transition. Such a model agrees well with the observed changes in cell parameters during I -II phase transition: the cell parameter a decreases with increasing pressure till the transition point, at which it instantly restores its initial value (Fig. S5 in the ESI †). At the same time the jump-wise change of the parameter c axis must be secondary. This also correlates well with the higher value of the bulk modulus of L-serine I in comparison to L-serine II, showing polymorph I being more stiff at the macroscopic level which can be considered as an indirect proof of its H-bond overstrain ( Fig. S3 and Table S2, ESI †).
The large hysteresis observed experimentally for I -II phase transition can also be interpreted at the microscopic level by the overstrain of the OHÁ Á ÁO H-bond upon increasing the pressureone can hardly expect that the overstrained H-bond will be formed also during the pressure decrease. If we start from L-serine II and stretch the system (modeling the pressure decrease), it would be very hard to rearrange the system in such a way, so that one ''spring'' (the OHÁ Á ÁO hydrogen bond) is pressed to the limit.
The situation with the phase transition from L-serine II to L-serine III is different. No significant overstrain could be observed for any of the pair-wise interactions ( Table S3 in the ESI †). This phase transition is a result of the co-operative effect when many types of hydrogen bonds are being changed to a similar extent, the bifurcated O-HÁ Á ÁO bond giving the maximum contribution to the energy gain. No change in the c parameter is observed upon II -III transition (Fig. S5 in the ESI †). The compression of the structure of L-serine II along the a axis still resembles a compression of a spring (but to a smaller extent, as compared to the I -II phase transition), until this spring relaxes at a sudden to restore its original length (Fig. S5 in the ESI †). However, in this case the limit of compression is reached not because of the overstrain of a selected individual H-bond, but because the neighboring molecules come so closely to each other that many contacts become unfavorable. This can explain the relatively small difference between the crystal structures of L-serine II and L-serine III and a smaller hysteresis of the phase transition between these two phases. The molecules approach each other as pressure is increased, so that at a certain point the formation of new H-bond becomes possible in the absence of any significant overstrain. Overall similarity of these structures, both showing no significant overstrain of H-bonds, was also showed by very close values for the bulk modulus, minimal volumes and energies in energy-volume profiles ( Fig. S3 and Table S2, ESI †). This model also agrees with the very small  changes in the values of b and c cell parameters, as well as those of the volume accompanying the phase transition from polymorph II to III. From aforementioned calculations one can expect that the nucleation barrier for the phase transition from I to II must be significantly higher than that for the transition from II to III. One can also expect that the rate of increasing pressure must have an effect on the phase transition from I to II. Another highpressure phase can be formed before the overstrain of the OHÁ Á ÁO hydrogen bond #5 (necessary to trigger the transformation into L-serine II) is achieved if pressure is increased slowly. It is also clear that polymorph III can be formed easily from polymorph II, but is not likely to form equally easily from polymorph I directly, even if pressure is increased very sharply immediately to the point of the second phase transition. These predictions of the modeling have been confirmed experimentally in a recent independent study. 26 Summing up, in this work a new approach to improving our understanding of the reasons and mechanisms of pressureinduced phase transitions was suggested and tested using L-serine. The structures were successfully optimized also out of the range of their stability. This enabled us to explore the difference between the hypothetical structures of different polymorphs at the same pressure. Despite the entropy effects were not considered the occurrence of the two phase transitions was predicted in the correct order; though the computed transition points were lower than the experimentally observed values, the difference between the points of the two phase transitions was close to that observed in experiments. It was shown that the PV term plays the decisive role at least in the first, I -II phase transition. The modeling of the effect of pressure on individual hydrogen bonds made it possible to rationalize a radical difference between the I -II and II -III phase transitions. The first one is triggered by a large overstrain of a selected intermolecular hydrogen bond OHÁ Á ÁO, which can explain the experimentally observed changes in cell parameters, and the large hysteresis. The second one seems to result from multiple small changes in many interactions without a significant overstrain of any selected hydrogen bond. The extension of this work using these results and the suggested computational approach in conjunction with new experimental data for the L-serine system will provide more understanding of phase transition L-serine I -L-serine IV and its co-existence with phases II and III.
The present work illustrates that relatively simple calculations complementing detailed experimental data can in fact give an insight into the macro-and micro-driving forces of pressureinduced phase transitions in hydrogen-bonded molecular crystals. The approaches used for model calculations are quite general and can be applied to pressure-induced transitions in other organic molecular crystals with different types of intermolecular hydrogen bonds.