highly 
 charged ( trivalent and tetravalent ) metal ions are of great 
 interest in both supramolecular , biomolecular , and rare - earth chemistry . some of these serve as the coordination center that performs structural , 
 catalytic , or electron - transfer functions , while others are well - known 
 biotoxins . for example , iron  sulfur clusters have biochemical 
 functions involving the respiration and photosynthesis processes .
lanthanide series elements have attracted significant attention 
 due to their specific electromagnetic and optical characteristics .
their complexes serve as luminescent probes because 
 of their large stokes shifts and emission lifetimes . the actinide series elements , such as th , u , and pu , are 
 well known for their radioactivity . all 
 of these can form highly charged metal ions , which pose a major challenge 
 to computational modeling .
indeed , effective and accurate modeling 
 of these ions will give insight into separation and recycling process 
 aimed at reducing their harmfulness to the environment .
there 
 are several theoretical methods to model metal ions : quantum 
 mechanics ( qm ) , molecular mechanics ( mm ) , and the hybrid qm / mm method .
classical force fields , 
 which use an analytical function to represent the relationship between 
 the energy and configuration of a system , have a significant speed 
 advantage over the quantum - based methods .
it is the state - of - the - art 
 tool to study systems at the molecular level when combined with molecular 
 dynamics or monte carlo methods . for metal ions ,
there are several widely used models including the 
 bonded model , nonbonded model , and cationic dummy model .
the bonded model represents the 
 interaction between the ion and its surrounding residues via the bond , 
 angle , torsion , coulombic , and van der waals ( vdw ) terms . because 
 of the harmonic approximation used in the bonded model , it does not 
 simulate the processes involving chemical bond formation and dissociation .
the nonbonded model usually places an integer 
 charge on the metal ion and only uses the coulombic and vdw terms 
 to represent the intermolecular interactions between the metal ion 
 and surrounding particles .
this simplification can result in a notable 
 underestimation for modeling systems with strong covalent bonds .
the dummy cationic model places the charge between 
 the metal ions and the ligating atoms to mimic the covalent bond . besides the models previously discussed
, there 
 are also some polarizable force fields that have been developed for 
 metal ions in recent years . even though more accurate models exist , the 12 - 6 lennard - jones 
 nonbonded model is widely used due to its simple form , computational 
 efficiency , and excellent transferability characteristics .
however , in previous research , we found that the 12 - 6 lennard - jones 
 ( lj ) nonbonded model could not reproduce several experimental properties 
 across a series of divalent metal ions due to the neglect of the ion - induced 
 dipole interaction .
for the divalent 
 metal ions , on average , there is a 50 kcal / mol underestimation for 
 the hydration free energy ( hfe ) if we want to reproduce the experimental 
 ion 
oxygen distance ( iod ) values , while there is 0.27 
  reduction for the iod values if we reproduce the experimental 
 hfe values . in light of this
, we proposed a 12 - 6 - 4 lj - type nonbonded 
 model to account for the charge - induced dipole interaction . after picking suitable parameters ,
it was demonstrated 
 that it was possible to reproduce the experimental hfe , iod , and coordination 
 number ( cn ) values simultaneously for a series of divalent metal ions . 

furthermore , it was shown the new nonbonded model was readily transferable 
 to mixed systems such as salt solutions and nucleic acid systems . in the present work ,
we have estimated the 12 - 6 lj parameters for 
 24 highly charged metal ions ( 18 m(iii ) ions and 6 m(iv ) ions ) for 
 three widely used water models ( tip3p , spc / e , and tip4pew ) , respectively .
this illustrated that the 12 - 6 
 lj nonbonded model has a much larger underestimation of the ion 
water 
 interactions for the highly charged ions than for the mono and dications . 

we next parametrized the 12 - 6 - 4 lj - type nonbonded model for the 24 
 highly charged metal ions for the same three water models . in general , 
 these parameters simultaneously reproduce both the experimental hfe , 
 iod , and cn values with good accuracy .
this work opens up new opportunities to simulate m(iii ) and m(iv ) 
 ions in aqueous solution using classical models .
stable 
 trajectories were obtained with the metal binding site being well - conserved , 
 further supporting the excellent transferability of these parameters .
in present work , we employed the nonbonded 
 model in the amber force field1 in eq 1 , the uij(rij ) is the nonbonded interaction 
 potential between 
 atoms i and j at distance rij .
herein e represents the charge of 
 the proton , while qi and qj are the partial charge of 
 atoms i and j. the partial charge 
 of metal ions is always treated as an integer number in the nonbonded 
 model .
the vdw interaction part uses a 12 - 6 lennard - jones ( lj ) potential , 
 in which there are two parameters ( ij and rmin , ij ) 
 that need to be determined . using the geometric combining rule ,
