The valence force-field model (VFM), initially proposed by Born and co-workers more than a century ago,1 is among the most useful approaches to deal analytically with the short-range valence forces in tetrahedrally coordinated crystals. In the VFM, the interactions in such systems are decomposed into bond-stretching and bond-bending forces. For solids with strong covalent bonding, where atom pair bonds play an essential role, the model is able to give a concise description of elastic properties. For example, the VFM has been most successful in describing the interactive forces in diamond and diamond-like structures,2 and it has been used to derive explicit formulas for the elastic constants and further relate them to the triply degenerate F_2g optical mode,3

ωVFM=8(k1+4k2)3M1/2,(1)where ω_VFM is the frequency of the diamond F_2g mode, k₁ and k₂ are the bond stretching and bond bending force constants for the harmonic part of the potential for diamond vibration, and M is the average atomic mass. In this equation, the anharmonic part of the carbon vibrations is ignored. Neglecting the bond bending term 4k₂ and expanding the free-energy terms of force constant, Vogelgesang et al.3 simplified Eq. (1) as follows:ωVFM=8KaM1/2,(2)where K is the bulk modulus and a is the lattice parameter of the diamond cubic structure. Using the VFM together with results from both Brillouin scattering and Raman spectroscopy, it has been possible to synchronize the optical mode frequency [the left hand side of Eq. (2)] and the elastic properties [the right-hand side of Eq. (2)] of diamond, silicon, and germanium under ambient conditions.

The coupling of elastic and optical properties in the framework of the VFM is critical for establishing an ab initio rule for the pressure parameter, which is fundamental for high-pressure sciences. In its simplest form, pressure p is the amount of force F applied perpendicular to a surface area S:p=FS.(3)However, neither F nor S can be easily measured inside the confined space of a pressurized apparatus. In addition, a limited variety of probes are available for such measurements. One possible approach is to use the integral of the bulk modulus with respect to volume, which follows from the relationdpdV=−KV,(4)although this relies on the possibility of calculating K directly from the Raman frequency of a pressure standard (e.g., diamond) at the corresponding volumes.

Deriving a high-quality, absolute primary scale is critical for high-pressure studies but remains a challenging task.4,5 Very early pressure scales were provided by dynamic compression using shock waves together with the equation of state (EOS) of NaCl.4 However, although shock-wave experiments achieve volume compression, there are also changes in internal energy, in addition to isentropic compression work corresponding to the observed density, leading to the so-called thermal pressure. This technical challenge has recently been partially overcome by employing multishock compression and precisely controlled compression wavelets such that shocks are delayed until after the experiment, in which there is forced compression along an adiabatic path.6 For static compression, Zha et al.7 constructed a primary pressure scale on the basis of Eq. (4) using a combined high-pressure Brillouin scattering and Raman spectroscopy system for polycrystalline MgO up to 55 GPa. This pressure range was later extended to 120 GPa by Murakami and Takata,8 using single-crystal MgO. The Brillouin scattering method provided the key information about elastic properties that was used for integrating the pressure. However, this method encountered difficulties with increasing pressure, because of the signal attenuation from Brillouin scattering spectroscopy and possibly the limited phase stability of MgO, which transforms to the B2 phase in the multi-megabar range.9 There have also been attempts to use theoretical simulations to construct primary pressure scales based on cubic SiC,10 despite its tendency to become structurally unstable above 40 GPa.11 The abovementioned experimental limitations have motivated us to develop a robust primary pressure scale, preferably based on diamond, one of the hardest known materials,12,13 under multi-megabar pressures.14

This work evaluates the validity of the VFM in diamond up to the terapascal pressure level, with the aim of establishing a primary pressure scale through an analytical relation between elastic and optical properties, namely, the K–ω relation. Both properties are calculated by first-principles computations, rather than extrapolation. Specifically, the elastic constant matrix of diamond at various volumes is obtained from the classic strain–stress relation. The gamma-centered optical mode ω_LTO is calculated by perturbing carbon atoms, i.e., through density functional perturbation theory (DFPT). Using the VFM, K is also calculated by rearranging Eq. (2) toKVFM=MωLTO28 a.(5)

