Implement Gaussian-sampling Langevin thermostat

[J-Hizzle]

Whilst the most common and obvious implementation for the Langevin thermostat is with Gaussian sampling, the current Langevin thermostat in MRMD samples the uniform distribution centered around 0.

This is permitted due to the statistics of the two sampling schemes being equivalent in the limit of small time steps and more performant (see ref. [1]).

However, to remain comparable with Gromacs, a Gaussian-sampling Langevin thermostat would be desirable and would hence need to be implemented.

References:
[1] Dünweg, B., Paul, W., International Journal of Modern Physics C, Volume 2, Issue 03, pp. 817-827 (1991).

[XzzX]

Ah, good to know. Feel free to implement it. Perhaps switchable if possible.

[junghans]

I think we had this discussion for quite a bit for @espressomd, but ultimately one can never compare two MD codes due to all the other roundoff errors, it will always be stochastic and only averages matter. So I am really not sure I would give this a high priority.

[J-Hizzle]

I think we had this discussion for quite a bit for https://github.com/espressomd, but ultimately one can never compare two MD codes due to all the other roundoff errors, it will always be stochastic and only averages matter. So I am really not sure I would give this a high priority.

Thanks for the insight. Would this argument also hold if I were to use a Langevin thermostat in HY and CG regions of AdResS and no thermostat within AT?
I could imagine the Gaussian-sampling thermostat to be more stable in this case, as it would more faithfully enforce the correct velocity distribution at the AT/HY interface compared to uniform sampling.

[J-Hizzle]

Interestingly enough, the Gaussian-distribution-sampling Langevin thermostat appears to be already implemented within MRMD! It can be found in action/VelocityVerletLangevinThermostat.hpp.

Some comments on this matter:

• From what I can say, the implementation of the Gaussian-sampling thermostat in action/VelocityVerletLangevinThermostat.hpp corresponds to the BAOAB scheme used in many contemporary MD codes, which is very nicely explained in Allen and Tildesley's book (see Code 12.1 on p. 385 of ref. [1]) and Kieninger and Keller's paper (see eqs. (10a - d) and (11) in ref. [2]). This splitting scheme basically performs half a conservative velocity and position update (B + A) before applying the friction and noise terms in integrated form to the velocities directly (O), followed by another half a position and velocity (A + B) update.
• In contrast to that, the implementation of the uniform-sampling thermostat in action/LangevinThermostat.hpp follows the philosophy of applying the friction and noise terms directly to the forces and then applying those forces to the velocities. I suppose that this is just a different representation of the same scheme and more inspired from the equation given for the change in momentum at the bottom of p. 384 of ref. [1].
• I guess one may use either of these thermostats and it's great to have them both, but for ease of use it might be nice to unify them to one single thermostat module eventually.
• Also, having these two philosophies in parallel is confusing since I don't think they have anything to do with the sampled distributions per se. I think one could implement either distribution for sampling in either of these representations.
• I like the force-based philosophy of the uniform-sampling thermostat since it elegantly detaches the application of the thermostat from the rest of the integration.

[1] Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids: Second Edition, 002 edition.; OXFORD UNIV PR: Oxford, United Kingdom, 2017.
[2] Kieninger, S.; Keller, B. G. GROMACS Stochastic Dynamics and BAOAB Are Equivalent Configurational Sampling Algorithms. J. Chem. Theory Comput. 2022, 18 (10), 5792–5798.

[XzzX]

Thank you for all the references.
What do you suggest? We can decouple the distribution from the rest, that is no problem.
Do you want to be able to apply what is currently applied to the velocity in action/VelocityVerletLangevinThermostat.hpp to the force in action/LangevinThermostat.hpp?

[J-Hizzle]

Thank you for all the references. What do you suggest? We can decouple the distribution from the rest, that is no problem. Do you want to be able to apply what is currently applied to the velocity in action/VelocityVerletLangevinThermostat.hpp to the force in action/LangevinThermostat.hpp?

Yes, that would be very nice.