The ionic atmosphere around a nucleic acid regulates its stability in aqueous salt solutions. One major source of complexity in biological activities involving nucleic acids arises from the strong influence of the surrounding ions and water molecules on their structural and thermodynamic properties. Here, we implement a classical density functional theory for cylindrical polyelectrolytes embedded in aqueous electrolytes containing explicit (neutral hard sphere) water molecules at experimental solvent concentrations. Our approach allows us to include ion correlations as well as solvent and ion excluded volume effects for studying the structural and thermodynamic properties of highly charged cylindrical polyelectrolytes. Several models of size and charge asymmetric mixtures of aqueous electrolytes at physiological concentrations are studied. Our results are in good agreement with Monte Carlo simulations. Our numerical calculations display significant differences in the ion density profiles for the different aqueous electrolyte models studied. However, similar results regarding the excess number of ions adsorbed to the B-DNA molecule are predicted by our theoretical approach for different aqueous electrolyte models. These findings suggest that ion counting experimental data should not be used alone to validate the performance of aqueous DNA-electrolyte models.

Nucleic acids (e.g., DNA and RNA) are essential components of living cells that allow organisms to store and transfer information to direct the synthesis of new proteins. These highly charged polyelectrolytes carry a −1e charge per nucleotide in its backbone. Water molecules and small ions screen the strong repulsive Coulombic interactions between the phosphate groups of nucleic acids. Experimental and theoretical data suggest the existence of a steep local concentration of counterions (i.e., cations) and the depletion of co-ions (i.e., anions) near the surface of nucleic acids due to these large Coulombic interactions.1 The ionic atmosphere around nucleic acids, which is essential for its proper biological functioning, is the main theme of the present work.

While performing their biological activities nucleic acids can bend, fold, denature, condense, and additionally they can also bind to other biomolecules, such as drugs, peptides, and proteins.2 During such events the charge densities of nucleic acids change and the surrounding ions and water molecules are redistributed. For instance, in order for a protein to bind DNA or RNA it displaces some surrounding waters and ions.3,4 Theory tells us that changes in the number of associated ions to nucleic acids is related to the derivative of the related free energy with respect to salt concentration.1 Thus, a better description of the ionic atmosphere around nucleic acids is important for explaining numerous salt dependent processes involving nucleic acids, such as melting, folding, packaging, and binding.

Although standard experimental techniques, such as X-ray crystallography can reveal ordered water molecules and specific ion binding sites around nucleic acids, they cannot provide a full description of the ionic atmosphere.5 Moreover, the identification of cation binding sites can be ambiguous since the electron density pattern of ions is similar to that of water molecules. Thus, ion counting experimental techniques such as anomalous small-angle X-ray scattering (ASAXS)5,6 and buffer equilibration and atomic emission spectroscopy (BE-AES)7 have been developed to probe the ionic atmosphere and measure the number of ions (both cations and anions) around nucleic acids. In parallel to the development of ion counting experimental techniques, improved treatment of long-range electrostatic interactions, new ion parametrizations, and enhanced sampling techniques now allow a detailed characterization of the ion distribution around nucleic acids with all-atom molecular dynamics (MD) simulations in the nanosecond and microsecond time scales.8–11 However, such explicit solvent simulations of nucleic acids in dilute electrolytes can be computationally very expensive.

Experimental and theoretical studies on the monovalent salt dependent thermodynamic and kinetic properties of nucleic acids support the use of the Poisson-Boltzmann equation (PBE), which is an implicit solvent (water) model that assumes a mean-field description of the ion distribution.12–15 However, these monovalent salt properties are only correctly predicted by the nonlinear PBE (NLBPE), not the linear PBE.16 Additionally, NLPBE and its size modified forms predict that like-charged polyelectrolytes always repel each other at all salt conditions. In fact, in the PBE approach the discrete nature of the solvent and ion-ion correlations are completely ignored. When these ingredients are taken into account properly, it is possible to provide a more realistic description of the structural and thermodynamic properties of highly charged polyelectrolytes as has been shown in previous studies of strongly correlated colloidal systems.17–28 

More sophisticated approaches such as classical density functional theory (DFT)29–31 and integral equation theory (IET)32,33 of liquid solutions have been proposed to address the limitations of continuum macroscopic electrostatics in describing the ionic atmosphere around DNA. These theories provide a rigorous statistical mechanics approach to incorporate the microscopic nature of the solute and solvent in the description of solvation phenomena.34 Classical DFT and integral equation theories like three-dimensional reference interaction site model/hypernetted chain closure (3D-RISM/HNC)35 often succeed in reproducing the main features of solvation of the DNA duplexes in dense liquids in which the correlation and excluded volume effects become important.32,33 Studies based on these theories reveal charge inversion for divalent electrolyte mixtures at high concentrations for polyelectrolytes. The predictions of HNC36 and DFT in combination with Mean Spherical Approximation (MSA)37,38 for the integrated charge (IC) and the mean electrostatic potential (MEP) within the restricted primitive model electrolyte (RPM) show agreement with the corresponding MD/Monte Carlo (MC) simulation results and disagreement with PB theory.33,39,40 This shows the importance of ion-ion correlations and excluded volume effect for capturing the charge inversion around polyelectrolytes.

However, 3D-RISM predictions usually depend on approximate closures such as HNC, MSA, Kovalenko-Hirata (KH)41 to calculate density distribution profiles and involves time-consuming three-dimensional fast Fourier transform (3DFFT) calculations. For instance, IET has a more pronounced deviation from MC simulations than PB and DFT predictions for anions in monovalent electrolytes at low concentrations around DNA.33 However, for divalent mixtures, IET reproduces the ionic density profiles better than PB but less accurate than DFT when compared with the corresponding MC simulations.

To provide a more detailed and accurate characterization of the water and ionic atmosphere surrounding a highly charged cylindrical B-DNA model, we use a classical DFT, that extends the capabilities of classical mean-field theories. By using a combination between the MSA42,43 and the White Bear version II of the fundamental measure theory (FMTWBII),44 we are able to properly describe the interplay between ion-ion correlation and excluded volume effects. In addition to these effects, the solvent crowding effects are also taken into account by using an explicit solvent model, considering a neutral fluid of small hard spheres. We estimate the hard sphere excess free energy using FMTWBII. This theory provides one of the most accurate expressions for uncharged inhomogeneous multi-component hard sphere fluids that recovers a generalization of the Carnahan-Starling equation of state in the uniform-fluid limit.45 This expression plays a fundamental role in modeling highly asymmetric ion sizes and solvent excluded volume effects at experimental size and concentration. It also provides a more accurate description of the physics at solute-solvent interfaces.46 Additionally, we use the MSA for multi-component charged hard sphere fluids which accounts for the electrostatic ion correlations at low computational cost. The two aforementioned approximations are implemented along with a suitable decomposition of the excess free energy, which plays a key role in capturing the complex interplay between ion correlations and excluded volume effects. These features make our DFT approach particularly useful for studying the aqueous electrolyte dependent structural and thermodynamic properties of highly charged polyelectrolytes when a good balance between accuracy and efficiency is desired. Previously, we have successfully demonstrated the usefulness of our DFT approach for studying the structural and thermodynamic properties of electrical double layers around spherical polyelectrolyte in different aqueous polyelectrolytes.30 

In this work, we focus the analysis on the excluded volume and ion-ion correlation effects on the properties of the ionic atmosphere around a model B-DNA. We apply the proposed computational model to study an infinitely long, and homogeneously charged cylinder, having the charge density of B-DNA shown in Fig. 1 in various ionic conditions; from pure NaCl, to electrolyte mixtures containing multivalent cations (i.e., NaCl and MgCl2, NaCl and AlCl3). First, this approach is validated against Monte Carlo simulations. Then, we calculate the ion density distributions in order to examine the role of asymmetric ion size and solvent excluded volume effects on the ionic atmosphere around a model B-DNA. Additionally, we calculate the IC and MEP to study charge inversion. Also, the theoretical predictions of the number of excess ions around B-DNA in electrolyte mixtures containing NaCl and MgCl2 are compared against experimental ion counting data7 using several models for ions and solvent molecules. Finally, we study the ionic potential of mean force by comparing the normalized MEP energy with hard sphere and electrostatic residual correlations.

FIG. 1.

The electrostatic potential (in 1/(βe), where

$\beta =\frac{1}{k_{B}T}$
β=1kBT⁠) of a 24 base pair ideal Watson-Crick B-DNA helix is mapped on its solvent-excluded molecular surface (a). The 6M Na+ iso-concentration contour around the same B-DNA helix is shown in green (b). The adaptive Cartesian grid-based Poisson-Boltzmann solver (CPB)79,80 is used with the following settings: the interior and exterior dielectric constant are set to 1 and 80, respectively, the solution temperature is 298.15 K and the NaCl concentration is 0.1M. An ion exclusion thickness of 2 Åis employed. A simplified B-DNA charge model, where each non-bridging phosphate oxygen is assigned a charge of −0.5e and all other atomic charges are set equal to zero, is employed to showcase sequence independent electrostatic features of the classical B-DNA.

FIG. 1.

The electrostatic potential (in 1/(βe), where

$\beta =\frac{1}{k_{B}T}$
β=1kBT⁠) of a 24 base pair ideal Watson-Crick B-DNA helix is mapped on its solvent-excluded molecular surface (a). The 6M Na+ iso-concentration contour around the same B-DNA helix is shown in green (b). The adaptive Cartesian grid-based Poisson-Boltzmann solver (CPB)79,80 is used with the following settings: the interior and exterior dielectric constant are set to 1 and 80, respectively, the solution temperature is 298.15 K and the NaCl concentration is 0.1M. An ion exclusion thickness of 2 Åis employed. A simplified B-DNA charge model, where each non-bridging phosphate oxygen is assigned a charge of −0.5e and all other atomic charges are set equal to zero, is employed to showcase sequence independent electrostatic features of the classical B-DNA.

Close modal

The system considered in this study consists of a rigid charged cylindrical solute surrounded by an electrolyte solution comprised of m ionic species. Each ion of species i is described by a hard sphere of diameter σi, charge ezi, and bulk density

$\rho _{i}^{0}$
ρi0⁠, where e is the electron charge and zi the corresponding ionic valence (see Tables I and II). The point-like charge model and Coulomb's law is used to describe the electrostatic interaction between ions in a continuum dielectric media of dielectric constant ε. The charge of the ions is placed at the center of ions such that the ion-ion pair-potential is given by

\begin{eqnarray}&&u_{ij}(r)\nonumber\\&&\;\;\,=\!\left\lbrace \begin{array}{l@{\;\;}l}\infty , & \textnormal { }r\!=\!|\mathbf {r}_{i}-\mathbf {r}_{j}|<(\sigma _{i}+\sigma _{j})/2,\\[8pt]\frac{z_{i}z_{j}e^{2}}{\epsilon r}, & \textnormal { }r\!\ge\! (\sigma _{i}+\sigma _{j})/2, \end{array}\right.\;\; i,j=1\ldots m,\nonumber\\\end{eqnarray}
uij(r)=,r=|rirj|<(σi+σj)/2,zizje2εr,r(σi+σj)/2,i,j=1...m,
(1)

where ri, rj are the positions of ionic species i and j, respectively.

Table I.

Ion size models. The diameter of water molecules is set to 2.75 Å78 for all SPMs.

 Diameter (Å)
ModelNa+Mg+2ClAl+3
Hydrated77  7.16 8.56 6.64 9.60 
Symmetric 2.75 2.75 2.75 2.75 
Shannon68  2.04 1.44 3.62 1.08 
 Diameter (Å)
ModelNa+Mg+2ClAl+3
Hydrated77  7.16 8.56 6.64 9.60 
Symmetric 2.75 2.75 2.75 2.75 
Shannon68  2.04 1.44 3.62 1.08 
Table II.

Electrolyte mixtures and water concentrations.

  Water concentration (M)
  Ion sizes: Symmetric and Shannon, hydrated
 Electrolyte mixture concentrationSPMPMSPM
a 0.025M NaCl + 0.005M MgCl2 55.56 55.56 
b 0.1M NaCl + 0.005M MgCl2 55.56 52.33 
c 1M NaCl + 0.125M MgCl2 55.56 19.86 
d 1M NaCl + 0.25M MgCl2 55.56 13.10 
e 1M NaCl + 0.125M AlCl3 55.56 16.90 
f 1M NaCl + 0.25M AlCl3 55.56 10.00 
  Water concentration (M)
  Ion sizes: Symmetric and Shannon, hydrated
 Electrolyte mixture concentrationSPMPMSPM
a 0.025M NaCl + 0.005M MgCl2 55.56 55.56 
b 0.1M NaCl + 0.005M MgCl2 55.56 52.33 
c 1M NaCl + 0.125M MgCl2 55.56 19.86 
d 1M NaCl + 0.25M MgCl2 55.56 13.10 
e 1M NaCl + 0.125M AlCl3 55.56 16.90 
f 1M NaCl + 0.25M AlCl3 55.56 10.00 

In order to study the ionic and water structure around a B-DNA molecule in ionic solutions, we use four aqueous electrolyte models to represent ions and water molecules as shown in Fig. 2. In the first one, called primitive model (PM), the ionic solution is approximated by a structureless continuum medium characterized by a dielectric constant ε = 78.35, and the ions are represented as hard spheres with point charges located at their centers (see Fig. 2(a)). As shown in Figs. 2(b)–2(d), when the water molecules are represented by neutral hard spheres as well, the resulting model is known as the solvent primitive model (SPM) or “civilized” primitive model.47,48 Even though SPM is not able to consider solvent polarization and other solvent electrostatic correlations, it includes solvent excluded volume effects. It has been shown that SPM provides a more realistic description of aqueous electrolyte properties when compared with PM.49 As shown in Figs. 2(b)–2(d), several approximations can be used to model the size of the ionic species in a SPM electrolyte. The ionic species can be considered unhydrated (Fig. 2(b)), with the same size as the solvent hard spheres (Fig. 2(c)), or hydrated ions (Fig. 2(d)).

FIG. 2.

Schematic representation of the different aqueous electrolyte models employed in this work. (a) Primitive model (continuum solvent) with hydrated ions; (b) solvent primitive model (SPM) with unhydrated ions; (c) SPM with equally sized ions and water hard spheres; and (d) SPM with hydrated ions and water hard spheres. In all the SPM models, the size of the water hard spheres is the same.

FIG. 2.

Schematic representation of the different aqueous electrolyte models employed in this work. (a) Primitive model (continuum solvent) with hydrated ions; (b) solvent primitive model (SPM) with unhydrated ions; (c) SPM with equally sized ions and water hard spheres; and (d) SPM with hydrated ions and water hard spheres. In all the SPM models, the size of the water hard spheres is the same.

Close modal

The B-DNA molecule is modeled as an infinite hard cylinder of radius R = 8 Å and uniform axial charge density λ = −0.94 nC/m. It has been shown that this approximate cylindrical model is able to provide a good description of salt-dependent properties of B-DNA.47,48 However, due to lack of structural detail this cylindrical B-DNA model cannot account for sequence specific properties.1,16 It is also unable to capture the biomolecule polarization effects.

The presence of the polyelectrolyte in the bulk electrolyte is represented by an external potential ui(r) acting on each ionic species i. This external potential generates inhomogeneous ion density profiles ρi(r) representing the probability of finding ions around the polyelectrolyte. In this work, ρi(r) is calculated from pair correlation functions and functional expansions around the bulk ionic densities

$\rho _{1}^{0},\rho _{2}^{0},\ldots ,\rho _{m}^{0}(\equiv \lbrace \rho _{j}^{0}\rbrace )$
ρ10,ρ20,...,ρm0({ρj0}) using the classical density functional theory.29–31 A detailed description of the theoretical approach is given in the Appendix.

Under this framework, the explicit expression for the ion density profiles ρi(r) depends on the way in which the excess Helmholtz free energy of the system Fex[{ρj}] is decomposed and subsequently approximated. We use the following decomposition of the excess free energy Fex:34 

\begin{equation}F^{ex}[\lbrace \rho _{j}\rbrace ]=F_{Coul}^{ex}[\lbrace \rho _{j}\rbrace ]+F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]+F_{res}^{ex}[\lbrace \rho _{j}\rbrace ],\end{equation}
Fex[{ρj}]=FCoulex[{ρj}]+Fhsex[{ρj}]+Fresex[{ρj}],
(2)

where the first term,

\begin{equation}F_{Coul}^{ex}[\lbrace \rho _{j}\rbrace ]=\frac{1}{2}\sum _{i}\sum _{k}q_{i}q_{k}\int\!\!\!\int d\mathbf {r}d\mathbf {r}^{\prime }\frac{\rho _{i}(r)\rho _{k}(r^{\prime })}{|\mathbf {r}-\mathbf {r}^{\prime }|}\end{equation}
FCoulex[{ρj}]=12ikqiqkdrdrρi(r)ρk(r)|rr|
(3)

represents the pure Coulomb interactions, and is also present in mean-field theories such as the nonlinear Poisson-Boltzmann equation. The second and third terms,

$F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]$
Fhsex[{ρj}] and
$F_{res}^{ex}[\lbrace \rho _{j}\rbrace ]$
Fresex[{ρj}]
, represent the contributions corresponding to the pure hard sphere and residual Coulomb interactions, respectively. These terms provide corrections to the nonlinear Poisson-Boltzmann (NLPB) predictions, and capture charge inversion, charge-asymmetry dependent effects, and other important phenomena characterizing highly charged interacting systems.50–52 The excess free energy decomposition in Eq. (2) yields the following fundamental expression for the ion density distributions:

\begin{align}\rho _{i}(\mathbf {r}) & =\left\lbrace \begin{array}{l@{\quad}l}\rho _{i}^{0}exp\lbrace -\beta q_{i}\psi (\mathbf {r};\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )\rbrace , & \textnormal { }r>R+\sigma _{i}/2,\\[8pt]0, & \textnormal { }r\le R+\sigma _{i}/2, \end{array}\right.\\& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad i=1\ldots m,\nonumber\end{align}
ρi(r)=ρi0exp{βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj})},r>R+σi/2,0,rR+σi/2,i=1...m,
(4)

where ψ(r; {ρj}) is the solution of NLPB equation representing the MEP due to the polyelectrolyte plus the ion distribution,

$\Delta c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )\equiv c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )-\tilde{c}_{i}^{(1)}(\lbrace \rho _{j}^{0}\rbrace )$
Δci(1)(r;{ρj})ci(1)(r;{ρj})c̃i(1)({ρj0})⁠, and
$c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)hs(r;{ρj})
and
$c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)res(r;{ρj})
are the hard sphere and residual electrostatic one particle direct correlation functions, respectively. They are defined as the first functional derivatives of the excess free energies
$F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]$
Fhsex[{ρj}]
and
$F_{res}^{ex}[\lbrace \rho _{j}\rbrace ]$
Fresex[{ρj}]
, respectively. Note that
$\tilde{c}_{i}^{(1)}(\lbrace \rho _{j}^{0}\rbrace )$
c̃i(1)({ρj0})
represents one particle direct correlation function for uniform fluids. Another point to note is that each of the terms in the argument of the exponential in Eq. (4) are actually dependent on each other, because all of them are functionals of the ion density profiles in a self-consistent manner. In particular, this means that ion-ion correlation effects cannot be completely decoupled from excluded volume effects in the current approach.

In our theoretical formulation, the residual electrical one particle direct correlation function

$c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)res(r;{ρj}) is calculated perturbatively around uniform fluids (see Subsection A 3 of the Appendix). To reduce the computational cost, only the first order term in such expansion in powers of
$(\rho _{i}(r)-\rho _{i}^{0})$
(ρi(r)ρi0)
is considered here.34 To calculate the first order approximation, we implement the MSA for charged multi-component hard sphere fluids, since it provides a compact analytic expression for an accurate and efficient evaluation of residual two particle direct correlation functions.42,53

To evaluate

$c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)hs(r;{ρj})⁠, we use the expression for the free energy density of homogeneous fluids introduced by Hansen-Goos and Roth44 in the FMTWBII that yields a new generalization of the Carnahan-Starling equation of state for additive mixtures of hard spheres,

\begin{equation}\beta P_{GR}=\frac{\xi _{0}}{1-\xi _{3}}+\frac{\xi _{1}\xi _{2}\left[1+\frac{1}{3}\xi _{3}^{2}\right]}{(1-\xi _{3})^{2}}+\frac{\xi _{2}^{3}\left[1-\frac{2}{3}\xi _{3}+\frac{1}{3}\xi _{3}^{2}\right]}{12\pi (1-\xi _{3})^{3}},\end{equation}
βPGR=ξ01ξ3+ξ1ξ21+13ξ32(1ξ3)2+ξ23123ξ3+13ξ3212π(1ξ3)3,
(5)

where P is the pressure and the scaled-particle variables ξa are defined as follows:

\begin{eqnarray*}&&\xi _{a}=\sum _{i=1}^{m}\rho _{i}^{0}\mathcal {R}_{i}^{a},\qquad \mathcal {R}_{i}^{0}=1,\quad ,\mathcal {R}_{i}^{1}=\sigma _{i}/2,\\\&& \quad \mathcal {R}_{i}^{2}=\pi \sigma _{i}^{2},\quad \mathcal {R}_{i}^{3}=\pi \sigma _{i}^{3}/6.\end{eqnarray*}
ξa=i=1mρi0Ria,Ri0=1,,Ri1=σi/2,Ri2=πσi2,Ri3=πσi3/6.

The superiority of the generalization of the Carnahan-Starling equation of state, compared with former equations of state, produces differences mainly at high densities and high ion size asymmetry. This becomes clear when re-writing Eq. (5) in terms of the Percus Yevick (PY),54 the well-known Boublik-Mansoori-Carnahan-Starling-Leland (BMCSL),55 and the so-called extended Carnahan-Starling (eCS) equations of state as follows:

\begin{equation}\beta P_{GR}=\beta P_{PY}^{c}-\Delta \beta P_{BMCSL}+\Delta \beta P_{eCS}+\Delta \beta P_{GR}.\end{equation}
βPGR=βPPYcΔβPBMCSL+ΔβPeCS+ΔβPGR.
(6)

The first term on the right-hand side of Eq. (6) represents the compressibility pressure predicted by PY approximation, which deviates significantly from molecular dynamics results for uniform fluids mainly at moderate to high ion concentrations.44 The second and third terms represent the first and second correction in power of the total packing fraction ξ3 to the PY prediction obtained empirically by BMCSL (

$\Delta \beta P_{BMCSL}=\beta P_{BMCSL}-\beta P_{PY}^{c}=\mathfrak {\mathcal {O}}(\xi _{3})$
ΔβPBMCSL=βPBMCSLβPPYc=O(ξ3)⁠) and eCS (
$\Delta \beta P_{eCS}=\beta P_{eCS}-\beta P_{PY}^{c}=\mathfrak {\mathcal {O}}(\xi _{3}^{2})$
ΔβPeCS=βPeCSβPPYc=O(ξ32)
), respectively. The last term in Eq. (6) represents the Hansen-Goos and Roth residual contribution,

\begin{equation}\Delta \beta P_{GR}=-\frac{\xi _{1}\xi _{2}\xi _{3}^{3}}{3(1-\xi _{3})^{3}}=\mathfrak {\mathcal {O}}\big(\xi _{3}^{3}\big)\end{equation}
ΔβPGR=ξ1ξ2ξ333(1ξ3)3=Oξ33
(7)

which provides a higher order correction in the prediction of the equation of state. This correction comes from the contribution to the excess free energy density

$\Delta \Phi _{hs}^{ex}\break=\int d\xi _{3}\Delta \beta P_{GR}$
ΔΦhsex=dξ3ΔβPGR⁠,44 which produces a correction
$\Delta F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]$
ΔFhsex[{ρj}]
to the excess free energy (Eq. (A12) in Subsection A 2 of the Appendix) as well as
$\Delta c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj})
to the direct correlation function (Eq. (A15) in Subsection A 2 of the Appendix). By construction, the residual term in Eq. (7) also restores the thermodynamic consistency (up to third order in the parameter ξ3) missed by empirically based BMCSL and eCS approximations. As a result, Hansen-Goos and Roth's44 approach is appropriate for modeling water excluded volume effects at experimental size and concentration and high size ionic asymmetry. It also plays a relevant role in the description of the physics at solute-liquid interfaces.46 

In order to test the accuracy of our theoretical results for the implicit solvent model, we perform Monte Carlo simulations of the B-DNA model immersed in a 1M 1:1 PM electrolyte solution. In these simulations, ions are represented by hard spheres with point charges placed at their centers that interact as described in Eq. (8).

The B-DNA is modeled as an infinite hard cylinder of radius R and homogeneous surface charge density. The electric field produced by this cylinder is equivalent to that generated by an infinite line of charge of density λ, which is placed in the center of the cylinder along its symmetry axis. In order to compare simulation data and theoretical results, we place 10 equally spaced point-charges along the cylinder's symmetry axis of value qi = e/10 each 1.7 Å. The interaction between an ion i of diameter σi and a point-charge qj is given by