the 
 well depth of the lj potential is2 for the rmin , ij parameters , 
 there are two widely used combining rules3and4 equation 3 is employed in the amber and charmm force 
 fields , while eq 4 is used in the opls force 
 field . in the present work
, we employed 
 the lorentz  berthelot combining rules , which is the union of 
 eqs 2 ( lorentz combining rule ) and 3 ( berthelot combining rule ) .
for the 12 - 6 - 4 lj - type 
 potential , we employed the following expression5 an r term was added to describe 
 the ion - induced dipole interaction , which can not be overlooked for 
 highly charged systems .
the parametrization work concentrated on the 
 determination of the rmin , i , i for the different metal 
 ions , and the c4 term between the metal ion and oxygen atoms of the different 
 water models . in the present work ,
we have determined the parameters 
 for three widely used water models ( tip3p , spc / e , and tip4pew ) independently .
previous work demonstrated that 
 it is necessary to design different parameters for these water models 
 due to their different geometries , charge distributions , and vdw parameters .
the thermodynamic integration ( ti ) 
 method was used to simulate the hfe values .
ti calculates the free - energy 
 change between two different states of a system .
an initial - state / final - state 
 mixing potential was used during the simulation , in which v0 and v1 represent 
 the potential of the initial and final state , respectively , while 
  governs the mixing between the two states .
k determines whether the mixing is linear or of higher order ( k = 1 is linear).6 the ion solvation process 
 is modeled 
 as transfer of an isolated metal ion from the gas to liquid phase . 
 in the present work
, we employed the thermodynamic cycle depicted 
 in figure 1 . to avoid the  end - point 
 catastrophe 
, we employed the linear scaling soft - core ti method to obtain the gvdw term.7 in eq 7 , rij is the 
 distance between the dummy atom and the surrounding 
 particles , while  is the distance at which the two particles 
 vdw interaction is equal to zero . here
the  is the well depth 
 and  is a constant set to 0.5 .
problem is largely avoided because there is limited energy penalty 
 induced between the dummy atom and proximal particles when the vdw 
 potential of the dummy atom is turned on.89 as seen from eq 8 , the free - energy of each 
 process is obtained via the integration of the derivative of the potential 
 with respect to . herein we employed gaussian quadrature ( eq 9 ) to evaluate the 
 integral in an efficient way .
as illustrated in figure 1 , we obtained 
 the hfe values based on the free - energy changes associated with four 
 processes : gvdw , gele+pol , gele+pol , and gvdw . herein the hfe value is computed using hfe = 1/2  ( gvdw + gele+pol  ( gele+pol  
 gvdw ) ) .
the gele+pol and gele+pol are gele and gele respectively , when the 12 - 6 nonbonded model 
 is employed . at first
, a dummy atom was placed in the center of a 
 cubic water box ( with size 32   32   
 32  ) with the closest water molecule 1.5  away 
 from it .
there are in total 722 water molecules in the tip3p or spc / e 
 water boxes , while there are 732 water molecules in the tip4pew water box . afterward
1000 steps of steepest descent minimization 
 were performed , followed by 1000 steps conjugate gradient minimization . 
 then , a 500 ps heating procedure was performed to heat the system 
 from 0 to 300 k in the nvt ensemble .
next , we equilibrated the system 
 for 500 ps at 300 k and 1 atm using the npt ensemble .
the final snapshot 
 from the equilibration simulation was used as the initial structure 
 for the calculation of gvdw . to 
 balance accuracy and speed
, we used the four - window linear soft - core 
 scaling process to obtain the gvdw value with  values of 0.1127 , 0.5 , 0.88729 , and 0.98 , respectively . 

the simulation of last window was used to further equilibrate the 
 system but was not used in the free - energy evaluation .
the final snapshot 
 was used to initiate the determination of gele+pol . for the gele+pol calculation process
, a nine - window linear scaling scheme was utilized 
 with  values of 0 , 0.2544 , 0.12923 , 0.29707 , 0.5 , 0.70292 , 
 0.87076 , 0.97455 , and 1 , respectively .
afterward , the gele+pol calculation process was carried out 
 in a similar manner . in these simulations ,
the first and last windows 
 (  equal to 0 and 1 ) are not used in the final free - energy calculation 
 but to further equilibrate the system .
finally , the gvdw simulation procedure was carried out using 
 a three - window linear soft - core scaling process in which  was 
 set at 0.1127 , 0.5 , and 0.88729 , respectively .
for the determination 
 of gvdw and gvdw , each window was simulated for 300 ps , with 
 the last 200 ps used for data collection .
the gele+pol and gele+pol simulations covered 200 ps , with the last 150 ps used for data collection 
 in each window .
