Predicting finite-temperature properties of crystalline carbon dioxide from first principles with quantitative accuracy

The temperature-dependence of the crystalline carbon dioxide (phase I) structure, thermodynamics, and mechanical properties are predicted in excellent agreement with experiment over a 200 K temperature range using high-level electronic structure calculations.


Method for computing the thermal expansion
To determine the structure of phase I carbon dioxide at a given temperature T and pressure P , one minimizes the Gibbs free energy G with respect to both the atomic positions in the unit cell and the unit cell parameters (a, b, c, α, β, γ). G(T, P ) = U el + P V + F vib (T ) (1) where U el is the internal electronic energy, P V is the pressure-volume contribution, and F vib represents the Helmholtz vibrational free energy contribution. This latter contribution is expressed in terms of the harmonic phonons ω k,i evaluated at multiple k-points in reciprocal space: where n is the number of unit cells in the supercell approximation and N A is Avogadro's number. Both U el and F vib are computed using the fragment-based hybrid many-body interaction (HMBI) model, [1][2][3][4] which allows one to perform high-level MP2 and even coupled cluster calculations on periodic systems like molecular crystals with reasonable computational cost. HMBI decomposes the intermolecular interactions in a crystal according to a many-body expansion, The important intramolecular (1-body) and short-range pairwise (SR 2-body) interactions are treated with quantum mechanics (QM), while the generally weaker long-range pairwise (LR 2-body) and manybody contributions are approximated with a polarizable molecular mechanics (MM) force field. In practice, the short-range 2-body QM treatment includes interactions involving molecules in the unit cell and those in nearby periodic image cells, while the MM terms are used to capture the long-range periodicity of the crystal via Ewald summation. The phonons at a given k-point are evaluated by constructing and then diagonalizing the massweighted dynamical matrix using the supercell approach, D α,β l, l ; k = 1 √ M l M l κ ∂ 2 U El ∂α l(0) ∂β l (κ) exp (−2πik · δx(0, κ)) where ∂ 2 U El ∂α l(0) ∂β l (κ) are individual elements in the supercell Hessian 5,6 involving the α-coordinate of atom l in the central unit cell (index 0) and the β-coordinate of atom l in the unit cell with index κ, and δx(0, κ) is the displacement vector between these two atoms. The M 's are the atomic masses. The exponential introduces the phase shift in the harmonic motions of the periodic image atoms relative to the atoms in the central unit cell.
Unlike conventional periodic boundary condition models, the fragment approach used in HMBI enables these equations to be evaluated for the supercell Hessian with minimal additional effort beyond the normal unit cell Hessian. 6,7 All the necessary supercell QM contributions to the force constants in Eq 4 can be constructed from contributions already available from the standard unit cell Hessian. The only new contribution required is the full supercell Hessian at the MM level. This allows one to capture phonon dispersion by using a large supercell and sampling many k points with comparatively low additional computation cost.
The treatment of F vib is potentially very computational demanding, since it would normally require many cycles of geometry optimization and a phonon calculation for each update in the unit cell parameters. Instead, we approximate the phonons for a given crystal volume using the quasiharmonic approximation (QHA). The QHA relates the i-th phonon frequency ω i at a given volume V to a reference frequency ω i,ref obtained at some reference volume V ref via the Grüneisen parameter for that phonon mode γ i , where Overall, then, a reference crystal structure, unit cell volume, and phonon frequencies are obtained at zero temperature and pressure by minimizing U el , and the harmonic phonon modes are computed using lattice dynamics. Two additional rounds of geometry optimization and phonon calculation are performed with fixed unit cell parameters at cell volumes which are slightly larger and smaller than V ref . Here, the cell volume was expanded/contracted by 10Å 3 , or roughly 5%. The Grüneisen parameter for each of the 3N vibrational modes is computed via finite difference using a linear model. 8 The phonon frequencies at the different volumes were matched via maximum overlap of the phonon eigenmodes.
With the reference frequencies, reference volume, and Grüneisen parameters obtained from these three optimizations and lattice dynamics phonon calculations, one can evaluate F vib at any given temperature and unit cell volume. That allows one to find the crystal structure at any temperature and pressure by minimizing the Gibbs free energy (Eq 1). Other properties like the heat capacity, enthalpy, and entropy can be calculated via standard partition function expressions based on the optimized structures and quasiharmonic phonon frequencies.
As described in the main paper, the QM contributions in the HMBI model were computed with MP2 or CCSD(T). At the MP2 level, analytic gradients were used, while the Hessian elements were calculated via finite difference of the gradients. The lattice dynamics phonon were performed on a 3 × 3 × 3 supercell using a 3 × 3 × 3 Monkhorst-Pack grid of k points. At the CCSD(T) level, analytic gradients are computationally expensive and are not widely available. Accordingly, the gradients were evaluated via finite difference of the energies. However, computing CCSD(T) second nuclear derivatives via finite difference of the energies introduced too much numerical noise into the two-body Hessian contributions, which made it difficult to obtain useful phonon frequencies. Therefore, for the CCSD(T) calculations, F vib was evaluated using MP2/CBS frequencies and Grüneisen parameters.

