Cation Disorder and Local Structural Distortions in AgxBi1–xS2 Nanoparticles

By combining X-ray absorption fine structure and X-ray diffraction measurements with density functional and molecular dynamics simulations, we study the structure of a set of AgxBi1−xS2 nanoparticles, a materials system of considerable current interest for photovoltaics. An apparent contradiction between the evidence provided by X-ray absorption and diffraction measurements is solved by means of the simulations. We find that disorder in the cation sublattice induces strong local distortions, leading to the appearance of short Ag–S bonds, the overall lattice symmetry remaining close to hexagonal.


Introduction
Dimetal chalcogenides of general formula I−V−VI 2 , in which I = Cu, Ag, or an alkali metal, V = Sb or Bi, and VI = S, Se, or Te, have recently attracted attention as an interesting class of materials for energy conversion from light (photovoltaics) and heat (thermoelectrics). Solar cells employing AgBiS 2 nanoparticle (NP) photo-absorbing films of just 40 nm in thickness exhibit 6.4% certified photovoltaic efficiencies [1]. This unique feature is attributed to the material's high optical extinction coefficient (>10 5 cm −1 ), and to a nearly 100% internal quantum efficiency across the VIS and NIR parts of the optical spectrum; among other attractive features of the material is its environmental friendliness. The photovoltaic performance of AgBiS 2 NPs is currently limited by the presence of electronic midgap trap states, which limit the maximum allowable thickness of films within a solar cell. The structural origin of these traps remains unknown and could be related to how Ag and Bi cations are distributed within the cation sub-lattice and/or to defects associated with structural distortions. Both AgBiS 2 and AgBiSe 2 have been demonstrated to have attractive thermoelectric properties due to the very low lattice thermal conductivity, a property related to bond anharmonicity [2][3][4].
The low-temperature phase of AgBiS 2 is that of the mineral matildite (hexagonal lattice, space group 164, P3m1), while above 179 • C, the stable phase is that of schapbachite (cubic lattice, space group 225, i.e., Fm3m) [5]. The two structures are related and differ in the relative arrangement of cations [6]. The matildite structure can be described in terms of a large cubic unit cell within which a smaller hexagonal one is found; the S atoms form a face-centered cubic (FCC) arrangement and Ag and Bi atoms are alternatively placed along the cubic (001) directions, forming a rocksalt structure. In schapbachite, the S atoms are also found in an FCC arrangement, but layers of BiS and AgS alternate in the (001) directions. Using density functional theory (DFT) simulations, Viñes et al. [6] predicted that the matildite structure is more stable than the schapbachite one and that the former is a semiconductor while the latter is a metal; the high photoactivity of the matildite structure has been attributed to different effective masses of electrons in the conduction band and holes in the valence band and a high optical absorption coefficient.
The relation between physical properties and the ordering of cations in dimetal chalcogenides has been an important theme of current research. An order-disorder phase transition in cubic AgBiS 2 nanocrystals, proposed to explain calorimetric measurements, has been suggested to be responsible for the high thermoelectric efficiency [2], but no direct structural proof was provided. The gradual transition from matildite to schapbachite by means of cation rearrangement has been simulated by means of DFT [7]. Cation rearrangement from matildite to schapbachite was found to induce a gradual decrease in the band gap as a result of conduction band states shifting downward toward the valence band; intermediate cases of ordering were suggested to result in significant local structural disruption, but this issue was never addressed experimentally. Additionally, small perturbations in the Ag to Bi ratio (that is, in the stoichiometry) originating from variations in the synthesis conditions may introduce further local distortions. A similar phase transition was also observed in the selenide analog of AgBiS 2 , namely, AgBiSe 2 , in which the cation arrangement related to a hexagonal-rhombohedral-cubic phase transition was associated with switching between different types of conductivity (p or n) and that a high thermoelectric performance was obtained for a disordered cation arrangement [3]. For a disordered arrangement in the cation sublattice, it can be expected to result in a similar effect on the properties of AgBiS 2 . Thus, to have an atomistic understanding of this material's properties and to guide the experimental efforts for improving the performance of AgBiS 2 -based technologies, it is crucial to explore the effects of cation disorder and stoichiometric perturbations in the Ag to Bi ratio.
In this paper, we address the relation between cation ordering and local structural distortions in a set of Ag x Bi 1−x S 2 NP thin films with purposely modified compositions, slightly below and above the stoichiometric 1:1 Ag to Bi ratio. We have used X-ray absorption fine structure (XAFS) and X-ray diffraction (XRD) aided by DFT ab initio and molecular dynamics (MD) simulations. Owing to its chemical selectivity and sensitivity, XAFS [8] is the premier experimental tool to measure interatomic distances in the first few coordination shells and to probe the relative atomic arrangement in condensed matter, in particular, in alloys; a further attractive feature is the applicability to ordered or disordered materials both in the bulk or as nanostructures. A limitation of XAFS is its insensitivity to lattice symmetry, which, along with the local atomic arrangement, also plays a fundamental role in determining physical properties. XRD is a standard characterization tool in this regard and its applicability to nanostructures has been recently greatly expanded, showing how refined information on symmetry, lattice strain, particle size, and shape can be obtained [9][10][11]. In this paper, the interpretation of the experimental data is greatly aided by DFT [12] and MD simulations, the efficiency and availability of which have constantly improved in recent years. We show that in Ag x Bi 1−x S 2 NPs, disorder in the cation distribution leads to significant local structural distortions. An apparent contradiction between XAFS and XRD is solved using DFT and MD simulations.

