GIPAW Pseudopotentials of d Elements for Solid-State NMR

Computational methods are increasingly used to support interpreting, assigning and predicting the solid-state nuclear resonance magnetic spectra of materials. Currently, density functional theory is seen to achieve a good balance between efficiency and accuracy in solid-state chemistry. To be specific, density functional theory allows the assignment of signals in nuclear resonance magnetic spectra to specific sites and can help identify overlapped or missing signals from experimental nuclear resonance magnetic spectra. To avoid the difficulties correlated to all-electron calculations, a gauge including the projected augmented wave method was introduced to calculate nuclear resonance magnetic parameters with great success in organic crystals in the last decades. Thus, we developed a gauge including projected augmented pseudopotentials of 21 d elements and tested them on, respectively, oxides or nitrides (semiconductors), calculating chemical shift and quadrupolar coupling constant. This work can be considered the first step to improving the ab initio prediction of nuclear magnetic resonance parameters, and leaves open the possibility for inorganic compounds to constitute an alternative standard compound, with respect to tetramethylsilane, to calculate the chemical shift. Furthermore, this work represents the possibility to obtain results from first-principles calculations, to train a machine-learning model to solve or refine structures using predicted nuclear magnetic resonance spectra.


Introduction
The application of NMR spectroscopy to rigid or semi-rigid solid samples allows the study a plan of systems as bio-molecules with high molecular weight, polymers, perovskites (e.g., solar cells absorbers) and cements in chemistry and chemical sciences. NMR is the oscillatory response of nuclei with non-zero nuclear spins (total angular momentum) immersed in an external field (B 0 ). The presence of B 0 removes the degeneracy of nuclear spins, leading to the energy difference: where γ is the gyromagnetic ratio; and σ is the chemical shielding around a nucleus, which is a characteristic of a specific isotope. Thus, the chemical structure can be revealed by NMR frequencies that are significantly affected by γ and σ. NMR frequencies are reported as a chemical shift (δ), which is the fractional difference between the frequency of a particular nucleus and a standard compound such as tetramethylsilane (TMS). If NMR seems to be, abstractly, the best way to determine chemical structure, the NMR frequencies are strictly anisotropic, being dependent on the relative orientation between B 0 and a sample, with the consequent generation of internuclear couplings and quadrupolar couplings. Actually, the quadrupolar coupling constant (C Q ) can be estimated as where eQ is the electric quadruple moment, V zz is the potential of electric field B 0 along z-axes, and h is the Planck constant. Thus, these anisotropic interactions need to be partially averaged through the molecular rotations, and measurement of motionally averaged NMR spectra and induced nuclear spin relaxation reveals the geometries and rates of motion. Furthermore, the nuclear magnetic resonance (NMR) signal is orders of magnitude lower in frequency than the microwave, infrared and ultraviolet frequencies employed in rotational, vibrational and electronic spectroscopes. This is due to the low population difference between nuclei with removed degeneracy and those in the ground state, causing low signal-to-noise ratios along the spectra. In solid-state NMR, low signal-to-noise ratio is accentuated by the presence of acoustic phonon deformation potential (ADP) scattering and optical phonon branches. They can be responsible for electron-phonon coupling, which can, alternatively, affect the population difference. Thus, the development of NMR is focused on the increasing of experimental sensitivity. This can be carried out through increasing the intensity of an applied magnetic field, with the consequent increase in ∆E, but it is very expensive. Or it can be carried out by recording NMR spectra in the domain following a radio-frequency pulse and obtaining the spectrum by Fourier transformation rather than by sweeping the frequency and measuring absorption or emission in classical spectra. Fourier transformation NMR spectroscopy increases by one order of magnitude and opens the door to multidimensional NMR spectroscopy. In spite of its sensitivity, the interpretation of NMR spectra can be less intuitive than microscopy or diffraction data, because structural information is encoded in frequency spectra rather than spatial density maps. The frequency peaks need to be assigned to individual atoms, which can be a significant challenge. However, the multitude of peaks in NMR spectra represent an exquisite chemical fingerprint of molecules, thus making NMR spectroscopy of great use to chemists. Thus, computational methods are increasingly used to support interpreting, assigning and predicting the solid-state NMR spectra of materials. Furthermore, density functional theory (DFT) gave excellent results with gauges including atomic orbitals [1] (GIAO) for soft-matter NMR; this approach cannot be applied to solid-state NMR, because all-electron calculations are not performable due to required computational resources and the necessity to preserve translation symmetry in solids. Thus, in the framework of planewaves DFT, a gauge including the projected augmented wave [2,3] (GIPAW) method was introduced to calculate nuclear resonance magnetic (NMR) parameters in solids, avoiding all-electron calculations. In the GIPAW approach [4,5], a uniform magnetic field is applied using boundary conditions, a periodic magnetic field with a finite wavelength r G is the gauge origin and is subsequently extrapolated in the limit r G → 0 to compute the chemical shielding. This formalism was seen to manage the numerical instabilities associated with the summation of two divergent terms and with the generalized gradient approximation exchange-correlation functional (the method of choice for condensed-matter simulations) to perform accurate results [6][7][8][9][10][11][12][13][14][15]. In this work we have developed GIPAW pseudopotentials for the elements of first, second and third rows of d elements excluding La, which is considered as a part of the Lanthanides. These pseudopotentials were tested on the oxides or nitrides optimizing the crystal structures and, subsequently, we calculated the NMR parameters. The developing of these pseudopotentials was carried out to increase the number of compounds for which NMR parameters can be calculated.

