Widom excess chemical potential calculation and RxMC performance improvement #4374

pm-blanco · 2021-10-20T17:47:53Z

Description of changes:

WidomInsertion.measure_excess_chemical_potential() was replaced by WidomInsertion.calculate_excess_chemical_potential() which returns the instantaneous value of the excess chemical potential instead of the accumulated mean and std. The mean value and std of the excess chemical potential must be now be calculated by the user elsewhere. An example of how to do it is given in samples/widom_insertion.py.
ReactionAlgorithm::do_reaction() now only computes the potential energy of the old state at the begininng of the reaction loop. This value is updated after each successful reaction trial move.

RudolfWeeber · 2021-10-21T07:21:02Z

Contentwise this looks good. However, we use semantic versioning. I.e., as long as the major version (4 in 4.1.2) is the same, we promise not to change the outcome of a simulation script silently.
To make your change non-silent, please rename the python function, i.g., calc_instantaneous_chemical_potential()

samples/widom_insertion.py

jonaslandsgesell · 2021-10-21T08:30:34Z

samples/widom_insertion.py

+        np.mean(gamma_samples[block * size_block:(block + 1) * size_block]))
+
+mu_ex_mean = -1 * np.log(np.mean(gamma_block)) * temperature
+mu_ex_std = np.sqrt(np.var(gamma_block) / np.mean(gamma_block)) * temperature


