Interaction between two polyelectrolytes in monovalent aqueous salt solutions

We use the recently developed soft-potential-enhanced Poisson-Boltzmann (SPB) theory to study the interaction between two parallel polyelectrolytes (PEs) in monovalent ionic solutions in the weak-coupling regime. The SPB theory is fitted to ion distributions from coarse-grained molecular dynamics (MD) simulations and benchmarked against all-atom MD modelling for poly(diallyldimethylammonium) (PDADMA). We show that the SPB theory is able to accurately capture the interactions between two PEs at distances beyond the PE radius. For PDADMA positional correlations between the charged groups lead to locally asymmetric PE charge and ion distributions. This gives rise to small deviations from the SPB prediction that appear as short-range oscillations in the potential of mean force. Our results suggest that the SPB theory can be an efficient way to model interactions in chemically specific complex PE systems.


I. INTRODUCTION
Polyelectrolytes (PEs) are polymers containing electrolyte groups, which dissociate in aqueous solutions into solvated counterions and charged polymers. PEs exhibit physical properties very different from those of the uncharged polymers, in particular in terms of their aqueous solubility and tunability with salt [1, 2]. Consequently, PE solutions and PE assemblies strongly react to changes in solvent environment, temperature, pH, and salt conditions [3][4][5]. When their molecular persistence length is long enough, PEs can be considered, to some extent, as charged rods. Examples of such PEs include, e.g., functionalized cellulose fibrils, tobacco mosaic virus, DNA, and actin. The interactions between rod-like PEs are widely addressed in chemistry and biology, leading to both rich functional and assembly behavior, such as in complex assembly of DNA and RNA [6,7], cellulose nanocrystal based advanced self-assembly [8,9], and complexation based synthetic PE materials [3,10]. The extensive applications of PEs in the field of nanostructured materials span from responsive materials [3], surface modification [11] to layer-by-layer assembly [12,13].
The fundamentals to understanding self-assembly, as well as structural and dynamical properties of dense PE systems, to a first approximation, come from the PE-PE pairwise interactions. These interactions in ionic solutions are however complex due to effects rising from both charge correlations and solution conditions. Replacing atomistically detailed models with lower resolution, coarse-grained (CG) counterparts have paved the way to efficiently study and simulate largescale processes at time scales inaccessible to all-atom models [14,15]. To this end, it is beneficial to consider simplified * tapio.ala-nissila@aalto.fi geometries, such as assuming rod-like, rigid PEs, axially symmetric charge distributions, and simplified descriptions of ions in the solution around the charged species.
The mean-field Poisson-Boltzmann (PB) theory has proven to be efficient in describing the condensation of monovalent ions around single low-charge rods [16][17][18][19][20][21]. This theory is able to predict various properties of PE solutions including electrophoretic migration, surface adsorption and osmotic pressure, for a wide range of concentrations [22][23][24]. However, PB theory fails when charge correlations become strong, such as for high ion valency, high surface-charge density, and low temperatures. There exists an extensive literature on the effects of charge correlations on ion condensation for PEs beyond the mean-field PB theory [25][26][27][28]. Perhaps the most striking effect of such correlations is the reversal of the effective PE charge for multivalent counterions, which also reverses the field-driven PE mobility [26,29].
Another shortcoming of these simplified models is that PEs, such as DNA or charged polypeptides, are atomistically structured, exhibiting a local spatially inhomogeneous charge distribution. The details of the local structure are however important and affect both the persistence length [30,31], the charge distribution around the PE, and the PE interactions [32,33]. Ion condensation and salt solution response of PEs are also salt species and PE dependent [32,33]. Such dependencies can, to some degree, be captured by an empirical modifications of the PB theory [21,34].
In systems consisting of two rod-like PEs, PB theory has also been often applied, in particular in the case of monovalent salt solutions. The interaction between two charged cylinders has been calculated in many theoretical works with different boundary conditions and approximations [35][36][37][38]. The linear version of PB theory with the Debye-Hückel approximation (LPB) provides an explicit analytical solution of the interaction between two parallel cylinders [39][40][41]. However, LPB theory fails in dealing with highly charged molecules or small cylinder radii (as comparable to the Debye length), such as, e.g., in the case of DNA [35]. Charge correlations in the strong-coupling regime lead to like-charge attraction [42][43][44][45][46][47][48]. For dense PE systems, the attraction induces, e.g., bundle formation of F-actin and toroidal aggregates of concentrated DNA [12,49].
Recently, we have developed a soft-potential-enhanced PB (SPB) theory to efficiently and accurately compute ion density distributions in the weak-coupling regime [50]. The soft potential in the SPB theory contains only one (physical) parameter, which can be fixed to give a good description for monovalent ion densities and electrostatic potentials around single rodlike PE for a wide range of salt and ion sizes. The SPB approach was shown to work well for monovalent salt concentrations up to 1 M and ion sizes ranging from those corresponding to small, hard monovalent ions, such as Na + and Cl − , to almost an order of magnitude larger ions with diameter comparable to the PE diameter [50].
Encouraged by the success of the SPB theory, in the present work we apply it to compute the ion distributions and corresponding PE-PE interactions via the potential-of-mean-force calculation for two rodlike PEs in the weak-coupling regime. We compare the SPB modelling outcome against a chemically specific PE. To this end, we choose a relatively highcharge-per-length cationic and very common polyelectrolyte poly(diallyldimethylammonium) (PDADMA) as the PE focus of the study. PDADMA is chosen because of its technological relevance in, e.g., water purification, flocculation, and paper industry [51]. Furthermore, the strong amide charges can be expected to cause localization, i.e. deviations from the mean field predictions. The SPB theory is benchmarked against molecular dynamics (MD) simulations of a CG model of PDADMA PE in salt solution, which allows the SPB fitting parameter to be optimized and fixed. We obtain good agreement between the mean-field and coarse-grained models in capturing the interactions between two PEs. To further elaborate the influence of the atomic structure of the PE, we show by atomistic-detail MD simulations of PDADMA how positional dependencies between the charged groups induce correlation in the ion distribution. For the simulation setup in which the atomistic-detail PE is stretched straight, these emerge as asymmetric, helical configurations. Such positional correlations in the ion distribution lead to small deviations from the SPB and CG model predictions in interaction strength and show as short-range deviation from the PB smooth prediction in the potential of the mean force.
The manuscript is organized as follows: Section II introduces the different methodologies used in this work. Specifically, in II-A we describe our all-atom MD simulations, in II-B the CG-MD simulations, and in II-C the SPB theory. In Section III-A we present the results regarding the SPB model parametrized against CG-MD simulations and compare the different models against the all-atom MD simulations. We discuss the ion density distributions predicted by the different approaches, underpinning the effect of the atomistic structure of the PE. Furthermore, in Section III-B we show good agreement between PE-PE interactions, as captured by the potential of mean force between two parallel PEs corresponding to PDADMA in monovalent salt solution obtained from the different description levels. Finally, in Section IV we summarize our findings and provide prospects for the application of SPB to model chemically specific complex PE systems.

