Chemically realistic coarse-grained models for polyelectrolyte solutions

Polymers bearing ionizable groups have numerous important technological applications, including ion-exchange membranes, drug delivery vehicles, and viscosity modifiers.^1–3 They are also of scientific interest. The mixing of polyions of opposite charge in the presence of salt results in an interesting phenomenon called complex coacervation,^4,5 which is a phase separation into a polyion rich and a polyion poor phase and is of fundamental and practical importance.^6–8 

A molecular understanding of polyelectrolytes in solutions has been the focus of active research. Theoretical and computational work has focused primarily on bead-spring models for polyelectrolyte solutions.^9–19 While these studies shed light on the qualitative behavior of charged polymers, they lack the chemical detail that will allow for a direct comparison with experiment, beyond scaling laws. Another deficiency of many of these models is that they use an implicit solvent. Chang and Yethiraj²⁰ found that an incorporation of explicit molecular solvent made qualitative differences to the polyelectrolyte behavior, as has also been emphasized in more recent investigations.^20–22 Atomistic models can provide a molecular level understanding but are computationally intensive. Coarse-grained (CG) models, where many atoms are grouped into a single site, are, therefore, an attractive middle ground for the computational and theoretical student of polyelectrolyte solutions.²³ Early CG models for polyelectrolytes,^24,25 however, did not incorporate the solvent explicitly.

The most widely studied polyelectrolyte, from a polymer physics perspective, is sodium polystyrene sulfonate (NaPSS). A number of experimental studies have focused on the behavior of NaPSS in semi-dilute solutions (salt-free and with added salt).^26–32 The excess scattering observed in the x-ray scattering studies of NaPSS semi-dilute solutions point toward the formation of aggregates in solution. SAXS (small angle x-ray scattering) and SANS (small angle neutron scattering) studies show that the addition of a good solvent (such as THF) in a partially sulfonated NaPSS solution causes the low wavevector scattering peak to disappear.^30,32 The influence of charge density on the structure and dynamics of NaPSS in aqueous solutions has been investigated.^28,29 Neutron scattering measurements revealed the shift of the high wavevector polyelectrolyte peak toward smaller wavevectors in the presence of divalent counterions.²⁸ Dynamics light scattering experiment of NaPSS in the presence of different concentrations of salt revealed the collapse from a coil to a globule.³³ 

There have been theoretical ideas and conjectures regarding these and other phenomena, but no molecular level understanding from a computational perspective. The development of chemically realistic coarse-grained models, therefore, is of interest and significance. Note also that water is a poor solvent for the backbone of PSS, whereas most previous theoretical and computational work considered good solvents for the backbone.

There are several questions CG models can address. What are the ion correlations in the presence of an explicit solvent? Is there ion-pairing or solvent separated ion pairing? Is there counterion condensation and are there long lived dipoles in the system? How do the dynamics differ from neutral polymer solutions? What is the balance between electrostatic and non-electrostatic correlations and how does this depend on the chain length? In this work, we develop chemically realistic CG models for polyanion and polycation.

A popular and successful approach to the development of coarse-grained models in biophysics^34,35 and soft matter^36,37 is the MARTINI framework. Mantha and Yethiraj³⁸ developed a CG model for NaPSS with two models for the explicit water, denoted by BMW and POL. They found that the POL model could be parameterized to reproduce the end-to-end distance potential of mean force (PMF) found in atomistic simulations. The model also reproduced the solvation energy of monomers in addition to the conformational properties of the polymer. The simulations showed that NaPSS collapsed into a cylindrical rather than spherical shape and that the conformational properties were sensitive to the sequence of the sulfonated (charged) sites on the polymer. Another CG model for NaPSS was proposed by Vögele et al.³⁷ to study polyelectrolyte complexes. The two aforementioned models are quite different in the choice of topology and partial charge and have never been previously compared.

In this work, we develop a CG model for polystyrene sulfonate (PSS) and polyvinyl(benzyltrimethyl) ammonium (PVBTMA) along with their respective counterions, i.e., Na and Cl (see Fig. 1). From a computational standpoint, the only difference in the chemical structure of monomeric unit of these two polyions is the charged pendant group. This is, therefore, a good model polycation–polyanion system in order to understand the differences in overall charge and complexation without additional complications due to chemical, e.g., hydorophobicity, differences. From a model building standpoint, most of the groups are similar, and therefore, few force field parameters have to be fit for PVBTMA, once a model for PSS is available. Following previous work,^37,38 we use the POL model for water.

