A generalized grand-reaction method for modeling the exchange of weak (polyprotic) acids between a solution and a weak polyelectrolyte phase.

We introduce a Monte-Carlo method that allows for the simulation of a polymeric phase containing a weak polyelectrolyte, which is coupled to a reservoir at a fixed pH, salt concentration, and total concentration of a weak polyprotic acid. The method generalizes the established grand-reaction method by Landsgesell et al. [Macromolecules 53, 3007-3020 (2020)] and, thus, allows for the simulation of polyelectrolyte systems coupled to reservoirs with a more complex chemical composition. In order to set the required input parameters that correspond to a desired reservoir composition, we propose a generalization of the recently published chemical potential tuning algorithm of Miles et al. [Phys. Rev. E 105, 045311 (2022)]. To test the proposed tuning procedure, we perform extensive numerical tests for both ideal and interacting systems. Finally, as a showcase, we apply the method to a simple test system that consists of a weak polybase solution that is coupled to a reservoir containing a small diprotic acid. The complex interplay of the ionization of various species, the electrostatic interactions, and the partitioning of small ions leads to a non-monotonous, stepwise swelling behavior of the weak polybase chains.