rmc-discord: reverse Monte Carlo refinement of diffuse scattering and correlated disorder from single crystals

Morgan, Z.J.; Zhou, H.D.; Chakoumakos, B.C.; Ye, F.

doi:10.1107/S1600576721010141

computer programs

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 54| Part 6| December 2021| Pages 1867-1885

https://doi.org/10.1107/S1600576721010141

Open

access

rmc-discord: reverse Monte Carlo refinement of diffuse scattering and correlated disorder from single crystals

Zachary J. Morgan,^a ^* Haidong D. Zhou,^b Bryan C. Chakoumakos ^a and Feng Ye ^a

^aNeutron Scattering Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA, and ^bDepartment of Physics and Astronomy, University of Tennessee, Knoxville, TN 37996, USA
^*Correspondence e-mail: morganzj@ornl.gov

Edited by T. J. Sato, Tohoku University, Japan (Received 10 June 2021; accepted 30 September 2021; online 23 November 2021)

A user-friendly program has been developed to analyze diffuse scattering from single crystals with the reverse Monte Carlo method. The approach allows for refinement of correlated disorder from atomistic supercells with magnetic or structural (occupational and/or displacive) disorder. The program is written in Python and optimized for performance and efficiency. Refinements of two user cases obtained with legacy neutron-scattering data demonstrate the effectiveness of the approach and the developed program. It is shown with bixbyite, a naturally occurring magnetic mineral, that the calculated three-dimensional spin-pair correlations are resolved with finer real-space resolution compared with the pair distribution function calculated directly from the reciprocal-space pattern. With the triangular lattice Ba₃Co₂O₆(CO₃)_0.7, refinements of occupational and displacive disorder are combined to extract the one-dimensional intra-chain correlations of carbonate molecules that move toward neighboring vacant sites to accommodate strain induced by electrostatic interactions. The program is packaged with a graphical user interface and extensible to serve the needs of single-crystal diffractometer instruments that collect diffuse-scattering data.

Keywords: diffuse scattering; correlated disorder; single crystals; rmc-discord; Monte Carlo refinement; supercells.

1. Introduction

Crystalline materials with disorder giving rise to unusual properties can be studied with diffuse-scattering measurements (Frey, 1995 ). These behaviors involve a wide range of phenomena from magnetically to structurally disordered crystals. Quantifying and understanding the diffuse scattering is crucial to developing property–structure relationships of the disorder in crystalline materials. Unlike conventional crystallography and diffraction in which the `average' structure is obtained, diffuse-scattering studies reveal the local `deviations' away from the average (Welberry & Butler, 1995 ). Often, the disorder is locally correlated, involving one or more of the three general disorder types. These include correlations of magnetic moments, site occupancies and/or atomic displacements.

Correlated disorder is essentially ubiquitous in all functional materials, and the nature of that disorder is closely related to its particular property functionalities (Keen & Goodwin, 2015 ). Diffuse scattering can be measured in both powder and single-crystal diffraction experiments; however, characterization of the three-dimensional reciprocal space is achievable by instrumentation advancements that allow complete volumes of single-crystal diffraction data to be collected (Welberry & Goossens, 2014 ). A recent technique has emerged to obtain the three-dimensional `difference' pair distribution function (3D-ΔPDF) from single-crystal diffuse-scattering data. This autocorrelation function describing the deviations away from the average structure can be obtained for occupational/displacive disorder (Weber & Simonov, 2012 ) and magnetic disorder (Roth et al., 2018 ) through the Fourier transform. The refinement program Yell (Simonov et al., 2014 ) is capable of refining the disordered structure from the full-volume diffuse-scattering data set interpreting with the 3D-ΔPDF.

Although powerful, the 3D-ΔPDF requires special care to separate the Bragg and diffuse scattering from the total scattering signal. For example, Bragg reflections are removed and the missing data are filled back in by interpolation to resemble the surrounding diffuse intensity data. However, this process may introduce artefacts in correlations with larger separation vectors. More advanced methods to deal with Bragg peaks include K-space algorithmic reconstruction (KAREN), which uses statistical methods to separate the Bragg from the diffuse scattering (Weng et al., 2020 ). Another challenge with the 3D-ΔPDF is the requirement of a complete volume of diffraction data. Often, it is impossible to obtain complete reciprocal-space volumes so data symmetrization is utilized (Michels-Clark et al., 2016 ). In many instances, the number of available symmetry operators is reduced by the disorder itself. Similarly, sample environments on instruments limit the detector coverage, further impinging the volume mapping. In these cases, alternative approaches to model-free methods are considered to obtain the three-dimensional correlations.

Atomistic Monte Carlo methods are one such technique, where discrete crystals are constructed with atoms located at their average positions within a unit cell (Proffen & Welberry, 1998b ). A supercell is generated by extending the unit cell in three-dimensional space, and disorder is introduced by decorating atoms with either magnetic moments, vacancies or atomic displacements. Changes are made by rotating moments, replacing atoms with vacancies and displacing atoms until the disordered supercell configuration agrees with experiment through the calculated diffraction intensity. This relatively straightforward approach is a robust method for quantifying the disorder from experimental single-crystal diffraction data, allowing for direct calculation of three-dimensional pair correlations.

Using forward Monte Carlo (FMC) modeling, a Hamiltonian is constructed that describes the system energy and involves coupling constants between the neighboring atoms (Welberry & Butler, 1994 ). For example, a general Heisenberg model (Ashcroft & Mermin, 1976 ) is used in magnetically disordered systems where each spin vector pair has an (isotropic or anisotropic) exchange coefficient describing the interaction strength between bonds. To simplify the Hamiltonian, often pairs up to a finite number of nearest neighbors are considered. In addition, long-range dipole interactions are greatly simplified assuming Ising spin constraints along easy directions such as in spin ice (Bramwell & Gingras, 2001 ). Similar simple energy models are used to capture occupational disorder using binary random variables and displacive disorder using harmonic potentials (Proffen, 2000 ). As these simplified models generate many of the diffuse-scattering features (Welberry, 2004 ), additional effort is required for more complicated crystal structures where the models become more difficult to adapt and the disorder is not well known or established.

The other atomistic method is reverse Monte Carlo (RMC). RMC is characterized not by an energy function with a model but rather by a goodness-of-fit parameter used to `fit' the supercell disorder to experimentally obtained scattering data (McGreevy & Pusztai, 1988 ; McGreevy, 2001 ). RMC does not use a Hamiltonian and is a good tool for systems with unknown disorder. This is generally the case for most complicated material systems that are difficult to model. Several RMC programs exist for analysis of diffuse scattering, including RMCProfile for polycrystalline and powder diffraction (Tucker et al., 2007 ), Spinvert for powder diffraction of magnetic materials (Paddison et al., 2013 ), and DISCUS for powder and single-crystal diffraction (Proffen & Neder, 1997 ).

The analysis of single-crystal diffuse scattering from diverse materials systems with complicated behaviors over different conditions is a serious challenge even with RMC (Nield et al., 1995 ). The resources required to perform refinements are computationally intensive since the diffuse-scattering pattern has to be recalculated for each proposed move of the RMC refinement. A user-friendly and efficient single-crystal RMC program that integrates into data-reduction workflows is highly desired at X-ray and neutron-scattering user facilities, since it avoids construction of complicated models and works on partial volumes of the reciprocal space. The new program rmc-discord is introduced to address this challenge for each of the chemical, structural and magnetic disorder types, optimizing algorithms for performance and providing a user-friendly experience.

2. Reverse Monte Carlo approach

The RMC method uses the Metropolis algorithm, similarly to FMC; details are available elsewhere (Proffen & Welberry, 1997 , 1998a ; Welberry & Proffen, 1998 ). For any candidate move that gives an improved fit χ², defined here as

$[\chi^{2} = \sum_{\bf{Q}}\bigg[{{I_{\rm {calc}}({\bf Q})-{I}_ {\rm {expt}}(\bf{Q})} \over {\sigma_{\rm {expt}}(\bf{Q})}}\bigg]^{ 2}, \eqno (1)]$

the disturbance is always accepted. Here I_calc(Q) and I_expt(Q) are the calculated and experimental intensity, respectively. The experimental error σ_expt(Q) is obtained from the data normalization (Michels-Clark et al., 2016). If instead the fit worsens, the move is either accepted or rejected with some probability according to an acceptance ratio that is a function of the current system temperature and magnitude of the change in fit. As the refinement progresses and overcomes local minima, fewer and fewer candidates that worsen the fit are accepted as it converges to a global minimum.

2.1. Single-crystal diffuse-scattering intensity calculations

A major advancement in the calculation of single-crystal diffuse-scattering patterns from atomistic models utilizes the fast Fourier transform (FFT), giving a speedup of at least a factor of one hundred over traditional algorithms in the program Scatty (Paddison, 2019 ). Although an RMC implementation of this algorithm is tempting to achieve faster refinements, further improvements can be made recognizing that the diffuse-scattering pattern is calculated repeatedly throughout the refinement. By storing intermediate results in memory, it is possible to update the diffuse-scattering pattern with order corresponding to the chosen sizes of the supercell and reciprocal-space volume data set. Therefore, the recalculation of a diffuse-scattering pattern after a given move in the refinement can be achieved in nearly linear complexity. The Scatty algorithm using FFT is summarized below for magnetic and structural disorder types (Paddison, 2019), and then we discuss the modification of storing intermediate results in memory.

