Dear all,
Due to a mistake in my code, I have learned that there is a small bug in the constant pH algorithm. I have observed that this bug happens in the actual development version of espresso. Consider the typical acid/base equilibrium of a weak acid
HA <-> A- + H+
HA: type 0
A-: type 1
H+: type 2
Consider that, by mistake, the user does not create any particle of types 0 and 1, but it creates particles of type 2. In these conditions, it is of course impossible to have the reaction in either direction. When attempting to do a reaction step, espresso enters in an infinite loop-like behaviour that needs to be killed from the terminal. Below I attach a code with a minimal setup that reproduces this bug. I have also noticed that if the user does not create any particle at all the program returns a segmentation fault error. These errors do not happen in using the reaction ensemble algorithm, which simply does not do the reaction if it is not possible.
In my opinion, this bug should be easily solved by checking the number of particles of each reactant at the beginning of the constant pH algorithm and exiting the subroutine if the reaction is impossible. I will work on a bugfix myself and submit it for your consideration.
Best,
Pablo M. Blanco
Test code with the described bug:
import espressomd
from espressomd import reaction_ensemble
system = espressomd.System(box_l=[10] * 3)
charge={0:0, 1:1, 2:-1}
typ=2
system.part.add(id=[0], pos=[[0,0,0]], type=[typ], q=[-1])
RE = reaction_ensemble.ConstantpHEnsemble(temperature=1, exclusion_radius=1, seed=1)
RE.constant_pH = 7
RE.add_reaction(gamma=1, reactant_types=[0], reactant_coefficients=[1],
product_types=[1,2],
product_coefficients=[1, 1], default_charges=charge)
print("Before reaction")
RE.reaction(1)
print("After reaction")
--
Dr. Pablo M. Blanco
Department of Physical and Macromolecular Chemistry;
Faculty of Science
Charles University
Faculty of Science
Albertov 6, 128 00 Praha 2