Sequential MC

[wiederm]

Description

Outline of the basic Sequential Importance Sampling (SIS) algorithm as outlined in 10.1021/acs.jctc.1c01198

Goal

Calculate free energy estimate between reference distribution \pi(\lamb=0; x) and target distribution \pi(\lamb=1; x) with resampling steps

Motivation

SMC might help focus the non-equilibrium sampling on the region of interest; it might help with RBFE calculations in analogy to the work of 10.1101/2020.07.29.227959 and avoid sampling of the target distribution \pi(\lamb; x) for bidirectional NEQ protocols

Todos

• [x] SMC algorithm
• [x] Unit tests
• [ ] Integration tests
• [ ] Update script

Questions

• [ ] Naive implementation, can be made much faster in vacuum. Should we spend effort on this?

Status

• [ ] Ready to go

[codecov-commenter]

Codecov Report

Merging #77 (82fa329) into main (08e6556) will decrease coverage by 1.37%. The diff coverage is 50.00%.

:exclamation: Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

[xiki-tempula]

@wiederm Thanks for inviting me to review this PR. I have one question, I think the importance sampling based on energy makes sense if the potential energy of lambda 0 and lambda 1 make sense in the absolute scale. For endstate state correction, this doesn't make sense as the energy difference between lambda 0 and lambda 1 includes the energy difference between the real lambda 0 and lambda 1 but is dominated by the offset between the ML model and the MM model, which is the reason that one needs the vacuum simulation to correct the offset. So if one samples the lambda bases on the offset-included energy difference, this would not give you the real distribution.

[wiederm]

Thanks for the comment @xiki-tempula ! I don't agree with your conclusion though :-)

To explain why let's consider the following example: $U(x; lamb=1)$ is the target potential (a NNP) and $U(x; lamb=0)$ is the reference potential (representing a MM force field). Then the importance weight for $x$ is calculated with some variation of $w_x \propto exp(U(x; lamb=0) - U(x; lamb=1))$. Let's consider what happens to the distribution of $w$ if we add a constant $\alpha$ to the above equation: $w_x \propto exp(U(x; lamb=0) - (U(x; lamb=1) + \alpha))$ This will change the value of $w_x$, but it won't change the distribution of $w$ for a set of $x$. The importance weights are not influenced by some arbitrary offset of any of the potential energy functions --- and we can subtract the offset without changing its distribution.

[wiederm]

Updated SMC implementation after discussions with @msuruzhon and @maxentile

[wiederm]

This PR also contains some rearrangement of functions. The relevant file changes are in smc.py and test_smc.py

[wiederm]

thanks @msuruzhon for all your very helpful comments! I have started with the trivial changes and will add the stratified resampling tomorrow.