Method for computing thermodynamic parameters
Using the crystal structures optimized at a given temperature according to Section 1 and the quasiharmonic estimates for the phonon frequencies, one can compute the various thermodynamic properties reported in the paper as follows: Enthalpy of sublimation: The sublimation enthalpy is given as the difference between the enthalpy of the gas and 1/4 the enthalpy of the crystal (which has 4 molecules in the unit cell).
The enthalpy of the crystal is computed as the sum of the electronic energy of the crystal U el,crystal plus a pressure-volume term P V and a vibrational internal energy contribution E vib,crystal : where the vibrational internal energy contribution is given as: The gas phase was modeled using standard ideal gas partition functions. Of course, the gas may exhibit deviations from ideality at the low temperatures considered here, but they are hopefully not too large. The enthalpy of the gas is given as the sum of electronic (U el,gas ), translational (3/2RT ), rotational (RT ), vibrational (E vib ) terms plus an additional factor of RT from the P V term.
where E vib is defined as: Entropy of sublimation: The entropy of sublimation was computed as, Once again, the gas was modeled using standard ideal gas partition function contributions from translation, rotation, and vibration: S gas = S trans,gas + S rot,gas + S vib,gas (13) where the terms are defined as: Note that the symmetry factor σ in S rot is 2 for carbon dioxide.
Heat capacity: The isochoric heat capacity of the crystal was computed using the standard harmonic oscillator expression, 3 Procedure for obtaining empirical entropy data The empirical entropies were derived from experimental data according to: This expression relates the sublimation entropy at a given temperature to the experimental sublimation entropy at the sublimation point (i.e. ∆S sub (194.7 K) = 129.62 J/mol K), as reported by Giauque and Egan. 9 By construction, Eq 19 exactly reproduces this experimental value at 194.7 K.
The crystal contribution was computed and integration of the isobaric heat capacity data of Giauque and Egan. 9 The highest temperature C p reported by Giauque and Egan is at 189.78 K. The C p at 194.7 K was linearly extrapolated to 55.51 J/mol K from the five highest temperature values reported by Giauque and Egan. A smooth integrable function for the heat capacities was created using a cubic spline in Mathematica 9.

