Gibbs ensemble MD moves instead of MC moves

[jonaslandsgesell]

Would you mind investigating using MD moves instead of MC moves to propagate the internal states of each system. Just instead of moving 1 particle and computing exp(-beta Delta E_pot) in an MC step, run x MD steps (e.g. 10-100) and compute exp(-beta Delta E_pot) from before and after running MD. This would allow to simulate polymer systems and should have much faster decaying autocorrelations of observables.

Running MD will require that you initialize to the velocities with a Maxwellian distribution like e.g. in https://github.com/espressomd/espresso/blob/ae3094f47cf78cfafaddb5cdc8c0e7edb5630b4f/src/core/reaction_methods/ReactionAlgorithm.cpp#L503 after creation or when you insert new particles etc.

You will also have to avoid inserting particles closer than e.g. 0.9 sigma to avoid instabilities in the MD (see https://github.com/espressomd/espresso/blob/ae3094f47cf78cfafaddb5cdc8c0e7edb5630b4f/src/core/reaction_methods/ReactionAlgorithm.cpp#L591 ) for fulfilling detailed balance you may then also not remove particles which were closer than this exclusion radius. You can just use the logic from within the reaction ensemble here.

Originally posted by @jonaslandsgesell in https://github.com/espressomd/espresso/issues/4243#issuecomment-837910703

[jhossbach]

Would it be necessary to check the energies before and after? Because we are running MD simulations using a thermostat we wouldn't need to check if the particle movement is correct, we could simply accept it as such. Or am I missing a point here?

[jonaslandsgesell]

You are right, if MD sampling is already correct (i.e. particles have Maxwell-Boltzmann distributed velocities etc.), then the additional MC criterion is not needed.

If you want to have guaranteed properties from the MC algorithm, you need to check energies before and after running MD (this is frequently done in hybrid Monte-Carlo methods https://de.wikipedia.org/wiki/Hybrid-Monte-Carlo-Algorithmus). The additional MC criterion helps that you obtain correct results (i.e. have correct sampling) even if the MD thermostat fails in sampling the correct pyhsical statistics if e.g. velocities are incorrectly assigned: for example with the additional MC criterion you can choose much higher timesteps in the MD part.

[jhossbach]

My algorithm seems to be working now, but I have problems with the particle exchange move. I use the logic from the snippet you proposed with d_min=0.9 sigma, it seems that almost every time a random particle is chosen or a random position is proposed, at least one of them will be in close range to another particle resulting in no particle exchanges being done in the overall simulation. Even for small initial densities this seems to be a problem, but even if it wasn't this would be a huge drawback for this algorithm. Is it possible to use different approaches to avoid instabilities, like performing a few MD steps in both system boxes to "shake it off" or use force caps?

[jonaslandsgesell]

Hey, I see your point! The Gibbs-Ensemble is applied here in the context of systems with high densities, where employing the exclusion radius based rejection strategy can become quite inefficient. I like your idea of trying to employ force capped, thermalized MD to propagte the state of the systems and then accept/reject the move with exp(-beta Delta E_pot). In the case of a rejection of the MD moves, the complete state prior to the attempted move needs to be restored. Potentially, one also needs to be careful to not introduce asymmetric proposal probabilities due to the force-capping, maybe there is literature already?

PS: In low density systems (e.g. concentrations below 1 mol/l), the above exlusion radius based strategy was successfully used for investigating the Donnan partitioning of ions between a polymeric gel and the supernatant solution.

[jhossbach]

I found that this is a common problem with Gibbs Ensemble simulations. This paper describes two ways to implement MD into the simulation: One uses pure MD simulations but introduces an extra degree of freedom describing whether the particle is in either one of the boxes or in a transition state inbetween to simulate the Gibbs Ensemble behaviour. But as of yet I don't see exactly how to implement this in ESPResSo. The other one, the so called Hybrid Monte Carlo/Molecular Dynamics seems to do the same as you proposed. For now I will concentrate one the latter algorithm.

[jonaslandsgesell]

Thanks for the paper :) A key concept there is summarized like this : "we implemented a short- range cutoff to the pairwise potential to remove instabilities resulting from the divergence of most intermolecular potentials at low separation distances. " This is promising for dense systems.

[jhossbach]

I opened a new issue dedicated to this short range cutoff in #4254, meanwhile I will use tabulated potentials and a LJ potential to mimick this behaviour.