For validation, the result is then compared with the bulk modulus obtained from the lattice stress–strain relation, K_LATT, at various volumes:KLATT=1+σKVFM,(6)where σ is the theoretical uncertainty of the bulk modulus using VFM and will be evaluated as a function of volume. We also perform an integration to calculate the pressure and estimate the total uncertainty at the end of the study. The theoretical K–ω relation obtained here is a preliminary result, but it provides a basis on which opportunities can be explored to establish a practical primary pressure scale through static compression experiments.

Our computation was implemented on the basis of density functional theory (DFT) and the projector augmented wave (PAW) method15 in the Vienna ab initio Simulation Package (VASP), version 5.4.4.16 The exchange-correlation energy was parameterized by a revised Perdew–Burke–Ernzerhof functional for solids (PBEsol).17 The core radius for carbon was taken as 1.128 Å. The Brillouin zone was sampled at a rate of 0.2 Å⁻¹, and the plane-wave basis set cutoff energy was taken as 1000 eV. The choices of a hard pseudopotential and a large energy cutoff were found to be sufficient to converge the interatomic force to within 5 meV/Å. We performed calculations for both unit cell and supercell (2 × 2 × 2 and 3 × 3 × 3) diamond to check the size effect. Calculations were conducted at various volumes with minimum V/V₀ = 0.507. Besides PBEsol, we also employed the PBE,15 Perdew–Wang 91,18 and local density approximation19 approaches to describe the exchange-correlation energy for the same calculation, and they achieved good agreement in the K–ω relation, as shown below.

We calculated the 0 K isothermal elastic constant matrix using the strain–stress relation20 and the bulk modulus of diamond, K_LATT, using the Voigt–Reuss–Hill averaging scheme.21 On the other hand, we calculated the optical mode of diamond through the DFPT method, perturbing the carbon atoms and constructing a Hessian matrix consisting of the second derivatives of energies. This was on the basis of a quasi-harmonic approximation. We then used the VFM [Eq. (5)] to calculate the K_VFM of diamond. We should note that calculation of K_VFM by DFPT and the VFM is fundamentally different from the method of strain–stress relation (giving K_LATT), the inputs to which are the responses to external strains, not vibrational properties.

It is worth mentioning that each diamond structure was optimized for atomic position and cell variables at fixed volume (or V/V₀). Results used for compression experiments should consider pressure correction of the elastic constant, which has been well documented.22,23

The diamond structure is optimized at various volumes. For example, V/V₀ = 0.507 corresponds to 1 TPa according to the Pulay stress. However, we intentionally discarded the parameter p as an input throughout our calculation, because our specific goal was to construct a primary pressure scale. The size of the supercell is known to affect simulation results,24 and we performed the same set of simulations for a selection of supercell sizes. Before proceeding to compressed volumes, we checked our calculated bulk modulus and shear modulus G at ambient pressure. The results of K = 450.6 GPa and G = 534.1 GPa matched well with those from the experimental EOS [K = 446(1) GPa],25 the Brillouin scattering method [K = 445(1) GPa and G = 538(1) GPa],3 ultrasonic determination [K = 442(5) GPa and G = 538(6) GPa],26 and a previous first-principles simulation,27 and thus our theoretical approach is capable of producing good-quality results. In Fig. 1, we plot the value of the uncertainty σ as defined in Eq. (6) vs the change in volume and compare the data with experimental results from Occelli et al.25

In principle, our first-principles simulation yielded a reasonable uncertainty σ. The size of the supercell has a moderate effect on the result. For example, in a unit cell of diamond, the σ values converged to ∼0.13 at more compressed volumes, whereas the calculations using supercells exhibited a turning point of σ at V/V₀ around 0.73. The inconsistency may come from the contribution of long-wavelength acoustic phonons, which was overcome by choosing a larger simulation supercell.28 In general, supercells greater than 2 × 2 × 2 only resulted in minor improvement. In the rest of this paper, without specifying the supercell size, we only show results from the 3 × 3 × 3 supercell.

