Low-pressure phase diagram of crystalline benzene from quantum Monte Carlo

Low-pressure phase diagram of crystalline benzene from quantum Monte Carlo

Sam Azadi s.azadi@ucl.ac.uk Department of Physics and Astronomy, University College London, Thomas Young Center, London Centre for Nanotechnology, London WC1E 6BT, United Kingdom    R. E. Cohen Extreme Materials Initiative, Geophysical Laboratory, Carnegie Institution for Science, Washington DC, 20015, USA; Department of Earth- and Environmental Sciences, Ludwig Maximilians Universität, Munich 80333, Germany; and Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
July 4, 2019

We study the low-pressure (0 to 10 GPa) phase diagram of crystalline benzene using quantum Monte Carlo (QMC) and density functional theory (DFT) methods. We consider the , , and structures as the best candidates for phase I and phase II. We perform diffusion quantum Monte Carlo (DMC) calculations to obtain accurate static phase diagrams as benchmarks for modern van der Waals density functionals. We use density functional perturbation theory to compute phonon contribution in the free-energy calculations. Our DFT enthalpy-pressure phase diagram indicates that the and structures are the most stable phases within the studied pressure range. The DMC Gibbs free-energy calculations predict that the room temperature to phase transition occurs at 2.1(1) GPa. This prediction is consistent with available experimental results at room temperature. Our DMC calculations show an estimate of 50.60.5 kJ/mol for crystalline benzene lattice energy.

I Introduction

Molecular crystals, including organic and inorganic, are vital in understanding the physics and chemistry of the Earth and planets. They are also of considerable technological interest. Low-Z molecular systems are among the most abundant in the solar system, as represented by planetary gases and ices. Their behaviour at high pressures is crucial in modelling the structure, dynamic, and evolution of the large planets. Moreover, compression of molecular systems provides the opportunities to form new materials, possibly with novel properties, such as high-temperature superconductivity and disordered and amorphous materials. One of the simplest organic molecular solids is crystalline benzene with aromatic van der Waals (vdW) interactions. Given its simplicity, high symmetric, and rigid molecular structure, crystalline benzene has become the model structure for calculating the lattice model vibrations in molecular crystals. Benzene has been extensively studied theoretically and experimentally Thièry and Lèger (1988); Wen et al. (2011); Piermarini et al. (1969); Raiteri et al. (2005). However, the phase transitions and intermolecular interactions are still controversial. The main goal of this paper is to present a comprehensive study of the phase transition of crystalline benzene at low pressures.

Early experiments by BridgmanBridgman (1941) revealed that liquid benzene crystallises at 68 MPa with space group symmetry and closest intermolecular distance of 3.5 . This structure, also, was confirmed at zero pressure and 270 KCox and Smith (1954); Cox et al. (1958). This phase I is also stable at lower temperatures of 218 and 138 KBacon et al. (1964). Since then, two experimental phase diagrams have been proposed for crystalline benzene. First, based on the phase diagram suggested by Thiéry and LégerThièry and Lèger (1988), liquid benzene crystallises at room temperature and pressure 700 bar within an orthorhombic structure , which is labeled as phase I. Phase II was suggested to exist between 1.4 and 4 GPa. Phases I and II primitive unit cells contain four benzene molecules (Z=4). Phase III is stable between 4 and 11 GPa. The symmetry of phase III is with two benzene molecules per monoclinic primitive unit cell (Z=2). Second, the phase diagram developed by Ciabini et al.Ciabini et al. (2005, 2007) which based on it phase I is orthorhombic Z=4 and phase II is monoclinic Z=2 Piermarini et al. (1969). Their results are obtained by means of infrared spectroscopy and X-ray analysis under high pressure. The phase is stable up to pressures 2025 GPa. This phase diagram only consists of two phases (I and II), and this same result has been reported by other experiments Katrusiak et al. (2010). Katrusiak et alKatrusiak et al. (2010), have determined the crystal structures of phases I and II at 295 K. The results of their study confirm Ciabini et al.’s phase diagram and show that the structures of phases I and II are Z=4 and Z=2, respectively. The results also indicate the absence of other benzene phases in the pressure range up to 5 GPa.

The crystalline benzene phase diagram is a challenge for first-principles theory because the energy differences are insignificant, and they are governed by vdW interactions. The energy difference between crystalline benzene and its low-energy polymorphs under pressure is less than few kJ/mol. Metadynamics calculations predict seven phasesRaiteri et al. (2005) as phases I ( Z=4), I ( Z=4), II ( Z=4), III ( Z=2), III ( Z=4), IV ( Z=4), and V ( Z=2). In their calculations, they have used numerous randomly generated metastable crystal structures as starting points for the metadynamics. A few metadynamics steps are often sufficient to obtain a more stable structure, which most of the time is similar to Z=4 or Z=2. Density functional theory (DFT) has also been used to compute the lattice energy of crystalline benzenePodeszwa et al. (2008); Grimme et al. (2010). Wen et al. employed DFT formalism and used Perdew-Burke-Ernzerhof (PBE)Perdew et al. (1996) exchange correlation functionals to study the phase diagram of crystalline benzene up to 300 GPaWen et al. (2011). They explained the complexities observed in benzene at high pressure. In the moderate pressure regime ( P 20 GPa), they found that the structure is stable up to 4 GPa, the phase is preferred in the pressure range of 47 GPa, and the structure shows the lowest enthalpy at higher pressures. Therefore, they labelled the , , and structures as phases I, II, and III, respectively. The present study shows that the structure is unstable in the pressure range of 010 GPa. Thus, the and structures are labelled as phases I and II, respectively.

