The properties of materials under extreme conditions are of fundamental importance in many areas of science and engineering.1–4Ab initio methods based on density functional theory (DFT)5,6 have played a key role in the study of matter under high pressure. Among these, the projector augmented wave (PAW) method7,8 developed by Blöchl has been widely used for high-throughput ab initio DFT studies, and has been implemented in many codes.9–12 In the PAW method, core electrons are “frozen” and replaced by a pseudopotential to obtain a PAW potential. There are various versions of PAW potentials, differing in the ways in which they treat valence electrons. We take as the standard PAW potential the one that treats only the outermost electrons as valence electrons. Other types of PAW potential can be regarded as extended versions of this standard potential.

The PAW potential is crucial to the accuracy of calculations,13 and it therefore needs to be selected with care. There have been a number of detailed comparative studies, including error analyses, of multiple methods, including a PAW method for ground-state elemental crystals.14–16 Some studies have also tested and optimized PAW datasets for rare-earth elements17 and transition metals.18 In these studies, the tests on the PAW potentials have been carried out at zero temperature or ambient pressure, thereby providing guidance for calculations under nonextreme conditions. However, under high pressure, the electronic structure of a material may change significantly compared with that under ambient pressure. One effect that has been recognized is that the transfer of conduction electrons between the s–p states and d states depends on the degree of compression and varies among crystal structures.19–21 The conduction electrons are closely related to the structural stability of a material. For most 3d transition metals, although the standard potential can give reasonable results at low pressure, it should be used with the uttermost care in high-pressure calculations.22 Using a PAW potential that cannot correctly describe the distribution of electrons when a material is compressed will predict erroneous thermodynamic properties under high pressure. Therefore, choosing a version with semicore electrons also included in the valence states is much safer.

However, in practical applications, besides accuracy, computational efficiency is also a crucial issue that must be taken into account. Especially when faced with large-scale computing tasks, such as ab initio molecular dynamics (AIMD) simulations of the melting of transition materials, we must strike a balance between accuracy and efficiency, and sometimes have to choose the cheapest method within a certain tolerance of errors. Taking the melting simulation of vanadium as an example, it is much cheaper to choose the PAW potential containing only five valence electrons (labeled V) than the version containing 13 valence electrons (labeled V_sv). Thus, the PAW potential V should also be tested as a candidate. Up to now, V_sv has been consistently used in DFT simulations of vanadium under high pressures.23–27 However, there have been few systematic tests of the performance of the V potential or discussions of the influence of the different PAW potentials for vanadium on simulations under extreme conditions. The effects of semicore electrons on the calculation of high-pressure thermodynamic properties of vanadium remain unclear, and there is as yet no definitive standard choice of PAW potential for simulation under high pressure.

In this work, the effects of different PAW potentials on the simulation of thermodynamic properties of vanadium under high pressure are investigated. Comparisons of various properties calculated using different PAW potentials are performed. In particular, the capability of the potential V at increased pressures is systematically tested by comparing its results with those from V_sv. Our purpose is to elucidate the origin of the differences between the results obtained with the different PAW potentials and finally to recommend a procedure for testing these potentials to guarantee reasonable results for thermodynamic properties of materials under extreme conditions.

The remainder of the paper is organized as follows. A brief introduction to the methods and techniques used in this work is given in Sec. II. Then, in Sec. III, the effects of different PAW potentials on the properties of body-centered cubic (bcc) vanadium simulated at zero temperature (zero-T) (e.g., cell parameters, bulk modulus, zero-T compression isotherm, energy band, and electron state density) and finite temperature (finite-T) (e.g., melting point and vibrational density of states) are comprehensively compared and discussed. Finally, our conclusions are summarized in Sec. IV.

In this study, various zero-T and finite-T tests are carried out employing the Vienna Ab initio Simulation Package (VASP).11,28–30 The standard VASP-PAW potential of vanadium, labeled V, contains only five electrons (4s², 3d³) as valence electrons. If the inner electrons of 3p and 3s3p are also treated as valence electrons, this gives versions with semicore states included, labeled V_pv (3p⁶, 3d³, and 4s²) and V_sv (3s², 3p⁶, 3d³, and 4s²). The cutoff radii of V, V_pv, and V_sv are 1.428, 1.217, and 1.217 Å, respectively. To obtain the all-electrons (AE) results as reference, some calculations are repeated with the full-potential (linearized) augmented plane-wave and local-orbitals (FP-LAPW+lo) method31,32 using the WIEN2k package.33,34