The underlying Scatty approach is to perform the FFT of the disordered magnetic or structural parameters over the discrete average position vectors of the supercell lattice points for each unique atom site labeled j = 0, 1,…. Through vector addition, the average location of an atom in a crystal can be defined in terms of a lattice vector R and its basis vector r_j. For a supercell with finite size N₁ × N₂ × N₃, these positions can be written in terms of the direct lattice vectors a, b and c, and integer cell coordinates (R₁, R₂, R₃):

$[{\bf R} = R_{1}{\bf a}+R_{2}{\bf b}+R_{3} {\bf c},\quad\quad R_{i} = 0,1,\ldots,N_{i}-1, \eqno (2)]$

which produces the uniform grid spacing of the discrete Fourier transform (DFT) with (N₁N₂N₃)² operations (Briggs & Henson, 1995 ). The FFT algorithm efficiently evaluates this using of the order of N₁N₂N₃log₂(N₁N₂N₃) operations. The average position vector r_j of a unique site j within the unit cell is written in terms of the direct lattice vectors and its fractional coordinates u_j, v_j and w_j:

$[{\bf r}_{j} = u_{j}{\bf a}+v_{j}{\bf b}+w_{j}{\bf c}. \eqno (3)]$

The true position of any atom in the supercell is the sum of its average position and its displacement u_j.

Given a real-space field g_j(R) that captures the disorder among magnetic moments, site occupancies and/or atomic displacements at sites j within a cell of the supercell R, its DFT is defined by

$[g_{j}({\bf k}) = \textstyle\sum\limits_{{\bf R}}g_{j}({\bf R})\exp({i {\bf k}\cdot{\bf R}}). \eqno (4)]$

Here, g_j(k) is the sequence of Fourier coefficients in the spatial frequency domain k. This uniformly spaced wavevector is given by

$[{{{\bf k}} \over {2\pi}} = {{k_{1}} \over {N_{1}}}{\bf a}^{ \ast}+{{k_{2}} \over {N_{2}}}{\bf b}^{\ast}+{{k_{3}} \over {N_{3}}}{\bf c}^{\ast},\quad\quad k_{i} = 0,1,\ldots,N_{i}-1, \eqno (5)]$

where a*, b* and c* are the primitive reciprocal lattice vectors, and k₁, k₂ and k₃ are grid points. This wavevector k is defined within the reciprocal lattice unit cell with resolution determined by the supercell size. To obtain the reciprocal-space wavevector, the reciprocal lattice vector K defined in terms of integer Miller indices h, k and l,

$[{{{\bf K}} \over {2\pi}} = h{\bf a}^{\ast}+k{\bf b}^{\ast}+l {\bf c}^{\ast}, \eqno (6)]$

is added to the uniformly spaced wavevector k.

The experimental reciprocal-space data set is composed of equally spaced bin sizes and can be written as

$[{{{\bf Q}} \over {2\pi}} = Q_{1}{\bf a}^{\ast}+Q_{2}{\bf b}^{ \ast}+Q_{3}{\bf c}^{\ast}, \eqno (7)]$

where Q₁, Q₂ and Q₃ are of size n₁ × n₂ × n₃, respectively. Although the supercell size gives the maximum reciprocal-space resolution from k defined in equation (5), a finer (or coarser) grid of reciprocal-space points Q corresponding to experimentally obtained scattering data can be implemented by using interpolation (or resampling). First, the phase factor in the scattering equations for both magnetic and structural scattering can be written as $[\exp[{i{\bf Q}\cdot({\bf R}+{\bf r}_{j}+{\bf u}_{j})}] ]$ = $[\exp({i{\bf Q}\cdot{\bf R}})\exp({i{\bf Q}\cdot{\bf r}_{ j}})\exp({i{\bf Q}\cdot{\bf u}_{j}})]$ . Recognizing the identity exp(iK · R) = 1, the term exp(iQ · R) ≃ exp[i(K + k) · R] = exp(ik · R). The resulting phase factor $[\exp({i{\bf k}\cdot{\bf R}})]$ × $[\exp({i{\bf Q}\cdot{\bf r}_{j }})\exp({i{\bf Q}\cdot{\bf u}_{j}})]$ allows for the DFT to be incorporated into the diffuse-scattering calculations that are efficiently evaluated using the FFT algorithm.

The procedure for calculating the diffuse-scattering intensities from a supercell of moments located at their average position (u_j = 0) is summarized below. For magnetic scattering, the DFT of each magnetic moment vector M_j(R) is computed according to

$[{\bf M}_{j}({\bf k}) = \textstyle\sum\limits_{{\bf R}}{\bf M}_{j}({\bf R})\exp({i{\bf k}\cdot{\bf R}}). \eqno (8)]$

Taking the transformed moments M_j(k) over each site, the magnetic structure-factor vector F(Q) is calculated by summing over each unique site j within the unit cell:

$[{\bf F}({\bf Q}) = \textstyle\sum\limits_{j}f_{j}(Q){\bf M}_{j}({\bf k})\exp({i{\bf Q}\cdot{\bf r}_{j}}), \eqno (9)]$

where f_j(Q) is the magnetic form factor of the ion. The structure factor is then projected into the scattering plane defined by the unit vector $[{\hat{\bf Q}}]$ of the reciprocal-space wavevector Q using

$[{\bf F}_{\perp} = {\hat{\bf Q}}\times{\bf F}\times {\hat{\bf Q}}. \eqno (10)]$

Finally, the magnetic intensity is calculated according to

$[I({\bf Q}) = {{C} \over {N}}\big{\langle}|{\bf F}_{\perp}({\bf Q})|^{2}\big{\rangle}, \eqno (11)]$

where C is a constant (0.07265 b) and N is the total number of atoms in the supercell. The angle brackets indicate that multiple supercell refinements can be averaged together (Paddison, 2019) to reduce the noise of a recalculated pattern. For each individual refinement, however, the refined intensity is calculated using only the refined supercell.

A similar approach is followed for occupational and displacive disorder. Assuming that each site can be occupied by either an atom or a vacancy, a binary site-occupancy parameter δ_j(R) can be defined according to

$[\delta_{j}({\bf R}) = \cases{0 &\rm unoccupied\cr 1 &\rm occupied} \eqno (12)]$

which can be incorporated into the structure-factor calculation. A site that is occupied may also take on some atomic displacement u_j(R) away from its average position. The structure factor accounting for both occupancy and displacement is given by

$[F({\bf Q}) = \textstyle\sum\limits_{j}\sum\limits_{{\bf R}}\delta_{j}({\bf R})b_{j} \exp({i{\bf k}\cdot{\bf R}})\exp({i{\bf Q}\cdot{\bf r}_{j }})\exp[{i{\bf Q}\cdot{\bf u}_{j}({\bf R})}], \eqno (13)]$

where b_j is the neutron-scattering length of atom j. For X-rays, the scattering length is replaced by form factor f_j(Q). Following the Scatty approach, it is necessary to take the Taylor expansion of $[\exp[{i{\bf Q}\cdot{\bf u}_{j}({\bf R})}]]$ truncated to order n about Q · u_j(R) = 0. A second-order expansion is typically sufficient, but higher orders are necessary to capture some diffuse features (Butler & Welberry, 1993 ). The trinomial expansion for the powers of the dot product in Cartesian coordinates for the mth term of the Taylor series can be written in a condensed notation as $[[{\bf Q}\cdot{\bf u}_{j}({\bf R})]^{m}]$ = $[[Q_{1}u_{1}^{j}({\bf R})+Q_{2}u_{2}^{j}({\bf R})+Q_{3}u_{3}^{j}({\bf R})] ^{m}]$ = $[\textstyle\sum_{|\alpha| = m}Q^{\alpha_{l}}U_{j}^{\alpha_{l}}({\bf R})]$ , where α indicates three-dimensional multi-index notation, a three-tuple of non-negative integers α = (α₁, α₂, α₃), and |α| = α₁ + α₂ + α₃. The Taylor expansion products Q^α and $[U_{j}^{\alpha}]$ are defined by

$[ Q^{\alpha} = \textstyle\prod\limits_{l = 1}^{3}Q_{l}^{\alpha_{l}} {\quad \rm and \quad} U_{j}^{\alpha}({\bf R}) = \prod\limits_{l = 1}^{3}{u}_{l}^{j}({\bf R})^{\alpha_{l}}. \eqno (14)]$

By definition, $[Q^{\alpha} = U_{j}^{\alpha}({\bf R})\equiv 1]$ when m = 0.

The occupancy parameter and displacement products can be combined to give a new structural parameter. To enforce the overall composition c_j for each atom j, the relative occupancy parameter is defined as

$[a_{j}({\bf R}) = {{\delta_{j}({\bf R})} \over {c_{j}}}-1, \eqno (15)]$

which takes the values a_j(R) = −1 and a_j(R) = 1/c_j − 1 when it is unoccupied and occupied, respectively. Multiplying the site-occupancy parameter by the displacement products gives a new structural parameter,

