choderalab / ensembler

Automated omics-scale protein modeling and simulation setup.
http://ensembler.readthedocs.io/
GNU General Public License v2.0
52 stars 21 forks source link

Default Collision Frequency #5

Closed kyleabeauchamp closed 9 years ago

kyleabeauchamp commented 10 years ago

So it looks like the default Langevin collision rate is set to 5.0 per picosecond (https://github.com/choderalab/msmseeder/blob/master/MSMSeeder/packaging.py#L114). According to Michael's paper (http://pubs.acs.org/doi/abs/10.1021/ct400109a), I think this might be slightly too strong. I think a better value might be around 0.5 per picosecond.

kyleabeauchamp commented 10 years ago

Or even 0.25 per picosecond. See figure 1 in that paper.

jchodera commented 10 years ago

Could you email a copy of the paper?

kyleabeauchamp commented 10 years ago

Done

jchodera commented 10 years ago

Look at Fig 5b. If I'm reading this right, it looks like weak thermostats actually can greatly reduce the decorrelation rate (the inverse of the correlation time).

I agree that our default of 5/ps is essentially an arbitrarily chosen value that could be optimized, and it may be too large. But we should actually do a domain-specific test on, say, hydration free energies or host-guest binding affinities computed using a variety of collision rates (if we're evaluating what to use for YANK, for example).

I'm not sure if Michael's paper mentions this, but there is actually a nice physically-motivated way to select the collision rate for systems in explicit solvent that Hans Andersen described in his 1980 thermostat paper---see Eq. 4.9. Given the thermal conductivity of the solvent (water) and number of particles in the system, we can derive (up to a shape-dependent constant that need only be computed once for each box shape) the appropriate collision rate to preserve the thermal conductivity. This provides a useful physical estimate to obtain appropriate kinetics.

jchodera commented 10 years ago

If we want realistic kinase kinetics, we should probably do an actual test for our water model and various box sizes to determine the Andersen parameter a.

kyleabeauchamp commented 10 years ago

IMHO I think you have it backwards: weaker coupling increases the decorrelation rate. Note that Fig. 5a reports timescales (stronger coupling = larger timescales) while Fig. 5b reports rates (weaker coupling = larger rates).

By my interpretation, Figure 5 suggests that by choosing a value closer to (10 ps) ^-1 = 0.1 per picosecond would give us both benefits:

  1. Closer to NVE results
  2. Faster decorrelation
kyleabeauchamp commented 10 years ago

FWIW, picosecond-timescale coupling is also in line with Gromacs "best practices", although that doesn't really mean anything: http://manual.gromacs.org/current/online/mdp_opt.html#run

jchodera commented 10 years ago

FWIW, picosecond-timescale coupling is also in line with Gromacs "best practices", although that doesn't really mean anything: http://manual.gromacs.org/current/online/mdp_opt.html#run

I frankly don't even trust the Langevin implementation they have in gromacs, so we should be sure to look at only results that are consistent between all the stochastic thermal control schemes they present

jchodera commented 10 years ago

IMHO I think you have it backwards: weaker coupling increases the decorrelation rate. Note that Fig. 5a reports timescales (stronger coupling = larger timescales) while Fig. 5b reports rates (weaker coupling = larger rates).

Fig 5b shows the "decorrelation rate" in 1/ps, which is defined in the lower-left of pg 2893 as 1/\tau_B. Larger correlation times \tau_B mean there is more simulation time required to generate an uncorrelated sample; hence larger decorrelation rates 1/tau mean that less simulation time is required to generate an uncorrelated sample. If you zoom in on Fig 5b, the largest decorrelation time is for the double-hatched lines, which uses an inverse collision rate of 10 ps, the slowest collision rate.

Unless Michael's figure is fucked up, I'm pretty sure this is the correct interpretation.

kyleabeauchamp commented 10 years ago

If you zoom in on Fig 5b, the largest decorrelation time is for the double-hatched lines, which uses an inverse collision rate of 10 ps, the slowest collision rate.

I disagree with this statement in particular. In Fig. 5b, the largest decorrelation rate is for the double-hatched lines. The decorrelation time is the inverse of the rate.

Or, another way to look at this is to just look at Fig. 5a. There, you can clearly see from the correlation functions that the NVE (black) decorrelates faster.

jchodera commented 10 years ago

I disagree with this statement in particular. In Fig. 5b, the largest decorrelation rate is for the double-hatched lines. The decorrelation time is the inverse of the rate.

There is no such thing as a "decorrelation time". It is a correlation time, and the larger the correlation time is, the longer the simulation must be run to generate an uncorrelated sample.

jchodera commented 10 years ago

Or, another way to look at this is to just look at Fig. 5a. There, you can clearly see from the correlation functions that the NVE (black) decorrelates faster.

You're right. Fig 5b is completely fucked up. Fig 5a is correct.

Someone should let Michael know he has to correct this paper.

jchodera commented 10 years ago

Sorry---Fig 5b is likely correct. The thermal control parameter that has the fastest decorrelation time (shortest correlation time) is \tau_T = 10 ps, or collision rate of 1/(10 ps) = 0.1/ps.

So yes, using very weak collision rates does give faster decorrelation and less perturbed diffusivities compared to NVE.

kyleabeauchamp commented 10 years ago

You're right though, there's a discrepancy between the caption text ("The rate of decorrelation, as measured by the time constant of a stretched exponential fit to S(k,t), is shown for all temperature control schemes tested") and the figure, which shows rates [ps ^ -1].

Basically, they say "time constant" but mean "rate constant".

kyleabeauchamp commented 10 years ago

That discrepancy probably explains why we couldn't agree faster.

kyleabeauchamp commented 10 years ago

Going back to the practical question: are there any objections to making 0.1 ps^-1 the default value? That would bring us closer to NVE and also give us better sampling.

jchodera commented 10 years ago

I'd be curious to investigate the thermal conductivity based approach. It would basically just require tuning the a parameter from the equation

nu = (2/3)*a*(kappa/k)*(rho**(1/3))*(N**(2/3))

where kappa is the thermal conductivity of water (0.58 W/m/K), rho is the density N/V, and k is the Boltzmann constant, and N is the number of atoms. a would be selected to minimally perturb dynamics. We could even compute this from the results in Michael's paper, presuming he used a cubic box.

jchodera commented 10 years ago

Going back to the practical question: are there any objections to making 0.1 ps^-1 the default value? That would bring us closer to NVE and also give us better sampling.

We should definitely not use this for equilibration with MSMSeeder, but it might be OK to use it for production dynamics (in the FAH packaging stage).

Do you think the OpenMM energy conservation is good enough for this really small collision rate?

kyleabeauchamp commented 10 years ago

I'm not sure how good it is, particularly on the crazy machines running on FAH. For example, I have a strong feeling that Gromacs would NOT be OK at 0.1 ps^-1

Let's do 0.25 ps ^-1 as a compromise? It should be very close to NVE but still "safe".

jchodera commented 10 years ago

Whatever the choice, I think we should test it out. 0.25/ps sounds like a reasonable choice for starting to test.

jchodera commented 10 years ago

We should probably take this opportunity to tune PME parameters too.

jchodera commented 9 years ago

It looks like the default collision rate for implicit solvent refinement is 20/ps. Is this OK for the paper?

We can certainly mention in the paper that lower collision rates may help, and cite the work above.

kyleabeauchamp commented 9 years ago

I'm happy. I brought this primarily for the eventual explicit solvent FAH simulations, which inherited values from here.