Crystal Structure Prediction and Lattice Dynamical Calculations for the Rare Platinum-Group Mineral Zaccariniite (RhNiAs)

: The crystal structures of newly found minerals are routinely determined using single-crystal techniques. However, many rare minerals usually form micrometer-sized aggregates that are difficult to study with conventional structural methods. This is the case for numerous platinum-group minerals (PGMs) such as, for instance, zaccariniite (RhNiAs), the crystal structure of which was first obtained by studying synthetic samples. The aim of the present work is to explore the usefulness of USPEX, a powerful crystal structure prediction method, as an alternative means of determining the crystal structure of minerals such as zaccariniite, with a relatively simple crystal structure and chemical formula. We show that fixed composition USPEX searches with a variable number of formula units, using the ideal formula of the mineral as the only starting point, successfully predict the tetragonal structure of a mineral. Density functional theory (DFT) calculations can then be performed in order to more tightly relax the structure of the mineral and calculate different fundamental properties, such as the frequency of zone-center Raman-active phonons, or even their pressure behavior. These theoretical data can be subsequently compared to experimental results, which, in the case of newly found minerals, would allow one to confirm the correctness of the crystal structure predicted by the USPEX code.


Introduction
Platinum-group minerals (PGMs) attract a great deal of research attention, as this group of interesting minerals comprises the following scarce, but technologically important, platinum-group elements (PGEs): Pt, Os, Ir, Ru, Rh, and Pd. The identification, study, and characterization of new PGMs is rarely straightforward because they are typically found as small grains, with mean diameters much smaller than 100 µm [1]. While microprobe analyses readily allow one to detect PGEs in mineral assemblages, it is not usually possible to employ conventional techniques, such as, single-crystal X-ray diffraction (XRD), to obtain structural information about the identified compounds.
Alternative methods are thus required in order to obtain crystallographic data from PGMs. For instance, it is possible to combine the study of a new PGM specimen with that of its synthetic analogue, thus allowing one to infer the crystal structure of the natural counterpart [2,3]. More involved techniques based on electron diffraction can also be employed to resolve the crystal structure of PGMs.
Zaccariniite (RhNiAs) provides a good example of the complexity associated with the structural characterization of PGMs. This mineral was first reported forming a PGM aggregate of ~40 µm, intergrown with other PGE-containing compounds [4]. However, it was only approved by the Commission on New Minerals, Nomenclature and Classification (CNMNC) of the International Mineralogy Association (IMA) after the structural identity between the mineral and synthetic RhNiAs was established [3]. For this purpose, the crystal structure of RhNiAs was refined from powder XRD data obtained from synthetic samples, whereas Raman scattering measurements were used to show the equivalence between the naturally occurring mineral and the synthetic analogue [3]. More recently, transition electron microscopy (TEM), with variable angle precession electron diffraction (PED) techniques, has been used to independently resolve the crystal structure of natural zaccariniite [5,6].
Crystal structure prediction methods provide an alternative means to infer the structural properties of materials. In particular, codes such as USPEX [7][8][9], based on evolutionary algorithms that provide a highly efficient means of predicting stable structures solely from their chemical composition, are capable of successfully predicting the existence of undiscovered phases. This has been particularly relevant over the last few years, in the case of materials exhibiting remarkable properties under high-pressure conditions, such as ultra-hardness [10,11] or superconductivity [12,13]. The discovery of novel, unexpected phases with unusual stoichiometries at high pressure [14][15][16], some of which may be relevant to understand planetary interiors [17], demonstrates the pivotal role of USPEX in high-pressure crystallography.
In the present work, we explore the usefulness of evolutionary structure prediction methods to resolve the crystal structure of PGMs. For this purpose, zaccariniite emerges as a particularly interesting case, due to the knowledge gained in previous works regarding its structural properties [3,5,6], and also due to its relatively simple structure and stoichiometry, thus involving relatively low computational costs. The aim of the present work is to demonstrate that crystal structure prediction algorithms provide a complementary means to experimental methods for this type of problem. In the first part of this work, USPEX is employed to successfully generate the crystal structure of zaccariniite. Subsequent density functional theory (DFT) calculations allow one to predict important properties of the mineral, such as zone-center phonon frequencies, which can then be used in combination with experimental techniques, such as Raman scattering, in order to confirm the crystal structure predicted by USPEX.