all ti simulations were performed in the npt ensemble . we have employed two different methods to evaluate the uncertainty 
 of the simulated hfe values in the present work .
the first method ( set 1 ) divided each sampling segment ( gvdw , gele+pol , gele+pol , and gvdw ) into two even portions and estimated the 
 uncertainties for the vdw and electrostatic plus polarization free - energy 
 determinations .
for example , we have 200 ps of sampling for each window 
 for the determination of gvdw and 
 gvdw .
we used the first 
 100 ps of sampling to calculate the gvdw - part  1 and gvdw - part  1 values , while we used the later 
 100 ps of sampling to obtain the gvdw - part  2 and gvdw - part  2 .
then , we assessed the uncertainty of the vdw free - energy determination 
 by calculating the standard deviation based on these four values . 

the uncertainty of the electrostatic plus polarization free - energy 
 determination was obtained in a similar manner with each fragment 
 using 75 ps of sampling .
subsequently , we obtain the total uncertainty 
 in the hfe value by adding the uncertainties of the vdw and electrostatic 
 plus polarization free - energy determinations .
the second method 
 ( set 2 ) uses eq 10 , in 
 which the a is the autocorrelation 
 time of observable a while ( ac ) is 
 the standard deviation of a. t is 
 the sampling time of the simulation and a is 
 the final uncertainty of the observable a. the final 
 hfe values given in the spreadsheets given in the si are depicted as a  
 . again , the uncertainty of the vdw and electrostatic plus 
 polarization free - energy determinations were evaluated separately , 
 and the final uncertainty was treated as their sum .
herein we used 
 a a value equal to 500 fs for the 
 vdw free - energy determinations , while we used 250 fs for the electrostatic 
 plus polarization scaling part .
we obtained these values based on 
 test simulations , and they are consistent with previous work.10 the set 1 approach yields uncertainties in the range of 0.17.5 
 kcal / mol with an average of 1.2 kcal / mol , while set 2 gives 
 uncertainties in the range of 1.1 to 2.0 kcal / mol with an average 
 of 1.4 kcal / mol .
on the basis of these analyses , we estimate 
 that the hfe uncertainty is in the 2.0 kcal / mol range , which 
 is quite small given the magnitude of hfes we are computing .
a metal ion ( with an integer 
 + 3 or + 4 partial charge ) was solvated in the center of a cubic water 
 box ( with the same size as described in the hfe simulation part ) . 

then , 1000 steps of steepest descent and 1000 steps of conjugated 
 gradient minimization were carried out to relax the initial structure . 

simulation was performed 
 in the nvt ensemble that took the system from 0 to 300 k. next , 500 
 ps of equilibration , followed by 2 ns of sampling were performed at 
 300 k and 1 atm .
snapshots were stored every 0.5 ps ( every 500 steps 
 for 4000 snapshots in total ) for the subsequent iod and cn analysis . 

the radial distribution function ( rdf ) of the ion and water oxygen 
 atom was then calculated based on the average volume of the entire 
 trajectory in the range of 05.0  with a grid resolution 
 of 0.01  .
the first quadratic fit was performed using the points 
 within 0.1  of the first peak of rdf . in this way
the second 
 quadratic fitting was done based on the points within 0.1  
 of the apex obtained from the first fitting . in total , 21 points were 
 used for each fit .
the maximum given by the second fitting was treated 
 as the final iod value to two decimal places .
the cn value was obtained 
 by integrating from the origin to the first minimum of the rdf .
the amber 12 suite of programs was used 
 to perform the simulations , while the amber tools suite of programs was utilized to carry out the data analysis . 

the particle mesh ewald ( pme ) method and periodic boundary 
 condition ( pbc ) were employed throughout .
the time - step was 1 fs , 
 while the cutoff was set to 10  .
for the temperature control , 
 the langevin algorithm was utilized with a collision frequency equal 
 to 5.0 ps .
the isotropic pressure algorithm was 
 used to control the pressure . the pressure relaxation time was set 
 to 10 and 1 ps in the ti and standard md simulations , respectively . 

the hfe values 
 for all m(iii ) and 
 m(iv ) ions investigated were taken from marcus .
it is one of the most complete databases regarding the 
 thermodynamic properties of ions .
the iod and cn values for m(iii ) 
 ions were taken from marcus review , while the iod and cn values for the m(iv ) ions were taken from a 
 number of sources .
the experimental effective ionic radii were obtained from shannon . on the basis of the iod values and effective 
 ionic radii , we estimated the effective radii of the coordinated water 
 and display the data in table 1 .
some highly 
 charged ions that readily hydrolyze water such as as , 
 sn , and pb ions were not considered in the present work . referenced from marcus . referenced or calculated from 
 marcus . from 
 ankudinov