\begin{equation}u_{ij}^{\prime }(r^{\prime },r)=\left\lbrace \begin{array}{l@{\quad}l}\infty , & \textnormal { }r^{\prime }<R+\sigma _{i}/2,\\[8pt]\frac{ez_{i}q_{j}}{\epsilon r}, & \textnormal { }r^{\prime }\ge R+\sigma _{i}/2, \end{array}\right.\end{equation}
uij(r,r)=,r<R+σi/2,eziqjεr,rR+σi/2,
(8)

where r is the perpendicular distance of an ion i to the line of charge, r is the distance between an ion i (of charge ezi) and a point-charge qj located along the line of charge.

A cubic simulation box of length L with periodic boundary conditions is used to perform NVT Monte Carlo simulations via the Metropolis algorithm.56,57 In order to avoid border effects, the length L is several times the Debye length for each ionic mixture. In typical runs, L is varied from 15 nm to 30 nm, involving 4000-7000 particles. The total number of ions is also adjusted to obtain the desired bulk concentration far away from the cylinder surface, but always fulfilling the electroneutrality condition,

\begin{equation}\sum _{i=1}^{N_{cyl}}q_{i}+N_{c}z_{c}e+\sum _{k=1}^{m}\sum _{j=1}^{N_{k}}ez_{k}=0,\end{equation}
i=1Ncylqi+Nczce+k=1mj=1Nkezk=0,
(9)

where Ncyl is the number of discrete charges placed on the line of charge, Nc and zc are the number and valence of free counterions of the line of charge, m is the total number of ionic species, and Nk and zk are the number of ions and valence corresponding to the ionic species k, respectively.

Recently, several computational schemes have been proposed to properly handle long-ranged Coulombic interactions.58–62 In this work, we use the classical Ewald sums approach57,58 with a damping constant α = 5/L (where L is the length of the cubic simulation box) and 725 vectors in the k-space to compute the reciprocal space contribution to the total electrostatic energy (details of the implementation are discussed elsewhere63–65).

In all simulations, at least 100 000 Monte Carlo cycles are used for thermalization, and 200 000-1 000 000 Monte Carlo cycles are performed to calculate the canonical average.

Additionally, we test the accuracy of our theoretical results in the solvent primitive model using the Goel et al. Monte Carlo simulation results.66 We analyze the ion and water density distributions for electrolytes containing mixed 1:2:1 (NaCl/MgCl2) having 1M 1:1 (NaCl) and various concentrations of 2:1 (MgCl2) in order to validate the approximations introduced in our proposed computational DFT approach including excluded volume effects of small neutral hard sphere water molecules.

We calculate the normalized ion density profiles around the model B-DNA cylinder immersed in the aqueous electrolytes listed in Table II using both the FMTWBII + MSA and MC simulations. The ions are modeled with hydrated diameters listed in Table I. Fig. 3 shows that the theoretical predictions are in good agreement with the corresponding MC results for the ion profiles in the 1M 1:1 NaCl electrolyte surrounding the B-DNA. In Fig. 3, the first peak of the solid curve reproduces the accumulation of Na+ ions near the surface of the B-DNA, where the electrostatic attraction between Na+ ions and the surface charge of the B-DNA dominates the steric repulsion of accumulated Na+ ions from each other. Additionally, the curve reproduces Cl depletion by the polyelectrolyte surface due to the electrostatic repulsion of Cl ions by the surface charge of B-DNA.

FIG. 3.

Normalized ion density profiles (ρ/ρo) around B-DNA (radius of 8 Å and linear charge density −0.94 nC/m, the same B-DNA model is used in all the sections) embedded in a 1M NaCl solution as a function of radial distance (r) from the B-DNA helical axis. The Na+ and Cl diameters are set to 7.16 Å and 6.64 Å, respectively. Solid and short-dashed lines represent our theoretical predictions for Na+ and Cl ion density distributions, respectively, and symbols represent the MC simulation results. Triangles represent Na+, and circles correspond to Cl.

FIG. 3.

Normalized ion density profiles (ρ/ρo) around B-DNA (radius of 8 Å and linear charge density −0.94 nC/m, the same B-DNA model is used in all the sections) embedded in a 1M NaCl solution as a function of radial distance (r) from the B-DNA helical axis. The Na+ and Cl diameters are set to 7.16 Å and 6.64 Å, respectively. Solid and short-dashed lines represent our theoretical predictions for Na+ and Cl ion density distributions, respectively, and symbols represent the MC simulation results. Triangles represent Na+, and circles correspond to Cl.

Close modal

Fig. 4 shows the effects of divalent and trivalent ionic species in the primitive model at several ionic concentrations. At lower ionic concentrations (Figs. 4(a) and 4(b)), the normalized ionic profile near the polyelectrolyte goes to the bulk monotonically. When the concentration of 1:1 electrolyte changes from 0.025M to 0.1M, the ionic strength increases. This in turn reduces the Debye screening length, which depends inversely on the ionic strength and is used as a measure of the thickness of the ionic atmosphere. As a result, the thickness of the ionic atmosphere reduces with increased electrolyte concentration and the ion density distributions reach their bulk values at faster rates as shown in Fig. 4(b). At strong ionic strength regimes (around 1M), layering effects are observed in the ion profiles of the electrolyte mixtures (see Figs. 4(c)–4(f)). Also, the addition of MgCl2 or AlCl3 to the pure 1M NaCl solution causes a displacement of sodium ions by the multivalent ions leading to a larger concentration of Mg+2 or Al+3 than the concentration of Na+ at the surface of B-DNA, and decrease of Na+ concentration at the surface compared to pure 1M NaCl. Similar results have been observed in other works.66,67 The comparison between 1:1 NaCl/2:1 MgCl2 (see Figs. 4(c) and 4(d)) and 1:1 NaCl/3:1 AlCl3 (see Figs. 4(e) and 4(f)) electrolyte mixture of equal concentrations shows that trivalent ion species are more strongly attracted to the polyelectrolyte than divalent ion species. This can be observed in Figs. 4(d) and 4(e) as the height of the first peak of Al+3 is much taller than that of the Mg+2, and consequently, a larger adsorption of counterions effectively screens the original B-DNA surface charge and promotes the adsorption of co-ions. In Figs. 4(c)–4(f), there is a hump in the density profiles of co-ions at around 20 Å. Fig. 4(d) shows a stronger hump with the increase in MgCl2 concentration compared to the hump at a lower concentration shown in Fig. 4(c). In the electrolyte mixture containing AlCl3, the hump in the density profiles becomes even stronger as we increase the electrolyte concentrations (see Figs. 4(e) and 4(f)). The behavior of the hump in the ionic density profiles is attributable to the ion-ion correlations.66 For further validation of our theoretical predictions, we compare the normalized ion density predictions for SPM with the corresponding MC simulation results of Goel et al.66 shown in Figs. 5 and 6. A good agreement is found between the MC simulation data and the DFT predictions.

FIG. 4.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in three component asymmetric electrolyte solutions. The hydrated ion sizes are given in Table I and the water molecules are not considered explicitly (PM). Plots (a)-(f) correspond to the electrolyte mixture concentrations given in Table II from top to bottom, respectively. Dashed, solid, and short-dashed lines correspond to our theoretical predictions for Mg+2 or Al+3, Na+, and Cl ion density distributions, respectively, whereas squares, triangles, and circles represent the corresponding MC simulation results, respectively.

FIG. 4.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in three component asymmetric electrolyte solutions. The hydrated ion sizes are given in Table I and the water molecules are not considered explicitly (PM). Plots (a)-(f) correspond to the electrolyte mixture concentrations given in Table II from top to bottom, respectively. Dashed, solid, and short-dashed lines correspond to our theoretical predictions for Mg+2 or Al+3, Na+, and Cl ion density distributions, respectively, whereas squares, triangles, and circles represent the corresponding MC simulation results, respectively.

Close modal
FIG. 6.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in three component symmetric electrolyte solutions. All ion sizes are equal to 4.0 Å. The water molecules are explicitly considered (SPM) at 10.3M concentration and size set to 4.0 Å. The plots (a) correspond to 1M NaCl + 0.125M MgCl2 + 10.3M H2O electrolyte mixture, whereas (b) corresponds to 1M NaCl + 0.25M MgCl2 + 10.3M H2O electrolyte mixture. Dashed, solid, dotted-dashed, and short-dashed lines correspond to our theoretical predictions for Mg+2 , Na+, solvent, and Cl ion density distributions, respectively, whereas solid squares, triangles, empty and solid circles represent the corresponding MC simulation results, respectively.66 

FIG. 6.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in three component symmetric electrolyte solutions. All ion sizes are equal to 4.0 Å. The water molecules are explicitly considered (SPM) at 10.3M concentration and size set to 4.0 Å. The plots (a) correspond to 1M NaCl + 0.125M MgCl2 + 10.3M H2O electrolyte mixture, whereas (b) corresponds to 1M NaCl + 0.25M MgCl2 + 10.3M H2O electrolyte mixture. Dashed, solid, dotted-dashed, and short-dashed lines correspond to our theoretical predictions for Mg+2 , Na+, solvent, and Cl ion density distributions, respectively, whereas solid squares, triangles, empty and solid circles represent the corresponding MC simulation results, respectively.66 

Close modal
FIG. 5.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in a 1M NaCl solution. The water molecules are explicitly considered (SPM) at 10.3M concentration. The Na+, Cl, and water molecule radii are set to 4.0 Å. Solid, dotted-dashed, and short-dashed lines represent our theoretical predictions for Na+, solvent, and Cl ion density distributions, respectively, whereas solid triangles, empty and solid circles represent the corresponding MC simulation results, respectively.66 

FIG. 5.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in a 1M NaCl solution. The water molecules are explicitly considered (SPM) at 10.3M concentration. The Na+, Cl, and water molecule radii are set to 4.0 Å. Solid, dotted-dashed, and short-dashed lines represent our theoretical predictions for Na+, solvent, and Cl ion density distributions, respectively, whereas solid triangles, empty and solid circles represent the corresponding MC simulation results, respectively.66 

Close modal

1. Solvent excluded volume effects on the ion density distributions

To understand the solvent excluded volume effects on the ionic atmosphere surrounding B-DNA, the ion profiles of SPM aqueous electrolyte are compared with those of the corresponding PM electrolyte. We also investigate the influence of water molar concentration and molecule size on the ion density distributions by fixing the partial packing fraction of water molecules around 0.4 and changing the water molar concentration accordingly for each electrolyte listed in Table II. It should be noted that even for PM, the solvent excluded volume effect is partially incorporated because of the use of large hydrated ionic diameters. At low NaCl concentration, the PM shows an exponential and monotonic behavior of the ion density profiles (see red lines in Figs. 7(a) and 7(b)). This behavior is similar to that of PB theory. By contrast, the explicit presence of solvent introduces layering in the ion density distributions (see black lines in Figs. 7(a)–7(f)). For SPM, the increased accumulation of ions at the contact point is clearly attributed to the solvent excluded volume effects, which push the ions toward the biomolecular surface. Additionally, at higher solvent concentrations (Figs. 7(a) and 7(b)), the oscillatory behavior within the ionic atmosphere is more pronounced, causing stronger layering of the ions.

FIG. 7.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in different solvent models and electrolyte mixtures. The electrolyte mixture concentrations for (a)-(f) plots are given in Table II from top to bottom, respectively. Dashed, solid, and short-dashed lines correspond to our theoretical predictions for Mg+2 or Al+3, Na+, and Cl ion density distributions with ion bulk concentrations described in Fig. 4. Red lines correspond to PM (without explicit water molecules) and black lines represent SPM (with explicit water molecules and hydrated ion sizes) ion density distributions at six different water concentrations from Table II.

FIG. 7.

Normalized ion density profiles (ρ/ρo) around B-DNA embedded in different solvent models and electrolyte mixtures. The electrolyte mixture concentrations for (a)-(f) plots are given in Table II from top to bottom, respectively. Dashed, solid, and short-dashed lines correspond to our theoretical predictions for Mg+2 or Al+3, Na+, and Cl ion density distributions with ion bulk concentrations described in Fig. 4. Red lines correspond to PM (without explicit water molecules) and black lines represent SPM (with explicit water molecules and hydrated ion sizes) ion density distributions at six different water concentrations from Table II.

Close modal

2. Ionic size asymmetry effects on the ion and water density distributions

To investigate the effect of asymmetrically sized ions on the ion density distribution around B-DNA, we study three electrolyte mixtures: a, d, and f listed in Table II, using the three different ionic size models listed in Table I. These three models cover the cases where the co-ions are smaller, equal, or greater than the counterions. In the first model, where hydrated ion sizes are used, counterions are bigger than co-ions. Using this model for different electrolyte mixtures requires varying water molar concentrations, whose values are listed in Table II. The second model employs symmetric ions with equal counterion and coion sizes. The third model uses Shannon effective ionic radii,68 in which the counterions are smaller than the co-ions. For the second and third models, the molar concentration of water is 55.56M.

