ions such as calcium , magnesium , sodium , potassium , and 
 chloride take part in vital bodily functions such as the generation 
 of action potentials , catalytic actions of enzymes , signal transduction , 
 and structural stabilization of nucleic acids . furthermore , 
 a large number of drugs contain metal ions , which allow those drugs 
 to interact with endogenous proteins . metals 
 interact with the biological molecule , either directly or through 
 the water present in their first solvation shell . divalent 
 ions , especially calcium ( ca ) and magnesium ( mg ) , play an important role in physiology as both act as cofactors 
 in many enzymatic reactions . mg drives the folding of 
 rna into stable states and 
 also participates in energy metabolism by binding to adenosine triphosphate 
 ( atp ) . most 
 researchers study physiological activity of ion - related mechanisms 
 by estimating the thermodynamic quantities of ions . they find their estimates by performing experiments on 
 mixtures of cations and anions in aqueous solution . to tease out 
 the properties of an individual ion , researchers have to make several 
 assumptions , leaving the values of an individual ion to reside within 
 an approximate range . however , they observed that accurate calculations 
 of thermodynamic quantities require careful selection of parameters 
 so that their theoretical values would match with the experimentally 
 derived values . in simulations , ions are represented using three basic 
 parameters : a charge and two lennard - jones ( lj ) potential parameters . the charge represents the valency of the ion , and the lj parameters 
 determine the size . the experimental measure most widely used for 
 validation is the hydration free energy or the negative work required 
 to insert a single ion in bulk water . despite careful selection of parameters and perfect match with these 
 experimental quantities , researchers have faced challenges in reproducing 
 accurate physiological behavior of biomolecules interacting with divalent 
 ions . historically , simulations of divalent ions with biomolecules 
 have suffered from two main problems . first , an accurate structure 
 of the coordination of ions with the biomolecules can not be obtained , 
 and second , the introduction of higher concentration of ions in a 
 solution leads to excessive aggregations and formation of ion clusters . the classical model of ions contains the entire charge of the ion 
 at the center . this property prevents the ion from responding directly 
 to the molecular environment and thus limits their usage with bigger 
 and complex biomolecules in classical simulations . researchers observed 
 that when they fit the parameters after precise matching with the 
 experimental hydration free energy , they were able to describe the 
 structure of the first hydration shell of ions with considerable accuracy . 
 for mg and ca , the ion - oxygen distances 
 in water have been reported as 2.12.5  and 2.43.0 
  , which lie close to the experimental values of 2.1 and 2.4  , respectively . the behavior of the metal cation when it interacts 
 with a biomolecule , however , involves substantial charge - transfer 
 effects , which are hard to capture with the localized , fixed , charged 
 model mentioned above . an external constraint placed between a ca ion and oxygen atoms in a theoretical study of ca binding with parvalbumin provided evidence of how the charge in 
 the center of the ion by itself is not sufficient to keep the ca ion stable in its binding site . another theoretical study of ca binding in calbindin 
 reported that a specialized treatment of electrostatics was necessary 
 to keep the ca ions stably bound to the two calcium binding 
 sites of the protein . the accurate structure 
 of direct coordination of a cation with biomolecules becomes much 
 more important in studies involving the description of conformational 
 changes occurring after ion binding . mg binding to unstructured 
 rna molecules is known to shift the folding versus unfolding equilibrium 
 toward folded states . similarly , the activation 
 mechanism of the calmodulin family of proteins also involve large 
 conformational changes occurring after ca binding . theoretical studies of this phenomenon have been extremely challenging 
 with the ion models in which the charge is located at the center of 
 the ion . various approximations for the interaction 
 of ions with biomolecules 
 also result in high - energy barriers for change of coordinating partners , 
 a process that occurs frequently in biological systems . for example , 
 during passage of ions through ion channels , the ions either partially or completely dehydrate , changing their 
 coordination partners from water oxygens to the protein oxygens . also , 
 in the presence of ca ions , the ef - hand family of proteins 
 replaces the already bound mg ion with the more preferred 
 ca ion . researchers have proposed several solutions to 
 treat this problem , including reparameterizing the charges of the 
 coordinating atoms or using polarizable 
 parameters for all atoms . these solutions are 
 either too specific for the system , rendering the parameters nontransferable , 
 or require the use of extensive computational resources . in a previous 
 work , saxena and sept offered a workable 
 solution with the development of a new model for the ions . a multisite 
 model was developed for divalent ions mg and ca , where the localized charge of the ion was split into multiple sites . this model was able to accurately represent the structure of the first 
 coordination shell of the ion in water , proteins , and nucleic acids . saxena and sept demonstrated that the model could be successfully 
 used to find conformational changes and thermodynamical changes associated 
 with ion binding in biomolecules , such as ef - hand proteins and rna . the second major problem faced by researchers in theoretical studies 
 of ion - mediated mechanisms is ion aggregation . ions at low concentrations 
 behave as ideal solutions , but as the ion concentration increases , 
 the solutions move away from ideality . the study of biological functions 
 of proteins requires the experiments to be performed at physiological 
 ion concentrations . at such concentrations , , at these concentrations , 
 the cation and anion interact with each other through a water molecule 
 forming a solvent separated pair ( ssp ) . simulation studies of ion solutions at finite concentrations 
 have found untimely aggregation of ion pairs . in - depth analysis 
 of the ion clusters showed the that at conditions where the oppositely 
 charged ions should have formed an ssp , they instead formed direct 
 contact pairs . this behavior is due to the inaccurate representation of ion the aforementioned artifact thus 
 leads to excessive ion pairing , poor solubility , and spontaneous crystallization 
 of ion solutions in theoretical calculations . the multisite 
 model developed by saxena and sept also suffers with 
 the excessive aggregation problem in concentrated mgcl2 and cacl2 solutions . in this work , we refine the parameters 
 of the multisite model such that it is devoid of both artifacts . the 
 multisite model is already able to predict accurate inner shell coordination 
 with biomolecules . to eliminate the ion aggregation problem , we modify 
 the interactions of the divalent cations with the monovalent anion 
 to reproduce the physical behavior of concentrated ion solutions , 
 while maintaining the parameters for the interaction of the ions with 
 water and biomolecules . experimentally , the physical behavior of regular 
 ion solutions can be studied by measuring properties such as osmotic 
 pressure and freezing point , but the 
 calculation of these properties through simulations has been difficult . osmotic pressure 
 calculations offer a method to measure the strength of solute interactions . 
 at a similar concentration of solutes , aggregating molecules apply 
 less pressure on the walls and results in reduced osmotic pressure , 
 while completely dissociated solute molecules move around in water 
 and exert more pressure on the walls , resulting in increased osmotic 
 pressure . by matching the calculated osmotic pressure to the experimental 
 values , the solute interaction can be improved to accurately reproduce 
 the behavior of concentrated electrolytes . this method has been successfully 
 used to obtain accurate behavior of aqueous solutions of nacl and 
 kcl and to obtain accurate behavior of 
 trimethylamine n - oxide . we use this method here to calculate the osmotic pressure of mgcl2 solutions at various concentrations and match them to the 
 experimental values to refine the parameters of the multisite model . anion radial distribution functions obtained after 
 the parameter refinement show that the solutions do not form excess 
 direct ion - pairs , matching the experimentally observed behavior . subsequently , 
 we use the refined parameters of the multisite model to study the 
 behavior of cacl2 solutions and find that the calculations 
 reproduce the experimental osmotic pressure of cacl2 with 
 high accuracy , without further refinement . structural characterizations 
 confirm that the divalent ions show accurate inner shell coordination 
 and are devoid of unusual ion  ion pairing at high salt concentrations 
 in both solution types ( mgcl2 or cacl2 ) . the 
 model captures the exchange of inner shell waters and better coordination 
 geometries with its binding partners . further , in comparison with 
 the work by saxena and sept , we show 
 that the model retains the ion solution properties seen in that study , 
 indicating that the accuracy of the model is not compromised after 
 the changes made in this work . through this study , we present an improved 
 multisite model , which can not only be used for prediction of free 
 energies and selectivity of divalent cation ( mg and ca ) binding sites in proteins and nucleic acids , but can now 
 be used over a wide range of physiological salt concentrations . additionally , 
 the model can be easily used with standard fixed charged force fields 
 with no extra computational cost . the multisite model consists of a central 
 atom m and n dummy charge centers , d , where n denotes the coordination number 
 of the ion . the dummy atoms are located along the vectors connecting 
 the ion center and the coordinating atoms , determined from the coordination 
 structure of the ion , octahedral for mg and pentagonal 
 bipyramidal for ca . the equilibrium bond length between 
 the ion center and the charged sites was fixed at 0.9  . each 
 dummy atom in a given ion has the same atom type and carries identical 
 charge and lennard - jones ( lj ) parameters . while the lj parameters of 
 the dummy atoms remain the same among different ions , the lj parameters 
 of the central atom differ and are based on their respective match 
 to the hydration free energy . the behavior of ions in concentrated 
 electrolytes causes the solutions to deviate from the vant 
 hoff ideal solution . osmotic pressure 
 has been widely used to estimate thermodynamics quantities , as it 
 is relatively simple to measure . osmosis is 
 the process where the solvent spontaneously passes through a semipermeable 
 membrane separating the solute the pressure applied by the 
 solute on the membrane is related to the concentration of the solute 
 by morse equation : where i is the dimensionless 
 vant hoff factor , m is the molarity , r is the gas constant , and t is the thermodynamic 
 ( absolute ) temperature . recently , luo and roux developed 
 a method for calculating osmotic pressure in simulations . the methodology involves creating virtual walls 
 in a simulation box , which act as semipermeable membranes allowing 
 only water to pass through . the average force required to keep the 
 solute molecules inside the simulation box is used to calculate the 
 osmotic pressure . the osmotic pressure also offers a measure of the 
 strength of interactions between solutes . at a similar concentration 
 of solutes , aggregating molecules apply less pressure on the walls 
 and result in reduced osmotic pressure , while completely dissociated 
 solute molecules move around in water and exert more pressure on the 
 walls , resulting in increased osmotic pressure . due to these characteristics , the osmotic pressure calculations provide a means to improve the parameters 
 of the solutes in simulations by matching with the experimental values . 
 here , we follow the same methodology as luo and roux to determine 
 accurate interaction between mg and cl ions and ca and cl ions . we systematically 
 modify the van der waals parameters of the divalent ions such that 
 the experimental osmotic pressure values of these solutions are matched . in particular , we modify the interaction of the cation dummy sites 
 with cl ions . all other parameters are kept the 
 same as in saxena and sept , such that the interactions of the divalent 
 cations with water and biomolecules do not change . we calculate the osmotic pressure 
 of mgcl2 and cacl2 solutions at various molalities 
 using molecular dynamics simulations at constant temperature and pressure . the simulation system consists of a central box of the ionic solution 
 flanked by a pure water box of equal volume . the two boxes were 
 generated individually at the required size , equilibrated using berendsen 
 pressure coupling and then concatenated along the z - direction . periodic boundary conditions ensured that another pure 
 water box flanked the center box from the other side . the resulting 
 box was energy minimized and equilibrated for 1 ns to allow equilibration 
 of water and the proper mixing of the ions . the final state obtained 
 from the previous step is then used to start three independent simulations , 
 each 10 ns long . simulations were conducted at constant pressure and 
 constant cross sectional area for each desired molality . this process 
 was followed individually for each ion concentration of the central 
 box . the wall restraints were applied along the virtual 
 walls by modifying gromacs 4.0.7 to enforce a flat - bottom , half - harmonic 
 potential . the osmotic pressure was calculated from the force exerted 
 by the wall to keep the solute inside the simulation box ,  
 = ( fwall/a ) , where a is the cross - sectional area of the virtual 
 wall obtained from the box length of the central box for each ion 
 concentration ( table 2 ) . fwall = ( k / np)npi(|zi|  |zwall| ) , with |zi| > |zwall| , k as the force constant ( 10 kcal / mol/ ) , np as the total number of configurations 
 in the production stage of each simulation ( 5 ns ) , and i as the index of ions . the force is averaged over two half - harmonic 
 walls to get the osmotic pressure . errors in the averages are estimated 
 using 2.5 ns block averages over the production stage of the three 
 independent simulations . the multisite model was used for mg and ca parameters , tip3p for water molecules , and cl parameters were obtained 
 from luo and roux s work . the 
 parameters for the cl anions with the cation dummy 
 sites were adjusted iteratively to match experimental osmotic pressure . ( a ) osmotic pressure 
 obtained for different concentrations of mgcl2 in water . the black line is the osmotic pressure data obtained 
 from experiments , and the magenta line 
 is the osmotic pressure obtained from simulations using the multisite 
 mg model . ( b ) snapshot of a simulation configuration 
 of mg and cl ions in water showing 
 aggregation ( mg as red spheres , cl as green spheres , water molecules occupying the first solvation 
 shell of mg shown in licorice ) . ( c ) radial distribution 
 function between mg and cl ions at 
 different concentrations of mgcl2 ( m represents 
 molality = mol / kg ) . ( d ) total number of cl ions 
 binding to mg ions at different concentrations of mgcl2 . we calculate 
 the osmotic pressure of mg in solution with cl ions at different concentrations . when comparing with the experimental 
 values , we find that at lower concentrations ( < 0.5 m ) the calculated 
 values of osmotic pressure match well with the experiments , but at 
 higher concentration the calculated values drop significantly and 
 the deviation grows larger as the concentration increased ( figure 1a ) . the decrease in osmotic pressure indicates that 
 the solute ( ions ) might have a strong attraction between them . to 
 verify this the mgcl2 mixtures , at almost all concentrations , show spontaneous 
 aggregation of ions , which occurs early in the simulations and remained 
 unchanged for the entire length of the simulation ( 100 ns , figure 1b ) . we found several instances where cl ions occupied the first coordination shell of mg ions 
 and vice versa . the radial distribution function ( rdf ) between mg and cl ions is proportional to the probability 
 of finding cl ions within a certain radius around 
 mg ions . the rdfs show that there is a large population 
 of cl ions located at a distance of 1.95  
 from the mg ions ( figure 1c ) , 
 indicating that the two ions form direct contact pairs . this behavior 
 does not match x - ray experiments of mgcl2 solutions , where 
 the minimum distance between mg and cl ions , below 3 m concentration , is close to 4.5  . structurally , this value corresponds to a solvent - separated 
 pair , which we do not see in our simulations . further , the x - ray and 
 neutron scattering experiments also reported that the total number 
 of cl ions coordinating with mg ions 
 does not exceed 2 below 3 m concentrations . from the simulation results , however , we find the 
 number sof cl ions coordinating with the mg ions are 7 , 5 , 3.2 , and 3 for the four different concentrations 
 studied ( figure 1d ) , all being larger than 
 2 . thus , from these observations it becomes clear that the strength 
 of interactions between ions is overestimated in the simulations , 
 leading to aggregation . we observe similar behavior in the simulations 
 of cacl2 solutions ( figure 2 ) . osmotic 
 pressure values start to deviate for concentrations > 0.5 m ( figure 2a ) . we find a direct correspondence between the 
 occurrence of ion clusters and the increase in deviation of the calculated 
 osmotic pressure from the experimental values . visual inspection shows 
 aggregate formation for all concentrations ( figure 2c ) . through rdf calculations we find that the probability 
 of finding cl around ca ions within 
 a minimum radius of 1.98  is very high , indicating the presence 
 of direct ion pairs ( figure 2b ) . as seen in 
 mgcl2 simulations , the cacl2 simulations show 
 coordination numbers of cl with ca more than 2 for all concentrations ( figure 2d ) . the black line is the osmotic pressure data obtained 
 from experiments , and the magenta line 
 is the osmotic pressure obtained from simulations after using the 
 multisite ca model . ( b ) snapshot of a simulation of ca and cl ions in water showing aggregation 
 of the ions ( ca as magenta spheres , cl as green spheres , water molecules are not shown for clarity ) . ( c ) radial distribution 
 function between ca and cl ions at 
 different concentrations of cacl2 ( m represents 
 molality = mol / kg ) . ( d ) total number of cl ions 
 binding to ca ions at different concentrations of cacl2 . the results described above clearly 
 indicate that the reduction 
 in osmotic pressure at higher concentration occurred due to clustering 
 and aggregation of ions . aggregation of ions leads to fewer collisions 
 of the solute particles with the membrane and thus leads to a reduction 
 in the osmotic pressure . based on the experimental results , however , 
 such aggregation only starts to occur at solute concentrations nearing 
 the crystallization conditions . simulation studies 
 by callahan and co - workers have also shown that the existence of direct 
 ion pairing in mgcl2 solutions is not energetically favorable . this suggests that the ion aggregation that 
 we observe in solute concentrations ranging below 3 m is certainly 
 an artifact of the model . such an artifact can be attributed to the 
 discrepancies in the interaction parameters of the atoms . as mentioned 
 before , the parameters for the multisite model were developed and 
 optimized to reproduce accurate coordination properties of a single 
 ion in water and are insufficient to correct for the excessive cation we focus 
 on refining the parameters of the multisite model such that an accurate 
 physical behavior of concentrated ionic solutions can be reproduced 
 through simulations . to reduce 
 ion aggregation , we focus on changing the cation  anion interaction 
 parameters , which affect the electrostatic and the van der waals interactions . we keep the same charge distribution and lj parameters for the central 
 atoms , as these choices were found to be critical for calculation 
 of sensitive thermodynamic parameters in our previous work . in the calculations described above , the lj 
 parameters between the cation and anions were determined by the lorentz following the work of luo and roux , we move away from the conventional mixing rules 
 and allowed the two unique atoms of the cation to interact differently 
 with the anion . we modify the force field by scaling the repulsive 
 component of the lj interaction ( c12ij , where i = dummy atom and j = cl ) until the calculated osmotic 
 pressure reached the experimental value . a systematic survey of c12ij values was first carried 
 out at a concentration of 1.0 m to get the upper and lower bounds 
 for c12ij values ( figure 3a ) . we started with log10[c12ij ] 
 = 6 and continued the variation until we obtained an optimal 
 value of c12ij that could 
 reproduce experimental osmotic pressure for a wide range of concentrations . osmotic pressures were plotted for intermediate sets of lj parameters 
 to find the best fit ( intermediate and final profiles in figure 3b ) . a final value of log10[c12ij ] = 5.5 
 is able to match the experimental osmotic pressure for all concentrations 
 ( figure 3b ) . ( a ) linear regression fit of osmotic pressure 
 values for 1 m mgcl2 at different lj parameters for interaction 
 between dummy 
 atoms and cl ions ( c12 ) ; ( b ) osmotic pressure profiles against different concentrations 
 of mgcl2 for 3 different van der waals parameter sets . the results for the original ( saxena and sept ) are shown in magenta , 
 intermediate values are shown in yellow , and the final fit is shown 
 in blue . ( c ) radial distribution function between all atom pairs in 
 1 m mgcl2 solution after using the optimized parameters 
 ( solid red line ) . the profile is compared with the data obtained from 
 x - ray experiments ( dotted red line ) . after 
 obtaining the refined lj parameters for the dummy atoms interaction 
 with cl , we test the structural characteristics of the mgcl2 solutions by running simulations with the modified parameters . we 
 find that the new simulations of mgcl2 solutions do not 
 show any signs of spontaneous ion aggregations ( figure 4a ) . the first coordination shell of both mg and 
 cl ions contained water molecules , showing evidence 
 of solvent separated ion pair formation . the first peak in the radial 
 distribution function profile of mg and cl for 1 m mgcl2 solution is now observed at 4.5  
 ( figure 3c ) , in contrast to the previous value 
 of 1.95  ( figure 1c ) . anion 
 pair at a separation of 4.5  corresponds to an ssp and matches 
 the results from x - ray crystallization experiments ( figure 3c ) . we also find that the 
 new simulations show a total of 2 cl ions coordinating 
 with mg ions , for all concentrations ( below 3 m ) , matching 
 the observations from neutron scattering and x - ray experiments on 
 mgcl2 . thus , the structural characteristics of mg and cl ions in the simulations of mgcl2 with modified 
 parameters confirm that the improved parameters for the multisite 
 model are able to eliminate the ion aggregation artifact and can now 
 reproduce accurate physical behavior of the ions at physiological 
 conditions . to test the transferability of the new multisite 
 model parameters 
 to other ions , we calculate the osmotic pressure of the cacl2 solutions using the new c12ij values for dummy atoms and cl ions , obtained 
 from the refinement of mgcl2 solutions . keeping everything 
 else the same , we replaced the old lj parameters of the dummy atom 
 with the improved values , thereby improving the interaction of the 
 dummy atoms with cl ions . the osmotic pressure 
 values , obtained from the simulations of cacl2 solutions 
 with the new parameters , match well with the experimental values , 
 without the need of further refinement ( figure 4c ) . this suggests that the changes made in the interaction of the 
 dummy atoms with the anions were sufficient to reproduce accurate 
 physiological behavior of ca ions with cl ions . we further verified this result by calculating the structural 
 properties of cacl2 in solution . the rdf profile between 
 ca and cl ions shows the first peak 
 at an ion  ion separation distance of 5  . this value corresponds 
 to a solvent - separated ion pair formation in cacl2 solutions . also , the number of cl ions coordinating with ca ions is 2 for all simulated 
 concentrations thus , this result suggests that 
 the new parameters for the dummy atoms are sufficient to improve the 
 interaction of another divalent cation with the anion , without the 
 need for further customization . the unique parameters of the dummy 
 atoms could be universally used whether they are part of the mg or ca multisite model . ( a ) snapshot of a simulation 
 of mg and cl ions in water after refinement 
 of the mutisite model ( mg as red spheres , cl as green spheres , water occupying 
 first solvation shell of mg are shown in licorice ) . ( b ) radial distribution functions between 
 mg and cl ions at different mgcl2 concentrations , after parameter refinement . ( c ) osmotic pressure 
 obtained for different concentration of cacl2 after using 
 2 sets of van der waals parameters ( original ( magenta ) and modified 
 parameters ( blue ) ) . ( d ) radial distribution function between ca and cl ions at different cacl2 concentrations , after parameter refinement . next , we calculate how the modeling of concentrated ionic solutions 
 affect the water coordination and exchange rate for the cations . the 
 original parameters of the multisite model allowed accurate prediction 
 of the structure of the coordination shell , with the ion - oxygen distances , 
 and the coordination number matching the experiments . in addition , 
 the model also gave accurate values for water exchange rates , indicating 
 that they are able to accurately represent long distance interactions . in this section , we used the parameters of the 
 new mg and ca multisite model and performed 
 simulations in water with and without cl ions . the calculations without cl ions use identical 
 parameters as in the work of saxena and sept . the resulting ion - oxygen distances , 
 coordination number , and water 
 exchange rates are then calculated from these simulations ( table 3 ) . a value of 2.1  has been reported for mg  o 
 interaction and 2.4  for ca  o 
 interactions . these values were observed as 
 2.1 and 2.4  , respectively , in previous work . as expected , these values remain unchanged for dilute solutions 
 ( table 3 ) . the refined model further allows 
 us to calculate cation - oxygen interactions in the presence of counterions , 
 which was not possible in our earlier simulations due to the ion aggregation 
 artifact . we observed that the cation  oxygen interactions remain 
 unchanged in solutions with higher anion concentrations ( 1 m ) . the 
 second property we look at is the coordination number of mg and ca ions . these values remain unchanged at 6 and 
 7 , respectively , before and after the parameter refinement . next we 
 calculate the water exchange rates for different cations in concentrated 
 cacl2 and mgcl2 aqueous solutions . the exchange 
 rates are difficult to determine from experiments : therefore , the 
 available data covers a wide range of values or bounds . for mg ions , for ca ions , the rates are much 
 faster ( less than a nanosecond ) . saxena and sept showed that the multisite model significantly enhanced 
 the water exchange properties of the cations from the previous models . they obtained a value of 20 ps for ca ions , while for 
 mg ions the simulations were not long enough to capture 
 the microsecond exchange rates . in this work , we calculate the exchange 
 rates from simulations of mg and ca ions 
 in the presence of cl ions . we find that the exchange 
 rates for ca ions still remain in the ps range ( 50 
 ps ) but with a slight increase . the increase in the exchange rate 
 can be attributed to the presence of counterions in the simulations . for mg ions , the inner shell waters remain unchanged 
 over a 100 ns simulation , indicating that the exchange rates would 
 be greater than 100 ns . thus , these results suggest that the refined 
 multisite model preserved its interaction with water . the properties 
 of the single ions remained unchanged from the original model and 
 the refined model acquired an additional ability to interact with 
 the counterions at any concentration . the 
 multisite model for divalent ions ( mg and ca ) developed previously achieved 
 better solvation , thermodynamics , and coordination properties of divalent 
 cations in molecular simulations . however , it suffered with artifacts 
 of cation  anion pairing and unusual ion aggregation in highly 
 concentrated ion solutions . in this work , we refined the parameters 
 of the multisite model such that the physical behavior of highly concentrated 
 solutions is reproduced accurately . the final parameters produced 
 a good fit with the experimental osmotic pressure values for various 
 concentrations of mgcl2 solutions , and the common artifact 
 of ion aggregation was eliminated . further , the changes made were 
 directly transferable to other ions and when used in cacl2 solutions , and experimental values of osmotic pressure were successfully 
 reproduced . the original multisite model had presented a number of 
 advantages over standard spherical ion models in molecular simulations 
 of cation binding biomolecules . the refined model went one step 
 further in extending the range of calculations from neutralizing salt 
 concentrations to excess salt . taken together , the multisite cation 
 model now presents a complete model for the use of divalent ions in 
 classical simulations . the model allows for precise predictions of 
 conformational and thermodynamic changes occurring in divalent cation 
 binding biomolecules ( proteins and nucleic acids ) at physiologically 
 relevant concentrations , an environment that was challenging to obtain 
 previously in simulations of divalent ions . since we did not change 
 the interaction parameters for the cation with water and the biomolecule , 
 the properties calculated for biomolecules by saxena and sept will 
 still hold for the modified model described here . additionally , the 
 model is compatible with currently available standard force fields 
 and has been implemented in gromacs . the 
 application of this model will enable numerous future studies involving 
 cation binding to biomolecules without compromising accuracy .
accurate 
 force field parameters for ions are essential for meaningful 
 simulation studies of proteins and nucleic acids . 
 currently accepted 
 models of ions , especially for divalent ions , do not necessarily reproduce 
 the right physiological behavior of ca2 + and mg2 + ions . saxena and sept ( j. chem . 
 theor . 
 comput.2013 , 9 , 35383542 ) described a model , 
 called the multisite - ion model , where instead of treating the ions 
 as an isolated sphere , the charge was split into multiple sites with 
 partial charge . 
 this model provided accurate inner shell coordination 
 of the ion with biomolecules and predicted better free energies for 
 proteins and nucleic acids . here , 
 we expand and refine the multisite 
 model to describe the behavior of divalent ions in concentrated mgcl2 and cacl2 electrolyte solutions , eliminating the 
 unusual ion  ion pairing and clustering of ions which occurred 
 in the original model . 
 we calibrate and improve the parameters of 
 the multisite model by matching the osmotic pressure of concentrated 
 solutions of mgcl2 to the experimental values and then 
 use these parameters to test the behavior of cacl2 solutions . 
 we find that the concentrated solutions of both divalent ions exhibit 
 the experimentally observed behavior with correct osmotic pressure , 
 the presence of solvent separated ion pairs instead of direct ion 
 pairs , and no aggregation of ions . the improved multisite model for 
 ( mg2 + and ca2 + ) can be used in classical simulations 
 of biomolecules at physiologically relevant salt concentrations .