Recently, quantum chemistry methods have been applied to benzene to obtain sub-kilojoule/mole accuracy in the lattice energy for crystalline benzene Yang et al. (2014). Tremendous measures are necessary to obtain such accuracy. In this work, we will show that QMC is an alternative efficient approach to achieve or surpass such accuracy in benzene crystals, as we previously demonstrated for the benzene dimerAzadi and Cohen (2015).

Quantum Monte Carlo (QMC), which approximately solves the electronic Schrödinger equation stochasticallyFoulkes et al. (2001), can yield highly accurate energies for atomsMarchi et al. (2009); Brown et al. (2007), moleculesGurtubay and Needs (2007); Trail and Needs (2008); Azadi et al. (2015), and crystalsMarchi et al. (2011); Mostaani et al. (2015); Azadi and Foulkes (2015). Previous studies have shown that diffusion quantum Monte Carlo (DMC) can provide accurate energies for vdW systemsCox et al. (2014); Ma et al. (2011); Al-Hamdani et al. (2014); Benedek et al. (2006). DMC can also produce an accurate description of the phase diagram of materials under pressureDrummond et al. (2015); Azadi et al. (2014, 2013). In general, QMC-based methods are faster than the most accurate post-Hartree-Fock schemes for large number of particles N. The computational cost of QMC methods scales usually as - depending on the method.

We have demonstrated that QMC can provide chemical accuracy for the benzene dimer systemAzadi and Cohen (2015). We have found optimal variational quantum Monte Carlo (VMC) and DMC binding energies of 2.3(4) and 2.7(3) kcal/mol. The best estimate of the CCSD(T)/CBS limit is 2.65(2) kcal/molMiliordos et al. (2014). The consistency among our results, experiments, and quantum chemistry methods, is an important sign of the capability of the QMC-based methods to provide an accurate description of weak intermolecular interactions based on vdW dispersive forces.

In this study, we examine the Z=4 to Z=2 phase transition of crystalline benzene at low pressures. We consider the and structures as best candidates for Z=4 and the structure for Z=2. We study pressures below 10 GPa. We obtain static and dynamic phase diagrams where the phonon contribution to the free energy is included. We employ different vdW functionalsKlimeš and Michaelides (2012) and compare them with conventional DFT functionals. We perform QMC calculations to obtain the static enthalpy-pressure phase diagram of crystalline benzene. We will show that DMC provides accurate results for the phase diagram of crystalline benzene.

Ii Computational Details

Given that the energy differences between crystalline benzene structures are small, the calculations must be performed with the highest possible numerical precision. Our DFT calculations were carried out within the pseudopotential and plane-wave approach using the Quantum ESPRESSO suite of programset al. (2009). All DFT calculations used ultrasoft pseudopotentialsRappe et al. (1990). Pseudopotentials were obtained by PBEPerdew et al. (1996) exchange correlation functionals. We used a basis set of plane waves with an energy cutoff 100 Ry. Geometry and cell optimisations employed a dense -point mesh. The quasi-Newton algorithm was used for cell and geometry optimisation, with convergence thresholds on the total energy and forces of 0.01 mRy and 0.1 mRy/Bohr, respectively, to guarantee convergence of the total energy to less than 1 meV/proton and the pressure to less than 0.1 GPa/proton.

To include the effects of zero point energy (ZPE), vibrational frequencies were calculated using density-functional perturbation theory as implemented in Quantum ESPRESSOet al. (2009). The ZPE per proton at a specific cell volume was estimated within the quasi-harmonic approximation: , where . and are the numbers of vibrational modes in the simulation cell and phonon wave vectors , respectively, and the summation over includes all -points on a grid in the Brillouin zone.

The thermodynamic properties are determined by the Helmholtz free energy . The free energy can be written as the sum of an electronic and a vibrational term. The electronic entropy is negligible for insulators: . In our calculations, the electronic part is obtained using the DMC method. Thus, the main quantity to calculate for obtaining the thermal properties and finite temperature phase diagram is the vibrational free energy . We use quasi-harmonic approximation to calculate the vibrational free energyBaroni et al. (2001):


where , , and are Boltzmann constant, unit cell volume, and eigenvalue of the phonon Hamiltonian, respectively. The pressures are calculated from the Helmholtz free energies by

We used the casino codeNeeds et al. (2010) to perform fixed-node DMC simulations with a trial wave function of the Slater-Jastrow (SJ) form:


where is a -dimensional vector of the positions of the electrons, is the position of the ’th spin-up electron, is the position of the ’th spin-down electron, is a Jastrow factor, and and are Slater determinants of spin-up and spin-down one-electron orbitals. These orbitals were obtained from DFT calculations performed with the plane-wave-based Quantum ESPRESSO codeet al. (2009), employing Trail-NeedsTrail and Needs (2005a, b) Hartree-Fock pseudopotentials. For the QMC study of C and CH-based systems, the Hartree-Fock description of the core is more accurateGreeff and Jr. (1998). A detailed study of silicon also showedZuo et al. (1997) that Hartree-Fock provides the most accurate description of the core density compared with generalised gradient approximation and local density approximation (LDA).