Materials and Methods
Colloidal Ag x Bi 1−x S 2 NP thin films were synthesized by a hot-injection low-temperature approach [1]. In this method, 0.8 mmol of silver acetate and 1 mmol of bismuth acetate are reacted with Nanomaterials 2020, 10, 316 3 of 13 oleic acid; then, 1 mmol of bis (trimethylsilyl) sulfide (TMS) diluted with octadecene is injected into silver-bismuth oleate solution at 100 • C to form the AgBiS 2 nanocrystals. To synthesize Ag x Bi 1−x S 2 NPs with four different x values, we adjusted the amount of silver precursor introduced into the reaction flask, keeping the other parameters of the reaction constant. Energy-dispersive X-ray (EDX) spectra acquisition and analysis were performed with a Helios Nanolab 600 (FEI Company) microscope combined with an X-Max detector and INCA ® system (Oxford Instruments). High-resolution transmission electron microscopy (TEM) measurements were performed at the Scientific and Technological Centers of the University of Barcelona (CCiT-UB). TEM micrographs were obtained using a JEOL 2100 microscope operating at an accelerating voltage of 200 kV. The samples were prepared by drop-casting a toluene-based colloidal NP solution on a holey carbon-coated grid. TEM ( Figure 1) indicated an average size of 3.5-4 nm; the particles tended to show oval projections, not perfectly rounded, with a short axis of about 3 nm and a long axis of 4-5 nm; the thickness for all samples was~100 nm. The correspondence between sample code and Ag concentration is summarized in Table 1, in which we report the reaction precursor ratio and values measured by EDX. Nanomaterials 2020, 10, x FOR PEER REVIEW 3 of 14 reaction flask, keeping the other parameters of the reaction constant. Energy-dispersive X-ray (EDX) spectra acquisition and analysis were performed with a Helios Nanolab 600 (FEI Company) microscope combined with an X-Max detector and INCA ® system (Oxford Instruments). Highresolution transmission electron microscopy (TEM) measurements were performed at the Scientific and Technological Centers of the University of Barcelona (CCiT-UB). TEM micrographs were obtained using a JEOL 2100 microscope operating at an accelerating voltage of 200 kV. The samples were prepared by drop-casting a toluene-based colloidal NP solution on a holey carbon-coated grid. TEM ( Figure 1) indicated an average size of 3.5-4 nm; the particles tended to show oval projections, not perfectly rounded, with a short axis of about 3 nm and a long axis of 4-5 nm; the thickness for all samples was ~100 nm. The correspondence between sample code and Ag concentration is summarized in Table 1, in which we report the reaction precursor ratio and values measured by EDX.  XAFS measurements at the Ag K (25.514 keV) and Bi LIII (13.409 keV) edges were performed at the BM08 (LISA) beamline of the European Synchrotron Radiation Facility (ESRF), Grenoble, France [13]. All measurements were carried out using a dynamically sagittal-focusing Si (311) monochromator [14], and the spectra were recorded in the fluorescence yield mode using a 12 element hyper-pure Ge detector [15]. Ag2S and Bi2S3 polycrystalline powders were measured in the transmission mode as references. The structure of Ag2S is that of the mineral acanthite (monoclinic lattice, space group 14, i.e., P21/c) with two inequivalent Ag sites present, see Figure 2, left panel. XAFS data analysis in the extended energy range (EXAFS) was performed using the DEMETER package [16] and theoretical signals based on FEFF [17]. Simulations of the X-ray absorption nearedge structure (XANES) part of the XAFS spectra were performed using FDMNES [18].
X-ray diffraction measurements were performed using a PANalytical X'Pert PRO MPD Alpha1 powder diffractometer with Cu Kα radiation (λ = 1.5406 Å , 45 kV-40 mA) in the Scientific and Technological Centers of the University of Barcelona (CCiT-UB). The samples were prepared by drop-casting a concentrated dispersion of AgxBi1−xS2 NPs onto a clean glass slide. The XRD patterns were analyzed with the TOPAS software for structural refinement [19], version 7, using the functions  XAFS measurements at the Ag K (25.514 keV) and Bi L III (13.409 keV) edges were performed at the BM08 (LISA) beamline of the European Synchrotron Radiation Facility (ESRF), Grenoble, France [13]. All measurements were carried out using a dynamically sagittal-focusing Si (311) monochromator [14], and the spectra were recorded in the fluorescence yield mode using a 12 element hyper-pure Ge detector [15]. Ag 2 S and Bi 2 S 3 polycrystalline powders were measured in the transmission mode as references. The structure of Ag 2 S is that of the mineral acanthite (monoclinic lattice, space group 14, i.e., P2 1 /c) with two inequivalent Ag sites present, see Figure 2, left panel. XAFS data analysis in the extended energy range (EXAFS) was performed using the DEMETER package [16] and theoretical signals based on FEFF [17]. Simulations of the X-ray absorption near-edge structure (XANES) part of the XAFS spectra were performed using FDMNES [18]. X-ray diffraction measurements were performed using a PANalytical X'Pert PRO MPD Alpha1 powder diffractometer with Cu K α radiation (λ = 1.5406 Å, 45 kV-40 mA) in the Scientific and Technological Centers of the University of Barcelona (CCiT-UB). The samples were prepared by drop-casting a concentrated dispersion of Ag x Bi 1−x S 2 NPs onto a clean glass slide. The XRD patterns were analyzed with the TOPAS software for structural refinement [19], version 7, using the functions recently developed for the analysis of nanocrystalline powders [20]. Different structures were implemented in the refinement to assess the degree of matching to the experimental data. Structural information on standard phases was taken from the "materials project" web interface [21]. In addition, to fit quality and structural information, data analysis provides the crystalline domain size, assuming spherical shape with a lognormal distribution of diameters, and inhomogeneous strain (also referred to as microstrain). Details on the procedure are reported elsewhere [20,22].
Tentative structures for the clusters were simulated using the DFT method, as implemented in the VASP code (version 5.2) [12]. The structures considered were the following: Bi-substituted acanthite (AcB), schapbachite (Sch), matildite (Mat), and matildite with random cation site occupation (MaR). Pure acanthite (Aca), Bi 2 S 3 (BiS), and sulfur (S) were also studied for the determination of the formation energies. The AcB cell was derived from that of Aca by substituting a Ag-II + ion with a Bi 3+ one and eliminating two Ag + ions, which were located as far as possible from Bi 3+ and from each other to ensure charge balance. A calculation carried out starting from Ag-I + led to a similar local structure and total energy, so these two structures will be considered as equivalent. For this compound, the first structural relaxation was carried out on a small cell (10 atoms) followed by a final relaxation on an 80 atom cell (subsequently also used for the MD run). Schapbachite was obtained by filling the cation sites randomly with Ag and Bi. A similar process was applied for random matildite: The cation sites were considered all equal and randomly filled with either Bi or Ag. Table 2 summarizes the main details of the cells used for the structural simulation. For each case, the geometry of the energy-minimized (relaxed) structure was determined, as well as the formation energy per atom and the simulated XAFS spectrum. For the structural relaxation, the PBE exchange-correlation functional [27] was used with a plane wave cut-off energy of 650 eV. For the cell energy calculation, increasingly dense k meshes were used until the energy difference between consecutive steps was less than 1/1000 of the total cell energy. The k point meshes were all centered on the Γ point and ranged from k = {3 × 3 × 3} for sulphur to k = {5 × 5 × 5} for the other compounds and k = {8 × 8 × 4} for hexagonal matildite. A Gaussian distribution with a smearing of 0.05 eV was used to populate the electronic states. The electronic density was considered converged when the difference between successive steps was less than 1 µeV. The structural relaxation process was considered converged when Feynman−Hellman forces were below 2 × 10 −4 eV/Å; every 20 ionic steps, the calculation was stopped and restarted with a new plane wave base recalculated with the new structure.
Formation energies of the cells (and related energies per atom) were calculated starting from chemical potentials of each element µ α (α = Ag, Bi, S). These potentials were obtained from the cell Nanomaterials 2020, 10, 316 5 of 13 energies E y of acanthite, Bi 2 S 3 , and S, and the number of each species in the cell N w y (see column 2, Table 2), solving the linear problem [28]: This means that the chemical potentials µ α are calculated in a condition in which acanthite, Bi 2 S 3 , and S are at equilibrium.
The MD simulations of the various materials, necessary for the simulation of the XAFS spectra, were carried out in the framework of the DFT method, as implemented in VASP. The starting points were the relaxed cells obtained previously. The runs lasted about 10 ps (longer than the typical optical phonon periods, expected to be <1 ps) and were carried out in steps of 1 fs. Although no structural evolution was expected in the MD runs, we have verified that the total energy in the last portion of the run was stabilized and that the temperature distribution in the frames peaked at the desired value. The temperature was stabilized using a Nosé thermostat set at 300 K and monitored along the full run.
The structures from the last 2.9 ps of the run were used for the simulation of EXAFS spectra at the Ag K-edge [29]. In particular, from 290 single structures (1 taken every 10 fs), a single EXAFS spectrum was calculated with the FEFF 8.4 code [30], each Ag atom in the cell contributing with its own spectrum. Single and multiple scattering paths were considered up to a maximum path length of 10 Å and a maximum scattering order of 4, and the data range in k space was 0-20 Å −1 . The data were then averaged and the final residual (i.e., the difference between the spectra taken with 289 and 290 frames) was about 3 × 10 −4 . The EXAFS spectra so obtained were successively filtered in the interval R = [0.4-6.8] Å in order to eliminate signals from atomic replicas generated by the periodic boundary conditions. Nanomaterials 2020, 10, x FOR PEER REVIEW 5 of 14 The MD simulations of the various materials, necessary for the simulation of the XAFS spectra, were carried out in the framework of the DFT method, as implemented in VASP. The starting points were the relaxed cells obtained previously. The runs lasted about 10 ps (longer than the typical optical phonon periods, expected to be <1 ps) and were carried out in steps of 1 fs. Although no structural evolution was expected in the MD runs, we have verified that the total energy in the last portion of the run was stabilized and that the temperature distribution in the frames peaked at the desired value. The temperature was stabilized using a Nosé thermostat set at 300 K and monitored along the full run.
The structures from the last 2.9 ps of the run were used for the simulation of EXAFS spectra at the Ag K-edge [29]. In particular, from 290 single structures (1 taken every 10 fs), a single EXAFS spectrum was calculated with the FEFF 8.4 code [30], each Ag atom in the cell contributing with its own spectrum. Single and multiple scattering paths were considered up to a maximum path length of 10 Å and a maximum scattering order of 4, and the data range in space was 0-20 Å −1 . The data were then averaged and the final residual (i.e., the difference between the spectra taken with 289 and 290 frames) was about 3 × 10 −4 . The EXAFS spectra so obtained were successively filtered in the interval R = [0.4-6.8] Å in order to eliminate signals from atomic replicas generated by the periodic boundary conditions.