one is 
 for rmin/2 from 0.9 to 2.3  with 
 0.1  intervals without the c4 term 
 ( namely using the 12 - 6 lj nonbonded model ) , while the other is for 
 the same rmin/2 value sets with a constant c4 term equal to 500 kcal / mol .
the  values are obtained for each rmin/2 value based on the noble gas curve ( ngc ) , which 
 was previously developed .
the hfe , iod , 
 and cn values for each parameter point are collected in tables si.1 
 and si.2 in the supporting information ( si ) . 

these data points might be useful to those who want to parametrize 
 the 12 - 6 lj or 12 - 6 - 4 lj - type nonbonded model with different target 
 values from the present work .
on the basis of the quadratic 
 fitting of the data points from the parameter scans without the c4 term ( see the si ) , we estimated the hfe and iod parameter sets .
the two parameter 
 sets are shown in tables 2 and 3 , while the estimated absolute and percent errors are shown 
 in table si.3 .
similar to our parametrization 
 of divalent metal ions using the 12 - 6 lj nonbonded model , the rmin/2 parameters 
 in the hfe parameter set are in excellent agreement with the vdw radii 
 calculated using the quantum - mechanical scaling principle ( qmsp ) method . for the al , y , and 
 la ions , the estimated rmin/2 values in the estimated hfe parameter set for the tip3p water 
 model are 0.981 , 1.454 , and 1.628  respectively .
the calculated 
 vdw radii are 1.046 , 1.481 , and 1.642  , respectively , based 
 on the qmsp method .
there is only a 6.2 , 1.8 , and 0.9% difference 
 between these two sets of values .
moreover , the estimated rmin/2 values used in the 12 - 6 lj parameter sets 
 could be used as the vdw radii for resp charge - fitting procedures . 

for example , in the work of kuznetsov et al . , they used 1.4  as the vdw radius for the resp charge fitting 
 for both the fe and fe ions , while for the 
 iod parameter set the fe radius was determined to be 
 1.409  in our previous research .
herein , the iod parameter sets for the fe ion estimated 
 the radius as 1.386 , 1.386 , and 1.375 for the tip3p , spc / e , and tip4pew water models , respectively . combining the data with 
 previous work on m(ii ) metal ions , we summarized the absolute and 
 percent errors for the 12 - 6 lj nonbonded model for m(ii ) , m(iii ) , 
 and m(iv ) ions in tables 4 and 5 .
these results , taken as a whole , show that the underestimation 
 of the 12 - 6 lj nonbonded model increases dramatically as the charge 
 on the metal ion increases .
for example , for the tip3p water model , 
 the average absolute error goes from 50 kcal / mol for m(ii ) 
 ions to 80 kcal / mol for m(iii ) and 240 kcal / mol for 
 m(iv ) ions for the iod parameter set , while the average absolute error 
 of the iod values for the hfe parameter set increases from 0.27 
  for m(ii ) ions to 0.29  for m(iii ) ions and 0.58  
 for m(iv ) ions .
for some of the monovalent ions , it is possible to 
 reproduce both the experimental hfe and iod values at the same time ; this underestimation of the 12 - 6 lj nonbonded 
 model is pretty small and can almost be neglected . because there are 
 significant errors associated with the 12 - 6 lj nonbonded model for 
 highly charged ions , we did not carry out further refinement work 
 on the estimated 12 - 6 lj parameters because the resultant parameters 
 would be of limited usefulness .
after initial parameter 
 selection and subsequent fine - tuning , the final 12 - 6 - 4 parameters 
 were determined .
the final optimized 12 - 6 - 4 parameters are given in 
 table 6 while the simulated hfe , iod , and cn 
 values are shown in table si.4 .
these parameters 
 reproduce the experimental hfe values by 1 kcal / mol and the 
 iod values by 0.01  for the m(iii ) ions , while they reproduce 
 the hfe values by 2 kcal / mol and the iod values by 0.01 
  for the m(iv ) ions .
just as in the 12 - 6 - 4 parameter sets for 
 divalent metal ions , the rmin/2 terms 
 are similar between the three water models , while the c4 term for tip4pew water is generally larger 
 than for the other two water models for the same metal ion . this may 
 due to the smaller dipole of the tip4pew water model ( 2.32 
 d ) relative to the tip3p ( 2.35 d ) and spc / e ( 2.35 d ) water models . 

figure 2 shows the accuracy comparison between 
 the 12 - 6 - 4 parameter set and the 12 - 6 parameter sets for divalent , 
 trivalent , and tetravalent metal ions .
we can see that there is significant 
 improvement in the accuracy using the 12 - 6 - 4 parameter set , which 
 is able to reproduce the experimental hfe and iod values simultaneously . 
 while for the 12 - 6 lj nonbonded model , if you want to reproduce the 
 experimental hfe values , the error in the simulated iod values would 
 increase along with the formal charge of the metal ions .
