Reaction method refactoring wish list

[RudolfWeeber]

actoring wish list (Disclaimer: I'm not terribly familiar with that part of the code)

This is a WIP

ReactionAlgorithm::make_reaction_attempt()

• Are the non const vector arguments actually used only as outputs? If so, they should be returned rather than being arguments, to make the dependencies clear.

ReactionAlgorithm::generic_oneway_reaction()

• This should take a reference to a SingleReaction to decouple the reaction from the storage of the reactions.
• Looks like some Wang Landau stuff is mixed in
• The "exponentials" accumulating accumulates something else than the Boltzmann factor used in the acceptance probability. Is that intended?
• The particle (property) hiding/deleting/restoring is not clear to me. It needs to be explained somewhere, or if it is, the explanation needs to be referenced. Why is hiding preferred to outright deletion?
• Rather than calling the on_.*() methods, the function should communicate the outcome via return value (a typed enum). This would allow to avoid mixing in code for a special case (Wang Landau) in the implementation of the general case. I guess, outcomes are: early rejection, rejection, acceptance.

ReactionAlgorithm::delete_particle() and ReactionAlgorithm::create_particle(int)

• The functionality of managing the max seen particle counter and obtaining an unused particle id for particle creation should be pulled out of the reaction ensemble. Other algorithms need that too (and have their own private implementations...)
• The expensive partCfg()and distto() stuff can be skipped if the exclusion radius is 0 (several instances throughout the code)

general

• The task of finding/checking for a close-by particle (currently done with distto()) can be accelerated considerably using the domain decomposition cell system. O(1) vs O(N)

ReactionAlgorithm::do_global_mc_move_for_particles_of_type()

• The task of finding ids and positions of particles to move should be outsourced to a function returning the two vectors. The logic can probably be simplified. Starting L563. Then, it could be tested independently.
• rather than immediately moving the particles, generate all new positions/velocities and store in an array. When setting positions, the particle data caches are invalidated. Or even better: Collect the particle masses when collecting the old positions and ids. then, no additional get_particle_data() is needed at all.
• there is an other instance of the expensive distto() logic for the exclusion radius
• Creating a dummy reaction and then using its acceptance probability method is an indication that the acceptance probability should be a free function.

ConstantpHEnsemble::do_reaction()

• The reaction choice method should be moved to a free function returning a reference to a SingleReaction. This would make the choosing logic testable.
• The logic for reaction choice can probably also be made more readable using STL algorithms (find or the like)

ConstantpHEnsemble::calculate_acceptance_probability

• What are the dummy variables for?

Wang Landau on_.*() hooks

• methods should be removed. Instead, return values from ReactionAlgorithm::generic_on_way_reaction() should be used to act on outcomes of a reaction.

WangLandauReactionEnsemble::get_flattened_index_wang_landau()

• This should make use of the functions in utils/include/utils/index.hpp, if possible.
• All the loops can/should probably be replaced by std::transform or the like

Histogram

• Make use of the histogram implementation in src/utils/iinclude/utils/Histogram.hpp
  This might remove the need for any index flattening in the first place.

WidomInsertion::measure_excess_chemical_potential()

• This should take a reference to a SingleReaction's argument to decouple from the SingleReaction storage.

Other points

From #4251 (comment)

• Optimizing the exclusion radius logic, make the search for particles close to the test position O(1). This will likely also solve the performance issue with tutorial 12
• Disentangling the Wang-Landau stuff from the ReactionAlgorithm base class, e.g., by communicating outcomes from reaction attempt and generic one way reactions via return value. Also, the different requirements for the acceptance probability need to be sorted out.
• Making use of Utils::Histogram and/or Utils::(un)_ravl_index() in Wang Landau. This will reduce the code and increase clarity considerably.

From #4251 (comment)

• reaction_ensemble.ConstantpHEnsemble takes as an argument 'temperature' but actually kT is expected as input. We think that for consistency this variable should be renamed as 'kT'.
• reaction_ensemble.ConstantpHEnsemble has as an argument reactant_coefficients and product_coefficients. However, for the constant pH method the only possible inputs for these variables are reactant_coefficients=[1] and product_coefficients=[1, 1]. Thus, we think that these variables should be internally hard coded and not be arguments of this function.

[jonaslandsgesell]

These are only some short answers. Here it is really late night...

ReactionAlgorithm::make_reaction_attempt()

• There are loops over the difference of number of reactants and products. Can this be written using the nu_bar property of the SingleReaction?

The value nubar to which you refer is an accumulated value: nu_bar = \sum_i \nu_i, the stochiometric coefficients nu_i are the changes in particle numbers and contain the necessary information for performing the reaction: \nu_1=+1,\nu_2=-1,\nu_3=+1 has \nu_bar=1, but \nu_1=+1,\nu_2=0,\nu_3=0 also has \nu_bar=1. The value of nu_bar alone is not enough .(see https://www.tandfonline.com/doi/abs/10.1080/08927020801986564)

• Looks like some Wang Landau stuff is mixed in

The ReactionEnsemble, cpHEnsemble, WangLandauReactionEnsemble currently inherit from ReactionAlgorithm, bceause each of them "is a ReactionAlgorithm". The idea to inject modifing behavior (e.g. in the WangLandauReactionEnsemble) makes use of the template method pattern: Hooks are employed which are e.g. empty in the superclass or are otherwise overridden.

• Line 324: if (particle_inside_exclusion_raidus_touched): I don't see where this can be set further up in the function (except for the initializaiton to false at the top)

In the superclass there is additional implementation details

espresso/src/core/reaction_methods/ReactionAlgorithm.cpp

Line 526 in 2d88a51

particle_inside_exclusion_radius_touched = true;

,

espresso/src/core/reaction_methods/ReactionAlgorithm.cpp

Line 592 in 2d88a51

particle_inside_exclusion_radius_touched = true;

Note that particle_inside_exclusion_raidus_touched turned out to be very important to ensure detailed balance (and correct sampling of the partition functions).

• The "exponentials" accumulating accumulates something else than the Boltzmann factor used in the acceptance probablity. Is that intended?

This nomenclature is left over from Tobias. Feel free to rename to something like acceptanceProbability.

• The particle (property) hiding/deleting/restoring is not clear to me. It needs to be explained somewhere, or if it is, the explanation needs to be referenced. Why is hiding preferred to outright deleiton?

The hiding is done in order to not have to take care of the "bonding topology" in polymeric systems when only changing some particle properties like the type, the charge,...

ReactionAlgorithm::delete_particle() and ReactionAlgorithm::create_particle(int

• the functionallity of managing the max seen particle counter and obtaining an unused particle id for particle creation should be pulled out of the reacion ensemble. Other algorithms need that too (and have their own privat implmeentations...)
• The expensive part_cfg()and dist_to stuff can be skipped if the exclusion radius is 0 (don't know, how common that is) (several instances throughout the code)

the exclusion radius can only be 0 if:

1. a non interacting system is simulated (can be solved analytically in many cases) or
2. an interacting system is simulated with displacment MC moves only (this is the main target)

One will always need the exclusion radius if MD is employed in an interacting system. Particles are placed close to each other in MC and resulting in force induced instabilities in MD with LJ-potentials.

general

• Considering the expensive check, why is the exclusoin radius needed? Wouldn't the move be rejected anyway?
  (see above)

• rather than immediately moving the particles, generate all new positions/velocities and store in an array. When setting positions, the particle data caches are invalidated. Or even better: Collect the particle masses when collecting the old positions and ids. then, no additional get_particle_data() is needed at all.

Tobi identified that it was easier to change charge & type instead of reconstructing potentially arbitrary properties of particles (like bondings, masses, etc.) This is the "hiding strategy".

ConstantpHEnsemble::do_reaction()

• The reaction choice method should be moved to a free function returning a reference to a SingeReaction. This would make the choosing logic testable.
• The logic for reaction choice can probably alo be made more readable using stl algorithms (find or the like)
  the logic in

This function can be completely deleted: https://github.com/espressomd/espresso/pull/4207/files#diff-31488bc3f9a69b6dd9b0349f71bdfba05fb8a1d20d4c640d837fd15bd43eb7feL41
Please merge this PR.

ConstantpHEnsemble::calculate_acceptance_probability

• What are the dummy variables for?

Currently, the variables int dummy_old_state_index, int dummy_new_state_index, bool dummy_only_make_configuration_changing_move exist because an identical signature is needed for overriding.

[RudolfWeeber]

As far as I can tell from reading through the code, there are three big items + a few minor ones.
The main points are:

1. Optimizing the exclusion radius logic, make the search for particles close to the test position O(1). This will liekly also solve the performance issue with tutorial 12
2. Disentangling the Wang-Landau stuff from the ReactionAglorithm base class, e.g., by communicating outcomes from reaction atempt and generic one way reactions via return value. Also, the different requirements for the acceptance probability need to be sorted out.
3. Making use of Utils::Histogram and/or Utils::(un)_ravl_index in Wang Landau. This will reduce the code and increase clarity considerably.

[RudolfWeeber]

I removed the items that were sorted out in Zoom the discussions

[pm-blanco]

@kosovan and I want to add a couple of things to the wish list:

• reaction_ensemble.ConstantpHEnsemble takes as an argument 'temperature' but actually kT is expected as input. We think that for consistency this variable should be renamed as 'kT'.
• reaction_ensemble.ConstantpHEnsemble has as an argument reactant_coefficients and product_coefficients. However, for the constant pH method the only possible inputs for these variables are reactant_coefficients=[1] and product_coefficients=[1, 1]. Thus, we think that these variables should be internally hard coded and not be arguments of this function.

[jngrad]

• reaction_ensemble.ConstantpHEnsemble takes as an argument 'temperature' but actually kT is expected as input. We think that for consistency this variable should be renamed as 'kT'.

I have no strong opinion on this proposed change. In the ESPResSo C++ code, variables called temperature in both the MD and Reaction Methods code is in units of kT. This is documented in the user guide for the MD code:

espresso/doc/sphinx/introduction.rst

Lines 391 to 400 in 0327d9a

	A special note is due regarding the temperature, which is coupled to the
	energy scale by Boltzmann's constant. However, since |es| does not enforce a
	particular unit system, we also don't know the numerical value of the
	Boltzmann constant in the current unit system. Therefore, when
	specifying the temperature of a thermostat, you actually do not define
	the temperature, but the value of the thermal energy :math:`k_B T` in
	the current unit system. For example, if you measure energy in units of
	:math:`\mathrm{kJ/mol}` and your real temperature should be 300K, then
	you need to set the thermostat's effective temperature to
	:math:`k_B 300\, K \mathrm{mol / kJ} = 2.494`.

• reaction_ensemble.ConstantpHEnsemble has as an argument reactant_coefficients and product_coefficients. However, for the constant pH method the only possible inputs for these variables are reactant_coefficients=[1] and product_coefficients=[1, 1]. Thus, we think that these variables should be internally hard coded and not be arguments of this function.

Maybe @jonaslandsgesell can tell us about this design decision. I guess it is meant to have the same python interface between all RE methods, but could also be due to planned extensions of the constant pH method.

[kosovan]

To my understanding, this parameter is called temperature because by the time when it was implemented, parameters of the thermostat were called the same. Later, the thermostat parameters were renamed to kT but the RE parameter was not renamed.

Therefore, we should rename it also in the reaction ensemble to ensure consistency. Then, the documentation could also become much simpler. Instead of repeating the generic discussion of units, that is now in the user guide, it would suffice to state that we provide the numerical value of kT in the reduced units and link to the user guide.

I recall a discussion that there should be just one global value of temperature or kT in Espresso, shared by all algorithms, because physically it does not make sense to use different temperatures within one system. I am not sure if there is a plan to implement this idea. Perhaps @RudolfWeeber can comment on this.

Requiring the user to provide parameters that have only one possible value is not a good practice. As far as I know, no other values of stoichiometric coefficients are possible for the constant-pH ensemble, and there is no such extension in the outlook.

[pm-blanco]

I can continue with the pending tasks on this issue tomorrow!

[pm-blanco]

@RudolfWeeber @jngrad I have been trying to improve the searching of a nearby particle for the exclusion radius using the cell system, but I was not able to find the way of getting the particles that are in the same cell than a given particle (or something similar). Is there somewhere in the code where I could get an example on how to do it? Thank you in advance!

[RudolfWeeber]

Can this be closed?

[pm-blanco]

There is still the issue of using the cell system for finding the neighboring particles. But maybe we can address that specifically in a separate issue and close this one. Otherwise it can be closed in my opinion.

[jngrad]

@pm-blanco could you please open a separate ticket outlining the neighbor search? Then we can close this ticket.