The comparison between the symmetric and Shannon size models (see green solid and blue dashed lines in Fig. 8) shows that the ion asymmetry influences the position and the height of the ion density peak around the biopolyelectrolyte surface. In the case of Shannon ion size model, the counterions being smaller in size than the co-ion come closer to the polyelectrolyte surface, whereas in the symmetric ion model, the counterions and co-ion come to the same distance to the B-DNA surface. On the other hand, in the hydrated ion size model (see black short-dashed lines in Fig. 8) the first peak in the density distribution of the counterions and co-ion are located at further distances from the B-DNA surface, compared with those predicted by the Shannon or symmetric ion sizes models, due to the bigger ionic sizes. The local concentrations of the ions are also affected by the ion sizes, since the height of the peaks in the ion density distributions is different for each ion size model (see Fig. 8). Overall, the results presented in Fig. 8 indicate that ion sizes as well as the asymmetry in ion sizes have a strong impact on the ionic atmosphere surround the B-DNA cylinder.

FIG. 8.

Normalized ion density profiles (ρ/ρo) around B-DNA for SPM electrolytes with different ionic diameters. Green solid lines correspond to symmetric ion sizes, blue dashed lines represent Shannon's effective ion sizes, and black short-dashed lines show the hydrated ion sizes as given in Table I. The (a), (b), (c) plots represent Na+, Mg+2, and Cl for a electrolyte mixture from Table II, the (d), (e), (f) plots represent Na+, Mg+2, and Cl for d electrolyte mixture, and the (g), (h), (i) plots represent Na+, Al+3, and Cl for f electrolyte mixture.

FIG. 8.

Normalized ion density profiles (ρ/ρo) around B-DNA for SPM electrolytes with different ionic diameters. Green solid lines correspond to symmetric ion sizes, blue dashed lines represent Shannon's effective ion sizes, and black short-dashed lines show the hydrated ion sizes as given in Table I. The (a), (b), (c) plots represent Na+, Mg+2, and Cl for a electrolyte mixture from Table II, the (d), (e), (f) plots represent Na+, Mg+2, and Cl for d electrolyte mixture, and the (g), (h), (i) plots represent Na+, Al+3, and Cl for f electrolyte mixture.

Close modal

To provide a more complete picture on layering formation and characterization of an electrolyte near its surface, we also study the ionic size asymmetry and solvent excluded volume effects on the hydration shells around B-DNA. In order to study water density distribution for different aqueous electrolyte models, we use the SPM with hydrated ion sizes and the a, d, and f electrolyte mixtures discussed above. We find that oscillatory behavior of the water density distributions depends strongly on the solvent and electrolyte mixture concentrations (see Fig. 9). Specifically, in the cases of a, d, and f, the height of the contact peaks and amplitude of oscillations of the solvent particles are consistent with the solvent volume fraction and the concentration imposed to keep the total volume fraction approximately at the same value. These features cannot be captured by PM since the water density distribution is absent. This behavior for PM is represented by the red line in Fig. 9.

FIG. 9.

Normalized water density profiles (ρ/ρo) around B-DNA. The red line corresponds to the PM, and the black dashed, blue short-dashed, and green dotted-dashed lines correspond to the SPM with hydrated ion sizes for a, d, and f electrolyte mixtures, respectively, introduced in Table II.

FIG. 9.

Normalized water density profiles (ρ/ρo) around B-DNA. The red line corresponds to the PM, and the black dashed, blue short-dashed, and green dotted-dashed lines correspond to the SPM with hydrated ion sizes for a, d, and f electrolyte mixtures, respectively, introduced in Table II.

Close modal

3. Ionic size and valence effects on the integrated charge and mean electrostatic potential

We use the ion density profiles to elucidate the dependence of polyelectrolyte electrostatic properties on ionic bulk concentrations and charge valences. We calculate the IC distribution surrounding the polyelectrolyte in order to investigate the charge inversion, which occurs when the sign of polyelectrolyte's net charge changes in aqueous salt mixture due to the accumulation of ions near the surface of a polyelectrolyte. The expression used to calculate the IC is given by

\begin{equation}{Q(r)=\lambda +2\pi \int _{R}^{r}\sum _{i}z_{i}e\rho _{i}(r^{\prime })r^{\prime }dr^{\prime },\quad i=1\ldots m,}\end{equation}
Q(r)=λ+2πRrizieρi(r)rdr,i=1...m,
(10)

where zi and e represent the charge valence of ith ion and electron charge, respectively. We investigate the effects of ionic size, electrolyte concentration, and ionic charge valence on the sign of the net charge of B-DNA cylinder (see Fig. 10) using the same three electrolyte mixtures a, d, and f from Table II discussed in Sec. III B 2. We observe that SPM with symmetric ion sizes predicts charge inversion for trivalent ionic mixture but fails to predict it for divalent ionic mixture at the same concentration (see Figs. 10(b) and 10(c)). This clearly demonstrates the contribution of ionic charge valence in predicting the charge inversion.

FIG. 10.

Integrated charge profiles corresponding to the SPM for (a) a, (b) d, and (c) f aqueous electrolyte mixtures. Red solid and green dashed lines represent our theoretical SPM predictions with symmetric and Shannon's ion sizes with corresponding water concentrations from Table II.

FIG. 10.

Integrated charge profiles corresponding to the SPM for (a) a, (b) d, and (c) f aqueous electrolyte mixtures. Red solid and green dashed lines represent our theoretical SPM predictions with symmetric and Shannon's ion sizes with corresponding water concentrations from Table II.

Close modal

We also calculate the MEP using Eq. (A11). The MEP is closely related to the zeta potential, which provides key information about the stability (aggregation/dispersion) of the polyelectrolyte suspension. We compare calculated MEP and IC predictions (PM using hydrated ion sizes) with MC simulations and NLPB predictions (see Fig. 11) in order to capture ionic size effects on the sign of the net charge of the B-DNA. Figs. 11(a) and 11(d) show that PM results for MgCl mixture at low concentration are in good agreement with MC simulations. Additionally, Figs. 11(b) and 11(e) show that MgCl2 results are in better agreement with MC simulations compared to AlCl3 results (see Figs. 11(c) and 11(f)). We only observe charge inversion in MEP and IC for divalent and trivalent mixtures at high concentrations. However, the NLPB does not predict the charge inversion for all studied concentrations. These observations expose the role of electrolyte concentration, ionic size, and valence effects in the charge inversion phenomenon. In particular, for large ion sizes and multivalent ions, the MEP and IC predictions reveal charge inversion at higher electrolyte concentrations.33 

FIG. 11.

Plots (a), (b), and (c) correspond to integrated charge, whereas (d), (e), and (f) represent mean electrostatic potential. We use a, d, and f aqueous electrolyte mixtures in plots (a), (b), (c) and (d), (e), (f), respectively. Red short-dashed lines represent our NLPB predictions. Black solid lines represent our PM predictions with hydrated ionic sizes. Symbols correspond to the MC simulation results.

FIG. 11.

Plots (a), (b), and (c) correspond to integrated charge, whereas (d), (e), and (f) represent mean electrostatic potential. We use a, d, and f aqueous electrolyte mixtures in plots (a), (b), (c) and (d), (e), (f), respectively. Red short-dashed lines represent our NLPB predictions. Black solid lines represent our PM predictions with hydrated ionic sizes. Symbols correspond to the MC simulation results.

Close modal

4. Number of excess ions around B-DNA: A comparison between theory and experiment

The validation of our DFT formalism against Monte Carlo simulations demonstrates that the approximations introduced in our computational model are adequate for the description of the ionic atmosphere surrounding cylindrical and highly charged polyelectrolytes, such as the ideal B-DNA model used here. The next step entails the comparison of the theoretical predictions against experimental data. This strategy is useful for verifying which physical approximations and model parameters better describe the experimental observations.7 

In an effort to assess the aqueous electrolyte models used in this work, we compare the number of excess monovalent and multivalent ions (Na+, Mg+2, and Cl) around B-DNA against experimental data. In particular, the number of excess ions of type i, Ni, can be calculated as follows:

\begin{equation}N_{i}(r)=\int _{R}^{r}\lbrace \rho _{i}(r^{\prime })-\rho _{i}(r_{max})\rbrace 2\pi hr^{\prime }dr^{\prime },\quad i=1\ldots m,\end{equation}
Ni(r)=Rr{ρi(r)ρi(rmax)}2πhrdr,i=1...m,
(11)

where [ρi(r) − ρi(rmax)] represents the excess ion concentration, r is the separation distance from the center of B-DNA, rmax is the maximum radial length (cutoff) at which ρi(rmax) is essentially indistinguishable from the bulk concentration, and h is the length of B-DNA.

Using Eq. (11) we compare the experimental data with the theoretical results obtained with both PM (Fig. 12(a)) and SPM (Figs. 12(b)–12(d)) for a 24 base pair B-DNA (e.g., h = 81.6 Å) immersed in three component electrolyte mixtures consisting of Na+, Mg+2, and Cl at different electrolyte mixture concentrations.

FIG. 12.

Number of excess ions around a B-DNA duplex as a function of the Na+ concentration. Different concentrations of NaCl are added into a background solution of 5mM Mg+2. Circles represent our theoretical predictions and triangular symbols represent experimental ion counting data.7 Red dashed lines correspond to Mg+2, blue solid lines correspond to Na+, and green short-dashed lines correspond to Cl. The plots represent the comparison between our theoretical predictions against experimental data for (a) PM with hydrated ion sizes, (b) SPM with Shannon's ion sizes and (c) same ion sizes, and (d) SPM with hydrated ion sizes.

FIG. 12.

Number of excess ions around a B-DNA duplex as a function of the Na+ concentration. Different concentrations of NaCl are added into a background solution of 5mM Mg+2. Circles represent our theoretical predictions and triangular symbols represent experimental ion counting data.7 Red dashed lines correspond to Mg+2, blue solid lines correspond to Na+, and green short-dashed lines correspond to Cl. The plots represent the comparison between our theoretical predictions against experimental data for (a) PM with hydrated ion sizes, (b) SPM with Shannon's ion sizes and (c) same ion sizes, and (d) SPM with hydrated ion sizes.

Close modal

At low Na+ concentration, all aqueous electrolyte models reproduce experimental results, which indicates that the ion size asymmetry and excluded volume effects are not significant in this regime. Figs. 12(b) and 12(c) show that at high Na+ concentrations, the SPM with Shannon's and symmetric ion sizes overestimate the experimental data for Na+ ions by 3-4 excess ion counts. In contrast, the PM and the SPM with hydrated ion sizes (Figs. 12(a) and 12(d)) underestimate the number of excess Na+ and Cl ions by 6-7 counts. Overall, our results show small changes in the number of excess ions predicted by different ion size models in comparison to the experimental data. However, it is possible that other physical modeling parameters such as detailed surface charge distribution on the B-DNA and dielectric discontinuity at the boundary between the B-DNA and aqueous electrolyte can change the number of excess ions by amounts similar to changes in ion size.69 

5. Electrostatic versus steric contributions on the ion density profiles

In Secs. III B 1–III B 4, we focus the analysis on the comparison between the ion density distributions, electrostatic properties, and experimental excess number of ion results predicted by different aqueous electrolyte models and electrolyte mixtures without identifying the driving force governing their behaviors. In this section, we study the ionic potential of mean force (PMF) per unit of thermal energy

$\triangle E_{i} =-\beta q_{i}\psi (\mathbf {r};\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
Ei=βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj})⁠, which appears in the expression for the ion density distributions (4), and whose negative gradient gives the corresponding average force acting on the ion species i due to the other ions and B-DNA. We compare the normalized MEP energy (−βqiψ(r; {ρj})) with hard sphere (
$\Delta c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj})
) and electrostatic residual (
$\Delta c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
Δci(1)res(r;{ρj})
) correlations. We analyze electrolytes at low and high NaCl concentrations using PM with hydrated ion sizes (see Fig. 13) and SPM with Shannon's ion sizes and experimental water concentration (see Fig. 14). We have chosen to study these two models because they have shown the major differences in the predictions of the number of excess ions around the B-DNA model at high electrolyte concentration. The results shown in Figs. 13 and 14 provide insight into the effects of the water crowding, hard sphere, and electrostatic residual correlation contributions on the PMF. Figs. 13(a), 13(b), 13(d), and 13(e) show a strong competition between hydrated hard sphere (green short-dashed lines) and electrostatic residual (blue dashed lines) correlation contributions for counterions using PM due to the opposite sign of each contribution. Figs. 13(c) and 13(f) show that both contributions have the same sign for co-ions. Additionally, ion hard sphere, electrostatic residual correlation contributions, and electrostatic potential energy contribution (red dotted-dashed lines) enhance the PMF (black solid lines). In Fig. 13(e), the electrostatic residual correlation provides the main contribution to the ionic PMF for divalent counterions at high electrolyte mixture concentration. In all the cases presented in Fig. 13, the magnitude of the electrostatic residual correlation is larger than the magnitude of the hard sphere contribution to the ionic potential of mean force.