We use molecular dynamics (MD) simulations and enhanced sampling methods to investigate the conformational behavior, static structure, and dynamics of dilute polyelectrolyte solutions. A previous simulation study on atomistic models for NaPSS has shown that simple molecular dynamics does not adequately sample the configurational space.³⁹ Park et al.³⁹ used Hamiltonian replica exchange molecular dynamics simulations. Other enhanced sampling methods have been used effectively to calculate the properties of CG models.^38,40 We used the most effective enhanced sampling method available to study the properties of atomistic and CG models. We find that specific choices in the coarse-graining methodology, such as the partial charge distribution and topology, make a dramatic difference to the observed properties. The polycation and polyanion, however, display similar properties in dilute salt-free solutions.

The rest of this paper is organized as follows: The simulation methodology is presented in Sec. II, the results are presented in Sec. III, and some conclusions are presented in Sec. IV.

The Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field^41,42 with extended simple point charge (SPC/E) water is used for the atomistic model. The polarizable MARTINI (POL)⁴³ water model is used in the CG simulations. The other CG bead types are similar to what has been used previously, and the interaction parameters are developed in this work. A pictorial representation of a single 16-mer chain of PSS and PVBTMA is depicted in Fig. 1.

The simulation cell is a cube with periodic boundary conditions in all directions. The linear simulation cell dimension (L) is taken to be very large as compared to the size of the polyion, i.e., R_g $<$ L to make sure that the system is in the dilute regime. For the largest chain size studied, R_g = 3.56 nm and L = 15.5 nm. All the systems are, therefore, in the dilute regime, i.e., the concentration c = N/L³ is smaller than the overlap threshold concentration $c*=N/Rg3$⁠.

Initial configurations for atomistic simulations are generated by placing a single polyion in a cubic simulation cell. The cell is solvated with water and then the counterions are added. The energy is minimized using a steepest decent algorithm, followed by two equilibration steps of 20 ns each in the isobaric–isothermal (number of molecules, pressure, P, and temperature, T, fixed) ensemble at 1 bar pressure and 350 and 300 K, respectively. The temperature and pressure at these steps are maintained using a Berendsen thermostat and barostat.^44,45 The results are averaged over a production run carried out for 300 ns at 300 K and 1 bar with the temperature and pressure maintained using a Nose–Hoover thermostat and a Parrinello–Rahman barostat.^46,47 Configurations are saved at every 10 ps for the calculation of properties. We verify that the radius of gyration and end to end vector autocorrelation functions decay during the course of the simulation.

Non-bonded interactions are calculated as follows: A cutoff of 1.2 nm is used for nonbonded Lennard-Jones interactions, and Coulomb interactions are handled with the particle mesh Ewald (PME) method.⁴⁸ The PME parameters are as follows: real space cutoff distance of 1.2 nm and interpolation order of 6 with a maximum fast Fourier transform grid spacing of 0.12 nm. Bond lengths in the solute molecule are constrained using the linear constraint solver (LINCS) algorithm and the water molecules are kept rigid using the SETTLE algorithm. All the atomistic simulations are performed using a leap frog integrator with an integration time step of 1 fs.

The protocol for the CG simulations is similar to that of the atomistic simulations and the differences are noted below. The time step is 20 fs. The two equilibration steps are for 100 ns at the higher temperature and 300 ns at the lower temperature. The production run is for 500 ns in the canonical ensemble.

All simulations are performed using the graphics processing unit (GPU) variant of GROMACS 2019.4.^49–51 

We calculate the potential of mean force (PMF) as a function of collective variables (CVs) using enhanced sampling methods. The two collective variables we consider are the root-mean-square radius of gyration, R_g, and the root mean square end-to-end distance, R_e, of the beads along the chain backbone. These are defined as

and

where N is the degree of polymerization, $R⃗i$ is the position of the ith backbone particle, and $R⃗c$ is the position of the center of mass of the backbone particles.

The PMF is computed along the CVs with replica exchange with collective variable tempering.⁵² This method integrates the well-tempered metadynamics approach with a Hamiltonian replica-exchange scheme. In this approach, a ladder of replicas with increasing values of ΔT = (γ − 1)T, where γ is the bias factor and T is the temperature of the system, are provided to the system. These values range from 0 to a value large enough to cross all relevant energy barriers. The replica having ΔT = 0 and γ = 1 can be used to obtain the unbiased statistics; thus, this method does not require re-weighting of the results. Further details for this methodology are illustrated by Gil-Ley and Bussi.⁵² 

