Calculation of the structure factor in computer simulations

Formal definition of the structure factor

We will follow here roughly the derivation by Egami and Billinge in Ref.¹, although it can be easily found in many textbooks. We begin with the sample scattering amplitude

with \(\mathbf{Q}\) the diffraction vector or momentum transfer. \(\mathbf{R}_i\) is the position of the particle \(i\), and \(\langle b\rangle\) is the average of the scattering amplitude of each particle in the vacuum \(b_i\). From this moment on, we will consider that all of the atoms are of the same species, \(b_i = b\).

From \(\Psi(\mathbf{Q})\) we define the structure factor \(S(\mathbf{Q})\) as²:

What follows immediately from this expression is that the structure function must be always positive for every value of \(\mathbf{Q}\). We can expand the scattering amplitude and use \(|z| = z\cdot z^*\) and, if all the atoms are of the same type,

Powder average: Debye formula

Usually we are interested in the powder average of the structure factor. This is the structure factor averaged for every possible orientation of the diffraction vector - because in a powder we have a lot of structures randomly oriented. We calculate therefore

This integral can be performed easily if we put the \(z\) axis along with the direction of \(\mathbf{Q}\) and perform the integration rotating the distances \(\mathbf{R}_{ij}\)

This is the famous Debye formula and, since its the average of an always positive quantity, it must be always positive.

Computer simulation

One of the most usual problems when we model and study systems in computer simulations is that we don’t have actual infinite systems. We do, however, use the periodic boundary conditions (PBC) usually to emulate the behavior of infinite systems. With the periodic boundary conditions we usually use the minimum image convention: for the distance between particle \(i\) and \(j\), we use whichever is the closest, considering all the possible positions through the boundaries. We calculate for a very simple test case (a simple cubic 3D lattice with 4x4x4=64 atoms) the structure factor with that prescription

A surprising result! We insisted several times that the structure factor should be always positive, yet we get, using the very same definition, a negative structure factor for wavenumbers near 10. Where did these negative values come from? From the construction that used to help us a lot, the minimum image convention. The pair distance now isn’t always \(r_{ij} = r_j - r_i\), but depends on whether we use the original particles or their images. Therefore, this “new” structure factor isn’t the product of two conjugate complex numbers³. What if we avoid the periodic boundary conditions? We have a comparision of the structure factor with and without boundary conditions (i. e., with the 64 atoms in a void):

This shows that the structure factor, when we use its definition without minimum image convention, is (as expected) always positive.

A bit further

We can also use the pair distribution function and calculate the structure factor as the Fourier Transform. But keep in mind that if you calculate the pair distribution function with PBC, when you get the structure factor related to it you might get negative numbers.

How can we simulate an infinite medium?

The question then, remains: how can we simulate an infinite medium when calculating structure factor? The first answer is that it is not that obvious that we would actually need this infinite medium, since the periodic images of the cell would be aligned in a crystal that might interfere with the structure within the cell — the one we actually do want to study. However, a couple of replicas should be enough to smear out some of the finite size effects. One of the possibilities is to replicate explicitly the box, creating the particles in the neighboring cells by duplication of the original ones. This, though, implies a calculation much harder, since the sum is over \(N^2\) particles, and replicating only one cell right and left in each direction would imply a computational time of \((3^3\cdot N)^2 \approx 700\cdot N^2\). In general, the complexity \(\mathcal{O}(N^2)\) makes structure factor calculation very expensive for large systems.

There is an alternative to add the boundary conditions. We begin with the definition of the sample scattering amplitude, but writing explicitly the periodic boundary images we want to consider:

where \(\mathbf{\Delta L}_j\) is the distance between a particle and its \(j\)-th periodic replica. Since the sums are independent, we can write:

Multiplying by the conjugate gives us the structure factor

The advantage of this calculation is that it is linear in the sum of the number of particles \(N\) and the number of replicas \(M\) consider, \(\mathcal{O}(N+M)\), much lower than the previous \(\mathcal{O}(N^2M^2)\). Consequently, if we want to focus in a region of \(\mathbf{Q}\), this new approach will be useful⁴. We are left with only one detail, respecting to the powder average. It is not trivial how to calculate this integral, since we need to give proper weights to each angle. One of the alternatives are to use the Lebedev quadrature⁵, although other methods like Importance Sampling Montecarlo can be useful in this situation.

Code

Here is the code with which we generated the figures above:

import numpy as np
import itertools as it

def ssf(x, size, q, pbc=False):
  """
  From a series of positions x in a cubic box of length size we get
  the structure factor for momentum q
  """

  natoms = np.shape(x)[0]
  sf = 0.0
  for i in range(natoms):
    x1 = x[i]
    for j in range(i+1, natoms):
      x2 = x[j]
      dx = x2 - x1
      if pbc:
        for i in range(3):
          if dx[i] >  size/2: dx[i] -= size
          if dx[i] < -size/2: dx[i] += size
      r = np.linalg.norm(dx)
      sf += 2*np.sin(q*r)/(q*r)
  sf /= natoms
  sf += 1
  return sf

def generate_sc(size, n):
  """
  Generate the positions of a simple cubic crystal in a box of
  length size with n atoms in each direction (order parameter =
  size/n)
  """

  natoms = n**3
  pos = range(n)
  x = np.zeros((natoms, 3))
  i = 0
  for px, py, pz in it.product(pos, pos, pos):
    x[i] = (px, py, pz)
    x[i] = x[i] * (size/n)
    i += 1
  return x

if __name__ == '__main__':
  import pylab as pl
  size = 1.0
  x = generate_sc(size, 4)
  q = np.linspace(0, 20.0, 101)[1:]
  sf_pbc = [ssf(x, size, _, True) for _ in q]
  sf_sing = [ssf(x, size, _, False) for _ in q]
  fig, ax  = pl.subplots()
  ax.plot(q, sf_pbc)
  ax.axhline(y=0, c='k', ls='--')
  ax.set_xlabel('Wavenumber')
  ax.set_ylabel('Structure factor')
  fig.tight_layout()
  fig.savefig('ssf_pbc.png')
  pl.close()

  fig, ax  = pl.subplots()
  ax.plot(q, sf_pbc, label='with PBC')
  ax.plot(q, sf_sing, label='without PBC')
  ax.axhline(y=0, c='k', ls='--')
  ax.set_xlabel('Wavenumber')
  ax.set_ylabel('Structure factor')
  ax.legend()
  fig.tight_layout()
  fig.savefig('ssf_comp.png')
  pl.close()