Theoretical Background
In plane-waves DFT the all-electron potential of an atom is substituted by a mathematical object, the so-called pseudopotential. The all-electron wave functions are substituted by pseudo-wavefunctions that eliminate the core states and describe only the chosen valence pseudo-wavefunctions. This limits the generation of a pseudopotential with a specific configuration. Actually, there are different types of pseudopotential: norm-conserving [16], ultrasoft [17] and projected augmented wave (PAW) [18]. The last introduces a linear operator that converts the pseudo-wavefunctions |φ R,n ⟩ to all-electron wavefunctions |ϕ R,n ⟩. In addition, ⟨p R,n | is a set of projectors such that ⟨p R,n |φ R,n ⟩ = δ R,R' δ n,m . Thus, it is possible to introduce in the Blöchl formalism [19] a field dependent transformation operator Γ B that restores the translational invariance.
Such re-formalism is called gauge including projected augmented wave (GIPAW) [2,3] and satisfies the translation relation |Ψ⟩ = Γ B |Ψ⟩. The GIPAWpseudo-operator O = Γ + B OΓ B , corresponding to a local or a semilocal operator O, is given by The GIPAW [2,3] method allows the calculation of an induced magnetic field at the nucleus position B ind (r) to the applied external magnetic field B 0 (r), according to: where σ is the magnetic shielding tensor that proceeds through the calculation of the first-order induced current density j (1) (r), which reads: The summation runs over the occupied states and n 0 (r) is the unperturbed charge density. In addition, ϕ 0 are the unperturbed Kohn-Sham orbitals and their first-order ϕ 1 counterpart, perturbed due to the external magnetic field . In Equation (7) is given the decomposition of the induced current into the so-called paramagnetic term j (1) p (r), which involves the first-order perturbed orbitals, and the diamagnetic term j (1) d (r), which depends on the unperturbed charge density. A(r) is the vector potential connected to B 0 through where r G is the so-called gauge origin. B ind (r) is finally obtained from the Biot-Savart law: Actually, there are connections between our GIPAW [2,3] approach and the GIAO [1] method widely used in the quantum-chemical community for the all-electron calculation NMR of molecules. However, it should be recognized that, in GIPAW [2,3], the phase required to maintain the translational invariance is carried by the operators, whereas in the GIAO [1] approach, the field-dependent phase is attached to the basis functions and to the occupied electronic orbitals, respectively.