Another important source of error may stem from the pseudopotential used in the simulations, which is intrinsic. This can be cross-checked against the experimental optical mode frequency at V₀, which has been reported by various studies.3,29 This error applies to the values calculated from both the VFM and DFPT. For example, at ambient pressure, the F_2g mode frequency of natural diamond is 1332.4 cm⁻¹, while our first-principles calculation gave 1349 cm⁻¹ with the VFM and 1325 cm⁻¹ with DFPT. These calculated values are usually tolerable for comparing experimental and simulation results. However, the error needs to be taken into consideration in the construction of a primary pressure scale, since it will increase the uncertainty in bulk modulus at more compressed volumes. Therefore, we shifted ω_VFM to be equal to ω_LTO at V₀ such that at ambient pressure, the uncertainty σ was normalized to 0 (the red and purple curves in Fig. 1). By shifting the zero-point frequency, the maximum uncertainty in bulk modulus using the VFM is 5.4%. Our results suggest that even with the bond bending term neglected, the VFM can still describe the elasticity and bond vibration of diamond with reasonably small uncertainty up to terapascal pressures.

It has not escaped our notice that the σ from the fitted experimental data of Occelli et al.25 diverged and σ blew up to greater than 0.24 at relatively high V/V₀ ratio (the orange curve in Fig. 1). In deriving the VFM equation, Vogelgesang et al.3 expanded the free energy F in terms of volume [Eq. (7) in Ref. 3]:ΔF=12(ΔV)2∂2F∂V2+⋯≈(ΔV)2B2V,(7)where ΔF is the change in Helmholtz free energy due to a perturbation ΔV of the volume. This expansion is no longer valid under stress, because the second-order elastic constants should be corrected for the pV term, and the expression for the bulk modulus will consequently change.22,23 The simplest example for cubic diamond under hydrostatic pressure is as follows:c11=1V∂2G∂e12T=0,c12=1V∂2G∂e1∂e2T=0+p,K=c11+2c123,(8)where G is the Gibbs free energy and e_i is the strain. Using the raw data points (volume and bulk modulus) from Occelli et al.25 and a revised formula for K_VFM, the VFM is able to reproduce the K_LATT from the EOS (the blue circles in Fig. 1) with an uncertainty of 4% or less. The experimental performance is even better than that of simulation within the same volume range. It is therefore important to process the bulk modulus data with caution when interpreting experimental results.

The most challenging part of establishing a primary scale is to obtain the elasticity of diamond under high pressure. Instead of directly determining the absolute K, the VFM can, through Eq. (2), link elastic properties with optical properties determined by Raman spectroscopy and x-ray diffraction, which have become routine probes for samples at multi-megabar pressures in diamond anvil cell (DACs). We plot the F_2g mode frequency of diamond vs K in Fig. 2(a), which is an important step before performing the pressure integral [Eq. (4)]. It is then necessary to perform error analysis to examine the validity of the K–ω relation. Figure 2(a) shows that ω reaches its maximum theoretical uncertainty (3.7%) when K_LATT ≈ 1100 GPa, equivalent to a volume of ∼4 ų/atom. This is due to the crossover of KLATT′ and KVFM′, indicating that according to the VFM, diamond is more compressible than predicted by DFT. In other words, the simplified VFM starts to underestimate C–C interactions beyond this critical point, although this simplification is still able to predict values of the bulk modulus with reasonable uncertainties. We fit our DFT results using a second-order polynomialω=a K2+bK+c,(9)with a = −0.591(1) × 10⁻⁵ (cm GPa²)⁻¹, b = 0.667(2) (cm GPa)⁻¹, and c = 1051(1) cm⁻¹. The fitting yields an R² of 0.9997. This is a robust relation, which is within reasonable error of those obtained by the VFM and is almost unbiased by the choice of pseudopotential [Fig. 2(b)].

We integrate our data according to Eqs. (4) and (8) and plot the diamond Raman peak vs accumulated pressure in Fig. 3. The results are compared with those obtained by Fratanduono et al.6 and Akahama and Kawamura29 using diamond Raman edge gauges and with those obtained by Dubrovinskaia et al.30 using a polycrystalline diamond gauge. We note that both of the diamond Raman edge gauges6,29 divided the pressure into two regimes and employed different functions for each of these. With the diamond gauge used by Dubrovinskaia et al.,30 the measurements were performed on a diamond chip inside the sample chamber, rather than on the diamond anvil itself. Our results show excellent agreement with those from the diamond chip Raman gauge within the experimental pressure range (0–270 GPa) and good agreement (<3% difference) with both of the diamond Raman edge gauges below megabar pressures. However, at higher pressures, our scale gives lower frequency values. For instance, extrapolation from the diamond Raman edge gauges overestimates the pressure at above 150 GPa and completely diverges, with >50% difference, when the pressure is above 500 GPa.

