ISISNeutronMuon / MDANSE

MDANSE: Molecular Dynamics Analysis for Neutron Scattering Experiments
https://www.isis.stfc.ac.uk/Pages/MDANSEproject.aspx
GNU General Public License v3.0
18 stars 4 forks source link

[BUG] DCSF assumes f(q,t)_AB = f(q,t)_BA is this correct? #464

Closed ChiCheng45 closed 2 weeks ago

ChiCheng45 commented 3 weeks ago

Description of the error In the DSCF job the following correlation is done.

    def combine(self, index, x):
        """
        Combines returned results of run_step.\n
        :Parameters:
            #. index (int): The index of the step.\n
            #. x (any): The returned result(s) of run_step
        """

        if x is not None:
            for pair in self._elementsPairs:
                corr = correlation(x[pair[0]], x[pair[1]], average=1)
                self._outputData["f(q,t)_%s%s" % pair][index, :] += corr

The pairs in self._elementPairs are unique pairs. In the weight section of the code it uses a symmetric setting so that it weights it according to the fact that the combination AB and BA occur. However, cross-correlations are not symmetric in general.

Additional details Here I calculate the DCSF with water trajectory and the grid qvector generator. I use this generator so that the qvectors aren't random.

image

Here I switch the pairs around in the code e.g.

corr = correlation(x[pair[1]], x[pair[0]], average=1)

image

lucas-wilkins commented 3 weeks ago

I think things are calculated correctly, The structure factor shouldn't be symmetric in general, which you can see from the definition of the van Hove function:

https://en.wikipedia.org/wiki/Dynamic_structure_factor#The_van_Hove_function

If you exchange the r, and r+r' terms you don't have the same equation, you kind of get one going backwards in time. I guess at equilibrium we might expect the two to match, but that should be something that happens in the data, not the calculation.

gonzalezma commented 3 weeks ago

Lucas is right: F_ab(Q,t) = F_ba(Q,t). As every time you run a calculation the number of vectors (if > Nmax) is chosen randomly, nothing ensures that you will get exactly the same result. I suppose this could explain why you got two different curves. Otherwise, some checking would be needed.

ChiCheng45 commented 2 weeks ago

I used the grid qvector generator which is not random. I think the difference is due to the small number of configurations that are average over at large t, see #435.