Add widget for VACF

[PardhavMaradani]

Changes made in this Pull Request:

• Add a new SlidingWindowVACF class that is O(n) for each window
• Support to compute and display particle VACFs
• Support for displaying running integral of VACF
• Support for configuring selection, dimension type and normalization
• Add new test trajectory data file that supports velocities

PR Checklist

• Tests?
• Docs?
• CHANGELOG updated?
• Issue raised/referenced?

---

Hi @orbeckst, @jeremyleung521, @amruthesht, @HeydenLabASU,

This PR adds a new 'VACF' widget as part of issue #34.

The code for VACF follows exactly the same lines as that of MSD support added in PR #35. I used this reference from MDAnalysis/transport-analysis for the VACF computations.

Here is an example showing VACF for 'all' and 'protein'. There is support to show normalized values and also particle VACFs if required as shown below:

Based on the reference, I also added an option to show the running integral if configured as shown below:

With PR #35 and this, we will have support for both the lag-time based computations mentioned in issue #34 .

Thanks

[read-the-docs-community]

Documentation build overview

📚 mdadash | 🛠️ Build #33637561 | 📁 Comparing ed82877 against latest (c213877)

  🔍 Preview build  

9 files changed · + 2 added · ± 7 modified

+ Added

• autosummary/mdadash.backend.analyses.vacf.html
• _modules/mdadash/backend/analyses/vacf.html

± Modified

• genindex.html
• py-modindex.html
• _modules/index.html
• autosummary/mdadash.backend.analyses.html
• autosummary/mdadash.backend.analyses.msd.html
• autosummary/mdadash.backend.kernel.core.html
• _modules/mdadash/backend/analyses/msd.html

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (c213877) to head (ed82877).

Additional details and impacted files

Components	Coverage Δ
frontend	100.00% <ø> (ø)
backend	100.00% <100.00%> (ø)

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[PardhavMaradani]

Given we have an option to display the plot of running integral of VACF, I added an option in MSD to display the diffusion coefficient plot from MSD values (using local derivatives).

Here is an example showing how the diffusion coefficient value is in agreement between both these plots:

[amruthesht]

This looks great, @PardhavMaradani. If it is not too much, I would suggest you add some of the following features, if feasible:

1. Make this a generic autocorrelation-based analysis and have support to choose the physical property being analyzed. For example, you can have VACF, a positional autocorrelation, or even a force autocorrelation, all of which follow the same mathematical formalism.

2. It might be nice to offer the option to do this computation with the mean subtracted and the ACF being normalized in the final representation/results.
   i. For mean subtraction, you have

