Predicting the phase diagram of solid carbon dioxide at high pressure from first principles

The physics of solid carbon dioxide and its different polymorphs are not only of great practical and fundamental interest but also of considerable importance to terrestrial and planetary chemistry. Despite decades of computer simulations, the atomic-level structures of solid carbon dioxide polymorphs are still far from well understood and the phase diagrams of solid carbon dioxide predicted by traditional empirical force fields or density-functional theory are still challenged by their accuracies in describing the hydrogen bonding and van-der-Waals interactions. Especially the “intermediate state” solid carbon dioxide phase II, separating the most stable molecular phases from the intermediate forms, has not been demonstrated accurately and is the matter of a long standing debate. Here, we introduce a general ab initio electron-correlated method that can predict the Gibbs free energies and thus the phase diagrams of carbon dioxide phases I, II and III, using the high-level second-order Møller-Plesset perturbation (MP2) theory at high pressures and finite temperatures. The predicted crystal structures, phase transitions, and Raman spectra are in excellent agreement with the experiments. The proposed model not only reestablishes the position of solid carbon dioxide in phase diagram but also holds exceptional promise in assisting experimental studies of exploring new phases of molecular crystals with potentially important applications.


INTRODUCTION
Carbon dioxide (CO 2 ), the major ingredient of the atmospheres of terrestrial planets, such as Mars and Venus, 1,2 is commonly found in ice form in planets and asteroids. 3 Solid CO 2 has the simple structure (with one carbon atom and two oxygen atoms), but it is more likely to transforms into a series of various stable, instable or metastable polymorphs at high pressures and temperatures. 4,5 Such polymorphs, occupying different positions in the phase diagram, have different symmetries, kinetic energies, intermolecular interactions, and crystal structures. [6][7][8] The experimental phase diagrams of CO 2 were proposed by Litasov et al., 2 Giordano et al., 9 and Iota et al., 10 where phase I occupies the low pressure and low temperature region, phase II and phase III locate at the high pressure, high temperature area and high pressure, low temperature area, respectively. Although intense experimental efforts have been devoted to identify the crystal structures and phase diagram for decades, the transition boundaries for CO 2 phases as well as their vibrational spectra are rather complicated and not yet well understood. For example, CO 2 -cubic phase I (Pa3) transforms to the orthorhombic phase III (Cmca) at [10][11][12][13][14][15][16][17][18][19][20][21][22] at temperatures ranging from 40 K to room temperature. The reported transition pressures range from 10 GPa, 14 12 GPa 13 to 18 GPa 11 at room temperature. Further subsequent heating at 20 GPa, CO 2 phase III transforms to a pseudotetragonal phase II (P4 2 /mnm). 16 The crystal structures of phases I, II, and III of CO 2 are shown in Fig. 1, where the phase I and phase III structures are closely related to each other and their phase transition can occur ideally by the molecular rotation from the cubic structure until the molecules become parallel to one of the base planes and the unit cell is deformed from a cube to a rectangular parallelepiped. 17 The unit cells of phase I and phase III contain four molecules, respectively, while the unit cell of phase II has two molecules.
However, solid CO 2 phase II, one of the most controversial structures, is far from the cognition and exploration. Prior to 2002, Yoo et al. 16 have determined the crystal structure of CO 2 phase II in situ at high pressures and high temperatures by angle resolved X-ray diffraction. Their data indicate that the carbon atoms in CO 2 phase II are pseudo-six-fold coordinated with oxygen atoms, whereas oxygen atoms are threefold coordinated with carbon atoms. They reported that CO 2 phase II had the P4 2 /mnm symmetry with the intermolecular distance almost less than twice the intramolecular distance. Based on the elongated intramolecular bond (C=O) distance (over 1.30 Å), and collapsed intermolecular (C-O) distance (below 2.38 Å), Yoo et al. 16 concluded that CO 2 -II was an intermediate between molecular and nonmolecular solids. However, Datchi et al. 17 and Bonev et al. 18 rejected this conclusion and measured a different intramolecular C=O bond length (1.14 Å), which suggested that phase II was a molecular crystal rather than an intermediate state. Another argument is about the transformation boundary between phase II and phase III in CO 2 phase diagram. Phase II was found by heating phase III above 16 GPa and 500 K, but the III-to-II transition was suggested to be irreversible, which seems to suggest that CO 2 -III is a metastable phase and the observed III-to-II transformation boundary is controversial. 19 There are two main controversial hypotheses, one is that the III-to-II transformation boundary is a kinetic barrier rather than a phase boundary and the other supports the III-to-II transformation boundary as a true phase boundary and assumes the existence of two distinct triple points separated by less than the experimental precision (5 K and 1 GPa). 20 Therefore, further additional work is needed to clarify above problematic aspects on solid carbon dioxide phase II on experimental or theoretical basis. 21,22 The large hysteresis of the phase diagram and the intermediate state of CO 2 -phase II make it difficult to characterize them by experimental means alone; the assistance of accurate computational calculations is essential. Pioneering ab initio methods have been previously performed on the calculations of CO 2 crystal structures and phase transitions. [23][24][25] The widely used density functional theory (DFT) 26 can give reasonable results within the valid period of prediction but is still challenged by the lack of dispersion interactions and the difficulty in improving the accuracy systematically. For example, Bonev et al. 18 calculated the phase transitions of CO 2 phase I, phase II, and phase III with DFT theory, which ruled out the experimentally derived hypothesis that both phase II and phase III were high-strength materials, and the positions of phase II and phase III were reversely compared with the experiment in the phase diagram. The most effective systematic approximation of ab initio molecular orbital theory is the second-order Møller-Plesset perturbation (MP2) theory, whose strengths spans three orders of magnitude and is able to capture the covalent, ionic, hydrogen-bond, and dispersion interactions accurately and can reproduce the crystal parameters quantitatively. In our previous work, we used the MP2 theory along with the embedded fragment method to predict the phase transition of CO 2 from phase I to phase III, 27 where the predicted transition took place at 13 GPa and 0 K and agreed well with the experiment. The embedded fragment method, dividing the system into fragments (including monomers and dimers), can reduce the tremendous computational cost of MP2 theory and be utilized for the free energy calculation of a molecular crystal and thus the phase transitions accurately.
In this article, we report ab initio calculations of crystal structures, Raman spectra, free energies of CO 2 phase I, phase II, and phase III at high pressures and high temperatures and thus of phase diagram from first principles. We use MP2 with the Møller-Plesset partitioning of the exact electronic Hamiltonian, which is able to describe dispersion interactions accurately when predicting the crystal structures of solid CO 2 . Furthermore, Raman spectrum is a distinct fingerprint for a particular molecule, and can be used to rapidly identify the materials, or distinguish them from others. The proposed MP2 theory along with the embedded fragment method reproduces the Raman spectra of CO 2 phase II, which was observed by Datchi et al. 17 and has not been supported by any theory. We demonstrate that CO 2 phase II is a molecular structure rather than an intermediate state. Based on the electron-correlated calculation, we accurately reproduce the volume-pressure relationships, enthalpies, phonon density of states, Raman spectra and then predict the phase transitions of solid CO 2 phase I, II, and III via the calculations of Gibbs free energies from 0 to 1000 K and 0 to 20 GPa. We support the hypothesis that the phase III-to-II transformation boundary is a true phase boundary in the carbon dioxide phase diagram.