Enhanced sampling simulations are performed using GROMACS 2019 patched with PLUMED 2.5.⁵³ Four replicas for each collective variable are constructed using initial configurations equilibrated in the canonical ensemble. For the metadynamics parameters, we use a Gaussian bias with the height of 0.05 kJ/mol and the width of 0.05 nm. For each replica, a different bias factor (γ) is chosen varying from 1 to 10. A Nose–Hoover thermostat is used for maintaining the temperature at 300 K. The simulation parameters are the same as described in the molecular dynamics simulation method. The simulations are performed for 500 ns for the atomistic models and 1 µs for the CG models.

The important aspects in the development of a CG model are (a) the water model, (b) the topology of the polyion, (c) the charge distribution on the sites of the polyion, and (d) the non-bonded interactions. We incorporate water explicitly using the POL-MARTINI water model. The model maps four water molecules into a single site that has a polarizable charge distribution. The performance of the model with polyions has been previously tested³⁸ and shown to give good agreement with atomistic simulations.

A schematic representation of the topology is shown in Fig. 2. We adopt the mapping scheme proposed by Rossi et al.⁴⁰ for polystyrene. Each monomer (styrene sulfonate) group consists of five CG sites. The CG bead for the backbone is assigned as type SC1 (labeled BB) such that its center of mass is in between two consecutive styrene groups. SC4 is assigned for the ring bead (labeled RB), Qa for the sulfonate group (labeled AB), and Qd for the sodium counterion (the SC1, SC4, Qa, and Qd types are the standard MARTINI particle types.). The partial charge distribution on the CG sites is taken such that it agrees with the atomistic partial charge distribution, where the atoms of the phenyl ring carry partial negative charge. For the non-bonded interactions between the different CG sites, we adopt the parameters proposed by Mantha et al., which were developed to reproduce the interactions of monomers in water when compared to atomistic simulations.³⁸ The CG model developed in this work is not able to represent tacticity because the two backbone atoms are collapsed into a single CG site. A snapshot of equilibrated conformation of the PSS chain along with explicit counterions and water is shown in Fig. 3.

The CG model is in good agreement with atomistic simulations. Figure 4 compares the simulation results from the CG model above with atomistic simulations for the PMF as a function of the two CVs. Near the minimum in the PMF, which is the equilibrium value of the CV, the CG model is in good agreement with atomistic simulations for both R_g and R_e. At lower values of the CV, the atomistic simulations predict a lower PMF and there are signs of a metastable compact configuration for R_g ≈ 0.6 nm or R_e ≈ 1 nm. The CG model also predicts a metastable compact configuration, but the value of the PMF is much higher (⁠$>$20 kJ/mol). Figure S3 shows good agreement between the form factor from the CG model and atomistic simulations. The counterion–polymer pair correlation function is not in good agreement, as expected, because the size of the counterion in the CG model is large, and therefore, the first peak in this correlation function is missed (see Fig. S3). Previous CG models are not in agreement with atomistic simulations for the PMF. In addition, the results shown in Fig. 4 are from previous CG models of Vögele et al.³⁷ and Mantha and Yethiraj,³⁸ denoted by CG_mod1 and CG_mod2, respectively. The models are described in the supplementary material, and Fig. S1 compares the topology and charge distribution. The CG_mod1 model, which has a different partial charge distribution than the atomistic model, has a global minimum at a much smaller value of R_g or R_e than the atomistic simulations. The difference in charge distribution leads to more compact conformations, which leads to a collapse for longer chains. The CG_mod2 model, which has a different topology from the CG model but the same charge distribution, is in good agreement with atomistic simulations near the global minimum. However, as the concentration is increased, the chains collapse in this model, which is unphysical (a movie of this collapse is provided in the supplementary material). We, therefore, conclude that the previous models have shortcomings that are rectified in the model of this work. The potential of mean force is interesting because it also shows metastable local minima at higher values of the free energy. From a practical perspective, however, the probability distribution function is more intuitive. The probability distribution functions of R_g for all the models are shown in Fig. S2 of the supplementary material. We have also validated the CG model for higher N, which shows that the behavior of the model does not deteriorate for longer chains. Figure 5 compares the PMF as a function of the rms radius of gyration (R_g) for the 30-mer NaPSS system. The minimum in the free energy profile of atomistic simulation is in close agreement with that in CG simulation.

