mean - field 
 methods provide a way to coarse grain the electrostatic 
 interactions between the solvent and the biomolecules of interest . they are commonly used to generate the electrostatic 
 potentials needed for biomolecular diffusion simulations and are used in calculations of the solvation free energies of biomolecules . the hundred - year - old poisson  boltzmann equation ( pbe ) 
 is 
 currently the flagship of the mean - field methods . despite the 
 success of the pbe , it is built upon several approximations ; 
 for example , the ions in the pbe are considered infinitesimally small . this approximation runs into trouble near highly charged biomolecular 
 surfaces , such as those of typical enzyme active sites . the reason is that , under the strong electrostatic 
 attractions near highly charged surfaces , the concentration of the 
 counterion can exceed its maximally allowed value without the restriction 
 from steric hindrance ; this will then lead to an overestimation of 
 electrostatic screening and an underestimation of the electrostatic 
 potential exerted by the biomolecule . to remedy 
 the lack of ion steric hindrance , or ion - size effects , in the pbe , 
 a stern layer , or ion - exclusion layer , outside of the molecular surface 
 has been used , within which no ion is allowed . to improve the pbe , the size - modified poisson boltzmann 
 equation ( smpbe ) , which incorporates the finite ion sizes through 
 a more physical lattice gas formulation , was developed and has been 
 applied to biomolecular systems . the computational cost of the smpbe is comparable to that of the 
 pbe when there are only two ion sizes . however , in its generalized form that can handle an arbitrary number 
 of different ion sizes , the smpbe takes much longer to solve . it is known that mean - field methods are parameter - dependent 
 and 
 that there is no standard way of choosing input parameters that guarantee 
 the optimum accuracy of prediction . the smpbe is reported to reproduce 
 experimental and molecular dynamics simulation results better than the pbe , but there still 
 lacks a comprehensive comparison between these two mean - field methods 
 using a variety of input parameters to determine whether the smpbe 
 generally outperforms the pbe . this paper discusses the accuracy 
 of prediction and the cost - efficiency 
 of the smpbe with respect to the pbe using several different parameter 
 sets . given that the mean - field methods are parameter - dependent , this 
 work also discusses the existence of a single parameter set with which 
 the mean - field methods give consistent high accuracy of prediction 
 across different ionic strengths and surface charge densities . specifically , 
 the equilibrium ion distributions outside of neutral and charged lipid 
 bilayers are calculated by the pbe and the smpbe using eight different 
 parameter sets and are compared to the molecular dynamics ( md ) simulation 
 results . the long - term goal of this work is to identify appropriate 
 tools for modeling diffusion dynamics of molecules that are of biological 
 or medical interest . in these diffusion processes , the system evolves from a nonequilibrium state to the equilibrium 
 state , and this work is an initial step to test the ability of the 
 smpbe in predicting ion distributions in the equilibrium state . lipid 
 bilayers are considered here because of their profound importance 
 in the activity of membrane proteins and cross - membrane signaling 
 pathways . the smpbe 
 formulation used here is capable of calculating the electrostatic 
 potential and the ion distributions for molecular systems containing 
 an arbitrary number of ion species with different ion sizes . a brief description of the formulation is included in the supporting information ( si ) section i. this formulation 
 is an extension of the earlier foundational smpbe theory by borukhov 
 et al . , which , in its original form , 
 can only handle one same - size cation  anion pair . boltzmann 
 solver ( apbs ) and is freely available 
 upon request . a brief explanation of the smpbe routine in apbs is 
 included in section ii of the si . molecular 
 dynamics ( md ) simulations can generate the equilibrium ion distributions 
 in a molecular system at higher accuracy and resolution than the mean - field 
 methods but also at significantly greater computational expense . in this work , the md simulation results are used 
 as a standard to assess the accuracy of the pbe and the smpbe predictions . the dimensions of the 
 lipid bilayer systems used in the md simulations 
 and the mean - field calculations are illustrated in figure s1 of the si . in this work , the bulk salt concentration 
 is the averaged salt ( nacl , kcl ) concentration at |z| = 60  ( z , the perpendicular distance to 
 the bilayer center , is defined in figure s1 of the si ) over the converged portion of the md trajectories . 1-palmitoyl-2-oleoylphosphatidylcholine 
 ( popc ) is a partially charged , 
 but overall neutral , lipid molecule , and each 1-palmitoyl-2-oleoylphosphatidylserine 
 ( pops ) lipid molecule carries one negative charge . all of the systems 
 in table 1 are built by the online application 
 charmm - gui membrane builder . the md simulations 
 are performed by namd2.9 software package with charmm36 lipid force field using mostly the default system settings offered 
 by the charmm - gui membrane builder ( see section iii of the si ) . each system in table 1 is simulated for 100150 ns . all of the popc simulations 
 are 150 ns long and reach convergence after 20 ns . all of the pops 
 simulations are 100 ns long and converge after 30 ns . the converged 
 md trajectories are utilized for all of the data analysis ( see section 
 iii of the si for the convergence criteria ) . to test 
 the accuracy of a mean - field method in predicting the ion distributions 
 in an md system listed in table 1 , three md 
 frames are evenly sampled from the converged md trajectory , and a 
 mean - field calculation is performed for each frame ; the final mean - field 
 results are averaged over all three calculations . note that apbs has 
 no periodic boundary conditions in its finite difference routine . 
 to avoid boundary effects , the lipid bilayer structures obtained from 
 the md simulations are manually extended 0.25 times along the positive 
 and negative x , y - axis , while the 
 mean - field results are only collected above and below the original 
 size of the lipid bilayer as in the md simulations , see figure s1 
 of the si . a 161(x ) by 
 161(y ) by 193(z ) finite difference 
 grid is used for all of the mean - field calculations following the 
 apbs electrostatic focusing scheme , and 
 the grid spacing is 0.43 (x ) by 0.43 (y ) by 0.47 (z ) , a reasonable resolution 
 for calculating ion distributions . the lipid atomic charges and radii 
 are taken from charmm36 charges and charmm36 van der waals radii . the 
 majority of the input parameters for the mean - field methods follow 
 the default settings of the apbs input preparation program pdb2pqr except for the ones mentioned in the next section . see section iv of the si for an example 
 of the apbs input file . the use of different mean - field method input parameters 
 is key 
 to the discussion of the accuracy of predictions of the smpbe with 
 respect to those of the pbe . among all of the mean - field calculation first , the ion size is chosen mainly because it is a particularly 
 important input parameter for the smpbe , and it is also used to define 
 the size of the stern layer . second , the definition of molecular surface 
 is chosen for its direct relationship to ion accessibility , thus the 
 ion distributions . specifically , four definitions of molecular surface , 
 figure 1 , and two sets of ion sizes , table 2 , are tested , yielding a permutation of eight different 
 sets of input parameters . vdws : the 
 van der waals surface , 
 drawn around the lipid atom centers with the charmm36 vdw radii , i.e. , 
 the lennard - jones rmin/2 values . fps : 
 the first - peak surface , drawn according to the md cation  lipid 
 carbonyl oxygen radial distribution function ( rdf ) first - peak positions . sas : the solvent - accessible surface , drawn by the center of the water 
 probe of radius 1.4  rolling across 
 vdws . ias : the ion - accessible surface , drawn in the same way as the 
 sas , only that the probe radius is that of the counterion . the dielectric 
 boundaries in all of the mean - field calculations are the same and 
 are defined by the solid light - blue line ( or equivalently the solvent - exclusion 
 surface ) . the van der waals ( vdw ) radii 
 are taken from the charmm36 ion vdw radii . the radial distribution 
 function ( rdf ) radii are the first - peak positions of the ion first , despite the multiple 
 definitions of the molecular surface in figure 1 , only the van der waals surface ( vdws ) is commonly recognized as 
 the molecular surface . second , the reason for the flexible use of 
 the term molecular surface in this work is as follows . supposedly , 
 the smpbe , without the help of the stern layer , is able to predict 
 solvation properties comparably accurately or more accurately than 
 the pbe with a stern layer . however , as it is shown later in the results section , unexpectedly , using the smpbe outside 
 of the vdws without a stern layer leads to large ion concentration 
 overestimation . to carry out further investigations , several alternative 
 molecular surfaces that are further away from the lipid molecules 
 than the vdws are tested . finally , the seemingly peculiar use of the 
 first - peak surface ( fps ) is inspired by the fact that the md radial 
 distribution function ( rdf ) first - peak positions are often smaller 
 than the sum of the vdw radii of the cation and its association partner , 
 the lipid carbonyl oxygen . the vdws is expanded by 
 0.6  for the na fps and 0.9  for the k fps . using the fps in the mean - field methods however , this turns out not to be 
 the case as shown in the results section . the rdf radius is chosen because it was used in a previous work , and an alternative set of ion radii is needed 
 to show how the pbe and the smpbe results change with respect to the 
 change of ion sizes . figure 2 shows an array of calculations on the number of 
 na ions bound per popc molecule . the corresponding results 
 for the k ion are shown in the si ( figure s2 ) . the number of ions bound is calculated by integrating 
 the concentration of ions within 25  to the lipid bilayer center . an example of the detailed position - dependent ion concentration is 
 shown in the si ( figure s3 ) . number of na bound per popc molecule at three nacl 
 bulk concentrations ( systems 13 in table 1 ) . boltzmann 
 equation ( without stern layer ) ; pbes : pbe with stern layer ; smpbe : 
 size - modified poisson  boltzmann equation ( without stern layer ) . subfigures ( a)(h ) use eight different parameter sets ; each 
 is a combination of a molecular surface ( vdws , fps , sas , ias , see 
 figure 1 ) and an ion radius set ( vdw radius , 
 rdf radius , see table 2 ) . ( g ) turns out to 
 be identical to ( e ) because the na vdw radius happens 
 to be the same as the radius of water , 1.4  , and so the sas 
 in ( e ) is exactly the same as the ias in ( g ) . by a crude observation , the pbe and the smpbe reproduce md 
 results 
 the best in figure 2 g , while the pbes ( the 
 pbe with stern layer ; the size of the stern layer is always set to 
 be the radius of the cation ) performs the best in figure 2a . however , the pbes calculation in figure 2a is actually the same as the pbe calculation in 
 figure 2 g . this is because the ion - accessible 
 surface ( ias ) used in figure 2 g exactly overlaps 
 with the outer surface of the stern layer in figure 2a . also , it should be reemphasized that the vdws is the physical 
 molecular surface , and it is only fair to compare the smpbe results 
 without the use of a stern layer to the pbe results with the use of 
 a stern layer . with this said , it can be concluded from figure 2a and b that when the smpbe is used outside of the 
 physical molecular surface without a stern layer , it significantly 
 overestimates the md ion concentrations compared to the pbe calculation 
 with the aid of a stern layer ; the size - modified correction of the 
 pbe does not seem to be as effective as the stern - layer - modified correction . this is a surprising result considering the fact that the popc bilayer 
 has zero net surface charge and yet the capability of the smpbe is 
 not enough to limit the ion concentrations to the correct range . furthermore , 
 caution should be taken to distinguish the real source of the corrections 
 or errors . given the poor performance of the smpbe in figure 2a , it is clear that the reproduction of the md results 
 by the smpbe in figure 2 g is mainly due to 
 the introduction of an effective stern layer through the use of the 
 ias instead of the smpbe formulation itself . the reason for the underestimation 
 of the md results by the pbes in figure 2 g 
 is the redundant use of a second stern layer on top of the ias . fortunately , as expected , at least the use of the smpbe consistently 
 alleviates the overestimation of ion concentrations by the pbe . the 
 strength of smpbe correction increases as the molecular surface shrinks 
 ( ias to vdws ) , which causes ions to experience a stronger electrostatic 
 potential , and as the ion radii increase ( vdw radius to rdf radius ) , 
 which results in more steric exclusion between ions in the smpbe formulation . however , in general , the smpbe gives similar results to the pbe ( figure 2a , c , e , f , g , and h ) . a previous study has shown 
 that the electrostatic free energies predicted by the pbe and the 
 smpbe can be quite different , but in 
 this work , the pbe and the smpbe electrostatic potentials are shown 
 to also be very similar under most parameter sets , see figure s4 of 
 the si . as mentioned in the introduction , the smpbe has a comparable computational cost 
 to the pbe when two ion sizes are used , but the cost significantly increases when the smpbe is generalized 
 to handle an arbitrary number of different ion sizes because the smpbe 
 can no longer be expressed analytically . if the simple newton s method is used to solve the generalized 
 smpbe , then the smpbe is up to hundreds of times slower than the pbe . 
 without developing a more efficient smpbe solver , only in the two - ion - size 
 case is smpbe a cost - efficient choice for its degree of correction 
 on the pbe . however , even with a faster smpbe solver , the stern layer 
 remains a more effective correction than the current smpbe formulation 
 at no extra computational cost . from figure 2 , it can be said that the vdws 
 combined with a stern layer of the size of the vdw radius of the cation 
 is the parameter set that gives the best pbe prediction of ion distributions . however , figure 3 shows that this does not 
 hold true for a different lipid bilayer surface charge density ( the 
 corresponding figure for the k ion is figure s5 of the si ) . number of na bound per lipid as 
 a function of the effective 
 stern layer thickness , the minimum distance between the vdws and the 
 molecular surface used for each calculation . ( a ) nacl of bulk concentration 
 0.13 m outside of the popc bilayer , system 2 in table 1 . ( b ) nacl of bulk concentration 0.17 m outside of the pops 
 bilayer , system 7 in table 1 ; vdw r and rdf 
 r are shorthand for vdw radius and rdf radius . gray shades mark where 
 the mean - field calculations using the vdw radius agree with md . the surface charge density of 
 the simulated pops bilayer is about 
 0.28 c / m or 1.73 e / nm ( for 
 reference , the surface charge density of dna is about 1 e / nm ) , a fairly high value . the gray shade positions in 
 figure 3a and b show that the thickness of 
 the stern layer for the pops bilayer needs to be smaller than what 
 is used for the popc bilayer ( the vdw radius of na ) for 
 the pbe to reproduce the md ion concentrations . figure 3 reveals that the same parameter set does not give consistent 
 accuracy of prediction when the surface charge density changes . this 
 is also true , although not as significantly , when the bulk ionic strength , 
 or the bulk ion concentration , varies . when looked at closely , the 
 slopes of the pbe , the pbes , and the smpbe plots in figure 2 are always equal to or smaller than that of the 
 md , leading to inconsistent performances of a parameter set across 
 different ionic strengths additionally , the same parameter set does 
 not give the same performance for different ions . for example , for 
 k , the ias with vdw radius does not perform as well as 
 it does for na , that is , the gray shade in figure s5 of 
 the si is not exactly overlapping with 
 the ias + vdw radius data point like it does in figure 3 . the inconsistency of the performance of a parameter set 
 is partly because mean - field methods are much more simplified than 
 the md simulations and do not describe several ionic effects such 
 as ion  ion correlations and fluctuations ; they also ignore the ion lipid specific interactions 
 modeled in the md force fields . the parameter 
 inconsistency can be circumvented if empirical guidelines 
 are available on how to choose parameters on the basis of the characteristics 
 of a system . for example , one can summarize a relationship between 
 the surface charge density and the effective stern layer thickness . two data points in this relationship are generated here : at zero surface 
 charge density , the vdw radius can be used as the size of the stern 
 layer , and at surface charge density 1.73 e / nm , roughly 0.8  for na and 0.6  for 
 k should be used as the size of the stern layer . more 
 data points need to be generated for a more complete knowledge of 
 this relationship , and the lipid - type specific ion absorption effects 
 should be considered . last but not least , even though theoretically 
 the pbe and the smpbe 
 become more distinguishable under a stronger electrostatic potential , 
 the strong electrostatic potential exerted by the higher surface charge 
 density of the pops bilayer is not yet enough to significantly separate 
 the pbe and the smpbe results when the common vdw ion radii are used , 
 figure 3b , putting the capability of the smpbe 
 under further doubt . first , 
 when it comes to studying ion equilibrium distributions , 
 the smpbe only offers a slight advantage over the pbe for the systems 
 considered here . it is neither faster nor more accurate than the simple 
 stern layer ion - size correction of the pbe . there may , however , be 
 other physical regimes where the smpbe has a strong advantage . further 
 improvements are needed for the smpbe to be an impactful correction 
 to the pbe . second , no single parameter set gives consistent mean - field 
 method accuracy of prediction across systems of varying surface charge 
 densities and ionic strengths . empirical rules may be able to be summarized 
 to guide the choice of parameters of the mean - field calculations .
mean - field 
 methods , such as the poisson  boltzmann equation 
 ( pbe ) , are often used to calculate the electrostatic properties of 
 molecular systems . in the past two decades , an enhancement of the 
 pbe , the size - modified poisson 
 boltzmann equation ( smpbe ) , 
 has been reported . here 
 , the pbe and the smpbe are reevaluated for 
 realistic molecular systems , namely , lipid bilayers , under eight different 
 sets of input parameters . 
 the smpbe appears to reproduce the molecular 
 dynamics simulation results better than the pbe only under specific 
 parameter sets , but in general , it performs no better than the stern 
 layer correction of the pbe . 
 these results emphasize the need for 
 careful discussions of the accuracy of mean - field calculations on 
 realistic systems with respect to the choice of parameters and call 
 for reconsideration of the cost - efficiency and the significance of 
 the current smpbe formulation .