FIG. 13.

Hard sphere contribution

$\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj}) (green short-dashed), electrostatic correlation contribution
$\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)res(r;{ρj})
(blue dashed), normalized electrostatic potential energy contribution (−βqiψ(r; {ρj})) (red dotted-dashed), and the PMF (sum of these three contributions
$(-\beta q_{i}\psi (r;\lbrace \rho _{j}\rbrace )\break +\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace ))$
(βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj}))
(black solid) (ρ/ρo) around B-DNA for the PM. (a), (b), (c) plots correspond to Na+, Mg+2, and Cl ions in 0.005M NaCl + 0.005M MgCl2 electrolyte mixture, respectively. (d), (e), (f) plots correspond to Na+, Mg+2, and Cl ions in 0.5M NaCl + 0.005M MgCl2 electrolyte mixture, respectively.

FIG. 13.

Hard sphere contribution

$\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj}) (green short-dashed), electrostatic correlation contribution
$\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)res(r;{ρj})
(blue dashed), normalized electrostatic potential energy contribution (−βqiψ(r; {ρj})) (red dotted-dashed), and the PMF (sum of these three contributions
$(-\beta q_{i}\psi (r;\lbrace \rho _{j}\rbrace )\break +\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace ))$
(βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj}))
(black solid) (ρ/ρo) around B-DNA for the PM. (a), (b), (c) plots correspond to Na+, Mg+2, and Cl ions in 0.005M NaCl + 0.005M MgCl2 electrolyte mixture, respectively. (d), (e), (f) plots correspond to Na+, Mg+2, and Cl ions in 0.5M NaCl + 0.005M MgCl2 electrolyte mixture, respectively.

Close modal
FIG. 14.

Hard sphere contribution

$\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj}) (green short-dashed), electrostatic correlation contribution
$\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)res(r;{ρj})
(blue dashed), normalized electrostatic potential energy contribution (−βqiψ(r; {ρj})) (red dotted-dashed), and the PMF (sum of these three contributions
$(-\beta q_{i}\psi (r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
(βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj})
) (black solid) for the SPM at 55.56M solvent concentration. Plots (a)-(f) represent the same electrolyte mixtures described in Fig. 13.

FIG. 14.

Hard sphere contribution

$\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)hs(r;{ρj}) (green short-dashed), electrostatic correlation contribution
$\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
Δci(1)res(r;{ρj})
(blue dashed), normalized electrostatic potential energy contribution (−βqiψ(r; {ρj})) (red dotted-dashed), and the PMF (sum of these three contributions
$(-\beta q_{i}\psi (r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)hs}(r;\lbrace \rho _{j}\rbrace )+\Delta c_{i}^{(1)res}(r;\lbrace \rho _{j}\rbrace )$
(βqiψ(r;{ρj})+Δci(1)hs(r;{ρj})+Δci(1)res(r;{ρj})
) (black solid) for the SPM at 55.56M solvent concentration. Plots (a)-(f) represent the same electrolyte mixtures described in Fig. 13.

Close modal

A different scenario is observed for explicit water molecules at experimental concentration values and Shannon's ion sizes. In Fig. 14, layering formation is observed in the PMF arising from the water crowding effect. Additionally, in all these cases the hard sphere correlation contribution is larger than the electrostatic residual correlation contribution at the contact point. Very near the B-DNA surface, Figs. 14(a), 14(b), 14(d), and 14(e) show the (ionic plus water) hard sphere and electrostatic residual correlation effects contribute with the same sign for counterions, whereas Figs. 14(c) and 14(f) show a significant competition between both contributions due to the opposite sign of each contribution for the co-ion.

In summary, we find that the contributions coming from beyond mean-field approximations play a dominant role in the behavior of the corresponding ion density profiles, depending on the aqueous electrolyte model used to describe the ionic atmosphere around the B-DNA.

In this work, we have introduced a classical density functional theory that is able to describe the influence of surrounding ions and water molecules on the structural and thermodynamic properties of polyelectrolytes in solution at infinite dilution. This approach goes beyond the nonlinear classical Poisson-Boltzmann theory by accounting for electrostatic ion correlations, size asymmetry, and excluded volume effects of ions and solvent particles. We have implemented a combination of the MSA and FMTWBII theories to study the influence of electrolyte mixtures of monovalent and multivalent ions and explicit water molecules at different densities and particle's sizes. The water crowding effects are taken into account by considering a neutral ion species of hard spheres representing explicit water molecules. The theoretical ionic density profiles display a good agreement with Monte Carlo simulations in the primitive model. The IC and MEP predicted by our DFT theory shows charge inversion at high electrolyte mixture concentrations.

Several models have been used to study the distribution of electrolyte mixtures surrounding a cylindrical rigid polyelectrolyte mimicking a B-DNA molecule at various ionic conditions. At microscopic level, we observe significant discrepancies in the ionic density profiles predicted by the different aqueous electrolyte models used in this work. However, the corresponding excess number of ions adsorbed to DNA are very similar in a wide range of electrolyte concentrations. Therefore, experimental data on the number of excess ions should not be considered alone, in general, as a conclusive test to determine the accuracy of DNA-electrolyte molecular models.10,70

On the other hand, the ionic potential of mean force has been analyzed and decomposed theoretically using our DFT approach in different contributions that depend on both electrostatic ion correlations and ionic excluded volume effects. At high electrolyte concentrations, the largest contribution to the ionic potential of mean force can be attributed to the electrostatic residual component in the presence of hydrated divalent counterions for the primitive model. At analogous conditions for the SPM model with Shannon's ionic sizes, it has been shown that the oscillatory behavior of the ionic potential of mean force associated to co-ions is promoted mainly by the hard sphere contribution.

Overall, the proposed theoretical DFT approach can be useful for studying ion and solvent dependent structural and thermodynamic properties of highly charged polyelectrolytes at low computation cost in comparison with explicit solvent simulations.

A better understanding of the ionic atmosphere around polymers, nanoparticles, and biological molecules requires the consideration of inhomogeneous charge distributions, irregular shapes, solvent polarization effects, and hydrogen bonding. Work along these directions is currently in progress.

We thank the National Science Foundation (NSF)-PREM Grant No. DMR-0934218, G.I.G.G. and M.O.C. thank the support of the NSF MRSEC Award No DMR-1121262. The computer cluster where the simulations were performed was funded by the office of the Director of Defense Research and Engineering (DDR&E) and the (U.S.) Air Force Office of Scientific Research (USAFOSR) under Award No. FA9550-10-1-0167. One of us (M.O.F.) would like to acknowledge support from National Institutes of Health (NIH)-NIGMS Grant No. 5R44GM073391 and thank Mr. Travis Mackoy for his help in preparing Fig. 1 and Dr. Robert C. Harris for providing some comments about the paper.

1. Density functional theory for inhomogeneous fluid mixtures

The density functional theory for inhomogeneous fluid mixtures is based on the definition of the exact Grand Potential functional of the ion density profiles Ω[{ρj}] from which the structural and thermodynamic properties of the system can be obtained.34 Specifically, the Legendre transform of the Helmholtz Free energy of the system F[{ρj}] defines the Grand Potential functional as follows:71 

\begin{equation}\Omega [\lbrace \rho _{j}\rbrace ]=F[\lbrace \rho _{j}\rbrace ]+\sum _{i=1}^{m}\int d\mathbf {r}\rho _{i}(r)\lbrace u_{i}(\mathbf {r})-\mu _{i}\rbrace ,\end{equation}
Ω[{ρj}]=F[{ρj}]+i=1mdrρi(r){ui(r)μi},
(A1)

where μi represents the chemical potential of the ith ionic component. The Helmholtz free energy functional is generally expressed as the sum of the ideal gas Fid and the excess Fex free energies

\begin{equation}F[\lbrace \rho _{j}\rbrace ]=F^{id}[\lbrace \rho _{j}\rbrace ]+F^{ex}[\lbrace \rho _{j}\rbrace ].\end{equation}
F[{ρj}]=Fid[{ρj}]+Fex[{ρj}].
(A2)

The ideal gas free energy is known exactly but approximations must be developed for the excess free energy. At equilibrium the Grand Potential functional is minimal with respect to variations in the density profiles yielding the following formal exact expression for the singlet ion density profiles {ρj(r)}:

\begin{eqnarray}\hspace*{-10pt}\rho _{i}(\mathbf {r})&=&\rho _{i}^{0}exp\lbrace \Delta c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )\nonumber\\&&-\beta q_{i}\psi ^{Solute}(\mathbf {r})-\beta u_{i}^{hs}(\mathbf {r})\rbrace,\quad i=1\ldots m,\hspace*{10pt}\end{eqnarray}
ρi(r)=ρi0exp{Δci(1)(r;{ρj})βqiψSolute(r)βuihs(r)},i=1...m,
(A3)

where β is the inverse of the thermal energy and the one particle direct correlation function for inhomogeneous fluids

$c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)(r;{ρj}) is defined as

\begin{equation}c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )=-\beta \frac{\delta F^{ex}[\lbrace \rho _{j}\rbrace ]}{\delta \rho _{i}(\mathbf {r})},\qquad i=1\ldots m,\end{equation}
ci(1)(r;{ρj})=βδFex[{ρj}]δρi(r),i=1...m,
(A4)

and

$u_{i}(\mathbf {r})\equiv q_{i}\psi ^{Solute}(\mathbf {r})+u_{i}^{hs}(\mathbf {r})$
ui(r)qiψSolute(r)+uihs(r)⁠. Note that ψSolute(r) represents the MEP generated by the polyelectrolyte surface charge. As a result, qiψSolute(r), where r = ‖r‖ determines the electric potential energy required by the fixed polyelectrolyte surface charge to bring an ion of species i from infinity to a distance r from the center of the polyelectrolyte

\begin{equation}\begin{array}{c}\textrm { }\\[3pt]q_{i}\psi ^{Solute}(r)=\left\lbrace \begin{array}{ll}-q_{i}\frac{2\lambda }{\epsilon }ln(r), & \textnormal { }\quad r>R,\\[6pt]\infty ,\textnormal { }\textrm { } & \quad \textrm { otherwise,} \end{array}\right.\quad i=1\ldots m. \end{array}\end{equation}
qiψSolute(r)=qi2λεln(r),r>R,, otherwise ,i=1...m.
(A5)

On the other hand,

$u_{i}^{hs}(r)$
uihs(r) represents the excluded volume potential energy generated by the cylindrical polyelectrolyte on an ion of species i located at a distance r from the center of the polyelectrolyte

\begin{equation}\begin{array}{c}\textrm { }\\[3pt]u_{i}^{hs}(r)=\left\lbrace \begin{array}{ll}0,\textnormal { }\quad r>R+\sigma _{i}/2, & \textnormal { }\\[6pt]\infty ,\textnormal { }\quad \textrm { otherwise,} \end{array}\right.\quad i=1\ldots m. \end{array}\end{equation}
uihs(r)=0,r>R+σi/2,, otherwise ,i=1...m.
(A6)

By replacing the latter expression into Eq. (A3) we obtain the following new expression for the ion (Boltzmann) distributions:

\begin{equation}\begin{array}{c}\textrm { }\\\rho _{i}(\mathbf {r})=\left\lbrace \begin{array}{ll}\rho _{i}^{0}exp\big\lbrace \Delta c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )-\beta q_{i}\psi ^{Solute}(r)\big\rbrace ,\textnormal { }\quad r>R+\sigma _{i}/2, & \textnormal { }\\[10pt]0,\textnormal { }\qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad \hspace{10.0pt}\,\,\, r\le R+\sigma _{i}/2, \end{array}\right.\quad i=1\ldots m. \end{array}\end{equation}
ρi(r)=ρi0expΔci(1)(r;{ρj})βqiψSolute(r),r>R+σi/2,0,rR+σi/2,i=1...m.
(A7)

Using the decomposition for the excess Helmholtz Free energy of the system Fex[{ρj}] given by the expression (2) and setting

$F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]=F_{res}^{ex}[\lbrace \rho _{j}\rbrace ]=0$
Fhsex[{ρj}]=Fresex[{ρj}]=0⁠, the first functional derivative of the excess free energy yields