In the zero-T tests, different PAW potentials are tested in two commonly used approximations for the electronic exchange-correlation (XC) functional: the Ceperley–Alder parametrization local density approximation (LDA)35–37 and the Perdew–Burke–Ernzerhof (PBE)38 generalized gradient approximation (GGA).39–41 Each calculation in this part is performed on a primitive bcc cell. The k-point grid is 17 × 17 × 17 for the equilibrium structure, and the number of k-points increases when the spacing between the two atoms is decreased. The cutoff energy is 260 eV for the standard PAW potential V and 360 eV for the semicore potentials V_pv and V_sv. These settings ensure energy convergence to better than 1 meV/atom.

In the finite-T tests, the thermal electronic excitations are included using Mermin’s DFT scheme.42 The melting properties are calculated based on the PBE approach. Some experimental43 and theoretical44,45 studies on the solid–solid phase transition of vanadium show that from ambient pressure to about 300 GPa, the bcc phase is the only stable solid phase of vanadium before melting occurs. Therefore, our melting simulations are performed on bcc cells. Each simulation cell, arranged by conventional cells into 5 × 5 × 5, contains 250 atoms. Only the Γ point in reciprocal space is used, and energy cutoffs are correspondingly the same as those of zero-T calculations.

The AIMD simulations of melting are implemented using the Z method.46,47 To achieve high accuracy, a standard simulation by the Z method needs thousands of atoms and is very time-consuming, especially at high pressures. Since the subject of this paper is the effect on calculations of thermodynamic properties of the use of different PAW potentials, it is sufficient to obtain some rough melting points rather than the exact ones on the phase diagram. A simulation cell containing 250 atoms and a simulation time of 3 ps for each melting simulation are used in this work. For more accurate results on the melting curve of vanadium, see our previous work24 on solid–liquid coexistence simulations and the results of Errandonea et al.25 using the standard Z method.

In this subsection, we show a series of zero-T properties of vanadium.

First, different PAW potentials and XC functionals are tested by calculating the zero-pressure properties of bcc vanadium, such as the equilibrium volume V₀ and bulk modulus B₀. Our results and reference data26,48–52 are listed in Appendix A in the supplementary material. The results show that at the level of the XC functional, the PBE functional is a better choice, since it gives a much more accurate result for V₀. For comparison, we perform an AE calculation with the FP-LAPW+lo method, in which the semicore electrons (3s3p electrons) are treated as valence electrons. The parameter RK_max, which determines the number of basis functions, is set to be 7. Unsurprisingly, the results using the PAW potential V_sv are almost the same as the AE results, and both are the closest to the experimental values. Meanwhile the variance of the results with the potential V is not notable, since the differences between these results V and those from V_sv are within 3%.

The zero-T compression isotherms of vanadium are calculated to test the performance of the PAW potentials under pressure. The results are plotted in Fig. 1(a) together with the FP-LAPW results and the data from the reference thermodynamic table of Young et al.53 As can be seen from Fig. 1(a), the LDA results differ from those of Young et al. in the low-pressure region. PBE gives better results than LDA. It indicates that the electron distribution of bcc vanadium has a degree of anisotropy and deviates from a uniform electron gas. We therefore focus our attention on the PBE results. To make the contrast more obvious, we define a ratio at each volume, which reflects the relative deviation of the PAW result from the FP-LAPW result:$$R_i=\frac{\mathrm{\Delta }{P}_i}{P_{\text{FP}-\text{LAPW}}},$$with$$\mathrm{\Delta }{P}_i=P_{{\text{PAW}}_i}-P_{\text{FP}-\text{LAPW}},$$where P is the value of the cold pressure, and the subscripts FP-LAPW and PAW_i indicate the results from FP-LAPW and the different PAW potentials. Figure 1(b) shows the ratio R_i as a function of volume. We find that when the volume is compressed to about 13 ų (a volume compression factor of about 0.96), the deviation caused by using the PAW potential V is almost 20%, which implies that V is unreliable even at low pressure. When the system is compressed to a volume below 8 ų (a volume compression factor of about 0.6), the results from the potential V clearly diverge from the AE reference results. As the volume decreases, the relative difference increases rapidly. Meanwhile, the potentials V_sv and V_pv still maintain a reasonable deviation of within 10% when the volume is compressed to less than 6 ų per atom.

