rethinking thermostats

[XzzX]

Can there be a more general interface for thermostats and integrators?

[XzzX]

Are these two implementations equivalent @J-Hizzle ?

mrmd/mrmd/action/LangevinThermostat.cpp

Lines 41 to 43 in fd144c2

	force(idx, 0) += p1 * vel(idx, 0) * m + p2 * (randGen.drand() - 0.5_r) * mSqrt;
	force(idx, 1) += p1 * vel(idx, 1) * m + p2 * (randGen.drand() - 0.5_r) * mSqrt;
	force(idx, 2) += p1 * vel(idx, 2) * m + p2 * (randGen.drand() - 0.5_r) * mSqrt;

mrmd/mrmd/action/VelocityVerletLangevinThermostat.cpp

Lines 52 to 63 in fd144c2

	auto damping = std::exp(-zeta * dtm);
	vel(idx, 0) *= damping;
	vel(idx, 1) *= damping;
	vel(idx, 2) *= damping;
	auto sigma = std::sqrt(temperature / mass(idx) * (1_r - std::exp(-2_r * zeta * dtm)));
	// Get a random number state from the pool for the active thread
	auto randGen = RNG.get_state();
	vel(idx, 0) += sigma * randGen.normal();
	vel(idx, 1) += sigma * randGen.normal();
	vel(idx, 2) += sigma * randGen.normal();

[XzzX]

When using velocity verlet integrator, do you always apply the thermostat at that position in the integration phase? @J-Hizzle

mrmd/mrmd/action/VelocityVerletLangevinThermostat.cpp

Lines 52 to 63 in fd144c2

	auto damping = std::exp(-zeta * dtm);
	vel(idx, 0) *= damping;
	vel(idx, 1) *= damping;
	vel(idx, 2) *= damping;
	auto sigma = std::sqrt(temperature / mass(idx) * (1_r - std::exp(-2_r * zeta * dtm)));
	// Get a random number state from the pool for the active thread
	auto randGen = RNG.get_state();
	vel(idx, 0) += sigma * randGen.normal();
	vel(idx, 1) += sigma * randGen.normal();
	vel(idx, 2) += sigma * randGen.normal();

[J-Hizzle]

Can there be a more general interface for thermostats and integrators?

I think it's definitely worthwhile to unify the two available thermostats while maintaining the choice whether to use the Gaussian or the uniform distribution for the velocity sampling if that's what you mean. Apart from that, I'm don't have an immediate idea for further generalization.

[J-Hizzle]

Are these two implementations equivalent @J-Hizzle ?

I can't say for certain. The second implementation is definitely the BAOAB scheme of Leimkuhler and Matthews (see p. 385 of ref. [1]), which is the generally preferred method. I think it would indeed be nice to decouple it from the integrator, then.

[1] Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids: Second Edition, 002 edition.; OXFORD UNIV PR: Oxford, United Kingdom, 2017.

[J-Hizzle]

As a side note: The implementation in MRMD assumes a Boltzmann constant of 1. This should be documented.

[J-Hizzle]

When using velocity verlet integrator, do you always apply the thermostat at that position in the integration phase? @J-Hizzle

yes

[J-Hizzle]

Are these two implementations equivalent @J-Hizzle ?

I can't say for certain. The second implementation is definitely the BAOAB scheme of Leimkuhler and Matthews (see p. 385 of ref. [1]), which is the generally preferred method. I think it would indeed be nice to decouple it from the integrator, then.

[1] Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids: Second Edition, 002 edition.; OXFORD UNIV PR: Oxford, United Kingdom, 2017.

I would write the update scheme implemented in LangevinThermostat.cpp as:
$v_{n + 1/2} = v_n + \frac{dt}{2 m}\left(F_n - \gamma m v_{n - 1/2} + \sqrt{\frac{24 T \gamma m}{dt}} R_n\right)$
$r_{n + 1} = r_n + \frac{dt}{m} v_{n + 1/2}$
$v_{n + 1} = v_{n + 1/2} + \frac{dt}{2 m} \left(F_{n + 1} - \gamma m v_{n + 1/2} + \sqrt{\frac{24 T \gamma m}{dt}} R_{n + 1}\right)$

This closely resembles the BBK update scheme (see appendix in ref. [1]), but not exactly: mind the factor of $(1 + \delta t \gamma/2)^{-1}$ in the third line of ref. [1] or the different scaling factor within the square root.

There is some criticism in the literature concerning the BBK update scheme (see ref. [1]), so it might be preferable to opt for BAOAB anyways and get rid of this implementation for now.

At the moment I conclude that we should not use this thermostat due to (1) the apparent differences to BBK that might yield disturbed statistics and (2) the preferable statistics that BAOAB seems to yield in the first place.

[1] Leimkuhler, B.; Matthews, C. Rational Construction of Stochastic Numerical Methods for Molecular Sampling. Appl Math Res Express 2013, 2013 (1), 34–56. https://doi.org/10.1093/amrx/abs010.