\begin{eqnarray*}c_{i}^{(1)}(\mathbf {r};\lbrace \rho _{j}\rbrace )&=&-\beta \frac{\delta F^{ex}[\lbrace \rho _{j}\rbrace ]}{\delta \rho _{i}(\mathbf {r})}=-\beta q_{i}\sum _{k}q_{k}\int d\mathbf {r}^{\prime }\frac{\rho _{k}(\mathbf {r}^{\prime })}{|\mathbf {r}-\mathbf {r}^{\prime }|}\nonumber\\&=&-\beta q_{i}\psi ^{ion}(\mathbf {r};\lbrace \rho _{j}\rbrace ),\,\,\,\, i=1\ldots m,\end{eqnarray*}
ci(1)(r;{ρj})=βδFex[{ρj}]δρi(r)=βqikqkdrρk(r)|rr|=βqiψion(r;{ρj}),i=1...m,

where ψion(r; {ρj}) represents the MEP due to the ionic distributions. By replacing the latter expression into Eq. (A7) we recover the conventional Boltzmann distribution approximation for the ionic species i

\begin{eqnarray}&&\hspace{-5pt}\rho _{i}(\mathbf {r})\nonumber\\&&\hspace{-5pt}\;\;=\!\left\lbrace\! \begin{array}{ll}\rho _{i}^{0}\text{exp}\lbrace -\beta q_{i}\psi (\mathbf {r};\lbrace \rho _{j}\rbrace \rbrace ,\textnormal { }\;\; r>R+\sigma _{i}/2, & \textnormal { }\\[10pt]0,\textnormal { }\qquad \qquad \qquad {\;\;\;\; \hspace{10.0pt}\quad r\le R+\sigma _{i}/2,} \end{array}\right.\;\;\! i\!=\!1\ldots m,\nonumber\\\end{eqnarray}
ρi(r)=ρi0exp{βqiψ(r;{ρj}},r>R+σi/2,0,rR+σi/2,i=1...m,
(A8)

where ψ(r; {ρj}) ≡ ψSolute(r; {ρj}) + ψion(r; {ρj}) represents the MEP due to the polyelectrolyte surface charge plus the ion distribution. The fundamental expressions for the ion distributions given in Eq. (4) are obtained by substituting Eq. (2) into Eq. (A4) and then substituting the resulting expression into Eq. (A7).

Under this framework, most of the relevant thermodynamic properties of ionic atmosphere (electrical double layers) can be obtained from the expressions for the ion distributions given in Eq. (4). For instance, the coordination number ni(r) gives the average ion number of species i that will be found in an imaginary cylinder of radius r centered on the polyelectrolyte axis. It is calculated in terms of the ion density profiles by the following expression:

\begin{equation}n_{i}(r)=2\pi \int _{R}^{r}\rho _{i}(r^{\prime })r^{\prime }dr^{\prime },\quad i=1\ldots m.\end{equation}
ni(r)=2πRrρi(r)rdr,i=1...m.
(A9)

Equation (10) and Gauss's law can be used to calculate the MEP ψ(r; {ρj}) due to the polyelectrolyte surface and the ion distributions in a solvent with isotropic and constant relative permittivity ε. Alternatively, it can be obtained by solving the Poisson equation,

\begin{equation}\nabla ^{2}\psi (\mathbf {r};\lbrace \rho _{j}\rbrace )=-\frac{1}{\epsilon }\sum _{i=1}^{m}\rho _{i}(r)z_{i},\quad i=1\ldots m,\end{equation}
2ψ(r;{ρj})=1εi=1mρi(r)zi,i=1...m,
(A10)

and using the boundary conditions ψ(r; {ρj})|r → ∞ → 0, ε∂ψ(r; {ρj})/∂r|r = R = −σ = −λ/(2πR), and the electroneutrality condition

\begin{equation*}\lambda +2\pi \int _{R}^{\infty }\sum _{i}z_{i}e\rho _{i}(r^{\prime })r^{\prime }dr^{\prime }=0,\quad i=1\ldots m.\end{equation*}
λ+2πRizieρi(r)rdr=0,i=1...m.

In particular, the formal solution of Eq. (A10) in cylindrical coordinates reads

\begin{eqnarray}&&\psi (r;\lbrace \rho _{j}\rbrace )=\frac{4\pi e}{\epsilon }\int _{r}^{\infty }\sum _{i}z_{i}\rho _{i}(r^{\prime })r^{\prime }\ln (r/r^{\prime })dr^{\prime }\nonumber\\&&\qquad\qquad r>R,\quad i=1\ldots m.\end{eqnarray}
ψ(r;{ρj})=4πeεriziρi(r)rln(r/r)drr>R,i=1...m.
(A11)

When the latter equation is substituted into Eq. (4) it yields an implicit equation for the ion profiles {ρi(r)}. The expression corresponding to the hard sphere one particle direct correlation functions for inhomogeneous fluids

$c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )$
ci(1)hs(r;{ρj}) is provided in Subsection A 2 of the Appendix. It plays a key role in the formulation of the approach introduced in this work.

2. The White Bear version II of FMT approach to calculate hard sphere direct correlation functions

The key assumption of FMT is that the hard sphere excess free energy functional for inhomogeneous has the form

\begin{equation}\beta F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]=\int \Phi _{hs}^{ex}[\lbrace \bar{\rho }_{a}(\mathbf {r})\rbrace ,\lbrace \bar{\bm {\rho }}_{b}(\mathbf {r})\rbrace ]d\mathbf {r},\end{equation}
βFhsex[{ρj}]=Φhsex[{ρ¯a(r)},{ρ¯b(r)}]dr,
(A12)

where the hard sphere free energy density

$\Phi _{hs}^{ex}$
Φhsex (per thermal energy unit) is a function of a set of scalar and vector weighted densities, each one defined in the form of Eq. (A12), i.e.,

\begin{equation*}\bar{\rho }_{a}(\mathbf {r})=\sum _{j=1}^{m}\int d\mathbf {r}^{\prime }w_{j}^{a}(\mathbf {r},{\bf r^{\prime }})\rho _{j}({\bf r}^{\prime }),\quad a=0,1,2,3\end{equation*}
ρ¯a(r)=j=1mdrwja(r,r)ρj(r),a=0,1,2,3

and

\begin{equation*}\bar{\bm {\rho }}_{b}(\mathbf {r})=\sum _{j=1}^{m}\int d\mathbf {r}^{\prime }\mathbf {w}_{j}^{b}(\mathbf {r},{\bf r^{\prime }})\rho _{j}({\bf r}^{\prime }),\quad b=1,2.\end{equation*}
ρ¯b(r)=j=1mdrwjb(r,r)ρj(r),b=1,2.

Under this framework, the weighted density functions are chosen to be equal to those characteristic functions present in the solution for the two particle direct correlation functions in the MSA, e.g.,

\begin{eqnarray*}w_{j}^{0}(\mathbf {r},{\bf r^{\prime }})=\frac{\delta (\bm{|}\mathbf {r}-{\bf r^{\prime }\bm{|}}-\sigma _{j}/2)}{\pi \sigma _{j}^{2}},&& \quad w_{j}^{1}(\mathbf {r},{\bf r^{\prime }})=\frac{\delta (\bm{|}\mathbf {r}-{\bf r^{\prime }\bm{|}}-\sigma _{j}/2)}{2\pi \sigma _{j}},\\\end{eqnarray*}
wj0(r,r)=δ(|rr|σj/2)πσj2,wj1(r,r)=δ(|rr|σj/2)2πσj,
\begin{eqnarray*}&&\\[-56pt]w_{j}^{2}(\mathbf {r},{\bf r^{\prime }})\!=\!\delta (\bm{|}\mathbf {r}\!-\!{\bf r^{\prime }\bm{|}}-\sigma _{j}/2),&&\quad w_{j}^{3}(\mathbf {r},{\bf r^{\prime }})\!=\!H(\bm{|}\mathbf {r}\!-\!{\bf r^{\prime }\bm{|}}\!-\!\sigma _{j}/2),\\\mathbf {w}_{j}^{1}(\mathbf {r},{\bf r^{\prime }})=\frac{\mathbf {r}-{\bf r^{\prime }}}{2\pi \sigma _{j}\bm{|}\mathbf {r}-{\bf r^{\prime }}\bm{|}}\delta \left(\bm{|}\mathbf {r}-{\bf r^{\prime }\bm{|}}-\sigma _{j}/2\right),&&\quad \mathbf {w}_{j}^{2}(\mathbf {r},{\bf r^{\prime }})=\frac{\mathbf {r}-{\bf r^{\prime }}}{\bm{|}\mathbf {r}-{\bf r^{\prime }}\bm{|}}\delta \left(\bm{|}\mathbf {r}-{\bf r^{\prime }}\bm{|}-\sigma _{j}/2\right).\end{eqnarray*}
wj2(r,r)=δ(|rr|σj/2),wj3(r,r)=H(|rr|σj/2),wj1(r,r)=rr2πσj|rr|δ|rr|σj/2,wj2(r,r)=rr|rr|δ|rr|σj/2.

The hard sphere correlation function is obtained when the expressions for the weighted densities are substituted into the expression for the excess free energy density for inhomogeneous fluids. The expression (5) is obtained from a third order expansion in the dimensionless scaled parameter ξ3 of the following free energy density for homogeneous fluids:

\begin{eqnarray}\Phi &=&\Phi _{SPT}+\phi _{4}(\xi _{3})\big[-c_{1}\xi _{2}^{3}\xi _{3}+b_{2}\xi _{1}\xi _{2}\xi _{3}^{2}+a_{3}\xi _{0}\xi _{3}^{3}\big]\nonumber\\&&+\phi _{5}(\xi _{3})c_{2}\big[\xi _{2}^{3}\xi _{3}^{2}-12\pi \xi _{1}\xi _{2}\xi _{3}^{2}\big],\end{eqnarray}
Φ=ΦSPT+ϕ4(ξ3)c1ξ23ξ3+b2ξ1ξ2ξ32+a3ξ0ξ33+ϕ5(ξ3)c2ξ23ξ3212πξ1ξ2ξ32,
(A13)

where

\begin{equation*}\phi _{4}(\xi _{3})=\frac{\frac{3}{2}\xi _{3}^{2}-\xi _{3}-(1-\xi _{3})^{2}\ln (1-\xi _{3})}{\xi _{3}^{3}(1-\xi _{3})^{2}}\end{equation*}
ϕ4(ξ3)=32ξ32ξ3(1ξ3)2ln(1ξ3)ξ33(1ξ3)2

and

\begin{equation*}\phi _{5}(\xi _{3})=\frac{\frac{9}{2}\xi _{3}^{2}-\xi _{3}^{3}-3\xi _{3}-3(1-\xi _{3})^{2}\ln (1-\xi _{3})}{\xi _{3}^{4}(1-\xi _{3})^{2}},\end{equation*}
ϕ5(ξ3)=92ξ32ξ333ξ33(1ξ3)2ln(1ξ3)ξ34(1ξ3)2,

and c1 = 1/(18π), b2 = 1/3, c2 = 1/(36π), a3 = c3 = 0, and b3 = −1/3. The first term in Eq. (A13) corresponds to the free energy density obtained with the scaled particle theory72,73 and the other terms introduce corrections in such a way that the corresponding equation of state (5) recovers the quasi-exact Carnahan-Starling result in the case of a one component fluid.45 Using the standard homogeneous-to-inhomogeneous generalization method, the expression for the approximate excess free energy for inhomogeneous fluids reads