We selected a very large basis-set energy cutoff of 200 Ry to approach the complete basis-set limitAzadi et al. (2010). The plane-wave orbitals were transformed into a localised “blip” polynomial basisAlfè and Gillan (2004). Our Jastrow factor consists of polynomial one-body electron-nucleus, two-body electron-electron, and three-body electron-electron-nucleus terms, the parameters of which were optimised by minimising the variance of the local energy at the VMC levelUmrigar et al. (1988); Drummond and Needs (2005). Our DMC calculations were performed at two different time steps 0.01 and 0.02 a.u. The target population control is two times larger for time step 0.02 a.u. We extrapolated our DMC energies to zero time step using a linear fitting. The time step error is linear in the time step. The population control error also is linear as function of reciprocal of the target population. Therefore, it is possible to remove both time step and population control errors simultaneously by linearly extrapolation to zero-time step.

Iii Results and discussion

iii.1 Geometry Analysis

In this section we discuss the results of our geometry optimization. We study the evolution of benzene molecule distances by increasing the pressure. The structure optimization results are compared with experiments.

The primitive unit cells of the , , and structures of solid benzene contain four, four, and two benzene molecules, respectively, as shown in figure 1. The and structures have orthorhombic and tetragonal primitive unit cells, respectively, whereas the primitive unit cell is monoclinic.

Figure 1: (colour online) Primitive unit cells of the (left), (middle), and (right) structures of solid benzene at low-pressure range. The and primitive unit cells contain four benzene molecules (Z=4), whereas the structure has two benzene molecules (Z=2) in a monoclinic primitive unit cell.

For geometry analysis of Z=4 and Z=2 structures, we focus on the and . We will show in the next section that these two structures are the best candidates for the phases I (Z=4) and II (Z=2). Our structure optimization indicates that the molecular orientations do not change significantly within the studied pressure range. We calculated the distances between C atoms on nearest-neighbour (nn) benzene molecules. The nn CC distances between molecules as function of pressure are reported in figure 2. The nn CC distances for and structures are calculated using vdWDion et al. (2004); Lee et al. (2010) and conventional DFT functionals.

Figure 2: (colour online) Pressure evolution of nearest-neighbour (nn) CC distances for and . The results are obtained by vdW and conventional functionals.

The vdW functionals, particularly vdW-DF2Lee et al. (2010), nn CC distances are in good agreement with experimentCiabini et al. (2005, 2007). The differences between vdW-DF1Dion et al. (2004) and vdW-DF2Lee et al. (2010) nn CC distances reduce with increasing pressure. The PBE nn CC distances are close to vdW functional results at lower pressures, whereas the differences between PBE and vdW results increases with increasing pressure. The PBE nn CC distances at higher pressures are close to LDA results. The BLYP nn CC distances are the largest at low pressures. However, BLYP nn CC distances are more similar to vdW results at pressures larger than 5 GPa.

Figure 3: (colour online) Pressure evolution of the difference between nearest-neighbour (nn) CC distances and the sum of van der Waals radii of C atoms (). The results are obtained for and structures using vdW and conventional functionals.

The van der Waals radius of Carbon atom is 1.7 . In crystalline benzene the benzene molecules are held together by van der Waals forces. The nearest that two atoms belonging to different benzene molecules can approach each other can be estimated by the sum of . We calculated the difference between nearest-neighbour (nn) CC distances and the sum of (). Figure  3 illustrates for and structures. The results are obtained by vdW and conventional functionals. At the same pressure, all the functionals give larger for structure. Our EOS calculations, which are presented in figure 7, indicate that at the same pressure molecular density of is larger than . LDA and BLYP provide the smallest and largest . Consequently they yield the smallest and largest vdW radii for C atom. Unlike the other functionals, the BLYP decline rapidly with increasing the pressure. According to LDA, results, benzene molecules are strongly bonded at pressures larger than 0.2 GPa. In lower pressures PBE is close to obtained by vdW functionals. With increasing the pressure PBE results become closer to LDA. According to vdW-DF1 results, the benzene molecules in structure is bonded above 3 GPa. Whereas vdW-DF2 results indicate that bonding between benzene molecules in phase could happen around 2 GPa. Based on the experimental phase digramCiabini et al. (2005, 2007), the phase is stable at pressures below 1.4 GPa. Our vdW results show that there are no strong bonds between benzene molecules in phase. In the structure the benzene molecules only interact through weak dispersive forces.

iii.2 Ground State DFT Phase Diagram

We begin our phase diagram study by DFT enthalpy-pressure calculations at zero temperature. We first present our static phase diagram results where the Born–Oppenheimer (BO) approximation is used. According to BO approximation the electronic and nuclear wave functions can be separated. It is also assumed that the nuclei are infinitely massive and the total nuclear momentum contribution in the Hamiltonian is zero. To find out the best candidate for Z=4 at the studied pressures, we used the PBEPerdew et al. (2008) and vdW-DF2Lee et al. (2010) functionals to calculate the enthalpy difference between the , , and structures. We performed calculations at six different volumes corresponding to DFT pressures of 0, 2, 4, 6, 8, and 10 GPa (Figure 7). Based on the linear fitting of the PBE results on two enthalpy-pressure points at P = 0 and 10 GPa, the structure is stable up to 3.6 GPa, whereas is stable in the pressure range of 3.66.8 GPa, and finally the structure has lowest enthalpy in pressures higher than 6.8 GPa.

A line between these two enthalpy-pressure points gives excellent agreement with the previous PBE computations by Wen et al.Wen et al. (2011) (Fig. 2(a)). However we find this result to be inaccurate, and a denser set of points in this pressure range is needed.