A. Atomistic molecular dynamics simulations
The system setup is shown in Figs. 1(a) and (b) where the atomistic-detail MD simulations include either one or two linear PE chains set to span the cuboid simulation box along the z axis as periodic chains (terminal group connected over the periodic boundary condition). The initial configuration preparation protocol follows Ref. 21 and leads to infinite, straight PE chains (covalently bonded over the periodic boundary). The stretched PE configuration is used as a simplification, enabling direct comparison with PB theory for cylindrical objects. In comparison, a free-ended PDADMA chain fluctuates in backbone conformations, which influences both ion distribution and complexation with oppositely charged PEs, see e.g. Refs. 52-55. PDADMA is chosen due to the relatively high charge per unit length and the charged nitrogen residing away from the backbone axis which induces some spatial asymmetry. The PE is fully charged (all monomers have dissociated electrolyte groups) and Cl − ions are considered as the counterions to neutralize the polyelectrolyte charge. For excess salt, NaCl as dissociated ions is introduced. PDADMA is a relatively symmetric molecule in its charge distribution. However, when set as an axially straight, periodic chain, the charged nitrogens (N + ) adopt a helical configuration around the PE main axis as their equilibrium configuration. This results from the charge groups repelling each other and steric effects. Additionally, methyl groups shield the nitrogens. PDADMA structure is shown in Fig. 1(e).
The atomistic simulations were performed using the GRO-MACS 2019.6 package [56,57]. OPLS-aa force field [58] was used to describe the PE, whereas the explicit TIP4P model [59] was employed for water molecules. For the Na + and Cl − ions in the simulations the parameters originate from Refs. 60 and 61, respectively.
The GROMACS solvate tool is used to solvate the PE chains. Excess salt (NaCl) in the concentration range 0.13 − 1 M is introduced by replacing random water molecules by the ions. Finally, the equilibrated simulation box has dimensions L x = L y = 7.9 nm and L z = 5.66 nm. The equilibrated straight PE chain length dictates L z . For PDADMA, 10 monomers long chain is used based on our earlier work [21] mapping this to be a sufficient length for the ion distribution to converge.
Van der Waals interactions are modelled using the Lennard-Jones potential. The particle-mesh Ewald (PME) method [62] was applied for the long-range electrostatic interactions with a 0.16 nm grid spacing and fourth-order spline. Both the van der Waals and real-space electrostatics employ a 1.0 nm direct cutoff (no shift). All the bonds in the PE and water molecules were controlled by LINCS [63] and SETTLE [64] algorithms, respectively. The stochastic V -rescale thermostat [65] with a coupling constant of 0.1 ps and reference temperature T = 300 K is used for temperature control. On the other hand, pressure control is semi-isotropic with Parrinello-Rahman barostat [66,67] alongside a coupling constant of 1 ps and reference pressure 1 bar. Following Ref. 21, the system is set to be incompressible along the z axis. A 2 fs integration time step within the leap-frog scheme was applied in the N pT simulations.
The initial configuration was first energy minimized via steepest descent algorithm until the largest force in the system was smaller than 1000 kJ mol −1 nm −1 . Next, a 2 ns semiisotropic N pT simulation in which the positions of the heavy (i.e. non-hydrogen) PE atoms were kept fixed was performed to adjust the x and the y dimensions of the system as well as the distribution of water and ions around the PE. Finally, the PE positional constraints were released and a 100 ns semiisotropic N pT simulation of a single PE was performed, of which the first 2 ns were disregarded from the analysis.
In order to examine the interaction between two PEs, a second PE was also placed parallel to the z axis into an identical in size simulation box. The initial center-to-center distance of the z-axial PEs was set to be d = 3.5 nm (cf. Fig. 1(b)). Solvation and initialization of the system followed the same procedure reported for the single-PE case.
An estimate of the potential of mean force (PMF) for the two PE interaction was obtained by umbrella sampling. To generate the umbrella sampling configurations, the two PEs were pulled together at a fixed rate of 10 nm/ns using a stiff spring (k = 3200 k B T /nm 2 ). We sampled 21 configurations uniformly distributed at distances d between 3.5 − 0.5 nm. For each umbrella sampling window, a 20 ns simulation with a harmonic tether potential with a spring constant of k = 320 k B T /nm 2 , but otherwise following the semi-isotropic N pT setup described above, were performed. The PMF was extracted using the Weighted Histogram Analysis Method (WHAM) of GROMACS [68]. The VMD package [69] was used for molecular visualizations.

B. Coarse-grained molecular dynamics simulations
Following Ref. 50, the PE is modeled by a linear series of CG spherical beads. The CG PE is confined along the z axis of a simulation box of size 20 × 20 × 20 nm 3 with implicit solvent and periodic boundary conditions in all directions. The initial configuration is shown in Fig. 1(c) using the VMD software package [69]. Each bead carries a charge of ζ e where e the elementary charge. This is equivalent to a line charge density λ = ζ e/b where b is the bead separation distance. In total, 74 consecutive beads, each with charge ζ e = 0.473e, are set at a distance b = 0.27 nm from one another. This leads to a line charge density λ = ζ e/b = 1.75 e/nm, matching with the atomistic-detail PDADMA. The ions were set initially at random positions in the simulation box but avoiding overlapping using the Packmol package [70].
The pair interactions between all the particles in the system are modelled via a standard Weeks-Chandlers-Andersen (WCA) [71] potential of the form The labels i, j denote either the CG ionic species (Na + , Cl − CG equivalents) or the polymer beads (B), and r is the distance between the pair i j; σ i and ε i denote the diameter and the depth of the potential well for species i. We use the Lorentz-Berthelot (LB) mixing rule σ i j = (σ i + σ j )/2 and ε i j = √ ε i ε j . r i j c = 2 1/6 σ i j is the cutoff radii of pair interactions. We set ε B = 0.1 kcal/mol, similar to the value of the ionic counterparts ε Na = 0.13 kcal/mol and ε Cl = 0.124 kcal/mol [72][73][74]. The diameter of the polymer beads σ B = 0.54 nm is determined by minimising the root-mean-squared error (RMSE) between ρ Cl and ρ Cl SPB (see Supplementary Material (SM)-1). The diameters of the CG ions σ Na = 0.234 nm and σ Cl = 0.378 nm are taken from Refs. 72, 74. The PE charge of 35e is neutralized by 35 monovalent counter-anions. In analogy to atomistic-detail simulations, we examine the system response to added monovalent salt concentrations of 0.13 M, 0.26 M, 0.52 M, and 1 M.
The electrostatic interactions were modeled via Coulombic potentials φ i j C . The pairwise interaction between two ionic species i and j with charges ζ i e and ζ j e is given by where β = 1/k B T and the Bjerrum length l B = β e 2 /(4πε) is the distance at which two unit charges have interaction energy on the order of k B T [18,75]. In water the permittivity is ε = 78ε 0 [76], where ε 0 is the vacuum permittivity. The implicit solvent is modeled by a continuum dielectric medium. All CG MD simulations employ the LAMMPS Jan2020 package [77,78]. The long-range electrostatic interactions were calculated using the particle-particle particle-mesh method (PPPM) [79]. Coulombic pairwise interactions were calculated in real space up to a cut-off of r i c = 13σ i , where i is any charged species, beyond which the contributions were obtained in reciprocal space. All simulations were performed in the NV T ensemble with a PPPM accuracy of 10 −5 . The reference temperature was set to 300K, controlled by the Nose-Hoover thermostat [80,81].
Once the initial configuration was prepared, we performed an energy minimization of the system. After this, a 0.2 ns NV T simulation in which the equations of motion were integrated with the velocity Verlet algorithm and 1 fs time step was carried out as initial equilibration of the system. Finally, a 20 ns NV T simulation was run for data analysis, out of which the first 1 ns is disregarded. Here, a 2 fs time step was used.
In the two-PE case, both PEs are also parallel to the z axis. The first PE is set at the center of the box with (x, y) = (0, 0), while the second PE is set at (x, y) = (d, 0), as shown in Fig. 1(d). With the exception of 70 CG Cl ions for neutralization, all initial settings for the two-PEs case are the same as in the single-PE case.
The equilibration and production run are performed analogous to the single PE case. For determining the PMF between the two CG PEs, 15 two-PE configurations with PE-PE axial separation distance d between 1 nm and 5 nm at equal intervals were generated. The PEs were kept fixed in position and for each fixed separation distance an MD run of 20 ns was performed. The PMF vs distance d was calculated by integrating the mean force f CG over the separation distance as Because of its success in describing monovalent ionic environments, the Poisson-Boltzmann (PB) equation is often used to describe the ion distribution and electrostatic potential in electrolyte solutions. The PE is modelled as an impenetrable and rigid cylinder with surface charge density λ /2πa, where a is the cylinder radius. We consider here the case of positively charged PEs which have to be neutralized by adding negative counterions with number density ρ ci . After this, salt is added to the system. Thus, the number density of cations (ρ Na 0 ) equals that of the added salt ρ 0 . The number density of anions (ρ Cl 0 ) then equals ρ 0 +ρ ci . The full Poisson-Boltzmann (PB) equation for the electrostatic potential φ PB (r) surrounding such a cylinder can be written as where i = Na, Cl, ζ i is the ion valency and ρ i 0 the number density of the i th ion species. The ion number densities from the PB theory can be obtained from In some cases we will also use data obtained from the linearized version of the Poisson-Boltzmann (LPB) equation, for which a simple closed-form analytic solution exists for two cylindrical polymers as explained in SM-2.
In a recent work [50], the prediction of ion number concentration from PB theory was enhanced with a cylindrical soft-potential correction (SPB) of the form which replaces the impenetrable cylinder in the standard PB theory with a soft cylinder potentialṼ i (r). The effective po-tentialṼ i (r) felt by the ions i is defined as a cylindrically symmetric WCA potential, i.ẽ where the cutoff radius is defined byr Bi c = 2 1/6σ Bi . Similar to Ref. 50, we useσ Bi = σ Bi ,ε BNa = 0.107 kcal/mol andε BCl = 0.11 kcal/mol, respectively.
To obtain accurate results from the SPB theory, the corresponding cylinder radius a SPB has to be adjusted in Eq. (5) for the electrostatic potential φ PB (r). This effective radius may change with system conditions, such as, e.g., ion concentration. The optimization is done by comparing Cl − ion number densities from the CG models with number density of the SPB theory prediction (Eq. (5)). From different values of a, Eq. (3) allows us to calculate φ PB (r), followed by Eq. (5) to calculate ρ Cl SPB . We identify the optimal radii a * SPB giving rise to the optimal potentials φ * SPB (r) in Eq. (5) based on minimizing the RMSE between ρ Cl CG and ρ Cl SPB . Additional information regarding soft-potential enhanced linear PB (SLPB) theory is provided in SM-2. The SPB optimized radii (as well as those from SLPB) at different salt concentrations are summarized in Table S1 of the SM. It is interesting to note that for different salt concentrations, the radius of the SPB theory only changes marginally (≈ ±13%). We thus use the average value of a * SPB = 0.44 nm throughout.

A. Single-PE case
We first assess the accuracy of the CG-MD simulations and the SPB and SLPB theories with respect to ion condensation around a single PE. We start by setting up all-atom simulations of a single PDADMA at added monovalent salt concentrations ρ 0 = 0.13 M, 0.26 M, 0.52 M and 1.0 M, respectively. To calculate the radially symmetric (perpendicular to the z axis) ion number density profiles ρ i (r), the x and y coordinates of the center of mass (CM) of the PE have been taken as the reference point. The density profiles ρ Na (r) and ρ Cl (r) at different salt concentrations are shown in Fig. 2. ρ Cl obtained from allatom MD simulations exhibit two peaks at r ≈ 0.45 nm and r ≈ 0.7 nm, respectively. We note that the second peak increases as a function of salt concentration. These peaks are a result of the atomic level structure of PDADMA. Notably, the charged N + that are usually found at roughly 0.25 nm from the center of the PE are responsible for the enhanced localization of the Cl ions. Details of the atomic configurations are presented in SM-3. To assess sampling equilibration, the orientational self-correlation function is calculated in SM-4.
In the CG model and in SPB theory, we modelled PDADMA as a cylinder with effective radius σ Bi and a * SPB , respectively. These are the only parameters to be optimized, and are obtained by minimizing the RMSE between the different models. Specifically, σ BCl = 0.46 nm is obtained from averaging the minimised RMSE between all-atom MD and CG-MD simulations at different added salt concentrations (see the CG-MD simulations part of Methods). On the other hand, the value of a * SPB = 0.44 nm is obtained by averaging the minimised RMSE between CG-MD simulations and SPB predictions at different added salt concentrations, as reported in the SM-1. Note that the latter value is also used for SLPB. The comparison between SPB and SLPB theories is shown in SM-5. In Fig. 2 we report the ion number density profiles from all the different models addressed here. We can see a satisfactory agreement between the different curves, with the exception of the small-scale oscillatory structure in the all-atom MD simulation results, mainly due to local asymmetry of PDADMA which cannot be captured by any of the CG techniques under the current assumptions.