Methods
The USPEX code, based on the evolutionary algorithm developed by Oganov and coworkers [7][8][9], was employed in order to determine the crystal structure of zaccariniite. For this purpose, fixed composition searches with a variable number of formula units were carried out. In this case, the formula unit, RhNiAs, corresponded to the ideal formula of the mineral obtained from microprobe analyses [3]. Calculations were performed starting from an initial population of 20 randomly generated structures. No seeds or antiseeds were considered. Subsequent generations were produced by applying the following variation operators (in parentheses are the typical values used in this work to produce newly generated structures): heredity (40-50%), soft mutation (5-10%), lattice mutation (5-10%), and permutation (5-10%). The rest of the newly generated lattices (25-40%) were randomly produced, half of which were produced with the random symmetric structure generator and the other half with the topological random generator. Trial structures generated by USPEX with 1 ≤ Z ≤ 4 formula units were fully relaxed in a three-step process by means of first-principles DFT calculations, using the Quantum Espresso open-source distribution [18][19][20]. The first and last steps only involved relaxation of atomic positions, while the intermediate step involved variable-cell relaxation. As is usual in USPEX calculations, convergence criteria and precision settings were increased after each relaxation step. The reciprocal space resolutions for k-point generation were set to 0.16, 0.12 and 0.08 2π/Å for the three calculation steps, respectively. All DFT calculations were performed using the plane-wave pseudopotential method, with ultrasoft pseudopotentials. For the calculation of the exchange and correlation functionals, generalized gradient approximation (GGA) combined with Perdew-Burke-Ernzerhof (PBE) parametrization [21] was used. No van der Waals (vdW) effects were included in the crystal structure prediction calculations. Due to the metallic character of RhNiAs predicted by the GGA calculations, ordinary Gaussian spreading was used to perform Brillouin-zone integrations. It should be noted, however, that preliminary hybrid functional calculations suggested that RhNiAs might be a narrow-band semiconductor, with a band gap of ~0.15 eV. More work should be performed in order to properly characterize the electronic band structure of this material.
Once the best structure for RhNiAs was found using USPEX, additional DFT calculations using the Quantum Espresso package were carried out in order to perform a tighter relaxation of the structure and determine different properties of the material, both in room conditions and, for completeness, as a function of pressure. As in the case of the USPEX predictions, these calculations were performed using the plane-wave pseudopotential method, with ultrasoft pseudopotentials. For the sake of comparison, the following two different parametrizations were considered for the exchange and correlation functionals: the local density approximation (LDA) and GGA-PBE. In both cases, the following two different sets of calculations were performed: (i) a first set in which long-range vdW interactions were neglected, and (ii) a second set including vdW interactions using Grimme's D2 correction [22]. In all cases, the crystal structure of RhNiAs was fully relaxed as a function of pressure until interatomic forces were lower than 10 −6 Ry/Bohr, with a plane-wave basis cut-off of 100 Ry. Monkhorst-Pack k-point grids were set to 8 × 8 × 6.
Lattice-dynamical calculations were performed using the density-functional perturbation theory (DFPT) approach implemented in the Quantum Espresso package [18,19], which allowed us to obtain Γ-point phonon frequencies at zero pressure from single qvector calculations on the primitive cell, and also to obtain their corresponding pressure coefficients by studying the pressure dependence of the phonon modes. The acoustic sum rule was applied to obtain the correct acoustic phonon frequencies at zero wave vector (q = 0). Additional calculations were carried out using the finite displacement (FD) approach, using Phonopy software [23], in order to obtain the phonon dispersion and the one-phonon density of states (1-PDOS). For this purpose, atomic forces arising from different FDs were obtained with Quantum Espresso DFT calculations using 2 × 2 × 2 supercells.

