Factors affecting the calculation of mean ionic activities in Monte Carlo simulations of primitive electrolytes

Abstract

This paper reports canonical ensemble Monte Carlo (CEMC) simulations of restricted primitive model (RPM) electrolytes (charged hard sphere ions of equal charge and radius in a dielectric continuum) in 2D and 3D non-Euclidean geometries. The effect of the maximum possible displacement for particle movements, Δx, was investigated systematically. As the maximum possible displacement increased, the energies of particle deletion and re-insertion became uncorrelated. The Widom deletion formula was insensitive to the magnitude of particle displacement but for the Widom insertion method the values of ln γ_± depended significantly upon the maximum possible displacement. The values of ln γ_± calculated from the Widom insertion formula were also profoundly sensitive to the inclusion or exclusion of cases of hard sphere overlap (HSO) in the ensemble averages. When the particle to be deleted or re-inserted interacted with a charged system, the ln γ_± values were affected. Attractive interactions (i.e. the remainder of the system had an opposite charge) decreased the values of ln γ_±, while repulsive interactions increased the ln γ_± values. Electroneutrality was maintained by the use of a notional or ‘‘ghost’’ point charge of the same charge as the ion to be deleted or inserted, placed temporarily at the opposite pole of the (hyper)sphere during the calculation of the energies of interaction of an ion with the rest of the ensemble during deletion and insertion. Based upon the above results, a new method for the calculation of ln γ_± values was proposed: (i) use of large random displacements during re-insertion; (ii) exclusion of cases of HSO from the ensemble averages; and (iii) use of a notional or ‘‘ghost ’’ point charge to maintain electroneutrality. When these conditions were used, values of ln γ_± calculated from the deletion and insertion formulae agreed within the simulation uncertainties for all system sizes (N) and the values did not depend upon N. The ln γ_± values calculated from the insertion and deletion formulae also agreed within the statistical uncertainties at all concentrations and they agreed exceptionally well with the values obtained from the ln(f(ϕ)/g(ϕ)) vs. βϕ plots. This proposed method was suitable for the calculation of chemical potentials of electrolytes up to volume fractions of 0.2, although it may be suitable for higher concentrations and volume fractions.