A Short Review of Current Computational Concepts for High-Pressure Phase Transition Studies in Molecular Crystals

: High-pressure chemistry of organic compounds is a hot topic of modern chemistry. In this work, basic computational concepts for high-pressure phase transition studies in molecular crystals are described, showing their advantages and disadvantages. The interconnection of experimental and computational methods is highlighted, showing the importance of energy calculations in this field. Based on our deep understanding of methods’ limitations, we suggested the most convenient scheme for the computational study of high-pressure crystal structure changes. Finally, challenges and possible ways for progress in high-pressure phase transitions research of organic compounds are briefly discussed.


Introduction
High-pressure chemistry and particularly crystallography is developing fast in recent decades [1][2][3]. High-pressure effects are known to be studied originally and mainly by physicists and geologists [4,5]. Such significant interest in various minerals at high pressure and extreme temperatures is caused by questions that usually arise from geology-how substances act in Earth's crust and mantle and how do they transform further when temperature and pressure decreases. To describe mineral behavior, many theoretical and experimental works have been done, resulting in a deep understanding of the formation of our and other planets [6]. Among more practical results, one should point out equations of state (EoS) many of which were originally developed for minerals and combined in special software [7][8][9][10] but now used frequently in different fields [11][12][13].
Construction and functional materials are also known for both precise and vast investigations [14][15][16]. One can easily understand that many of the materials we use every day are exposed to relatively high pressures and sometimes high temperatures and, for sure, should be studied for possible undesirable phase transitions. Another application of high pressures is discovering new forms of different elements and substances, which may exhibit new valuable properties [17][18][19][20][21][22]. Stabilizing new forms obtained at high pressures at ambient conditions is an ambitious aim for applied science and industry-graphite transforming to diamonds is the most popular example of this possible process. The importance of finding and stabilizing new inorganic and metal-organic phases that may arise at high pressure is evident. These materials are almost in every part of different devices and constructions surrounding us.
Nevertheless, new forms of organic substances are also in demand by industry [23][24][25]. Production of different forms of active pharmaceutical ingredients is a hot topic of crystallography of organic substances and can be used in practice [25][26][27][28]. Among others, one can find the application of high-pressure being a non-trivial way for obtaining new polymorphs of desired substances [29][30][31][32]. Several groups are doing extensive work looking over many organic substances for new phases at high pressure. Nowadays, there are more than 1000 structures at high pressure according to the CCDC database, and its number is growing steadily [33].
Thus, one can understand that a high-pressure phase transition study in molecular crystals is an important field of modern chemistry and as any advanced direction should be studied both experimentally and computationally.

The Main Research Directions
High-pressure phase transitions can be studied mainly using two different concepts: 1. Experimental techniques usually provided in diamond anvil cells (DAC) by X-Ray diffraction and Raman spectroscopy [1] but not limited to these methods [34,35]. Significant progress in engineering made these experiments possible and even routine in some sense, but still, very time consuming and complicated. These experiments give atomic coordinates and unit cell parameters of molecular crystal and rarely some information on atom/functional groups motion, which can be interpreted in terms of energetic characteristics.  45], CPMD [46], CASTEP [47] and CRYSTAL [48] codes, but researcher's choice is not limited to abovementioned software. Reasonable choice of level of theory (LOT), which combines many parameters, gives very accurate energies of investigated systems but requires much more resources and time in comparison to FF methods. Here can be also mentioned many specific methods that can shed some light on different aspects of structures' nature or phase transitions, including electron/charge density analysis [49,50]. Nevertheless, these methods are relatively rarely used for molecular crystals at high-pressure conditions and have a specific discussion in literature [51,52], thus would not be described further in this work. We also do not describe molecular dynamics (MD) methods which are very useful [53][54][55], having not much experience in these methods. Modern machine learning and big data approaches are also not specified here, showing significant progress in the crystal structure and properties prediction, but not being applied to high-pressure phase transitions in molecular crystals [56][57][58][59][60].
Important to understand that both computational DFT and FF methods can barely be used without experiments being done at all because atom coordinates are needed for calculations [13,61]. Nevertheless, computational methods should be used when possible to understand and explain the nature of phase transitions of molecular crystals at high pressure because crystallographic and spectroscopic data rarely gives unambiguous reasons for phase transitions and can be interpreted differently by different scientific groups [61][62][63][64][65][66][67].

