Reaction methods bug fixes #4589
Conversation
jngrad
force-pushed
the
reaction_bugfixes
branch
from
aaabe63 to
8400107
Compare
Restore original velocities when a displacement move is rejected. When one move in the sequence is rejected, skip the remaining move attempts. Fix infinite loop when not enough particles are available.
When no user-defined type is available, explicitly assign type 0. This is necessary for tracking particles ids in reaction methods.
jngrad
force-pushed
the
reaction_bugfixes
branch
from
8400107 to
292bf32
Compare
|
Thanks for involving me in this PR! There are some things that I am not convinced about the MC displacement move method:
@davidbbeyer do you think that such an implementation would still be compatible for your project? If yes and all of you agree with this alternative implementation I volunteer to change the code accordingly |
|
Offline discussion with @pm-blanco: the Monte Carlo displacement move method in its current state could benefit from a few improvements, but that work might take a couple of weeks to complete. There is also a newly discovered bug when combining constant-pH with Widom insertion, which may be relevant for Monte Carlo displacement moves too, since displacement moves can be carried out from any reaction method class instance. Since this PR fixes important bugs and is blocking the parallel Monte Carlo project, we agreed that this PR should be merged now. The other points can be taken care of in an upcoming PR that will be orthogonal to the parallel MC project. |
| m_tried_configurational_MC_moves += 1; | ||
| particle_inside_exclusion_range_touched = false; | ||
|
|
||
| auto const particle_number_of_type = number_of_particles_with_type(type); | ||
| if (particle_number_of_type == 0 or | ||
| particle_number_of_type_to_be_changed == 0) { | ||
| if (n_part == 0) { | ||
| // reject | ||
| return false; | ||
| } |
There was a problem hiding this comment.
Is this really intentional that we increase the counter of m_tried_configurational_MC_moves when n_part is zero? Did one really try a configurational move if one doesn't change anything?
There was a problem hiding this comment.
That's an oversight. Calling this method with 0 particles doesn't make sense, and it isn't meaningful to return true or false because no move was actually attempted. It's even worse if n_part is a negative number... I'll fix that.
| if (n_part <= 0) { | ||
| throw std::domain_error( | ||
| "Parameter 'particle_number_to_be_changed' must be >= 1"); | ||
| } | ||
| auto const n_particles_of_type = number_of_particles_with_type(type); | ||
| if (n_part > n_particles_of_type) { | ||
| throw std::domain_error( | ||
| "Parameter 'particle_number_to_be_changed' must be less or equal to " | ||
| "the number of particles in the system"); | ||
| } |
There was a problem hiding this comment.
We should spend a second thought on this change, because it will probably break old espresso scripts.
What I was thinking of is a reactive system (e.g. GCMC) where you have initially setup a species that will react with time. At some point it might happen that this species has completely reacted, but one might have used displacement moves in the script to change the configuration. In such a case one would have just rejected the move in the old implementation, but now an error is thrown. An additional, explicit check before the displacement moves would be necessary to prevent this error.
The very same happens when providing zero, which might happen when using the current particle number as the number of displacement moves.
Disclaimer: I'm not claiming any of this is good practice, but it used to work and I have written scripts where this change would also break my code.
jngrad
force-pushed
the
reaction_bugfixes
branch
from
cb58c0e to
7893235
Compare
Description of changes: - track particles with default-constructed type 0 (fixes espressomd#4588) - restore the original particle velocity when a Monte Carlo displacement move is rejected (fixes espressomd#4587)
Description of changes: