I am trying to reproduce the 2D SAXS interference pattern of bulk water via ab initio MD methods as shown in the figure below:

enter image description here

The main aim is to recreate figure (a) via averaging MD configurations of bulk water. To do this. I ran QM/MM simulations of water in CP2K to get the electron density of bulk water. The water molecules were selected so as to create a spherical volume in the middle of the simulation box.

The structure factor $S(\vec{G})$ as a function of the scattering vector $\vec{G}$ can be written as the fourier transform of the electron density:

$$ n(\vec{G}) = \frac{1}{V}\sum\limits_j e^{-i\vec{G}\cdot\vec{r}_j}\int n_j\left(\vec{r}'\right)e^{-i\vec{G}\cdot\vec{r'}}d\vec{r}'\tag{1}$$

\begin{equation} S(\vec{G}) =Vn(\vec{G}) = \sum\limits_j f_j\left(G\right)e^{-i\vec{G}\cdot\vec{r}_j}, \tag{2}\end{equation}

so the simplest way to reproduce the above figure seem to be to perform the fourier transform of the electron density cube files obtained from the MD simulations and average them over multiple configurations. However in my case, I am getting a spike in the q-space corresponding to the zero scattering vector in the center.

enter image description here

Does anyone here has any experience with this type of scenario or has tried to replicate the SXRS spectra of water? Any suggestions would be helpful.

  • $\begingroup$ The g(r) graph looks like the gyration radii that it is a discrete quantity (depend on the position of each particle). I m not seen how to calculate it from electronic density (continuous function of x,y,x). $\endgroup$
    – Camps
    Commented Mar 9, 2021 at 15:46
  • 1
    $\begingroup$ @Camps The g(r) graph is the radial distribution function of water which can be experimentally determined by X-ray (hence based on electronic structure of the system). For ex: A.K. Soper's 2013 paper on water RDF (doi.org/10.1155/2013/279463) $\endgroup$
    – mykd
    Commented Mar 11, 2021 at 19:19

1 Answer 1


In short, there is nothing wrong with the spike that you have obtained. The low-wavevector spike is a general feature of small-angle diffraction. For example, see Figure 1 of this paper where the scattered intensity is plotted towards small scattering vectors. The scattered intensity increases exponentially towards zero scattering, which makes for the spike you see from your simulation.

Experimentally, small-angle scattering experiments are performed on relatively thin samples. This is done to avoid multiple scattering, which will increase with the thickness of sample being measured. And as scattering is generally rare, the majority of the incident beam is transmitted through the sample in the first place, making the transmitted beam the zero limit of scattering vector.

If your goal is to re-create the radial distribution function $g(r)$, it is perfectly okay to cut off the scattering vector at some low magnitude and eliminate the zero-wavevector spike for doing the Fourier transform back to real space. This is okay for two reasons. First, because the size of your simulation is finite, and the simulation only contains information down to a smallest vector, not to zero (that would require a macroscopically large simulation). The information contained in the simulation is limited to a minimum wavevector of $2\pi$ divided by the simulation cell size. Second, computationally the zero-wavevector point in the Fourier transform of real-space densities is the constant term. It arises when the strictly positive radial distirbution functions (or electron densities) are Fourier transformed, and the zero-wavevector term is the average value of the real-space functions.

As a very brief pointer, I would recommend dynasor for convenience of calculating structure factors (static or dynamic) from molecular dynamics. Hope this helps the community compute away.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .