Molecular Dynamics simulations of concentrated aqueous electrolyte solutions

Transport properties of concentrated electrolytes have been analyzed using classical molecular dynamics simulations with the algorithms and parameters typical of simulations describing complex electrokinetic phenomena. The electrical conductivity and transport numbers of electrolytes containing monovalent (KCl), divalent (MgCl$_2$), a mixture of both (KCl + MgCl$_2$), and trivalent (LaCl$_3$) cations have been obtained from simulations of the electrolytes in electric fields of different magnitude. The results obtained for different simulation parameters have been discussed and compared with experimental measurements of our own and from the literature. The electroosmotic flow of water molecules induced by the ionic current in the different cases has been calculated and interpreted with the help of the hydration properties extracted from the simulations.


I. INTRODUCTION
Understanding the transport properties of electrolytes in aqueous solution is important in a wide range of electrokinetic phenomena such as streaming current experiments 1 , ionic transport in biological and syntetic nanochannels 2,3 or colloidal electrophoresis 4,5 .
Traditionally, electrokinetic and ionic transport phenomena have been described using primitive models in which the solvent is approximated as a continuum of dielectric constant ǫ. Although such models provide an accurate description of a wide variety of phenomena, they fail for cases in which the discrete nature of the solvent plays a fundamental role, or to describe phenomena related to specific hydration of ions.
An alternative theoretical approach to understand electrokinetic phenomena which circumvents the limitations of primitive models is the use of molecular dynamics simulations of ions in explicit solvent. Due to recent improvements in algorithms and computer power, all-atom molecular dynamics simulatios studies of electrokinetic phenomena in some realistic systems are now possible 6,7 . From such atomistic descriptions, one can elucidate the microscopic mechanisms responsible for macroscopic measurable phenomena. For example, Aksimentiev and Schulten studied at atomic detail the permeation of individual ions through the transmembrane channel α-Hemolysin with the help of molecular dynamics simulations 6 . This approach has also been employed to understand from a microscopic perspective the electrophoresis of DNA immersed in multivalent ionic solutions when an external electric field is applied 8 , the transport of monovalent and divalent ions through polymeric and silica nanopores 9,10 , and the ionic selectivity of the OmpF porin biological nanochannel 11 .
The realistic atomic description of such electrokinetic phenomena involves complex systems containing large number of particles. To be able to cope with such systems using molecular dynamics simulations, the use of algorithms which enhance the computer performance are required. A molecular dynamics simulation package which has been proven to be very successful in the description of biological molecules and elecrokinetic phenomena is NAMD2 12 . To accelerate the calculations of such large systems, in NAMD2 the temperature T of the system is controlled by using the Langevin thermostat instead of the more rigorous but much more computationally demanding Nosé-Hoover thermostat. Also to speed up the simulations, NAMD2 is usually employed with a multiple time step: a basic timestep for short range interactions and a longer one for the evaluation of k-space contribution to the long range electrostatic forces in the PME technique.
In this paper we study, under different conditions, the transport properties obtained from Molecular Dynamics simulations of electrolyte solutions using the same algorithms and conditions employed in simulations of complex systems. In particular, we focus our analysis on the transport properties of a monovalent electrolyte (KCl), a divalent electrolyte (MgCl 2 ), the mixture of both (KCl + MgCl 2 ) and a trivalent electrolyte (LaCl 3 ) when an external electric field is applied. While the use of the KCl electrolyte is standard in molecular dynamics studies and its transport properties have been briefly analyzed from this perspective before 6 , the reliability of such simulations methods in dealing with transport properties of multivalent electrolytes have not been tested in detail in spite of the rich phenomenology that they originate 1,2,13,14 and their use in the simulation of electrokinetic phenomena involving divalent and trivalent ions 8,15-17 .

