Exclusively Relativistic: Periodic Trends in the Melting and Boiling Points of Group 12

Abstract First‐principles simulations can advance our understanding of phase transitions but are often too costly for the heavier elements, which require a relativistic treatment. Addressing this challenge, we recently composed an indirect approach: A precise incremental calculation of absolute Gibbs energies for the solid and liquid with a relativistic Hamiltonian that enables an accurate determination of melting and boiling points (MPs and BPs). Here, we apply this approach to the Group 12 elements Zn, Cd, Hg, and Cn, whose MPs and BPs we calculate with a mean absolute deviation of only 5 % and 1 %, respectively, while we confirm the previously predicted liquid aggregate state of Cn. At a non‐relativistic level of theory, we obtain surprisingly similar MPs and BPs of 650±30 K and 1250±20 K, suggesting that periodic trends in this group are exclusively relativistic in nature. Ultimately, we discuss these results and their implication for Groups 11 and 14.

TABLE I. Valence-projector energies and ENMAX value (in eV) for the default relativistic and non-relativistic PBE-PAW potentials used in and made for this study. NR PAWs for Zn, Cd, and Hg use the exact same parameters as the default ones (12 valence electrons). The PAWs for Cn are taken from a previous study and include the s and p semi-core levels to a total of 20 valence electrons. Liquid and solid phases are represented by 64-atom configurations, except for liquid Hg, where we additionally consider a 121-atom configuration with very similar results. Previous studies have confirmed that these atom numbers provide numerically converged results. [6,8] The influence of electronic entropy is accounted for through Fermi-smearing of the populations. [9] Two reference levels of theory are used for each element. One with reduced accuracy that allows for efficient sampling of the configuration space, and another one to obtain the final numerically converged energies. The lower reference level of theory exploits reduced global precision (PREC = NORMAL) and an energy-convergence of 10 −4 eV (EDIFF = 1E-04) to accelerate the simulations and in in general scalar/non-relativistic (LSORBIT = FALSE). Liquid simulations conducted at this level use the Nose-Hover thermostat with an SMASS parameter of 2-4, solid simulations employ a Langevin thermostat with a LANGEVIN GAMMA value of 3. The higher reference level uses (PREC = ACCURATE), an energy-convergence of 10 −6 eV (EDIFF = 1E-06), and explicitly includes spinorbit coupling in the valence space (LSORBIT = TRUE, if applicable). The element specific parameters like kineticenergy cut-offs, k-point grids, and timesteps for the simulations are collected in Tab. II.   TABLE II. Element-specific parameters for the calculations (identical for NR, SR and SOR levels). The lower reference level used in the simulations (TDI) and the higher reference level for single-point energies (TPT) are given in the form (TDI→TPT). The timestep used in the normal simulations is given last, followed by the shorter timestep used in the simulations near the non-interacting limit (liquid only). element (phase Tsim) cut-off [eV] k-grid timestep [fs] Zn (solid 700 K) 300 → 600 2 3 → 4 3 3 Zn (liquid 1000 K) 300 → 600 2 3 → 4 3 3 (1) Cd (solid 500 K) 300 → 600 2 3 → 4 3 4 Cd (liquid 1000 K) 300 → 600 2 3 → 4 3 4 (2) Hg (solid 250/500 K) 250 → 500 2 3 → 5 3 6 Hg (liquid 300/700 K) 250 → 500 2 3 → 5 3 6 (2) Cn (solid 200 K) 400 → 600 2 3 → 5 3 8 Cn (liquid 250/300 K) 400 → 600 2 3 → 5 3 8 (2)

