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, charge regulation (CR), 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 CR 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 large-scale Atomic/Molecular Massively Parallel Simulator 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 CR on the nature of the polyelectrolyte coil–globule transition and on the effective interaction between oppositely charged nanoparticles.

The formation of coacervates through phase separation of oppositely charged polyelectrolytes (PEs) is critical for understanding biological condensates and developing responsive materials. Traditionally, coacervates are viewed as spherical droplets with growth dynamics resembling liquid-liquid phase separation. However, our fluid particle dynamics simulations incorporating hydrodynamic and electrostatic interactions challenge this perspective. Here, we find that oppositely charged PEs form a percolated network even in semi-dilute solutions, coarsening with a unique growth law, ℓ ∝ t^1/2. This self-similarity, absent for neutral polymers in poor solvents, arises because PEs in good solvents exhibit weaker, longer-range attractions due to spatial charge inhomogeneity under global charge neutrality. This results in a lower density of the PEs-rich phase and reduced interfacial tension. Increased charge asymmetry further slows network coarsening. Additionally, coacervate droplets initially display irregular shapes due to weak interfacial tension, transitioning slowly to spherical forms. Our research provides new insights into coacervate morphology and coarsening dynamics.

We have developed an accurate and efficient method for molecular dynamics simulations of charged particles confined by planar dielectric interfaces. The algorithm combines the image-charge method for near field with the harmonic surface mapping, which converts the contribution of infinite far-field charges into a finite number of charges on an auxiliary spherical surface. We approximate the electrostatic potential of far-field charges via spherical harmonic expansion and determine the coefficients by fitting the Dirichlet-to-Neumann boundary condition, which only requires the potential within the simulation cell. Instead of performing the direct evaluation of spherical harmonic series expansion, we use Green’s second identity to transform the series expansion into a spherical integral, which can be accurately represented by discrete charges on the sphere. Therefore, the fast multipole method can be readily employed to sum over all charges within and on the sphere, achieving truly linear $\mathcal{O}\left(N\right)$ complexity. Our algorithm can be applied to a broad range of charged complex fluids under dielectric confinement.

Ionic distributions near charged interfaces control processes from colloidal aggregation to solvent flow in nanodevices. Such interfaces are often characterized by a jump in the permittivity, which gives rise to the surface polarization charge. This induced charge may significantly affect the ionic distributions so that efficient methods for modeling spatially varying dielectrics are needed. We formulate a method with $\mathcal{O}\left(N\log N\right)$ scaling for electrolytes between charged planar interfaces with asymmetric dielectric contrasts. Our approach, which builds on earlier work, is based on combining image charges with the particle–particle particle–mesh algorithm and representing uniform surface charges via an electric field. This enables simulations of complex dielectric interactions that outperform most alternative methods in speed and accuracy. To make the method practically useful, we provide guidelines—based upon careful tests—for choosing optimal simulation parameters. Explicit expressions for the electrostatic forces are given to facilitate the implementation of our algorithm in standard molecular dynamics packages.