Figure 4: (colour online) Enthalpy difference between the , , and structures as function of applied pressure. (a) The results are calculated using PBE and linear fitting on two enthalpy-pressure points at P = 0, and 10 GPa. (b) The phase diagram is simulated using DFT-PBE and vdW density functional of vdW-DF2. We used the Vinet EOS and six enthalpy-pressure points at P = 0, 2, 4, 6, 8, and 10 GPa.

Using the VinetVinet et al. (1986) equation of state (EOS) we found that the structure is not stable in the pressure range of 010 GPa. The results of our EOS calculations are presented in figure 7. The enthalpy difference between the , , and structures versus pressure is calculated using PBE and vdW functionals (Figure 4(b)). We find that instability of is independent of employed functional. Our results indicate that and are the most stable structures in the studied pressure ranges. These results are consistent with the experimental phase diagram proposed by Ciabini et al.Ciabini et al. (2005, 2007); Piermarini et al. (1969); Katrusiak et al. (2010). Therefore, in the rest of this paper, we label and as phases I and II, respectively.

To study the importance of dispersion effects, we calculated the phase diagram of crystalline benzene using different functionals (Figure 5). We employed vdW-DF1Dion et al. (2004), vdW-DF2Lee et al. (2010), vdW-DF-obk8, vdW-DF-ob86, vdW-DF2-B86RKlimeš et al. (2012); Klimeš and Michaelides (2012), vdW-DF-C09, vdW-DF2-C09Cooper (2010), vdW-DF-cxBerland and Hyldgaard (2014), and vdW-rVVSabatini et al. (2013); Vydrova and Voorhis (2010) vdW functionals. Except rVV functional, the nonlocal term in the other vdW functionals is either vdW-DF1Dion et al. (2004) or vdW-DF2Lee et al. (2010). Employing various gradient corrections to the exchange energy results in a variety of vdW functionals. We also determined the phase diagram using conventional DFT functionals, including PBEPerdew et al. (2008), LDAPerdew and Zunger (1981), and BLYPLee et al. (1988).

Figure 5: (colour online) Enthalpy difference between (phase I) and (phase II) as function of pressure obtained with vdW and conventional DFT funcationals. The left panel shows the results of vdW-DF1Dion et al. (2004), vdW-DF2Lee et al. (2010), rVVSabatini et al. (2013); Vydrova and Voorhis (2010), obk8, B86R, ob86Klimeš et al. (2012); Klimeš and Michaelides (2012), DFCxBerland and Hyldgaard (2014), DFC09, and DF2C09Cooper (2010) vdW functionals. The right panel illustrates the results of conventional DFT including PBEPerdew et al. (2008), LDAPerdew and Zunger (1981), and BLYPLee et al. (1988).

The vdW functionals yield different III phase transition pressure. Figure 6 illustrates to phase transition pressures which are obtained by different vdW functionals. The CxBerland and Hyldgaard (2014) and DF2C09Cooper (2010) functionals show the lowest and highest phase transition pressures, respectively. The difference between largest and smallest phase transition pressure is about 1.1 GPa. This value corresponds to inaccuracy in prediction of phase transition pressure by vdW functionals. It should be noted that the experimental to phase transition occurs within 1.4 GPa pressure window. The results of the PBE, LDA, and BLYP functionals predict that the phase III transition occurs at 5.2, 5.6, and 3.5 GPa, respectively.

Figure 6: (colour online) to phase transition pressure. The results are calculated using different vdW functionals as explained in text.

Our phase diagram calculations indicate that vdW results are in better agreement with experiments than the conventional functionals. Between the PBE and BLYP functionals, the PBE results are closer to vdW-DF1 and vdW-DF2 at low pressures. The difference between PBE phase transition pressure and vdW-DF1 and vdW-DF2 phase transition pressures are 4.55 and 3.9 GPa, respectively. However the difference between BLYP phase transition pressure and vdW-DF1 and vdW-DF2 phase transition pressures are 2.85 and 2.2 GPa, respectively. As we discussed in the previous section, vdW-DF1 and vdW-DF2 are positive for both and structures below 2 GPa where the phase transition between them happens. Therefore the phase III transition in crystalline benzene occurs without any intermolecular contacts. This transition occurs only due to dispersion effects.

Figure 7: (colour online) EOS of and structures obtained by vdW ((a) and (b)) and conventional ((c) and (d)) DFT functionals. The results are compared with experimental data which are reported in References Ciabini et al., 2005; Katrusiak et al., 2010.

Using our DFT results we compute the EOS of and structures. Figure 7 represents the results which are obtained by vdW and conventional functionals. We compare our DFT results with experiments which are reported in References Ciabini et al., 2005; Katrusiak et al., 2010. The experimental results in Ref. Ciabini et al., 2005 are data for crystalline benzene at 540 K that have been fitted by the Vinet EOS. The second experimental results Katrusiak et al. (2010) belong to crystalline benzene at lower pressures and 295 K. Among DFT conventional functionals used in this study only the PBE results are close to experiments. The BLYP and LDA curves lie far above and below experimental curves, respectively. In general, the vdW results are in good agreement with experiments. At lower pressures vdW-DF1Dion et al. (2004), vdW-DF2Lee et al. (2010), and DFcxBerland and Hyldgaard (2014) points for phase are close to experiments. With increasing the pressure, the curves computed with vdW-DF2Lee et al. (2010), obk8Klimeš et al. (2012); Klimeš and Michaelides (2012), rVVSabatini et al. (2013); Vydrova and Voorhis (2010), and B86RKlimeš et al. (2012); Klimeš and Michaelides (2012) are close to experimental points. The rVV functional has a different nonlocal correlation kernel, whereas other vdW functionals are the modified versions of vdW-DF1 or vdW-DF2. Our EOS calculations indicate that the modifications bring the vdW-DF1 and vdW-DF2 curves below experimental ones. It is hard to conclude whether these modifications improve the accuracy of vdW-DF1 and vdW-DF2 fucntionals, especially in the case of vdW-DF2 fucntional, which overall gives the most accurate results. Our ground state EOS calculations indicate that at fixed pressure the volume per benzene molecule for phase is larger than . This is in agreement with finite temperature experimental measurements. This conclusion is also independent of used DFT functionals.