FREE-ENERGY CALCULATIONS
Solid -The free energy of the solids is calculated in four main steps, in the first of which the electronic energy and ionic frequencies of the super-cells are calculated at their respective equilibrium volumes either non-or scalarrelativistically. For the ionic frequencies, the numerical precision is increased by setting Precision=Accurate and LREAL=False in the INCAR file to avoid aliasing errors with the finite-differences method. To make these calculations more accurate, four displacements are considered for each atom in each direction (NFREE=4). Secondly, using the previously determined forceconstants, the vibrational contributions to the free energy are calculated in the harmonic approximation. These calculations are carried out with the Phonopy program using a very fine k-point mesh (16 3 ). [11] Thirdly, thermodynamic integration from the har- monic to the (non-or scalar-relativistic) DFT solid is carried out by interpolating between harmonic and DFT forces. The integral is evaluated using a Gauss-Legendre three-point rule with at least 5000 steps time at each point (1000 steps equilibration, 4000+ steps productive). This provides the anharmonic correction to the free energy of the harmonic crystal, which is the main source of errors for the free energy of the solid with a statistical uncertainty of 0.1 − 0.4 meV.
Finally, the difference between the low and high reference levels (spin-orbit coupling, k-point mesh, kineticenergy cut-off, and increased numerical precision, see Tab. II) is accounted for using TPT, i.e., by recalculating the electronic energy for at least 10 frames at both  where the index after the angle bracket indicates that the difference ∆U 1−2 is evaluated for configurations generated by H 1 . Instead of the exact equation, we use the second-order cumulant expansion which is sufficiently accurate since already the secondorder term is 1 meV/atom in all cases, indicating a parallel run of the respective potential-energy surfaces.
Statistical errors are < 0.2 meV in all cases. The contributions and final Gibbs free energies of the solids are collected in Tabs. III-VI.
Liquid -The Gibbs energy of the liquid is calculated through TDI from a non-interacting reference. For this purpose the difference of the internal energies is integrated along the coupling parameter λ The free energy of the non-interacting reference at the Here, N is the number of atoms, m their mass, V the atomic volume, and Θ the electronic degeneracy. Note that the pV term is negligible at the volume of the liquid ( 1 meV/atom) and only taken into account for the gas phase below. Since the kinetic energy part of U 1 and U 0 is identical it cancels, and the potential part vanishes at zero interaction stength (U 0 ), the value of the integrand is the average internal potential energy calculated at full interaction strength U pot 1 (R) for configurations R generated with reduced interaction strength (at the respective λ). This integral is evaluated using numerical quadrature in the form of a n-point Gauss-Lobatto rule, in principle requiring one NVT simulation for each value of λ. We use n = 6 for Zn and Cd, n = 8 for Hg and Cn. While most of these simulations are straightforward, the ones very close to the ideal-gas limit (λ 0.01 or < 1% of the DFT forces) are tedious, whereas the simulation for the end point λ = 0 is not possible with a PAW+DFT methodology. This is because close-encounters between the (almost) non-interacting atoms lead to singularity in the energy resulting in numerical instabilities in errors in the simulations, partly resulting from overlapping coreelectrons. An approach to circumvent these issues was devised and implemented by Kresse and coworkers and will be used here with slight modifications. [8] The approach is based on substituting λ in eq. (3) with λ(x) = ( x+1 2 ) 1/(1−κ) , which yields This introduces an explicit dependency on λ in the integrand, which not only dampens the impact of the technically challenging calculations near the non-interacting limit, and eliminates the point for λ = 0. The parameter κ guides the mapping of the quadrature points between the domains. While a value close to 0 retains the original (equidistant) spacing of the Gauss-Lobatto quadrature, κ close to unity increases the density of quadrature points in the λ domain in the region close to λ = 0, where the slope of f (λ) is the largest. Kresse and coworkers suggest κ > 0.8. In a recent work, [12] we demonstrated that smaller values can suffice and lead to much more manageable simulations. We use κ = 0.7 for all elements, having tested 0.8 for Cn (300 K) with very similar results. In analogy to the solid, TPT is used to correct and the final energies with a desired precision of ≈ 1 meV/atom, which requires about 10 snapshots from each simulation. Gas Phase -The Gibbs energy of the gas phase is calculated for the non-interacting (ideal) gas at its equilibrium volume and ambient pressure. For a given atomic degeneracy Θ, volume V = nkT /p, temperature T , particle number N and mass m this is (7) where Z is defines in eq. (5). For the gas phase, the logarithm in Z is solved using the Stirling approximation, which is sufficiently accurate since we are considering an arbitrary number of particles.
To check if the non-interacting model is sufficient, we evaluate the first virial (two-body) correction for each of the examples assuming a Lennard-Jones (12,6) potential with the parameters for the dimers provided in the introduction. This leads to the following integral which can be evaluated as described in ref. 13. Due to the weak interaction between isolated Zn, Cd and Hg atoms and their large R e s, the correction is repulsive and negligibly small at the respective BPs (≤ 0.1 meV/atom). At the distinctly lower BP of Cn, the virtial term is attractive and somewhat larger at −0.2 meV/atom, yet still too small to affect the calculated BP.

CALCULATIONS OF ERRORS
Statistical errors for averages taken from simulations are calculated using block-averaging. The typical block size for the solid configuration is about 50, and distinctly larger with 50-150 for the liquid depending on time-step and λ. Based on these data, we have chosen the length of the simulations to achieve a statistical error smaller than 1 meV/atom. The other significant source of errors is the linear extrapolation of G(T ) for the solid and liquid (cf. Fig. 3 in the main article). To mitigate these errors, we have conducted the free-energy calculations near the respective MP/BP. Altogether, this leads to an uncertainty of about 15 K in the predicted MPs and about 5 K in the BP, the latter of which is less sensitive due to the steeper intersection of the respective free-energy curves.

ORIGIN OF THE DIFFERENCE TO THE PREVIOUS ESTIMATE OF STEENBERGEN
We note that the interface-pinning method employed in ref. 14 requires much larger cells of 250 atoms (125 liquid and 125 solid atoms), which was only possible (despite using a large supercomputer) at a scalar-relativistic level with a rather coarse k-point grid (2 × 2 × 1 or effectively 2 × 2 × 2 points for the solid/liquid part as these are joined third axis). In contrast, our approach allows to use smaller 64-atom solid (4 × 4 × 4 super-cell) and 64-or 121-atom liquid configurations, which in combination with the incremental corrections via TPT allows to include much finer k-point grids of up to 5 × 5 × 5 as well as explicit spin-orbit coupling. Inspection of the respective TPT corrections from the 2 × 2 × 2 grid employed in our simulations, which corresponds to the final level of ref. 14, to the final 5 × 5 × 5 grid used here (cf. Table V) reveals a strong influence and large differences in the respective corrections for the liquid (< 1 meV/atom), relativistic solid (8 meV/atom), and non-relativistic solid phases (13 meV/atom). Although a detailed analysis of the influence this has on the predicted MP is not possible, the large differential corrections strongly suggest this as the origin of the deviation.