Force Field Methods
Force Field methods represent the functional form and parameter sets used to calculate the potential energy of a system of atoms or coarse-grained particles in molecular mechanics and molecular dynamics simulations. They are known widely due to their implementation in molecular dynamics for different systems, including biological ones [68,69]. Much less known but still frequently used "static" FF methods can evaluate lattice energies and some properties of molecular crystals [70][71][72][73][74]. The first successful attempt (accepted by the scientific community) was made by A. Gavezzotti, showing AA-CLP method for very fast evaluation of lattice energy and enthalpy of molecular crystals, based on atomic charges [36]. In this work [36], it has been postulated that "Structural papers with discussions of crystal packing involving user-selected atom-atom distances or geometries only should no longer be allowed in the literature" and gave the right direction of further research on molecular crystals and their forms. Later this method was substituted by more precise and sophisticated PIXEL [37], which used electron density for energy evaluation and Gaussian software [75] as a backend for this type of calculations. The energy of pair-wise interactions was described as a sum of different terms (Etot = Eele + Epol + Edis + Erep), which are given in analytical form [37]. Summation of these pair-wise interactions over a cluster of molecules in the crystal structure of radius equal to 20-50 A 3 (depending on the system) gives total lattice energy (Etot). This leads to a very fast estimation of energies and, what is more important, to the "chemical" partition of different energy terms. Definitely, this very clear energy concept together with DFT was applied to high-pressure research to find the reasons for phase transition at extreme conditions [67,76,77]. One of the pioneer works in this direction was done by Wood et al. [67] who combined experimental study with extensive computational work, where both DFT and FF methods were used. Relatively good convergence with DFT methods has been shown, taking into account some limitations applied to DFT optimization procedure. Moreover, energies of structures of L-serine polymorphs were reasonable and close for both DFT and PIXEL calculations. The most valuable feature is that the energy of intermolecular interactions can be monitored for experimental structures at different pressures, and it is possible to divide into chemically sensible terms and visualize using new CE17 implementation [40]. Another important feature is the possibility to monitor the energy of any pair-wise interaction in the crystal structure and thus follow phase transitions at high pressure at a molecular level. If one needs to understand deeply what the reason for a phase transition is and why some interactions or H-bonds are broken and new appeared, energy calculations are mandatory and should be done if atom coordinates are known for different crystal phases at high pressures.
Nevertheless, one should understand that despite FF methods are very fast and relatively easy to use, they should be checked carefully for specific tasks, one of which is high-pressure research. Energy is written as a summation of different terms with coefficients Etot = keleEele + kpolEpol + kdispEdisp + krepErep, where ki is parametrized on DFT energies or experimental enthalpies (depending on implementation) for molecular crystals at ambient conditions. H-bond length is also parametrized for ambient condition calculations, and cannot be switched off due to the very low precision of X-ray experimental measurements in defining hydrogen bond lengths. Summing up, an extensive benchmark on high-pressure data is needed to use FF methods unambiguously for high-pressure research in molecular crystals or should be critically checked by DFT calculations.

DFT Methods
In contrast to FF methods, DFT calculations for systems with periodic boundary conditions prove to give very accurate energies if multiple parameters are used correctly [78][79][80][81][82]. Discussing any phase transition, we assume that Gibbs energy change should be negative for spontaneous phase transition, which can be estimated from enthalpy and entropy terms. Entropy term is frequently neglected (especially if space group preserves after phase transition), supposing to have a small impact on Gibbs energy, while enthalpy can be calculated relatively easy [61,67,80,81,83]. Nevertheless, a lack of entropy calculations can lead to significant mistakes in the prediction of phase stability, which was shown for many inorganic materials [84,85]. The introduction of thermal effects in phase stability (DFT-QHA) increases computational costs drastically, but phonon effects are crucial to define accurate thermodynamic properties and Gibbs energies [86]. Calculation of all terms of Gibbs energy is strongly recommended to correctly predict stability and (subtle) phase transitions at high-pressure and temperature conditions [85,87,88]. Thus, the enthalpy diagram should be used as the first step of an investigation. Enthalpy is a sum of internal energy and PV term, where internal energy can be described as lattice energy-a sum of intermolecular and intramolecular (conformational) interactions. The thing that can be calculated using FF methods is intermolecular interactions only, while conformational (relaxation) energies should be calculated using DFT methods [80,81]. Constant-pressure (or, more precisely, fixed-stress) geometry optimization has to be carried out in DFT energy calculations to accurately define static pressure. Then, P-V-T EoS of the substance must be known via phonon dispersion calculations to define total pressure as a sum of static, zero-point and thermal contributions. Correct assessment of PV term at high-pressure and temperature conditions could be rather complex and expensive from a computational point of view and could strongly affect P-T location of phase transition boundaries [88]. Despite the abovementioned statements were supported mostly by examples from inorganic materials, we do not see any significant difference when these calculations applied to molecular crystals. Following Gibbs energy contributions can lead to a new understanding of reasons for phase transitions for any crystalline material. In some cases, it is possible to claim volume change (PV), energy (Einter or Eintra), enthalpy (H) or entropy (S) as a driving force for phase transition [61,67,81,89] DFT methods can be applied to experimental structures, obtained at some pressure, and all the above-mentioned parameters can be calculated. If one needs thermodynamic parameters for crystal structure at a pressure where no experimental data available, there are two main possibilities. The first is to calculate the equation of state for the phase and estimate the volume of the system, which will be fixed during ionic relaxation (structure optimization) to obtain energies. It is important to note that all calculations, in this case, should be performed for fixed unit cell volume. This technique is dependent on the accuracy of the equation of state, which needs at least 5-7 pressure-energy points within a pressure interval of a couple of GPa [90]. In this case, the quality of experimental data is very important because it influences the whole procedure of calculations. More convenient and independent from the number and quality of experimental data is the procedure with programmed pressure, e.g., as a starting structure can be taken one at ambient conditions and optimized to programmed pressure of, e.g., 3 GPa (e.g. PSTRESS option at VASP package). It was shown before that this kind of optimization to programmed pressure works accurately [61,91]. In this case, the volume is not fixed. In the case of DFT calculations enthalpy, the sum of inter-and intramolecular interactions (Ucrystal), entropy and PV term are obtained, which may give unambiguous answers with reasons for a pressure-induced phase transition. It is recommended to verify computational parameters (functional, dispersion correction scheme, k-points, and Ecutoff, etc.) on the structure, which is used as a guess structure for further calculations, usually ambient-pressure experimental data. One of the most important parameters is dispersion correction, which is crucial for accurate simulation of molecular crystals at ambient and non-ambient conditions. Nowadays, highly sophisticated density functionals are available to account for dispersive interactions in DFT (e.g., DFT-D3, B3LYP-D*, HF-3c, DFT-TS, M06, etc.) [92][93][94][95][96][97][98][99].
In relation to FF methods no "chemical" (electrostatic, polarization, dispersion, repulsion) terms are obtained in case of DFT calculations, but in common case gives more reliable energies. The advantages and disadvantages of DFT and FF methods applied to high-pressure researches are summarized in Table 1. 1 not mandatory if programmed pressure is used during optimization procedure but preferable for selection of parameters of optimization [61]. 2 Not implemented in PIXEL and CE17 software.