B. Two-PE case
In the present case where the system has monovalent ions and is in the weak-coupling regime, there is no charge inversion and the two PDADMA molecules will always repel each other. We quantified the interaction between two PEs via numerical computations of the PMF. Our results from all the different models and different salts here are summarized in Fig. 3, with the all-atom MD being the ultimate benchmark for comparison. The first general observation is that the range of the PMF curves decreases with increasing salt concentration. This is expected and is due to increasing screening effect mediated by the condensing anions. We find good agreement between most of the methods, with the exception of smallamplitude oscillations at d ≈ 1.2 nm that only show up in the all-atom MD data. This feature, as well as the differences ris- ing from the different models, are addressed in more detail in the following sections.

The CG model and SPB theory
In the case of the CG-MD model, the value σ B = 0.54 nm is inherited from the single-PDADMA case. The CG model qualitatively captures the PMF curves when benchmarked against all-atom MD simulations at different salt concentrations up to 1 M (see Fig. 3). We find that due to the second PE, the condensation peak maximum of the Cl ions corresponds to 3 − 4 times higher number density than that of the single PE case at short distances. The detailed results on this are shown in Fig. S3 and S4 of the SM-6.
In order to obtain φ SPB (abbreviated with φ in this section), we employ a finite-element package (COMSOL5.2) to solve Eq. (3) in the xy plane with the von Neumann boundary conditions at the surfaces of the two disks, i.e. ∇φ ·n| x 2 +y 2 =a 2 = λ /2πaε and ∇φ ·n| (x−d) 2 +y 2 =a 2 = λ /2πaε. To do so, we use a = a * SPB = 0.44 nm as obtained in the single-PE case. The same is applied to the SLPB theory. Heren is the unit normal vector of each surface, as sketched in Fig. 1(g). The mean electrostatic force between the two cylinders per unit length f (d) at distance d was calculated via the stress integral [35] f (d) = ε where κ = e 2 (ρ Na 0 + ρ Cl 0 )/εk B T is the screening parameter. The mean force can be integrated to get the potential of mean force between two cylinders as V PMF (d) = ∞ d f (x )dx . The SPB prediction of the PMF is in good agreement with the CG model and the all-atom MD simulations, as shown in Fig. 3. In contrast, the PMFs from linear SPB are inaccurate at low salt concentrations. The linear SPB predictions are discussed in more depth in SM-5.

