Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

introduce symmetric proposal probability in cpH method #4207

Conversation

jonaslandsgesell
Copy link
Member

@jonaslandsgesell jonaslandsgesell commented Apr 6, 2021

Fixes #4197

Description of changes:

  • Cleanup constant pH code via removal of asymmetric proposal probabilities, now implements constant pH method with symmetric proposal probability in Landsgesell et al. 2017 (doi:10.1140/epjst/e2016-60324-3), equation 15
  • Fix infinite loop when the H+ and A- types are accidentally switched when defining the acid-base reaction

@jonaslandsgesell jonaslandsgesell changed the title introduce symmetric proposal probability introduce symmetric proposal probability in cpH method Apr 6, 2021
@jonaslandsgesell jonaslandsgesell force-pushed the simplify_cph_ensemble_with_symmetric_proposal_probability branch from 2604465 to 29a614a Compare April 6, 2021 20:48
@jngrad
Copy link
Member

jngrad commented Apr 6, 2021

@jngrad what is the difference between the tests https://github.com/espressomd/espresso/blob/python/testsuite/python/constant_pH.py and https://github.com/espressomd/espresso/blob/python/testsuite/python/constant_pH_stats.py ? Can one be removed?

The former checks the interface using a chemical reaction with cherry-picked parameters to run in a few seconds in coverage builds. The latter is a statistical test that checks convergence. They are otherwise almost identical, and constant_pH_stats.py will probably be removed in the future.

@jonaslandsgesell
Copy link
Member Author

@jngrad I get a message about some cpH core tests failing which you introduced, could you please have a look? I guess the values are hardcored and do not match anymore.

@jngrad
Copy link
Member

jngrad commented Apr 6, 2021

@jngrad I get a message about some cpH core tests failing which you introduced, could you please have a look? I guess the values are hardcored and do not match anymore.

You will need to adapt the acceptance rate in the unit test with the result of auto const g = calculate_factorial_expression_cpH(...).

But before you spend more time on this, the core team should first decide whether we really want to introduce this silent breaking change to the Constant pH algorithm.

@jonaslandsgesell
Copy link
Member Author

The statistical properties of the algorithm before and after the change are shown to be the same in the referred paper. I would be in favour of delivering this change.

@kosovan
Copy link
Contributor

kosovan commented Apr 7, 2021

@jonaslandsgesell can you please explain how the asymmetric probabilities fix the bug when the reaction is not possible in either direction, as described in #4197 As far as I understand, moving from symmetric to asymmetric proposal probabilities still ensures that the sum of probabilities is equal to 1.0.

@kosovan
Copy link
Contributor

kosovan commented Apr 7, 2021

The statistical properties of the algorithm before and after the change are shown to be the same in the referred paper. I would be in favour of delivering this change.

I am convinced that it works on paper but I think that the implementation should be tested on various non-trivial cases before it is merged. I am afraid that the standard test cases might not be sufficient for that. We should also verify its impact on performance.

@RudolfWeeber proposed to have a discussion on the future of reaction ensemble. This is a big change that affects the core of the reaction ensemble, so I suggest to wait with this merge until this discussion takes place.

@jonaslandsgesell
Copy link
Member Author

@kosovan all_reactant_particles_exist() is checked in generic_oneway_reaction. If it returns false no reaction is attempted because a reaction step removing the non-existing reactant particles is not possible. You can also see this via constructing any situation you like and attaching a debugger. The problem in the previous code was the way the asymmetric proposal probability was implemented.

I would not defer this merge because it
A) removes some bogous code which made the proposal probabilities asymmetric.
B) If you want to make other changes please do not introduce a coupling with regards to this PR. The changes you want to make should be regarded independently.

What non-trivial test you have in mind?

@kosovan
Copy link
Contributor

kosovan commented Apr 7, 2021

A) removes some bogous code which made the proposal probabilities asymmetric.

I am not objecting against these changes, I just want to be sure that they are thoroughly tested.

B) If you want to make other changes please do not introduce a coupling with regards to this PR. The changes you want to make should be regarded independently.

I do not understand. If I just mention another issue, does github automatically treat them as dependent?

What non-trivial test you have in mind?

I need some time to think of that. I was not aware that these changes were planned, so I have not thought of it before.