vice versa , 
 if you simulate the iod values using the 12 - 6 lj nonbonded model , 
 the error of the calculated hfe would increase markedly with an increase 
 in the oxidation state of the metal ion in question .
( a ) hfe errors for the 
 12 - 6 iod and 12 - 6 - 4 parameter sets for m(ii ) , 
 m(iii ) , and m(iv ) metal ions .
( b ) iod errors for the 12 - 6 hfe and 
 12 - 6 - 4 parameter sets for the m(ii ) , m(iii ) , and m(iv ) metal ions .
these metal ions 
 have much stronger ion  water interactions than the ln ions .
they form a stable 
 octahedral structure with water molecules in the first solvation shell . 

data in table 1 indicate that the average effective 
 radius of the first solvation shell water is 1.35  for 
 the first several metal ions , which is consistent with strong interactions 
 between the coordinated water molecules and these metal ions .
these 
 values are close to previously proposed coordinated water radius ( 1.34 
  ) .
the corresponding average values 
 are 1.49  and 1.46  for the ln and the m(iv ) metal ions , respectively , which implies a smaller 
 electronic cloud overlap between the metal ion center and each of 
 the coordinated water molecules .
al , in , and tl are group iv ions . for the c4 parameters derived herein , we obtained a sequence of 
 tl > al > in .
the al ion is the smallest m(iii ) ion , resulting in a relatively 
 larger c4 term due to its strong covalent 
 interaction 
 with coordinated water molecules .
tl has two oxidation states , + 1 
 and + 3 , and the hfe values of the tl and k are almost the same in marcus hfe set .
the reduced electric potential of m(iii ) 
 + 3e = m is 1.67 , 0.3382 , and + 0.72 
 ev for al , in , and tl , respectively .
it readily obtains electrons from 
 its surroundings , which may be the reason for a strong charge - transfer 
 effect between the tl ion and the surrounding water molecules . 

the 12 - 6 - 4 parameters of in and tl ions 
 gave an excellent prediction for the hfe and iod values but overestimated 
 the cn value ( 8 instead of 6 ) , and this is mainly due to the lack 
 of a correction for the water  water interactions in the first 
 solvation shell during the simulations .
water interactions 
 were parametrized to reproduce the pure liquid water properties in 
 the original parameter design .
however , the first solvation shell 
 water molecules of the highly charged metal ions should more strongly 
 repel one another due to their bigger charge separations .
this effect 
 is smaller for m(i ) and m(ii ) metal ions , but it dramatically increases 
 for the highly charged ions .
meanwhile , this kind of effect may decrease 
 in protein systems due to the preorganization of the metal ion binding 
 sites .
fe has a larger c4 term 
 than al , the smallest m(iii ) ion , which suggests that 
 fe has a stronger interaction with its surrounding water 
 molecules .
this is consistent with quantum - mechanical charge - field 
 molecular dynamics ( qmcf - md ) simulations , which shows that the force 
 constant between the ion and the oxygen of first solvation shell water 
 molecules ( kion - o ) is 198 n / m for 
 fe compared with 185 n / m for al .
the fe ion has an average 
 1.85e charge ( from a mulliken analysis ) 
 in the qmcf simulation , while al ion has a corresponding 
 value of 2.5 e , which implies that there 
 is a stronger charge - transfer effect for the fe ion and 
 its surrounding water molecules than for the al ion . 

also , as previously discussed , this may be due to the underestimation 
 of the interactions between the first solvation water molecules . while 
 this effect is operative in aqueous solution , it will likely be less 
 of an issue in protein systems ( see discussion below ) .
the + 3 oxidation 
 state is the typical oxidation state of the ln elements , with the 
 exception that eu and ce could also be observed . 

this is because eu has a half - filled 4f orbital while 
 ce has the same electronic configuration as xe .
the interaction 
 of the ln ions with surrounding water molecules would 
 be expected to have more ionic character than the m(iii ) ions previously 
 discussed .
for example , the c4 terms between 
 the ln ions and water molecules are between 152 and 282 
 kcal / mol , which is smaller than the 258456 
 kcal / mol range seen for the other + 3 metal 
 ions previously discussed .
previous simulations found that the kion - o values are much smaller for the 
 ln ions ; for example , la , ce , lu , and er have kion - o values 110 n / m , while the values for the 
 al and fe ions are 185 and 203 n / m , respectively .
the ln ions have effective ionic 
 radii in the range of 0.86 to 1.03  and iod values in the range 
 of 2.34 to 2.55  .
these values are similar to that of the ca ion ( whose effective ionic radius and iod value is 1.00 and 2.46  , respectively ) .
therefore , they have been used as probes to investigate 
 the role of ca ions in biological systems . from table 6
, we 
 observe that the la and gd ions have the 
 smallest c4 terms among the ln ions
. this may be because they have either totally empty or half - filled 
 4f orbitals , making them more likely to form isolated ions , which 
 reduces the covalent character of their bonds with coordinated water 
 molecules .
the effective ion radius 
 decreases monotonically with an increase in the metal ions 
 atomic number due to the poor shielding of the 4f electrons toward 
 5s and 5p orbitals .
a similar tendency 
 can also be seen for the hfe and iod values along the series .
previous work reached the conclusion 
 that the lighter ln ions ( la to nd ) prefer a cn of 9 and the heavier ions ( gd to 
 tb ) prefer a cn of 8 , while the middle ions such 
 as sm and eu have cn between these two values .
it was proposed that the former ln ions have a tricapped 
 trigonal prism structure , which then shifts to a distorted bicapped 
 trigonal prism structure for the heavier elements as one of the two 
 capping water molecules leaves the first solvation shell .
generally speaking , there is a good agreement 
 between the 12 - 6 - 4 parameters for the hfe , iod , and cn values with 
 experiment , with the exception that some of the cn values were slightly 
 overestimated .
the final parameters gave a cn in the range of 9 to 
 10 for ln ions rather than the range of 8 to 9 reported 
 in the literature .
as previously discussed , this may be due to the 
 fact that there is no water  water interaction correction term 
 in the present parametrization process .
moreover , the cn values given 
 by marcus ( as shown in table 1 ) likely also 
 vary under different experimental conditions ( counterions used solute 
 concentration , etc . ) . among all ions ,
al , y , and la have the same electron configurations as the 
 noble gas atoms ne , kr , and xe , respectively . using the tip3p water 
 model as an example
, we can see their c4 values decrease from 399 to 216 and 152 kcal / mol , respectively .
there are only 
 a few m(iv ) ions that exist in aqueous solution , while the others 
 are readily hydrolyzed into polynuclear complexes in water .
the cn values 
 of these metal ions are greater than 8 , with some of them being 10 
 according to experiment .
previous work found that pu(iv ) , th(iv ) , and u(iv ) 
 could strongly bind to transferrin , an iron - transport protein .
hence , the parameters developed herein might 
 facilitate theoretical research on the biotoxicity of these m(iv ) 
 ions .
the zr and hf are in the ivb 
 group . even though hf has a larger atomic number than 
 zr , due to the  lanthanide contraction  

effect , it has a smaller effective radius , smaller iod , smaller hfe , 
 and a bigger c4 term than zr .
in contrast , ce and th are 
 in the same group where the larger atomic number ( th ) 
 has the bigger ionic radius , bigger hfe , and smaller c4 terms . this may be because they share the same electronic 
 structure as xe and rn , respectively . besides these ions , zr is another m(iv ) ion that has a noble gas atom s electronic 
 structure ( the same as kr ) .
there is also a single trend for the c4 terms of zr , ce , 
 and th ions for each specific water model .
for instance , 
 the c4 terms for the tip3p water model 
 are 761 , 706 , and 512 kcal / mol , respectively .
th , u , and pu are in the an series and are the largest elements 
 investigated in the present work .
canaval et al . investigated th in aqueous solutions using the qmcf - md method .
they found 
 a stable nine - coordinate complex , and even the third layer of water 
 molecules has a bigger mean residence time than that of pure water , 
 implying they are stabilized by the highly charged th ion .
u fluoresces due 
 to the electron transition between the 6d5f and 5f electronic configurations .
the u ion has the largest c4 term of all of the m(iv ) metal ions investigated .
frick 
 et al . investigated the u ion in aqueous solution using 
 the qmcf - md method , and the cn value was characterized as 9 , while 
 the average charge of u was predicted to be + 2.68 from 
 mulliken population analyses .
simulated the pu , pu , puo2 , and puo2 ions in water solution 
 using the car  parrinello molecular dynamics method .
they predicted that the pka value for the first hydrolysis step for the pu , pu , puo2 , and puo2 ions is 6.65 , 0.17 , 9.51 , and 5.70 , respectively , showing 
 a general tendency that the larger the charge of the metal center , 
 the lower the pka value of the first hydrolysis 
 reaction .
hf , zr , and pu have 
 relatively smaller iod values among the tetravalent ions , where they 
 all have experimental cns of 8 .
ce was determined to have an experimental cn of 9 , while u and th have 
 cns between 9 and 11 .
proposed that counterions also play a key role in the first solvation 
 shell structure , while the 9- , 10- , or 11-coordinated th have very small energy differences and are in a dynamic equilibration .
the simulated hfe and iod values of the 12 - 6 - 4 
 parameter set are in excellent agreement with the experiment .