The influence of the atomistic structure
When the two atomistic detail PDADMA chains in the allatom MD simulations are separated by a distance comparable to the size of a monomer, the orientations of the PE side chains will affect their interaction. As shown in Fig. 1(e), the 10 charged N + atoms of each PDADMA chain form in the elongated, axially straight configuration a right-handed cylindrical helix. In cylindrical coordinates, with the longitudinal axis centered at the backbone of each PDADMA chain (see the two black dashed lines in left panel of Fig. 4(a)), the N + atom of the l th monomer of the k th PE has coordinates (r l k , θ l k , z l k ), with k = 1, 2 and l = 1, 2, ..., 10. We take the positive x direction to correspond to θ = 0. We set a plane at z = 0, where the orientations of the two helices are defined by the angles θ k (k = 1, 2), respectively, as shown in the right panel of Fig. 4(a). Because the pitch of the helices is L z /2, the angles θ k can be calculated from the location of the first monomer (i.e. l = 1) as θ k = θ 1 k − 4πz 1 k /L z . We then define the variable ξ ≡ |θ 1 − θ 2 | to represent the relative orientation between two PDADMA chains. In Fig. 4(b) we show the average value of ξ /π as a function of the distance d between the two PEs. At long distances, the positioning of the PDADMA charges remains uniformly random with ξ /π = 1/2. In contrast, at small PE-PE separations of d < 1 nm, the charge group positioning is highly correlated which shows as two orientations coupling strongly. The coupling rises from minimization of electrostatic and steric repulsion by rotation and axial displacement of the PDADMA chains. Nevertheless, the outcome is a strong interchain repulsion evident in the corresponding PMFs.

IV. SUMMARY AND CONCLUSIONS
In this work we have addressed the interactions between two PEs by using the recently developed soft-potential-enhanced Poisson-Boltzmann (SPB) theory. We have here focused on the case of PDADMA molecules in monovalent salt solutions in the weak-coupling regime, where the interactions are expected to be fully repulsive at all distances. To this end, we have first addressed the ion number density around a single PDADMA. In the elongated, axially straight PE configuration, the charged N + atoms of PDADMA adopt a helical configuration as their minimum energy configuration by electrostatics and sterics. This leads to the two peaks in the ion distribution around the atomistically detailed PE (see Fig. 2). Although the isotropic SPB theory cannot reproduce such features, it can qualitatively capture the atomistic ion distributions at salt concentrations up to 1 M. By minimizing the RMSE between the number density distribution of the CG model which has been fitted on all-atom MD simulations and the one obtained from SPB theory, we find an optimal effective radius for PDADMA (a * SPB = 0.44 nm), which is the only (physical) fitting parameter in the SPB theory.
To study the interactions between two PEs we have considered two parallel rodlike PDADMA chains and numerically computed the PMF curves at different salt concentrations. We find that the predictions from the SPB method are accurate when the distance between the two PDADMAs is larger than their radius. There are small oscillations in the PMFs at short distances as revealed by all-atom MD simulations. These oscillations, rising from the chemical structure of the PE, are due to the coupling between charged atoms, steric packing of the side chains, and correlations in the relative orientations of the chains at distances smaller than about 1.3 nm.
While our work suggests that the SPB theory is a simple and accurate way to model interactions in complex PE systems for a large range of system parameters, the limitations of the theory should be kept in mind. We naturally expect the theory to fail when the approximations in the weak-coupling theory are no longer valid. This would happen in cases where the PE line charge density is high, there are electrostatic correlations arising from strong Coulombic interactions, or the assumption of a straightened PE as in our simulation setup leads to a large entropy loss that destabilizes such configurations (see SM-7 for some additional discussion).

