Simulated tempering

[Peter Eastman]

I've written a script for doing simulated tempering. You can get it from the add-on repository at http://wiki.simtk.org/openmm/VirtualRepository. If you decide to try it out I'd really appreciate your feedback, both about the interface and about how it works.

Peter

[Kyle Kihn]

Hi Peter,
As you have probably guessed by some of my other forum posts, I have been running some systems using this simulated tempering script. So far it has worked really well for our lab, and seems to be way less finicky than the dual amd boost method. Once we get to a point where we can publish some of this data, what is the best way to site you and this code?
Thank you,
Kyle

[Peter Eastman]

Thanks! I'm glad it's working for you. The best way to cite OpenMM is to cite the first publication listed at https://simtk.org/plugins/publications/ ... blications (the 2017 one).

[lewis martin]

Thanks for this. I've had great initial success using just the default settings.

How would you recommend restarts are handled? Can the final weights from simulation_1 be used as the initial for simulation_2, for example? Cheers
lewis

[Peter Eastman]

Assuming the weights are well converged, and assuming you're simulating an identical system in a similar range of configurations, you can reuse them. If you've changed your system, or if you're exploring a very different region of configuration space, they might not be appropriate. In that case, the simulation will spend a lot more time at some temperatures than others.

[lewis martin]

Hi Peter,
Would you mind explaining the difference between `self.weights` and `self._weights` ? Looking at the code it seems that self.weights is just relative to the bottom temperature weight. Any time I do a restart using the last line of weights that were printed to the reporter file, my simulations spend all their time at the top temperature (450). The system has no conformational change, it's just protein and ligand in bulk at this time, and used 40ns equilibration with no initial weights set and thus had weights determined by the Wang-Landau algorithm which appeared converged. So I looked at `simulatedtempering.py` and found out that it is self.weights that gets printed to file, but perhaps self._weights is the appropriate one to use as weights for a restart.

I don't yet understand the difference in usage between the two, but for an example of what I'm looking at, i added these two lines to `_attemptTemperatureChange`:

Code: Select all

        print('weights:', self.weights)
        print('weights2:', self._weights)

and these print something like:

Code: Select all

weights: [0.0, 4095.0, 8191.0, 12287.0, 16383.0, 18431.0, 22527.0, 26623.0, 30719.0, 34815.0, 36863.0, 40959.0, 43007.0, 47103.0, 49151.0]
weights2: [-6144.0, -2049.0, 2047.0, 6143.0, 10239.0, 12287.0, 16383.0, 20479.0, 24575.0, 28671.0, 30719.0, 34815.0, 36863.0, 40959.0, 43007.0]

Thanks for your time.
lewis

[lewis martin]

I'll post this separately because it's a different question:

Would it be appropriate to turn on weight updates, even when the weights are provided?
For example, the __init__ method currently has:

Code: Select all

        if weights is None:
            self._weights = [0.0]*numTemperatures
            self._updateWeights = True
            self._weightUpdateFactor = 1.0
            self._histogram = [0]*numTemperatures
            self._hasMadeTransition = False
        else:
            self._weights = weights
            self._updateWeights = False

whereas I would change it to:

Code: Select all

        if weights is None:
            self._weights = [0.0]*numTemperatures
            self._hasMadeTransition = False
        else:
            self._weights = weights
            
        self._updateWeights = True
        self.weightUpdateFactor = 1.0
        self._histogram = [0]*numTemperatures

From the perspective of the concatenated trajectory, this would have the effect of resetting the weightUpdateFactor to 1 at each restart but otherwise I think it's the same.

[lewis martin]

Last one !
The Wang Landau paper states that the weight update factor 'f' is > 1, and in that paper it is updated as (f)^0.5, which would converge to 1 in the limit. What's the reason this script uses f*0.5 instead, which converges to 0? I'm not sure what effect it has. That paper also says 'The MC steps needed for a given f generally increase as we refine the modification factor', and if f is converging to zero instead of 1 then this increase is going to be a lot larger.
Thanks
lewis

[Peter Eastman]

Going through your questions in order:

1. The values returned by self.weights just have a constant subtracted to make the first value equal to zero. This shouldn't make any difference. The weights can be interpreted as free energies, so only differences should matter. If you have a case where it does matter, I'm not sure what's going on!

2. Technically speaking, simulated tempering only produces correct distributions when the weights are constant. Ideally you should know what weights you want in advance, specify them at the start, and then never change them. But of course, you often have no idea what values will work well, so in that case you can let it work out weights automatically. The Wang-Landau algorithm is designed to converge asymptotically, so once it finds good values they pretty much stop changing. And it's best to throw out the initial part of the simulation when they were changing quickly from your analysis.

So basically the class is intended to be used in one of two modes. If you know good weights to start with, you specify them and they should never change. If you don't know good weights, you let it figure them out and you throw out the initial part of the simulation.

3. This script operates on free energies, while the paper works with a density of states. The former depends on the log of the latter, so instead of multiplying by f we add log(f). And taking the square root of f is equivalent to multiplying log(f) by 0.5.

[lewis martin]

Thankyou that makes a lot of sense.