Equations of state
We use our in-house parallel-execution program (based on the MP2 theory and the embedded-fragment method) to perform the structural optimization, vibrational analysis, and Gibbs free energy calculations. The computed pressure-volume curves at T = 0 for three proposed structures of CO 2 phases I, II, and III are plotted in Fig. 2, along with the available experimental data (solid dots). At the given pressure range, phase II has a very close volume with phase III, both of which are obviously smaller than that of phase I. The MP2 calculated P-V curves, agreeing well with the experimental data, 7,16 record ca. 2 and 0.5% volume reductions upon transitions from phase I to phase III and phase III to phase II, respectively. The visible volume collapse shown in the shade of Fig. 2 is also consistent with the experimental data. 16 The decreases of volumes of phases I, II, and III imply that CO 2 crystal will become progressively stiffer with the increasing of the pressure. As a result, electrons localized within intramolecular bonds become less stable as pressure increases and the intermolecular potential becomes repulsive. The volume relations for the three phases are: phase I > phase III > phase II. From previous experimental work, the transition from phase I to phase III is more favorable than to phase II, and the transitions from phase II to phase III or to phase I is prohibited during the compressing. Such results suggest that what is considered phase III in some experiments is actually phase II, when compressing phase I. 16,18,19,22 Figure 3 shows the MP2/aug-cc-pVDZ calculated lattice constants (a = b and c) of carbon dioxide phase II compared with experimental data. The calculated lattice constants are consistent with the experimental results.
The big debate of CO 2 phase II is whether it is an intermediate state or a molecular crystal, which depends on the intramolecular bond (C=O) distance and the intermolecular bond (C…O) distance. Yoo et al. 16 measured the intermolecular distance almost less than twice the intramolecular distance (over 1.30 Å), while Datchi et al. 17 reported the intermolecular distance twice more than the intramolecular distance (1.14 Å). The former result tends to prove that the CO 2 phase II is an intermediate state, while the later leads to a molecular crystal of phase II. In the present study, the C=O bond lengths are calculated by the high-level MP2 theory, which is on the order of 1.1 Å and consistent with the work given by Datchi et al. 17 As shown in Fig. 4, the bond lengths of all three phases decrease nearly linearly as the pressure increases. From 10 to 30 GPa, the intramolecular C=O bond lengths of phase II, presenting about 15% drop after full structural optimization, start from 1.174 to 1.164 Å, which is consistent with the work of Datchi et al. 17 Therefore, the MP2 calculation, matching the work proposed by Datchi et al., 17 indicates that CO 2 phase II is a molecular crystal rather than an intermediate state. Besides, with the increasing of pressure, the C=O bond lengths of phases II and III have the similar decreasing tendency and are both obviously larger than that of phase I. Such relation is consistent with the pressure-volume figure in Fig. 2.
Vibrational spectra Raman spectrum is a distinct chemical fingerprint for a particular molecule or material, and can be used to rapidly identify the materials, or to distinguish them from others. Our MP2 calculations reproduce the Raman spectra of carbon dioxide phases with remarkable accuracy. Figure 5 compares the calculated and observed Raman spectra in the librational region of carbon dioxide phase I at 11.7 GPa (Fig. 5a), phase II at 26 GPa (Fig. 5b), and phase III at 26 GPa (Fig. 5c), respectively. The observed curves in Fig. 5a-c are taken from the work done by Olijnyk et al. 15 (phase I and phase III) and Datchi et al. 17 (phase II). The calculated Raman spectra of phase I has three Raman bands (at 127 cm −1 , 179 cm −1 and 273 cm −1 ), while phase II (at 265 cm −1 and 350 cm −1 ) and phase III (at 305 cm −1 and 367 cm −1 ) have two Raman bands in the librational regions, respectively. The calculated Raman spectra are consistent with Datchi et al.'s 17 experimental result in both the numbers and positions of Raman peaks, which establishes the correctness of our CO 2 phase calculation at finite pressures. Figure 6 summarizes the calculated and observed pressuredependence of the frequencies of the Raman bands of carbon dioxide phase II. The calculated vibrational frequencies (red curves) of phase II are in line with the observed data (black dots), which are taken from Iota et al. 10 The MP2 calculated Raman bands of carbon dioxide phase I and phase III can be found in previous work by Li et al. 27 Gibbs free energy Figure 7 shows the MP2/aug-cc-pVDZ calculated Gibbs free energy differences between carbon dioxides phases II-III and phases I-III at different pressures (Fig. 7a) and temperatures (Fig.  7b), where the positive value means that phase III is more stable. As shown in Fig. 7a, the transition temperature from phases III to II increases slightly with the increasing of pressure, while in Fig. 7b the transition pressure from phase III to I decreases slightly with the increasing of temperature. Such small pressure and temperature differences in Fig. 7a, b are in line with the experiments, 13 and expose the limitations of the DFT and empirical calculations. 12,18 Figure 8 plots the MP2/aug-cc-pVDZ calculated threedimensional Gibbs free energy surfaces as functions of pressure  Fig. 2 Pressure-volume relationship of solid CO 2 , for phase I (black curve), phase II (red curve), and phase III (blue curve) calculated by MP2/aug-cc-pVDZ. The black dots are the experimental data from refs. 7,16 The phase III and phase II exhibit the volume collapse between 11-13 GPa (shaded area) as indication of CO 2 phase transitions and temperature for carbon dioxide phases I, II and III, respectively, in which the intersections of two energy surfaces indicate the transition boundaries of the two phases. The lower the free energy surface, the more stable the structure.
Phase diagram Given the accurate structures and Gibbs free energies of carbon dioxide phases I, II and III, predicted by MP2/aug-cc-pVDZ, Fig. 9 shows the calculated phase diagram, along with the experimental phase boundaries, in the pressure range of 0 to 20 GPa and the temperature range of 0 to 1000 K. MP2 calculation predicts that the phase I-III transition pressure is 13 GPa at 0 K (or a few GPa less at higher temperatures), and the phase II-III transition temperature is slightly increased with the increasing of pressure, such as 12 GPa at 515 K, 13 GPa at 531 K, 15 GPa at 552 K, etc. The transition curve between phases II and III is considered as a kinetic line. As shown in Fig. 9, phase I is a low temperature and low pressure structure, while phase III is more stable than phase II at lower temperature and high pressure, which means that phase III can be reached by the low-temperature compression of phase I, and phase II can be obtained by the high-pressure heating of phase III. The present MP2/aug-cc-pVDZ calculation match the experimental phase diagram of carbon dioxide very well and expose the limitations of the DFT calculations and empirical calculations, where the positions of phase II and phase III were predicted reversed or dislocated in the phase diagram. 12,18,28,29 In conclusion, we have examined the CO 2 phase transitions at high pressure and finite temperature from first principles based on an embedded fragment MP2 method, which, therefore, can calculate the enthalpy, Raman frequencies and Gibbs free energies of solid CO 2 phases I, II and III, thereby allowing an ab initio Fig. 4 The pressure dependence of the MP2/aug-cc-pVDZ calculated intramolecular C=O bond lengths of CO 2 phases I (black curve), II (red curve) and III (blue curve) determination of these phase transitions. With this method, we have reproduced quantitatively the lattice constants, equation of state, and the Raman bands of phases I, II and III, which agree well with experiments. We conclude that the carbon dioxide phase II is a molecular crystal rather than an intermediate state, and the phase III is more favorable at low temperature and high pressure region while the phase II occupies the high temperature and high pressure area in the phase diagram. The predicted phase transition boundaries of solid carbon dioxide phases I, II, and III are reasonably in line with experimental work, especially the transitions between phase III and phase II. The proposed study can not only be very helpful to reestablish the phase diagram of solid carbon dioxide but also motivate the possibility of exploration of new phases with important applications.