A CG model for PVBTMA is obtained following the protocol used for PSS. The topology of the chain is taken to be the same as in PSS, the charge distribution is obtained from the CHARMM force field for the atomistic model, and the non-bonded interactions are the same as in PSS, except for the positively charged group labeled TB in Fig. 6. There are different interactions defined within the MARTINI framework, and we assign the charged site as Q0 and adapt its non-bonded interactions with SC1 and SC4. We use the Qa particle type for the Cl ion. Figure 6 depicts the topology and charge distribution in the polycation model. A snapshot of the equilibrated configuration of a PVBTMA chain along with explicit counterions and water is shown in Fig. 7. Figure S4 and Table S1 in the supplementary material depict the structure and partial charges proposed by the CHARMM atomistic force field for a monomeric unit of PVBTMA. The positive charge seems to be localized on the side group attached to the phenyl ring; thus, here we assign the charge on the pendant bead (labeled TB) = +1.0.

The CG model is in good agreement with the atomistic simulations for the PMF as a function of R_g and R_e. Figure 8 shows that both models have a minimum in R_g ≈ 0.9 nm and R_e ≈ 2.5–3 nm. Both models also have a higher free energy metastable compact conformation. Overall, the agreement between CG and atomistic simulations is similar to what was seen for NaPSS. In fact, the overall behavior of the isolated polycation is remarkably similar to that of the isolated polyanion. The probability distribution of R_g for the CG and AA models, shown in Fig. S5, also shows the validity of the CG model.

The choice of topology in the CG model plays an important role. In Fig. 9, the properties of the present CG model are compared to those of CG₁, which has the same charge distribution as the CG model, but a topology similar to CG_mod2 [Fig. S1(b)]. Figure 9(a) compares the probability distribution function of R_g for CG₁ with the CG model developed here and the atomistic model. For CG₁, the maximum of this distribution is at ∼0.45 nm, which corresponds to a collapsed state for the polycation chain in dilute solution, in contrast to the results from atomistic simulations or the CG model of this work. Consistent with this, the Kratky plot of the form factor [Fig. 9(b)] shows a peak corresponding to a collapsed chain.

The results of the CG model are in good agreement with atomistic simulations for the conformational properties, counterion distribution, and microscopic dynamics. Figure 10 compares the CG and atomistic model simulations for the form factor defined as

where q is the momentum transfer variable and r_ij = r_i − r_j, where r_i is the position of bead i. To facilitate a comparison between models, only the backbone beads connected to the phenyl rings are considered in the calculation. There is good agreement between the models within uncertainties. In the scaling regime, both models predict P(q) ∼ q⁻¹, which corresponds to rod-like scaling, as expected in dilute solution.

The CG model is in qualitative agreement for the distribution of counterions around the polyion and the residence time of the counterion near the charged site. Figure 11(a) compares CG predictions for the pair correlation function, g(r), between the charged site and the counterion. In the atomistic model, the charged site is taken to be the nitrogen atom, and in the CG model, it is the TMA bead (labeled TB). The position of the peak is the same in both models, but the CG model has a higher first peak as well as stronger oscillations at larger separations. This is a consequence of the larger diameter of the sites in the CG model and is a known drawback of coarse-graining the van der Waals interaction of four heavy atoms into a single Lennard-Jones site. The peak intensities in the charge–charge interactions do not overlap because the radial distribution functions (RDFs) are dependent on the volume of the system. In addition, the species for which RDF is reported is the N atom and the Cl ion in the AA model and the TMA site and the Cl ion in the CG model.

There is no counterion localization, and the residence time of counterions near the charged site is similar in the two models. Counterion localization occurs when the counterions remain in the vicinity of polyion as compared to the bulk solution. However, here, we see that these counterions are exchanged from the bulk, and the average lifetime of these counterions around the charge sites of polyion is around 0.3 ns. The lifetime of counterion–polyion complex is calculated using the function

where n(t) is 1 if the counterion is within a distance of the first minimum in g(r) and zero otherwise. This correlation function is plotted in Fig. 11(b), which shows that the two models show a similar behavior for the counterion dynamics near the polyion. The residence time, τ, is defined as the time when ϕ(t) has decayed to 1/e of its initial value. The residence time takes on values of 0.31 ns for the atomistic model and 0.35 ns for the CG model. Although there is no localization of counterions, here we do observe some degree of counterion condensation, which is also referred to as the solvent separated ion pair. Further details are presented in Sec. III D.