We should note that the divergence stems from the methodology rather than the experimental procedure, although all of the above-mentioned methods use diamond as the calibrant. A diamond Raman edge gauge measures the first derivative of the diamond Raman peak at the center of the diamond anvil, the tail shape of whose spectrum is generated by the pressure gradient over the anvil culet. The cupping and even double cupping behavior of a beveled diamond anvil will generate a complex pressure distribution on the culet surface and thus distort the spectrum of the diamond peak at more compressed volumes.31 In the case of a primary pressure scale, the diamond will be put into the sample chamber,30 and the observed Raman spectrum will have a single peak corresponding to the triply degenerate F_2g mode, similar to that found by Dubrovinskaia et al.30 Differences in the ways in which experiments are performed have led to substantial deviations when the pressure is close to or beyond the deformation limit of single-crystal diamond.31

The plot of pressure residual vs absolute (integrated) pressure in Fig. 4 reflects the excellent match between the integrated pressure and the standard pressure p_std. In this work, p_std is regarded as an independent checkpoint for the VFM. For example, in our computations, p_std is taken to be the Pulay stress, which is the diagonal component of the stress tensor. In the work by Occelli et al.,25 the secondary ruby gauge pressure was chosen as p_std. The validity of the VFM and our proposed K–ω relation is verified by the experimental results of Occelli et al.25 (inset of Fig. 4) and theoretical calculations. Up to ∼1 TPa pressure, the pressure residual is no greater than 10 GPa, which is particularly useful for calibrations in the multi-megabar pressure range.

A primary scale is rarely used in high-pressure experiments. Instead, ruby, Raman, and x-ray gauges have been used over greatly extended pressure ranges by extrapolating literature data or achieving pressure self-consistency,5,32 however, even so, they must ultimately be tied to a primary scale. This leaves room for systematic discrepancies when different standards are used, particularly at elevated pressures and temperatures.33 Therefore, a high-precision primary scale up to multi-megabar pressures is an urgent requirement for ultrahigh-pressure experiments. Diamond has a number of advantages as a basis on which such a scale can be established. It features very strong C–C covalent bonding, which is ideal for implementing the VFM. It is also thermodynamically stable throughout the pressure range achievable in current static compression experiments.14,34 The diamond Raman edge gauge has been de facto among the most popular calibrants for high-pressure experiments because of the popularity of DAC.35,36 A primary pressure scale for diamond will not only satisfy the requirements for ultrahigh-pressure experimentation but also streamline pressure calibrations for all DAC users.

The theoretical pressure scale using diamond and the VFM agrees well with previous primary scales7,8 and also with the latest shock compression experiments6 up to megabar pressures, but exhibits substantial differences at higher pressures (Fig. 3). These sizable inconsistencies in pressure will not reduce the importance of recent highlights in ultrahigh-pressure work, including hot superconductors37,38 and dense metal-like hydrogen,39 although they may lead to lower values of their onset pressures and further constrain the phase stability fields. However, large overestimates of pressure will be a hurdle in interpreting experimental data in the earth sciences. For example, pressures in Earth’s core cover a wide range of 140–360 GPa, and estimation of light element content relies heavily on the structure and the EOS of iron.40 Our results suggest that the overestimated pressures given by the most widely used pressure gauges will indicate a much “lighter” Earth’s core, and, as a consequence, existing models will suffer considerable density deficits. Therefore, an additional number of volatile elements may need to be incorporated into the conventional light-element-free inner core, leading to a revision of the currently accepted large-scale element budget and geochemical condition of the core.41

In short, a robust primary pressure scale that extends to multi-megabar pressures is vital for the next generation of experiments under extreme conditions. The classic VFM is a powerful tool to verify the K–ω relation (Fig. 2) at terapascal pressures, and theoretical uncertainties will not accumulate with increasing pressure. Our study lays the foundations for a path to establish a primary pressure scale using currently available structural and optical probes, which could contribute to many aspects of high-pressure sciences.

Acknowledgment. We thank Howard Sheng, Hanyu Liu, and Wenzhong Wang for helpful discussions. Q.H. is supported by the CAEP Research Project (Grant No. CX20210048) and a Tencent Xplorer prize. The work is also partially supported by the National Natural Science Foundation of China (NSFC) Grant No. U1930401.