Please add a note that this is due to linearized error propagation. In principle you can do better via not using linearization: mu_ex_Delta =abs( -1 * np.log(np.mean(gamma_block)+np.std(gamma_block, dof =1)) * temperature -(-1 * np.log(np.mean(gamma_block)-np.std(gamma_block, dof =1)) * temperature)

jonaslandsgesell · 2021-10-21T08:32:16Z

src/python/espressomd/reaction_ensemble.pyx

-        the excess chemical potential. The error estimate assumes that
-        your samples are uncorrelated.
-
+        Returns the instantaneous excess chemical potential


There is no instantaneous chemical potential (which is always an ensemble average, please find a better name), it would also miss the application of -kT log().

If you change the logic of how this function works, please describe the needed computations in the doc string

testsuite/python/widom_insertion.py

jonaslandsgesell · 2021-10-21T08:38:43Z

There is no such thing as an instanteous chemical potential. The chemical potential is per definition an ensemble average. Additionally, the value returned cannot be energy: The application of -kT log is missing (surely because log and the ensemble average do not commute). You could name it compute_exp_of_potential_energy_change_for_excess_chemical_potential_computatiom() which would at be a telling name.

src/core/reaction_methods/WidomInsertion.cpp

jonaslandsgesell

The current naming is misleading, please improve

src/core/reaction_methods/ReactionAlgorithm.cpp

jonaslandsgesell · 2021-10-21T15:47:17Z

samples/widom_insertion.py

@@ -123,7 +123,7 @@
 n_iterations = 5000
 for i in range(n_iterations):

-    mu_ex = widom.calc_instantaneous_chemical_potential(insertion_reaction_id)
+    E_ins = widom.calc_particle_insertion_potential_energy(insertion_reaction_id)


consider E_pot_insertion. The inserted particles also carry kinetic energy which is not reported in this observable

jonaslandsgesell · 2021-10-21T15:48:19Z

src/core/reaction_methods/WidomInsertion.cpp


-  const double mu = E_pot_new - E_pot_old;
+  const double E_ins = E_pot_new - E_pot_old;


I think we should avoid abbreviation whenever possible https://medium.com/ltunes/what-is-clean-code-naming-conventions-part-1-426d383eb85d

src/core/reaction_methods/WidomInsertion.hpp

pm-blanco · 2021-10-21T16:21:00Z

@jonaslandsgesell I belive I have resolved all your comments, please take a look and let me know if you find anything else

jonaslandsgesell · 2021-10-21T17:53:02Z

Thank you for taking care.

samples/widom_insertion.py

jonaslandsgesell · 2021-10-21T20:04:53Z

This PR would be a good check for your benchmark PR. What is the performance before/after your suggested change, when running e.g. 10 reactions in a row for 1000 times?

pm-blanco · 2021-10-22T09:22:46Z

@jonaslandsgesell I did a performance test with the benchmark script of my PR. Performing 1000 MC reactions steps prior my implementation takes ~ 1.5 s on my laptop. Performing the same number of MC steps with the implementation proposed in this PR takes ~ 0.9 s on my laptop. For a sufficiently large number of reaction steps you gain roughly a factor of 2 in computational time (only in the reaction side, of course).

jonaslandsgesell · 2021-10-28T10:44:14Z

Could you please add the unittests which checks if the energy is altered after it is passed in? Then this PR is fine from my side

jngrad · 2021-10-28T15:28:01Z

Once the comments by Jonas are resolved, it would be nice to rebase the PR. I find it challenging to follow the changes, because the 19 commits alternate between the Widom refactoring and the ReactionAlgorithm::do_reaction for loop optimization (plus a few random merge and formatting commits). Also most commits don't pass CI, so if I select a commit revision and the code doesn't build or doesn't pass the tests on my workstation, I don't know if the issue comes from my build environment or from the PR.

src/python/espressomd/reaction_ensemble.pyx

samples/widom_insertion.py

pm-blanco · 2021-10-29T07:16:30Z

@jonaslandsgesell @jngrad I will try to attend your comments ASAP. Honestly, I have not yet done the unit test because I am unsure on how to implement such test straightforwardly, but I will try to propose something.

jonaslandsgesell · 2021-10-29T11:45:49Z

@jonaslandsgesell @jngrad I will try to attend your comments ASAP. Honestly, I have not yet done the unit test because I am unsure on how to implement such test straightforwardly, but I will try to propose something.

@jngrad , could you please assist in the unittest? It might require refactoring to get something which is easily testable.

pm-blanco · 2021-10-31T18:29:21Z

@jonaslandsgesell it took me a while to figure it out how to do it properly, but I managed to implement one unit test in src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp. Basicaly I define a simple identity change reaction D <-> E with a reaction constant Gamma = infinity (i.e. completly shifted to products). Then I create one particle D in the system and I call ReactionAlgorithm::generic_oneway_reaction providing a fake energy of the system equal to 1. Since it is an ideal system with no interactions and the reaction always will be accepted, after this call the energy should be updated to 0. I checked and the test fails if in ReactionAlgorithm::generic_oneway_reaction the energy is provided by value and not by reference and it passes otherwise.
@jngrad I am sorry, but I did not rebase the git history this time. I prefer to practice and learn how to do it properly than mess this PR. Next time I will try to do a more clean PR. I applied all your other revisions.

pm-blanco · 2021-10-31T19:50:53Z

@jngrad it is now failing the CI due to some weird error due to an "unexpected indentation in the docstring" in \src\python\espressomd\reaction_ensemble.pyx, any idea why?

src/python/espressomd/reaction_ensemble.pyx

src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp

src/python/espressomd/reaction_ensemble.pyx

jngrad · 2021-11-01T14:53:58Z

For the rebase, it's as simple as

# backup current state of the branch
git branch Widom+RxMC_fix_backup
# squash all changes
git reset 5bab6bd51e7

Then do git add -p to select which lines belong to the performance improvement change for the reaction algorithm, commit it, and add the rest to the Widom commit. You can recover the original history at any time with the backed up branch.

Cache the potential energy value. This saves the cost of an energy calculation per loop when running multiple reactions.

Widom now returns the particle insertion potential energy and calculates the excess chemical potential afterwards with block analysis.

pm-blanco · 2021-11-02T14:23:58Z

@jngrad the PR is now blocked because it does not detect that I did the changes you requested, I belive

pm-blanco changed the title ~~Widom+rx mc fix~~ Widom+RxMC fix Oct 20, 2021