Since the result from V_pv is close to that from V_sv, in order to simplify the problem and reveal more clearly what is going on, the following comparison is only performed between V and V_sv. As we have mentioned, the main difference between V and V_sv is in the treatment of the valence electrons. To investigate the effect of the valence electrons in the high-pressure calculation, we implement another FP-LAPW calculation, in which the 3s3p electrons are treated as core electrons. As can be seen from Fig. 1(c), the treatment of the 3 s and 3p states has a significant influence on the FP-LAPW calculation: the FP-LAPW result without restriction of 3s3p is very close to the PAW result using V_sv: when the 3s3p electrons are restricted to the core during the FP-LAPW calculation, the slope of the cold pressure curve obtained becomes flatter and closer to the result from the PAW potential V. This indicates that the underestimation of the cold pressure by the PAW potential V is due to the improper definition of the valence electrons of vanadium under high pressure. The inner electrons of 3s3p must also be treated as valence electrons.

The electronic structures of systems under high pressure are also compared between the PAW potentials V and V_sv. The energy bands and the corresponding electronic density of states (DOS) of bcc vanadium at different volumes are calculated along the specific k-points (Γ–H–N–Γ–P) in the Brillouin zone. The results for two volumes, namely, 9.84 and 6.91 ų, are shown in Fig. 2. Appendix B in the supplementary material shows the orbital projection of energy bands and partial DOS. We must first clarify that the sp states in Fig. 2 do not refer to the 3s3p, which should lie in a much lower energy range. However, the difference induced by the treatment of 3s3p electrons can be observed from the spd states around the Fermi level, as shown below. At a volume of 9.84 ų, the band structures calculated using these two PAW potentials are almost indistinguishable. As the pressure increases, the valence electrons delocalize further, resulting in a lower electron density at the Fermi level. The spatial extension of the s–p orbitals is greater than that of the d orbitals, which leads to upward shifts of the s–p bands relative to the d bands and finally to the formation of a pattern of inversion of the s–p and d bands. Both V and V_sv describe the increasing band dispersion and the characteristic of band inversion caused by pressure-induced electron transfer between the s–p and d states, but there are obvious differences between the results from the two potentials at higher pressures. At a volume of 6.91 ų, the band structures described by V_sv show much wider dispersion at the valence band compared with those described by V. The enhancement of band dispersion even leads to some topological difference from the results with V. For example, the electron hole at the k-point along the Γ–H path (at about three-quarters of Γ–H) and the shallow electron pocket at Γ show the qualitative difference between the V and V_sv results in terms of metallic properties. This difference may result in a substantial divergence of thermal stability, melting points, or other important thermodynamic properties between V and V_sv. A possible reason for this is that the semicore electrons in V_sv provide an extra degree of freedom to describe the anisotropy of the potential in reciprocal space. With a more detailed description of the potential in reciprocal space, the energy states at different k-points are found to couple much more tightly.

Besides the different anisotropies in reciprocal space between the V and V_sv results, the difference in local potential polarization is also investigated here. The interatomic forces are calculated using V and V_sv. As can be seen in Fig. 3, below the equilibrium bond length, the potential curve u(r) and corresponding force F(r) obtained using V_sv have steeper slopes than those obtained using V. For systems with the same average interatomic distance r (compressed to the same volume), the probability of collisions between atoms calculated using V_sv is greater than that calculated using V, which results in a higher pressure of the system. This is consistent with the results for the zero-T compression isotherms. The weak force at short range observed in the results from V indicates that the polarization of electrons obtained using V is different from the strong bond-like structure obtained using V_sv.