Supplementary Material 1. ROOT-MEAN-SQUARE ERROR ANALYSIS
We optimize the parameters σ B and a * SPB by minimizing the root-mean-square error (RMSE) between the Cl − density profile obtained from the different models. The definition of RMSE between two density profiles ρ 1 (r) and ρ 2 (r) is defined as: where N is the number of sampling points. Specifically, in order to obtain the effective diameter of the PE in the CG model, σ B , ρ 1 is replaced by ρ Cl from all-atom MD simulations, whereas ρ 2 by ρ Cl CG , obtained from the CG model, which depends on σ B . The values of the optimized diameter σ B for different salt concentrations are shown in Table. SI. We choose the average over the four salt concentrations, i.e. σ B = 0.54 nm.
On the other hand, in the case of SPB theory, ρ 2 is replaced by ρ Cl SPB , which depends on a SPB . The optimal values from fitting a * SPB at different salt concentrations are shown in Table  SI. Also in this case we choose the average over the different salt concentrations, i.e. a * SPB = 0.44 nm, used throughout the work (also for the SLPB theory).

SOFT-POTENTIAL-ENHANCED LINEAR POISSON-BOLTZMANN THEORY
When β eφ (r) 1, which corresponds to small electrostatic potentials, we can linearize the full PB equation. This linear approximation is commonly referred to as the Debye-Hückel approximation. Equation (3) in the main text then becomes where 1/κ is the Debye or screening length [28]. The analytic solution of Eq. (S2) for r > a is [41] φ LPB (r) = λ 2πaκε where K 0 and K 1 are modified Bessel functions of order zero and one, correspondingly. The soft-potential enhancement can be used for both the full and linearized PB theories. Finally, for the two-PE case, the corresponding force from the linear theory f SLPB and the PMF are obtained by numerically solving Eq. (S2) similarly to the full PB equation, with the boundary conditions at the surface of the two disks with radii a = a SLPB (cf. Fig. 1(g) in the main text) ∇φ ·n| x 2 +y 2 =a 2 = λ /2πaε and ∇φ ·n| (x−d) 2 +y 2 =a 2 = λ /2πaε.

THE ATOMIC STRUCTURE OF PDADMA
Our atomistic PDADMA MD model consists of 10 monomers and the polymer axis is set along the z direction. Each monomer corresponds to a backbone length of L z /10. Since the charged groups of the monomers adopt a helixlike structure (see Fig. 4(a) in the main text) for the straightened PDADMA chain, we define the location of N + of the l th monomer by (r l , θ l , z l ), l = 1, 2, ..., 10 in cylindrical coordinates with longitudinal axis parallel to (and centered at) the backbone of the PE. The ion number density around the l th monomer (in the region z l − L z /20 < z < z l + L z /20) is projected into a two-dimensional density map ρ Cl,l (r, θ ) in polar coordinates. The ion number densities around different monomers are similar but the angular orientation of the density follows the relative orientation of the monomer θ l . To account for this, we rotate the ion number density map by −θ l degrees as ρ Cl,l (r, θ − θ l ), to overlay the orientations of the monomers such that superposed ion number density data in terms of orientation with respect to monomer is obtained. Then we take an average over ρ Cl,l to get the density ρ Cl = ρ Cl,l (r, θ − θ l ) l for different salt concentrations as shown in Fig. S1(b). This ion number density distribution around the PE monomers concretely shows the origins of the two-peak profile in the ion distribution in the main manuscript: ions form a semi-circular ring around the charged group with ion condensation localization visible both at the sides and at the tip of the charged group, leading to the binodal ion number density distribution. Similar findings were already reported in Ref. 21. In particular, the Cl − ions stack mainly on the directions ±π/2 and 0, corresponding to the first and the second peaks of Cl − in the radial density profiles, respectively. This configuration is caused by two main interactions on Cl − ions: attraction from positively charged group centered at N + atoms and steric repulsion from the surrounding neutral atoms. The positionally correlated configurations of the PE charge groups and its effect on the ion density span over different salt concentrations.