Method and Computational Details
GIPAW pseudopotentials were developed for 21 d, elements excluding La which is considered as a part of Lanthanides, using Quantum Espresso version 6.6 [20,21] and they are written in the universal pseudopotential format (UPF) version 2 (pseudopotentials will be available to everybody as UPF2 format, see available dataset). GIPAW were developed solving the scalar relativistic wave equation (Koelling-Harmon-like equation) [22] with Rappe-Rabe-Efthimios-Kaxiras-Joannopoulos [23] (RRKJ) form of pseudo-wavefunctions modeled by double projectors and semi-core states. A non-linear core correction (NLCC) [24] was employed for all developed GIPAWs. The NLCC allows to avoid the necessity to separate spin-up and spin-down ionic pesudopotentials, treating explicitly the nonlinear exchange and correlation interaction between the core and the valence charge densities. In particular, the spin-polarized configurations are well-described with a single potential. The analysis of logarithmic derivatives, i.e., derivatives of an l-state d(log(Ψ l (E))) dE computed for the exact atomic problem and with the GIPAW dataset, was computed to verify the presence of highly-localized negative energy ghosts that could affect the quality of pseudopotential. No element presented highly localized negative-energy ghosts and only highly localized positive-energy ghosts are seen in the d orbital of Nb Os and Ta, but they are located too high in energy (i.e., 1 Ry for Nb and Os, 4 Ry for Ta) to affect the quality of pseudopotential also, cases with the promotion of one electron in the d orbital will be described. All the logarithmic derivatives and NLCC for each pseudopotential are shown in the available data. Considered oxides and nitrides of d elements of first-, second-and third-period, which are semiconductors, were fully optimized in Quantum Espresso version 6.6 [20,21] with previously developed GIPAW [2,3] of d elements for PBE exchange-correlation density functional [25]. The crystal structures were taken from Material Project [26] [25]. The geometry optimization relies on Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [27][28][29][30], with force tolerance for the maximum net force on atoms fixed at 10 −6 Ry/Angstrom and the total energy at 10 −8 Ry. The kinetic energy plane-wave energy cutoff of 100 Ry, and the Gaussian smearing [31]

Results and Discussion
In metal transitions, the so-called semi-core states overlap with the d-valence states. Thus, we introduced the semi-core states in valence states with a small core region radius, improving the accuracy. For the 3d transition metals, the semi-states 3s 3p 3d were chosen as part of the valence partition. The GIPAWs for the 4d transition metals all contain the 4s 4p 4d semi-states in the valence. For the 5d transition metals, the 5s 5p 5d semi-states are contained in the valence of GIPAWs. In 5d transition metals, 4f states were frozen due to the complexity, to be rightly described lying at the same energy range as the 5s and 5p states. It is noteworthy that, for such elements, the ground-state properties can be described well enough with 4f frozen. Meanwhile, for optical properties and GW approximation [32,33], this may not be the case. This is the case even for elements such as Au, where the 4f electrons lie about 3 Ha below the Fermi level. Instead, for GIPAW for 3d and 4d elements, GW approximation [32,33] can also be performed with such pseudopotentials. The ground-state properties of 5d elements are not affected by 4f orbitals because they are contracted close to the core due to relativistic effects and consequent jj-spin orbit coupling. This means that, firstly, for each orbital its orbital angular momentum is coupled with the magnetic angular momentum of the electron located in such an orbital, obtaining an angular momentum j and, subsequently, all these j are coupled to make the total angular momentum of a heavy atom. We fully relaxed the crystal structures of the oxides and nitrides of the 21 d elements of which we developed GIPAW pseudopotentials using the PBE DFT functional. The PBE DFT functional is known in the literature to overestimate the unit volume by close to 4-4.5%, but the agreement between the calculated and experimental NMR parameters is generally found to be significantly improved after the DFT optimization of the structure geometry [34]. This motivated our choice to optimize the structure before calculating the NMR parameters. We are aware that all oxides and nitrides of d elements are considered strong correlated systems. This means that heuristic DFT alone seems to not be enough to describe correctly the electronic structure spreading the d electrons within the unit cell. Usually, DFT+U [35] typically works well for strongly correlated systems localizing the d electrons on metals, but previous works showed that the GIPAW formalism feasibly describes the transition metals with an heuristic DFT [36][37][38][39]. Thus, we have used a heuristic DFT to fully relax the chosen systems. The fully optimized lattice parameters of such systems agree with those that come from experiments and with the computed structures showed in Material Project [26], see Figure 1 and Tables S1 and S2 in Supporting Information. Such computed structures, reported in Material Project [26], were actually optimized using projector augmented-wave (PAW) pseudopotentials [19] in VASP code [40][41][42][43]. Confident of the quality of the obtained lattice parameters of the oxides and nitrides of the d elements, we were able to use the GIPAW approach to calculate the total NMR chemical shift (σ) by adopting the Simpson [44] convention for anisotropy and asymmetry and the average value of the d elements within the oxide or nitrides are reported in Table 1. We have reported, in Table 1, the average values because, in some structures, the atoms that are in non-equivalent for symmetry positions have different values and in Supporting Information from Tables S3-S23 we reported all atomic positions with corresponding (σ). It is noteworthy that we have obtained negative (σ) in some cases and this is due to the fact that some d elements are more shielded with respect to the chosen reference (i.e., TMS) compound. This make us to consider if the evolution of NMR should be focused to find an inorganic compound to set σ. Indeed, here we considered nitrides or oxides that are semiconductors. Thus, our results are not affected by Knight shift, as in metals, due to the high population of d electrons at the Fermi level. We also calculated C Q and its average values for the d elements within oxides or nitrides, see Table 1, while the value for each atom is reported in Supporting Information. Of note is the impossibility of using predicted NMR spectra from first-principles calculations to solve or refine structures due to the time and cost of the calculation, which poses challenges to a real-time automated solution. To address this problem, machine-learning approaches were introduced to calculate chemical shifts in molecular solids, which reduces computational cost by orders of magnitude while maintaining the accuracy of DFT [45]. All these machine-learning models are trained on first-principles calculations, making the latter of fundamental importance to improve the accuracy of machine-learning models.

