Specify insertion proximity in reaction ensemble

[pw0908]

Hi Espresso Developers,

I'm currently trying to run some simulations involving the reaction ensemble where I have a reaction:

$$ \mathrm{HA} \rightarrow \mathrm{A}^- + \mathrm{H}^+$$

To give a more realistic representation of the reactions, during the monte carlo steps, I would like the proton to be inserted / removed in some proximity to the $\mathrm{HA}$ / $\mathrm{A}^-$ particle with which it's reacting. I've taken a look at the source code and, as far as the insertion is concerned, this doesn't seem too challenging to implement. The deletion does appear to be more difficult.

Would you have any guidance for how to proceed?

Thanks!

[jngrad]

Hi!

Our MC experts @kosovan and @pm-blanco may have some idea of how to implement this functionality.

Before you invest time on modifying the code, please be aware we have rewritten the reaction ensemble in Python to make it easier to implement such changes. The original C++ code will be completely replaced. This refactoring project is currently on pause, presumably until after we release ESPResSo 5.0 (in a few weeks), and is kept in a different repository because it is not yet in a state suitable for non-developers.

[pw0908]

Thanks for your reply!

While waiting for the MC experts to reply, where is this python implementation being developed? I don't think this will be a big change and I'd rather do it sooner rather than later. Of course, I would much prefer it be a contribution that will be useful long-term, so if I can help set-up this implementation such that it will be used in the refactored version of espresso, I would rather do that.

[kosovan]

Hi,

It is possible to bias the reaction algorithm, such that reactions occur only between reactants within a given region. For details, see e.g. our recent paper with @pm-blanco https://doi.org/10.1063/5.0226122

However, such implementation will not give you a "more realistic representation of the reactions", as you probably expected. The whole point of using MC is that it allows you to take an arbitrary path through the configuration space and, as long as you preserve detailed balance, it always converges to the correct equilibrium distribution. The path you take affects the rate of convergence, so you can choose the one that converges fastest.

If you implement such bias for acid-base reactions, as suggested, then you need to introduce a correction for the bias in the acceptance rate. This comes at extra computational overhead and increased complexity of the code. So, unless a random reactant from any location in the box poses a technical problem or restricts the sampling efficiency, then it is usually the best choice. See the paper cited above for an example where it was not the best choice. Another example is already present in the code: we prevent the reactions for occurring below a certain threshold radius. In this case, the correction is straightforward, so there is no overhead.

If you are interested in modeling reaction kinetics, then it is a different story. You have to start with postulating the reaction mechanism which you want to represent by a stochastic process. Then, the main problem is postulating a realistic mechanism, which requires experimental verification.