Figure 8: (colour online) ZPE of and structures obtained by DFT. Geometries are accurately optimised by two functionals: (a) vdW-DF2Lee et al. (2010) and (b) vdW-DF1Dion et al. (2004).

To investigate the ZPE contribution in phase diagram calculations, we simulated the difference between the gas and crystal ZPEs. The ZPE of the and structures with respect to gas phase is shown as function of pressure (Figure 8). We used the vdW-DF2 and vdW-DF1 functionals to optimise the structures for phonon calculations. ZPE is obtained using quasi-harmonic approximation, as explained in the previous section. Within the studied pressure range, the difference between the ZPE of phases I and II is less than 2 meV/atom. The vdW-DF2 results indicate that the phase III ZPE transition happens at 0.6 GPa, whereas the vdW-DF1 results predict that the phase III ZPE transition occurs at 1.65 GPa (Figure  8). The difference between the ZPE of phases I and II increases with pressure.

Figure 9: (colour online) Static and dynamic phase transition of to obtained by (a) vdW-DF2Lee et al. (2010) and (b) vdW-DF1Dion et al. (2004) vdW functionals.

The ZPE correction to the cohesive energy of crystalline benzene was previously calculatedBludský and Rubes̆ (2008). They evaluated the ZPE using -point harmonic frequencies at the PBE level. They found that the ZPE of the structure is 44 meV/molecule. In their calculations, they employed experimentally reportedDavid et al. (1992) orthorhombic cell without full three-dimensional optimisation. Finite-temperature experimentsNakamura and Miyazawa (1969) show that the ZPE of crystalline benzene is 2.8 kJ/mol (29.02 meV/molecule). The ZPE experimental result is also employed to investigate the binding energy of benzene crystalLu et al. (2009). An estimate of 4.8 kJ/mol was obtained using DFT many-body dispersion methodReilly and Tkatchenko (2013). This ZPE is significantly larger than an estimate of 2.8 kJ/mol which is obtained by finite molecular cluster calculationsPodeszwa et al. (2008); Ringer and Sherrill (2008). Our ZPE results are close to PBC-DFT calculationsde-la Roza and Johnson (2012), where an estimate of 2.6 kJ/mol is obtained using the PBE functional.

The static phase diagrams in Figure 5 assume that the atoms are infinitely massive. We calculate the dynamic phase diagram by adding the ZPE to the static results. Figure 9 illustrates the dynamic phase diagrams of crystalline benzene at the DFT level. The vdW-DF2 results indicate that adding ZPE lowers the phase transition by 0.02 GPa, and the to phase transition pressure is 1.42 GPa. The vdW-DF1 results predict that the phase transition occurs at 0.66 GPa, which is 0.03 GPa higher than the static phase transition pressure. The results of comparing the static and dynamic phase diagrams indicates that the ZPE contribution is negligible.

iii.3 Finite Temperature DMC Phase Diagram

In this section we present our finite temperature phase diagram calculations. We use QMC based methods to calculate the electronic structure ground state energy. The inadequacy of mean-field-like DFT calculations of hydrogen-rich systems was demonstrated beforeAzadi and Foulkes (2013). To obtain reliable results, going beyond DFT-based methods and properly considering many-body effects are necessary. The DMC is generally considered as the most accurate first-principle method available in studying the phase diagram of hydrogen-dominant materialsDrummond et al. (2015); Azadi et al. (2014). In addition DMC is an effective method to study non-covalent systems. It can reach and go beyond the chemical accuracy which is desired for non-covalent systemsAzadi and Cohen (2015).

We perform DMC calculations to obtain the wave-function-based phase diagram for crystalline benzene at low-pressures. We use the vdW-DF2 optimized structure for our DMC calculations. As we demonstrated in our DFT calculations, vdW-DF2 functional gives the closest results to experiment. The DMC results for energies in the limit of infinite system size are obtained by extrapolation using DMC energy data at and simulation cells. Extrapolation is advantageous because it can approximately account for finite-size effects that are not considered in the other correction schemes, such as finite-size effects in the fixed-node error. In addition, it does not suffer from the reliance on stochastically optimised trial wave functions that affects the kinetic-energy correction, because it is purely based on SJ DMC energiesAzadi and Foulkes (2015); Drummond et al. (2008).

Vol E() E() E()
781.086 -1024.7464(4) -1022.7063(8) -1022.4140(8)
693.335 -1024.6376(5) -1022.5976(7) -1022.3052(7)
646.304 -1024.4880(5) -1022.4478(5) -1022.1556(5)
619.413 -1024.3588(4) -1022.3185(5) -1022.0264(5)
Table 1: DMC energies of the (phase I) structure. Energies are obtained at two simulation cells containing = 48 and = 384 atoms. Linear extrapolated energies are shown as E(). Energy (E) and volume (Vol) are in eV and Bohr per benzene molecule, respectively.