A. Simulation Methods
We have studied the transport properties of different electrolytes by performing molecular dynamics simulations under external electric fields of different magnitudes. The systems considered were ionic solutions of KCl, MgCl 2 , a mixture of KCl and MgCl 2 , and LaCl 3 , see Table I. To study the effect of the system size on the ionic transport properties of bulk electrolytes extracted from molecular dynamics simulations, two cubic simulation boxes of different size were used, L ≃ 4nm and L ≃ 8nm.
All simulations have been carried out with the molecular dynamics simulation package NAMD2 12 , since it is widely employed in the simulation of biological macromolecules. Water was described using the TIP3P water model as implemented by the CHARMM force field. The ions were modeled as charged Lennard-Jones particles with parameters given by the CHARMM force field for K + , Mg 2+ and for Cl − , while the Lennard-Jones parameters for La 3+ were taken from Ref. 18 .
The initial configuration of the simulation was constructed as follows. The ions (K + , Mg 2+ and Cl − ) were inserted at random positions by employing the AutoIonize plugin of the Visual Molecular Dynamics (VMD) software package 19 inside a cube of preequilibrated TIP3P water molecules, obtained with the help of the Solvate plugin of VMD (see Fig. 1 for a diagram of the system). The faces of the cube are parallel to the XY, XZ, and YZ planes. In all simulations we employed periodic boundary conditions in all directions. Lennard-Jones interactions were computed using a smooth (10-12Å) cutoff, as it is customary done in NAMD2 simulations 6,8 . The electrostatic interactions were calculated using the particlemesh Ewald (PME) method with a precision of 10 −6 using a 128 × 128 × 128 grid and a 12Åcutoff for the real space calculation. These are common parameters used to simulate complex and large biological macromolecules 6 . For each case, the equilibration procedure consisted of 50000 steps of energy minimization, a 1ns run in the NpT ensemble (with p = 1atm and T = 296K) followed by a 1ns run in the NVT ensemble (T = 296K). Production runs in the NVT ensemble were performed in the presence of different electric fields to induce electromigration of the ions in solution. The NpT simulation runs were performed employing a combination of the Nosé-Hoover constant pressure method with the piston fluctuation control implemented using Langevin dynamics as implemented in NAMD2 (p = 1atm, period of the oscillations of 0.1ps and relaxation constant of 0.05 ps). As mentioned above, to speed up calculations the NVT runs were carried out by applying Langevin dynamics, with parameters (also in the NpT run) T = 296K and a damping coefficient of 1ps −1 (using 0.2ps −1 to test its effect in the dynamics for some cases specified later). Langevin forces were applied to all atoms except for hydrogens, which thermalize through interactions with the rest of the system.
In all cases, the equations of motion were solved using a multiple time step in order to speed up the simulations. A basic time step of 2fs was used for the evaluation of short range interactions and a longer time step of 4fs for the calculations of the k-space contribution to the long range electrostatic forces in the PME technique.
In the production runs, the instantaneous current induced by the external electric field applied along the Z-direction is calculated with the help of 6 where z i and q i are the Z-coordinate and the charge of atom i, respectively. L is the size of the simulation box and ∆t is the time interval employed to record data, which was chosen to be ∆t = 10ps. The average current is computed by linearly fitting the cumulative current that is obtained by integration of the instantaneous current. To ensure constistency of the results, the current was also computed by counting the flux of ions crossing a plane perpendicular to the direction of the electric field and located in the center of the simulation box. The conductivity κ of the solution is defined by: where A = L 2 is the cross section area perpendicular to the electric field E. In order to obtain the conductivity of the electrolyte, we have performed simulations at different electric fields E and calculated the current I induced by them. In all cases, an ohmic behaviour has been observed (i.e. consistent with a linear relation between I and E), and a linear fit of the data to Eq. (2) gives the value of the conductivity κ.
The electromigration of ions in aqueous solutions as a result of the application of an external electric field is often accompanied by a net flow of water, the so-called electroosmotic flow. The electroosmotic flow has been evaluated by keeping track, every ∆t = 10ps, of the accumulated number of water molecules crossing a plane perpendicular to the electric field. Crossings of such plane in the direction of the electric field are counted as positive, whereas crossings in the opposite direction are counted as negative. Hence, the sign of the overall flow gives the direction of the electroosmotic flow with respect to the elecric field direction. As reported before 6 , the simulations might give a drift of the whole system which is unphysical since no net force is applied to the system (the electrolyte is globally electrically neutral). To avoid such spurious effects, the computation of the electroosmotic flow has  The electric current flowing in an electrolyte solution is caused by the motion of anions and cations moving in opposite directions under the applied field. The fraction of the total current induced by each ion type defines its transport number, which is in general a function of the electrolyte concentration. The fraction of electrical current carried by cations defines the cationic transport number, t + , and the fraction of the electrical current carried by anions the anionic transport number, t − . A completely equivalent quantity to transport numbers is the ratio of the current induced by anions (anionic current) over the current induced by cations (cationic current), which will be be refered as the ratio of currents throughout the paper. Due to the global electroneutrality of the system, the total electric current flowing through an electrolyte as a result of the application of an external electric field is independent of the frame of reference in which it is measured. The ratio of currents and transport numbers, however, are frame dependent, and their computation must be done carefully. A proper account of the contribution to the total current from every ion in electrolyte solutions is important to describe phenomena in other more complex systems in which not only the total electric current but also the flow of each type of ion is relevant (e.g. in the study of the selectivity of ionic channels). Experimentally, the relevant frame of reference in which data of transport numbers and ratio of currents is provided is the frame of reference of the moving fluid. Therefore, to facilitate the comparison with experimental values, transport numbers and ratio of currents will be given in the frame of reference comoving with the fluid.