the 
 simulated cn values of most of the m(iv ) ions are 10 , with 
 hf having a cn 8 for the tip4pew and 
 spc / e water models . herein , the tip3p model always predicted a larger 
 cn value than the other two water models , which may be because it 
 has a smaller c12 term ( 582.0  10 kcal / mol ) than that of the spc / e ( 629.4 
  10 kcal / mol ) and tip4pew ( 656.1  10 kcal / mol ) models .
later , we analyze several 
 redox pairs to explore the consistency of the c4 parameters we determined with respect to the behavior of 
 these pairs in aqueous solution .
we also calculated the relative hfe 
 between each redox pair for the tip3p water model ( see below ) .
the 
 divalent metal ions 12 - 6 - 4 parameters are from previous work . a nine - windows ti simulation ( 50 ps of equilibration 
 and 150 ps of sampling for each window ) was performed forward and 
 backward , respectively , to obtain the final results .
for example , the 
 simulated relative hfes of fe / fe , cr / cr , and ce / ce ion 
 pairs were 580.2 , 516.9 , and 698.2 kcal / mol , while the experimental 
 values are 579.6 , 516.2 , and 697.6 kcal / mol , respectively .
the experimental 
 iod values shrink 0.08  from 2.11  of fe to 2.03  of fe ion .
n / m for fe , which is almost twice as strong as that of 
 fe ( 93 n / m ) , while the effective charges ( from a mulliken 
 population analysis ) for the fe ion are in the range 
 of 1.25 to 1.45 ( with an average of 1.36 ) , and for the fe ion it is in the range of 1.70 to 1.95 ( with an average of 1.85 ) .
fe / fe redox pairs 
 exist broadly in biologically related system such as the fe  s 
 proteins and heme structures .
moreover , fe / fe redox pairs play fundamental roles in many electron - transfer 
 processes .
we have determined the 12 - 6 - 4 parameters for fe ion with three different water models in previous work . for the tip3p water model ,
the final optimized 
 parameters were rmin/2 = 1.457  , 
  = 0.02710805 kcal / mol , and c4 = 
 163 kcal / mol . from the 12 - 6 - 4 parameters determined for the fe ion 

herein , we find that the rmin/2 decreases 
 slightly as the outer shell electron number decreases , while the c4 term increases by 2.5 times relative 
 to fe .
this is consistent with a ratio between the c4 terms for a trivalent and divalent ion of 
 [ 3/2 ] = 2.25 , which is derived from the original ion - induced 
 dipole equation ( eq 11 in a prior publication ) .
the cr and cr ions 
 have [ ar]3d and [ ar]3d electronic structures , 
 respectively , where the cr(h2o)6 complex has a strong jahn  teller effect while the cr(h2o)6 molecule has a standard octahedral 
 configuration .
the iod values decreases from 2.08 to 1.96  for 
 the cr and cr ions . using the tip3p water 
 model as a representative example , we observe that the rmin/2 parameter decreases from 1.431 to1.405  , while 
 the c4 term increases from 137 
 to 258 kcal / mol for the cr and cr ions ( the cr parameters are reported 
 in reference ( 18 ) ) .
ce is the most stable 
 state because it shares the same electronic configuration with xe . 
 just like the two oxidation pairs previously discussed we find that 
 the rmin/2 value decreases while the c4 term increases significantly with increasing 
 charge .
for example , for the parameters determined for the tip3p water 
 model , the rmin/2 decreased 0.03 
  , while the c4 term increased by 
 480 kcal / mol .
pdb entry 4bv1 was used to obtain 
 the starting coordinates for this modeling exercise .
it is a protein tetramer with each monomer having a metal site 
 containing an fe ion .
the structure has been determined 
 by using x - ray crystallography to a resolution of 1.90  .
the 
 tetramer structure is shown in figure 3 , while 
 chain c with its metal site is shown in figure 4 .
the metal site contains four histidine groups , one cysteine group , 
 and one water molecule . by treating chain c as the initial structure
, 
 we performed three simulations with different parameter sets ( the 
 hfe , iod , and 12 - 6 - 4 parameter sets ) .
for the 12 - 6 - 4 parameter set , 
 the c4 terms between the fe ion and atom types other than water oxygen were evaluated using 
 eq 11 .
details of the simulation procedures and the polarizability 
 of each atom type are given in the si.11 a total of 10 ns of 
 sampling was performed during the simulation , and snapshots were stored 
 after each 500 fs .
the hfe parameter set prefers a smaller cn ( of 
 4 ) , and the metal ion moves out from the binding pocket , while stable 
 metal complex structures were obtained for the simulations using the 
 iod and 12 - 6 - 4 parameter sets .
an rmsd analysis was performed over 
 the heavy atoms of the backbone and the metal site for the simulations 
 by treating the initial structure ( experimental structure ) as reference . 

