Tables for
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 4.2, pp. 528-530   | 1 | 2 |

Section 4.2.7. Computer simulations and modelling

F. Frey,a H. Boysena and H. Jagodzinskia

aDepartment für Geo- und Umweltwissenschaften, Sektion Kristallographie, Ludwig-Maximilians Universität, Theresienstrasse 41, 80333 München, Germany

4.2.7. Computer simulations and modelling

| top | pdf | Introduction

| top | pdf |

The various analytical expressions given in Section 4.2.5[link] are mostly rather complex, so their application is often restricted to relatively simple systems. Even in these cases analytical solutions are often not available. For larger displacements the approximations that use an expansion of the exponentials [e.g. equation ([link]] are no longer valid. Hence there is a need for alternative approaches to tackling more complicated disorder models. In the past, optical transforms (e.g. Harburn et al., 1974[link], 1975[link]) and the videographic method (Rahman, 1993[link]) were developed for this purpose. In the first method, a 2D mask with holes with sizes that represent the scattering power of the atoms is generated, which is then subjected to coherent light (from a laser) to produce a diffraction image. Problems arise for strong scatterers requiring very large holes. These problems were overcome in the second method, where the mask is replaced by a computer image with intensities proportional to the scattering power for each pixel. With the advent of more and more powerful computers these methods are now replaced by complete computer simulations, both to set up the disorder model and to calculate the diffraction pattern. Simulation programs

| top | pdf |

Having established a disordered crystal with the types and positions of all atoms involved (a configuration), e.g. by using one of the methods described below, computer programs employing fast-Fourier-transform techniques can be used to calculate the diffraction pattern, which may be compared with the observation. It has to be borne in mind, however, that there is still a large gap between a real crystal with its ~1023 atoms and the one simulated by the computer with only several thousand atoms. This means that very long range correlations can not be included and have to be treated in an average manner. Furthermore, the limited size of the simulated crystal leads to termination effects, giving rise to considerable noise in the calculated diffraction pattern. Butler & Welberry (1992[link]) have introduced a technique to avoid this problem in their program DIFFUSE by dividing the simulated crystal into smaller `lots'. For each lot the intensity is calculated and then the intensities of all lots are summed up incoherently, which finally results in a smooth intensity distribution. Note, however, that in this case long-range correlations are restricted to even smaller values. To overcome this problem Boysen (1997[link]) has proposed a method for suppressing the subsidiary maxima by multiplying the scattering density of the model crystal by a suitably designed weighting function simulating the effect of the instrumental resolution function.

A very versatile computer program, DISCUS, which allows not only the calculation of the scattering intensities but also allows the model structures to be built up in various ways, has been designed by Proffen & Neder (1999[link]). It contains modules for reverse Monte Carlo (RMC, see below) simulation, the calculation of powder patterns, RMC-type refinement of pair-distribution functions (PDFs, see below) and many other useful tools for analysing disorder diffuse scattering.

Other computer programs have been developed to calculate diffuse scattering, some for specific tasks, such as SERENA (Micu & Smith, 1995[link]), which uses a collection of atomic configurations calculated from a molecular dynamics simulation of molecular crystals. Modelling procedures

| top | pdf |

Several well established methods can be used to create the simulated crystal on the computer. With such a crystal at hand, it is possible to calculate various thermodynamical properties and study the effect of specific parameters of the underlying model. By Fourier transformation, one obtains the diffraction pattern of the total scattering, i.e. the diffuse intensities and the Bragg peaks. The latter may lead to difficulties due to the termination effects mentioned in Section[link]. These may be circumvented by excluding regions around the Bragg peaks (note, however, that in this case valuable information about the disorder may be lost), by subtracting the average structure or by using the approximation of Boysen (1997[link]) mentioned above. A major advantage of such modelling procedures is that realistic physical models are introduced at the beginning, providing further insight into the pair interactions of the system, which can only be obtained a posteriori from the correlation parameters or fluctuation wave amplitudes derived from one of the methods described in Section 4.2.5[link]. Molecular dynamics

| top | pdf |

Molecular dynamics (MD) techniques have been developed to study the dynamics of a system. They may also be used to study static disorder problems (by taking time averages or snapshots), but they are particularly useful in the case of dynamic disorder, e.g. diffusing atoms in superionic conductors. The principle is to set up a certain configuration of atoms with assumed interatomic potentials Φij(rij) and subject them to Newton's equation of motion,[F_i(t)=m_i\, {{\rm d}^2r_i(t)\over{\rm d}t^2},\eqno(]where the force Fi(t) is calculated from the gradient of Φ. The equations are solved approximately by replacing the differential dt by a small but finite time step Δt to find new positions ri(t + Δt). This is repeated until an equilibrium configuration is found. MD techniques are quite useful if only short-range interactions are effective, even allowing the transfer of potential parameters between different systems, but are less reliable in the presence of significant long-range interactions. Monte Carlo calculations

| top | pdf |

Monte Carlo (MC) methods appear to be more suited to the study of static disorder and many examples of their application can be found in the literature. Different variants allowing simulations and refinements have been applied:

  • (1) Direct MC simulation. In this method, introduced by Metropolis et al. (1953[link]), a starting configuration is again set up in accordance with the known average structure and other crystal, chemical and physical information if available, and an appropriate interaction potential is chosen. The choice of this potential can be a quite delicate task. On the one hand, it should be as simple as possible to allow as large a simulation box size as can reasonably be handled by the computer's capacity. On the other hand, it must be detailed enough to include all the relevant interaction parameters of the system. Frequently used interaction potentials are a pseudo spin Ising Hamiltonian, where the `spins' can either be binary (e.g. atom types or molecular orientations) or continuous (displacements), and harmonic springs for the displacements.

    Having set up the starting configuration and defined the Hamiltonian, one proceeds by choosing a specific site at random and changing its parameters (occupancies and positions) by a random amount (a `move'). The energy of this new configuration is calculated and compared with the old one. If the difference ΔE is negative, the move is accepted. If it is positive it is accepted with a probability[P=\exp\{-\Delta E/k_{\rm B}T\}/[1+\exp\{-\Delta E/k_{\rm B}T\}], \eqno(]where T is a temperature and kB is Boltzmann's constant. Then another site is chosen at random and the process is repeated again and again until an equilibrium configuration is found, i.e. until the energies fluctuate around some average value. After this, the corresponding intensity distribution is calculated and compared with the observed one. The parameters of the Hamiltonian may then be modified to improve the agreement between the calculated and observed diffraction patterns. This rather cumbersome trial-and-error method may be used to study the influence of the various parameters of the model. Further details of this method together with some illustrative examples may be found in Welberry & Butler (1994[link]).

    In an attempt to automate the adjustment of the model parameters, Welberry et al. (1998[link]) have used a least-squares procedure to minimize[\chi^2=\textstyle \sum \limits w(\Delta I)^2, \eqno(]where ΔI are the differences between the observed and calculated intensities and w are the appropriate weights. The problem with this approach is that the necessary differentials δΔIpi must be calculated numerically by performing full MC simulations with parameters pi and pi + δpi at each iteration step, which presents a formidable task even for the fastest modern computers.

  • (2) Reverse MC calculations. Application of the direct MC method may be very time consuming, as many MC simulations are necessary to arrive at a final configuration. To overcome this problem, the so-called reverse Monte Carlo (RMC) method has been developed, originally for liquids and glasses (McGreevy & Pusztai, 1988[link]) and later for disordered single crystals (Nield et al., 1995[link]).

    RMC is a model-free approach, i.e. without the need to define a proper interaction Hamiltonian. Otherwise it proceeds in a very similar way to direct MC analysis. A starting configuration is set up and random `moves' of the atoms are carried out. The only difference is that acceptance or non-acceptance of a move is based on the agreement between observed and calculated diffraction intensities, i.e. equation ([link] is replaced by[P=\exp \{-\Delta \chi^2/2\}, \eqno(]where χ2 is defined in equation ([link] and [\Delta\chi^2 = \chi^2_{\rm new} - \chi^2_{\rm old}]. Usually only a single MC run is necessary to arrive at a configuration with a diffraction pattern consistent with the observation. A drawback of the RMC technique is that usually only one configuration is found satisfying the observed diffraction pattern. In fact, it is possible that different configurations may produce the same or very similar intensity distributions (see e.g. Welberry & Butler, 1994[link]). In this context, it has to be borne in mind that the diffuse scattering contains information about two-body correlations only, while the physical reason for a particular disordered structure may also be influenced by three- and many-body correlations. Such many-body correlations can easily be incorporated in the direct MC method. Moreover, RMC is susceptible to fitting artefacts (noise) in the diffraction pattern. To avoid unreasonable atomic configurations, restrictions such as a limiting nearest-neighbour approach can be built in. Weak constraints, e.g. ranges of interatomic distances or bond angles (i.e. three-body correlations) can be introduced by adding further terms in the definition of χ2 [equation ([link]]. In the same way, multiple data sets (e.g. neutron and X-ray data, EXAFS measurements etc.) may be incorporated. For more details and critical reviews of this method see e.g. McGreevy (2001[link]) and Welberry & Proffen (1998[link]). Many applications of this method can be found in the literature, mainly for powder diffraction, but also for single-crystal data.

  • (3) Simulated annealing and evolutionary algorithms. A general problem with MC methods is that they may easily converge to some local minima without having found the global one. One way to reduce such risks is to start with a high probability in ([link] by using a large `temperature' T, i.e. initially allowing many `false' moves before gradually reducing T during the course of the simulation cycles. This is called `simulated annealing', although it has nothing to do with real annealing. In the RMC method one may introduce a weighting parameter similar to T in ([link].

    Another effective way to find the global minimum and also to accelerate the optimization of the energy parameters of a given MC model is by using the principles of evolutionary theory: selection, recombination and mutation. Two such evolutionary (or genetic) algorithms have been proposed: the differential evolution (DE) algorithm (Weber & Bürgi, 2002[link]) and the cooperative evolution (CE) algorithm (Weber, 2005[link]). A single parameter of the model is called a gene and a set of genes is called a chromosome, p being the genotype of an individual. First a population consisting of several individuals is built up and the corresponding diffraction patterns are calculated and compared with the observation. The fitness of each individual is quantified by χ2. Children are then created by choosing one parent individual and calculating the second one from three randomly chosen individuals according to[p_c' = p_c +f_{\rm m}(p_a -p_b), \eqno(]where fm governs the mutations. The chromosome of the child is obtained by combining the genes of the two parents governed by a crossover (or recombination) constant fr. If its fitness is higher than that of the parent, it replaces it. This procedure is repeated until some convergence criterion is reached. This DE method is still rather time-consuming on the computer. Therefore, the CE technique was introduced, where only one crystal is built up during the refinement. Here a large population is created spanning a large but reasonable parameter space. Individuals are selected at random to decide upon acceptance or rejection of an MC move via their own energy criteria. Then χ2 is calculated and according to its positive or negative change a gratification or penalty weight is given to that individual. It may live as long as this weight is positive, otherwise it is replaced by a new individual calculated from ([link]. This way, useful individuals live longer to act on the same crystal, while unsuccessful ones are eliminated early. Note that the recombination operation is not used in this technique.

  • (4) The PDF method. The pair-distribution function (PDF) has long been used for the analysis of liquids, glasses and amorphous substances (Warren, 1969[link]), but has recently regained considerable interest for the analysis of crystalline substances as well (Egami, 1990[link]; Billinge & Egami, 1993[link]; Egami, 2004[link]). The PDF is obtained by a Fourier transformation of the total (Bragg plus diffuse) scattering in a powder pattern,[\rho_0G(r)= \rho_0 +\textstyle{1\over 2}\textstyle\int H[S(H) -1]\sin(2\pi H)\ {\rm d}H. \eqno(]This is nothing other than the van Hove correlation function ([link] or the related Patterson function ([link] averaged spherically and taken at t = 0 (a snapshot) and is given by[G^{\rm PDF}(r)=4\pi r[\rho(r) -\rho_0] =\left(\textstyle \sum \limits c_i b_i\right)^{-2}4\pi r \rho_0G(r), \eqno(]where ρ0 is the average number density and[\rho(r)= (1/N) \textstyle \sum \limits(b_ib_{i'}/\langle b\rangle^2)\delta(r-r_{ii'}).\eqno(]The δ functions are then convoluted with a normalized Gaussian to account for (harmonic) thermal motion. Parameter refinement proceeds in a similar way to the RMC technique. First a model is built, initially within just one unit cell and with periodic boundary conditions, then its PDF is calculated, compared with the observed one and further improved using, for example, MC simulated annealing. The model is then enlarged to include longer-range correlations. Owing to the small size of the models, this technique is much faster than the RMC method. It is essential, however, that data are measured up to very high H values to minimize truncation errors in ([link]. General remarks

| top | pdf |

All of the different modelling techniques mentioned in this section have their specific merits and limitations and have contributed much to our understanding of disorder in crystalline materials following the interpretation of the corresponding diffuse scattering. It should be borne in mind, however, that application of these methods is still far from being routine work and it requires a lot of intuition to ensure that the final model is physically and chemically reasonable. In particular, it must always be ensured that the average structure remains consistent with that derived from the Bragg reflections alone. This may be done by keeping the Bragg reflections, i.e. by analysing the total scattering, or by designing special algorithms, e.g. by swapping two atoms at the same time (Proffen & Welberry, 1997[link]). Moreover, possible traps and corrections like local minima, termination errors, instrumental resolution, statistical noise, inelasticity etc. must be carefully considered. All this means that the analytical methods outlined in Section 4.2.5[link] keep their value and should be preferred wherever possible.

The newly emerging technique of using full quantum-mechanical ab initio calculations for structure predictions along with MD simulations may also be applied to disordered systems. The limited currently available computer power, however, restricts this possibility to rather simple systems and small simulation box sizes, but, with the expected further increase of computer capacities, this may open up new perspectives for the future.


Billinge, S. J. L. & Egami, T. (1993). Short-range atomic structure of Nd2−xCexCuO4−y determined by real-space refinement of neutron-powder-diffraction data. Phys. Rev B, 47, 14386–14406.
Boysen, H. (1997). Suppression of subsidiary maxima in computer simulations of diffraction intensities. Z. Kristallogr. 212, 634–636.
Butler, B. D. & Welberry, T. R. (1992). Calculation of diffuse scattering from simulated disordered crystals: a comparison with optical transforms. J. Appl. Cryst. 25, 391–399.
Egami, T. (1990). Atomic correlations in non-periodic matter. Mater. Trans. JIM, 31, 163–176.
Egami, T. (2004). Local crystallography of crystals with disorder. Z. Kristallogr. 219, 122–129.
Harburn, G., Miller, J. S. & Welberry, T. R. (1974). Optical-diffraction screens containing large numbers of apertures. J. Appl. Cryst. 7, 36–38.
Harburn, G., Taylor, C. A. & Welberry, T. R. (1975). An Atlas of Optical Transforms. London: Bell.
McGreevy, R. L. (2001). Reverse Monte Carlo modelling. J. Phys. Condens. Matter, 13, R877–R913.
McGreevy, R. L. & Pusztai, L. (1988). Reverse Monte Carlo simulation: a new technique for the determination of disordered structures. Mol. Simul. 1, 359–367.
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. & Teller, E. (1953). Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092.
Micu, A. M. & Smith, J. C. (1995). SERENA: a program for calculating X-ray diffuse scattering intensities from molecular dynamics trajectories. Comput. Phys. Commun. 91, 331–338.
Nield, V. M., Keen, D. A. & McGreevy, R. L. (1995). The interpretation of single-crystal diffuse scattering using reverse Monte Carlo modelling. Acta Cryst. A51, 763–771.
Proffen, Th. & Neder, R. B. (1999). DISCUS, a program for diffuse scattering and defect-structure simulations – update. J. Appl. Cryst. 32, 838–839.
Proffen, Th. & Welberry, T. R. (1997). Analysis of diffuse scattering via reverse Monte Carlo technique: a systematic investigation. Acta Cryst. A53, 202–216.
Rahman, S. H. (1993). The local domain configuration in partially ordered AuCu3. Acta Cryst. A49, 68–79.
Warren, B. E. (1969). X-ray Diffraction. Reading: Addison-Wesley.
Weber, T. & Bürgi, H.-B. (2002). Determination and refinement of disordered crystal structures using evolutionary algorithms in combination with Monte Carlo methods. Acta Cryst. A58, 526–540.
Weber, Th. (2005). Cooperative evolution – a new algorithm for the investigation of disordered structures via Mont Carlo modelling. Z. Kristallogr. 220, 1099–1107.
Welberry, T. R. & Butler, B. D. (1994). Interpretation of diffuse X-ray scattering via models of disorder. J. Appl. Cryst. 27, 205–231.
Welberry, T. R. & Proffen, Th. (1998). Analysis of diffuse scattering from single crystals via the reverse Monte Carlo technique. I. Comparison with direct Monte Carlo. J. Appl. Cryst. 31, 309–317.
Welberry, T. R., Proffen, Th. & Bown, M. (1998). Analysis of single-crystal diffuse X-ray scattering via automatic refinement of a Monte Carlo model. Acta Cryst. A54, 661–674.

to end of page
to top of page