Results and Discussion
The background-subtracted Ag K-edge EXAFS spectra of all samples and of Ag2S are shown in Figure 3: They are all very similar and resemble the Ag2S spectrum. The main difference is the absence of a shoulder at ~5.1 Å −1 in the samples (indicated by the arrow); presumably, this is due to the presence of Bi in the local coordination around Ag. The Fourier Transforms (FTs) of the EXAFS spectra, performed in the range of 2.0-10.5 Å −1 of Ag2S and sample 2, which are also quite similar, are shown in Figure 4 as black continuous lines. The XANES spectra at the Ag K-edge is somewhat featureless because of the core-hole broadening. The experimental spectra of the samples and Ag2S are reported in the Supplemental Material ( Figure S1) and also show similarities. We also report ( Figure S2) a comparison between the experimental spectrum and FDMNES simulations for the schapbachite (cubic), acanthite (monoclinic), and matildite (hexagonal) bulk phases. There are similarities of the experimental spectrum with the hexagonal and monoclinic phases, but it is not possible to draw firm conclusions from these simulations, also because of the core-hole broadening.

Results and Discussion
The background-subtracted Ag K-edge EXAFS spectra of all samples and of Ag 2 S are shown in Figure 3: They are all very similar and resemble the Ag 2 S spectrum. The main difference is the absence of a shoulder at~5.1 Å −1 in the samples (indicated by the arrow); presumably, this is due to the presence of Bi in the local coordination around Ag. The Fourier Transforms (FTs) of the EXAFS spectra, performed in the range of 2.0-10.5 Å −1 of Ag 2 S and sample 2, which are also quite similar, are shown in Figure 4 as black continuous lines. The XANES spectra at the Ag K-edge is somewhat featureless because of the core-hole broadening. The experimental spectra of the samples and Ag 2 S are reported in the Supplemental Material ( Figure S1) and also show similarities. We also report ( Figure S2) a comparison between the experimental spectrum and FDMNES simulations for the schapbachite (cubic), acanthite (monoclinic), and matildite (hexagonal) bulk phases. There are similarities of the experimental spectrum with the hexagonal and monoclinic phases, but it is not possible to draw firm conclusions from these simulations, also because of the core-hole broadening.  The Ag2S spectrum was fitted using the monoclinic structure described above. As mentioned, Ag has two structurally inequivalent sites and, especially for site II, there is a rather broad distribution of interatomic distances around Ag. This fact complicates the data analysis. We, therefore, adopted a simplified approach, which, despite being an approximation, is able to reproduce well all spectra. Specifically, we fitted the spectra with the local structure of only the Ag-I site as the absorbing atom; the rationale is that the average of all scattering paths giving rise to the EXAFS signal relative to both Ag sites is equivalent to the sum of scattering paths for the more ordered Ag-I site. Clearly, we do not claim that the distances and Debye−Waller factors obtained in this way reproduce faithfully the full atomic distribution around Ag; however, the adopted approach is able to reproduce key features and trends in the local structural parameters and has the advantage of simplicity.
EXAFS scattering paths were calculated ab initio using FEFF, and the highest amplitude ones were included in the fit, which was performed in the ranges = 2-10.5 Å −1 and = 1-3.5 Å . The interatomic distances, an energy origin shift, and Debye−Waller factors were considered as fitting parameters; the many body amplitude reduction factor was fixed to the value obtained from a fit of a Ag metal reference spectrum (S0 2 = 0.834). We found that the Ag2S spectrum could be fitted with single scattering paths due to two S atoms at ~2.50 Å (Ag-S1), one S atom at ~3.06 Å (Ag-S2), seven Ag atoms at ~2.97 Å (Ag-Ag1), and two more Ag atoms at ~3.11 Å (Ag-Ag2). These contributions give rise to the first two peaks in the FT of the spectra reported in Figure 4A, in which the components are labelled as just described. There are nine Ag atoms at different distances in the cation coordination  The Ag2S spectrum was fitted using the monoclinic structure described above. As mentioned, Ag has two structurally inequivalent sites and, especially for site II, there is a rather broad distribution of interatomic distances around Ag. This fact complicates the data analysis. We, therefore, adopted a simplified approach, which, despite being an approximation, is able to reproduce well all spectra. Specifically, we fitted the spectra with the local structure of only the Ag-I site as the absorbing atom; the rationale is that the average of all scattering paths giving rise to the EXAFS signal relative to both Ag sites is equivalent to the sum of scattering paths for the more ordered Ag-I site. Clearly, we do not claim that the distances and Debye−Waller factors obtained in this way reproduce faithfully the full atomic distribution around Ag; however, the adopted approach is able to reproduce key features and trends in the local structural parameters and has the advantage of simplicity.
EXAFS scattering paths were calculated ab initio using FEFF, and the highest amplitude ones were included in the fit, which was performed in the ranges = 2-10.5 Å −1 and = 1-3.5 Å . The interatomic distances, an energy origin shift, and Debye−Waller factors were considered as fitting parameters; the many body amplitude reduction factor was fixed to the value obtained from a fit of a Ag metal reference spectrum (S0 2 = 0.834). We found that the Ag2S spectrum could be fitted with single scattering paths due to two S atoms at ~2.50 Å (Ag-S1), one S atom at ~3.06 Å (Ag-S2), seven Ag atoms at ~2.97 Å (Ag-Ag1), and two more Ag atoms at ~3.11 Å (Ag-Ag2). These contributions give rise to the first two peaks in the FT of the spectra reported in Figure 4A, in which the components are labelled as just described. There are nine Ag atoms at different distances in the cation coordination The Ag 2 S spectrum was fitted using the monoclinic structure described above. As mentioned, Ag has two structurally inequivalent sites and, especially for site II, there is a rather broad distribution of interatomic distances around Ag. This fact complicates the data analysis. We, therefore, adopted a simplified approach, which, despite being an approximation, is able to reproduce well all spectra. Specifically, we fitted the spectra with the local structure of only the Ag-I site as the absorbing atom; the rationale is that the average of all scattering paths giving rise to the EXAFS signal relative to both Ag sites is equivalent to the sum of scattering paths for the more ordered Ag-I site. Clearly, we do not claim that the distances and Debye−Waller factors obtained in this way reproduce faithfully the full atomic distribution around Ag; however, the adopted approach is able to reproduce key features and trends in the local structural parameters and has the advantage of simplicity.
EXAFS scattering paths were calculated ab initio using FEFF, and the highest amplitude ones were included in the fit, which was performed in the ranges k = 2-10.5 Å −1 and R = 1-3.5 Å. The interatomic distances, an energy origin shift, and Debye−Waller factors were considered as fitting parameters; the many body amplitude reduction factor was fixed to the value obtained from a fit of a Ag metal reference spectrum (S 0 2 = 0.834). We found that the Ag 2 S spectrum could be fitted with single scattering paths due to two S atoms at~2.50 Å (Ag-S1), one S atom at~3.06 Å (Ag-S2), seven Ag atoms at~2.97 Å (Ag-Ag1), and two more Ag atoms at~3.11 Å (Ag-Ag2). These contributions give rise to the first two peaks in the FT of the spectra reported in Figure 4A, in which the components are labelled as just described. There are nine Ag atoms at different distances in the cation coordination shell, and we found that slightly different groupings of nine Ag atoms (for example, six plus three instead of seven plus two) did not give any significant differences except small variations in the Debye−Waller factor; clearly, in this approach, we are approximating the real distribution with a sum of two Gaussian components with different coordination numbers. The best fit and the individual components are shown in Figure 4A. A model derived from the monoclinic structure of Ag 2 S was considered to fit the EXAFS spectra of the samples. As shown in Figure 2 (right panel), the Ag-II atom was replaced with Bi, and the corresponding theoretical scattering paths were generated and used in the fit. The spectra of the samples were fitted as for Ag 2 S with a modification regarding the Ag-Ag contributions. The cation contribution was thus considered as a linear combination of Ag-Ag and Ag-Bi scattering paths with coordination numbers 9y and 9(1 − y), respectively, and 0 ≤ y ≤ 1; in this way, the total coordination of cations around the central Ag atom in the alloy is 9, as in AgS 2 . We define normalized Ag-Ag and Ag-Bi coordination numbers as N Ag−Ag = y and N Ag−Bi = 1 − y, respectively. The best fit for sample 2 is shown in Figure 4B, including the individual components (note that Ag-Bi substitutes Ag-Ag2) and the fits for all samples, for all of which the R-factor was less than 0.01 are shown in Figure S3.
The most important numerical results of the fit are reported in Table 3. We note that the Ag-S bond length is in the range of 2.48-2.51 Å, which is compatible with a low S coordination number in the first shell as in monoclinic acanthite and not with a high coordination site as in cubic schapbachite or hexagonal matildite [31]. This is a further confirmation that the local structure has features similar to those of the monoclinic phase. Furthermore, as shown in the Supplementary Material, the Bi L III edge EXAFS spectra were fitted using the same model, albeit only in the first shell, and the results are reasonably compatible with the Ag K-edge results. Figure S4 reports the FT and fit for a Bi 2 S 3 reference sample and Figure S5 the FTs and fits for the samples. Table 3. Summary of quantitative analysis of Ag K-edge X-ray absorption fine structure in the extended energy range (EXAFS) spectra of AgBiS 2 samples. Uncertainties in the least significant figures are reported in brackets (italics). y is the Ag concentration in the cation coordination shell and α is the Cowley short-range order parameter.

