Static Structure Factor S(q) or S(k)

Reference: K. Zhang, On the Concept of Static Structure Factor, arXiv:1606.03610 (2016) (link)

see code at

https://github.com/statisticalmechanics/gr

Also check a more comprehensive review on this topic

Unifying the concepts of scattering and structure factor in ordered and disordered samples (link), Dingning Li and Kai Zhang, J. Appl. Cryst. 54, 644-660 (2021)

and source code for it

https://github.com/statisticalmechanics/scatter

LJ system

Lammps input script is here

configuration snapshots are stored in "movie.xyz" in xyz format

C code to calculate g(r) gr.c , which generates "gr.dat"

C code to calculate S(q) sq.c , which calculates S(q) with two methods as bellow.

*****************************************************************************************************************************************************

After being complied (with executive file e.g. sq), the "sq.c" code can be run in command line as

./sq

Some parameters need to be input in the window.

NumberOfParticles // total number N of particles in the system

L // simulation box size, assuming cubic shape. This is used to define wavevector increment dk=2 pi/L, and may allow to be scaled by fluctuating box size (not implemented in the current version)

NumberOfrBins // number of r-bins in the file "gr.dat" containing the radial distribution function g(r)

NumberOfqBins // number of q-bins in S(q) used in the discrete numerical integration of the Fourier transform

qmax // the maximum q value for S(q) in the Fourier transform method

NumberOfConfig // total number of configurations in your xyz movie file. put 1 if you have only one snapshot

NumberOfEqui // number of configurations considered as equilibration and to be skipped in your xyz movie file. put 0 if the sampling starts from the first snapshot

nkbound // number of k vectors in each direction (constrained by the allowed size of the data structure). So in total nkbound^3 wavevectors

dkmultiplier // put 1 if you want the increment dk to be 2 pi / L. Otherwise, it is dk = 2 pi / L * dkmultiplier

*****************************************************************************************************************************************************

Method 1. Fourier transform of g(r)

comment: Not defined at k=0. Also should not apply to crystals. Oscillation of S(k) for k<k* is due to oscillation of g(r) about 1 at r>r*.

Method 2. Direct calculation with Fourier density correlation.

comment: S(k) -> N as k->0, due to finite size and periodic boundary condition, where N is the total number of particles. S(k) scales with N at peak k's for crystals. The noise of the data can be reduced by choosing smaller increment dk (e.g. using larger box size L in the input)

1. rho = 0.8 at T = 2.0, liquid