Energy calculation
For the molecular crystal system, the internal energy of a unit cell can be calculated by: where n is the three-integer index of a unit cell, E i(n) is the calculated MP2 energy of the ith molecule in the nth unit cell, and E i(0)j(n) is the calculated MP2 energy of the dimer for the ith molecule in the central unit cell and the jth molecule in the nth unit cell. In this work, we take a 3 × 3 × 3 supercell of the crystal. The first term in Eq. (1) gives all the single molecule energy in the central unit cell. The second term gives the interaction of two-body system whose distance is shorter than a given cutoff distance λ in the central unit cell. The third term gives the interactions between two molecules in the central unit cell and the nth unit cell that has a shorter distance than the cutoff distance λ (λ was set to 5 Å in this study). The first three terms, referring to the short-range interactions, are calculated at the MP2/aug-cc-pVDZ level in the electrostatic field of the rest crystal represented by the electrostatic potential (ESP) charges fitted at the HF/aug-cc-pVDZ level. The long-range interactions of the dimer, in which the distance between two molecular is larger than λ, are approximately treated as charge-charge Coulomb interactions. We take into account the background charges in the 11 × 11 × 11 supercell. The last term E LR gives the long-range electrostatic interactions in a 41 × 41 × 41 supercell.
Considering the effect of external pressure, 30 the enthalpy H e per unit cell can be calculated by: where P is the external pressure and V is the unit cell volume. The Gibbs free energy dependence on temperature and pressure of a unit cell G e is calculated by: where U v and S v are the zero-point vibrational energy and entropy per unit cell at temperature T and pressure P. For molecular crystals, the zero-point vibrational energy U v and the entropy S v are obtained by Eqs. (4) and (5) with the harmonic approximation: where ω nk is the frequency of the phonon with lattice vector k, β = 1/k 0 T, and k 0 is the Boltzmann constant. The product over k must be taken over all K evenly spaced grid points of k in the reciprocal unit cell. In this study, the k-grid of 21 × 21 × 21 has been used (K = 9261).