Conclusions
In this work, we developed GIPAW pseudopotentials for 21 d elements and tested them on their oxides or nitrides. The developed pseudopotentials present semi-core states to increase their flexibility, to be employed in several compounds with feasible approximation, and are free of ghosts states, making them conceptually right. The obtained σ and C Q for the d elements seem to be coherent with expected values for such elements, respectively, in their oxides or nitrides. The GIPAW pseudopotentials increase the number of compounds for which it will be possible to calculate the NMR parameters with first-principle calculations and subsequently use them to develop machine learning models for real-time refining structure from NMR spectra.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ma15093347/s1, Table S1: Table of lattice parameters of 10 fully-relaxed unit cells with developed GIPAW pseudopotentials and those optimized with PAW pseudopotentials in VASP code and the experimental ones; Table S2: Table of lattice parameters of 11 fully-relaxed unit cells with developed GIPAW pseudopotentials and those optimized with PAW pseudopotentials in VASP code and the experimental ones; Table S3: Atomic crystallographic positions with vectors lattice in Angstrom of AgN 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S4: Atomic crystallographic positions with vectors lattice in Angstrom of Au 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S5: Atomic crystallographic positions with vectors lattice in Angstrom of CdO 2 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S6: Atomic crystallographic positions with vectors lattice in Angstrom of Cr 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S7: Atomic crystallographic positions with vectors lattice in Angstrom of HgO with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S8: Atomic crystallographic positions with vectors lattice in Angstrom of IrN 2 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S9: Atomic crystallographic positions with vectors lattice in Angstrom of Lu 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S10: Atomic crystallographic positions with vectors lattice in Angstrom of MoO 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S11: Atomic crystallographic positions with vectors lattice in Angstrom of Nb 2 O 5 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S12: Atomic crystallographic positions with vectors lattice in Angstrom of OsO 4 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S13: Atomic crystallographic positions with vectors lattice in Angstrom of PdN 2 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S14: Atomic crystallographic positions with vectors lattice in Angstrom of PtO 2 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S15: Atomic crystallographic positions with vectors lattice in Angstrom of Re 2 O 7 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S16: Atomic crystallographic positions with vectors lattice in Angstrom of Rh 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S17: Atomic crystallographic positions with vectors lattice in Angstrom of RuO 4 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S18: Atomic crystallographic positions with vectors lattice in Angstrom of Sc 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S19: Atomic crystallographic positions with vectors lattice in Angstrom of Ta 2 O 5 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S20: Atomic crystallographic positions with vectors lattice in Angstrom of Tc 2 O 7 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S21: Atomic crystallographic positions with vectors lattice in Angstrom of V 2 O 5 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S22: Atomic crystallographic positions with vectors lattice in Angstrom of WO 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element; Table S23: Atomic crystallographic positions with vectors lattice in Angstrom of Y 2 O 3 with total chemical shift σ in ppm and quadrupolar coupling constant C Q in MHz for each element.