\begin{eqnarray}&&\Phi _{hs}^{ex}[\lbrace \bar{\rho }_{a}(\mathbf {r})\rbrace ,\lbrace \bar{\bm {\rho }}_{b}(\mathbf {r})\rbrace ]\nonumber\\&&\quad\equiv \Phi _{SPT}[\lbrace \bar{\rho }_{a}(\mathbf {r})\rbrace ]+\phi _{4}(\bar{\rho }_{3}(\mathbf {r}))\nonumber\\&&\qquad\times[-c_{1}\lbrace \bar{\rho }_{2}(\mathbf {r})^{3}-3\bar{\rho }_{2}(\mathbf {r})\bar{\bm {\rho }}_{2}(\mathbf {r})\cdot \bar{\bm {\rho }}_{2}(\mathbf {r})\rbrace \bar{\rho }_{3}(\mathbf {r})\nonumber\\&&\qquad+b_{2}\lbrace \bar{\rho }_{1}(\mathbf {r})\bar{\rho }_{2}(\mathbf {r})-\bar{\bm {\rho }}_{1}(\mathbf {r})\cdot \bar{\bm {\rho }}_{2}(\mathbf {r})\rbrace \bar{\rho }_{3}(\mathbf {r})^{2}\nonumber\\&&\qquad+a_{3}\bar{\rho }_{0}(\mathbf {r})\bar{\rho }_{3}(\mathbf {r})^{3}]+\phi _{5}(\bar{\rho }_{3}(\mathbf {r}))c_{2}[\lbrace \bar{\rho }_{2}(\mathbf {r})^{3}\nonumber\\&&\qquad -3\bar{\rho }_{2}(\mathbf {r})\bar{\bm {\rho }}_{2}(\mathbf {r})\cdot \bar{\bm {\rho }}_{2}(\mathbf {r})\rbrace \bar{\rho }_{3}(\mathbf {r})^{2}\nonumber\\&&\qquad-12\pi \lbrace \bar{\rho }_{1}(\mathbf {r})\bar{\rho }_{2}(\mathbf {r})-\bar{\bm {\rho }}_{1}(\mathbf {r})\cdot \bar{\bm {\rho }}_{2}(\mathbf {r})\rbrace \bar{\rho }_{3}(\mathbf {r})^{2}].\qquad\quad\end{eqnarray}
Φhsex[{ρ¯a(r)},{ρ¯b(r)}]ΦSPT[{ρ¯a(r)}]+ϕ4(ρ¯3(r))×[c1{ρ¯2(r)33ρ¯2(r)ρ¯2(r)·ρ¯2(r)}ρ¯3(r)+b2{ρ¯1(r)ρ¯2(r)ρ¯1(r)·ρ¯2(r)}ρ¯3(r)2+a3ρ¯0(r)ρ¯3(r)3]+ϕ5(ρ¯3(r))c2[{ρ¯2(r)33ρ¯2(r)ρ¯2(r)·ρ¯2(r)}ρ¯3(r)212π{ρ¯1(r)ρ¯2(r)ρ¯1(r)·ρ¯2(r)}ρ¯3(r)2].
(A14)

Having the expression for the

$\Phi _{hs}^{ex}[\lbrace \bar{\rho }_{a}(\mathbf {r})\rbrace ,\lbrace \bar{\bm {\rho }}_{b}(\mathbf {r})\rbrace ]$
Φhsex[{ρ¯a(r)},{ρ¯b(r)}]⁠, the hard sphere one particle direct correlation function for inhomogeneous fluids is obtained by differentiation of the free energy as follows:

\begin{equation}c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}\rbrace )=-\beta \frac{\delta F_{hs}^{ex}[\lbrace \rho _{j}\rbrace ]}{\delta \rho _{i}(\mathbf {r})},\qquad i=1\ldots m.\end{equation}
ci(1)hs(r;{ρj})=βδFhsex[{ρj}]δρi(r),i=1...m.
(A15)

The corresponding expression for the homogeneous fluids

$c_{i}^{(1)hs}(\mathbf {r};\lbrace \rho _{j}^{0}\rbrace )$
ci(1)hs(r;{ρj0}) can be subsequently obtained from Eq. (A14) by replacing
$\lbrace \bar{\rho }_{a}(\mathbf {r})\rbrace \rightarrow \lbrace \xi _{a}\rbrace$
{ρ¯a(r)}{ξa}
and by setting
$\lbrace \bar{\bm {\rho }}_{b}(\mathbf {r})\rbrace =0$
{ρ¯b(r)}=0
.

3. Residual (electrostatic) correlation functions

Another important contribution to the ion density profiles appearing in the argument of the exponential in Eq. (4) is given by the remaining (residual) one particle correlation function for inhomogeneous fluids. It is usually approximated by using functional Taylor expansion at first order in power of the ion density profiles around uniform fluids as follows:74 

\begin{eqnarray}&&\triangle c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )=\!\sum _{j=1}^{m}\!\int\! d\mathbf {r}^{\prime }c_{ij}^{(2)res}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{j}^{0}\rbrace )\,\big[\rho _{j}({\bf r}^{\prime })-\rho _{j}^{o}\big],\nonumber\\&&\quad i=1\ldots m,\end{eqnarray}
ci(1)res(r;{ρj})=j=1mdrcij(2)res(r,r;{ρj0})ρj(r)ρjo,i=1...m,
(A16)

where the following relationship is used:

\begin{eqnarray}&&c_{ij}^{(2)res}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )\nonumber\\&&\;\;\,=\left[\frac{\delta c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{k}\rbrace )}{\delta \rho _{i}(\mathbf {r})}\right]_{\lbrace \rho _{k}\rbrace =\lbrace \rho _{k}^{o}\rbrace },\;\; i,j=1\ldots m.\qquad\end{eqnarray}
cij(2)res(r,r;{ρk0})=δci(1)res(r;{ρk})δρi(r){ρk}={ρko},i,j=1...m.
(A17)

The latter equation implies that

$c_{ij}^{(2)res}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )\break =\;c_{ij}^{(2)}(\mathbf {r},\;{\bf r^{\prime }};\;\lbrace \rho _{k}^{0}\rbrace )\,+\,z_{i}z_{j}e^{2}/(\epsilon |\mathbf {r}\;-\;\mathbf {r}^{\prime }|)$
cij(2)res(r,r;{ρk0})=cij(2)(r,r;{ρk0})+zizje2/(ε|rr|)
$\,-\,c_{ij}^{(2)hs}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )$
cij(2)hs(r,r;{ρk0})
where the expression for
$c_{ij}^{(2)hs}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )$
cij(2)hs(r,r;{ρk0})
is given by Eq. (A15). An advantageous approximate explicit expression for the two particle correlation functions for multicomponent homogeneous fluids
$c_{ij}^{(2)}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )\break \equiv c_{ij}^{(2)short-MSA}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )+c_{ij}^{(2)el-MSA}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )$
cij(2)(r,r;{ρk0})cij(2)shortMSA(r,r;{ρk0})+cij(2)elMSA(r,r;{ρk0})
is provided by Hiroike.42 

The two particle correlation function for homogeneous fluids predicted by the conventional MSA reduces to the PY correlation function when the ions are not charged, such that

$c_{ij}^{(2)hs}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )=c_{ij}^{(2)short-MSA}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )$
cij(2)hs(r,r;{ρk0})=cij(2)shortMSA(r,r;{ρk0})⁠.75 Here, it is assumed that the following difference for short-range two particle correlation functions between the MSA and Hansen-Goos and Roth FMT approximations:

\begin{eqnarray*}&&\sum _{j=1}^{m}\int d\mathbf {r}^{\prime }\big[c_{ij}^{(2)short-MSA}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{j}^{0}\rbrace )\nonumber\\&&\quad-c_{ij}^{(2)hs-Roth}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{j}^{0}\rbrace )\big]\big[\rho _{j}({\bf r}^{\prime })-\rho _{j}^{o}\big],\quad i=1\ldots m\end{eqnarray*}
j=1mdr[cij(2)shortMSA(r,r;{ρj0})cij(2)hsRoth(r,r;{ρj0})]ρj(r)ρjo,i=1...m

is negligible compared with the corresponding contribution coming from the electrostatic ion-ion correlations. Under this assumption, the residual one particle direct correlation function for inhomogeneous fluids is approximated as follows:

\begin{eqnarray}&&\triangle c_{i}^{(1)res}(\mathbf {r};\lbrace \rho _{j}\rbrace )\nonumber\\&&\quad =\sum _{j=1}^{m}\int d\mathbf {r}^{\prime }\big[c_{ij}^{(2)el-MSA}(\mathbf {r},{\bf r^{\prime }};\lbrace \rho _{k}^{0}\rbrace )\nonumber\\&&\qquad+z_{i}z_{j}e^{2}/(\epsilon |\mathbf {r}-\mathbf {r}^{\prime }|)\big][\rho _{j}({\bf r}^{\prime })-\rho _{j}^{o}],\quad i=1\ldots m.\nonumber\\\end{eqnarray}
ci(1)res(r;{ρj})=j=1mdr[cij(2)elMSA(r,r;{ρk0})+zizje2/(ε|rr|)][ρj(r)ρjo],i=1...m.
(A18)
4. MSA approach to calculate electrical direct correlation functions

We use the analytic expressions for c(2)res obtained by Hiroike using the MSA,42 when σi < σj and 0 ≦ r ≦ (σj − σi)/2

\begin{eqnarray}&&\nonumber\\[-40pt]c_{ij}^{(2)res}(r) &=&-\frac{2e^{2}}{\epsilon kT}[-z_{i}N_{j}+X_{i}(N_{i}+\Gamma X_{i})\nonumber\\&&-(\sigma _{i}/3)(N_{i}+\Gamma X_{i})^{2}]\quad i,j=1\ldots m.\qquad\end{eqnarray}
cij(2)res(r)=2e2εkT[ziNj+Xi(Ni+ΓXi)(σi/3)(Ni+ΓXi)2]i,j=1...m.
(A19)

When (σi − σj)/2≦r < (σi + σj)/2,

\begin{eqnarray}{c}c_{ij}^{(2)res}(r)&=&\frac{e^{2}}{\epsilon kT}\Bigg ((\sigma _{i}-\sigma _{j})\Big \lbrace \frac{(X_{i}+X_{j})}{4}[(N_{i}+\Gamma X_{i})-(N_{j}+\Gamma X_{j})]\nonumber\\&&-\frac{(\sigma _{i}-\sigma _{j})}{16}[(N_{i}+\Gamma X_{i}+N_{j}+\Gamma X_{j})^{2}-4N_{i}N_{j}]\Big \rbrace /r\nonumber\\&&-\Big \lbrace (X_{i}-X_{j})(N_{i}-N_{j})+\big(X_{i}^{2}+X_{j}^{2}\big)\Gamma +(\sigma _{i}+\sigma _{j})N_{i}N_{j}\nonumber\\&&-\frac{1}{3}[\sigma _{i}(N_{i}+\Gamma X_{i})^{2}+\sigma _{j}(N_{j}+\Gamma X_{j})^{2}]\Big \rbrace \qquad\qquad\qquad\qquad\qquad i,j=1\ldots m,\\&&+\Big \lbrace \frac{X_{i}}{\sigma _{i}}(N_{i}+\Gamma X_{i})+\frac{X_{j}}{\sigma _{j}}(N_{j}+\Gamma X_{j})+N_{i}N_{j}\nonumber\\&&-\frac{1}{2}[(N_{i}+\Gamma X_{i})^{2}+(N_{j}+\Gamma X_{j})^{2}]\Big \rbrace r \nonumber\\&&+\Big \lbrace \frac{(N_{i}+\Gamma X_{i})^{2}}{6\sigma _{i}^{2}}+\frac{(N_{j}+\Gamma X_{j})^{2}}{6\sigma _{j}^{2}}\Big \rbrace r^{3}\Bigg ), \nonumber\end{eqnarray}
ccij(2)res(r)=e2εkT((σiσj){(Xi+Xj)4[(Ni+ΓXi)(Nj+ΓXj)](σiσj)16[(Ni+ΓXi+Nj+ΓXj)24NiNj]}/r{(XiXj)(NiNj)+Xi2+Xj2Γ+(σi+σj)NiNj13[σi(Ni+ΓXi)2+σj(Nj+ΓXj)2]}i,j=1...m,+{Xiσi(Ni+ΓXi)+Xjσj(Nj+ΓXj)+NiNj12[(Ni+ΓXi)2+(Nj+ΓXj)2]}r+(Ni+ΓXi)26σi2+(Nj+ΓXj)26σj2r3),
(A20)

where Xi is defined as

\begin{eqnarray}X_{i}&\!=\!&\frac{z_{i}}{1+\Gamma \sigma _{i}}-\frac{c\sigma _{i}^{2}}{1+\Gamma \sigma _{i}}\nonumber\\&&\times\frac{\sum _{j=1}^{m}\rho _{j}\sigma _{j}z_{j}(1+\Gamma \sigma _{j})^{-1}}{1+c\sum _{j=1}^{m}\rho _{j}\sigma _{j}^{3}(1+\Gamma \sigma _{j})^{-1}},\quad i\!=\!1\ldots m\qquad\quad\end{eqnarray}
Xi=zi1+Γσicσi21+Γσi×j=1mρjσjzj(1+Γσj)11+cj=1mρjσj3(1+Γσj)1,i=1...m
(A21)

and

\begin{equation}N_{i}=\frac{(X_{i}-z_{i})}{\sigma _{i}},\qquad i=1\ldots m.\end{equation}
Ni=(Xizi)σi,i=1...m.
(A22)

The parameter Γ in Ref. 42 is obtained numerically by solving the set of nonlinear equations given by

\begin{equation}\Gamma ^{2}=\frac{\pi e^{2}}{\epsilon kT}\sum _{i=1}^{m}\rho _{i}^{0}X_{i}^{2}\end{equation}
Γ2=πe2εkTi=1mρi0Xi2
(A23)