Table  1 lists the DMC energies of the structure at four primitive unit-cell volumes. We consider two simulation cells for each density containing 48 and 384 atoms. DMC energy at thermodynamic limit is obtained by linear extrapolation in .

Vol E() E() E()
760.8398 -1024.8824(5) -1022.6683(8) -1022.3514(8)
670.2382 -1024.7736(6) -1022.5582(8) -1022.2426(8)
622.5712 -1024.6240(5) -1022.4099(7) -1022.0929(8)
589.6330 -1024.4948(5) -1022.2793(8) -1021.9637(8)
Table 2: DMC energies of the (phase II) structure. Energies are obtained at two simulation cells containing = 24 and = 192 atoms. Linear extrapolated energies are shown as E(). Energy (E) and volume (Vol) are in eV and Bohr per benzene molecule, respectively.

Table  2 shows the DMC energies of the structure at different primitive unit-cell volumes. We consider two simulation cells for each density containing 24 and 192 atoms. DMC energy at infinite system size limit is calculated by linear extrapolation in .

To identify enthalpy-pressure curves for the and structures, we fitted model equations of state to our finite-size-corrected DMC energy against volume . We used the Vinet EOSVinet et al. (1986) to fit our total energies and propagate errors using classical statistics. The pressure and the enthalpy is , where is DMC electronic structure energy of system.

Figure 10: (colour online) (left) DMC energy of the and structures as function of volume per benzene molecule. Energy error bars are included in point sizes and are of the order of meV. (right) Relative enthalpies of the and structures as function of pressure. The widths of the DMC lines indicate the estimated uncertainties in the enthalpies because of statistical and systematic errors.

Figure 10 (left) illustrates the DMC energy of phases I and II of crystalline benzene as function of volume per benzene molecule. With increasing density, phase II becomes favourable over phase I in the structure. Figure 10 (right) shows the relative enthalpies of the and structures. Based on our static enthalpy-pressure phase diagram, the to phase transition occurs at pressure 2.60.1 GPa. The use of the DMC method has significant consequences for the static-lattice relative enthalpies of the studied structures. Compared with vdW-DF2, the DMC enthalpy-pressure results predict that the phase III transition occurs at 1.2 GPa higher pressure. Among conventional DFT functionals, the BLYP results are closest to DMC. The difference between DMC and BLYP phase transition pressure is 0.9 GPa.

Figure 11: Phonon contribution to the Helmholtz free energies of the and structures of crystalline benzene. The geometries are optimized using vdW-DF2Lee et al. (2010) functionals.

To obtain the phase diagram at finite temperature, we used quasi-harmonic approximation to obtain lattice dynamic contribution to the free energies. Phonons have contributed to the Helmholtz free energies of crystalline benzene (Figure 11). We used vdW-DF2 functionals to optimise the and structures at different pressures. Vibrational free energies are calculated at different temperatures of 50, 100, 200, and 300 K. At room temperature and 0 GPa vibrational free energy of is higher than . Meanwhile, the vibrational free energies of become higher than by increasing the pressure. This results indicates the stability of the phase at ambient conditions, which is also observed experimentallyCiabini et al. (2005, 2007).

We calculated relative Gibbs free energies of the and structures at different temperatures (Figure 12). The static electronic structure results are obtained by DMC calculations. Our results predict that the room temperature to structure transformation happens at 2.1(1) GPa. Experiments indicate the transition to phase II occurs at room temperature and around 1.4 GPaPiermarini et al. (1969). The I to II phase transition was found to be extremely sluggish, and it can be speeded up by heating the sampleThièry and Lèger (1988). Keeping the low-pressure phase I, , in a metastable state at least up to 3 GPa is possible without heatingBridgman (1941). Experimentally achieving low-enough temperature results is extremely difficult. Our DMC phase diagram at low temperature predicts that the phase could be stable up to 2.6(1) GPa. The phase diagram that we obtained by combining DMC static-lattice energies and quasi-harmonic vibrational energies can be extended to higher pressures.

Figure 12: Relative Gibbs free energies of the and structures. (a) 0 K, and (d) 300 K. The Gibbs free energies are calculated using static-lattice DMC calculations together with DFT quasi-harmonic vibrational calculations.

As the final step of our study, we calculated the lattice energy of crystalline benzene at ambient conditions. The cohesive energy yields the strength of the vdW forces holding the crystalline benzene together. We used our DMC and ZP energies for structure. The cohesive energy is calculated using the difference between total energies of structure and its fragments. Cohesive energy calculation is a precise test of DMC method, since it has to accurately describe two different systems of benzene molecule and crystalline benzene. The electronic structure of these two systems are not similar. In our DMC lattice energy calculation, we used same time step of 0.01 a.u for both crystal and molecule. We found an estimate of 50.60.5 kJ/mol for lattice energy. Ab initio many-electron wave functions methods provide an estimate of 55.900.76 kJ/mol for benzene crystal lattice energy at zero temperatureYang et al. (2014). The experimental lattice energy at same condition is 55.32.2 kJ/molYang et al. (2014). We used conventional Jastrow factor in our DMC calculations. In principle, the DMC lattice energy can be systematically improved by accurately taking into account the correlation energy and also decreasing the fixed-node errors. These purposes can be fulfilled by adding additional terms in Jastrow factor and using backflow transformationsAzadi and Cohen (2015). However, improving the DMC lattice energy until it converges to exact results requires huge amount of computational time.