Sample Code
R Ag-S1 σ 2 Ag−S1 R Ag-S3 σ 2 Ag−S3 R Ag-Ag σ 2 Ag−Ag R Ag-Ag σ 2 Ag−Ag R Ag-Bi σ 2 In Figure 5, we report the normalized Ag-Ag and Ag-Bi coordination numbers, as defined above, with the Ag concentration x. It is evident that N Ag−Ag is always greater than x and N Ag−Bi is always smaller than 1 − x, indicating that there is a tendency for cation ordering (i.e., clustering of similar cations). In order to quantify the degree of cation ordering, we calculated the Cowley [32,33] short-range order parameter, defined as The value of α can vary in the interval between −1 and 1; α = 0 indicates a random distribution (y = x), negative values indicate an excess of Ag-Ag correlations with respect to the random case (like-atom clustering), and positive values indicate cation ordering (preference for Ag-Bi correlations); Nanomaterials 2020, 10, 316 8 of 13 values of α are reported in Table 3. The value of α for all samples indicates a tendency for like-atom clustering, suggesting the presence of Ag-rich and Bi-rich volumes within the NPs. The indication from XAFS that the local structure is similar to that of a monoclinic phase raises a serious issue as the XRD patterns exhibit peaks that, despite considerable broadening, are close to those of matildite and schapbachite and have no correspondence to those of monoclinic acanthite. In fact, in Figure 6, we report the experimental XRD patterns for sample 2 (the others are very similar) and simulations for the three structures. In all cases, the domain shape is assumed to be spherical, with a lognormal distribution of diameters (D), and with a small but non-negligible microstrain component of line broadening. Here, and in the following modelling, the background was represented by a Chebyshev polynomial of degree 12 to account for incoherent scattering, noise, and a fraction of non-crystalline phase. The acanthite pattern is so far from the experimental one that it cannot be refined (in the simulated pattern of Figure 6A, mean domain size was set to <D> = 5 nm, with a distribution standard deviation (s.d.) of 1.5 nm). Instead, the schapbachite and matildite simulations are close to the experiment pattern. Refined mean domain size <D> is 2.9 and 3.1 nm, respectively, for the cubic and hexagonal phases, with s.d. = 1 nm in both.  The indication from XAFS that the local structure is similar to that of a monoclinic phase raises a serious issue as the XRD patterns exhibit peaks that, despite considerable broadening, are close to those of matildite and schapbachite and have no correspondence to those of monoclinic acanthite. In fact, in Figure 6, we report the experimental XRD patterns for sample 2 (the others are very similar) and simulations for the three structures. In all cases, the domain shape is assumed to be spherical, with a lognormal distribution of diameters (D), and with a small but non-negligible microstrain component of line broadening. Here, and in the following modelling, the background was represented by a Chebyshev polynomial of degree 12 to account for incoherent scattering, noise, and a fraction of non-crystalline phase. The acanthite pattern is so far from the experimental one that it cannot be refined (in the simulated pattern of Figure 6A, mean domain size was set to <D> = 5 nm, with a distribution standard deviation (s.d.) of 1.5 nm). Instead, the schapbachite and matildite simulations are close to the experiment pattern. Refined mean domain size <D> is 2.9 and 3.1 nm, respectively, for the cubic and hexagonal phases, with s.d. = 1 nm in both. The indication from XAFS that the local structure is similar to that of a monoclinic phase raises a serious issue as the XRD patterns exhibit peaks that, despite considerable broadening, are close to those of matildite and schapbachite and have no correspondence to those of monoclinic acanthite. In fact, in Figure 6, we report the experimental XRD patterns for sample 2 (the others are very similar) and simulations for the three structures. In all cases, the domain shape is assumed to be spherical, with a lognormal distribution of diameters (D), and with a small but non-negligible microstrain component of line broadening. Here, and in the following modelling, the background was represented by a Chebyshev polynomial of degree 12 to account for incoherent scattering, noise, and a fraction of non-crystalline phase. The acanthite pattern is so far from the experimental one that it cannot be refined (in the simulated pattern of Figure 6A, mean domain size was set to <D> = 5 nm, with a distribution standard deviation (s.d.) of 1.5 nm). Instead, the schapbachite and matildite simulations are close to the experiment pattern. Refined mean domain size <D> is 2.9 and 3.1 nm, respectively, for the cubic and hexagonal phases, with s.d. = 1 nm in both.
(A) 10   An apparent serious contradiction is, therefore, present. XAFS indicates a local structure around Ag cations similar to monoclinic acanthite; the detection of short (~2.48 Å ) Ag-S bonds is particularly significant as EXAFS is a reliable method to measure interatomic distances in the first shell. On the other hand, the positions of XRD diffraction peaks are also, undeniably, a reliable fingerprint of the presence of particular crystallographic phases, even for NPs characterized by large particle-sizeinduced peak broadening and strain.
Ab initio DFT and MD simulations provide original insight to solve this contradiction. Analysis of the DFT-relaxed structures of the model compounds (acanthite, Bi2S3, and S) revealed an accuracy in the determination of the (Ag, Bi)-S distances of 1% with respect to experimental values and will not be discussed here. Table 4 reports some data retrieved from the static DFT (formation energy) and DFT−MD (structure) modelling of candidate structures to be considered for the interpretation of the XAFS data. An image of the structure resulting from the DFT−MD simulations is reported in Figure S6. The results of the formation energy confirm the observation of Viñes et al. [6] that matildite is more stable than schapbachite. However, these values should be considered just from a qualitative point of view as they are relative to a 0 K calculation and to a specific equilibrium configuration. Due to the complex environments around the cations involving many sites, structural data were derived from the first peak of the partial radial pair distribution functions (PRPDFs) obtained from the MD runs. In particular, the first peaks in the Ag−S and Ag−Ag PRPDFs were fitted with two Gaussians and the results are presented in Table 4. An apparent serious contradiction is, therefore, present. XAFS indicates a local structure around Ag cations similar to monoclinic acanthite; the detection of short (~2.48 Å) Ag-S bonds is particularly significant as EXAFS is a reliable method to measure interatomic distances in the first shell. On the other hand, the positions of XRD diffraction peaks are also, undeniably, a reliable fingerprint of the presence of particular crystallographic phases, even for NPs characterized by large particle-size-induced peak broadening and strain.
Ab initio DFT and MD simulations provide original insight to solve this contradiction. Analysis of the DFT-relaxed structures of the model compounds (acanthite, Bi 2 S 3 , and S) revealed an accuracy in the determination of the (Ag, Bi)-S distances of 1% with respect to experimental values and will not be discussed here. Table 4 reports some data retrieved from the static DFT (formation energy) and DFT−MD (structure) modelling of candidate structures to be considered for the interpretation of the XAFS data. An image of the structure resulting from the DFT−MD simulations is reported in Figure S6. The results of the formation energy confirm the observation of Viñes et al. [6] that matildite is more stable than schapbachite. However, these values should be considered just from a qualitative point of view as they are relative to a 0 K calculation and to a specific equilibrium configuration. Due to the complex environments around the cations involving many sites, structural data were derived from the first peak of the partial radial pair distribution functions (PRPDFs) obtained from the MD runs. In particular, the first peaks in the Ag−S and Ag−Ag PRPDFs were fitted with two Gaussians and the results are presented in Table 4. It must be noted that in the random matildite structure, "short" Ag-S bonds appear. This is illustrated in Figure 7 in which the PRPDFs of matildite and random matildite are compared. EXAFS spectra were calculated based on the MD simulations, as described in the previous section; in Figure 8, we report the comparison between the experimental spectrum of sample 2 and the simulations for the four candidate structures: Bi-substituted acanthite (AcB), matildite (Mat), matildite with random cation occupation (MaR), and schapbachite (Sch). It is clear that the only candidate structure for which there is a good comparison in terms of overall lineshape is the random matildite one. As a final check, atomic coordinates from the DFT-MD simulations were used to simulate the XRD patterns, and the result is reported in Figure 9. Here, the hexagonal structure of Ag 12 Bi 12 S 24 random matildite was mapped in the P1 space group, and the lattice parameters of the 2a, 2b, c supercell were refined starting from values for the standard hexagonal phase (a = b c, α = β = 90 • , γ = 120 • ). The best fit shows a visible improvement with respect to the standard matildite of Figure 6C, with the statistical quality index Rwp decreasing from 3.05 to 2.55, and a corresponding improvement in the goodness of fit (GoF) from 2.01 to 1.72 [34]. Unit cell parameters change slightly from the starting hexagonal values to: a = 0.8119(5) nm, b = 0.8127 (5) Figure 6 and TEM observations (cfr. Figure 1). The inhomogeneous strain is small but measurable (~0.0003), and is likely due to the surface relaxation commonly observed in nanocrystals [22,35]. The fact that random matildite is the preferential structure in the synthesized AgBiS 2 nanocrystals, as revealed by our results, may have consequences in key properties relevant to the photovoltaic performance of this material, such as the electronic structure and the energy band gap. Revision of previously reported hybrid DFT calculations (which were done with the pure matildite structure) [6] would bring further insight to evaluate the impact of the found distortions. Although dramatic changes are not expected, such variations could explain the moderate mismatch between theoretical and experimental data.
Summarizing, DFT-MD simulations show that a random distribution of Ag and Bi in the cation sites of the hexagonal matildite structure gives rise to "short" Ag-S bonds detected experimentally by EXAFS; at the same time, the lattice remains hexagonal, providing a good comparison to the XRD data. The EXAFS analysis indicated that the distribution of cations around Ag is not actually random, but a preference for Ag-Ag correlations is present. However, more than 10% Bi is always found in the local cation environment of Ag, enough to induce strong local distortions in the Ag-S, producing short bond lengths that are not present in the bulk matildite structure. It is possible that these local distortions are related to electronic defects that affect the photovoltaic properties, a subject worthy of further investigation. Finally, these results illustrate the power of combined experimental-simulation studies of complex structures found in nanostructured materials.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: Ag K-edge XANES spectra of the samples and Ag2S. Figure S2: Comparison of the Ag K-edge XANES spectrum of sample 2 with FDMNES simulations for schapbachite (cubic), acanthite (monoclinic), and matildite (hexagonal) bulk phases. Figure S3: Best fits of the Ag K-edge EXAFS spectra for all samples. Figure S4: FT of Bi LIII edge EXAFS spectrum of a Bi2S3 reference sample with its best fit. Figure S5: FT of Bi LIII edge EXAFS spectra of the samples with their best fits. Figure S6: Image of the structure resulting from the MD-DFT simulations. Table S1: Summary of quantitative analysis of Bi LIII-edge EXAFS spectra of the samples. Hexagonal P1 "random" matildite Ag12Bi12S24 Figure 9. Comparison between the experimental XRD pattern of sample 2 (red points) and the simulation based on atomic coordinates obtained from the density functional theory (DFT)-MD simulations (blue line), including an amorphous-like background (grey). The difference between the experimental data and model (residual) is shown below (black line).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2079-4991/10/2/316/s1, Figure S1: Ag K-edge XANES spectra of the samples and Ag 2 S. Figure S2: Comparison of the Ag K-edge XANES spectrum of sample 2 with FDMNES simulations for schapbachite (cubic), acanthite (monoclinic), and matildite (hexagonal) bulk phases. Figure S3: Best fits of the Ag K-edge EXAFS spectra for all samples. Figure S4: FT of Bi L III edge EXAFS spectrum of a Bi 2 S 3 reference sample with its best fit. Figure S5: FT of Bi L III edge EXAFS spectra of the samples with their best fits. Figure S6: Image of the structure resulting from the MD-DFT simulations. Table S1: Summary of quantitative analysis of Bi L III -edge EXAFS spectra of the samples.