use correct Maxwell-Boltzmann velocity distribution in reaction ensemble

[thepith]

Fixes #1770 and #2377

The velocities of the particles inserted by the reaction ensemble did not follow the Maxwell-Boltzmann velocity distribution [1]. Even worse, with increasing temperature the mean squared velocity was decreasing. The effect of this PR on the velocity distribution function are shown in the following figure:

As you can see, now the measured velocity distribution from the simulation (marked as Sim in the plots).becomes wider with increasing temperature and follows the Maxwell-Boltzmann distribution (marked as Theo in the plots).

[1]: Atkins, P. W., and J. de Paula. Atkins’ Physical Chemistry. Oxford, UK: Oxford University Press, 2006.

[codecov]

Codecov Report

Merging #2377 into python will increase coverage by <1%.
The diff coverage is 100%.

@@           Coverage Diff           @@
##           python   #2377    +/-   ##
=======================================
+ Coverage      72%     72%   +<1%     
=======================================
  Files         397     397            
  Lines       18773   18778     +5     
=======================================
+ Hits        13580   13585     +5     
  Misses       5193    5193

Impacted Files	Coverage Δ
src/core/reaction_ensemble.cpp	72% <100%> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4ba6b58...1146eca. Read the comment docs.

[fweik]

Thank you for your contribution. This looks more correct than before to me, but finally @jonaslandsgesell has to decide.

[jonaslandsgesell]

This PR fixes #1770 and #2377 for particles of mass 1. For particles with different masses the Maxwell-Boltzmann distribution needs to take into account the particle mass which has to be provided to the reaction ensemble algorithm as input then. This is not yet fixed with this PR and remains an issue. I opened issue #2388 to document this.

[fweik]

I don't get it, ReactionAlgorithm::create_particle only creates particles of mass 1, no?

[jonaslandsgesell]

Yes but this should be changed if one wants to be as correct as possible. One example would be if one wants to sample the (mass dependent) velocity distribution in a reacting system. Then the masses of all types need to be correct. One can also construct exotic cases where masses become a problem. This is obvious when you think about only the particle moving Monte Carlo moves (which are also implemented in Espresso in the reaction ensemble module). Consider a case where a user created particles with mass 0.5 but has not assigned velocities. If there is no other thermostat running the velocities will stay 0 if we only use MC moves to sample the configuration space. The configuration space is sampled correctly. However the kinetic part of the virial pressure will be 0 which is not how a termalized system would behave. In order to solve this problem the correct velocity distribution (which is mass dependent) needs to be assigned to the particles... All in all this is not a problem if the system becomes correctly thermalized using e.g. a Langevin thermostat before taking a sample. Am Mo., 26. Nov. 2018 um 18:53 Uhr schrieb Florian Weik < [email protected]>:

…

I don't get it, ReactionAlgorithm::create_particle only creates particles off mass 1, no? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2377 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AE8wPf9NnjMYYkjEQi1CWjp0AZseVNKTks5uzCqTgaJpZM4YgCuP> .

[fweik]

I mean you could just draw the velocity from the correct MBD for the particle, and this would be correct for ever right? Why can't you just put in the mass?

[jonaslandsgesell]

Yes exactly. This will be the way to go. Pascal suggested to postpone this though since it requires the Reaction Ensemble to be aware of masses (User input). Am Mo., 26. Nov. 2018, 19:43 hat Florian Weik <[email protected]> geschrieben:

…

I mean you could just draw the velocity from the correct MBD for the particle, and this would be correct for ever right? Why can't you just put in the mass? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2377 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AE8wPR8M1WD6Rr1HaiEpKeXIcZFrHhnAks5uzDZsgaJpZM4YgCuP> .

[thepith]

Or perhaps to rephrase myself: I currently do not have the capacity to do it. If @jonaslandsgesell wants to do it, we can incorporate it in this pull request.

Maintainers are allowed to edit, so I believe Jonas can push to my branch...