B. Conductivity measurements
The electrolyte conductivity measurements were performed using a MeterLab CDM210 (R) conductivity meter, Radiometer Analytical SAS (France). The solutions were prepared using water from a Water Purification System Millipore Simplicity 185. Magnesium Chloride 6-hydrated (MgCl2) and Potassium Chloride (KCl) from Panreac in all cases following ACS specifications. Weighting of the compounds were done with a Mettler Toledo AB104-S, in the quantities necessary to achieve a 1M concentration. The conductivity measurement included a stirring of the solution with a magnetic stirrer and heater JPSelecta Agimatic-E, with temperature monitoring using a laboratory thermometer at 296.0 K.

A. A test case: Transport properties of 1M KCl
The results obtained from molecular dynamics simulations for the electrical properties of the 1M KCl electrolyte are given in Fig. 2 and Table.II. We have used two different sizes for the simulation box, L = 3.88nm and L = 7.82nm, to test the dependence of the results on the simulation's box size. As can be seen, the values for the ionic conductivity for the two different box sizes do not differ significantly, κ = 13.4 S/m for the small box and κ = 12.6 S/m for the big one. Both results are in good agreement with our measured experimental value κ exp = 11.24 ± 0.01 S/m, being the value of the larger system closer to the experimental result as it should be expected. Note that in previous studies of 1M KCl bulk electrolyte (see Ref. 6 ), it was argued that a 0.2 ps −1 damping constant is necessary in order to reproduce correctly transport properties. However, our present results were obtained using a Langevin thermostat with a damping constant of 1 ps −1 , which is the typical choice for simulations of complex systems in contact with electrolyte (such as protein channels or silica nanochannels 6,10,11 ). Our results show that with the standard simulation parameters employed in complex systems the conductivity of the electrolyte KCl is correctly reproduced.
We have also calculated the ratio of the different contributions to the total current from anions and cations (equivalent to the ratio of transport numbers). As explained in the previous section, it is defined as the ratio of the anionic and cationic currents in the frame of reference of the moving fluid. The results for such ratio of currents are given in Table II. As shown in the table, the values calculated from molecular dynamics simulations are in good agreement with the experimental value of I Cl /I K = 1.048 (equivalent to a value of the cation transport number t + = 0.488 20 ) for both sizes of the cubic simulation box.

Electoosmotic flow
In addition to ionic currents, we also observe a net flow of water molecules. In Fig. 3, the accumulated number of water molecules per unit area crossing a plane perpendicular to the electric field is represented as a function of simulation time. A linear fit of these magnitudes provides the flux of water molecules for every value of the electric field, as shown in Table III. As shown in Fig. 3 and Table III, the electroosmotic flow for the 1M KCl electrolyte is in the direction opposite to the external field, that is, in the direction of the flow of chloride ions.