$[V_{j}^{\alpha}({\bf R}) = c_{j}[1+a_{j}({\bf R})]U_{j}^{\alpha}({\bf R}), \eqno (16)]$

which can be converted to the wavevector domain by DFT,

$[V_{j}^{\alpha}({\bf k}) = \textstyle\sum\limits_{{\bf R}}V_{j}^{\alpha}( {\bf R})\exp({i{\bf k}\cdot{\bf R}}). \eqno (17)]$

The structure factor in equation (13) can then be approximated by

$[F({\bf Q}) = \sum\limits_{j}b_{j}\exp({i{\bf Q}\cdot{ {\bf r}_{j}}})\sum\limits_{m = 0}^{n}i^{m}\sum_{|\alpha| = m}{{1} \over {\alpha!}}V_{j} ^{\alpha}({\bf k})Q^{\alpha}, \eqno (18)]$

where α! = α₁!α₂!α₃!. Lastly, the diffuse-scattering intensity is determined by

$[I({\bf Q}) = {{1} \over {N}}\Big\langle|F({\bf Q})-\langle F({\bf Q})\rangle|^{2}\Big\rangle, \eqno (19)]$

where 〈F(Q)〉 is the Bragg structure factor (Frey, 1995; Paddison, 2019) and the outer angle brackets indicate averaging over multiple refinements. Again, this averaging over supercells is only carried out for the recalculated pattern. In the limit of no displacements, u_j(R) = 0 and only the m = 0 term remains in the Taylor expansion products such that the two right-most sums in equation (18) simplify to δ_j(R). Similarly, in the limit of no occupational disorder, δ_j(R) = c_j[1 + a_j(R)] = 1, and the parameter simplifies to $[V_{j}^{\alpha} = U_{j}^{\alpha}]$ .

2.2. Storing of intermediate results in memory

For each RMC move, one atom in the supercell is changed, requiring the complete recalculation of the scattering pattern. A straightforward implementation for recalculating the intensity from the structure factor is to recompute the FFT in N₁N₂N₃log₂(N₁N₂N₃) operations for parameters affected by the change at location (R₁, R₂, R₃, j). Next, recalculate the structure factor by multiplying the FFT result with phase and form factors summing over each site j. Finally, obtain the new intensity from the magnitude of the structure factor. However, the intermediate results for each one of these operations can be computed before the refinement and stored into memory. Significant speedups in recalculating the intensity can be made by taking advantage of these intermediate results, reusing and updating them as needed according to whether a move is accepted or rejected. Using the definition of the DFT, the transformed results of a given field computed initially by FFT can be updated with order corresponding to the supercell size. The structure factors can be updated using a similar procedure with order corresponding to the size of the experimental data set. The two procedures for updating the DFT result and structure factor are outlined below.

Beginning with the DFT update procedure, every proposed move is accompanied by a change in only one atom of the supercell with coordinates (R₁, R₂, R₃, j). All other atoms remain unchanged. Considering the field represented by the arbitrary function g_j(R) = g_j(R₁, R₂, R₃) in equation (4), its DFT g_j(k) = g_j(k₁, k₂, k₃) is affected for all k_i = 0, 1,…, N_i − 1. Using the FFT, recalculation of equation (4) reduces the order of operations from (N₁N₂N₃)² to N₁N₂N₃log₂(N₁N₂ N₃). However, intermediate results from the DFT calculation itself can be stored in memory to update g_j(k) on the order of N₁N₂N₃.

First, the exponential factors exp(ik · R) are calculated and stored in an array exp_fact[k1,k2,k3,R1,R2,R3] of size (N₁N₂N₃)² at the beginning of the refinement. For a chosen atom with coordinates (R₁, R₂, R₃, j), the N₁N₂N₃ values of g_j(k₁, k₂, k₃) are copied into a smaller `original' array for site j. A single `candidate' value is generated at (R₁, R₂, R₃, j) for g_j(R₁, R₂, R₃). Next, a Fourier components candidate array of size N₁N₂N₃ corresponding to site j is computed by adding to the original array the product of the exponential factors with the difference between the candidate and original value of g_j(R₁, R₂, R₃). This is carried out for all k_i = 0, 1,…, N_i − 1. Pseudocode of this procedure is shown in Fig. 1 for updating the Fourier components of the magnetic moment vector M_j(k) and the structural parameter $[V_{j}^{\alpha}({\bf k})]$ . Once the candidate array is computed, it is then used to update the structure factor.

Figure 1
Pseudocode for updating the magnetic and structure Fourier components using the intermediate results. For a given change in the supercell at location (R₁, R₂, R₃, j), a nested for loop of size N₁ × N₂ × N₃ is needed to update the Fourier components. By comparison, recalculation of the DFT by FFT requires on the order of N₁N₂N₃log₂(N₁N₂N₃) operations.

A similar procedure is utilized to update the structure factors. At the beginning of the refinement, two copies of the structure factor are made. For every proposed update, one is used to calculate the candidate structure factor while the other is used to store its original values. The structure-factor parameters that never change are also stored in an array fact[q1,q2,q3,j]. Here, q₁, q₂ and q₃ are integer indices corresponding to Q₁, Q₂ and Q₃ in equation (7). For magnetic disorder, the factors are $[f_{j}(Q)\exp({i{\bf Q}\cdot{\bf r}_{j}})]$ , and for occupational/displacive disorder, the factors are $[b_{j}\exp({i{\bf Q}\cdot{\bf r}_{j}})]$ . The products of these factors with all of the corresponding Fourier components are saved in a product array: these products include one for each of the three vector components of the magnetic moment vector M_j(k) or one for the structural parameter $[V_{j}^{\alpha}({\bf k})]$ . These arrays have size n₁n₂n₃ multiplied by the number of unit-cell atoms. For each update of a Fourier component corresponding to site j, the n₁n₂n₃ values in the product array for site j are copied into an original array before a candidate array is computed by multiplying the Fourier components with the factors array. The candidate structure-factor array is computed for all Q_i by adding to the original array the difference between the candidate and original product arrays. This removes the need to sum over all j as in equations (9) and (18). All structure-factor updates can be accomplished in the order of reciprocal-space size. This update procedure is illustrated in Fig. 2, which displays the pseudocode for the magnetic structure-factor vector F(Q) and nonmagnetic structure factor F(Q).

Figure 2
Pseudocode for updating the structure factors using the intermediate results and a nested for loop of size n₁ × n₂ × n₃. Multiplying the updated Fourier components by precomputed fixed factors into a product array, the structure factor is updated by simply adding the difference between the candidate and original arrays. A mapping function is used to account for nearest-neighbor interpolation due to differences in reciprocal-space and Fourier-space resolution.

After the candidate structure factors are calculated, the intensities are calculated using equations (11) and (19). The typical RMC refinement continues. If the proposed move is accepted, the candidate Fourier component array is copied into the larger array of components. The candidate products array is also copied into its corresponding larger array. If the move is rejected, no copies are made into the larger arrays and a new move is proposed. This procedure of saving intermediate results allows for fast recalculation of the intensity. It also benefits from multi-threaded parallelism.

2.3. Implementation and libraries