Some points that come to my mind immediately: The standard test cases are not exactly reproducible, they have to rely on some arbitrarily chosen tolerance. This change affects the acceptance probability, and if there is a bug, it might not be immediately visible with a tolerance that is set so that the test completes quickly. Therefore, I would try to generate high-quality data using the old and the new implementation, and check that they produce identical results. Definitely, I would focus on various corner cases, with alpha approx. 0 and alpha approx. 1, and I would like to test if/how the performance has changed with the new code.

@KaiSzuttor
Copy link
Member

Some points that come to my mind immediately: The standard test cases are not exactly reproducible, they have to rely on some arbitrarily chosen tolerance. This change affects the acceptance probability, and if there is a bug, it might not be immediately visible with a tolerance that is set so that the test completes quickly. Therefore, I would try to generate high-quality data using the old and the new implementation, and check that they produce identical results. Definitely, I would focus on various corner cases, with alpha approx. 0 and alpha approx. 1, and I would like to test if/how the performance has changed with the new code.

If you don't trust the current test setup it's impossible to safely change any part of the reaction ensemble...

@kosovan
Copy link
Contributor

kosovan commented Apr 7, 2021

If you don't trust the current test setup it's impossible to safely change any part of the reaction ensemble...

The current cpH test completes on my laptop in 2.807s after it runs 1000 samples for one specific set of conditions. As far as I remember, it has been chosen as a compromise between accuracy and speed, with a preference for speed. This is sufficient in many cases, if you want to check that the code does not crash due to a programming error but it is not sufficient to check that the physics has not changed.

I think that one needs different levels of cautiousness, depending on the possible impact a particular change has. The change in this pull request has similar impact on the cpH calculations as if you would replace the integrator scheme or a random number generator. It is definitely worth more testing before it comes to production.

For example, the test case only tests an ideal system (no interactions), so it cannot catch even such a simple mistake as a wrong sign before Delta U in the acceptance probability.

@jonaslandsgesell
Copy link
Member Author

For example, the test case only tests an ideal system (no interactions), so it cannot catch even such a simple mistake as a wrong sign before Delta U in the acceptance probability.

I know you just made an example, but the changes introduced do not affect interactions because the factor involving interactions is untouched...

@kosovan
Copy link
Contributor

kosovan commented Apr 7, 2021

I know you just made an example, but the changes introduced do not affect interactions because the factor involving interactions is untouched...

It should not affect interactions. But maybe it affects something else that you have not considered. Just give us some time to test it on a non-trivial problem.

@pm-blanco
Copy link
Contributor

pm-blanco commented Apr 7, 2021

Just to be sure, I can perform simulations for some interacting systems in some non-trivial conditions and check that the results I get with your new implementation are the same that the ones I was getting with the actual code. Also could be a good chance to check the performance :)

@jonaslandsgesell jonaslandsgesell force-pushed the simplify_cph_ensemble_with_symmetric_proposal_probability branch from 5cda7b4 to 3d142c3 Compare April 7, 2021 20:53
@pm-blanco
Copy link
Contributor

pm-blanco commented Apr 10, 2021

I did some calculations of the titration curve of a weak polyacid (with WCA and electrostatic interactions). Attached you can find the results with the new implementation (label "new") and the current development code (label "devel"). In all the pH range both algorithms yield to the same average degree of ionization (\alpha) within the estimated error. (Note: the results are rounded to the third decimal).
I also tested that this new implementation does indeed solve the infinite loop bug that I mentioned in #4197, by returning a no-op. I also noticed that the new implementation seems to be slightly faster than the current develop version.
In my opinion, the branch can be safely merged.
alpha-comparison
alpha-abs

@jonaslandsgesell
Copy link
Member Author

jonaslandsgesell commented Apr 10, 2021

PS: The cpp cpH unittest also checks the impact of the input energy.

@jngrad jngrad added this to the Espresso 4.2 milestone May 18, 2021
@jngrad jngrad added the BugFix label May 18, 2021
@jngrad jngrad added the Core label May 18, 2021
@jngrad jngrad added the automerge Merge with kodiak label May 18, 2021
@kodiakhq kodiakhq bot merged commit d7d9b49 into espressomd:python May 18, 2021
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Constant pH method can enter an infinite loop if some particles are not created initially
6 participants