To further clarify the essential cause of the discrepancies in the results obtained using different PAW potentials, the differences in charge densities are compared, and the results are shown in Appendix C in the supplementary material. Since the contrast in direct charge density is not obvious, the charge density differences (CDDs) of systems with different degrees of compression described by the PAW potentials V and V_sv are calculated and shown in Fig. 4. The CDD is given by the difference between the charge density of a vanadium crystals and the superposition of the charge densities of the neutral vanadium atoms, and clearly reveals the characteristics of the electrons when vanadium atoms form a vanadium crystal, such as the charge transfer properties and the direction of bonding polarization in the process of bonding and coupling of bonding electrons. However, the CDDs shown here do not represent the actual charge distribution in the atomic sphere, since any PAW potential will always give the wrong electron density inside the core region. Comparing the CDDs calculated using the potentials V and V_sv, we find that the electrons in the bcc vanadium described by V tend to be more discretely distributed between atoms, and their spatial distribution is less oriented. The more homogeneous charge distribution in the results from V may lead to weaker bonding between atoms compared with that obtained from V_sv. This phenomenon is probably caused by the isotropic potential of the semicore electrons in the pseudopotential V. In contrast to V, in V_sv, the 3s3p electrons are treated as valence electrons, which makes the potential fully polarized, with an anisotropic crystal environment. As can be seen from our results, this polarization becomes more obvious and more significant as the pressure grows. Therefore, we speculate that under high pressure, the insufficient bonding strength described by the potential V will give rise to premature deformations of the system before it reaches its true melting temperature, which will result in an inaccurate calculation of the melting properties of vanadium. This will be verified in the following finite-T simulations.

In this subsection, several finite-T properties of bcc vanadium at high pressure are calculated using the PAW potentials V and V_sv, and the results are compared.

The high-pressure melting curve is one of the most important properties of vanadium. The results obtained using the PAW potential V are compared with those from previous calculations24,25 using V_sv. Taking the results for volumes of 12.84, 9.84, and 6.91 ų as examples, it can be seen from Figs. 5(a)–5(c) that the neighborhood where the isochoric curve bends represents the approximate melting point. The melting pressures and temperatures calculated using V_sv are consistent with the Z-method results from Errandonea et al.25 and our coexistence results.24 However, the melting points calculated using V are very different from these results.

In the case of 12.84 ų (a volume compression factor of 0.95) [Fig. 5(a)], the melting temperature T_m and pressure P_m calculated using the PAW potential V are already known to be underestimated by about 20%. Hence, it can be seen that, even at low pressures, using V will lead to misleading results for the melting point of vanadium. As the compression ratio continues to increase, the underestimations of P_m and T_m become more pronounced. When the volume is compressed to 6.91 ų [Fig. 5(c)], the pressure and temperature at the melting point calculated using V are respectively about 140 GPa and 4000 K lower than those calculated using V_sv. More alarmingly, in this case, the potential V gives an abnormal relationship in which the pressure decreases with increasing temperature, resulting in a melting curve with a zero or even negative slope dT_m/dP_m at high pressure. According to the Clausius–Clapeyron relation dT_m/dP_m = V^ls/S^ls (where V^ls is the volume and S^ls is the entropy of fusion, which is always positive), the negative slope of the melting curve indicates that the volume of the system is reduced after melting, which contradicts experimental results.25

The calculated thermodynamic quantities are also compared between the two PAW potentials. It can be seen from Fig. 5(d) that for systems of the same volume and temperature, the equilibrium pressures calculated using the potential V are always lower than those calculated using V_sv. This is consistent with the calculated zero-T compression isotherms and is due to the different descriptions of forces in V and V_sv (as shown in Fig. 3).

The vibrational density of states (VDOS) is also calculated. The VDOS is defined by the distribution of vibrational normal modes of a system.54 It is the frequency Fourier transform of the velocity autocorrelation function, which can be obtained from temporal trajectories of particles from MD simulations. Figure 6 shows a comparison of the VDOS obtained using the PAW potentials V and V_sv for the same simulation temperatures and volumes. At a temperature of 4000 K, for both cases with volumes of 9.84 and 6.91 ų, the VDOS given by V_sv still corresponds to the solid state, while the system described by V has already melted. This supports the results from the melting simulation, in which the potential V underestimates the melting point of vanadium under high pressure.