Crystal Structure Determination
Zaccariniite provides a good example for illustrating the usefulness of USPEX in order to determine the crystal structures of newly found PGMs, many of which exhibit relatively simple chemical formulas and/or crystal lattices (i.e., eminently inorganic with no H2O). The present work can be envisaged as a practical example of how one would proceed in order to determine and experimentally confirm the crystal structure of new PGMs or similar minerals, assuming that only the chemical formula is known (for instance, from microprobe measurements, as in the case of zaccariniite [3]).
Indeed, USPEX correctly predicts the tetragonal, ambient-pressure crystal structure of zaccariniite (space group P4/nmm, No. 129, Z = 2), after just a few computed generations. Figure 1 shows a selected example of a typical USPEX run, in which the correct structure was found after 16 generations. The figure shows the enthalpy of the best individual, per formula unit, after each generation; the minimum of the enthalpy corresponds to the tetragonal structure of zaccariniite. Repeated calculations with slightly modified variation operators provided, in all cases, the correct result. Remarkably, in one case, the correct structure was achieved after only four generations. The projection onto the c-a plane of the resulting structure, with double layers of Rh and As atoms, can be observed in Figure 1.
USPEX runs were also performed at 10 GPa and 20 GPa to assess the stability of the P4/nmm structure under high-pressure conditions. These runs provided the same structure as for the zero-pressure case, thus suggesting a large stability range for the ambientpressure structure of RhNiAs. A tighter relaxation of the zero-pressure P4/nmm structure of RhNiAs was subsequently performed with a full variable-cell relaxation, setting the convergence criterion to achieve interatomic forces lower than 10 −6 Ry/Bohr. However, the presence of double Rh and As layers in the structure of the mineral suggests that vdW interactions could be relevant. Therefore, in order to evaluate the importance of dispersion effects in RhNiAs, variable-cell relaxations were performed by including, or neglecting, vdW forces (using Grimme's D2 correction [22]). For completeness, calculations were performed with LDA or PBE GGA functionals, and as a function of hydrostatic pressure, which allowed further comparison of the resulting compression behavior. The results of these calculations are shown in Figure 2 and summarized in Table 1. Bulk modulus values were obtained by fitting the calculated data to a third-order Birch-Murnaghan equation of state.   It can be observed that both the type of functional and the inclusion of vdW effects had a sizable bearing on the calculated equation of state (see Figure 2). As expected, PBE calculations yielded slightly larger unit-cell volumes, while the inclusion of vdW in LDA produced somewhat disrupted results. This can be attributed to the fact that, in LDA, there is a compensating effect between the overestimated covalent part of the interlayer bonds and the neglected vdW interactions [24]. In the case of LDA-vdW, such compensating effects no longer take place. On the other hand, the zero-pressure volume and lattice parameters of PBE approach those of LDA upon the inclusion of vdW effects.
As no experimental data on the compression properties of natural or synthetic RhNiAs are available, it is not possible to conclude, with the present results alone, which functional (PBE-vdW vs. LDA) better predicts the zero-pressure lattice parameters and the P-V equation of state of this compound. However, it is well known that LDA fails to reproduce the pressure behavior of materials whose fundamental properties are strongly determined by vdW interactions. This is the case, for instance, for the vibrational properties of layered transition metal dichalcogenides (TMDCs), such as MoS2, HfS2, and ReS2 [25,26]. For these materials, it is found that PBE calculations do properly describe their high-pressure lattice dynamical properties, provided that vdW corrections are considered. In the case of RhNiAs, as will be discussed below, only PBE-vdW provided reliable results for the phonon modes, even at zero pressure. Hence, we have compiled, in Table  2, the best calculated results for RhNiAs, obtained with PBE parametrization and vdW corrections. Note the excellent agreement between the theoretical data (lattice parameters and atomic positions) and the experimental values given in Ref. [3]. Table 2. Calculated structural data for zaccariniite (RhNiAs) obtained from DFT calculations within the generalized gradient approximation (GGA), with Perdew-Burke-Ernzerhof (PBE) parametrization, including vdW corrections (Grimme's D2 method).

Lattice Dynamical Calculations
At this point, it is important to recall that in the present work we are concerned with the problem of determining, from ab initio calculations alone, the crystal structure of zaccariniite. Therefore, the following natural step after having constrained the structural data for the mineral (Table 2) would be a comparison with experimental results. Assuming that one has no access to involved analytical techniques, such as TEM, or to the synthetic material, it is clear that Raman scattering is a very appropriate tool to confirm the correctness of the crystal structure data obtained so far. Indeed, micro-Raman measurements provide a fast, non-destructive means of characterizing materials at the micron scale. Here, we present the results of DFT lattice-dynamical calculations, and compare them with the Raman results reported in Ref. [3].
The crystal structure predicted for zaccariniite, summarized in Table 2, has six atoms in the unit cell and, therefore, 18 normal modes of vibration, 15 of which correspond to optic modes (Γoptic = 2A1g + 2A2u + B1g + 2Eu + 3Eg). Among these, the Raman-active modes are ΓRaman = 2A1g + B1g + 3Eg, where only Eg is doubly degenerate. For this structure, with a center of inversion symmetry, odd modes (u, ungerade) are infrared (IR) active, while even modes (g, gerade) are Raman active. Information about the Raman or IR activity of the optical modes of RhNiAs can be found in Table 3. Table 3. Theoretical zero-pressure frequencies (ω0) for the zone-center phonons of RhNiAs and their corresponding linear and quadratic pressure coefficients, a and b, as obtained by fitting the PBE-vdW DFPT calculations of Figure 3 as a function of pressure (P) to the equation ω(P) = ω0 + aP + bP 2 .

Symmetry (Activity)
Frequency  Our results of the lattice-dynamical calculations of the zone-center optical modes of RhNiAs, performed within the framework of DFPT, using PBE parametrization and including vdW corrections using Grimme's D2 method, are shown in the second column of Table 3. Unfortunately, no detailed information about the experimental Raman peaks of natural and synthetic RhNiAs is given in Ref. [3]. However, a direct comparison of the spectra published in that work with the theoretical data of Table 3 is still possible, and allows one to find very good agreement between both works. First, it can be observed that the Raman spectra reported in Ref. [3] exhibit a sharp peak just below 300 cm −1 that can be readily assigned to first-order Raman scattering by A1g and Eg modes, which are expected to overlap around 295 cm −1 , according to the present calculations (Table 3). On the other hand, the four features that appear in the Raman spectra of RhNiAs within the 100-200 cm −1 range can be attributed to Eg, B1g, Eg, and A1g modes. According to these attributions, the sharpest peak that dominates the Raman spectrum of RhNiAs would correspond to an A1g mode located around 180 cm −1 .
It is important to note that imaginary frequencies for the low-frequency B1g mode were obtained when the DFPT phonon calculations were performed without considering vdW interactions, regardless of the functional parametrization (LDA or PBE). This is due to the fact that, as can be observed after diagonalization of the dynamical matrix for the zone-center phonons, the low-frequency B1g mode involves atomic displacements along the x-y plane, with the Ni and As layers oscillating in phase, and 180° out of phase, relative to the Rh layers. As occurs in layered materials, such as TMDCs, these types of vibrations are strongly dependent on vdW forces [25,26], and, as a consequence of this, no reliable results were obtained for this mode when the lattice dynamical calculations were performed neglecting vdW effects. In other words, it is essential to include vdW corrections to calculate the phonon properties of layered materials, such as RhNiAs. In Figure S1 of the Supplementary Material, the calculated ambient-pressure phonon band structure and the corresponding phonon density of states, obtained using the FD approach, are shown. As can be observed in the figure, no imaginary (negative) frequencies were found in any of the high-symmetry directions considered. This is in contrast to what occurred when vdW corrections were not taken into account (not shown). Similar results were obtained at 10 and 20 GPa, which further confirms the stability of the P4/nmm structure of RhNiAs under high pressure, as predicted by USPEX.
Finally, it should be noted that the experimental Raman spectrum of RhNiAs might also exhibit second-order features and/or disorder-induced modes that could be stronger than some of the first-order Raman-active modes. In order to unambiguously assign the different features in the Raman spectrum of zaccariniite, a detailed study combining theoretical data and Raman measurements with different scattering geometries would be required. However, such experiments would be difficult due to the lack of information regarding the orientation of the available mineral aggregates.

Zone-Center Phonons of RhNiAs as a Function of Pressure
An alternative means of identifying the nature of zone-center phonon modes in materials has been recently outlined by Manjón et al. [27]. These authors have stressed that, as in the cases of many rare minerals, such as zaccariniite, a large number of materials cannot be obtained in a single-crystal form, and, therefore, it is not possible to correctly identify their Raman spectra by polarization-dependent Raman measurements. In the case of understudied materials, conventional zero-pressure DFT-based lattice dynamical calculations do not suffice for this purpose, due to large associated uncertainties in the calculated frequencies, which are found to strongly depend on functional parametrization. Manjón and co-workers have shown, with numerous examples, that such identification is possible by combining experimental and theoretical results of high-pressure Raman scattering measurements and lattice-dynamical calculations. These authors demonstrate that careful monitoring of the experimental pressure dependence of the Raman modes, with the support of lattice-dynamical calculations, allows one to correctly assign the symmetry of most Raman-active features. The important point is that such an assignment can then be extrapolated to the conventional ambient-pressure spectrum, thus providing unambiguous identification of the experimental Raman spectrum of a material.
In the case of rare minerals, such as zaccariniite, once the crystal structure has been inferred from the USPEX computations, it is clear that the joint theoretical and experimental method proposed by Manjón et al. [27] would allow one to further support the results of the crystal structure predictions. With regards to this, it should be mentioned that high-pressure Raman measurements involve microsamples similar to those found in PGMs, which should be loaded into a diamond anvil cell (DAC) and measured using a micro-Raman setup. This type of experiment would be totally feasible in the particular case of zaccariniite.
Hence, for the sake of completeness, we show, in Figure 3, the pressure dependence of DFPT-calculated, zone-center phonon modes of RhNiAs up to 25 GPa, which is a typical pressure value achieved with DACs. Note also that, as previously mentioned, the USPEX calculations predicted that the P4/nmm structure was stable, at least up to 20 GPa. Linear and quadratic pressure coefficients for all the zone-center modes (Raman and IR active) are shown in the last two columns of Table 3. It can be observed, in Figure 3, that the lowfrequency Eg mode exhibits a small pressure coefficient, which could be easily observed in a typical high-pressure Raman experiment. On the other hand, according to these results, progressive overlapping of the Eg and A1g modes at intermediate frequencies would be expected. Finally, the high-frequency Eg and A1g peaks would be found to progressively split with increasing pressure, giving rise to peak broadening at low pressures, and, presumably, yielding two separate peaks at higher pressures. All these observations would be fingerprints for the RhNiAs compound, further confirming the correctness of the predicted crystal structure.

Conclusions
The crystal structure of eminently inorganic, relatively simple minerals, such as zaccariniite (RhNiAs), which are found in nature as micro-sized aggregates that are difficult to analyze using standard single-crystal techniques, can be reliably determined with crystal structure prediction methods, such as USPEX. For this purpose, only information about the chemical composition of the mineral is required. Subsequent DFT calculations can be performed in order to accurately determine different fundamental properties of the mineral, such as, for instance, lattice parameters or Raman-active phonon frequencies. The latter can then be used to confirm the correctness of the predicted structure by means of conventional Raman scattering experiments, which are typically less involved than other analytical tools. In addition, following the recent Gedankenexperimente of Manjón et al. [27], high-pressure Raman measurements might be used to properly assign the phonon modes of the mineral, just by comparing the experimental results with theoretical phonon frequencies as a function of pressure. This method, applied by Manjón and collaborators over the last few years to many different material systems [27], might provide further confirmation regarding the correctness of the crystal structure predictions of USPEX.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/article/10.3390/min12010098/s1, Figure S1: Calculated phonon band structure of RhNiAs using the finite displacement (FD) method with GGA-PBE parametrization and including van der Waals corrections (Grimme's D2). The right panel shows the corresponding phonon density of states (PDOS).
Funding: This research received no external funding. Data Availability Statement: Not applicable.