THE ORIENTATIONAL SELF-CORRELATION FUNCTION
The atomistic-detail PDADMA molecule spans the z axial direction of the simulation box as a periodic molecule. However, it can both translate and rotate along the z axis during the MD simulations. The orientation of the chain on the z = 0 plane is given by the angle θ 1 as defined in Fig. S1 (l = 1). During a 100 ns simulation we measured the time dependence of θ 1 (t) and its time averageθ 1 = θ 1 (t) t . We then define the normalized orientational (fluctuation) self-correlation function C(t) as where ∆θ 1 = θ 1 −θ 1 . As an example, for a single PDADMA chain at salt concentration ρ 0 = 0.26 M, the self-correlation function decreases exponentially as exp (−t/t c ) as shown in Fig. S2. From this data we can estimate the self-rotational correlation time as t c ≈ 15 ns. We addressed the differences in the electrostatic potential φ SPB (r) and φ SLPB (r) with the same linear charge density λ , effective radii a SLPB = a SPB = a and monovalent salt concentration ρ 0 . A theoretical comparison between the two approximations is reported in Ref. 82. However, here we focus on PDADMA. The difference is expressed as the relative error ∆ φ ≡ |φ SLPB (r) − φ SPB (r)| /φ SPB (r) r , which is averaged over a < r < L x /2. Using the Debye length κ −1 as a unit of length, the only characteristic length scale in the system is κa. Thus, ∆ φ solely depends on λ and κa. We show this FIG. S5. The Cl − and Na + number density profiles along the x axis connecting the centers of the two PEs from the CG-MD model (as sketched in the inset). The center-to-center distance d changes from 1.1 nm to 3.0 nm. in Fig. S3. As expected, for large values of λ , small values of a and low salt concentration ρ 0 (corresponding to low κ), we find the highest-error region (red). We can see that in the case of PDADMA (circles), the PB and LPB are almost identical at all salt concentrations considered here. This is due to a relatively small linear charge density λ = 1.767 e/nm, giving rise to small electrostatic potential φ . We conclude that in the case of PDADMA the linear approximation is successful in capturing the ion distribution.
For the two-PDADMA case, the relative error of mean force and electrostatic potential at x = d/2 are expressed as ∆ f ≡ | f SLPB − f SPB | / f SPB and ∆ φ ≡ |φ SLPB (d/2, y) − φ SPB (d/2, y)| /φ SPB (d/2, y) y , respectively; the latter is averaged over −L y /2 < y < L y /2. Their comparison is shown in Fig. S4. In contrast to the single PE case, the mean force from SLPB is very different from that of the SPB, especially when the two PEs are close.

ION DENSITY AROUND TWO CG-MD PDADMAS
In the CG model, we can examine the ion number densities around two PEs. As shown in Fig. S6, a significant fraction of the Cl − ions accumulate in the region between the two PEs when the distance between the two rods is relatively short (d ≤ 1.5 nm). In Fig. S6 we report the two-dimensional number density for d = 1.2 nm and d = 2.0 nm.
To better quantify this, we obtained the one-dimensional cross section of the ion-number density profiles between the two PEs at 0 < x < d, y = 0 (see the inset of Fig. S5) for different values of d, and compared them to the single-PE case, as shown in Fig. S5. When the two rods are close (d = 1.1 nm), the peak value of the Cl − ion-number density is 3 to 4 times larger than that of the single-PE case. On the other hand, the Na + ion number density is much lower than that of a single rod. This is due to the stronger electrostatic repulsion, as well as the Cl − ions occupying more space in between the rods.

LIMITATION OF THE SPB APPROACH
In the original work of Vahid et al. [50] where the SPB theory was introduced it was applied to model monovalent ion distributions around single polystyrene sulphonate (PSS) molecules. The corresponding ion densities were accurately reproduced for a wide range of system parameters, including salt and ion sizes. We performed additional testing for the case of two such PEs with atomistic-level MD simulations and the same setup as for the PDADMA molecules here. However, we found that there was an attractive force induced between two PSS molecules at short distances. Such attraction is not expected in the weak-coupling regime and here it may be caused by entropy loss due to the axially straight, infinite PE chain