The CG model is in reasonable agreement with the atomistic model for the self-diffusion coefficients. The self-diffusion coefficients for the center of mass of the polyion and counterion are computed from the (mean square displacement) MSD(t),

where r_i(t) is the position of the center of mass of the ith species at time t. Table I compares the self-diffusion coefficients in the two models. For the counterions, the agreement is excellent and for the polyions, the two models are within an order of magnitude. The CG model is not expected to be in quantitative agreement with the AA model for these dynamic properties because four water molecules are grouped into one site and features on smaller length and timescales are not reliable. It is, therefore, heartening that the two models predict self-diffusion coefficients that are of the same order of magnitude.

The conformational, structural, and dynamic properties of dilute solutions of PVBTMACl are very similar to those of NaPSS. In both cases, the polyions display rod-like behavior, as expected. Figure 12 depicts the variation of R_g as a function of N on a logarithmic plot. We find R_g ∼ N^ν with ν = 0.95 ± 0.005 and 0.97 ± 0.004 for PSS and PVBTMA, respectively, consistent within errors to a rod-like behavior⁹ and with experiments on NaPSS.^54,55 We are not aware of any experimental data for PVBTMACl. The variation of R_e divided by the chain length as a function of N is shown in Fig. S7. We can see that for both PSS and PVBTMA, the value decreases as we increase N. There is no stiffening of the chains as N increases.

The spatial correlations between the polyion charged site and counterions are similar in the two solutions. Figure 13(a) (top two panels) compares the g(r) between the charged site on the polyion and the counterion for various values of N. Interestingly, the oscillations in NaPSS g(r) are more pronounced, of lower period, and persist to larger distances than in PVBTMACl. This is a consequence of the higher charge on the TMA site compared to the sulfonate site. If the polycation charge distribution is altered so that the charge on the TMA is decreased, the g(r) is very similar to that in NaPSS (see Fig. S6 of the supplementary material).

The second peak in g(r) is slightly more prominent for NaPSS solutions, and this is reflected in the larger number of counterions in the vicinity of NaPSS compared to PVBTMACl. Figure 13(b) depicts the fraction of counterions within a cutoff distance of r_cut = 0.62 nm around the charged site of the polyion (n_ion is the number of counterions within the cutoff distance and n is the total number of counterions). This fraction is slightly higher for NaPSS than PVBTMACl, and in both cases, it decreases with increasing N over the range investigated. This ratio is slightly more pronounced in the case of polyanion, which is due to more counterions near the charged site as depicted earlier. However, these trends are contrary to that observed for freely jointed bead spring chains where the fraction of counterions near the polyion increases with increasing N.^18,56,57 We have made sure that the concentration of the polymer for each system is same irrespective of the system size.

We do not observe any counterion localization in these systems, and contacts of the polyion with the counterions decay within 1 ns. There is, therefore, no evidence of any long-lived dipole between the polyion and counterion, contrary to what is invoked to explain the so-called extra-ordinary transition in the dynamics.^58,59 When a counterion is within a distance corresponding to the first minimum in the g(r) (0.62 nm) between the charged site and counterion, we consider it a contact. The residence time is the duration of a contact. Figure 14(a) depicts the probability distribution of the residence time of counterions on the charged site of the polyion in PSS and PVBTMA. There is a peak at very short times and the distribution decays within 1 ns. There is a slightly higher residence time around PSS than PVBTMA, consistent with a higher second peak in g(r).

The short residence times of the counterions are also evident in the time variation in contacts made by the 16 counterions on a specific monomeric charged site (in a 16-mer PVBTMACl) over a 1 μ s trajectory [see Fig. 14(b)]. In the figure, the y axis is the index of the counterion. The x axis is the time of the simulation, and if the counterion is in contact with the monomer (in this case, the middle monomer on the chain), then a red dot is noted for the times during which it is in contact. As shown in the figure, there are no long lived contacts in the entire duration of the simulation. A similar behavior is seen for all the monomers on the chain.

We find that the correlations between the backbone and side group of the chain with water are similar for the polyanion and polycation (see Fig. S8). This suggests that there is no difference in hydration and the counterions do not penetrate beyond the charged polymer sites. The pair correlation function between the counterions is qualitatively similar for the polycation and polyanion solution (see Fig. S9), although the peaks are more prominent in NaPSS solutions.