Embedded-fragment method
Although the ab initio methods perform well on property calculation, they cannot be efficiently applied to large systems. The embedded fragment method 27 treats a large system as many small subsystems known as fragments. The calculations will only need to be done on these subsystems, and then we can obtain properties of the entire large system by bringing the subsystems' results together. Varying in the way to form fragments and deal with the interactions between fragments, there are many schemes of the embedded-fragment method. The fragment-based quantum mechanical method can only be applied to electron-localized systems, such as liquid water, salt solutions, proteins, DNAs and ionic liquids, but not to electron-delocalized systems, such as metals. [31][32][33][34][35] Here, we utilize the binary interaction method, which was developed by Hirata and coworkers 30,36,37 for molecular systems and periodic crystalline systems. For molecular crystal systems, each individual molecular is treated as a fragment. The interactions between these fragments are handled by a many-body expansion. Since the many-body expansion can converge easily, the binary interaction method neglects higher-order terms and only includes one-body and two-body contributions. The interaction energy between two fragments, whose distance is shorter than the cutoff distance, is calculated by quantum mechanics (QM). The interactions between two long-range fragments are treated as pairwise charge-charge Coulomb interactions. [32][33][34][35][38][39][40] To deal with the environmental effects, all QM calculations are embedded in an electrostatic field of point charges.