Combined Techniques
Recently we suggested the scheme where different phases are simulated not only in the intervals of their stability according to experimental data but also out of these intervals [61], and used it later for another system [91]. To the best of our knowledge, these are the first examples of such calculations for organic systems, studying intrinsic vs. extrinsic phase stability of polymorphs for molecular crystals. Nevertheless, it is widely used for inorganic systems, e.g., [85,100]. This helps to understand what would happen to the polymorph structure if no phase transition occurs. Such simulation provides structural and energy data to direct comparison of all phases at the same conditions, even if some phases are not found at these pressures experimentally (Figure 1).

Figure 1.
Schematic representation of calculated energy parameters (∆Ucrystal, ∆H, etc.) for different phases in the whole pressure range, even at pressures where a specific phase is not found experimentally. Solid dots filled with color, starting experimental structure for DFT optimization; solid empty dots, DFT optimized structure at fixed pressure as observed experimentally; dashed dots, DFT optimized structure at fixed pressure out of experimental pressure range. Blue, phase I; green, phase II; red, phase III. Dashed squares show pressure-energy conditions close to high-pressure phase transitions.
In previous work, we also simulated pair-wise interactions using DFT gas-phase calculations with fixed geometries after solid state optimization. This simulation shows energy changes when distance decreases at pressure, building energy well for each H-bond [61]. This kind of energy well simulation is definitely not perfect due to the absence of molecular conformation change because of the molecular surrounding. This possibly could be also simulated using FF methods as described before and was done in other works [67,80,81,101].
Summing up, we suggest that the following scheme should be applied for high-pressure phase transition research of molecular crystals (Scheme 1).

Scheme 1.
A suggested scheme for the computational study of high-pressure phase transitions of molecular crystals. All steps are calculated for structures both in and out of structural stability regions.
Finally, we would like to point out that scrupulous checks for the validity of all parameters and methods should be done for every particular system. If this is done correctly, one can assume this computational scheme is a numerical experiment aimed at understanding the reasons and mechanisms of phase transitions of molecular crystals at high pressure.
In this review, we put aside many questions related to high-pressure polymorphism of organic compounds: kinetic barriers for phase transitions which play an important role [102][103][104][105], temperature corrections for DFT methods [106][107][108], prediction of new phases (which in future can significantly decrease amount of experimental work) [109][110][111], etc. Nevertheless, as it is shown above, current computational methods can bring out many more answers than pure experimental results.

Prospects
To the best of our knowledge, there is no specific method or software for computationally studying high-pressure phase transitions in molecular crystals. Still, wise usage of already developed methods and their implementations with appropriate validation of different parameters can unveil reasons and even some mechanisms for phase transitions. We hope that in the future all experimental works be complemented with computational techniques. This may bring out true reasons and a deep understanding of already observed phase transitions and the prediction of new polymorphs. Finally, taking into account progress in crystal structure prediction [109] [112], one can expect the prediction of new structures of molecular crystals at extreme conditions. Nevertheless, we are sure that the concept shown in this work gives the most convenient and all-round scheme for high-pressure research of organic compounds using computational methods and can be improved by obtaining different more specific properties: mechanical, electronic, optical, etc. Finally, we do not see any issues that can prevent usage of the proposed scheme for metal-organic and inorganic materials if proper parameters of the simulation are chosen.