and Eq. (A21).

5. Algorithm

The theoretical model described above yields a set of implicit nonlinear equations that are used to calculate the ion density profiles. The corresponding numerical solution of the coupled Eqs. (4), (A10), (A12), (A14), (A15), and (A18) are obtained using iterative techniques based on the Picard algorithm. The numerical solution for the ion density profiles is obtained with an accuracy of 6 digits of precision, and satisfies the electroneutrality condition.

The nonlinear equations of MSA are solved by tensor methods using Algorithm 768 of ACM: TENSOLVE software package, which is suitable for small-medium problems with up to 100 unknowns.76 The three-dimensional integrals involved in the computation of correlation functionals are reduced to one-dimensional integrals along the radial direction due to the cylindrical symmetry of the polyelectrolyte. These one-dimensional integrals are evaluated using the trapezoidal algorithm with a linear mesh resolution of the order 0.05-0.1 Å. The cutoff (rmax) for the grid used in the theoretical calculations is varied between 60 Å and 100 Å depending on the ionic solution conditions.

1.
J. D.
Ballin
,
D. M.
Wilson
, and
E. M.
Tadashi
,
Role and Applications of Electrostatic Effects on Nucleic Acid Conformational Transitions and Binding Processes: In Application of Thermodynamics to Biological and Materials Science
(
InTech
,
Rijeka
,
2011
).
2.
G. C. L.
Wong
and
L.
Pollack
,
Annu. Rev. Phys. Chem.
61
,
171
(
2010
).
3.
M. T.
Record
,
T. M.
Lohman
, and
P. L.
de Haseth
,
J. Mol. Biol.
107
,
145
(
1976
).
4.
G. S.
Manning
,
Q. Rev. Biophys.
11
,
179
(
1978
).
6.
S. A.
Pabit
,
S. P.
Meisburger
,
L.
Li
,
J. M.
Blose
,
C. D.
Jones
, and
L.
Pollack
,
J. Am. Chem. Soc.
132
,
16334
(
2010
).
7.
Y.
Bai
,
M.
Greenfeld
,
K. J.
Travers
,
S. D. V. B.
Chu
,
J.
Lipfert
, and
D.
Herschlag
,
J. Am. Chem. Soc.
129
,
14981
(
2007
).
8.
A. A.
Chen
,
M.
Marucho
,
N. A.
Baker
, and
R. V.
Pappu
,
Methods Enzymol.
469
,
411
(
2009
).
9.
S.
Kirmizialtin
and
R.
Elber
,
J. Phys. Chem. B
114
,
8207
(
2010
).
10.
J.
Yoo
and
A.
Aksimentiev
,
J. Phys. Chem. B
116
,
12946
(
2012
).
11.
M. A.
Young
,
B.
Jayaram
, and
D. L.
Beveridge
,
J. Am. Chem. Soc.
119
,
59
(
1997
).
12.
M. O.
Fenley
,
C.
Russo
, and
G. S.
Manning
,
J. Phys. Chem. B
115
,
9864
(
2011
).
13.
P. L.
Privalov
,
A. I.
Dragan
, and
C.
Crane-Robinson
,
Nucl. Acids Res.
39
,
2483
(
2011
).
14.
C.
García-García
and
D. E.
Draper
,
J. Mol. Biol.
331
,
75
(
2003
).
15.
K. A.
Sharp
,
R. A.
Friedman
,
V.
Misra
,
J.
Hecht
, and
B.
Honig
,
Biopolymers
36
,
245
(
1995
).
16.
M. O.
Fenley
,
R. C.
Harris
,
B.
Jayaram
, and
A. H.
Boschitsch
,
Biophys. J.
99
,
879
(
2010
).
17.
E.
Gonzáles-Tovar
,
M.
Lozada-Cassou
, and
D.
Henderson
,
J. Chem. Phys.
83
,
361
(
1985
).
18.
N.
Bagatella-Flores
and
P.
González-Mozuelos
,
J. Chem. Phys.
117
,
6133
(
2002
).
19.
G. I.
Guerrero-García
,
P.
González-Mozuelos
, and
M. O.
de la Cruz
,
J. Chem. Phys.
135
,
164705
(
2011
).
20.
A.
Martín-Molina
,
J. G.
Ibarra-Armenta
,
E.
González-Tovar
,
R.
Hidalgo-Álvarez
, and
M.
Quesada-Párez
,
Soft Matter
7
,
1441
(
2011
).
21.
D.
Frydel
and
Y.
Levin
,
J. Chem. Phys.
137
,
164703
(
2012
).
22.
M.
Nedyalkova
,
S.
Madurga
,
S.
Pisov
,
E. V. I.
Pastor
, and
F.
Mas
,
J. Chem. Phys.
137
,
174701
(
2012
).
23.
A.
Bandopadhyay
and
S.
Chakraborty
,
Electrophoresis
34
,
2193
(
2013
).
24.
A.
Moncho-Jordá
,
J. A.
Anta
, and
J.
Callejas-Fernández
,
J. Chem. Phys.
138
,
134902
(
2013
).
25.
F.-H.
Wang
,
Y.-Y.
Wu
, and
Z.-J.
Tan
,
Biopolymers
99
,
370
(
2013
).
26.
G. I.
Guerrero-García
and
M. O.
de la Cruz
,
J. Chem. Theory Comput.
9
,
1
(
2013
).
27.
G. I.
Guerrero-García
,
Y.
Ying
, and
M. O.
de la Cruz
,
Soft Matter
9
,
6046
(
2013
).
28.
P.
González-Mozuelos
,
G. I.
Guerrero-García
, and
M. O.
de la Cruz
,
J. Chem. Phys.
139
,
064709
(
2013
).
29.
R.
Ramirez
and
D.
Borgis
,
J. Phys. Chem. B
109
,
6754
(
2005
).
30.
B.
Medasani
,
Z.
Ovanesyan
,
D. G.
Thomas
,
M. L.
Sushko
, and
M.
Marucho
,
J. Chem. Phys.
140
,
204510
(
2014
).
31.
M. G.
Knepley
,
D. A.
Karpeev
,
S.
Davidovits
,
R. S.
Eisenberg
, and
D.
Gillespie
,
J. Chem. Phys.
132
,
124101
(
2010
).
32.
J. J.
Howard
,
G. C.
Lynch
, and
B. M.
Pettitt
,
J. Phys. Chem. B
115
,
547
(
2011
).
33.
C. N.
Patra
and
A.
Yethiraj
,
J. Phys. Chem. B
103
,
6080
(
1999
).
34.
J.-P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
, 3rd ed. (
Academic Press
,
2006
).
35.
H. C.
Andersen
and
D.
Chandler
,
J. Chem. Phys.
53
,
547
(
1970
).
36.
J. S.
Rowlinson
,
Rep. Prog. Phys.
28
,
169
(
1965
).
37.
J. S.
Hoye
,
G.
Stell
, and
E.
Waisman
,
Mol. Phys.
52
,
1071
(
1976
).
38.
J. S.
Hoye
and
G.
Stell
,
J. Chem. Phys.
67
,
524
(
1977
).
39.
M.
Deserno
,
F.
Jiménez-Ángeles
,
C.
Holm
, and
M.
Lozada-Cassou
,
J. Phys. Chem. B
105
,
10983
(
2001
).
40.
M.
Quesada-Pérez
,
E.
González-Tovar
,
A.
Martín-Molina
,
M.
Lozada-Cassou
, and
R. H.
Álvarez
,
ChemPhysChem
4
,
234
(
2003
).
41.
A.
Kovalenko
and
F.
Hirata
,
J. Chem. Phys.
110
,
10095
(
1999
).
44.
H.
Hansen-Goos
and
R.
Roth
,
J. Chem. Phys.
124
,
154506
(
2006
).
45.
N. F.
Carnahan
and
K. E.
Starling
,
J. Chem. Phys.
51
,
635
(
1969
).
46.
J.
Hughes
,
E. J.
Krebs
, and
D.
Roundy
,
J. Chem. Phys.
138
,
024509
(
2013
).
47.
S. L.
Carnie
and
D. Y. C.
Chan
,
J. Chem. Phys.
73
,
2949
(
1980
).
48.
T.
Goel
,
C. N.
Patra
,
S. K.
Gosh
, and
T.
Mukherjee
,
J. Chem. Phys.
129
,
154707
(
2008
).
49.
T.
Goel
,
C. N.
Patra
,
S. K.
Gosh
, and
T.
Mukherjee
,
J. Phys. Chem. B
115
,
10903
(
2011
).
50.
C. N.
Patra
and
A.
Yethiraj
,
Biophys J.
78
,
699
(
2000
).
51.
K.
Wang
,
Y.-X.
Yu
, and
G.-H.
Gao
,
Phys. Rev. E
70
,
011912
(
2004
).
52.
K.
Wang
,
Y.-X.
Yu
,
G.-H.
Gao
, and
G.-S.
Luo
,
J. Chem. Phys.
123
,
234904
(
2005
).
53.
G.
Stell
and
S. F.
Sun
,
J. Chem. Phys.
63
,
5333
(
1975
).
54.
J. K.
Percus
and
G. J.
Yevick
,
Phys. Rev.
110
,
1
(
1958
).
55.
T.
Boublík
,
J. Chem. Phys.
53
,
471
(
1970
).
56.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
New York, USA
,
1989
).
57.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation
(
Academic
,
London
,
2002
).
58.
I.
Fukuda
,
Y.
Yonezawa
, and
H.
Nakamura
,
J. Chem. Phys.
134
,
164107
(
2011
).
59.
G.
Vernizzi
,
G. I.
Guerrero-García
, and
M. O.
de la Cruz
,
Phys. Rev. E
84
,
016707
(
2011
).
60.
I.
Fukuda
and
H.
Nakamura
,
Biophys. Rev.
4
,
161
(
2012
).
61.
P. X.
Viveros-Mendez
and
A.
Gil-Villegas
,
J. Chem. Phys.
136
,
154507
(
2012
).
62.
N.
Kamiya
,
I.
Fukuda
, and
H.
Nakamura
,
Chem. Phys. Lett.
568
,
26
(
2013
).
63.
G. I.
Guerrero-García
,
E.
González-Tovar
, and
M.
Chavez-Páez
,
Phys. Rev. E
80
,
021503
(
2009
).
64.
G. I.
Guerrero-García
,
E.
González-Tovar
, and
M. O.
de la Cruz
,
Soft Matter
6
,
2056
(
2010
).
65.
G. I.
Guerrero-García
,
E.
González-Tovar
, and
M. O.
de la Cruz
,
J. Chem. Phys.
135
,
054701
(
2011
).
66.
T.
Goel
,
C. N.
Patra
,
S. K.
Gosh
, and
T.
Mukherjee
,
J. Chem. Phys.
132
,
194706
(
2010
).
67.
R.
Messina
,
C.
Holm
, and
K.
Kremer
,
Eur. Phys. J. E.
4
,
363
(
2001
).
68.
R. D.
Shannon
,
Acta Crystallogr.
32
,
751
(
1976
).
69.
R. C.
Harris
,
A. H.
Boschitsch
, and
M. O.
Fenley
,
J. Chem. Phys.
140
,
075102
(
2014
).
70.
G. M.
Giambaşu
,
T.
Luchko
,
D.
Herschlag
,
D. M.
York
, and
D. A.
Case
,
Biophys. J.
106
,
883
(
2014
).
72.
H.
Reiss
,
H. L.
Frish
, and
J. L.
Lebowitz
,
J. Chem. Phys.
31
,
369
(
1959
).
73.
J. A.
Barker
and
D.
Henderson
,
Rev. Mod. Phys.
48
,
587
(
1976
).
74.
Y.-X.
Yu
,
J.
Wu
, and
G.-H.
Gao
,
J. Chem. Phys.
120
,
7223
(
2004
).
75.
E.
Waisman
and
J. L.
Lebowitz
,
J. Chem. Phys.
56
,
3086
(
1972
).
76.
A.
Bouaricha
and
R. B.
Schnable
,
ACM Trans. Math. Softw.
23
,
174
(
1997
).
77.
A. G.
Volkov
,
S.
Paula
, and
D. W.
Deamer
,
Bioelectrochem. Bioenerg.
42
,
153
(
1997
).
78.
Y.
Zhang
and
Z.
Xu
,
Am. Mineral.
80
,
670
(
1995
).
79.
A. H.
Boschitsch
and
M. O.
Fenley
,
J. Chem. Theory Comput.
7
,
1524
(
2011
).
80.
A. H.
Boschitsch
and
M. O.
Fenley
,
J. Comput. Chem.
28
,
909
(
2007
).