Refactored MC displacement move of particles

[pm-blanco]

Description of changes:

• Substituted bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type, int n_part) for a new method to perform MC displacement moves of particles void ReactionAlgorithm::do_particle_displacement_MC_move( int mc_steps, std::vector<int> particle_types_to_move)
• The new method takes the number of trial MC steps mc_steps as input instead of the number of particles to move, avoiding confusion on the actual number of MC steps performed. If no value is given performs one MC steps by default.
• Only one particle is moved in each MC step.
• The new method also takes an optional argument particle_types_to_move, which is the list of particle types from which particles will be selected. If no value is given all particle can be moved by default.
• The method ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) has been removed since it was only used by ReactionAlgorithm::do_particle_displacement_MC_move
• I renamed the variables m_tried_configurational_MC_moves and m_accepted_configurational_MC_moves to N_trial_particle_displacement_MC_moves and N_accepted_particle_displacement_MC_moves for clarity.
• I adapted all tests according these changes.

I provide below a minimal working script to test the new implementation.

import espressomd
from espressomd import reaction_methods
import numpy as np
import espressomd.visualization
import threading


BOX_LENGTH=10

system = espressomd.System(box_l=[BOX_LENGTH,]*3)

types = {"A": 0, "B": 1, "C": 2}
N_types= {"A": 1, "B": 1, "C": 2}
particle_types_to_move=[2]
N_MC_STEPS=1

system.part.add(type=[types["A"]]*N_types["A"], pos=np.random.rand(N_types["A"], 3) * BOX_LENGTH)
system.part.add(type=[types["B"]]*N_types["B"], pos=np.random.rand(N_types["B"], 3) * BOX_LENGTH)
system.part.add(type=[types["C"]]*N_types["C"], pos=np.random.rand(N_types["C"], 3) * BOX_LENGTH)

RE = espressomd.reaction_methods.ReactionEnsemble(kT=1, exclusion_range=0, seed=12)


visualizer = espressomd.visualization.openGLLive(system)

def main_thread():
    while True:
        RE.do_particle_displacement_MC_move(mc_steps=N_MC_STEPS, particle_types_to_move=particle_types_to_move)
        visualizer.update()

t = threading.Thread(target=main_thread)
t.daemon = True
t.start()
visualizer.start()

[jonaslandsgesell]

How many particles are moved per MC-step? I would make it a parameter of the function.

[pm-blanco]

@jonaslandsgesell only one particle is moved in each MC step. You can attempt to move as many particles as you want within your sampling cycle by just doing more MC steps:

for sample in range(N_samples):
      RE.do_particle_displacement_MC_move(mc_steps=N_particles)

Optionally, one could add an extra parameter to select how many particles you move in each new configuration before the Metropoli acceptance criteria but I do not see any practical use of it since basically would reduce the efficiency of the MC sampling of any system that I can think of.

[jonaslandsgesell]

@pm-blanco thanks for your answer :) The more particles are moved per accepted MC step the more the system decorrelates per accepted MC step. I expect a system dependent sweet spot in the number of uncorrelated samples which can be obtained per (realtime) hour. This sweet spot will depend on the force field, the particle density etc. and the a) number of particles to be moved per single MC step and b) the needed compute time per MC step. Since the compute time for the new electric energy computation (in Espresso) is essentially independent from moving one or 10 particles there will be a system dependent sweet spot. The explicit loop over 1-particle MC moves will decorrelate the system slower (i.e. needing more accepted MC steps) while potentially having a higher acceptance rate (measured per MC step).
If (on average) the compute time per accepted N-particle MC-step is lower than the compute time per N accepted 1-particle MC steps, then the N-particle MC step is favourable with regards to compute time. This can happen because in the one N-particle MC step there is 2 energy computations involved vs. In the case of N accepted 1-particle moves there are of the order N energy computations needed.

[pm-blanco]

