Probing the Local Atomic Structure of In and Cu in Sphalerite by XAS Spectroscopy Enhanced by Reverse Monte Carlo Algorithm

: The distortion of atomic structure around In and Cu dopants in sphalerite ZnS was explored by extended X-ray absorption ﬁne structure (EXAFS) spectroscopy enhanced by reverse Monte Carlo (RMC) simulation (RMC-EXAFS method). These data were complemented with quantum chemical Density Functional Theory (DFT) calculations and theoretical modeling of X-ray absorption near edge spectroscopy (XANES) spectra. The RMC-EXAFS method showed that in the absence of Cu, the In-bearing solid solution is formed via the charge compensation scheme 3Zn 2 + ↔ 2In 3 + + (cid:3) , where (cid:3) is a Zn vacancy. The coordination spheres of In remain undistorted. Formation of the solid solution in the case of (In, Cu)-bearing sphalerites follows the charge compensation scheme 2Zn 2 + ↔ Cu + + In 3 + . In the solid solution, splitting of the interatomic distances in the 2nd and 3rd coordination spheres of In and Cu is observed. The dopants’ local atomic structure is slightly distorted around In but highly distorted around Cu. The DFT calculations showed that the geometries with close arrangement (clustering) of the impurities—In and Cu atoms, or the In atom and a vacancy—are energetically more favorable than the random distribution of the defects. However, as no heavy In atoms were detected in the 2nd shell of Cu by means of EXAFS, and the 2nd shell of In was only slightly distorted, we conclude that the defects are distributed randomly (or at least, not close to each other). The disagreement of the RMC-EXAFS ﬁttings with the results of the DFT calculations, according to which the closest arrangement of dopants is the most stable conﬁguration, can be explained by the presence of other defects of the sphalerite crystal lattice, which were not considered in the DFT calculations.


Introduction
In our previous works [1,2], we studied the local atomic environment and electronic state of dopants (In, Cu, Ag, Au, Cd) in natural and synthetic ZnS crystals by means of X-ray absorption spectroscopy (XAS). Analysis of the X-ray absorption near edge structure (XANES)/Extended X-ray absorption fine structure (EXAFS) spectra revealed that in substitutes for Zn in the structure of The methods of synthesis of the doped sphalerites and reference materials, analytical methods, and the methods of XAS measurements and preliminary spectra treatment (normalization, background removing) are described in detail in our previous publication (Section 2 in Trofimov et al. [2], this issue). The sample IDs are identical throughout both the contributions.

Sample Preparation and Characterization
In the present work, we investigated three samples of synthetic sphalerite. The compositions of the samples are given in Table 1 in Trofimov et al. [2]. One sample was doped with one element-In (Sample 3757-0.7 wt.% In), and the other two samples were doped with two elements-In + Cu (Sample 4108-0.1 wt. Cu and 0.09 wt.% In; 4186-2.7 wt.% Cu and 4.9 wt.% In). The reference materials include pure sphalerite ZnS and roquesite CuInS 2 .

