Introduced a bugfix to solve the infinite loop in constant pH algorithm #4200

pm-blanco · 2021-04-06T09:37:15Z

When attempting to do a reaction step using the constant pH algorithm, if there are not enough reactants nor product particles, the program simply not does the operation. The apparition of an infinite loop is prevented by existing the reaction subroutine in case that neither direction of the reaction is possible.

Description of changes:

The program checks if at least one direction of the reaction is possible and otherwise exits the subroutine returning a 1. This value should be later processed so the user gets a warning message informing of this error (work in progress).

…n algorithm when any direction of the acid-base reaction is possible. The program checks if at least one reaction is possible and otherwise exits the subrutine returning a 1. This value should be later processed so the user gets a warning message informing of this error.

KaiSzuttor · 2021-04-06T10:14:44Z

src/core/reaction_ensemble.cpp

+  bool reaction_possible = false;
+
+  for (int reaction_id = 0; reaction_id < reactions.size(); reaction_id++) {
+    if (all_reactant_particles_exist(reaction_id)) {
+      reaction_possible = true;
+    }
+  }
+
+  if (reaction_possible == false) {
+    return 1;
+  }
+


isn't it possible to get this kind of logic implicitly? If no particle of the reactant type is found nothing should happen automatically without manually checking?

The problem that causes the infinite loop has its source in the while loop in lines 1521-1546. This while loop searches the particle ids that have a type that matches the reactant/product types. Of course, if you do not create such particle, then this loop runs for ever.
As you suggest, another possible solution would be to do this task without using a while statement so that if there is no reactant then the program does nothing. However, in my opinion, manually checking it has the advantage that we can integrate an exception error preventing the user from doing such simulation (since I cannot think where one would attempt to run the algorithm without reactants if it is not by mistake).

The point is that such situation should not occur. If it occurs, then it is an error in the system setup, so it cannot be a simple no-op. It should raise an error.

in my mind the reaction algorithm itself should be independent of the state of the system. I can also define a non bonded interaction of two non existing particle types without having an infinite loop. The reaction of reactants is also some kind of particle interaction, right? Why should it be special here?

All methods default to a no-op when not enough particles are present, except for WidomInsertion which throws an error. I don't immediately see why it should be an issue to not have enough particles of type HA if there is at least one particle of type A-. Could you please clarify this for us?
Also if you reverse the reactant types with the adduct types, you will get a bug in the other direction, if I understand e4104d2 correctly.

@jngrad thank you for your clarification, now I understand your point. Now I agree with you that the desired output would be a no-op so one could couple the acid/base reaction with other reactions

Indeed, @jngrad is right that all reactants could be consumed if different reactions are competing for the same product or reactant. However, I would discourage anybody from intentionally using RxMC in this regime. It suffers from serious finite-size artifacts if the number of (any type of) reacting particles is on the order of unity. It is related to the fact that particle numbers in a simulation have to be integers whereas most existing theories work with concentrations as continuous variables.

Note: SN2 refers to a specific reaction mechanism. A reaction of the type A- + BX <--> AB + X- is not compatible with the way currently reactions are implemented in Espresso. The current implementation allows to create or destroy particles at random positions but it does not allow to transfer a particle between two binding partners. It is also one of the points on our wishlist.

Be aware it definitely depends on what you mean:
Of course it is totally valid to create this reaction
A+B <-> C+D
Whit artbitrary A, B, C, D. I.e. also :
A- + BX <--> AB + X-
However all the types will be "primitive" beads i.e. not have bonds in the bead "BX". But it is completely fine for the "BX" type to have different lenard jones parameters (i.e. appear bigger) etc...
Maybe this would require having a type dependent exclusion radius for some more extreme variations...

A reaction of the type A- + BX <--> AB + X- is not compatible with the way currently reactions are implemented in Espresso. The current implementation allows to create or destroy particles at random positions but it does not allow to transfer a particle between two binding partners.

Good point. I also thought about the case of the titration curve of a buffer solution, i.e. a mixture of a weak acid and a strong acid. Titration by a strong base can trigger the bug once the buffer capacity is used up. But that's an edge case, one probably wouldn't want to add more base than there is buffer.

… occur because there are not enough reactants

…ticles

RudolfWeeber · 2021-04-06T13:53:50Z

In my view, there are two issues:

The function ConstantPHEnsemble::get_random_valid_p_id() is incorrect, because it contains an empty loop if it has no solution. This should be fixed by changing the return type to std::optional<int> and returning {} if the system is empty. The quickest fix is to check that the total number of particles is >0. Long-term, the function can be replaced by something more efficient for the case that there are many holes in the particle id list.
The question of what should happen if no adducts are present, more more technically, what should happen in do_reaction() if get_random_valid_pid() returns {}. In my view, no-op is the correct behavior. This is analogous to running

system.cell_system.skin = 0
system.time_step = 1
system.integrator.run(1)

without particles (no-op, as-well)

RudolfWeeber · 2021-04-06T13:59:15Z

Espresso also allows definition of non-bonded interactions for types for which there are (at that point) no particles.
The interaction will become used as soon as relevant particles are there.
So, there as well no-op would be the consistent behavior.

The Python call re.reaction() might return the number of attempted and performed reactions, just like system.integrator.run() does.
In this way, the user could handle the case of no educts present, if they wish to.

kosovan · 2021-04-06T14:17:39Z

The Python call re.reaction() might return the number of attempted and performed reactions, just like system.integrator.run() does.
In this way, the user could handle the case of no educts present, if they wish to.

I like the solution proposed by @RudolfWeeber Indeed, the user might want to first create a reaction and then particles.

However, if we change the behaviour and return value of re.reaction(), then it would require more extensive changes, beyond a simple bugfix. I have some other suggestions on my list when it comes to modifying the reaction ensemble behaviour. Therefore, I propose to implement a simple bugfix now (throw an exception) and start creating a wishlist of deeper changes for the next coding day.

jonaslandsgesell · 2021-04-06T20:44:58Z

The problem #4197 is fixed via circumventing the asymetric proposal probability in #4207. The PR effectively shortens the code and makes it simpler. There is no need for a more complicated solution or further discussion of symptomatic fixes ;)

jngrad · 2021-05-18T17:00:27Z

PR superseded by #4207.

pm-blanco added 2 commits April 6, 2021 11:18

Solved formatting problems

e275731

KaiSzuttor reviewed Apr 6, 2021

View reviewed changes

pm-blanco added 5 commits April 6, 2021 13:05

Avoided the use of redundant bool variables and simplified the code

2c0615b

Exception error introduced that warns the user if the reaction cannot…

e4104d2

… occur because there are not enough reactants

formatting errors corrected

f808fbb

fixed bug preventing the user to do reactions when only having HA par…

c5b0c50

…ticles

formmating corrected

9c47f03

jngrad closed this May 18, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Introduced a bugfix to solve the infinite loop in constant pH algorithm #4200

Introduced a bugfix to solve the infinite loop in constant pH algorithm #4200

pm-blanco commented Apr 6, 2021 •

edited by jngrad

Loading

KaiSzuttor Apr 6, 2021 •

edited

Loading

pm-blanco Apr 6, 2021

kosovan Apr 6, 2021

KaiSzuttor Apr 6, 2021

jngrad Apr 6, 2021

pm-blanco Apr 6, 2021

kosovan Apr 6, 2021 •

edited

Loading

kosovan Apr 6, 2021

jonaslandsgesell Apr 6, 2021

jngrad Apr 7, 2021

RudolfWeeber commented Apr 6, 2021 •

edited by jngrad

Loading

RudolfWeeber commented Apr 6, 2021

kosovan commented Apr 6, 2021 •

edited

Loading

jonaslandsgesell commented Apr 6, 2021 •

edited

Loading

jngrad commented May 18, 2021

Introduced a bugfix to solve the infinite loop in constant pH algorithm #4200

Introduced a bugfix to solve the infinite loop in constant pH algorithm #4200

Conversation

pm-blanco commented Apr 6, 2021 • edited by jngrad Loading

KaiSzuttor Apr 6, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kosovan Apr 6, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

RudolfWeeber commented Apr 6, 2021 • edited by jngrad Loading

RudolfWeeber commented Apr 6, 2021

kosovan commented Apr 6, 2021 • edited Loading

jonaslandsgesell commented Apr 6, 2021 • edited Loading

jngrad commented May 18, 2021

pm-blanco commented Apr 6, 2021 •

edited by jngrad

Loading

KaiSzuttor Apr 6, 2021 •

edited

Loading

kosovan Apr 6, 2021 •

edited

Loading

RudolfWeeber commented Apr 6, 2021 •

edited by jngrad

Loading

kosovan commented Apr 6, 2021 •

edited

Loading

jonaslandsgesell commented Apr 6, 2021 •

edited

Loading