This seems like a reasonable feature to add to freud. This function has a time-dependence to it. We have only one other analysis, MSD, in freud which computes a function with explicit time dependence, so the API should mirror that and take the box as a parameter to the constructor, and then give a 3-dimensional array of positions to the compute method.

In addition, it will be important to average over k-vectors properly. A good example to work from is in the StaticStructureFactorDirect class.

With all of those things in mind, here's how I would propose the new class to look, along with all the properties that make sense to expose:

isf = freud.diffraction.IntermediateScattering(box, k_max, k_min, num_sampled_k_points)
isf.compute(points, query_points, N_total, reset)

# access properties
isf.k_max
isf.k_min
isf.k_points
isf.nbins
isf.num_sampled_k_points
isf.min_valid_k
isf.bounds
isf.bin_edges
isf.bin_centers

isf.self_function
isf.distinct_function

The query_points and N_total in the compute method are only needed when computing partial scattering functions. It would be completely okay to ignore that and just compute full scattering functions in a PR that creates the class, then come back later and implement partials only if necessary.

@Roy-Kid How comfortable would you be trying to implement this yourself? It's not a simple calculation to write, but I would be able to assist to get past roadblocks in implementation.

Thanks for your concern.
I have started to read the source code of StaticStructureFactorDirect and msd and hope to combine those together. I think the difficulty is to understand the code about how to contracture reciprocal space.

I am try to do this work, but i'm not sure I can complete or give a robust code.

Okay. The sampling of reciprocal space is done in StaticStructureFactorDirect::reciprocal_isotropic, that code should actually be reused for this new feature, because the sampling method will be the same.

When I have the time, I can assist with getting some of the boilerplate code and class structure in place first and then we can work more on implementation from after that.

@tommy-waltmann
Hi waltmann,

I start a PR as a start for intermediate scattering calculation. I write boilerplate code for both python and c++. Here is my plan:

1. We derive/extend the StaticStructureFactorDirect class to use reciprocal_isotropic compute_F_k etc.
2. Since the isf is time-dependent and reply on r(t=0), so we record the position of first frame, and construct a neighborlist by AABQUERY((box, point)).query(r(t=0)).
3. The neighborlist computes the all distance between all particles and there should be a method to separate "self-distance" and "distinct-distance" (I think it can be done by the construct a index filter).
4. When we get the r_i_t - r_i_t0 and r_i_t - r_j_t0, we can reuse the compute_F_k to compute exponent factor.
5. To compute F(k, t), we need to pass all frames to compute method in python like what MSD does. The inner for loop accumulate the isf result one-by-one frame.
6. Finally we get a (n_frame, n_k) array for both self and distinct function.

What do you think? Please indicate any physical errors or code-style differences.