B. Transport properties of 1M MgCl2
A similar procedure has been followed to obtain the transport properties of a 1M MgCl 2 electrolyte. The electric current induced by the electromigration of ions was obtained for different values of the electric field using two different sizes of  The ratio between the different contributions of anions and cations to the total current has been evaluated. The results for such ratio of currents are given in Table IV. These results emphasize the need for using a large enough simulation box. The values for the ratio of currents and transport numbers obtained for the cubic simulation box of size L = 3.82nm exhibit a spurious dependence on the applied electric field. For the larger simulation box, the results are not field dependent, as should  be expected. Such results show that the anionic contribution to the current is significantly larger than the cationic contribution (being a factor of 1.4 between them). In our simulations these differences in the anionic and cationic contributions to the current can be attributed, to a large extent, to differences in diffusion coefficients between both ions. In order to disentangle the diffusional and correlation contributions to the transport number, we have evaluated the translational diffusion coefficient of each ion by computing the mean square displacement of each ion in a 2ns long NVE simulation run with no external field applied. The results of D Mg = 0.95 × 10 −5 cm 2 /s and D Cl = 1.69 × 10 −5 cm 2 /s for the diffusion coefficients of Mg 2+ and Cl − , respectively, lead to a ratio between the diffusion coefficients of D Cl /D Mg = 1.8. Experimentally, both diffusion coefficients and transport numbers can be obtained, so we can compare both simulation results with experimental data. Using NMR, Struis et al. 23 ob- The simulations reported here reproduce with good accuracy the diffusion coefficient for Mg 2+ but the diffusion coefficient for Cl − is greatly underestimated. Concerning transport numbers, electrochemical measurements 21,22 give a transport number for cations around 0.3, so the experimental ratio between anionic and cationic currents is ≃ 2.3. The difference between the ratio of transport numbers obtained by electrochemical methods and the ratio between diffusion coefficients obtained by NMR can be interpreted in several ways. First of all, thermodynamic arguments 24 show that some experimental procedures mixed up different reference frames, so the electrochemical results have to be interpreted with caution. If the difference between both ratios is indeed physical, the difference can be attributed to electrokinetic processes (not accounted by diffusion coefficients) which are supposed to affect the transport numbers of each ion 20 . In any case, these electrokinetic processes were predicted in the framework of classical, continuum theory of  electrolytes and are not clearly observed in our simulations, so its molecular origin remains unclear. The effect of the damping constant of the Langevin thermostat was investigated by performing molecular dynamics simulations on the same system of size L = 7.73nm and same conditions but using a Langevin thermostat with a damping constant of τ Lan = 0.2 ps −1 (instead of the damping constant of 1 ps −1 employed in all other simulations). The electrolyte exhibits an ohmic behaviour, with conductivity κ = 13.5 S/m, which slightly differs from the value obtained before and from the experimental value. Nevertheless, the major effect of the damping constant is in the different ionic contributions to the total current, as can be seen in Fig. V. In comparison to the results of the simulations with τ Lan = ps −1 , it is evident that there is a spurious dependence of the ratio of the anionic current over the cationic current on the magnitude of the applied electric field.