The program is written in Python (https://www.python.org/), allowing for efficient development and interfacing with many widely available software libraries including NumPy (https://numpy.org/) for array manipulation, SciPy (https://www.scipy.org/) for advanced scientific functions like kD trees, and Matplotlib (https://matplotlib.org/) for plotting. Performance-critical functions and subroutines are written in Cython (https://cython.org/), providing extensions with the performance of compiled C code with support for parallelism through OpenMP (https://www.openmp.org/) with parallel `for loops'. To read crystallographic information framework (CIF) files, PyCifRW (Hester, 2006 ) is utilized, allowing for the creation of crystals. In addition, it allows for the exportation of supercell data for visualization in external viewers like VESTA (Momma & Izumi, 2008 ). Diffuse-scattering patterns saved in the NeXus data format (Könnecke et al., 2015 ) are read using the Python package nexusformat (https://github.com/nexpy/nexusformat). A graphical user interface (GUI) is implemented using PyQt5 (https://riverbankcomputing.com/software/pyqt/). The package is currently available in Windows and Linux distributions through the Python Package Index (PyPI) (https://pypi.org/). A Mac OS/X is planned in a future release.

2.4. Capabilities and user requirements

The RMC program is designed to be simple yet robust to help facilitate the analysis of diffuse-scattering measurements. The only major requirements are a CIF file of the average structure of the crystal and a diffuse-scattering measurement data set for refinement. The program currently supports the reading of NeXus files that are built on top of the HDF5 file format. Alternative data readers will also be implemented and supported. After performing the refinement, the user can calculate and visualize correlations. All data sets can be exported for further analysis in external programs.

3. Scientific user cases

To illustrate the effectiveness of the RMC method in extracting the correlated disorder from more complicated systems, two user cases are presented below: one for magnetic disorder and a second for structural disorder. The first is from the mineral bixbyite and the other is from a triangular lattice system Ba₃Co₂O₆(CO₃)_0.7. The magnetic case is presented first since it has been previously analyzed using the magnetic 3D-ΔPDF method (Roth et al., 2018, 2019 ). The second nonmagnetic structural case has not been previously analyzed; the RMC program is one tool capable of extracting the correlated disorder. Both data sets were collected from the elastic diffuse-scattering spectrometer CORELLI (Ye et al., 2018 ).

3.1. Magnetic diffuse scattering of frustrated magnet bixbyite

The first user case is from the diffuse scattering of single-crystal bixbyite, (Fe, Mn)₂O₃, a natural mineral with magnetic iron and manganese ions (Roth et al., 2019). It is a cubic crystal with space group 206 $[Ia\bar{3}]$ and lattice parameter a = 9.409 Å (Pauling & Shappell, 1930 ). Shown in Fig. 3 is the polyhedral model of bixbyite and its network of (Mn, Fe)₂O₃ polyhedra. The magnetic Mn³⁺ and Fe³⁺ ions have nearly identical magnetic form factors (Lisher & Forsyth, 1971 ) and a similar magnitude of magnetic moments. The two metal sites in the structure are occupied by both manganese and iron: one is Fe rich and the other is Mn rich. The Fe-rich site resides at the origin (0, 0, 0), while the Mn-rich site is located at (0, 0.25, 0.285): that is, nearly (0, 0.25, 0.25). The structure can be thought of as a deviation from a higher-symmetry one with the position of the Mn-rich sites slightly altered. This reduced cubic structure has space group 221 $[Pm\bar{3}m]$ with a lattice parameter that is half of the original. Considering identical ions, the structure corresponds to a face-centered cubic crystal, which itself is a common frustrating magnetic lattice (Ramirez, 1994 ).

Figure 3
Bixbyite. (a) A polyhedral model of (Fe, Mn)₂O₃ and (b) a network of frustrated hexagons with nearest-neighbor Mn-rich sites surrounding central Fe-rich sites. Below its transition temperature T = 32.5 K, bixbyite is highly geometrically frustrated.

For this $[Pm\bar{3}m]$ structure, the nearest-neighbor pairs of Fe–Mn and Mn–Mn have identical distances. forming a cuboctahedral network of manganese pairs surrounding a central iron atom. The deviation of Mn-rich sites from the high-symmetry position creates a nearest-neighbor network consisting of a hexagonal arrangement of nearly coplanar manganese pairs surrounding an iron atom, as shown in Fig. 3(b). Bulk property measurements of field-cooled and zero-field-cooled magnetization suggest spin glass behavior with inverse susceptibility indicating strong antiferromagnetic interactions (Roth et al., 2018). The single-crystal diffuse scattering of bixbyite (Roth et al., 2018) is displayed in Fig. 4, which shows Bragg and non-Bragg layers in (a) and (b), respectively.

Figure 4
Single-crystal neutron scattering of bixbyite at (hkl) planes (a) l = 0 and (b) l = 0.5, obtained by subtracting high-temperature T = 300 K from low-temperature T = 7 K data. The complete reciprocal-space volume covers −10 to 10 in each h, k and l dimension using 501 × 501 × 501 bins. The diffuse-scattering features repeat every four reciprocal lattice units in h, k and l.

A strategy for fast refinement is to focus on a region of the diffuse scattering with minimal overlapping symmetry that covers a primitive segment of the repeating pattern. For bixbyite, the h, k and l range 0–4 satisfies this requirement. In general, a two-dimensional slice is not appropriate for refinements with scattering that show three-dimensional symmetry. Some liberty can be taken to restrict the range of a particular direction while maintaining a good refinement. It is found that the range of one dimension (along l) can be reduced by half in bixbyite. However, high and low intensities are both needed for refinement. For any volume not refined against, the program has no way of calculating goodness of fit in those missing regions so the intensity values are unconstrained there. Hence, the primitive features are the `minimum' requirement to obtain a reasonable refinement. A wide coverage with overlapping symmetry is preferred to obtain the best refinement.

Using the ranges described above, Fig. 5 shows a comparison between the experimental and refinement data. The Bragg peaks are removed and, unlike for the 3D-ΔPDF, the missing data do not have to be filled back in. Figs. 5(a) and 5(b) shows the comparison between the experiment and refined data, respectively, for a slice taken from the plane l = 0. Figs. 5(c) and 5(d) show the corresponding data at slice l = 0.5. Both sets of slices are in good agreement. Spin–spin correlations are calculated using the refined supercell structure by taking the dot product between spin vectors of every possible ion pair and averaging the results over common distances and ion pairings within a cutoff distance. With RMC refinement, the subtle differences in pair lengths can be easily resolved in the correlation calculations. As the magnetic 3D-ΔPDF has no underlying information on the average positions of ions within the crystal, it is generally more difficult to distinguish pairs that are closely spaced, especially like they are in bixbyite when considering first- and second-nearest-neighbor manganese pairs.

Figure 5
Diffuse-scattering intensity from bixbyite showing selected (hkl) planes of the [(a), (c)] cropped and rebinned experimental data compared with [(b), (d)] RMC refinement values. The region of interest has h and k ranging from 0 to 4 and an l range from 0 to 2 with a grid size of 0.08. A 4 × 4 × 4 supercell is used. The match between experiment and refinement is excellent. The comparison for plane l = 0 with strongest diffuse-scattering features is shown in (a) for the experiment data and (b) for the refinement values, while (c) and (d) show the corresponding data for plane l = 0.5. The RMC refinement captures all of the subtle features of the experimentally obtained pattern with the nuclear Bragg peaks punched out.

The calculated spin–spin correlations for the refinement are plotted in a variety of ways. Fig. 6(a) shows the spherically averaged correlations where the closest nearest-neighbor pairs have strong antiferromagnetic preference. The calculated 3D correlations in the plane z = 0 are displayed in (b) and compared with the magnetic 3D-ΔPDF in (c) using the preprocessed intensity data corresponding to Fig. 4 by removing Bragg peaks and interpolating the missing data. One advantage of RMC refinement is that the 3D correlations are easily mapped to the proposed frustrated structure. Fig. 6(d) displays the 3D spin-pair correlations corresponding to the nearest-neighbor hexagonal network shown in Fig. 3(b). The calculated spin correlations of Fig. 6(a) are labeled by their corresponding ion pair. It is clear that the six manganese atoms surrounding the central iron atom have antiferromagnetic preference.

Figure 6
Magnetic spin-pair correlations for bixbyite. (a) Spherically averaged radial spin correlations show antiferromagnetic first- and second-nearest-neighbor Fe–Mn (orange fill) and Mn–Mn (blue fill) pairs with the color indicating the bond type. (b) The corresponding 3D spin correlations in the xy plane with z = 0 reveal that the second-nearest-neighbor Fe–Mn antiferromagnetic (blue fill) correlations are stronger than the first. (c) The 3D correlations resolve the subtle differences in bond length between neighboring pairs (upper right) better than the magnetic 3D-ΔPDF (lower left). (d) A portion of the nearest-neighbor hexagonal network along a principal axis with manganese surrounding iron. Here, iron is placed at the origin (red fill).

The calculated correlations offer additional detail of pairs with nearly equal separation vectors. This is clear comparing the RMC correlations with the 3D-ΔPDF in Fig. 6(c). For example, the first and second Fe–Mn pairs have slightly different separation vectors deviating from the aforementioned $[Pm\bar{3}m]$ structure. From the refinement, the antiferromagnetic interactions between second-nearest-neighbor pairs are stronger than the first Fe–Mn pairs, the latter of which result in the frustrated hexagonal network. This is also observed in the radial spin correlations plotted in Fig. 6(a), which shows the difference in strength of the antiferromagnetic first and second Fe–Mn pairs. The nearest-neighbor Mn–Mn pairs show similar preference, where the second are stronger than the first. The frustration extends beyond the hexagonal network.

Ion pairs that share the same separation can be individually calculated from the RMC refinement. The first Fe–Fe pairs and third Mn–Mn pairs both share common separation vectors, each apart by half the length of the unit-cell edge. Both of these correlations are strongly ferromagnetic, as indicated in Fig. 6(a). In the context of the hexagonal network, the magnetic moments that reside at the center tend to align with the other central ions of nearby hexagons. Similarly, the magnetic moments along the vertices of the hexagons have alignment preference with the next closest corresponding ion in the neighboring hexagon. The ability to distinguish different ion types allows complicated frustrated structures to be analyzed through the construction of 3D networks, like the one shown in Fig. 6(d), to help guide the building models of the magnetic interactions.

3.2. Structural diffuse scattering of disordered Ba₃Co₂O₆(CO₃)_0.7

The next example demonstrates RMC refinement of structural disorder from the triangular lattice system Ba₃Co₂O₆(CO₃)_0.7. It is a hexagonal crystal with space group 174 $[P\bar{6}]$ and lattice parameters a = 9.683 Å and c = 9.518 Å (Boulahya et al., 2000 ), composed of CoO₆ octahedra and carbonate CO₃ molecular chains along its c axis (Iwasaki et al., 2009 ), as shown in Figs. 7(a) and 7(b). These chains are visualized in Fig. 7(c), where the average occupancy of the polyatomic ion CO₃²⁻ molecule is 0.7. The magnetic cobalt ions form honeycomb layers, which is one lattice that can be geometrically frustrated. A single layer is shown in Fig. 7(d).

Figure 7
A polyhedral model of Ba₃Co₂O₆(CO₃)_0.7 showing chains of CoO₆ and carbonate CO₃ in (a) the standard orientation and (b) along the c axis. (c) Isolated views of chains of polyatomic CO₃²⁻ ions with occupancy 0.7. Single-crystal measurements of inverse susceptibility reveal inter-layer ferromagnetic and intra-layer antiferromagnetic interactions that suggest a spin liquid candidate of Ising spins on the (d) cobalt honeycomb lattice.

The absence of long-range magnetic ordering at low temperatures suggests a spin liquid candidate (Igarashi et al., 2012 ); however, the diffuse-neutron-scattering pattern does not show any signs of magnetic correlation even at low temperature. The characteristic drop off in intensity at larger scattering vectors indicative of the magnetic form factor is not present. Instead, diffuse features are observed within different non-integer l planes, as shown in Fig. 8, that persist all the way to room temperature. The Bragg layer at (a) l = 0 shows no diffuse features, but such features are clearly observed in non-integer layers (b) l = 0.8 and (c) l = 2.8. The correlated disorder is therefore along the c axis of the crystal and its origin is structural rather than magnetic.

Figure 8
Single-crystal neutron scattering of Ba₃Co₂O₆(CO₃)_0.7 at (hkl) planes (a) l = 0.0, (b) l = 0.8 and (c) l = 2.8 obtained at T = 50 K. The complete reciprocal-space volume covers −10 to 10 in each h, k and l dimension using 501 × 501 × 501 bins. The corresponding single-crystal X-ray scattering data at (hkl) planes (d) l = 0.0, (e) l = 0.8 and (f) l = 2.8 obtained at T = 150 K from a Rigaku laboratory source do not exhibit any diffuse-scattering features, suggesting that the correlated disorder originates from the light elements which are more sensitive to neutrons than X-rays. For the neutron data, the diffuse scattering remains at temperatures above T = 150 K.

Single-crystal X-ray diffraction measurements with a Rigaku laboratory source reveal no diffuse scattering at these planes, as indicated in Figs. 8(d)–8(f). Because X-ray scattering length is proportional to atomic number (Seltzer, 1995 ), the absence of diffuse scattering in the X-ray data would suggest that the heavier elements are not the origin as the neutron-scattering lengths of oxygen and carbon are similar to those of barium and cobalt (Sears, 1992 ). Synchrotron measurements could help confirm this. Since the average occupancy of CO₃ is 0.7, occupational disorder of the molecule itself is one possible origin. Calculating the structure factor of CO₃ only, it is observed that the Bragg scattering intensity at l = 0 and l = 1 in Figs. 9(a) and 9(b) resembles that of the diffuse scattering at the l = 2.8 and l = 0.8 planes in Figs. 8(b) and 8(c), respectively. This indicates that the CO₃ groups are responsible for the disorder. Neutron-diffraction data collected using the TOPAZ instrument (Schultz et al., 2014 ) result in atomic displacement factors with prolate oxygen spheroids generally oriented along the chain, as indicated in Fig. 9(c), suggesting the presence of displacive disorder as well. For these reasons, a detailed analysis of CO₃ disorder is key to understanding the origin of the diffuse scattering in Ba₃Co₂O₆(CO₃)_0.7.

Figure 9
The calculated CO₃-only Bragg scattering intensity at (hkl) planes (a) l = 0 and (b) l = 1 resembles the diffuse neutron scattering at (hkl) planes l = 2.8 and l = 0.8, respectively, suggesting that the diffuse pattern originates from CO₃ disorder along the c axis. Additional neutron-scattering data of (c) atomic displacement ellipsoids (99% probability) for the carbon and oxygen atoms on the CO₃ molecule elongated along the c axis further hint at correlated displacements of the CO₃ molecules. The orientation variants are labeled as A, B and A′.

In each unit cell, the three CO₃ molecules take on one of two orientation variants that are ∼180° from one another. The pattern repeats as A–B–A′–A–B–A′–A–B–A′⋯, where A and A′ have the same orientation distinguished from B with opposite orientation. After RMC refinement, it is straightforward to separate calculated correlations, placing different variants at the origin to investigate the correlated disorder. For molecular-disorder refinement, CO₃ is kept as a rigid stoichiometric molecule and allowed to be either vacant on a site or present with finite displacement and rotation from its average position and orientation. To further speed up the refinement, the structure factor of the molecule is precalculated and stored as a prefactor. Considering only the displacement of the center of mass of the molecule (no rotation), a speedup corresponding to the number of atoms per molecule can be achieved. Such an approach is a general method for handling rigid assemblies of polyhedra.

Refinement of rigid CO₃ displacement and occupancy is performed to analyze the structural disorder. Starting with a random distribution of molecules with 70% occupied and no initial displacement, occupational disorder is refined for the first half of the refinement. The remaining molecules are allowed to displace away from their average positions in the second half. Although it is possible to refine both disorder types simultaneously, faster convergence is achieved when the occupational refinement is given more weight earlier in the refinement, since removing a molecule affects the structure factor much more than displacing it. For the displacive refinements, the molecules are allowed to displace anywhere uniformly within a sphere of radius 1.5 Å, which is about half the distance between molecules along the chain. The experimental data set is compared with the resulting refinement pattern obtained with RMC in Fig. 10: (a) and (b) compare the plane l = 0.8, while (c) and (d) compare the plane l = 2.8. This refinement provides a good fit of the experimental data. Incorporating rotation and distortion of the molecule would help improve the refined pattern. The refined features are also broader due to the finite size of the supercell. Extending its size will improve the sharpness of the refined intensity at the cost of speed. Correlations can then be obtained considering different variants as the origin.

Figure 10
Diffuse-scattering intensity from Ba₃Co₂O₆(CO₃)_0.7, showing selected (hkl) planes of the [(a), (c)] cropped and rebinned experimental data compared with [(b), (d)] RMC refinement values. The region of interest has an h range from −10 to 0, a k range from 0 to 10, an l range from 0 to 4 and a grid size of 0.08. A 12 × 12 × 12 supercell is used. The match between experiment and refinement is good. The comparison for plane l = 0.8 is shown in (a) for the experiment data and (b) for the refinement values, while (c) and (d) show the corresponding data for plane l = 2.8 with the strongest diffuse-scattering features. There are some additional cloudy features in the refined pattern due to the simplified model of rigid molecular displacements. The extra correlations generated by this constraint introduce these artefacts. The RMC refinement is also insensitive to the powder rings associated with the cryomagnet observed in the experimentally obtained pattern as they do not overlap with the main diffuse-scattering features.

Having obtained the occupation and position of the CO₃, various occupancy, displacement and vacant–displacement pair correlations can be obtained along the molecule chains. For occupational disorder, occupied and vacant molecules are mapped to a binary parameter that is either positive or negative unity, respectively (Welberry & Goossens, 2008 ). The calculated occupancy correlations along the chain in Fig. 11(a) are essentially featureless, ruling out any clustering of vacant sites where some characteristic distance by which the correlations become negative would be observed. Considering the average occupancy is 70%, two to three molecules on average are between every pair of vacant sites. The displacement correlations are analogous to magnetic spin-pair correlations or static displacement variables (Welberry & Goossens, 2014); inspecting them in Fig. 11(b), their modulation corresponds to a periodicity of four unit cells along the chain in the c-axis direction. This indicates that, within the first unit cell, neighboring occupied molecules tend to move away from one another.

Figure 11
Correlations along the CO₃ chains of Ba₃Co₂O₆(CO₃)_0.7 with the A′ molecule variant located at the origin: (a) occupancy, (b) displacement and (c) vacant–displacement pairs. The chain pair correlations follow the A′–A–B order illustrated in (d), the relative unit-cell fractional coordinates and the distance between variants. Representing the molecules as atoms, (e) an illustration of chain disorder is shown with displaced molecules (hatch fill) moving away from their average position (blue fill) toward nearby vacant sites (pink fill) along the chain.

The interactions between occupational and displacive disorder are quantified by identifying all vacant–occupied pairs. The dot product of the unit separation vector between pairs with the direction vector of the displaced molecule gives a measure of the vacant–displacement correlation. This metric is analogous to an atomic size effect parameter (Welberry, 1986 ; Butler et al., 1992 ) that quantifies the preference of an atom to move toward or away from a vacancy. In Fig. 11(c), the vacant–displacement correlations have the same periodicity as the displacement correlations but are out of phase by one unit cell. Illustrations of the intra-chain disorder are shown in Figs. 11(d) and 11(e). The spacings between the average positions of the molecules along the chain are nearly equal. As no clustering occurs along the chain, two to three molecules will be present between every pair of vacant sites, and those molecules displace away from their neighbors toward the closest missing site. This reduces the overall spacing between molecules producing the displacement and vacant–displacement modulations. These correlated displacements accommodate the electrostatic interactions between the polyatomic CO₃ ions and generate the diffuse-scattering pattern observed by neutrons. This behavior of Coulomb field relaxation is similar to another hexagonal system with channeled molecules. For the case of urea inclusion compounds that form a hexagonal network, rigid alkane groups within its channels exhibit orientational disorder (Welberry & Mayo, 1996 ) with a repeat distance different from the spacing of the urea itself, which generates diffuse-scattering features in non-integer l layers similar to Ba₃Co₂O₆(CO₃)_0.7.

4. Graphical user interface

An important aspect of the program is its GUI. Although scripts written to perform an RMC refinement are included, they are intended for development of new modules and testing purposes. From the perspective of the user, a simple GUI guides much of the workflow and reduces the learning barrier to obtaining meaningful information from a diffuse-scattering pattern. Four main steps achieve this: (1) building a supercell from the average structure, (2) processing the intensity data to make it suitable for refinement, (3) executing the refinement, and (4) calculating and interpreting the real-space correlations. These major steps are incorporated as tabs of the GUI window.

4.1. `Crystal' tab

The `Crystal' tab allows the user to read a CIF file corresponding to the average structure of the crystal being analyzed. In many cases, a user will have a CIF file from preliminary experiments on the sample. If not, several structural databases are available with published CIF files, including the Bilbao Crystallographic Server (Aroyo, Perez-Mato et al., 2006 , 2011 ; Aroyo, Kirov et al., 2006 ), Cambridge Structural Database (Groom et al., 2016 ), Inorganic Crystal Structure Database (Belsky et al., 2002 ), Crystallography Open Database (Gražulis et al., 2009 ) and others. This generates the unit-cell data including atomic sites, atom positions, atom types and site occupancies without requiring knowledge about the symmetry operators or lattice parameters used to define it. The user can specify the refinement type as either magnetic or nonmagnetic (structural). The size of the supercell can be specified at any time. In addition, the `Crystal' tab allows the unit-cell data to be freely adjusted, including switching atoms on or off, changing occupancies, and manipulating atom types. The program contains many magnetic form factors of ions and neutron-scattering lengths for different isotopes that can be simply selected. Future support for X-ray form factors for structural disorder is expected.

An example of the GUI for bixbyite is given in Fig. 12. In 12(a), the `Crystal' tab with finalized data used for refinement is shown. The overall unit-cell data are given on the left; a portion of the 32 iron and manganese magnetic ions are shown. The right shows the three atom sites read from the original CIF file. In this case the oxygen atom is unselected as it is not magnetic. Figs. 12(b)–12(d) illustrate this crystal-building process for bixbyite. Loading in the unit-cell data, the oxygen atoms are removed and the supercell is generated by repeating the modified unit cell according to its chosen size. The final supercell may be exported as another CIF file for visualization in an external viewer like VESTA.

Figure 12
(a) The GUI `Crystal' tab for building the supercell for refinement. Illustrations of (b) the original unit cell for bixbyite, (c) the unit cell with the oxygen atoms removed and (d) the final supercell visualized with VESTA.

4.2. `Intensity' tab

Once the supercell is generated, the experimental intensity is loaded from a supported data format. Current support is for the NeXus files output by many synchrotron and neutron sources. Future releases will have support for other data formats including NumPy arrays. Support for commercial single-crystal X-ray diffractometers is also possible. Once the data set is loaded, the user has options to clean up, crop and rebin the data for refinement. First, the data are rebinned to a coarser bin size using available predetermined sizes from a dropdown menu. The binning sizes available are factors of the original data-set size. The rebinning can be selected for each of the h, k and l dimensions. Next, the data set is cropped to a smaller region of interest by selecting the range for h, k and l.

Bragg peaks may be removed using a punch method. The user selects the cell centering with the default set from the space group within the CIF file. A box or ellipsoid punch size in h, k and l space is specified. Finally, an outlier parameter is defined giving the range of data to be excluded based on the interquartile range of the data within the punch. A default value of 1.5 is the standard definition of an outlier: the parameter (1.5) times the interquartile range excludes the data above and below the third and first quartile, respectively. This is a simple approach to removing Bragg peaks. The punch parameters are adjustable and multiple passes may be performed to give a clean data set for refinement. The punch may be reset at any time. Future implementations of the program may include more advanced algorithms, such as KAREN (Weng et al., 2020) which identifies and removes Bragg peaks from the total scattering without the need for the unit-cell data. In addition, it can remove other anomalies due to sample environment.

The steps for generating clean volume data for refinement are shown in Fig. 13 for bixbyite: the finalized data set is produced by rebinning, cropping and Bragg-peak removal. The original data set of size 501 × 501 × 501 with h, k and l range −10 to 10 is reduced to a smaller bin size of 0.8 in each dimension, and cropped to a range of 0–4 in h and k and 0–2 in l. The Bragg peaks are then removed in one pass using the default parameters and a box size of five voxels. The resulting data set may be further investigated at different slices and the displayed figure may be saved. At any time during this process, the original data set can be recovered using the reset button.

Figure 13
(a) The GUI `Intensity' tab for generating raw volume data showing (b) the original NeXus file for bixbyite. The figure also shows the data set (c) rebinned with a larger bin size, (d) cropped to a smaller region and (e) with Bragg peaks removed.

4.3. `Refinement' tab

Having built the supercell and preprocessed the intensity data set, a refinement is now set up. The interface for refinement is flexible with options for both magnetic and structural refinements. The Gaussian filter size may be specified for each of the h, k and l dimensions to reduce the noise in the refined intensity and/or mimic the effects of the instrumental resolution broadening. The number of cycles used for the refinement can be specified. One cycle is defined as the number of proposed moves which corresponds to the number of atoms in the supercell. It is typically best to start with a short refinement (1–10 cycles) to optimize the annealing temperature prefactor and decay constant. In the first 10–20% of the refinement, the temperature should be high enough that all of the proposed moves are accepted. This allows the system to overcome any local minima as it progresses toward a global minimum. During this part of the refinement, the temperature should cool gradually such that some of the bad moves are rejected. By the final 50% of the long refinement (100–200 cycles), nearly all of the bad moves should be rejected as the system is moving toward the global minimum.

These accepted and rejected moves can be monitored in a plot window that is automatically updated throughout the refinement and gives the goodness-of-fit value (χ²) of accepted and rejected moves. This window offers several other diagnostic plots, including overall χ², energy, temperature and scale factor. In addition, the refined intensity plot can be compared with the experimental data set. Other parameters relevant to each refinement type (magnetic, occupational, displacive) can be specified before running the refinement. A batch job can be specified to give multiple sequential refinements that can be used to improve the refinement statistics by averaging. An example of a completed refinement is shown in Fig. 14, which gives the final screenshot in (a) and an illustration of the refinement progress over time in (b)–(e). The tab allows for user interaction during a refinement and the parameters may be updated as it progresses. At any time, the refinement can be stopped or reset with new parameters as needed. Further capabilities to control the refinement are expected additions, including combined occupational/displacive disorder of polyhedra.

Figure 14
(a) The GUI `Refinement' tab showing progress of the calculated intensity during the refinement (b)–(e) from 0–100%. The moves that are accepted and rejected are displayed in auxiliary monitor plots that update automatically along with the intensity throughout the refinement.

4.4. `Correlations' tab

The `Correlations' tab provides the ability to calculate both spherically averaged and three-dimensional correlations for occupational, displacement and magnetic disorder. A cutoff distance is employed to limit the search of pairs using the underlying kD tree. The distance tolerance is also employed to account for rounding error in distances. Too tight a tolerance may identify more pairs than actually exist by subtle differences in decimal places introduced by rounding error. Too loose a tolerance may not distinguish real differences in bond lengths. Support for size-effect correlations will also be implemented in future releases.

In the case of spherically averaged correlations, a line plot up to a specified radial cutoff distance provides a simple way to identify the strongest correlations. Three-dimensional correlations are visualized as two-dimensional slices through the specification of lattice planes. Future capabilities will allow one-dimensional plots along particular directions. For any pair correlations, particular atom-pair types may be separated out to help visualize the strength of their interaction. This is especially important if two or more different pairs share the same separation vectors. Without accounting for pair identity, they cannot be distinguished if they have the same separation. However, the correlations may be recalculated with these pairs that are the same distance apart averaged together. In bixbyite, several Fe–Fe, Fe–Mn and Mn–Mn pairs share the same separation vector. Shown in Fig. 15(a) are the calculated three-dimensional correlations at z/a = 0 with all pairs averaged together. The calculated spherical average and three-dimensional correlations may be exported as a comma-separated values (CSV) file or as visualization toolkit (VTK) files, respectively, for further processing. For example, Figs. 15(b) and 15(c) show the three-dimensional correlations visualized with the external viewer ParaView (Ayachit, 2015 ) along two different orientations. Further correlation calculations are expected additions, including the occupancy–displacement pair correlations. The capability to map the calculated correlations to the average structure and visualize them is also planned.

Figure 15
(a) The GUI `Correlations' tab showing calculated three-dimensional magnetic spin-pair correlations. The figure also displays visualizations of nearest neighbors in ParaView viewed (b) along the [001] axis, showing those observed in the z/a = 0 plane corresponding to the 3D-ΔPDF, and (c) along the [111] axis, showing those corresponding to the nearly hexagonal frustrated network.

4.5. `Recalculation' tab

Having obtained the refined structure, it is possible to recalculate the diffuse-scattering pattern back over the original (or any) size with resolution higher than that used for refinement. It is then possible to examine and compare with the original experimental data set. During this recalculation, the symmetry in reciprocal space may optionally be restricted. The Laue symmetry can be selected or inferred from the loaded CIF file in the `Crystal' tab. Of course, the diffuse-scattering pattern is not required to have the same Laue symmetry of the average structure but instead may have symmetry that is lower. For the case of bixbyite as shown in Fig. 16, the symmetry of the diffuse scattering appears to have the same Laue symmetry as the average structure $[m\bar{3}]$ . The recalculated scattering pattern can be exported to a VTK file for external visualization in ParaView.

Figure 16
The GUI `Recalculation' tab showing the recalculated scattering intensity. The recalculation is symmetrized using the underlying $[m\bar{3}]$ Laue symmetry of bixbyite. The resulting pattern compares favorably with the full-volume experimental data set despite using a smaller region for refinement.

5. Program availability

The program is built to accommodate new features based on user needs and evolving data approaches. The source code is located on GitHub for download. An accompanying web site containing tutorials and installation instructions is also available: https://zjmorgan.github.io/rmc-discord/. The content includes example refinements of classical disordered systems with data and instructions to help new users get started. This includes disorder on triangular, honeycomb, kagome and pyrochlore lattices.

6. Limitations and guidance

This program is only suited for crystalline materials and not amorphous or powder materials. Programs like RMCProfile (Tucker et al., 2007) can handle those cases. The maximum correlation length that can be resolved in the refinement along any supercell dimension is half the size of the corresponding supercell length. This also determines the sharpness of the features in reciprocal space. A three-dimensional volume data set is required, which appropriately weights measurement scans together and normalizes for Lorentz and spectrum corrections. To refine the diffuse scattering only, Bragg peaks can be removed, but information about the absolute scaling of magnetic moments, site occupancies and atomic displacements is lost in the overall scale factor. Diffuse magnetic scattering data should be isolated by subtracting from the low-temperature data a higher-temperature data set (Paddison et al., 2013; Roth et al., 2018). Alternatively, background subtractions should be performed from empty instrument measurements scaled to the same incident flux as the normalization data (Michels-Clark et al., 2016). Although powder lines that do not overlap with diffuse-scattering features do not affect the overall refinement as demonstrated with the Ba₃Co₂O₆(CO₃)_0.7 example, it is generally advised to exclude them in areas where there is severe contamination with the diffuse-scattering features.

7. Future development

An opportunity for the single-crystal RMC approach is a co-refinement of the diffuse and Bragg scattering to obtain both the short- and long-range ordering (Dove et al., 2002 ; Proffen & Kim, 2009 ). By removing Bragg peaks, average structure information is lost that may contain critical information about the origin of the disorder. Future implementations will include refining the integrated Bragg intensities which would give average structure constraints, guiding the RMC refinement and ensuring consistency over the full range of the collected data. It would also allow the absolute values for moments, occupancies and displacements to be obtained. For complete volume data sets, 3D-PDF analysis could also be incorporated prior to performing RMC refinements of the total scattering, like that implemented by the Yell program. Previous studies demonstrated this approach using PDF analysis before RMC refinement with powder data (Whitfield et al., 2016 ). Also, the data-preprocessing step would be simplified by eliminating the need for Bragg-peak removal. This approach may benefit magnetic structure analysis for magnetically dilute systems and may also allow for analysis of structures with twinning.

Future work is planned to improve upon the molecular-structure-factor work, extending the displacement and rotation of rigid structures to account for networks of corner-linked polyhedra with constraints. This will enable the diffuse-scattering refinement of networks with rigid units. Similarly, the use of Z matrices for molecules (Goossens et al., 2011 ) is another capability that could be implemented, allowing for the flexibility and rigidity of various portions of a molecule to be specified. An interface to interactively build and visualize polyhedra within the GUI is also planned.

Recently, a formalism for interpreting diffuse scattering based on a disordered superspace approach has been introduced for systems with substitutional disorder (Schmidt & Neder, 2019 ). Physical-space structures are realized that generate diffuse-scattering patterns with specified width and location in reciprocal space. Compared with RMC where the interpretation of diffuse scattering is in terms of correlation parameters, this method has the potential to yield a structure solution with more physical meaning at a reduced computational cost, allowing direct interpretation of the diffuse-scattering features in terms of the modulation functions. An innovation for RMC that considers the symmetry of disorder is highly desired. Symmetry analysis is key to understanding distortion modes that are responsible for unusual properties, for example in layered perovskites that undergo switching from positive to negative thermal expansion (Senn et al., 2016 ). One opportunity is to utilize ISODISTORT to identify distortion symmetries from the input CIF file of the parent structure (Campbell et al., 2006 ). This could be used to restrict the possible distortions of individual atoms in disordered crystals that undergo phase transformations, to speed up RMC refinements and improve the analysis of these systems.

8. Conclusions

A complete and user-friendly RMC refinement program has been developed to analyze single-crystal diffuse-scattering data. The program addresses each of the three disorder types: occupational, displacive and magnetic. It is optimized for performance and is accompanied by a user-friendly graphical interface that allows a user to load in collected data, process them efficiently, perform a refinement and extract correlations that describe the disorder within the crystal. It effectively demonstrates the extraction of correlations from single-crystal diffuse-scattering data in two different scientific user cases. In mineral bixbyite, RMC offers richer detail of the spin-pair antiferromagnetic correlations than the magnetic 3D-ΔPDF by distinguishing the subtle differences in separation vectors between first- and second-nearest neighbors of Fe–Mn and Mn–Mn pairs. In the triangular lattice system Ba₃Co₂O₆(CO₃)_0.7, RMC reveals no clustering of vacant CO₃ sites along c-axis chains, and the occupied CO₃ sites move toward their closest neighboring vacant sites along the chain to relax strain induced by electrostatic interactions.

Acknowledgements

Acknowledgments are given to Thomas Proffen (ORNL) for modeling suggestions of CO₃ in Ba₃Co₂O₆(CO₃)_0.7, Xiaoping Wang (ORNL) for help with the average structure refinements of Ba₃Co₂O₆(CO₃)_0.7 and Christina Hoffmann (ORNL) for helpful discussions regarding the development of the software.

Funding information

Research at Oak Ridge National Laboratory (ORNL) Spallation Neutron Source (SNS) was sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, US Department of Energy (DOE). ZJM acknowledges support from the Higher Education Research Experience at ORNL. The work of HDZ at the University of Tennessee [crystal growth of Ba₃Co₂O₆(CO₃)_0.7] was supported by the US DOE, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under award DE-SC-0020254.

References

Aroyo, M. I., Kirov, A., Capillas, C., Perez-Mato, J. M. & Wondratschek, H. (2006). Acta Cryst. A62, 115–128. CrossRef CAS IUCr Journals Google Scholar
Aroyo, M. I., Perez-Mato, J. M., Capillas, C., Kroumova, E., Ivantchev, S., Madariaga, G., Kirov, A. & Wondratschek, H. (2006). Z. Kristallogr. 221, 15–27. CrossRef CAS Google Scholar
Aroyo, M. I., Perez-Mato, J. M., Orobengoa, D., Tasci, E., de la Flor, G. & Kirov, A. (2011). Bulg. Chem. Commun. 43, 183–197. CAS Google Scholar
Ashcroft, N. W. & Mermin, N. D. (1976). Solid State Physics. New York: Holt, Rinehart and Winston. Google Scholar
Ayachit, U. (2015). The Paraview Guide (Full Color Version): a Parallel Visualization Application. Kitware. Google Scholar
Belsky, A., Hellenbrandt, M., Karen, V. L. & Luksch, P. (2002). Acta Cryst. B58, 364–369. CrossRef CAS IUCr Journals Google Scholar
Boulahya, K., Amador, U., Parras, M. & González-Calbet, J. M. (2000). Chem. Mater. 12, 966–972. CrossRef ICSD CAS Google Scholar
Bramwell, S. T. & Gingras, M. J. P. (2001). Science, 294, 1495–1501. Web of Science CrossRef PubMed CAS Google Scholar
Briggs, W. L. & Henson, V. E. (1995). The DFT: an Owner's Manual for the Discrete Fourier Transform. Philadelphia: Society for Industrial and Applied Mathematics. Google Scholar
Butler, B. D. & Welberry, T. R. (1993). Acta Cryst. A49, 736–743. CrossRef CAS Web of Science IUCr Journals Google Scholar
Butler, B. D., Withers, R. L. & Welberry, T. R. (1992). Acta Cryst. A48, 737–746. CrossRef CAS IUCr Journals Google Scholar
Campbell, B. J., Stokes, H. T., Tanner, D. E. & Hatch, D. M. (2006). J. Appl. Cryst. 39, 607–614. Web of Science CrossRef CAS IUCr Journals Google Scholar
Dove, M. T., Tucker, M. G. & Keen, D. A. (2002). Eur. J. Mineral. 14, 331–348. Web of Science CrossRef CAS Google Scholar
Frey, F. (1995). Acta Cryst. B51, 592–603. CrossRef CAS Web of Science IUCr Journals Google Scholar
Goossens, D. J., Heerdegen, A. P., Chan, E. J. & Welberry, T. R. (2011). Metall. Mater. Trans. A, 42, 23–31. Web of Science CrossRef CAS Google Scholar
Gražulis, S., Chateigner, D., Downs, R. T., Yokochi, A. F. T., Quirós, M., Lutterotti, L., Manakova, E., Butkus, J., Moeck, P. & Le Bail, A. (2009). J. Appl. Cryst. 42, 726–729. CrossRef IUCr Journals Google Scholar
Groom, C. R., Bruno, I. J., Lightfoot, M. P. & Ward, S. C. (2016). Acta Cryst. B72, 171–179. Web of Science CrossRef IUCr Journals Google Scholar
Hester, J. R. (2006). J. Appl. Cryst. 39, 621–625. Web of Science CrossRef CAS IUCr Journals Google Scholar
Igarashi, K., Shimizu, Y., Satomi, E., Kobayashi, Y., Takami, T. & Itoh, M. (2012). J. Phys. Conf. Ser. 400, 032024. CrossRef Google Scholar
Iwasaki, K., Yamamoto, T., Yamane, H., Takeda, T., Arai, S., Miyazaki, H., Tatsumi, K., Yoshino, M., Ito, T., Arita, Y., Muto, S., Nagasaki, T. & Matsui, T. (2009). J. Appl. Phys. 106, 034905. CrossRef Google Scholar
Keen, D. A. & Goodwin, A. L. (2015). Nature, 521, 303–309. Web of Science CrossRef CAS PubMed Google Scholar
Könnecke, M., Akeroyd, F. A., Bernstein, H. J., Brewster, A. S., Campbell, S. I., Clausen, B., Cottrell, S., Hoffmann, J. U., Jemian, P. R., Männicke, D., Osborn, R., Peterson, P. F., Richter, T., Suzuki, J., Watts, B., Wintersberger, E. & Wuttke, J. (2015). J. Appl. Cryst. 48, 301–305. Web of Science CrossRef IUCr Journals Google Scholar
Lisher, E. J. & Forsyth, J. B. (1971). Acta Cryst. A27, 545–549. CrossRef IUCr Journals Web of Science Google Scholar
McGreevy, R. L. (2001). J. Phys. Condens. Matter, 13, R877–R913. Web of Science CrossRef CAS Google Scholar
McGreevy, R. L. & Pusztai, L. (1988). Mol. Simul. 1, 359–367. Web of Science CrossRef Google Scholar
Michels-Clark, T. M., Savici, A. T., Lynch, V. E., Wang, X. P. & Hoffmann, C. M. (2016). J. Appl. Cryst. 49, 497–506. Web of Science CrossRef CAS IUCr Journals Google Scholar
Momma, K. & Izumi, F. (2008). J. Appl. Cryst. 41, 653–658. Web of Science CrossRef CAS IUCr Journals Google Scholar
Nield, V. M., Keen, D. A. & McGreevy, R. L. (1995). Acta Cryst. A51, 763–771. CrossRef CAS Web of Science IUCr Journals Google Scholar
Paddison, J. A. M. (2019). Acta Cryst. A75, 14–24. Web of Science CrossRef IUCr Journals Google Scholar
Paddison, J. A. M., Ross Stewart, J. & Goodwin, A. L. (2013). J. Phys. Condens. Matter, 25, 454220. Web of Science CrossRef PubMed Google Scholar
Pauling, L. & Shappell, M. (1930). Z. Kristallogr. Cryst. Mater. 75, 128–142. CrossRef CAS Google Scholar
Proffen, T. (2000). Z. Kristallogr. Cryst. Mater. 215, 661–668 CrossRef CAS Google Scholar
Proffen, T. & Kim, H. (2009). J. Mater. Chem. 19, 5078–5088. Web of Science CrossRef CAS Google Scholar
Proffen, Th. & Neder, R. B. (1997). J. Appl. Cryst. 30, 171–175. CrossRef CAS Web of Science IUCr Journals Google Scholar
Proffen, Th. & Welberry, T. R. (1997). Acta Cryst. A53, 202–216. CrossRef CAS Web of Science IUCr Journals Google Scholar
Proffen, Th. & Welberry, T. R. (1998a). J. Appl. Cryst. 31, 318–326. CrossRef CAS IUCr Journals Google Scholar
Proffen, T. & Welberry, T. R. (1998b). Phase Transit. 67, 373–397. CrossRef CAS Google Scholar
Ramirez, A. P. (1994). Annu. Rev. Mater. Sci. 24, 453–480. CrossRef CAS Web of Science Google Scholar
Roth, N., May, A. F., Ye, F., Chakoumakos, B. C. & Iversen, B. B. (2018). IUCrJ, 5, 410–416. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Roth, N., Ye, F., May, A. F., Chakoumakos, B. C. & Iversen, B. B. (2019). Phys. Rev. B, 100, 144404. CrossRef Google Scholar
Schmidt, E. M. & Neder, R. B. (2019). Phys. Rev. B, 100, 054201. Web of Science CrossRef Google Scholar
Schultz, A. J., Jørgensen, M. R. V., Wang, X., Mikkelson, R. L., Mikkelson, D. J., Lynch, V. E., Peterson, P. F., Green, M. L. & Hoffmann, C. M. (2014). J. Appl. Cryst. 47, 915–921. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sears, V. F. (1992). Neutron News, 3(3), 26–37. CrossRef Google Scholar
Seltzer, S. (1995). X-ray Form Factor, Attenuation and Scattering Tables. NIST Standard Reference Database 66. National Institute of Standards and Technology, Gaithersburg, MD, USA. Google Scholar
Senn, M. S., Murray, C. A., Luo, X., Wang, L., Huang, F.-T., Cheong, S.-W., Bombardi, A., Ablitt, C., Mostofi, A. A. & Bristowe, N. C. (2016). J. Am. Chem. Soc. 138, 5479–5482. CrossRef ICSD CAS PubMed Google Scholar
Simonov, A., Weber, T. & Steurer, W. (2014). J. Appl. Cryst. 47, 1146–1152. Web of Science CrossRef CAS IUCr Journals Google Scholar
Tucker, M. G., Keen, D. A., Dove, M. T., Goodwin, A. L. & Hui, Q. (2007). J. Phys. Condens. Matter, 19, 335218. Web of Science CrossRef PubMed Google Scholar
Weber, T. & Simonov, A. (2012). Z. Kristallogr. 227, 238–247. Web of Science CrossRef CAS Google Scholar
Welberry, T. R. (1986). J. Appl. Cryst. 19, 382–389. CrossRef CAS Web of Science IUCr Journals Google Scholar
Welberry, T. R. (2004). Diffuse X-ray Scattering and Models of Disorder. IUCr Monographs on Crystallography, No. 16. Oxford University Press. Google Scholar
Welberry, T. R. & Butler, B. D. (1994). J. Appl. Cryst. 27, 205–231. CrossRef CAS Web of Science IUCr Journals Google Scholar
Welberry, T. R. & Butler, B. D. (1995). Chem. Rev. 95, 2369–2403. CrossRef CAS Web of Science Google Scholar
Welberry, T. R. & Goossens, D. J. (2008). Acta Cryst. A64, 23–32. Web of Science CrossRef CAS IUCr Journals Google Scholar
Welberry, T. R. & Goossens, D. J. (2014). IUCrJ, 1, 550–562. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Welberry, T. R. & Mayo, S. C. (1996). J. Appl. Cryst. 29, 353–364. CrossRef CAS Web of Science IUCr Journals Google Scholar
Welberry, T. R. & Proffen, Th. (1998). J. Appl. Cryst. 31, 309–317. Web of Science CrossRef CAS IUCr Journals Google Scholar
Weng, J., Dill, E. D., Martin, J. D., Whitfield, R., Hoffmann, C. & Ye, F. (2020). J. Appl. Cryst. 53, 159–169. Web of Science CrossRef CAS IUCr Journals Google Scholar
Whitfield, R. E., Goossens, D. J. & Welberry, T. R. (2016). IUCrJ, 3, 20–31. Web of Science CrossRef CAS PubMed IUCr Journals Google Scholar
Ye, F., Liu, Y., Whitfield, R., Osborn, R. & Rosenkranz, S. (2018). J. Appl. Cryst. 51, 315–322. Web of Science CrossRef CAS IUCr Journals Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 54| Part 6| December 2021| Pages 1867-1885

https://doi.org/10.1107/S1600576721010141

Open

access

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

computer programs\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

rmc-discord: reverse Monte Carlo refinement of diffuse scattering and correlated disorder from single crystals

1. Introduction

2. Reverse Monte Carlo approach

2.1. Single-crystal diffuse-scattering intensity calculations

2.2. Storing of intermediate results in memory

2.3. Implementation and libraries

2.4. Capabilities and user requirements

3. Scientific user cases

3.1. Magnetic diffuse scattering of frustrated magnet bixbyite

3.2. Structural diffuse scattering of disordered Ba3Co2O6(CO3)0.7

4. Graphical user interface

4.1. `Crystal' tab

4.2. `Intensity' tab

4.3. `Refinement' tab

4.4. `Correlations' tab

4.5. `Recalculation' tab

5. Program availability

6. Limitations and guidance

7. Future development

8. Conclusions

Acknowledgements

Funding information

References

computer programs

3.2. Structural diffuse scattering of disordered Ba₃Co₂O₆(CO₃)_0.7