Remove Wang-Landau Reaction Ensemble algorithm (#4288)

The current Wang-Landau Reaction Ensemble (WLRE) implementation is tightly coupled to the other reaction methods in the C++ core. This accidental complexity slows down the development and improvement of the other reaction methods. Rewriting the WLRE method to improve modularity would be a significant project for the core team. This reaction method is not actively used by the community, and a recent call on the ESPResSo mailing list did not receive any response against the removal of the method:
https://lists.nongnu.org/archive/html/espressomd-users/2021-06/msg00001.html

Partial fix for #4251 (removing WLRE hooks)
Closes #4156

Description of changes:
- Remove the unused WLRE method. It is not actively maintained and not fully tested (86% code coverage).
- Fix undefined behavior in `ReactionAlgorithm` (member `particle_inside_exclusion_radius_touched` was uninitialized).
- Write additional unit tests.

doc/doxygen/bibliography.bib

@@ -589,17 +589,6 @@ @Article{tironi95a
  doi = {10.1063/1.469273},
 }
 
-@Article{troester05a,
- title = {{Wang-Landau sampling with self-adaptive range}},
- author = {Tr\"{o}ster, Andreas and Dellago, Christoph},
- journal = {Physical Review E},
- volume = {71},
- number = {6},
- pages = {066705},
- year = {2005},
- doi = {10.1103/PhysRevE.71.066705},
-}
-
 @Article{tyagi10a,
  author = {Tyagi, Sandeep and S\"{u}zen, Mehmet and Sega, Marcello and Barbosa, Marcia C. and Kantorovich, Sofia S. and Holm, Christian},
  title = {{An iterative, fast, linear-scaling method for computing induced charges on arbitrary dielectric boundaries}},
@@ -630,17 +619,6 @@ @Article{wang01a
  doi = {10.1063/1.1398588},
 }
 
-@Article{yan02b,
-author = {Yan, Qiliang and Faller, Roland and {de Pablo}, Juan J.},
-title = {{Density-of-states Monte Carlo method for simulation of fluids}},
-journal = {The Journal of Chemical Physics},
-volume = {116},
-number = {20},
-pages = {8745-8749},
-year = {2002},
-doi = {10.1063/1.1463055},
-}
-
 @Article{yeh99a,
  title = {{Ewald summation for systems with slab geometry}},
  author = {Yeh, In-Chul and Berkowitz, Max L.},

doc/sphinx/advanced_methods.rst

@@ -1770,27 +1770,6 @@ be converted to :math:`K_c` as
 where :math:`p^{\ominus}=1\,\mathrm{atm}` is the standard pressure.
 Consider using the python module pint for unit conversion.
 
-.. _Wang-Landau Reaction Ensemble:
-
-Wang-Landau Reaction Ensemble
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Combination of the Reaction Ensemble with the Wang-Landau algorithm
-:cite:`wang01a`. Allows for enhanced sampling of the reacting system
-and for the determination of the density of states with respect
-to the reaction coordinate or with respect to some other collective
-variable :cite:`landsgesell17a`. Here the 1/t Wang-Landau
-algorithm :cite:`belardinelli07a` is implemented since it
-does not suffer from systematic errors.
-
-Multiple reactions and multiple collective variables can be set.
-
-An example script can be found here:
-
-* `Wang-Landau reaction ensemble <https://github.com/espressomd/espresso/blob/python/samples/wang_landau_reaction_ensemble.py>`__
-
-For a description of the available methods, see :class:`espressomd.reaction_ensemble.ReactionEnsemble`.
-
 .. _Grand canonical ensemble simulation using the Reaction Ensemble:
 
 Grand canonical ensemble simulation

doc/sphinx/zrefs.bib

@@ -101,18 +101,6 @@ @article{beenakker86a
  doi = {10.1063/1.451199}
 }
 
-@article{belardinelli07a,
- title={Fast algorithm to calculate density of states},
- author={Belardinelli, R. E. and Pereyra, V. D.},
- journal={Phys. Rev. E},
- volume={75},
- number={4},
- pages={046701},
- year={2007},
- doi={10.1103/PhysRevE.75.046701},
- publisher={APS}
-}
-
 @ARTICLE{bereau15,
  author = {Tristan Bereau},
  title = {Multi-timestep integrator for the modified Andersen barostat},
@@ -513,17 +501,6 @@ @Book{kruger12a
 doi={10.1007/978-3-8348-2376-2},
 }
 
-@Article{landsgesell17a,
-title = {{W}ang-{L}andau Reaction Ensemble Method: {S}imulation of Weak Polyelectrolytes and General Acid-Base Reactions},
-author = {Landsgesell, Jonas and Holm, Christian and Smiatek, Jens},
-journal = {J. Chem. Theory Comput.},
-volume = {13},
-number = {2},
-pages = {852--862},
-year = {2017},
-doi = {10.1021/acs.jctc.6b00791}
-}
-
 @Article{landsgesell17b,
 title = {Simulation of weak polyelectrolytes: {A} comparison between the constant p{H} and the reaction ensemble method},
 author = {Landsgesell, Jonas and Holm, Christian and Smiatek, Jens},
@@ -929,17 +906,6 @@ @article{wagner02
  doi = {10.1023/A:1014595628808}
 }
 
-@article{wang01a,
- title={Efficient, multiple-range random walk algorithm to calculate the density of states},
- author={Wang, Fugao and Landau, David P},
- journal={Phys. Rev. Lett.},
- volume={86},
- number={10},
- pages={2050},
- year={2001},
- publisher={APS}
-}
-
 @ARTICLE{wolff04a,
  author = {Ulli Wolff},
  title = {{M}onte {C}arlo errors with less errors},

src/core/reaction_methods/CMakeLists.txt

@@ -22,7 +22,6 @@ target_sources(
  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp
  ${CMAKE_CURRENT_SOURCE_DIR}/ReactionEnsemble.cpp
  ${CMAKE_CURRENT_SOURCE_DIR}/ConstantpHEnsemble.cpp
- ${CMAKE_CURRENT_SOURCE_DIR}/WangLandauReactionEnsemble.cpp
  ${CMAKE_CURRENT_SOURCE_DIR}/WidomInsertion.cpp
  ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp)
 

src/core/reaction_methods/ConstantpHEnsemble.cpp

@@ -35,7 +35,7 @@ namespace ReactionMethods {
  */
 double ConstantpHEnsemble::calculate_acceptance_probability(
  SingleReaction const &current_reaction, double E_pot_old, double E_pot_new,
- std::map<int, int> const &old_particle_numbers, int, int, bool) const {
+ std::map<int, int> const &old_particle_numbers) const {
  auto const beta = 1.0 / temperature;
  auto const pKa = -current_reaction.nu_bar * log10(current_reaction.gamma);
  auto const ln_bf = (E_pot_new - E_pot_old) - current_reaction.nu_bar / beta *

src/core/reaction_methods/ConstantpHEnsemble.hpp

@@ -46,8 +46,8 @@ class ConstantpHEnsemble : public ReactionAlgorithm {
 protected:
  double calculate_acceptance_probability(
  SingleReaction const &current_reaction, double E_pot_old,
- double E_pot_new, std::map<int, int> const &old_particle_numbers, int,
- int, bool) const override;
+ double E_pot_new,
+ std::map<int, int> const &old_particle_numbers) const override;
 };
 
 } // namespace ReactionMethods

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -293,12 +293,9 @@ void ReactionAlgorithm::generic_oneway_reaction(
 
  current_reaction.tried_moves += 1;
  particle_inside_exclusion_radius_touched = false;
- int old_state_index = -1; // for Wang-Landau algorithm
- on_reaction_entry(old_state_index);
  if (!all_reactant_particles_exist(current_reaction)) {
  // makes sure, no incomplete reaction is performed -> only need to consider
  // rollback of complete reactions
- on_reaction_rejection_directly_after_entry(old_state_index);
  return;
  }
 
@@ -336,23 +333,15 @@ void ReactionAlgorithm::generic_oneway_reaction(
  ? std::numeric_limits<double>::max()
  : calculate_current_potential_energy_of_system();
 
- int new_state_index = -1; // save new_state_index for Wang-Landau algorithm
- int accepted_state = -1; // for Wang-Landau algorithm
- on_attempted_reaction(new_state_index);
-
- constexpr bool only_make_configuration_changing_move = false;
  auto const bf = calculate_acceptance_probability(
- current_reaction, E_pot_old, E_pot_new, old_particle_numbers,
- old_state_index, new_state_index, only_make_configuration_changing_move);
+ current_reaction, E_pot_old, E_pot_new, old_particle_numbers);
 
  std::vector<double> exponential = {
  exp(-1.0 / temperature * (E_pot_new - E_pot_old))};
  current_reaction.accumulator_exponentials(exponential);
 
  if (m_uniform_real_distribution(m_generator) < bf) {
  // accept
- accepted_state = new_state_index;
-
  // delete hidden reactant_particles (remark: don't delete changed particles)
  // extract ids of to be deleted particles
  auto len_hidden_particles_properties =
@@ -374,7 +363,6 @@ void ReactionAlgorithm::generic_oneway_reaction(
  current_reaction.accepted_moves += 1;
  } else {
  // reject
- accepted_state = old_state_index;
  // reverse reaction
  // 1) delete created product particles
  for (int p_ids_created_particle : p_ids_created_particles) {
@@ -386,7 +374,6 @@ void ReactionAlgorithm::generic_oneway_reaction(
  restore_properties(changed_particles_properties,
  number_of_saved_properties);
  }
- on_end_reaction(accepted_state);
 }
 
 /**
@@ -566,22 +553,14 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
  * Performs a global MC move for a particle of the provided type.
  */
 bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
- int type, int particle_number_of_type_to_be_changed, bool use_wang_landau) {
+ int type, int particle_number_of_type_to_be_changed) {
  m_tried_configurational_MC_moves += 1;
  particle_inside_exclusion_radius_touched = false;
 
- int old_state_index = -1;
- if (use_wang_landau) {
- on_reaction_entry(old_state_index);
- }
-
  int 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) {
  // reject
- if (use_wang_landau) {
- on_mc_rejection_directly_after_entry(old_state_index);
- }
  return false;
  }
 
@@ -596,20 +575,8 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
 
  double beta = 1.0 / temperature;
 
- int new_state_index = -1;
- double bf = 1.0;
- std::map<int, int> dummy_old_particle_numbers;
- SingleReaction temp_unimportant_arbitrary_reaction;
-
- if (use_wang_landau) {
- new_state_index = on_mc_use_WL_get_new_state();
- bf = calculate_acceptance_probability(
- temp_unimportant_arbitrary_reaction, E_pot_old, E_pot_new,
- dummy_old_particle_numbers, old_state_index, new_state_index, true);
- } else {
- // Metropolis algorithm since proposal density is symmetric
- bf = std::min(1.0, bf * exp(-beta * (E_pot_new - E_pot_old)));
- }
+ // Metropolis algorithm since proposal density is symmetric
+ auto const bf = std::min(1.0, exp(-beta * (E_pot_new - E_pot_old)));
 
  // // correct for enhanced proposal of small radii by using the
  // // Metropolis-Hastings algorithm for asymmetric proposal densities
@@ -618,24 +585,16 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
  // std::pow(particle_positions[0][1]-cyl_y,2));
  // double new_radius =
  // std::sqrt(std::pow(new_pos[0]-cyl_x,2)+std::pow(new_pos[1]-cyl_y,2));
- // bf = std::min(1.0,
- // bf*exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
+ // auto const bf = std::min(1.0,
+ // exp(-beta*(E_pot_new-E_pot_old))*new_radius/old_radius);
 
  // Metropolis-Hastings algorithm for asymmetric proposal density
  if (m_uniform_real_distribution(m_generator) < bf) {
  // accept
  m_accepted_configurational_MC_moves += 1;
- if (use_wang_landau) {
- on_mc_accept(new_state_index);
- }
  return true;
  }
- // reject
- // modify wang_landau histogram and potential
- if (use_wang_landau) {
- on_mc_reject(old_state_index);
- }
- // restore original particle positions
+ // reject: restore original particle positions
  for (auto const &item : original_positions)
  place_particle(std::get<0>(item), std::get<1>(item));
  return false;