  We consider in this paper random batch interacting particle methods for
solving the Poisson-Nernst-Planck (PNP) equations, and thus the
Poisson-Boltzmann (PB) equation as the equilibrium, in the external unbounded
domain. To justify the simulation in a truncated domain, an error estimate of
the truncation is proved in the symmetric cases for the PB equation. Then, the
random batch interacting particle methods are introduced which are $O(N)$ per
time step. The particle methods can not only be considered as a numerical
method for solving the PNP and PB equations, but also can be used as a direct
simulation approach for the dynamics of the charged particles in solution. The
particle methods are preferable due to their simplicity and adaptivity to
complicated geometry, and may be interesting in describing the dynamics of the
physical process. Moreover, it is feasible to incorporate more physical effects
and interactions in the particle methods and to describe phenomena beyond the
scope of the mean-field equations.