System E EF , eV
ZnS + Cu + In (In and Cu atoms are in neighboring sites) 0.18 ZnS + Cu + In (In and Cu atoms are placed far from each other) 0.28 ZnS + 2In, (

X-ray Absorption Spectroscopy (XAS) Measurements
The Zn, In, and Cu K-edge spectra were recorded at the Rossendorf Beamline BM20 of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France. Energy calibration was performed using the K-edge excitation energy of Zn (9659 eV), In (27,940 eV), and Cu (8979 eV) metal foil. The spectra of reference substances were collected in transmission mode, while the spectra of sphalerite samples were recorded in total fluorescence yield mode using a 13-element high-throughput Ge-detector. The detected intensity was normalized to the incident photon flux. Normalization and extraction of the χ(k) signal were performed using the ATHENA program (IFEFFIT software package [3], see details in Trofimov et al. [2]).

Density Functional Theory (DFT) Calculations
In our previous studies [1,2], it was demonstrated by the analysis of the X-ray absorption spectra of In-Cu-bearing sphalerites that the formation of the sphalerite solid solution takes place according to the charge compensation schemes: 3Zn 2+ ↔2In 3+ + and 2Zn 2+ ↔In 3+ + Cu + . In the present study, the formation energy of In-and In-Cu-bearing sphalerite solid solutions was calculated by means of Density Functional Theory (DFT). Besides, inthis way, the distortion of the local atomic geometry around dopants was explored.
The software package QUANTUM ESPRESSO (version 6.1) [9] was used for the calculations. We employed a projector-augmented wave description of the electron-ion interactions [10,11] with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional. The self-consistent field (SCF) method with a 70 Ry kinetic energy cut-off for the plane waves, a 1000 Ry charge density cut-off, and a SCF tolerance >10 −9 Ry was applied in the electronic structure calculations. The optimizations of the crystal structure and supercell parameters were performed using the Broyden-Fletcher-Goldfarb-Shanno algorithm for the atomic coordinates with convergence threshold 10 −3 Ry/a.u. for the forces and 10 −4 Ry for the energy. The relaxation of the atomic positions and cell parameters were applied for a 3 × 3 × 3 supercell, which contained 108 Zn and 108 S atoms, with periodic boundary conditions. In all cases, the large unit cell allowed gamma point approximation to be employed. Our previous investigations [1,12,13] demonstrated the high accuracy of this approach to DFT calculations in the case of sphalerite-type structures. The discrepancy between the calculated interatomic distances and EXAFS results was 0.02, 0.01, and 0.04 Å for the 1st, 2nd, and 3rd coordination shells, respectively.
The formation energies of sphalerite solid solutions (E FE , [12,14]) were obtained using the calculated total energies. We considered the following models: two In atoms substitute for two Zn atoms in ZnS crystal structure; two In atoms substitute for two Zn atoms in ZnS crystal structure with the formation of a vacancy at Zn position; one In and one Cu atom substitute for two Zn atoms in ZnS crystal structure. Accordingly, the E FE values per impurity atom after geometry relaxation were calculated as defined by the following equations: where E is the total energy of the ZnS, ZnS-In, or ZnS-Cu-In supercells; E(ZnS), E(CuInS 2 ), and E(In 2 S 3 ) are the energies calculated by DFT for pure sphalerite, roquesite, and In 2 S 3 crystal structures; E(S) is calculated for the isolated S atom placed in a vacuum box. The structures with lower E EF are more stable.

Reverse Monte Carlo (RMC) EXAFS Spectra Fitting
To explore distortion of the local atomic structure caused by the presence of dopants, we employed a reverse Monte Carlo (RMC) method [15,16]. The evolutionary algorithm (EA) of the RMC method was used for the EXAFS spectra fitting as implemented in EvAX code [8]. An important advantage of the RMC method is the automatic handling of multiple scattering paths (MS). The inclusion of MS paths in the conventional EXAFS fitting procedure is usually associated with a drastic increase in the number of fitted parameters, which often makes the obtained results unreliable.
The RMC fitting aimed to obtain the best agreement between the experimental and the calculated spectra by small variations of the atomic coordinates. The maximum displacement of atoms from their initial positions was fixed at 0.4 Å; this value is enough to describe the thermal movement of atoms and distortion of the sphalerite crystal structure. A comparison of experimental and theoretical spectra was made in the Wavelet transform (WT) space [17,18] (see Trofimov et al. [2] for details of the WT calculation). The experimental EXAFS spectra were fitted in the R-range between 1.3 and 4.5 Å; the k-range was chosen for every spectrum depending on the signal/noise ratio. The E 0 shift was determined for each EXAFS spectrum by fitting the first peak of Fourier transform magnitude.
Calculations of the scattering amplitude and phase functions were performed by the FEFF8 [19] code. We used atomic structures derived from ZnS crystal by substituting the central Zn atom by In or Cu as a starting geometry for EXAFS refinement. Such clusters contain a central absorbing atom (In or Cu) and three coordination shells with 4 S, 12 Zn, and 12 S atoms (see Figure S1). The radii of coordination shells around central Cu or In atoms were corrected to match the results of the conventional EXAFS spectra fitting [2]. To speed up the calculations, in the calculations of doped sphalerites, we used isolated clusters without periodic boundary conditions. EXAFS spectra at the Cu K-edge were fitted using a cluster of ZnS structures, with and without In in the second coordination shell, to examine the possibility of the close arrangement of Cu and In dopants. The In K-edge EXAFS spectra were fitted without Cu in the second coordination shell since Cu has the same scattering amplitudes and phase functions as Zn. As will be shown below, the 5 Å cluster radius is enough to fit the EXAFS spectra in the R range from 1.3 to 4.5 Å. It is important to have at least a few thousand absorbing atoms involved in the fitting procedure to ensure good quality of the calculated spectra, Radial Distribution Function (RDF) curves, and reliable results of the spectra approximation [20]. For the RMC procedure, we used 128 ensembles of 432 clusters placed at significant distance from each other to mimic uniform distribution of dopants in sphalerite.
To obtain the reference RDFs, we performed a fit of the EXAFS spectrum of pure ZnS measured at the Zn K-edge. The RMC-EXAFS procedure for pure sphalerite was applied for a 4 × 4 × 4 supercell with the periodic boundary condition. In all the calculations, the amplitude reduction factor S 0 2 was fixed at values estimated for the references in [1]: S 0 2 (Zn) = 0.85, S 0 2 (Cu) = 0.75, S 0 2 (In) = 0.95.
As a result of the RMC optimization, the atomic coordinates were derived. Based on these data, the RDFs were built for each pair of atoms, and the structural parameters (interatomic distances) were obtained. Gaussian functions were used to approximate the calculated RDF curves. Coordination numbers associated with Gaussian decomposition were calculated by integration of the RDFs. Analysis of the RDF shape provided information which cannot be derived by means of the conventional EXAFS fitting procedure, including detection of small distortion of the local atomic structure.

XANES Spectra Calculation
Two series of XANES spectra calculations were performed. The first series included In and Cu K-edge XANES spectra for which the models of the local atomic geometry around the absorbing atom were obtained by means of DFT geometry optimization. The second series (only Cu K-edge XANES) was based on the atomic cluster geometries derived by means of RMC-EXAFS fitting. To model the RMC-EXAFS spectra, about 100 random geometries obtained by RMC-EXAFS fitting were averaged.
Theoretical calculations of Cu and In K-edge XANES spectra were performed using the finite difference method (FDM) employing the FDMNES program [21,22]. Relativistic self-consistent field FDMNES calculations were carried out with the exchange-correlation part of potential in a local density approximation [23]. The final electronic states were calculated in a full core-hole screening. Atomic clusters inside the spheres with radii of 5 Å were used in the calculations. Many-body effects and core-hole lifetime broadening were accounted for in the frame of the arctangent convolution model [24]. Parameters of the XANES spectra convolution were chosen to obtain the best co-incidence of the experimental and calculated data.

DFT Calculations
The optimization of crystallographic parameters of the pure ZnS obtained for a 3 × 3 × 3 supercell yields the lattice parameter a DFT = 5.439 Å, which is slightly bigger than the experimental XRD value of a XRD = 5.410 Å, with relative deviation of 0.5%. This accuracy level is high enough to model the local atomic geometry and compare it to the structural parameters derived from EXAFS spectra.
In case of substitution of Zn atoms by dopants, the atomic structure of sphalerite is preserved, but the interatomic distances are relaxed to match the In and Cu ionic radii. The close arrangement of dopants and/or Zn vacancies results in the significant distortion of the coordination shells, which is translated to the splitting of the interatomic distances. Table S1 (Supplementary Materials) shows the interatomic distances of the optimized structure of ZnS around the In and Cu dopants. Localization of dopants (Me = In, Cu) and vacancy in the nearest cationic positions causes splitting of the Me-Zn distances up to 0.12 and 0.16 Å in the 2nd and 3rd coordination shells, respectively. This splitting is sufficiently high to be detected by EXAFS spectroscopy. The range of variation of the interatomic distances in the 1st coordination shell, which consists of four S atoms, is less pronounced (<0.05 Å). It is difficult to determine the contributions of individual atoms, which are characterized by small variations of interatomic distances within the 1st coordination shell, by means of EXAFS spectra analysis. Accordingly, the effect of the close arrangement of defects in sphalerite structure-In and a Zn vacancy, and In and Cu atoms-results in important splitting of interatomic distances in the 2nd and 3rd coordination shells, which can be detected and characterized by means of EXAFS.
Calculated values of the formation energy of sphalerite solid solutions are presented in Table 1. The positive sign of E EF obtained in all calculations implies endothermic character of the substitution reactions. Embedding two In atoms in ZnS without a vacancy creation does not obey the charge compensation scheme and is characterized by high value of E EF . Therefore, this configuration is not stable. Comparison of E EF calculated for different arrangements of In atom and a vacancy in the cationic sublattice shows that the location of In and a vacancy in the nearest cationic positions is energetically favorable. The same conclusion can be reached for Cu and In atoms: the close position of these dopants is more stable because of low values of E EF .
To summarize the results of the DFT calculations, the close arrangement of the impurities-In and Cu atoms, or In atom and a vacancy-is energetically more favorable than the random distribution of the defects. The close arrangement of the defects causes distortion of the coordination polyhedra of the 2nd and 3rd coordination shells which is high enough to be detected by means of EXAFS. These results will be compared with the data derived from the experimental EXAFS spectra treatment.

EXAFS Spectra Analysis (RMC-EXAFS)
The main structural information which is obtained from the RMC-EXAFS spectra fitting is the shape and position of RDF peaks, which are calculated using atomic coordinates produced by the fit. In the primary stage of the calculations, the Zn K-edge EXAFS spectrum of pure ZnS was fitted using the RMC method. Figure 1 compares the experimental and calculated spectra of ZnS. The RDF consists of three distinct peaks which correspond to the 1st, 2nd, and 3rd coordination shells. These peaks are accurately fitted by Gaussian functions (optimized Gaussian functions are shown by solid curves in Figure 1). The Gaussian curves that correspond to the 2nd coordination shell (Zn-Zn RDF, red peak in Figure 1) will be further compared with RDFs derived by the analysis of In and Cu K-edge EXAFS spectra fitting. The centroid positions and variances of Gaussian functions, which fit coordination shells around Zn in the ZnS structure, are presented in Table 2.  The next step consisted of fitting EXAFS spectra measured at K-edges of dopants in Samples 3757 (In K-edge), and 4186 and 4108 (In K-edge, Cu K-edge). The results of the fits are shown in Figure 2 (In K-edge spectra) and Figure 3 (Cu K-edge spectra). The positions of the Gaussian curves, which were used to approximate the RDFs, are listed in Table 2. Visual examination of the images in Figures 2 and 3 shows good agreement between the experimental data and model spectra. The RDFs are described by three distinct peaks which correspond to coordination shells of the probed cation with 4 S, 12 Zn, and 12 S atoms. Images of Cu K-edge EXAFS fitting with In in the 2nd coordination shell and corresponding RDFs are shown in Supplementary Materials ( Figures S2 and S3). The coordination shell radii derived from EXAFS fitting for different samples measured at the same edge are very close to each other, which demonstrate high stability of the RMC fitting procedure. Adding In to the 2nd coordination shell doesn't improve the description of the experimental EXAFS data (compare the fit results in Figure 3 and Figure S2), which is in agreement with the performed wavelet analysis applied to the collected EXAFS data [2].
The first peak in the RDF plots, which represents the mean bond length between the central atom (Cu or In) and its nearest neighbors (4 S atoms), is of asymmetrical and sharp shape. The peak asymmetry can be explained by strong distortion of the nearest atomic coordination or limited k-range of measured EXAFS spectra, especially in the case of the Cu-K edge. We believe that the limited k-range of the spectra is the more probable reason for the observed asymmetrical peak shape. The second peak of RDF functions at~3.9 Å (In K-edge) and 3.84 Å (Cu K-edge) corresponds to the 2nd coordination shell, which was approximated by two models-12 Zn or (only for Cu K-edge) 11 Zn + 1 In atoms. The last peak in RDFs corresponds to 12 S atoms located in the 3rd coordination shell. As follows from data presented in Table 2, the mean radii of the three coordination shells determined by RDF peak fittings with Gaussian functions are in good agreement with the results of conventional EXAFS spectra fits (Trofimov et al. [2]).
The most important information concerning the atomic arrangement around a dopant can be obtained by exploring the 2nd coordination shell. The peaks, which correspond to the metal atoms located in the 2nd coordination shells of Cu and In, have asymmetrical shape for all fitted spectra and cannot be accurately approximated by single Gaussian function (Figures 2 and 3). The asymmetry can be accounted for by decomposition of the peaks into two Gaussian components for the In-Zn RDFs of Samples 3757 and 4108, and three Gaussian components in the case of the In-Zn RDF of Sample 4186 and Cu-Zn RDFs. The results of the decomposition are given in Table 3 and are shown in Figure 4.  The first peak in the RDF plots, which represents the mean bond length between the central atom (Cu or In) and its nearest neighbors (4 S atoms), is of asymmetrical and sharp shape. The peak asymmetry can be explained by strong distortion of the nearest atomic coordination or limited krange of measured EXAFS spectra, especially in the case of the Cu-K edge. We believe that the limited k-range of the spectra is the more probable reason for the observed asymmetrical peak shape. The second peak of RDF functions at ~3.9 Å (In K-edge) and 3.84 Å (Cu K-edge) corresponds to the 2nd coordination shell, which was approximated by two models-12 Zn or (only for Cu Kedge) 11 Zn + 1 In atoms. The last peak in RDFs corresponds to 12 S atoms located in the 3rd coordination shell. As follows from data presented in Table 2, the mean radii of the three coordination shells determined by RDF peak fittings with Gaussian functions are in good agreement with the results of conventional EXAFS spectra fits (Trofimov et al. [2]).
The most important information concerning the atomic arrangement around a dopant can be obtained by exploring the 2nd coordination shell. The peaks, which correspond to the metal atoms located in the 2nd coordination shells of Cu and In, have asymmetrical shape for all fitted spectra and cannot be accurately approximated by single Gaussian function (Figures 2 and 3). The  The positions of the three Gaussian functions, which were used to fit the Cu-Zn RDFs in Samples 4108 and 4186 (Cu K-edge EXAFS), are very stable and only slightly change from one sample to another. The significant splitting of the 2nd coordination shell of Cu results in the formation of three groups of Zn atoms, two of which are shifted by 0.16-0.2 Å from the average position of 3.84 Å (the average Cu-Zn distance is close to the Zn-Zn distance of 3.83 Å in pure sphalerite, Table 2). All Cu K-edge EXAFS spectra were accurately fitted without In in the 2nd coordination shell. Addition of an In atom does not improve the agreement between the experimental and calculated spectra. The shapes of the RDFs derived from the EXAFS spectra fitting with and without In are close to each other ( Figure S4). Since the WTs ( Figure S2) do not reveal any heavy atoms in the 2nd coordination shell of Cu, we can conclude that the 2nd coordination shell consists mostly of Zn atoms.
The results of RMC-EXAFS spectra fittings demonstrate that the structures of the local atomic environments of In and Cu dopants are not identical. The local atomic symmetry of In depends on the concentration of Cu. The local atomic structure of In in sphalerite with only In impurities is close to one of Zn in pure sphalerite. An increase in Cu and In concentration causes distortion of the local atomic structure of In. In contrast, the interatomic distances in the 2nd coordination sphere of Cu exhibit considerable splitting, which doesn't depend on the amount of impurities in sphalerite. This splitting cannot be explained by the close arrangement of the dopants (In and Cu) as (i) no heavy In atoms were detected in the second shell of Cu (Trofimov et al. [2]), (ii) inclusion of In atoms in the model does not improve the quality of the RMC-EXAFS fitting, and (iii) the 2nd shell of In in Sample 4108, with strong distortion of atomic structure around Cu, is almost undistorted. We note that in the first two samples (3757, 4108), the coordination numbers associated with the distorted Zn atoms in the 2nd coordination shell of In do not exceed N = 1 (Table 3). This means that the distortion of the 2nd coordination shell of In is almost negligible. Sample 4186 contains a high amount of Cu and In impurities (nearly 2 mol.% of Cu and In [2]), and in this case, the 2nd coordination shell around In is much more distorted.
The positions of the three Gaussian functions, which were used to fit the Cu-Zn RDFs in Samples 4108 and 4186 (Cu K-edge EXAFS), are very stable and only slightly change from one sample to another. The significant splitting of the 2nd coordination shell of Cu results in the formation of three groups of Zn atoms, two of which are shifted by 0.16-0.2 Å from the average position of 3.84 Å (the average Cu-Zn distance is close to the Zn-Zn distance of 3.83 Å in pure sphalerite, Table 2).
All Cu K-edge EXAFS spectra were accurately fitted without In in the 2nd coordination shell. Addition of an In atom does not improve the agreement between the experimental and calculated spectra. The shapes of the RDFs derived from the EXAFS spectra fitting with and without In are close to each other ( Figure S4). Since the WTs ( Figure S2) do not reveal any heavy atoms in the 2nd coordination shell of Cu, we can conclude that the 2nd coordination shell consists mostly of Zn atoms.
The results of RMC-EXAFS spectra fittings demonstrate that the structures of the local atomic environments of In and Cu dopants are not identical. The local atomic symmetry of In depends on the concentration of Cu. The local atomic structure of In in sphalerite with only In impurities is close to one of Zn in pure sphalerite. An increase in Cu and In concentration causes distortion of the local atomic structure of In. In contrast, the interatomic distances in the 2nd coordination sphere of Cu exhibit considerable splitting, which doesn't depend on the amount of impurities in sphalerite. This splitting cannot be explained by the close arrangement of the dopants (In and Cu) as (i) no heavy In atoms were detected in the second shell of Cu (Trofimov et al. [2]), (ii) inclusion of In atoms in the model does not improve the quality of the RMC-EXAFS fitting, and (iii) the 2nd shell of In in Sample 4108, with strong distortion of atomic structure around Cu, is almost undistorted. Therefore, despite the splitting of interatomic distances in the 2nd shell of Cu, RMC-EXAFS results argue for random (or at least, not close) arrangement of the dopants-In and Cu. The disagreement with the results of DFT calculations, according to which the closest arrangement of dopants is the most stable configuration, can be explained by the presence of other defects of the sphalerite crystal lattice, which were omitted during the DFT calculations.

In K-Edge
Since XANES spectra of Samples 3757, 4108, and 4186 are identical [2], we used the experimental XANES spectrum of Sample 3757 (In-bearing sphalerite) as a reference to compare with calculated In K-edge XANES spectra. The simulations of In K-edge XANES spectra were performed using atomic configurations optimized by means of DFT. The In-S distance calculated by DFT differs slightly from the one derived by EXAFS analysis. Therefore, XANES calculations were performed with atomic coordinates corrected to match the In-S distance obtained by EXAFS. The calculation of XANES spectra were performed for (i) sphalerite crystal structure with two Zn atoms substituted by In and Cu atoms ( Figure 5a); (ii) sphalerite crystal structure with two Zn atoms substituted by two In atoms and one Zn vacancy ( Figure 5b); (iii) sphalerite crystal structure with two Zn atoms substituted by two In atoms (this structural model is beyond the charge compensation scheme) (Figure 5c). Calculations for the three models listed above were performed with dopants and vacancy placed in the nearest positions and far from each other to check the influence of the crystal structure defects arrangement on XANES spectra. As can be seen in Figure 5a-c, all calculated spectra reproduce correctly all features of the experimental spectrum. The difference between the calculations performed for the different structural models is negligible, which can be explained by the significant natural broadening of In K-edge spectra.

Cu K-Edge
The Cu K-edge XANES spectrum of Sample 4186 (In-Cu-bearing sphalerite) was used as a reference to compare with the calculated spectra. The results of the Cu K-edge XANES calculations are shown in Figure 6. The figure compares the results of simulations which are based on atomic geometries obtained by DFT optimization and RMC-EXAFS fitting. In both cases, calculations were performed on structures with and without In in the 2nd coordination shell. The experimental XANES spectrum contains spectral features A, B, C, and D. All the structures used and both the calculation methods (DFT and RMC-EXAFS) reproduce the overall shape of the experimental spectrum. However, the geometries obtained by RMC-EXAFS fitting provide better agreement with the experimental data regarding the energy position and intensity of spectral features B and C. Another useful conclusion follows from the comparison of calculated spectra of structures with and without In in the 2nd coordination shell: influence of one In atom on the shape of the XANES Figure 5. Comparison of experimental (Sample 3757) and calculated In K-edge XANES spectra. The calculation was performed for different models of In-bearing sphalerite obtained by DFT: (a) 2Zn 2+ ↔Cu 1+ + In 3+ substitution; (b) 2Zn 2+ ↔In 3+ + substitution; (c) 2Zn 2+ ↔2In 3+ substitution. Black lines are used to plot the experimental spectra of Sample 3757; red lines-calculated spectra for structures with impurities and/or vacancy placed far from each other; green lines-calculated spectra for the structures with impurities and vacancy placed in the nearest positions.

Cu K-Edge
The Cu K-edge XANES spectrum of Sample 4186 (In-Cu-bearing sphalerite) was used as a reference to compare with the calculated spectra. The results of the Cu K-edge XANES calculations are shown in Figure 6. The figure compares the results of simulations which are based on atomic geometries obtained by DFT optimization and RMC-EXAFS fitting. In both cases, calculations were performed on structures with and without In in the 2nd coordination shell. The experimental XANES spectrum contains spectral features A, B, C, and D. All the structures used and both the calculation methods (DFT and RMC-EXAFS) reproduce the overall shape of the experimental spectrum. However, the geometries obtained by RMC-EXAFS fitting provide better agreement with the experimental data regarding the energy position and intensity of spectral features B and C. Another useful conclusion follows from the comparison of calculated spectra of structures with and without In in the 2nd coordination shell: influence of one In atom on the shape of the XANES spectrum is negligible.

Discussion
In the present study, the local atomic geometry of In and Cu in sphalerite ZnS was determined by the reverse Monte Carlo (RMC) method applied to EXAFS spectra approximation. The EXAFS spectra of synthetic In-and In-Cu-bearing sphalerites were fitted by means of RMC with determination of the interatomic distances (for the 2nd shell atoms) in the first three coordination shells of dopants. The advantage of the RMC-EXAFS method in comparison with conventional EXAFS spectra fitting (realized, for example, in the IFEFFIT software) is the possibility to obtain reliable information about small distortion of coordination shells. In general, the results of the present study are in agreement with our data obtained previously by conventional EXAFS spectra fittings. Incorporation of In into the cationic position of ZnS crystal structure causes expansion of the nearest-to-In coordination polyhedra. The bond lengths increase (relatively to pure ZnS) from 2.34 to 2.45 Å in the 1st coordination shell (NS = 4), from 3.85 to 3.91 Å in the 2nd coordination shell (NZn = 12), and are close to the Me-S distance in the pure sphalerite in the 3rd coordination shell (NS = 12, RIn-S = 4.48 Å ). The In Radial Distribution Function (RDF) consists of three distinct peaks, which correspond to the first three coordination shells. We observe no splitting of the coordination shells of In into subshells when In is the only dopant. The presence of Cu in sphalerite distorts the local atomic surrounding of In, and the distortion extent of the In-Zn coordination shell depends on the concentration of impurity atoms. In Sample 4108 (concentration of In and Cu is less than 0.1 mol%), the shape of In-Zn RDF is close to the one of In-Zn for the sphalerite with only In dopants. An increase in concentrations of In and Cu causes the splitting of In-Zn RDF into three contributions.
The Cu impurity in the In-Cu-bearing sphalerite also presents in the solid solution state substituting for Zn. Since the ionic radii of Zn and Cu are close to each other, only small changes of the average coordination shell sizes take place around Cu compared to the pure sphalerite.

Discussion
In the present study, the local atomic geometry of In and Cu in sphalerite ZnS was determined by the reverse Monte Carlo (RMC) method applied to EXAFS spectra approximation. The EXAFS spectra of synthetic In-and In-Cu-bearing sphalerites were fitted by means of RMC with determination of the interatomic distances (for the 2nd shell atoms) in the first three coordination shells of dopants. The advantage of the RMC-EXAFS method in comparison with conventional EXAFS spectra fitting (realized, for example, in the IFEFFIT software) is the possibility to obtain reliable information about small distortion of coordination shells. In general, the results of the present study are in agreement with our data obtained previously by conventional EXAFS spectra fittings. Incorporation of In into the cationic position of ZnS crystal structure causes expansion of the nearest-to-In coordination polyhedra. The bond lengths increase (relatively to pure ZnS) from 2.34 to 2.45 Å in the 1st coordination shell (N S = 4), from 3.85 to 3.91 Å in the 2nd coordination shell (N Zn = 12), and are close to the Me-S distance in the pure sphalerite in the 3rd coordination shell (N S = 12, R In-S = 4.48 Å). The In Radial Distribution Function (RDF) consists of three distinct peaks, which correspond to the first three coordination shells. We observe no splitting of the coordination shells of In into subshells when In is the only dopant. The presence of Cu in sphalerite distorts the local atomic surrounding of In, and the distortion extent of the In-Zn coordination shell depends on the concentration of impurity atoms. In Sample 4108 (concentration of In and Cu is less than 0.1 mol%), the shape of In-Zn RDF is close to the one of In-Zn for the sphalerite with only In dopants. An increase in concentrations of In and Cu causes the splitting of In-Zn RDF into three contributions.
The Cu impurity in the In-Cu-bearing sphalerite also presents in the solid solution state substituting for Zn. Since the ionic radii of Zn and Cu are close to each other, only small changes of the average coordination shell sizes take place around Cu compared to the pure sphalerite. However, analysis of Cu RDF shows significant deformation of the 2nd coordination shell. The Cu-Zn RDF peak, which corresponds to the 2nd coordination shell of the cation, can be decomposed into three Gaussian contributions with centroids at~3.64,~3.80, and~3.99Å, which means considerable deformation of ZnS crystal structure near Cu impurity.
The In and Cu K-edge XANES spectra were simulated using structural models obtained by the DFT and RMC-EXAFS methods. It was found that the deformed coordination polyhedra around Cu are better reproduced by the RMC-EXAFS method, which provides better agreement between the simulated and the experimental spectra. The XANES spectra calculations confirm conserving of the of ZnS structure near In impurity and its deformation near Cu impurity.
Substitution In 3+ →Zn 2+ results in the formation of a positively charged site located on the dopant. On the contrary, the Cu + →Zn 2+ substitution results in the formation of a negatively charged site. The DFT calculations show that the close arrangement of In and Cu is more energetically favorable than random distribution of dopants, which can be explained by the local charge balance. The DFT shows significant distortion in ZnS crystal structure around the dopants (In and Cu). The splitting of interatomic distances, calculated by means of DFT, reaches 0.10 Å for the 2nd shell around Cu and 0.12 Å for the In 2nd shell. This coordination shell splitting should be observed by means of EXAFS, even in the samples with low concentrations of In and Cu. Analysis of the experimental EXAFS spectra shows considerable deformation of ZnS geometry around Cu, which doesnot depend on the concentration of Cu or In. In the case of In, the local atomic distortion depends on the concentration of impurities, and its extent correlates with the amount of dopants in sphalerite. This implies that Cu and In impurities in sphalerite with low concentration of dopants are placed at a significant distance from each other. The apparent disagreement between the results of the DFT calculations and experimental EXAFS spectra fittings can stem from the oversimplified model of doped sphalerite structure used in the DFT calculations. In the calculations, we employed the simplest models with defects localized at the atoms of In, Cu, or vacancy. The real structure of doped sphalerite can be more complicated and contain defects, which were not considered in our calculations.

Conclusions
The reverse Monte Carlo (RMC) simulation was applied to interpret the In and Cu K-edge EXAFS spectra of In-and In-Cu-bearing sphalerites. In the studied samples, both In and Cu occur in the solid solution state, substituting for Zn. The local atomic structure of In in ZnS without Cu is undistorted. The presence of both dopants in ZnS causes variations of the interatomic distances, which is clearly seen in the splitting of the 2nd coordination spheresof In and Cu. The degree of distortion of the 2nd coordination sphere of In correlates with the concentration of In + Cu, but is weak in general. The degree of distortion of the 2nd coordination sphere of Cu shows no correlation with the concentration of the dopants, but is much more pronounced than that of In. The results of RMC-EXAFS fittings argue for a random (or at least not close) arrangement of the dopants-In and Cu. In contrast to the RMC-EXAFS method, the Density Functional Theory (DFT) calculations showed that the close arrangement of In and Cu is energetically favorable. We attribute the disagreement between the RMC-EXAFS and DFT methods to the presence of defects or vacancies in ZnS structure in close vicinity of the dopant sites, which were not taken into account in the DFT calculations. The nature and types of such defects are the questions of further studies.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-163X/10/10/841/s1, Figure S1: Top: Coordination polyhedra of the cations in the structures of sphalerite ZnS [3] and roquesite CuInS 2 [4] (Adopted from Trofimov et al. [2]). Bottom: Crystal structures of sphalerite (doubled cell is shown to make comparison with CuInS 2 easier) and roquesite [4], Figure S2: The RMC fit results of Cu K-edge EXAFS spectra of samples 4108 and 4186 (the 2nd coordination shell of Cu is described by Zn + In atoms). Panel A: experimental WT image; panel B: fitted WT image; panel C: calculated RDF for first three coordination shells around absorbing atom; panel D: experimental (black dotted line) and fitted (red solid line) EXAFS signal χ(k)*k2; panel E: Fourier transform magnitudes of experimental (black dotted line) and fitted (red solid line) χ(k)*k2 function; panel F: close view of Cu-Zn and Cu-In RDF. The RDFs fitted by Gaussians are shown by solid lines, Figure S3: Comparison of Cu-Zn and In-Zn RDFs with Zn-Zn RDF derived from RMC EXAFS fits. The model assumes In atoms in the 2nd coordination shell of In and Cu, Figure S4: Comparison of Cu-Zn and Cu-(Zn + In) RDFs for Cu K-edge EXAFS fitting with and without In in the 2nd coordination shell of Cu. Note that the shapes of the curves representing the two models (11Zn + In and 12Zn) are very close to each other which means that incorporation of In atom in the 2nd sphere of Cu is not necessary for the accurate approximation of the experimental data, Table S1: Interatomic distances in sphalerite with and without dopants determined by DFT calculations. Literature data on the interatomic distances of unrelaxed pure sphalerite structure are given at the bottom of the table. Data for the systems ZnS + In was obtained in the present study, all other data are adopted from [1,2]. The method of the calculations is described in [1].