Iv Conclusion

We have comprehensively studied the crystalline benzene phase diagram at pressures below 10 GPa. We have used different vdW functionals and also three most used conventional functionals to obtain DFT energy of system. The vdW-DF2 results of our study indicated that the and structures are the best candidates for phases I and II, respectively. We have used the accurate DMC method to calculate the ground-state electronic structure energy of system. We have compared static enthalpy-pressure phase diagrams which are obtained by DFT and DMC methods. We used quasi-harmonic approximation and density functional perturbation theory to calculate the phonon contribution to the free energy of system. Our Gibbs free energy phase diagram predicts that at room temperature, the phase III transition occurs at 2.1(1) GPa, which is in good agreement with experiments. We have found DMC lattice energy of 50.60.5 kJ/mol for crystalline benzene at ambient conditions. The results of our study indicate the importance of many-body electronic structure calculation to obtain a reliable phase diagram for molecular crystals.

This study utilised computing facilities provided by ARCHER, the UK national super computing service, and by the University College London high-performance computing centre. S. Azadi acknowledges that the results of this research have been obtained using the PRACE-3IP project (FP7 RI-312763) resource ARCHER based in the UK. The authors acknowledge the financial support of the European Research Council under the Advanced Grant ToMCaT (Theory of Mantle, Core, and Technological Materials). R. E. Cohen acknowledges the support of the Carnegie Institution for Science


  • Thièry and Lèger (1988) M. M. Thièry and J. M. Lèger, J. Chem. Phys. 89, 4255 (1988).
  • Wen et al. (2011) X.-D. Wen, R. Hoffmann,  and N. W. Ashcroft, J. Am. Chem. Soc. 133, 9023 (2011).
  • Piermarini et al. (1969) G. J. Piermarini, A. D. Mighell, C. Weir,  and S. Block, Science 165, 3461 (1969).
  • Raiteri et al. (2005) P. Raiteri, R. Martonak,  and M. Parrinello, Angew. Chem. Int. Ed. 44, 3769 (2005).
  • Bridgman (1941) P. W. Bridgman, J. Chem. Phys 9, 794 (1941).
  • Cox and Smith (1954) E. G. Cox and J. A. S. Smith, Nature 173, 75 (1954).
  • Cox et al. (1958) E. G. Cox, D. W. Cruickshank,  and J. A. S. Smith, Proc. R. Soc. London A 247, 1 (1958).
  • Bacon et al. (1964) G. E. Bacon, N. A. Curry,  and S. A. Wilson, Proc. R. Soc. A 279, 98 (1964).
  • Ciabini et al. (2005) L. Ciabini, F. A. Gorelli, M. Santoro, R. Bini, V. Schettino,  and M. Mezouar, Phys. Rev. B 72, 094108 (2005).
  • Ciabini et al. (2007) L. Ciabini, M. Santoro, F. A. Gorelli, R. Bini, V. Schettino,  and S. Raugei, Nature Mater. 6, 39 (2007).
  • Katrusiak et al. (2010) A. Katrusiak, M. Podsiadzo,  and A. Budzianowski, Crystal Growth & Design 10, 3461 (2010).
  • Podeszwa et al. (2008) R. Podeszwa, B. M. Rice,  and K. Szalewicz, Phys. Rev. Lett. 101, 115503 (2008).
  • Grimme et al. (2010) S. Grimme, J. Antony, S. Ehrlich,  and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  • Perdew et al. (1996) J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  • Yang et al. (2014) J. Yang, W. Hu, D. Usvyat, D. Matthews, M. Schutz,  and G. K.-L. Chan, Science 345, 640 (2014).
  • Azadi and Cohen (2015) S. Azadi and R. E. Cohen, J. Chem. Phys. 143, 104301 (2015).
  • Foulkes et al. (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs,  and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  • Marchi et al. (2009) M. Marchi, S. Azadi, M. Casula,  and S. Sorella, J. Chem. Phys. 131, 154116 (2009).
  • Brown et al. (2007) M. D. Brown, J. R. Trail, P. L. Ríos,  and R. J. Needs, J. Chem. Phys. 126, 224110 (2007).
  • Gurtubay and Needs (2007) I. G. Gurtubay and R. J. Needs, J. Chem. Phys. 127, 124306 (2007).
  • Trail and Needs (2008) J. R. Trail and R. J. Needs, J. Chem. Phys. 128, 204103 (2008).
  • Azadi et al. (2015) S. Azadi, R. Singh,  and T. D. Kuhne, Int. J. Quantum Chem. 115, DOI: 10.1002/qua.25005 (2015).
  • Marchi et al. (2011) M. Marchi, S. Azadi,  and S. Sorella, Phys. Rev. Lett. 107, 086807 (2011).
  • Mostaani et al. (2015) E. Mostaani, N. Drummond,  and V. Falḱo, Phys. Rev. Lett. 115, 115501 (2015).
  • Azadi and Foulkes (2015) S. Azadi and W. M. C. Foulkes, J. Chem. Phys. 143, 102807 (2015).
  • Cox et al. (2014) S. J. Cox, M. D. Towler, D. Alfè,  and A. Michaelides, J. Chem. Phys. 140, 174703 (2014).
  • Ma et al. (2011) J. Ma, A. Michaelides,  and D. Alfè, J. Chem. Phys. 134, 134701 (2011).
  • Al-Hamdani et al. (2014) Y. S. Al-Hamdani, D. Alfè,  and O. A. von Lilienfeld, J. Chem. Phys. 141, 18C530 (2014).
  • Benedek et al. (2006) N. A. Benedek, I. K. Snook, M. D. Towler,  and R. J. Needs, J. Chem. Phys. 125, 104302 (2006).
  • Drummond et al. (2015) N. Drummond, B. Monserrat, J. Lloyd-Williams, P. L. Rìos,  and R. J. Needs, Nature Communications 6, 7749 (2015).
  • Azadi et al. (2014) S. Azadi, B. Monserrat, W. M. C. Foulkes,  and R. J. Needs, Phys. Rev. Lett. 112, 165501 (2014).
  • Azadi et al. (2013) S. Azadi, W. M. C. Foulkes,  and T. D. Kuhne, New Journal of Physics 15, 113005 (2013).
  • Miliordos et al. (2014) E. Miliordos, E. Aprà,  and S. S. Xantheas, J. Phys. Chem. A 118, 7568 (2014).
  • Klimeš and Michaelides (2012) J. Klimeš and A. Michaelides, J. Chem. Phys. 137, 120901 (2012).
  • et al. (2009) P. G. et al., J. Phys.: Condens. Matter 21, 395502 (2009).
  • Rappe et al. (1990) A. M. Rappe, K. M. Rabe, E. Kaxiras,  and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
  • Baroni et al. (2001) S. Baroni, S. de Gironcoli,  and A. D. Corso, Rev. Mod. Phys. 73, 515 (2001).
  • Needs et al. (2010) R. J. Needs, M. D. Towler, N. D. Drummond,  and P. L. Ríos, J. Phys.: Condens. Matter 22, 023201 (2010).
  • Trail and Needs (2005a) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 174109 (2005a).
  • Trail and Needs (2005b) J. R. Trail and R. J. Needs, J. Chem. Phys. 122, 014112 (2005b).
  • Greeff and Jr. (1998) C. W. Greeff and W. A. L. Jr., J. Chem. Phys. 109, 1607 (1998).
  • Zuo et al. (1997) J. M. Zuo, P. Blaha,  and K. Schwarz, J. Phys.: Condens. Matter 9, 7541 (1997).
  • Azadi et al. (2010) S. Azadi, C. Cavazzoni,  and S. Sorella, Phys. Rev. B 82, 125112 (2010).
  • Alfè and Gillan (2004) D. Alfè and M. J. Gillan, Phys. Rev. B 70, 161101 (2004).
  • Umrigar et al. (1988) C. J. Umrigar, K. G. Wilson,  and J. W. Wilkins, Phys. Rev. Lett. 60, 1719 (1988).
  • Drummond and Needs (2005) N. D. Drummond and R. J. Needs, Phys. Rev. B 72, 085124 (2005).
  • Perdew et al. (2008) J. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou,  and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
  • Lee et al. (2010) K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist,  and D. C. Langreth, Phys. Rev. B 82, 081101(R) (2010).
  • Vinet et al. (1986) P. Vinet, J. Ferrante, J. R. Smith,  and J. H. Rose, J. Phys. C Solid State 19, L467 (1986).
  • Dion et al. (2004) M. Dion, H. Rydberg, E. Schröder, D. C. Langreth,  and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
  • Klimeš et al. (2012) J. Klimeš, D. R. Bowler,  and A. Michaelides, Phys. Rev. B 83, 195131 (2012).
  • Cooper (2010) V. R. Cooper, Phys. Rev. B 89, 035412 (2010).
  • Berland and Hyldgaard (2014) K. Berland and P. Hyldgaard, Phys. Rev. B 89, 035412 (2014).
  • Sabatini et al. (2013) R. Sabatini, T. Gorni,  and S. de Gironcoli, Phys. Rev. B 87, 041108(R) (2013).
  • Vydrova and Voorhis (2010) O. A. Vydrova and T. V. Voorhis, J. Chem. Phys. 133, 244103 (2010).
  • Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  • Lee et al. (1988) C. Lee, W. Yang,  and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  • Bludský and Rubes̆ (2008) O. Bludský and M. Rubes̆, Phys. Rev. B 77, 092103 (2008).
  • David et al. (1992) W. I. F. David, R. M. Ibberson, G. A. Jeffrey,  and J. R. Ruble, Physica B 180, 597 (1992).
  • Nakamura and Miyazawa (1969) M. Nakamura and T. Miyazawa, J. Chem. Phys. 51, 3146 (1969).
  • Lu et al. (2009) D. Lu, Y. Li, D. Rocca,  and G. Galli, Phys. Rev. Lett. 102, 206411 (2009).
  • Reilly and Tkatchenko (2013) A. M. Reilly and A. Tkatchenko, J. Chem. Phys. 139, 024705 (2013).
  • Ringer and Sherrill (2008) A. L. Ringer and C. D. Sherrill, Chemistry 14, 2542 (2008).
  • de-la Roza and Johnson (2012) A. O. de-la Roza and E. R. Johnson, J. Chem. Phys. 137, 054103 (2012).
  • Azadi and Foulkes (2013) S. Azadi and W. M. C. Foulkes, Phys. Rev. B 88, 014115 (2013).
  • Drummond et al. (2008) N. D. Drummond, R. J. Needs, A. Sorouri,  and W. M. C. Foulkes, Phys. Rev. B 78, 125106 (2008).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description