@jonaslandsgesell For displacement MC moves, single particle moves are actually more efficient than N particle moves and they are specially recommended for simulations of condensed matter (See https://www.epfl.ch/labs/lcbc/wp-content/uploads/2019/03/MDMC_3.pdf, Section 3.4.1 pages 30 to 33). Of course, depending of the system, the efficiency of the MC conformational sampling can be improved by moving groups of particles using some smart scheme (for instance pivot moves for polymer chains). But this system specific MC moves should have their own function and they are, in my opinion, out of the scope of this general function.

[jonaslandsgesell]

@pm-blanco the argument presented in the book makes a lot of assumptions for arguing that the random N-particle displacement moves are less computationally efficient but agrees that exact performance has to be benchmarked.

To exactly decide if a random N-particle move can be more efficient than N 1-particle moves in Espresso I expect that we need to run benchmarks.

Let me make a theoretical example where my expectation is, that running an N particle move is more efficient.
In a system simulated with p3m and the amount |q| of charges per particle goes to 0 (i.e. vanishing electrostatic and other interactions also going to zero) we essentially simulate an ideal gas at the cost of running the p3m code and almost all particle moves will be accepted (be it the N particle move or the 1 particle moves). In this case the particle positions will decorrelate faster with the N particle move. I would be interested in benchmark results for "the number of moved particles per compute time" depending on the number of moved particles per MC step (1...10). In the case that the optimum is not at 1, but system dependent, we might want to document the result because this is relevant for big simulations.

I bring this topic up because the PR changes the existing behavior that the code can move multiple particles at once.

[RudolfWeeber]

My impression is that this could be modularized more:

there should be functions (or maybe functors as they can store parametesr and can be managed via script interface) to

• select a particle based on specific criteria

  • one that selects any particle
  • one that selects a particel from a list of given particle tyeps

• generate a new position for a particle

  • completely random (as is done now)
  • with an incremental position update by a step width ($x_1 =x_0 \detlta *r_u$, with $r_u$ a uniform random number)

• the check wehter the particle is with an exclusion radius should be a free funciton. Also, as I understnad the current function, the decision whether to use the efficient neighbor search is not automated. This can be decided by querying ghe max_range() of the cell system. In rhw lonf run, the entire functionality should go to the cell_system, I guess.

very rough sketch:

template <typename ParticleSelector, typename ParticleMoveGenerator>
single_parcile-move(int n_steps, const ParticleSelector& select, ParticleMovegenerator move_gen) {
auto last_energy = calc_potential_energy90
for (int i=0;i <mc_steps; i++)
   auto particle_id = select(); // reutrn a std::optioanl<size_t>
   if (! particle_id) return false; // no candidates
   old_pos = get_particle_position((*particle_id));
   auto new_pos = gen_move((*particle_id));
   move_particle((*particle_id), new_pos)

  if (! within_exclusion_radius(particle_id) {
     auto new_energy = calc_potential_energy()
     if (metropolis_accept(new_energy-old_energy, kT) {
       last_energy = new_energy;
       continue
     }
  }
   // else reject and unrolll move
   move_particle((*paritcle_id), old_pos);
}

Notes:

• as long as the code is not parallelized, the use of particle ids is not an issue, since the fetch cache will prevent unneeded communicaton ot other nodes.
• once the code is parallelized, only the node that has the selected particel will move it. After each step, it has to be checked via an mpi_reduce with a logical and whether a global resort is necessary (if the particle was moved beyond the node's spatial domain)
• Later on the full energy recalculaiton can be replacd by something more targeted
• the current id seleciton (to my understanding) loads all particle ids of the system into a vector. That doesn't scale. we can come up with sth better, once the code is modularized.

[davidbbeyer]

@pm-blanco @jonaslandsgesell In principle, I agree with Jonas that there may be situations where moving multiple particles might be more efficient. However, since I seem to be the only person using the displacement move feature (and I don't need to move multiple particles at once) and the choice was between removing it or fixing it, I would be in favor of leaving it as it is right now.