$ACF_{\text{centered}}(\tau) = \langle (v(t)-\mu) (v(t+\tau) - \mu) \rangle$
$ACF_{\text{centered}}(\tau) = \langle v(t) (v(t+\tau)\rangle + \mu^2 - 2 \langle v(t) \rangle \langle v(t+\tau) \rangle$
$ACF_{\text{centered}}(\tau) \approx ACF_{\text{raw}}(\tau) - \mu^2$ (when $N_{\text{samples}} \gg \tau$ - can be added as a warning that the number of data samples must be much greater than the lag-time window for this to be accurate)

Here, $\mu$ would simply be the running mean calculated as data arrives and can be updated on the plot accordingly.

ii. Normalization usually involves normalizing the ACF by its variance or its value at $\tau=0$. I do see you have a normalization parameter—is that accessible as a toggleable switch from the UI interface or something like that?

3. Finally, I see that you have particle-wise VACFs, which is very nice and useful. It would also be very useful to have cross-correlations between particles of, say, a certain selection. This would be a costly computation, so it would be better to have an optional flag to trigger this. As for the mean subtraction for cross correlations, it will involve a $\mu_i \mu_j$ term instead of $\mu_i^2$.
   The actual plotting or representation of this can be tricky (it can be too crowded without a proper legend) and is not always necessary. But it would be nice for the user to have access to these calculations as a part of the final result from the loop.

Finally, a related suggestion is to have an ACF function that takes in two selections and computes cross-correlations between them.

I know this might involve a lot of code additions, so take your time and let me know if you have any questions/need any help with these code implementations. I would also be glad to contribute to some of these if you need any help!

[PardhavMaradani]

Hi @amruthesht, thanks for reviewing the changes.

1. Make this a generic autocorrelation-based analysis and have support to choose the physical property being analyzed. For example, you can have VACF, a positional autocorrelation, or even a force autocorrelation, all of which follow the same mathematical formalism.

Yes, this is feasible. I quickly tried this out, please see next section for some results...

2. It might be nice to offer the option to do this computation with the mean subtracted and the ACF being normalized in the final representation/results.

$ACF_{\text{centered}}(\tau) = \langle (v(t)-\mu) (v(t+\tau) - \mu) \rangle$
Here, $\mu$ would simply be the running mean calculated as data arrives and can be updated on the plot accordingly.

I tried with this option by doing something like this:

...
        current = getattr(self.ag, self.physical_property)[:, self._dim]
        self.running_sum += current
        self.running_count += 1
        mu = self.running_sum / self.running_count
        for i in range(n):
            lag = n - 1 - i
            _ = self.u.trajectory[i]  # set trajectory to past frame
            previous = getattr(self.ag, self.physical_property)
            if self.centered:
                corr = (current - mu) * (previous[:, self._dim] - mu)
            else:
                corr = current * previous[:, self._dim]
            sum_corr = np.sum(corr, axis=-1)
            self.acf_sums[lag] += np.mean(sum_corr)
            self.acf_counts[lag] += 1
...

I got plots that looked like this:

The plots labelled 'Centered' are the ones with the mean subtracted and the 'Normalized' ones are the normalized ACFs.

All of them run through the same code pasted above. The VACF looks the same as before. But the ones for positions and forces don't seem right, but I don't know much and don't know if they are expected / correct. Could you please check and let me know.

ii. Normalization usually involves normalizing the ACF by its variance or its value at $\tau=0$. I do see you have a normalization parameter—is that accessible as a toggleable switch from the UI interface or something like that?

Yes, this is already available for customization from the UI. Any attribute that shows up in the _inputs array in the widget definition can be configured from the UI.

Here is an example (includes some of the changes suggested here):

We are currently dividing by the value at 0 as follows:

...
        if normalized:
            avg_acfs = avg_acfs / avg_acfs[0]
...

3. Finally, I see that you have particle-wise VACFs, which is very nice and useful. It would also be very useful to have cross-correlations between particles of, say, a certain selection. This would be a costly computation, so it would be better to have an optional flag to trigger this. As for the mean subtraction for cross correlations, it will involve a $\mu_i \mu_j$ term instead of $\mu_i^2$.
   The actual plotting or representation of this can be tricky (it can be too crowded without a proper legend) and is not always necessary. But it would be nice for the user to have access to these calculations as a part of the final result from the loop.

Finally, a related suggestion is to have an ACF function that takes in two selections and computes cross-correlations between them.

Could we add cross-correlations as a separate widget? The current one (if we rename as ACF instead of VACF) takes a single selection phrase. Adding another selection here might make it a bit confusing / complex code wise. A separate widget could make this a clean separation both in terms of what to expect and how we structure it?

Could you please let me know if I should go ahead with the generic ACF (for velocities, positions, forces) and mean subtracted (centred) option based on your review of the plots above? (I didn't commit the changes because I wasn't sure if the plots are correct / expected)

[amruthesht]

@PardhavMaradani - thanks for making the requested changes.

The plots look pretty good and reasonable.
The VACF remains similar, as velocities on average have zero or close to zero mean by construction in MD (Newton's Laws) (even with NVT, it shouldn't be too huge of a non-zero mean).
The Force ACF also looks good and has an expected exponential-like decay.
As for the positional ACF, could you run this for longer and for about $10 \ ps$ worth of lag times? It looks reasonable, but I would like to double-check just in case.

Could we add cross-correlations as a separate widget? The current one (if we rename as ACF instead of VACF) takes a single selection phrase. Adding another selection here might make it a bit confusing/complex code-wise. A separate widget could make this a clean separation, both in terms of what to expect and how we structure it?

I would leave that up to you.
I was thinking of having a generic correlation widget/function that can perform cross-correlations between any 2 selections (and, without a second selection, perform a correlation with itself by default). This would include all the correlation terms (as a rectangular or square matrix, depending on whether the two selections have different numbers of atoms). Depending on input flags, it can restrict itself to just the diagonal (self) terms, if needed (can be used for particle-wise ACFs). This function will also be able to normalize and subtract the mean as flags. The ACF widget/function and cross-correlation widgets can then call this correlation function as needed for relevant calculations.
Ideally, even the MSD widget can use the same correlation function in the background by decomposing the MSD calculation as follows:

$MSD(\tau) = < |r(t + \tau) - r(t)|^2 >$
$MSD(\tau) = < r(t + \tau)^2 > + < r(t)^2 > - 2 * < r(t + \tau) · r(t) >$
$MSD(\tau) = 2 * < r^2 > - 2 * ACF(\tau)$

Here, however, the convergence of both will slightly depend on how well the system is equilibrated, so maybe you can have it as an option as to how the user wants MSDs calculated (using correlation or explicit positional subtraction)

Apologies for making this progressively more complex and suggesting different implementation approaches. I would leave the final decision to you on how you implement this overall. Whether as separate widgets/functions or otherwise.

But great implementation and progress so far on all the correlation and lag-time-based analyses! Everything looks robust and well implemented to me, and I'm sure it will be super useful for real-world analysis use cases!