Electoosmotic flow
The electroosmotic flow induced by the ionic current was also computed. In Fig. 6, the accumulated number of water molecules per unit area crossing a plane perpendicular to the electric field is represented as a function of simulation time.
The flux of water molecules for each electric field is given in Table III. It is interesting to note the different direction and magnitude of the water flow obtained for 1M KCl and for 1M MgCl 2 electrolytes. Indeed, while in the presence of 1M KCl the direction of the net flow of water is opposite to the direction of the electric field, in the presence of 1M MgCl 2 the net flow of water goes in the direction of the applied field. The magnitude of the net flow of water also differs significantly in both cases, being much larger (almost an order of magnitude) for the 1M MgCl 2 electrolyte. These differences can be understood from the hydration properties of the ions involved. For each ion, we have computed its hydration, i.e., the average number of water molecules in its first coordination shell (as defined by the first minimum of the radial distribution function, see Figs. 7 and 8). We have also computed the average number of water molecules which have remained a time τ in the first coordination shell of each ion to test the robustness of the hydration layer, see Table VI. The values given in Table VI do not show any dependence on the value of the electric field applied in the simulations. The reported results for the hydration values of the different ions agree with results from previous molecular dynamics studies [25][26][27] and Ab initio calculations 28 . The results for the average number of water molecules that spend a time τ in the first coordination shell of the different ions are also in agreement with previous computational studies, in which residence times of water molecules in the first coordination shell of the order of ∼ 10ps were obtained for K + and Cl − and of the order of ∼ 10ns for Mg 2+ . Such difference in the robustness of the hydration layer of Mg 2+ and the other ions is a broadly accepted fact, established by ab initio and DFT calculations 28,29 as well as by NMR, Raman spectroscopy, and X-ray adsorption spectroscopy experiments 23,30,31 .
Hence, from the analysis of the hydration properties of the ions it is possible to understand the difference in the electroosmotic flow between the 1M KCl and 1M MgCl 2 electrolytes. For the 1M KCl electrolyte, the anionic current is slightly greater than the cationic current and the hydration of Cl − is higher than that of K + so it exhibits a net electroosmotic flow in the direction of the flow of Cl − . For the 1M MgCl 2 electrolyte, although the current of Cl − is greater than the current of Mg 2+ , the latter ion is much more effective dragging water molecules (its hydration layer is much more robust over time) so the net electroosmotic flow is in the direction of the flow of Mg 2+ along the applied electric field.

Concentration dependence of transport properties of MgCl2
For the MgCl 2 electrolyte we have also tested the dependency of the conduction properties of the electrolyte on its concentration, c. We have performed simulations at different concentrations of MgCl 2 for a fixed value of the applied electric field E = 0.0433V/nm in a cubic simulation box of side L ≃ 8nm (see Table I for precise values). The summary of the results are given in Table VII and Figures 9 and 10. Here, for each value of the electrolyte concentration the conductivity has been obtained from a single point (I, E), under the assumption that the electrolyte exhibits an ohmic behavior (suggested by the studies at 1M concentration, see Fig. 4). As expected, the electrical conductivity, κ, of the MgCl 2 solution increases with increasing concentration, as can be seen in Table VII and Fig. 9. In Fig. 9, the conductivity of the electrolyte obtained from molecular dynamics simulations for each concentration is compared with experimental results 21 and the Kohlrausch's limiting law. The results obtained from MD simulations reproduce well the dependency of the electrolyte conductivity on the concentration, as compared to the experimental tendency. Besides, the absolute values of the experimental and simulated conductivities agree well with each other, especially for high and low salt conditions. In Fig. 9 we represent the ratio of the anionic current over the cationic current as a function of the MgCl 2 concentration. The results from MD simulations indicate that the contribution of the anions becomes more prominent with increasing concentration of electrolyte. This tendency agrees with the experimental behaviour for 2 : 1 electrolytes 20 .

C. Transport properties of mixture composed of 1M KCl and 1M MgCl2
We have also studied the transport properties of an electrolyte composed of 1M KCl and 1M MgCl 2 . The addition of monovalent ions to a system immersed in a multivalent electrolyte is sometimes used experimentally to screen the electric charge of the multivalent ions 1 . This effect is also used in experiment and simulations to determine whether electrostatic correlations are responsible for a certain macroscopic effect 11,32,33 . The analysis is similar to the analysis done for the previous electrolytes. The electric currents caused by the drift of ions have been obtained for four different values of the applied electric field. In this case, a single cubic box of size L = 7.72 nm was employed. The results for the conductivity are summarized in Fig. 11, which results in a value for the conductivity κ = 14.8 S/m. Our own experimental measurements give κ = 16.83 ± 0.01 S/m. The different contribution to the total current of the different ions, as well as the cationic transport number are given in Table VIII for every value of the applied electric field.