There is another important phenomenon that needs attention in the VDOS results. It stands to reason that after the system has been compressed from 9.84 to 6.91 ų, the distribution of atomic vibrational frequencies in reciprocal space should move toward higher frequencies overall, owing to the shortening of the lattice. According to the results from V_sv, after the volume has been compressed, all the peaks move about 200 cm⁻¹ toward higher frequencies. However, the VDOS given by the potential V does not exhibit such behavior. By comparison with the results from V_sv, the interatomic force described by V is weaker, and the atomic vibrational frequency is softer. Thus, it is confirmed that the bonding strength given by V is too weak.

From the results of the zero-T and finite-T simulations, we have found the following.

First of all, the high-pressure, high-temperature properties of vanadium are sensitive to the choice of PAW potential. Using the PAW potential V will give a cold pressure deviation of over 10%, resulting in an underestimation of the melting point by more than 20% at a given compression ratio. The errors in thermodynamic properties at high temperature and high pressure calculated using different PAW potentials can be predicted by the deviations of the zero-T compression isotherms. We suggest that PAW potentials should first be evaluated by testing zero-T properties under high pressure in the following steps:1.The PAW radii in the PAW potential should be small enough for the problem of interest. The cutoff radii should be checked in advance to ensure that the volume fraction of overlap of PAW spheres is within an acceptable limit.⁵⁵2.The PAW method should be used with different PAW potentials and the AE method to calculate the zero-T compression isotherm of the system within the target pressure range.3.The AE result should be taken as the benchmark, and PAW potentials with an unacceptable deviation (e.g., >10% for vanadium) should be excluded.4.Among the remaining PAW potentials, the one that can reproduce the benchmark and at the same time has an acceptable computational cost for the subsequent calculation is selected.

Using this procedure, a pseudopotential suitable for simulations under high temperatures and high pressures can be selected.

Second, the energy bands and CDDs reveal that the discrepancies between the PAW potentials arise because of the differences in the definitions of the valence electrons in these potentials. A potential with semicore (3s3p) electrons as valence electrons gives much better results at high pressure from two aspects. On the one hand, at high pressure, the coupling between different electron states is much tighter, and this leads to a more obvious variation of the potential in reciprocal space. The additional semicore electrons provide an extra degree of freedom to enable a more detailed description of the anisotropy of the potential in reciprocal space. This is important for the band dispersion, which in turn affects the metallic properties of the system. On the other hand, the reduced distance between nuclei under high pressure leads to polarization of the distribution of electrons around atoms. The free semicore electrons can contribute to the formation of a local polarized potential, which gives a better description of bonding under high pressure.

The results for the VDOS corroborate the above conclusions. They verify that the insufficient bonding strength described by the potential V causes the melting point under high pressure to be underestimated.

We have reported simulations of the high-pressure properties of vanadium using different PAW potentials under both zero-T and finite-T conditions. We find that there is no significant distinction between the results from the potentials V_sv and V for equilibrium volume and bulk modulus at zero pressure, but there are striking differences between the results from these potentials for the zero-T compression isotherms and thermodynamic properties at high pressure. The calculations using the potential V badly underestimate the melting points of vanadium, and even give a qualitatively wrong characterization of the melting curve as having negative slope.

The high-pressure electronic structure, charge distribution, and interatomic force at zero temperature have been analyzed to explain the above phenomena. This analysis reveals that the discrepancies in thermodynamic properties calculated using different PAW potentials arise mainly from the differences in the electron dispersion and anisotropy of electronic bonding described by these pseudopotentials. These lead to a significant difference in lattice stability with increasing temperature.

Finally, and most importantly, we have proposed a procedure for testing PAW potentials to ensure that reasonable results are obtained for thermodynamic properties when simulations are performed under extreme conditions.

Acknowledgment. The authors thank Shijie Xiong for valuable suggestions. Our deepest gratitude goes to the editors and the anonymous reviewers for their careful work and thoughtful suggestions that have helped improve this paper substantially. This work is supported by the Science Challenge Project (Grant No. TZ2016001) and the National Natural Science Foundation of China (Grant Nos. U1930401, 51671033, and 12004048). We thank the Tianhe platforms at the National Supercomputer Center in Tianjin and the foundation of LCP.