The self-diffusion coefficients for polyion and counterions are similar in both solutions. Figure 15 depicts the self-diffusion coefficients for the polymers and counterions as a function of N. The scaling of D_self with N is consistent with the Rouse behavior, but the uncertainties in the exponent are quite large. For example, D_self ∼ N^−ν with ν = 0.8 ± 0.4 and 0.9 ± 0.2 for PSS and PVBTMA, respectively. The observed exponents with theoretical predictions for salt-free dilute solutions.⁶⁰ The self-diffusion coefficients of the counterions are less sensitive to N.

The electrophoretic mobilities of counterions and polyions are insensitive to N, and as a consequence, the transference number calculated from the mobilities is significantly different from that calculated from the self-diffusion coefficients. The electrophoretic mobility of species i, denoted by u_i, is calculated from the Onsager transport coefficients,

where z_i and c_i are the charge and concentration of species i, respectively, F is Faraday’s constant, and L_ij are the Onsager coefficients defined by

where $Δriα(t)≡riα(t)−riα(0)$⁠, r_iα(t) is the position of the center of mass of molecule α of species i at time t, n_i is the number of molecules of species i, k_B is Boltzmann’s constant, T is the temperature, and V is the volume.

Figure 16(a) shows the variation in the electrophoretic mobilities of the polyions and their counterions, respectively, as a function of N. The mobilities of the polyions are essentially independent of the chain length, consistent with theoretical predictions for rod-like chains.⁶⁰ 

Transference numbers can be calculated using the Nernst–Einstein relation from the self-diffusion coefficients, denoted by t_D, or from the mobilities, denoted by t_u,^61,62 which are defined as

and

respectively, where the subscripts c and p refer to the counterion and polyion, respectively, z_c is the charge on the counterion, and z_P is the charge on each monomer of the polyion. Figure 16(b) depicts the transference numbers as a function of N. In all cases, there is essentially no N dependence. For both solutions, t_D > 0.9 but t_u ≈ 0.4–0.5, suggesting that the Nernst–Einstein approximation is not valid because distinct correlations are not taken into account while computing the self-diffusion coefficients.

We have developed coarse-grained models for a polyanion (NaPSS) and a polycation (PVBTMACl) in the MARTINI framework with an explicit solvent and counterions and with the POL model for the water. We find that two factors, the topology of the chain and the distribution of partial charges in a monomer, play a crucial role in the robustness of the CG model. The CG model developed here is in very good agreement with atomistic simulations for the distribution of chain size, the static structure, pair correlations, and chain and counterion dynamics. Changes in the topology or charge distribution result in dramatically different properties include unphysical collapse of the charged chain. For dilute solutions, the polycations and polyanions display a similar behavior with a slight difference arising from a more diffuse counterion cloud in the case of the polycation. We see no evidence for counterion localization or long-lived dipoles with the longest counterion residence times of the order of a few ns and a peak in the residence time distribution function of the order of 10 ps. The success of the models, when compared to atomistic simulations, suggests that this is a good model system for the study of polyelectrolyte complexes and coacervates.

See the supplementary material for the mapping scheme for previously used PSS CG models, a comparison of P(q) and g(r) of AA and CG NaPSS, atomistic model partial charge distribution on the monomeric unit of PVBTMA, a comparison of the R_g probability distribution function of AA and CG NaPSS and PVBTMACl, variation in R_e as a function of N, and g(r) for hydration of polyions and counterion-self-interactions in dilute solutions. Force field files for PSS and PVBTMA MD simulations: https://github.com/KaurSupreet05/CG-PSS-PVBTMA.git. Movie for the collapse of the NaPSS model CG_mod2: https://github.com/KaurSupreet05/PSS-Na-Collapse.git.

This work was supported by the National Science Foundation under Grant No. CHE-1856595. All simulations presented here were performed using computational resources provided by UW–Madison Department of Chemistry HPC Cluster under NSF Grant No. CHE-0840494 and the UW–Madison Center for High Throughput Computing (CHTC). S.K. thanks Dr. Ajay Muralidharan for the technical help in performing free energy calculations and Dr. Tyler Lytle for valuable discussions. The authors extend their thanks to Dr. Martin Voegele for providing them raw data for comparison with the MARTINI NaPSS CG model and Dr. Sriteja Mantha for providing data of the previously proposed NaPSS CG model.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.