  The net charge of solvated entities, ranging from polyelectrolytes and
biomolecules to charged nanoparticles and membranes, depends on the local
dissociation equilibrium of individual ionizable groups. Incorporation of this
phenomenon, \emph{charge regulation}, in theoretical and computational models
requires dynamic, configuration-dependent recalculation of surface charges and
is therefore typically approximated by assuming constant net charge on
particles. Various computational methods exist that address this. We present an
alternative, particularly efficient charge regulation Monte Carlo method
(CR-MC), which explicitly models the redistribution of individual charges and
accurately samples the correct grand-canonical charge distribution. In
addition, we provide an open-source implementation in the LAMMPS molecular
dynamics (MD) simulation package, resulting in a hybrid MD/CR-MC simulation
method. This implementation is designed to handle a wide range of
implicit-solvent systems that model discreet ionizable groups or surface sites.
The computational cost of the method scales linearly with the number of
ionizable groups, thereby allowing accurate simulations of systems containing
thousands of individual ionizable sites. By matter of illustration, we use the
CR-MC method to quantify the effects of charge regulation on the nature of the
polyelectrolyte coil--globule transition and on the effective interaction
between oppositely charged nanoparticles.