Equation of state fitting procedures
Pressure versus volume curves were generated at several different temperatures. At a given temperature and pressure, the crystal structure was optimized by minimizing the Gibbs free energy to obtain the unit cell volume. This procedure was typically repeated for 10 (21 at 296 K) different pressures at each temperature. Sample MP2/CBS data at 296 K is plotted in Figure S1. The full set of predicted P -V data is provided in Table S3.  Figure S1: Comparison of the predicted and experimental pressure versus volume data at 296 K. Note that CO 2 transitions to phase III above ∼10 GPa, as seen by the sudden drop in the experimental volumes at higher pressures.
Next, the P -V data was fitted to the Vinet equation of state, 12 whereṼ = (V /V 0 ) 1/3 . The constants V 0 , B 0 , and B 0 were obtained via non-linear least squares fitting.
To explore the multitude of least squares solutions, ∼1000 fits with randomly sampled initial guess parameters spanning two-orders of magnitude were performed. For every method/basis set, the best P -V fits exhibited residuals orders of magnitude smaller than the other fits identified, and that best fit is reported here. Note that for the Birch-Murnaghan equation of state, however, the same was not true. The Birch-Murnaghan fits proved ill-constrained, with a large number of very different parameter sets (e.g. with B 0 and B 0 values varying by up to an order of magnitude) producing fits with nearly identical residuals. Fixing the volume and linearizing the Birch-Murnaghan equation helps, but the results were extremely sensitive to the V 0 used.
Exploratory fits to the experimental data proved less well-constrained than those to the theoretical predictions due to the increased noise in the P -V data, but the parameters extracted from Monte Carlo sampling of fits to the Vinet equation exhibit much less variability than those from the Birch-Murnaghan equation. For these reasons of numerical stability, we adopted the Vinet equation.
The zero-pressure volume V 0 was treated as an adjustable fitting parameter in the Vinet equation of state. However, one could instead extract V 0 directly from the zero-pressure geometry optimization of the crystal as the given temperature, and fit only B 0 and B 0 . Table S2 compares the results from both procedures. In most cases, the resulting V 0 , B 0 , and B 0 values are nearly identical, with variations in the three fitted parameters of ∼0.1 or less. The largest variations between the two fitting procedures occur at 296 K, but these variations are smaller than the typical experimental uncertainties 13 in the fit parameters at the same temperature.
We also predicted the bulk modulus at 190 K using CCSD(T)/CBS. As shown in Table S2, the resulting equation of state parameters differ only marginally from the MP2/CBS ones. Given the high computational expense of performing 10 separate CCSD(T) geometry optimizations at various pressures for each temperature and the small difference between the coupled cluster and MP2 results, CCSD(T) calculations were not performed at the other temperatures. Table S2: Comparison between Vinet EOS fits with and without constrained volumes. In the constrained volume case, the parameter V 0 was fixed to match the geometry optimized unit cell volume. In the unconstrained case, V 0 was treated as a fitting parameter.

Unconstrained Volume
Constrained Volume Temperature Method The crystal was not bound at zero pressure and 296 K with MP2/aDZ, so V 0 could not be obtained. b The crystal was not bound at zero pressure and 296 K with the force field, and the Vinet EOS fits proved numerically unstable. As discussed in the main article and Section 1 above, lattice dynamics calculations were performed on a 3 × 3 × 3 k-point grid to capture phonon dispersion. As shown in Figure S2, this proves important for capturing the proper slope of the thermal expansion. A Γ-point only phonon treatment underestimates the thermal expansion at higher temperatures. 6 Intermolecular Amoeba force field parameters for CO 2 The HMBI model requires a polarizable force field for the long-range pairwise and many-body intermolecular interactions. These parameters were generated for a CO 2 using version 1.1.3 of the Poltype utility. 14 This utility first optimizes the molecular geometry at the MP2/6-31G* level of theory. Initial atomic multipole moments are calculated using Gaussian distributed multipole analysis on an MP2/6-311G** electron density. These multipole moments are subsequently refined by fitting against electrostatic potentials at various points in space computed from MP2/6-311++G(2d,2p

Tables of predicted properties
The data tables below list the predicted properties used to generate the figures in the main paper.

Temperature-dependence of the crystal structure
The space group for phase I CO 2 is P a3. This space group is cubic so all lattice length are the same and all lattice angles are 90 • . There are 4 CO 2 molecules in the unit cell. The is one carbon atom and one oxygen atom in the asymmetric unit. The carbon in the asymmetric unit is fixed at the origin. The carbon-oxygen double bond(C=O) forms a 45 • angle to each of the lattice angles when projected into the sides of the unit cell. The only adjustable parameters are the length of the carbon-oxygen double bond (C=0) and the length of the unit cell (a). These values at each level of theory and basis set are presented in Table S4. The location of the remaining atoms can be determined using space group operators. The crystal was not bound at zero pressure and 296 K with MP2/aDZ or the force field. b The CCSD(T)/CBS structure at 296 K was not computed.

Bulk moduli
Tables S11 and S12 present the experimental and predicted bulk modulus values used to generate Figure 5 in the main article. Note that only DFT calculations which used a van der Waals dispersion correction (specifically, the D3(BJD) correction) were included in Figure 5.