Electroosmotic flow
The results for the net flow of water induced by the ionic current are given in Fig. 12, in which the accumulated number of water molecules per unit area crossing a plane perpendicular to the electric field is represented as a function of simulation time. The flux of water molecules for each electric field is given in Table III. Similarly to the case of the 1M MgCl 2 electrolyte, the net electroosmotic flow is in the direction of the flow of cations, in spite of being smaller than the net flow of anions (see Table VIII). We interpret this result along the same lines as with the case of 1M MgCl 2 electrolyte. As it is shown in Table VI, the hydration layer of Mg 2+ is much more robust than the hydration of Cl − and K + , so the electroosmotic flow induced by the flow of Mg 2+ dominates.

D. Transport properties of 1M LaCl3
The lanthanum cation La 3+ is a highly charged and polarizable ion which, for most problems, requires the use of Ab  33 . For such cases, La 3+ can be modeled as a charged Lennard-Jones particle and one can obtain the transport properties of a LaCl 3 ionic solution by using classical molecular dynamics calculations. The Lennard-Jones parameters used here to describe the lanthanum cation La 3+ were taken from Ref. 35 The electrical transport properties of the 1M LaCl 3 electrolyte in the presence of an external electric field are obtained using molecular dynamics simulations which employ a nonpolarizable force field (described above). The summary of the results is given in Fig.13 and Table IX. The electrolyte shows an ohmic behaviour (see Fig.13), with an electrical conductivity of κ = 12.8 S/m. Considering the limitations of the description, this value is quite satisfactory compared to the experimental value of κ exp = 15.3 S/m 36 . The ionic contribution, the ratio of currents and the cationic transport number are given in Table IX for each value of the external electric field.

Electroosmotic flow
The results for the net flow of water induced by the ionic current are given in Fig. 14, in which the accumulated number of water molecules per unit area crossing a plane perpendicular to the electric field is represented as a function of simulation time. The flux of water molecules for each electric field is given in Table III.
The hydration properties of the 1M LaCl 3 electrolyte are summarized in Table VI. In spite of the limitations of the non-polarizable model used to describe La 3+ , the values obtained for the hydration of La 3+ are close to the values given in the literature (hydration ∼ 8 − 9 and a residence time of hydrated water in the first coordination shell of ∼ 1ns 34 ). As showed in Fig. 14, the electroosmotic flow obtained for the 1M LaCl 3 electrolyte is in the direction of the flow of chloride ions, opposite to the direction of the applied electric field. In this case, the higher hydration of La 3+ is not enough to compensate the higher flow of chloride ions versus the flow of lanthanum cation La 3+ .

IV. CONCLUSIONS
Classical molecular dynamics simulations of various electrolytes (KCl, MgCl 2 , KCl + MgCl 2 , LaCl 3 ) in external electric fields were performed to study their transport properties. We have employed algorithms and parameters typically used in simulations of complex electrokinetic phenomena to test against experimental data their validity in the description of transport properties of electrolytes. The electrical conductivity and transport numbers (ionic contributions to the total current) of electrolytes containing monovalent (1M KCl), divalent (1M MgCl 2 , 1M MgCl 2 + 1M KCl ), and trivalent (1M LaCl 3 ) cations were computed from the simulated ion trajectories. It is shown that in all cases the electrical conductivities obtained from the simulations are in good agreement with our own experimental measurements and values from the literature. Also, for the MgCl 2 electrolyte, the dependence of the electrical conductivity on the electrolyte concentration was investigated and found that the MD results follow the experimental trend and give accurate results for the lowest (0.1M) and the highest concentrations (1M) studied. Transport numbers obtained for the 1M KCl electrolyte agree well with the experimental values found in the literature. For the divalent electrolyte 1M MgCl 2 , there is a significant discrepancy between the transport numbers obtained from MD simulations and the experimental values, probably due to the poor result for the self difusion coefficient for anions provided by the Lennard-Jones parameters used in the simulations. The effect of the simulation box size on the electrical transport properties was also explored. It is found that a big enough simulation box is needed to avoid spurious effects of the simulation. The electroosmotic flow of water induced by the ionic flow has also been computed from the MD trajectories. It is shown that the direction and magnitude of the water flux depends on the electrolyte: while the flux is in the direction of the cation flow for electrolytes containing Mg 2+ , it is in the opposite direction and weaker for all the other cases. These results are interpreted with the help of the hydration properties of the ions, which are calculated and successfully compared with previous studies on the subject.