Add widget for generic ACF

[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 #33672403 | 📁 Comparing 03162f5 against latest (c213877)

  🔍 Preview build  

12 files changed · + 2 added · ± 10 modified

+ Added

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

± Modified

• genindex.html
• getting_started.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
• _modules/mdadash/backend/kernel/core.html
• _modules/mdadash/backend/state/manager.html

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (50b0525) to head (03162f5).

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. But there are some issues that might need a little attention to detail.

As for the current ACF implementation with centering, there might be an issue with how it's implemented. It looks good to subtract the mean in place, but since these are running updating means, doing so will bake the current value into the product and thus the overall ACF average. It would be more accurate thus to add the mean contribution separately as suggested here,-Here%2C).

Now, as for the plots, could you share the input scripts of your test system so I can double-check the results—

• Btw 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 a non-zero mean).
• The Force ACF also looks good and has an expected exponential-like decay. But the decorrelation time is too short. I'll see if there is any reference data I can compare it against.
• 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. I hope you are using the no-jump transformation here as well.

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 * ACF(0) - 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!

[HeydenLabASU]

just in case, I performed a reference calculation with a stored trajectory for the position, velocity and force ACF (without subtraction of the mean):

https://www.dropbox.com/scl/fo/5eihvfegx1ebcm801q825/ANwykkNWJl2wvE79c-hoBFM?rlkey=ov9chkpnuhqvfv767d08mvuj0&dl=0

The simulated system should be the same as in your tests: lysozyme with 1960 atoms in water + 150mM salt (all files in MD sub-directory).

I stored coordinates, velocities and forces every 0.01 ps and the average atomic correlation functions for coordinates, velocities and forces (computed with 200 correlation steps, max correlation time = 1.99 ps) are in the files "acf_results.dat".

I cross-checked the python code in acf.ipynb against a C-code that I wrote a while ago.

The results looks pretty much like yours, which is great! Only the coordinate correlation functions look different, because I remove jumps across periodic boundary conditions.
In any case, I hope this may be helpful.

[PardhavMaradani]

Hi @amruthesht, @HeydenLabASU, thank you both for additional details. They are very helpful.

The artifact in the Postion ASF in my previous comment at the very beginning was because the running count was initialized to 1 (a carry over from the MSD code). I now fixed this and this looks much better.

I am now using the system provided by MH in the comment above. I made minor modifications to stream that simulation to test with our dashboard. (Just for reference, I was so far using the sample gromacs sim from imd workshop 2025) In the dashboard, I changed the total steps to 10000 and the buffer size to 200 to get the 2ps lag.

I got an exact match to the plots that MH had in acf.ipynb as shown below:

ACF's from acf.ipynb:

However, I got the above match only after removing the nojump transformation I had for the universe (u.trajectory.add_transformations(NoJump())) during it's creation. With the nojump transformation added, the position ACF graph was fluctuating quite a bit during the sim and eventually even after the full sim, didn't match what MH had. I thought that the nojump was needed here and would be equivalent to the additional handling for coords in MH's code (c_pbc)?

Update: After trying this out a bit more, I can confirm that the code in this PR is equivalent to MH's reference. The code in acf.ipynb still leaves out a bunch of lags uncomputed - for the frames left in the buffer after the main trajectory loop is done. Once that is accounted, I see exactly the same plots. After adding the nojump transformation to both cases, I get the same plots in all cases. So, overall, I think we are good here.

Here is the screenshot with all the three physical properties, with centering and normalization:

As for the current ACF implementation with centering, there might be an issue with how it's implemented. It looks good to subtract the mean in place, but since these are running updating means, doing so will bake the current value into the product and thus the overall ACF average. It would be more accurate thus to add the mean contribution separately as suggested here,-Here%2C).

I didn't change the way centering was implemented from the previous comment. Based on the plots above, should I retain this as is or change? (The plots looked fairly consistent all through this shorter sim - after the window filled up)

Now, as for the plots, could you share the input scripts of your test system so I can double-check the results—

I will use MH's sim going forward as it is something I can run till the end and also check with the trajectory file if needed for any future cases.

As for the positional ACF, could you run this for longer and for about $10 \ ps$ worth of lag times?

I am assuming this is no longer needed. If needed, please let me know and I can run this.

I hope you are using the 'no-jump transformation' here as well.

Please see above about what I ran into with 'no-jump' added/removed. Please let me know if there is anything you'd want me to check here.

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).

I've for now updated the PR with the generic ACF changes (with centering support, etc). I will try out how best to add cross-correlations and see how it goes - I think I understood what is required. The changes to MSD to add an additional option on how it is calculated is something I can explore after that. Hope that is ok.

Please let me know if there is anything else you want me to check.

Thanks

[amruthesht]

@PardhavMaradani, thank you for the new changes. They look great and are in the right direction in making the analysis more generic and robust!

As for the nuances,

1. I didn't change the way centering was implemented from the previous comment. Based on the plots above, should I retain this as is or change? (The plots looked fairly consistent all through this shorter sim - after the window filled up)

You should still change the centering to the modified mean subtraction, (with a note that the running-updated mean based on all data processed so far is used; and a small warning message on stats),-Here%2C) I suggested earlier. Here, the difference is not pronounced as the means are negligible in the simulation. However, in cases where it does matter, the difference might be more pronounced.

2. However, I got the above match only after removing the nojump transformation I had for the universe (u.trajectory.add_transformations(NoJump())) during it's creation. With the nojump transformation added, the position ACF graph was fluctuating quite a bit during the sim and eventually even after the full sim, didn't match what MH had. I thought that the nojump was needed here and would be equivalent to the additional handling for coords in MH's code (c_pbc)?

Let me get back to you on this. In my opinion, there might be something weird happening with the transformation.

3. The changes to MSD to add an additional option on how it is calculated is something I can explore after that. Hope that is ok.

This is completely alright!

Great work again in making the changes—this looks in good shape and should be ready once the finer details are done!

[PardhavMaradani]

You should still change the centering to the modified mean subtraction, (with a note that the running-updated mean based on all data processed so far is used; and a small warning message on stats) I suggested earlier.

I changed the centering to the modified mean subtraction as suggested and also added the note to the description of this input so that it shows up in the UI. Thanks