Binary mixtures of asymmetric continuous charge distributions: molecular dynamics simulations and integral equations.

Binary mixtures of asymmetric continuous charge distributions: Molecular dynamics simulations and integral equations D. M. Heyes and G. Rickayzen Citation: The Journal of Chemical Physics 142, 074904 (2015); doi: 10.1063/1.4908046 View online: http://dx.doi.org/10.1063/1.4908046 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/7?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Molecular Dynamics Simulation of Smaller Granular Particles Deposition on a Larger One Due to Velocity Sequence Dependent Electrical Charge Distribution AIP Conf. Proc. 1415, 209 (2011); 10.1063/1.3667258 A general formulation for the efficient evaluation of n-electron integrals over products of Gaussian charge distributions with Gaussian geminal functions J. Chem. Phys. 134, 244115 (2011); 10.1063/1.3600745 Electrostatic potential of mean force between charged bovine serum albumin molecules in aqueous NaCl solutions by hypernetted-chain integral equation J. Chem. Phys. 117, 407 (2002); 10.1063/1.1481380 Comparison between integral equation method and molecular dynamics simulation for three-body forces: Application to supercritical argon J. Chem. Phys. 114, 5674 (2001); 10.1063/1.1350643 Self-consistent integral equation theory for polyolefins: Comparison to molecular dynamics simulations and xray scattering J. Chem. Phys. 114, 2847 (2001); 10.1063/1.1338505

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 149.150.51.237 On: Thu, 02 Apr 2015 19:17:58

THE JOURNAL OF CHEMICAL PHYSICS 142, 074904 (2015)

Binary mixtures of asymmetric continuous charge distributions: Molecular dynamics simulations and integral equations D. M. Heyes1,a) and G. Rickayzen2,b)

1 2

Department of Physics, Royal Holloway, University of London, Egham, Surrey TW20 0EX, United Kingdom School of Physical Sciences, University of Kent, Canterbury, Kent CT2 7NH, United Kingdom

(Received 11 December 2014; accepted 2 February 2015; published online 20 February 2015) An investigation is carried out of the association and clustering of mixtures of Gaussian charge distributions (CDs) of the form ∼Q exp(−r 2/2α 2), where Q is the total charge, r is the separation between the centers of charge and α governs the extent of charge spreading (α → 0 is the point charge limit). The general case where α and Q are different for the positive and negatives charges is considered. The Ewald method is extended to treat these systems and it is used in Molecular Dynamics (MD) simulations of electrically neutral CD mixtures in the number ratios of 1:1 and 1:4 (or charge ratio 4:1). The MD simulations reveal increased clustering with decreasing temperature, which goes through a state in which each large CD is overlapped by four of the oppositely signed CD in the 1:4 case. At very low reduced temperatures, these mini-clusters progressively coalesce into much larger tightly bound clusters. This is different from the 1:1 mixture case, where the low temperature limit is a random distribution of neutral dimers. At higher temperatures, the MD radial distribution functions g(r) agree well with those from the hypernetted chain solution of the Ornstein-Zernike integral equation, and (at not too high densities) a previously introduced mean field approximation extended to these charge distribution systems. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4908046]

I. INTRODUCTION

Electrostatic interactions play an important role in all branches of the physical sciences.1 One such area is the behavior of large macroions (e.g., proteins and colloids) of the same charge sign in a compensating background of counterions, which forms an electroneutral system. For dilute solutions, the Debye-Hückel (DH) or screened coulomb effective pair potential between the macroions, which is repulsive at all separations, can be used. With increasing concentration screening effects arising from positional correlations between the two species which go beyond a mean-field level become important; it has been argued that these can lead to a long range attraction between macroions with the same charge sign,2–5 although it appears, this has still not been fully resolved (see, for example, Refs. 6 and 7 for a contrary view). More generally, it may be supposed that clustering of the macroions can arise through the many-body correlations between the spatial distribution of macroions and counterions. Clustering of model charged particle systems near the critical point has been investigated. For ionic solutions, as the critical temperature (Tc) is approached, the crossover from meanfield-like to Ising criticality is delayed compared to that for solutions of van der Waals molecules; this may be due to ion pairing or clustering in the subcritical state.8,9 Also Monte Carlo simulations of the Restricted Primitive Model (RPM) model ions reveal that at temperatures below Tc, the vapor phase consists almost exclusively of strongly held together a)Electronic mail: [email protected] b)Electronic mail: [email protected]

0021-9606/2015/142(7)/074904/10/$30.00

ion pairs.10 Significant ion pair dissociation occurs very close to Tc, which suggests that there are significant fluctuations in the number of clusters very close to Tc. This work is concerned with providing a computationally manageable description of macroion-counterion systems in which there can be a large asymmetry in size and net charge between the two species. As an example, a soft macroion such as a polyelectrolyte can be represented at a coarse grained level using a charge density if the molecule is represented by the Gaussian form ∼ exp(−r 2/2α 2), where r is the distance of a volume element from the center of charge and α determines the breadth of the charge distribution (CD) (α → 0 is the point charge limit). This analytic form gives rise to a bounded effective pair potential which has been modelled a number of times.11–16 For an overall neutral binary mixture of oppositely charged charge distributions, there is no need for an additional short range repulsion, as unlike in the point charge limit, these charge distribution systems are thermodynamically stable according to the Fisher-Ruelle stability criteria.17 Let the total charges of the positively and negatively charged distributions be Q+ and Q−, respectively. In the special case of Q+ = −Q− = Q, the classical ground state average potential energy per particle √ for α > 0 in the zero temperature limit is −V (0)/2 = −Q2/2α π, where V(0) is the value of the pair potential energy at zero separation. An effective reduced temperature can be defined as T ◦ = k BT/V (0) which is therefore proportional to α. If α (and therefore T∗) is small enough the attraction between the oppositely charged distributions overcomes the thermal energy, and as revealed by integral equation and molecular simulation studies, the system separates into almost neutral dimers.11–16,18

142, 074904-1

© 2015 AIP Publishing LLC


074904-2

D. M. Heyes and G. Rickayzen

J. Chem. Phys. 142, 074904 (2015)

We have recently reported preliminary results for the mean spherical approximation (MSA) and hypernetted chain (HNC) closure of the Ornstein-Zernike (OZ) equations and carried out Molecular Dynamics (MD) simulations of a 1:1 binary mixture of oppositely charged Gaussian CDs with the same value of α.21 The nature of the clustering and dynamical features of the dimers were investigated. The focus here is on mixtures of unequal charges with arbitrary charge spreading functions. In the more general case of a 1:m relative number density of charges ratio, where m > 1, the efficacy of the integral equation solutions and the nature of the clustering is less well established as m increases (the charge ratio is m:1). Although some of the final form equations derived herein were obtained by Coslovich et al.,11 their comparison with simulation in that work was carried out using only symmetric Gaussian distributions. More general 1: m cases are treated here using Ornstein-Zernike integral equations and molecular dynamics simulation. In addition, a Mean Field approximation (MFA) which is an alternative to the Ornstein-Zernike routes applied by us to Gaussian and bounded soft sphere systems in Ref. 19 is extended here to treat asymmetric charge systems. The Ewald formulas applicable to an arbitrary binary mixture of Gaussian charge distributions used in the MD simulations are derived here, in the Appendix.

the OZ equation for mixtures is   hα β (r) = cα β (r) + ργ d 3r ′cαγ (r′ − r′) hγ β (r′) ,

(2.7)

γ

where ργ is the number density of the γ-component. As a consequence of Eq. (2.6), there are just three independent O-Z equations. These equations must be supplemented by the HNC closure relations which in this case are gα β (r) = exp − βVα β (r) + hα β (r) − cα β (r) . (2.8) The 12 and 21 terms for h and c are equal. The HNC closure is used here as this has been found to be successful in representing soft bounded particle interaction for 1:1 systems.19–21 Some of the technical details relevant to this work are given in our previous article on a Gaussian charge distribution extension of the One Component Plasma (OCP).22 The energy per unit volume is23  ∞ u=π ρα ρ βVα β (r) gα β (r) r 2dr. (2.9) α, β

0

This expression excludes the self-energies of the individual particles.

II. THEORY OF UNEQUAL CHARGES

A. HNC

Consider a mixture of two species of Gaussian charge distributions; the different species carry a charge Qi, has a spread of αι and a number density ρi, for i = 1,2. Species i has a charge spherically symmetric density (relative to its center) of

Because of the long range nature of the potential in Eq. (2.4), the procedure of Springer et al. originally derived for the OCP24 is extended to apply to the HNC equations for these systems. In this scheme, the pair potential is split into a short-range and long-range part, i.e., V (r) = Vs (r) + Vl (r), where

ρi (r) ∝ e−r

2/2α 2 i

,

(2.1) Vs (r) = Q2

whose Fourier Transform (FT) is 2 2/2

ρˆ i (k) = Q i e−α i k

.

(2.2)

(2.3)

The interaction energy of two such unequal charge distributions is  sin (kr) 2Q1Q2 ∞ dk exp − α12 + α22 k 2/2 V12 (r) = π kr 0 ( ) 2 2 1/2 α + α2 Q1Q2 r + , = erf ; β=* 1 (2.4) r 2β , 2 with FT 4πQ1Q2 Vˆ12 (k) = exp − β 2 k 2 . (2.5) 2 k The effective pair potential in Eq. (2.4) is bounded so the charge distributions can interpenetrate and pass through each other. In general, there are three independent pair distribution functions, namely, g11 (r) ; g22 (r) ; g12 (r) = g21 (r) ;

(2.10)

and 1 − e−pr . (2.11) r The parameter p in Eqs. (2.10) and (2.11) is an arbitrary decay parameter and is chosen to be the order of inverse decay length of the radial distribution function. The FT of these functions is   exp(−k 2α 2) p2 (2.12) V˜s (k) = 4πQ2 − k2 k 2(k 2 + p2) Vl (r) = Q2

If the fluid is to be neutral overall, one must have ρ1Q1 + ρ2Q2 = 0.

erf (r/2α) − 1 + e−pr r

(2.6)

and V˜l (k) = 4πQ2 p2/k 2(k 2 + p2),

(2.13)

respectively. The FT of the expressions in Eq. (2.7) is h˜ α β = c˜α β +

2 

ρδ c˜αδ h˜ δ β ,

(2.14)

δ=1

where ρi is the number density of species i. Now define γ(r) = h(r) − c(r), without subscripts, where the indices 11, 12, or 22 are implied. Then if γ˜ s = γ˜ − V˜l β and c˜s = c˜ + V˜l β, where β = 1/k BT and


074904-3


J. Chem. Phys. 142, 074904 (2015)

s l γ˜11 = h˜ 11 − c˜11 − V˜11 β s l ˜ ρ1c˜11c˜11 − V11 β + ρ2c˜12 h˜ 12 = , (2.15) 1 − ρ1c˜11 where Eq. (2.14) has been used in the second line. Similarly, l s β = h˜ 12 − c˜12 − V˜12 γ˜12 s l ˜ ρ1c˜11c˜12 − V12 β + ρ2c˜12 h˜ 11 = , 1 − ρ1c˜11

exp(− βVisj (r)

As (2.16)

γisj (r)).

from which gi j (r) = + The iterative procedure for solving these equations is given in Ref. 25. B. MSA

The random phase or MSA is another closure of the OZ equation, which has been applied by Coslovich, Hansen, and Kahl to these systems.11 This amounts to assuming that the direct correlation functions have their asymptotic forms throughout space, and therefore,11,15 ( ) Q2i r ; kT cii (r) = −Vii (r) = − erf r 2α i ( ) Q1Q2 r kT c12 (r) = −V12 (r) = − erf ; r 2β (2.17) 4πQ2i kT c˜ii (k) = − 2 exp −α 2i k 2 ; k 4πQ1Q2 kT c˜12 (k) = − exp − β 2 k 2 , k2 where i is 1 or 2. The FT’s of the pair correlation functions are then obtained as solutions of three Ornstein-Zernike Eq. (2.7) with the FT of the total correlation functions given by (1 − ρ2c˜22 (k)) c˜11 (k) + ρ2c˜21 (k) c˜12 (k) h˜ 11 (k) = , D (k) c˜12 (k) h˜ 12 (k) = h˜ 21 (k) = , (2.18) D (k) (1 − ρ1c˜11 (k)) c˜22 (k) + ρ1c˜12 (k) c˜21 (k) h˜ 22 (k) = , D (k) where D (k) = (ρ1c˜11 (k) − 1) (ρ2c˜22 (k) − 1) − ρ1 ρ2c˜12(k)2,

(2.19)

and in the MSA, c˜α β (k) = − βV˜α β (k)

where Vi j (r) is the electrostatic potential between particles of species i and j, and vi j (r) is the corresponding self-consistent potential. The pair distribution functions are related to these potentials by gi j (r) = exp − βvi j (r) . (2.22)

(2.20)

as long as the determinant of coefficients D(k) is not zero. Because the direct correlation function is symmetric to swapping α and β in Eq. (2.20), it follows that h12 = h21 and the same equally applies to their FT. C. MFA

In the MFA,19 for mixtures, the closure relations of Eq. (2.8) are replaced by the self-consistent equations vi j (r) = Vi j (r)   + ρk dr′Vik (r − r′) exp − βvk j (r′) − 1 , k

(2.21)

Vi j = Vj i ,

(2.23)

it follows that vi j = v j i ;

gi j = g j i .

(2.24)

Equation (2.21) can be solved by Fourier transformation and using iteration. The Picard method was used to iterate to a converged solution, which is sufficient for these soft systems (see Ref. 23). The advantage of the MFA formulation is that it is physically more intuitive than the HNC approach and does not require a Springer-like decomposition of the potential terms in Eq. (2.21) as the self-consistent potential vij is much shorter ranged because of the charge screening, and there is an analytic solution for the FT of Vij. The solution of the electrostatics of these systems as Ewald formulas which are used in the molecular dynamics simulations is derived in the Appendix.

III. RESULTS AND DISCUSSION

Molecular dynamics simulations were carried out using the generalized Ewald summation method derived in the Appendix. The artificial charge screening parameter defined in the Appendix is set to γ = L/16, where L is the side length of the cubic simulation box. The upper limit on the magnitude of the reciprocal lattice vector, K, satisfied, K 2 L 2/2π 2 ≤ 99. Two sets of calculations were carried out, the first, S1, was for a 1:1 system, and the second, S2, for a 1:4 mixture. The unit of charge is taken to be the total charge on one “particle” of species 1. The charges were Q2 = −Q1 = −1 for S1 and Q2 = −Q1/4 = −1/4 for S2. The positive charge was taken to be Q1. The unit of distance, σ, is twice the Wigner-Seitz radius rs = (3/4πρ1)1/3 of species 1, where ρ1 is the number density of species 1. The unit of energy is ε = Q12/σ, and temperature is T ∗ = k BT/ε, as for our previous reports on the soft OCP20 and the 1:1 equimolar mixture of Gaussian charge distributions.21 The a charge decay distance parameter is taken to be proportional to the square root of Q, to be consistent with ideal chain statistics (Q is assumed to be proportional to the molar mass of a polyelectrolyte molecule). For species 1, the parameter a was σ/16 and for species 2 it was, σ/32, with L being typically ∼6σ. There were N1 charge distributions of species 1 and N2 of species 2 in the simulation cell. The reduced packing fraction is defined to be ζ1 = π N1/6L3 and ζ2 = π N2/6L3. The total charge of the system was zero. The Verlet leapfrog algorithm21 was used to integrate Newton’s equations of motion with a time step of 0.002 a1/2/T 1/2 to accommodate a decrease in the time step with decrease in the range of the charge (parameter a in the Appendix). The simulations were carried out for 105–106 time steps. We concentrate on 1:4 mixtures where the asymmetry is


074904-4


sufficiently large to exhibit noticeable differences in physical and structural properties from the previously studied 1:1 case without suffering too much from statistical noise in the properties involving the larger species (of which there are fewer than the other species) with the computer resources available to us. The state points covered here are at densities lower than the critical point value for the 1:1 system and of order or lower than Tc as given in Refs. 11 and 18. For the 1:4 case, one expects a lower Tc than for 1:1 as this system is approaching the soft OCP which does not have a critical point (i.e., Tc tends to zero in 1: ∞ limit). The simulations were carried out in the single phase domain which has been shown previously to manifest only the clustering transformation with decreasing temperature. Figure 1 shows the radial distribution functions, g(r), for the like and unlike sign charge distributions of the 1:1 system with ζ1 = ζ2 = 1/2 and T = 20, at which temperature there is no dimerization of the oppositely charged distributions. There were 150 CDs of each species. The HNC closure of the OZ equation passes through the MD symbols within their statistical uncertainty. The MSA solution underestimates the like charge g(r) (i.e., for ++ and −− combinations) and the +− g(r). The MSA solution of the OZ equation is symmetric about unity and therefore not able to reproduce the asymmetry about unity of the (++,−−) and (+−) radial distribution functions for strongly correlated CD state points. The MD statistics deteriorate in the r → 0 limit because of the small volume element there. The DH analytic solution for g(r) agrees well with MD from about r > 0.2. The MFA and HNC solutions are too close to distinguish between them on the scale of the figure. The corresponding functions for the 1:4 system (i.e., four negatively charged CD for every positively charged one) are shown in Fig. 2, where ζ1 = 1 and T = 20. There were 200 positively charged distributions and 800 of the negatively charged species in the MD simulation. The HNC curve again goes through the simulation data points for the three component g(r)s. The MSA solution for the ++ g(r)

FIG. 1. Radial distribution function for a 1:1 system at a reduced temperature, T = 20 and ζ = 1. (ζ1 = ζ2 = 1/2) Key: MD (symbols), the OrnsteinZernike equation with closures: HNC (blue lines) and MSA (red lines). The DH approximation is also shown (green line). The MFA curves are also given, but on the scale of the figure, the HNC and MFA curves are so close they are hard to distinguish.

J. Chem. Phys. 142, 074904 (2015)

underestimates the MD data, but is in excellent agreement with MD and HNC for the −− and +− curves. The deviation from the MSA solution of the + + component of g(r) is greatest because the pair potential for these species is the largest, and expanding the exponential is less accurate. The MSA radial distribution functions involving the negatively charged CD are in closer agreement with the MD data (i.e., the “high temperature” limit) because the charges are smaller. Figure 2 also shows that the MFA solution is hardly distinguishable from the HNC solution, as may be seen also in Fig. 1, for all species components of the radial distribution function. Figure 3 compares the absolute values of the selfconsistent potential vij and the bare pair potential Vij for the S2 data at T = 20 (note the lin-log scale). The upper three curves are for Vij and the lower three curves with symbols are for vij. The HNC prediction for vij, which in reduced units is −T ln(g(r)), is also shown on the figure as continuous lines, which are seen to be hardly distinguishable from the MFA data points. The relationship between the three closures can be brought out by using the OZ equation to rewrite the HNC closure as gi j (r) = exp − βVi j (r) + hi j (r) − ci j (r)  2      (r) = exp  − βV + ρ dr′cik (r − r′) hk j (r′)  i j k .  k=1   (3.1) If in this expression the direct correlation is approximated by its asymptotic form, ci j (r) = − βVi j (r) .

(3.2)

Equation (3.1) becomes the MFA closure and the OZ equation becomes the MSA closure. At high temperatures, the replacement Eq. (3.2) leads to the agreement of the 3 closures with each other and with the simulation. At lower temperatures, the MSA differs from MD, the other two closures still agree with each other, and MD except at very small values of r (see Fig. 2), which are probably due to statistical fluctuations in the MD data in that region (the volume element tending to zero in this limit). Also as density increases for ρ1 up to 15, differences

FIG. 2. As for Fig. 1, except that a 1:4 mixture is considered. The total charges of the two types of charge distribution are Q+ = Q1 = 1 and Q2 = Q− = −1/4 charges. The charge spreading lengths, a and b, for + and − charged species, respectively, are in the ratio: 1 to 1/2.


074904-5


FIG. 3. The effective pair potential between the various species combinations obtained by iteration from (a) the HNC solution of the OZ equation and (b) the MFA. The 1:4 system with ρ1 = 1 and ρ2 = 4, and T = 20 was used. The top three curves, labelled V++, V−−, and V+−, are the true pair potential (Vij in Eq. (2.21)). The corresponding mean field potentials, vij are the lowest three curves on the figure. Again the HNC and MFA curves are hardly distinguishable. Note the lin-log scales.

J. Chem. Phys. 142, 074904 (2015)

FIG. 5. Projections of the coordinates of the two species on one of the faces of the MD cell at one time step. The positively charged species are shown as blue circles online and the negatively charged species are red online. Key: A 1:1 mixture for T = 1, ρ1 = 1, and ρ2 = 1 is shown.

between the MFA and HNC become evident, with the HNC results being closer to those produced by MD. Figure 4 shows the component radial distribution functions for a 1:10 mixture, where ρ1 = 1.0 and T = 5 using the integral equations (this large asymmetric ratio would for us be impractical to employ in MD). The three components of g(r) are given for HNC and MFA, and again the agreement is very good. The soft OCP solution which treats the negatively charged counterions as forming a continuous background22 is also shown for the ++ pair and may be seen to agree well with the HNC and MFA data where the counter CDs are explicitly included. The MSA agrees well with HNC and MFA for the −− component, but poorly for ++, for reasons discussed above. In contrast, the function, g(r) = exp(gMSA(r) − 1) which is labelled, MSA-HNC on the figure, does agree well with the HNC and MFA results for all three species components. As temperature decreases below about 5, differences between the HNC and MFA radial distri-

bution functions in the soft OCP model become progressively larger. This is in a temperature range where convergence of the HNC radial distribution becomes problematic in the explicit CD mixture case, presumably reflecting the onset of clustering in the MD system. Figure 5 shows the coordinates of the two species projected onto one side of the MD cell at the end of a 1:1 mixture simulation using T = 1, with the two species in a different color online. Even at this temperature, there is some evidence of a change of structure of the fluid; the non-random appearance of discs of clustering of the two species in chainlike structures in the projection is evidence for clustering in the bulk. This is accentuated at T = 0.1 in Fig. 6, in which quite densely packed clusters of the two species are evident. At T = 0.01, shown in Fig. 7, virtually all of the CD are paired up with another of opposite sign, forming a distribution of dimers which presumably are weakly interacting.

FIG. 4. The component radial distribution functions for a 1:10 mixture, where ρ1 = 1.0 and T = 5. The three component g(r)s are given for HNC and MFA. The soft OCP solution is shown for the ++ pair. The MSA solution, g MSA(r ) and the function, g (r ) = exp(g MSA(r ) − 1) which is labelled, MSA-HNC are also presented on the figure.

FIG. 6. As for Fig. 5, except that T = 0.1.


074904-6


J. Chem. Phys. 142, 074904 (2015)



Figure 8 shows the coordinates of the two species projected onto one side of the MD cell at the end of a 1:4 mixture simulation using T = 1. The larger discs are for species 1 and the smaller ones are for species 2, and there is some evidence of the onset of clustering of species 1, perhaps mediated by species 2. Figure 9 shows the corresponding display for T = 0.1, where there is clear evidence of a partitioning of the particles into small clusters. The basic cluster consisting of one of species 1 and up to four of species 2 is a common occurrence as are apparently more connected combinations of many 1 and 2 CD clusters. Not all of the small clusters are electroneutral however, as some of the positively charged species have less than four “counterions” associated with them. The trend continues for T = 0.05 where at that temperature the “1 + 4” cluster is the dominant structural unit. This is a neutral cluster, which is the nearest equivalent to the neutral dimer limit of the 1:1 mixture case.21 It is not unreasonable for these small clusters to be electrically neutral, which appears to be the case for many from the images. In

contrast to the 1:1 dimer, the 1 + 4 (and 1 + 3) cluster can be the building block of a more extended cluster even at lower temperature. Figure 10 for T = 0.01 shows that the system has transformed into a mixture of a few compact clusters coexisting with smaller clusters. The large clusters have a blue boundary (online), which indicates accumulation of the counterions in the middle of the large cluster. Because of the repulsion between the negatively charged species, it is presumably more likely that the 1 + 4 cluster will have a residual dipole and higher multipoles which would cause them to cluster together in groups. One might have to go to physically unrealistically low temperatures to suppress this effect (and produce neutral 1 + 4 units akin to the 1:1 “dimers”). The end point observed in Fig. 10 for T = 0.01 appears to have more in common with the collapsed state implied in Fisher and Ruelle’s theory of thermodynamic stability of pair potentials,17 (charged interactions are also treated in that publication), which we have explored for other model interaction potentials in Refs. 26 and 27. If the charge spreading function lengths are doubled, the

FIG. 8. As for Fig. 5, except that a 1:4 mixture is shown. The blue circles are for the positively charged distributions and the (smaller) red circles are for the negatively charged species. ρ1 = 1 and ρ2 = 4, and T = 1.



074904-7


FIG. 11. The average potential energy per particle, u, as a function of number density of the large charge (species 1) along several isotherms. The continuous lines are HNC results and the larger open circles are the MD data. The lower set of data is for the 1:1 mixture system. For the upper set, of 1:4 data, the points are multiplied by 5 and shifted upwards by 3.0, for clarity.

state of many isolated “1 + 4” (and similar) clusters is less well defined, and the system more readily forms the very large clusters seen in Fig. 10. With increasing number density, lower temperature is required to achieve the compact isolated cluster state; for 1:4, this is at a temperature of ca. 0.1 for ρ1 equal to 5 and ca. 0.01 for ρ1 equal to 15. Coslovich et al.11 (see their Eq. (17)) showed, using a thermodynamic stability criterion, that low temperature clusters are not expected for charge ratios greater than ca. 22:1 (or number density ratios greater than 1:22). Presumably, in the large ratio limit, the system goes over to the charge distribution version of the OCP considered in Ref. 22, as the counterions can then be considered to form a continuum background which ensures electroneutrality of the whole system. Figure 11 shows the average potential energy per particle, u, as a function of number density of the large charge (species 1) along several isotherms. The continuous lines are HNC integral equation results and the larger open circles are the corresponding MD data. Both the 1:1 and 1:4 mixtures are represented. The agreement between the two sets of data is very good up to ρ1 equal to 15, the highest value used in both treatments. In the vicinity of T ≈ 2, there are convergence problems with the integral equations which coincide with the start of well-defined clustering in the corresponding MD systems.

IV. CONCLUSIONS

The results of MD calculations are reported here of an asymmetric binary mixture of Gaussian charge distributions, in which the magnitude and distribution width of the charge can be arbitrarily different. As a methodological step in achieving this, the Ewald electrostatic interaction summation method is extended to encompass this more general coarsegrained representation of charged systems. For a 1:4 mixture

J. Chem. Phys. 142, 074904 (2015)

it, is found from the MD simulations that at not too small temperatures, when the system displays insignificant clustering, the HNC closure of the Ornstein-Zernike equation reproduces the local structure as represented by the radial distribution function very well, taking the MD results to be essentially exact. Also the MSA solution is very good for all component g(r) except for that between the highly charged macroions (++ in this case). An analytically quite different MFA gives numerically the same radial distribution functions as those produced by HNC in the high temperature limit (i.e., ca. T > 5) for number densities up to ca. 1-5. At lower temperatures, where the first signs of clustering become noticeable, the HNC solution shows convergence problems and small differences from the MFA solution. The soft OCP solutions of g(r) for the same systems converge excellently for HNC and MFA at high temperatures but show progressive differences in these functions for T lower than about 5. Instantaneous configuration snapshots for all of the two species in the simulation cell show that with decrease in temperature the 1:4 system forms a series of cluster states which are qualitatively different to those of the 1:1 system. In the 1:1 system in the T → 0 limit, previous studies have shown that the charge distributions self-assemble into essentially neutral dimers composed of overlapping positive and negatively charges (note there are no excluded volume interactions in the pair potential of these systems). In the 1 + 4 cluster case, there is a competition between the mutual repulsion between the small charge distributions (“counterions”) and their attraction for the larger charge distribution of opposite sign. One might expect the “1 + 4” building block clusters to have residual multipole moments which would cause them to attract and assemble into larger agglomerates. It is not inconceivable that at extremely low temperature states (which might be difficult to achieve in practice in a simulation), these clusters would break up into distinct 1 + 4 clusters. The effective interactions between the clusters could in principle be determined from the species-resolved radial distribution functions as was performed for 1:1 clusters in Ref. 11. The “counter anions” here are in fact coarse-grained themselves, so they may be considered to represent regions of excess charge (of opposite sign to the macroion) in the solvent. The advantage of the MD simulations carried out here is that they automatically take into account positional coupling between the model macroions and the counterions. The coarse grained representation of the distribution of counterions employed makes these asymmetric charge distributions computationally tractable for both molecular simulation and integral equation treatments.

APPENDIX: APPLICATION OF EWALD’S METHOD TO UNEQUAL CHARGE DISTRIBUTIONS 1. Electrostatic interactions of unequal charge distributions in a periodic lattice

This appendix considers further the formulas for the energy of a lattice of continuous charge distributions.28 Because the charge is assumed to be periodic in space standard


074904-8


J. Chem. Phys. 142, 074904 (2015)

procedures may be used to expand its FT in the series  Nα       ρ˜ 0 (K) = d 3reiK.r ρα r − r m α + R n    n  α=1,2 m α  cube  N α       = d 3reiK.r  ρ r − r α m α   mα α=1,2   all space =

Nα  

ρα (K) e

iK.r m α

.

(A1)

α=1,2 m α =1

Here, Nα is the number of molecules of species α in a unit cell, the vectors, Rn , are the lattice vectors, the vectors, K, are the reciprocal lattice vectors, and rm α are the co-ordinates of the particles of species α in the unit cell. Note that because a cell is neutral overall, ρ˜ 0 (0) is zero. Similarly, the electrostatic potential has the same periodicity as the charge density and can also be expanded in a similar Fourier series  φ (r) = L −3 φ˜ (K) exp (iK.r) , (A2) K

with 4π ρ˜ 0 (k) (A3) φ˜ (k) = k2 except that φ˜ (0) is, in the case of overall neutrality, an indeterminate constant which has no effect on the force nor the energy and may be taken to be zero. Thus, when K , 0, Nα     4π  . (K) φ˜ (K) = 2  ρ ˜ exp −iK.r α m α   K  α=1,2 m α =1 

(A4)

The total electrostatic energy of one cell is  1 1  d 3r ρ (r) φ (r) = ρ˜ 0 (K) φ˜ (−K) W= 2 2L 3 K,0 cell

ρ 0 (−K) 2π  ρ˜ 0 (K) ⌢ = 3 L K,0 K2 ⌢ 2π   ρ˜ α (K) ρ β (−K) = 3 L K,0 α, β=1,2 K2   ( ) × exp iK. rm α − rn β .

(A5)

m α, n β

Because ρ˜ 0 (0) = 0, there is no contribution to the energy from the term with K = 0. When the charges are Gaussian, we have, say, ρ1 (K) = Q1 exp −a2 K 2/2 , (A6) ρ2 (K) = Q2 exp −b2 K 2/2 , where the Gaussians have different spread parameters and charges. Hence, for K , 0, ρ0 (K) = Q1 exp −a K /2 2

2

N1 

eiK.rm 1

m 1=1

+ Q2 exp −b2 K 2/2

N2  m 2=1

eiK.rm 2.

(A7)

Similarly, when K , 0, exp −a2K 2/2  exp −iK.rm1 K2 m1 2 2 exp −b K /2  + 4πQ2 exp −iK.rm2 . (A8) 2 K m

φ˜ (K) = 4πQ1

2

From Eq. (A5) for Gaussian distributions, the energy per unit volume is ⌢ 2π   ρ˜ α (K) ρ β (−K) w= 3 L K,0 α, β=1,2 K2   ( ) × exp iK. rm α − rn β m α, n β

=

⌢ 2π   ρ˜ α (K) ρ β (−K) L 3 K,0 α, β=1,2 K2   ( ) cos K. rm α − rn β , ×

(A9)

m α, n β

with the charge densities given by Eq. (A6). In the last step, the fact that the terms in K and −K are complex conjugates has been exploited to write the result in its usual real form suitable for computation. When the spread of the charges, i.e., defined through the parameters, a and b, is sufficiently large, the series in Eqs. (A8) and (A9) converge rapidly and the sums can be performed without Ewald’s method.17 When one or both of these distances are small, however, it is necessary to use his method. We return to this case in Subsection 2 of the Appendix. The expression for the energy in Eq. (A5) includes half the energy arising from the interaction of a particle with itself because the potential is due to all the particles. In the case of point particles, this is an infinite energy. In the case of spreadout particles, it is finite. Nevertheless, it is not included in the theoretical expression in Eq. (A9). To compare the two we need to subtract this self-energy. The self-energy of one particle of species 1 is  1 Q21 V (0) = , (A10) π a with a similar expression for the other particles. Hence, half the self-energy per unit volume is   N2 N1 1 Q21 1 Q22 + . (A11) w0 = 2L 3 π a 2L 3 π b To obtain similar expressions in Cartesian co-ordinates, note that the potential at r, due to a particle of species 1 at rm, is ( ) |r − rm | Q1 erf . (A12) √ |r − rm | 2a The total potential at r is N1   r − rm1 − nL Q1 * + erf φ (r) = √ r − rm1 − nL 2a m 1=1 n , N 2   r − rm2 − nL Q2 +. erf * + √ r − rm2 − nL 2b m 2=1 n , (A13) Similarly, the energy per unit volume is


074904-9


J. Chem. Phys. 142, 074904 (2015)

N1  N1     rm1 − rl1 − nL 1   2   * +   er f Q1       2a r − r − nL   m l 1 1   m 1=1 l 1=1 n ,       N N   1 2      r − r − nL 1   1  m1 l2 + * . + 2Q Q er f w= √ 1 2 3   2L  rm1 − rl2 − nL 2 (a2 + b2) -    n m =1 l =1 ,   1 2       N N 2 2        r − r − nL 1 m2 l2   2   * +   + Q er f   2   2b r − r − nL m2 l2   m =1 l =1 n , 2

2

This includes half the self-energy of the interacting particles which can again be subtracted at the end.

2. The electrostatic potential and energy when the charge is periodic and α, β are small

In this case, the Fourier Series converge slowly and it is better to follow Ewald and split the potential into a shortranged potential, φ s (r) and a slowly varying long-range one, φl (r). This is accomplished by introducing a third potential due to a Gaussian charge distribution with a greater spatial spread, γ. There is some flexibility in the choice of this potential, and it is reasonable to assume that the potential arising from a mixture of N1 particles each carrying a charge

φ s (r) =

N1   m 1=1 n

+

Q1 |r − rm − nL|

N2   m 2=1 n

(A14)

of Q1 and N2 particles each carrying a charge of Q2 and both having a spread γ would give the fastest convergence; we use that potential here. The potential of the original particles is, say, denoted by φα (r) and the new one by φγ (r). Then, φl (r) = φγ (r)

(A15)

is chosen and φ s (r) = φα (r) − φ β (r) ;

φα (r) = φl (r) + φ s (r) . (A16)

The long-range potential can be calculated using Eq. (A13) but with a and b both replaced by γ. The short ranged potential may be calculated in the same way. The result is

 r − rm1 − nL  r − rm1 − nL + − erf * + erf * √ √  , 2a 2γ - ,

Q2 |r − rm − nL|

 r − rm2 − nL  r − rm2 − nL + * + . erf * − erf √ √  , 2b 2γ , -

(A17)

Similarly, the energy per unit volume should be split into two parts w = wl + w s ,

(A18)

where wl =

1 2L 3

 d 3r ρ0 (r) φγ (r) = cell

1  ρ0 (K) φγ (K) 2L 6 K

 exp − a2 + γ 2 K 2/2    2   exp iK. rm1 − rn1 Q     1 2   K     m , n K,0   1 1     2 2 2       exp − b + γ K /2   2     + Q exp iK. r − r m n   2 2 2 2    K   2π  m2, n 2 K,0  . = 6 2 2 2    exp − a + γ K /2   L    + Q1Q2 exp iK. rm1 − rn2      2   K     m , n K,0 1 2      2 2 2      exp − b + γ K /2         exp iK. r − r m n  + Q1Q2  2 1 2 K   m , n K,0 2 1

(A19)


074904-10


Finally, ws =

1 2L 3

J. Chem. Phys. 142, 074904 (2015)

 d 3r ρ0 (r) φα (r) − φγ (r) .

(A20)

The first term in the equation above is the sum of all the

ws =

N1  N1  Q21  1 2L 3 m =1 l r − r m1 l 1 − nL n

 r − rl1 − nL  r − rl1 − nL  + − erf * m1 +   erf * m1   2a 2c 1 , ,   1 1=1 N N 2  2   Q2  r − rl2 − nL   * rm2 − rl2 − nL + 1 + * m2   + 23 erf − erf  2b 2c 2L m =1 l =1 n rm2 − rl2 − nL  2 , ,   2 2 N1  N2    r − rl2 − nL r − rl2 − nL   Q1Q2 1 + − erf * m1 +   erf * m1 +  2a 2c 2L 3 m =1 l =1 n rm1 − rl2 − nL  1 , -  , 1 2 N2  N1   rm2 − rl1 − nL rm2 − rl1 − nL   Q1Q2  1 + − erf * + ,  + erf * 3   2b 2c2 2L m =1 l =1 n rm2 − rl1 − nL  , ,  2 1

5Y.-G.

where 1 2

c1 = a2 + γ 2 /2 ,

1 2

c2 = b2 + γ 2 /2 .

(A22)

Levin, Rep. Prog. Phys. 65, 1577 (2002). Ray and G. S. Manning, Macromolecules 30, 5739 (1997). 3X. Chu and D. T. Wasan, J. Colloid Interface Sci. 184, 268 (1996). 4R. Messina, C. Holm, and K. Kremer, Phys. Rev. Lett. 85, 872 (2000).

2J.

(A21)

Chen and J. D. Weeks, PNAS 103, 7560 (2006).

6J. Baumgartl, J. L. Arauz-Lara, and C. Bechinger, Soft Matter 2, 631 (2006).

Thus, if γ is chosen suitably, wl can be evaluated from the FT of Eq. (A19) and w s from the real space lattice summation of Eq. (A21). The self-energy of the α-particles is contained in the terms with mα = l α and n = 0 in Eq. (A21). These terms are    1 1  N1Q21 N2Q22 N1Q21 N2Q22  w0 = + − − . (A23) b c1 c2  2L 3 π  a  Hence, to remove the self-energy of the α-particles, we need only replace the diagonal terms in the sum by  1 1 * N1Q21 N2Q22 + − 3 + . (A24) π , c1 c2 2L This last term needs to be retained as it is included in wl . Note that in the case where b > a and b sufficiently large, one may choose γ = b. 1Y.

energies resulting from the interactions of the original Gaussian charges with each other, and the second arises from the interactions of the γ charges with the original. Both can be derived from Eq. (A13) which leads to the result

7R.

Tehver, F. Ancilotto, F. Toigo, J. Koplik, and J. R. Banavar, Phys. Rev. E 59, R1335 (1999). 8G. Stell, J. Phys. Condens. Matter 8, 9329 (1996). 9Y. C. Kim and M. E. Fisher, Phys. Rev. Lett. 92, 185703 (2004). 10P. J. Camp and G. N. Patey, Phys. Rev. A 60, 1063 (1999). 11D. Coslovich, J.-P. Hansen, and G. Kahl, J. Chem. Phys. 134, 244514 (2011). 12D. Coslovich, J.-P. Hansen, and G. Kahl, Soft Matter 7, 1690 (2011). 13A. Nikoubashman, J.-P. Hansen, and G. Kahl, J. Chem. Phys. 137, 094905 (2012). 14A. Nikoubashman, J.-P. Hansen, and G. Kahl, J. Chem. Phys. 138, 109901 (2013). 15P. B. Warren and A. J. Masters, J. Chem. Phys. 138, 074901 (2013). 16P. B. Warren and A. Vlasov, J. Chem. Phys. 140, 084904 (2014). 17M. E. Fisher and D. Ruelle, J. Math. Phys. 7, 260 (1966). 18J.-M. Caillol and D. Levesque, J. Chem. Phys. 140, 214505 (2014). 19G. Rickayzen and D. M. Heyes, Mol. Phys. 109, 1373 (2011). 20A. Lang, C. N. Likos, M. Watzlawek, and H. Löwen, J. Phys.: Condens. Matter 12, 5087 (2000). 21D. M. Heyes and G. Rickayzen, Mol. Phys. 112, 1398 (2014). 22D. M. Heyes and G. Rickayzen, J. Chem. Phys. 140, 024506 (2014). 23J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 3rd ed. (Academic Press, Amsterdam, 2006). 24J. F. Springer, M. A. Pokrant, and F. A. Stevens, Jr., J. Chem. Phys. 58, 4863 (1973). 25D. M. Heyes, Phys. Chem. Liq. 20, 115 (1989). 26D. M. Heyes, J. Chem. Phys. 132, 064504 (2010). 27G. Rickayzen and D. M. Heyes, J. Phys.: Condens. Matter 19, 416101 (2007). 28D. M. Heyes and A. C. Bra´ nka, J. Chem. Phys. 138, 034504 (2013).


Asymmetric charge distributions in planar bilayer systems.

Vapor Interface Revisited by Path-Integral and Ab Initio Molecular Dynamics Simulations.

Hydrogen Bonding and Dielectric Spectra of Ethylene Glycol-Water Mixtures from Molecular Dynamics Simulations.

Mutual and Self-Diffusivities in Binary Mixtures of [EMIM][B(CN)4] with Dissolved Gases by Using Dynamic Light Scattering and Molecular Dynamics Simulations.

Stokes shift dynamics in (non-dipolar ionic liquid + dipolar solvent) binary mixtures: a semi-molecular theory.

Molecular dynamics simulations of complex mixtures aimed at the preparation of naproxen-imprinted xerogels.

alkali salt mixtures interface.

Hydration of porphyrin and Mg-porphyrin: ab initio quantum mechanical charge field molecular dynamics simulations.

Thermal conductivity of molten salt mixtures: Theoretical model supported by equilibrium molecular dynamics simulations.

Charge templates in aromatic plus ionic liquid systems revisited: NMR experiments and molecular dynamics simulations.

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations.

Properties of lipid electropores I: Molecular dynamics simulations of stabilized pores by constant charge imbalance.

Molecular dynamics simulations in biology.

Floating orbital molecular dynamics simulations.

Charge-Neutral Constant pH Molecular Dynamics Simulations Using a Parsimonious Proton Buffer.

A new approach to calculate charge carrier transport mobility in organic molecular crystals from imaginary time path integral simulations.

Impact of distributions and mixtures on the charge transfer properties of graphene nanoflakes.

Molecular dynamics simulations: advances and applications.

Fluid-solid coexistence from two-phase simulations: binary colloidal mixtures and square well systems.

Microsecond molecular dynamics simulations of lipid mixing.

Molecular dynamics simulations of classical stopping power.

Compressibility, thermal expansion coefficient and heat capacity of CH4 and CO2 hydrate mixtures using molecular dynamics simulations.

Molecular dynamics simulations of large macromolecular complexes.

Accelerated molecular dynamics simulations of protein folding.