the rmsd 
 of the heavy backbone atoms fluctuated around 1.2  , 
 while the rmsd of the metal site was 0.5  .
these values 
 illustrate that the metal binding site is stable during the course 
 of the simulations .
we have also performed an rmsf analysis 
 of the backbone heavy atoms 
 for each residue together with the oxygen atom in the metal site binding 
 water .
it can be seen that the metal site residues : residue his10 
 ( residue number 11 ) , his 35 ( residue number 36 ) , his 41 ( residue number 
 42 ) , cys 97 ( residue number 98 ) , and his 100 ( residue number 101 ) , 
 all have relatively small rmsf values ( 0.5  ) . the metal 
 site binding water ( residue number 115 , is not shown in the figure 
 since the protein ends at residue number 112 ) has a rmsf of 0.6 
  for both the simulations with the iod and the 12 - 6 - 4 parameter 
 sets .
these results further validated that the metal ion site is stable 
 over the course of the simulations .
water molecules are not shown in the 
 figure , the ferric ion is 
 shown as a sliver sphere .
chain c in pdb entry 4bv1 ( left ) and a close up of the metal site in chain c 
 ( right ) .
the ferric ion is represented as a silver sphere , and it 
 is coordinated by one cys , four his , and one water molecule .
rmsd of heavy atoms of backbone ( left ) and the metal site ( right , 
 including the binding water molecule ) for simulations with iod and 
 12 - 6 - 4 parameter sets using the initial structure ( experimental structure ) 
 as reference .
rmsf of heavy atoms of 
 the protein residues in the simulations 
 using iod and 12 - 6 - 4 parameter sets .
in this work , we have estimated the 12 - 6 
 lj parameters and determined 
 the 12 - 6 - 4 lj - type parameters for 24 highly charged metal ions ( 18 
 m(iii ) ions and 6 m(iv ) ions ) with three water models ( tip3p , spc / e 
 and tip4pew ) based on a parameter scanning protocol .
we 
 have shown that with the increasing charge of the metal ions there 
 is a notable decrease in the accuracy of the 12 - 6 lj nonbonded model . 
 using tip3p as an example
, the average underestimation of the hfe 
 values increases from 50 kcal / mol for m(ii ) ions to 80 
 kcal / mol for m(iii ) ions and to 240 kcal / mol for m(iv ) ions 
 when trying to reproduce the experimental iod values .
the average 
 underestimation of the iod values increases from 0.27 to 0.29 
  and 0.58  for the m(ii ) , m(iii ) , and m(iv ) ions , 
 respectively , when trying to reproduce the experimental hfe values .
the 12 - 6 - 4 lj - type nonbonded model , which we previously described , 
 addresses this problem in a consistent manner .
it improves the accuracy 
 of the 12 - 6 lj nonbonded model remarkably with just a slight increase 
 in computational cost .
this parameter set , derived in this work , reproduced 
 several experimental values ( hfe , iod , and cn ) with good accuracy . 

they reproduce the hfe within 1 kcal / mol for the m(iii ) ions 
 and 2 kcal / mol for m(iv ) ions while reproducing the experimental 
 iod values to within 0.01  .
moreover , excellent quantitative 
 and qualitative agreement with previous experimental and computational 
 work supports the validity of the 12 - 6 - 4 lj - type nonbonded model . 

testing in a protein system also revealed good transferability of 
 the parameters determined herein .
 highly charged metal ions act as 
 catalytic centers and structural 
 elements in a broad range of chemical complexes . the nonbonded model 
 for metal ions is extensively used in molecular simulations due to 
 its simple form , computational speed , and transferability . 
 we have 
 proposed and parametrized a 12 - 6 - 4 lj ( lennard - jones)-type nonbonded 
 model for divalent metal ions in previous work , which showed a marked 
 improvement over the 12 - 6 lj nonbonded model . in the present study , 
 by treating the experimental hydration free energies and ion  
 oxygen 
 distances of the first solvation shell as targets for our parametrization , 
 we evaluated 12 - 6 lj parameters for 18 m(iii ) and 6 m(iv ) metal ions 
 for three widely used water models ( tip3p , spc / e , and tip4pew ) . 
 as expected , the interaction energy underestimation of the 12 - 6 
 lj nonbonded model increases dramatically for the highly charged metal 
 ions . 
 we then parametrized the 12 - 6 - 4 lj - type nonbonded model for 
 these metal ions with the three water models . 
 the final parameters 
 reproduced the target values with good accuracy , which is consistent 
 with our previous experience using this potential . 
 finally , tests 
 were performed on a protein system , and the obtained results validate 
 the transferability of these nonbonded model parameters . 
Introduction
Method
Results and Discussion
Conclusions