Structure optimization and properties
The enthalpy H e and the internal energy of a unit cell E e can be analytically differentiated with atomic positions and lattice constant. 36,37 We can perform the crystal optimization to relax the molecular crystal and obtain the most stable structure at a specific pressure by minimizing the enthalpy. In this work, we introduce the quasi-Newton algorithm 41 for the crystal structure optimizations. The approximate Hessian matrix was updated by the BFGS procedure. 42 The convergence criterion for the maximum gradient was set to 0.001 Hartree/Bohr.

Raman spectra simulation
The dynamical force constant matrix of a periodic system can be expressed as where m A and m B denote the mass of atom A and atom B, respectively, k represents a given point in the Brillouin zone, and H(r A,0 , r B,n ) is the secondorder derivative of the total energy per unit cell with respect to atom A in the 0th cell and atom B in the nth cell at the equilibrium geometry. 39 The number of neighboring unit cells for QM treatment was truncated at S = 1, and the number of k-points was set to 21 in each of x, y, z dimensions. Once the dynamical matrix of a periodic system is obtained, one can calculate the vibrational frequencies as well as the corresponding normal modes by diagonalizing the matrix D(r A ,r B ,k). In simulating the Raman spectra, only the zone-center (k = 0) vibrations have nonzero intensities and thus are Raman-active. Therefore, the force-constant matrix D(0) was used for the vibrational frequency calculation of the Raman spectra. The Raman intensity (R n0 ) for the fundamental transition of the mode in the nth phonon branch with wave vector k = 0 can be obtained by the following derivatives 39 where Q n0 is the corresponding normal mode. The dipole moment derivative ∂μ p =∂Q n0 and the polarizability derivative ∂α pq =∂Q n0 in the central unit cell can be derived based on the embedded-fragment quantum mechanical method. 39

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. Fig. 8 Gibbs free energy surfaces of carbon dioxides phase I (black surface), phase II (red surface) and phase III (blue surface) as functions of temperature and pressure, calculated by MP2/aug-cc-pVDZ. The intersection of black and blue surfaces is the phase transition boundary between phase I and phase III, while the intersection of blue and red surfaces is the phase transition boundary of phase II and phase III Fig. 9 Phase diagram of carbon dioxide. The red dots and blue line are the calculated transition boundary between phases II-III and phases I-III, respectively. The experimental data are taken from Giordano et al., 9,43 Iota et al., 10 and Litasov et al. 2