Effective particle size from molecular dynamics simulations in fluids

We report molecular dynamics simulations designed to investigate the effective size of colloidal particles suspended in a fluid in the vicinity of a rigid wall where all interactions are defined by smooth atomic potential functions. These simulations are used to assess how the behavior of this system at the atomistic length scale compares to continuum mechanics models. In order to determine the effective size of the particles, we calculate the solvent forces on spherical particles of different radii as a function of different positions near and overlapping with the atomistically defined wall and compare them to continuum models. This procedure also then determines the effective position of the wall. Our analysis is based solely on forces that the particles sense, ensuring self-consistency of the method. The simulations were carried out using both Weeks–Chandler–Andersen and modified Lennard-Jones (LJ) potentials to identify the different contributions of simple repulsion and van der Waals attractive forces. Upon correction for behavior arising the discreteness of the atomic system, the underlying continuum physics analysis appeared to be correct down to much less than the particle radius. For both particle types, the effective radius was found to be $$\sim 0.75\sigma $$∼0.75σ, where $$\sigma $$σ defines the length scale of the force interaction (the LJ diameter). The effective “hydrodynamic” radii determined by this means are distinct from commonly assumed values of $$0.5\sigma $$0.5σ and $$1.0\sigma $$1.0σ, but agree with a value developed from the atomistic analysis of the viscosity of such systems.

are important in an increasing number of technological applications that involve flow in micrometer-and nanometer-scale channels [6][7][8].
Atomistic descriptions, such as molecular dynamics (MD) simulations, can characterize the fluid problems that are affected by these molecular aspects. While it remains difficult to use MD simulations to study macroscale hydrodynamics, confined regions, such as fluid-solid interfaces where most breakdowns of the continuum description occur, can be effectively characterized by these means. It is therefore important to understand details of how and at what length scales the breakdown of the continuum description occurs and how to best map the results of the atomistic simulations onto the continuum models. One aspect of this is defining where surfaces exist at an atomic scale, so that comparisons with continuum models can be accurately assessed.
One relevant physical situation that we are studying and where such detailed scrutiny can be accomplished is the behavior of a colloidal particle near a wall. For this case, analytic solutions of the Stokes equations exist and the problem is tractable for molecular dynamics simulations. In a seminal paper, Brenner obtained exact solutions of the Stokes equation in a viscous fluid for a spherical particle moving steadily toward and away from a planar surface of infinite extent [9]. For the case of a solid surface with the conventional continuum assumption of no-slip between the surface and the solvent, Brenner found that the drag force, F, on the particle is given by where μ is the solvent viscosity, b the radius of the sphere, U its velocity, c = 6 for no-slip boundary conditions, and λ is a correction factor which, for no-slip conditions, is given by with where h is the distance from the center of the sphere to the plane. At a large distance from the surface, h → ∞, λ tends to unity, so that Eq. (1) simplifies to Stokes' law for no-slip conditions (c = 6). At short separation distances from the plane surface (compared to the radius of the sphere), this continuum force diverges as 1/ h. At distances comparable to molecular dimensions, the continuum assumption should probably break down because of the discrete molecular nature of the fluid, and one would expect deviations from this ideal equation. Other assumptions implicit for Brenner's solution are a constant fluid density and infinitely smooth solid and colloid surfaces. In real atomistic systems, these assumptions are not completely valid at sufficiently small length scales, so discrepancies would be expected. Also, the full realization of no-slip conditions essentially requires that the solvent be irreversibly bonded to the colloidal particle, which blurs the definition of what should be part of the particle. Here, we note that Brenner's formulation has been extended to consider pure slip (c = 4) and partial slip (4 < c < 6) conditions, where this also requires modification in the definition of λ [10].
There have been a several previous studies employing molecular dynamics (MD) on the validity of Stokes' law in the case of a spherical particle moving toward a wall, with mixed opinions on the interpretation of the results. Here we note that in most MD studies, the somewhat indefinite sizes of atoms are represented by describing the forces between atoms using analytic functions that decay to zero over some range. Vergeles et al. studied the translational and rotational motion of a particle in a viscous Lennard-Jones fluid both near and away from a flat surface [11,12]. Their study, which included both simple spherical and atomistically structured particles, showed that the drag and torque on the sphere in an essentially unbounded fluid are in good (though not quantitative) agreement with the continuum hydrodynamics results captured by Stokes' law. Similar studies were also performed by Heyes et al. [13,14]. As the solute approaches the wall, Vergeles et al. concluded that the Brenner model was supported by their data, although some variations were observed at separation distances comparable to the solvent radius. These authors also concluded that the no-slip condition breaks down near the wall for their systems.
In a more recent study that applied repulsive Weeks-Chandler-Andersen (WCA) potentials, Challa and van Swol revisited parts of the Vergeles studies [15]. In their analysis, they found that as a sphere approaches a flat surface, the net drag force on the sphere can be represented by a superposition of two contributions: a static solvation force and the hydrodynamic force. They found that the static solvation force contributes most prominently at small distances between the sphere and the planar surface and at small velocities. This force oscillates between positive and negative values and reflects the discrete size of the solvent molecules, an effect that is not included in the continuum hydrodynamic result [1,[6][7][8]. After compensating for this, they found that the divergence in the force on approaching the wall was not nearly as strong as that predicted by the Brenner analysis.
There have also been several related MD studies assessing the validity of the Stokes law for particle drag forces, where the results demonstrate qualitative validity, but there is disagreement on the inferred value for c in Eq. (1) and how it should be interpreted [5,[16][17][18]. In order to make any of the above comparisons to the continuum models, it is required to define values for the radius b of the colloidal particle (or tagged solvent) particle and for the distance h between the particle and the wall. While critical for the quantitative evaluation of these relationships, these values have been treated in an ad hoc manner. For example, if the range of the Lennard-Jones (LJ) interaction is σ , one possible choice for the radius of the particle would be σ/2 for both the LJ and WCA models. However, the use of such values implies values for c that appeared to be unphysically high. Consequently, concepts of a larger effective "hydrodynamic" radius were developed that considered the volume from which solvent particles (or more specifically the solvent nuclear position) are excluded. In this case, the radius of the aforementioned LJ particle would be increased by about an additional σ/2 to a net value of ∼ σ . Such values tend to be more consistent with the expected values for c, though the physical basis for this argument is somewhat weak. For defining the position of the wall, both of these options have been used as well as the most simplistic assumption that the wall position coincides with the surface atom centers. These differences are not trivial adjustments for these analyses and may be contributing to the differing conclusions. They certainly are required if one wishes to do truly quantitative analysis of the simulations.
In this paper, we focus on investigating the wall-particle interactions, including solvation forces, for a single suspended, spherical particle in a viscous fluid. Ultimately our analysis should be applied to realistic nanoparticles composed of multiple atoms. However, this is not a good starting point, as it would require a vast parameter space of defining various clusters of atoms in near-spherical shapes, defining the interaction potentials and spacings between those atoms. These parameters would also affect the drag coefficients significantly. Moreover, all this work must hinge upon first defining the effective size of the atoms on a physical basis.
It is also important to remember that all high-quality atomic potentials have analytic representations with diffuse tails due to approximately exponential decay of electron density around atoms. This aspect makes defining the size of the atoms difficult, leading to finite forces of interaction between even neutral atoms at rather large atomic separations (i.e., well beyond what one might define as "contact"). The literature that we cite above illustrates different ad hoc assumptions for the effective radius of an atom that is described by the Lennard-Jones-type potentials. Typically, the atomic radius is taken to be either 0.5 or 1.0 times the σ parameter. Based on those assumptions, researchers obtain very different (and non-physical) values for drag coefficients and how a particle in a fluid interacts with the wall compared to Brenner's model. Consequently, in the following narrative we reduce the problem to a spherically shaped particle for which there is an analytic continuum solution that can be directly compared to. It is by no means the finish point of the work, but a necessary and unavoidable starting point.
Within the scope defined by the previous paragraph, our interest is to elucidate the effective sizes of the particles and the position of the wall by using several different metrics. In a subsequent paper, we will build on these results to carry out the analysis of the drag force on a colloid particle and to assess the validity of Brenner's expression as the suspended particle approaches the wall. The work here will focus on the analysis of equilibrium calculations with a static colloid particle. By separating the force contributions between the three species (colloid particle, fluid, and wall), and applying continuum analogies to the behavior of those resulting forces with respect to the radius and separation distances, the effective dimensions of the particles can be derived. This will rely on configurations that could be considered unphysical (e.g., strong overlap between the colloid and the wall), but these are required to demonstrate where particles effectively begin and end interactions with the other species. The results of this study are also relevant to the general study of fluid surface interactions and establish means to define the effective position of the wall in those studies [5,17,19]. Our emphasis is on analyzing the forces directly observed by the colloid particle and other particles of the system. This is a crucial aspect for developing a self-consistent model for the analysis.

Simulation approach
All simulations reported here were performed with the LAMMPS package for molecular dynamics [20]. The fluid particles have a mass m and interact via a pairwise potential. For the fluid-fluid interaction potentials, we used two different variants of the Lennard-Jones (LJ) potential. For the first set of simulations, we used the well-known repulsive WCA potential [21] that was used in many previous studies as shown in Eqs. 4 and 5.
The cutoff radius r c for this potential is chosen at the minimum of the LJ potential, r c = 2 1/6 σ ≈ 1.12σ , and r i j is the distance between the centers of the two particles, i and j. The parameter ε is the depth of the LJ potential and represents our fundamental energy unit, and similarly, σ is the fundamental length unit. The potential between the fluid particles and the suspended sphere is described using a shifted WCA potential, as was used in previous studies, shown in Eq. 6.
Here, b 0 is the radius of the hard-core region of the sphere (which allows its effective size to be arbitrarily adjusted) and r i j is the distance between the center of the sphere and the center of a fluid particle. These potentials are effectively continuous and smooth to arbitrary order and short-ranged. (There is a numerical divergence at the inner cutoff b 0 but because of the rapidly increasing repulsive force approaching that cutoff, this region is not accessible to the other particles under reasonable setup geometries and timestep controls.) Consequently, there are no complicating factors in their numerical evaluation and, for the standard time step of t = 0.005σ √ m/ , there is a very high degree of energy conservation. We used the same parameter values for WCA (and for modified LJ) as Challa and van Swol [15].
To explore the effect of attractive van der Waals forces under similar simulation conditions, we also utilized a modified version of the LJ potential. Because of the moderately long range of the LJ potential, many previous studies either truncate the interaction at an arbitrary distance (resulting in an expression that is neither continuous nor smooth at the cutoff distance), or employed a cut-and-shifted potential form that is continuous at the cutoff distance, but is not smooth there. Consequently, the simulations do not conserve energy during the dynamics runs, and the method of simulation and/or the analyses require corrections for this artifact. Other approaches to this problem are either to include a very large range of pairwise interactions (which is very computationally intensive) or to more smoothly truncate the tail of the LJ potential [22]. Here, we define a modified LJ form that smoothly (to second order) approaches zero at a distance r max . This potential consists of the standard LJ expression up to its inflection point at r spl = (26/7) 1/6 σ ≈ 1.24σ , as proposed earlier by Holian et al. [22]. There (and distinct from [22]), it switches to a quartic polynomial spline that is smoothly matched (to second order) at this point and goes smoothly to zero at a non-arbitrary value r max . Thus, this potential has the form shown in Eq. 7 below, where φ 12−6 r i j is given in Eq. 5.
The constants A and B are chosen so that the function φ spl and its two first derivatives are continuous at r spl and r max . The net result is a LJ-type potential with a tractable finite range (< 2σ ) that has a smooth second-order derivative and a continuous third-order derivative. Consequently, the simulations naturally conserve energy to a high level for the timestep size chosen above and require neither thermostatting during the simulations nor corrections for their analysis. For efficient calculations, we have implemented this as a tabular formulation using standard LAMMPS methods. We have verified a very high level of energy conservation with this potential and its application to long timescale simulations [23]. The fluid-sphere interactions are given by similar expressions where the sphere core radius b 0 is analogously introduced as in Eq. 6. The simulations were performed in a box of dimensions L x × L y × L z = 13.68σ × 13.68σ × 32.20σ with periodic boundary conditions applied to the x and y directions, while a pair of parallel walls was placed normal to the z direction near the two ends of the box. The atomistic wall that was used to explore the behavior of the suspended spherical particle consisted of two layers each of 200 atoms. These were placed in a facecentered cubic (fcc) arrangement with a lattice spacing of 1.094σ to form a wall that would be impervious to the solvent atoms. The inner most layer of wall atoms was positioned with their atomic centers at z = 0. These wall atoms interacted with the solvent atoms with the same potential as the solvent-solvent interaction, but their positions were fixed during the simulations because they were effectively given very large masses (10 4 m). Collisions with them were perfectly elastic and conserved energy, though their inclusion corrupts the manner by which the pressure is calculated within LAMMPS. The wall at the far end of the simulation box (z = 32.20σ ) was a simple smooth wall that interacted with the solvent atoms by a planar 9-3 potential as used in earlier work [15]. This wall is well removed from the colloidal particle and the atomistic wall so that it did not perturb any of the interactions there, as the distance was several times greater than the cutoff distances of the interaction potentials. As discussed below, this wall was included as means of directly measuring the pressure of the system. The potential energy between the suspended sphere and the wall particles was set to zero so that the wall particles exerted no direct force on the sphere. In all simulations described here, the position of the colloid with respect to the wall was fixed during the data acquisition periods, so the magnitude of the wall-colloid interaction has no significance on the simulation.
The approach described above does not take the thermal motion of the wall atoms into consideration. This is a deliberate assumption that we have made to limit the parameter space of our study. Let us elaborate on the physical reasons justifying this assumption. Generally, the walls are denser and the atoms are much more tightly bound to one another than the fluid-fluid or fluid-wall interactions. Consequently, the range of thermal motion of the atoms in the wall will be much smaller than those in the fluid, at least a factor of 10 for say a polymer plastic-based walls interacting with an organic fluid, and about a factor of 100 for a denser metal wall and an organic fluid. Again, there is a large range of choices to be made here, but in all cases, this should be a minor contribution.
The simulations were initialized by filling up the region between the confining walls with fluid particles. The fluid particles were initially placed in an fcc crystal structure with a lattice constant 1.094σ giving an initial solvent density of ρσ 3 = 0.8. The first layer of solvent was placed at z = 1.094σ . Then, a colloid sphere with a core radius of b 0 and mass M was placed at a specified point defined below. Any solvent particle whose center was within b 0 + σ of the colloid center was then removed from the simulation. This method generated fairly consistent solvent densities regardless of the size and initial position of the colloid. The system was then randomized for 40,000 time steps using a Nosé-Hoover thermostat at a high temperature, kT /ε = 3, to liquefy and randomize the solvent. The temperature was then lowered to a value of either kT /ε = 1.2 or kT /ε = 1.6, again using a Nosé-Hoover thermostat, and the system was allowed to equilibrate for 20,000 time steps. Subsequently, the simulation was continued in an NVE ensemble, and the z-component of the net solvent force on the sphere was evaluated a total of 1000 times, collecting the atomic forces every 200 time steps so that each data collection point should be uncorrelated with respect to the others.
This procedure was performed for 100 different static positions of the suspended sphere for each selected core radius, potential set, and temperature. The initial starting position was chosen to be 8σ from the position of the atomistic wall, which is then 20σ from the back wall. For the largest sphere, this leaves ∼ 3σ between the outside of the sphere and the surface of the wall (see results below). This is the amount of separation where perturbations from the discrete nature of the solvent have been observed to be suppressed from previous simulations [11][12][13][14][15][16][17]. The range of sampling positions was then incremented down from this starting location where the center of the sphere was well separated from the wall, and was repeated until the sphere had moved completely through the wall and there was no longer any perturbation to the solvent. Effectively, this tracks the equilibrium forces resulting from moving the colloid from the bulk fluid and through the wall. By reinitializing The pressure of the simulation system was calculated by integrating the total force from the solvent on the smooth wall at the far end of the simulation box. This was averaged over 100 snapshots that were collected every 200 timesteps. Since the area of the back wall is precisely known (13.68σ × 13.68σ ), the pressure is simply calculated by dividing the net average force by the area. Figure 1 illustrates the coordinate system used for the analysis and the expected results based on a continuum model. The position of the wall is assumed to be not known exactly and is defined to be at z = + 1 relative to the origin, where the origin is the position of the atomic centers that define the surface layer of the wall. The radius of the colloid is also not known exactly and is defined as b = b 0 + 2 . Here, b 0 can be associated with the core radius defined in the potential energy of Eq. 6, and 2 defines the additional radial extent of the effective colloid radius that will be determined from the derived value for b. The z position of the center of the colloidal particle is defined as being at r with respect to the origin. The (x, y) coordinates are defined as (0, 0) where this defines the centerline of the simulation box. The values of 1 and 2 will be determined by characterizing the net force on the colloid from the solvent as the colloid passes through the wall. The force between the wall and colloid will be neglected in this analysis, though all other interactions are included in the simulation.

Analysis method
In the continuum limit, the solvent force on the colloid comes from the hydrostatic pressure, p, that is exerted by the fluid normal to the surface of the suspended particle. In this coordinate system, the net force in the z direction is given by the difference of the force on the right hemisphere minus the force on the left hemisphere. In this continuum approximation, the net forces on the sphere will be zero for all positions Fig. 1a shows the case of r = b + 1 , which is where the sphere is just touching the wall surface. The net forces in the perpendicular directions (x, y) will be zero under all conditions because of the cylindrical symmetry of the system about the z axis.
As the sphere starts penetrating the wall (Fig. 1b), the intersection plane between the wall surface and the sphere can be defined by the angle θ a , defined with respect to the equatorial plane of the sphere. Our convention is to define θ as positive for the left half of the sphere, and it will be negative if the intersection is on the right half of the sphere. In this scenario, the net force along the z direction on the hemisphere on the right side can be calculated by integrating the z component of the force exerted by the fluid on the right hemisphere, f z = psinθ , over an annular element of area, 2πb 2 cosθ dθ , where the integration is over the half surface − π/2 ≤ θ ≤ 0: We note that f z as defined above is a negative quantity (pointing to the left) since the sign convention is that θ and sin θ are negative for the right side of the sphere, so the net force is negative as observed. The corresponding force on the left side of the sphere is: Here, as shown in Fig. 1b, r a is the distance from the surface of the wall to the center of the sphere, and r b is the radius of the disk defining the intersection of the wall and sphere. This force is positive and points away from the wall. At r = 1 , when the center of the sphere is at the wall surface (Fig. 1c), F left is seen to go to zero, and there is a maximum force on the sphere of − π pb 2 . As the sphere continues passing through the wall (Fig. 1d), F left remains zero, and F right diminishes as the integration limits are decreased to − π/2 ≤ θ ≤ θ a .(θ a is now negative.) That integral then equals −π p(b 2 − r 2 a ). More succinctly, one can write a general equation for the net force by summing these two terms: This then needs be recast into the reference frame for use in the analysis. Since r a defines the position of the center of the sphere relative to the wall surface, the absolute z coordinate for the center of the sphere is r = r a + 1 . As mentioned above, the colloid radius will also be more specifically defined as b = b 0 + 2 . Substituting those forms into Eq. 10 generates the equation form that will be used as a basis for the analysis below: In the ideal case, the results can be fit to a general parabolic expression (e.g., cr 2 +dr + e), and the values of 1 and 2 directly extracted from that fit, as well as a value for p, by using Eq. 11. A simple direct method of evaluating these quantities follows from the derivation above. The most extreme value for the force, F ex , should occur at r = 1 , and its value is given in Eq. 12.
The two roots of the quadratic, r 1 and r 2 , are where the net force becomes zero and these are given in Eq. 13.
Formally, the value of p can be defined from the second derivative of Eq. 11 as shown in Eq. 14.
We note that there are more accurate means for evaluating p from the simulations (e.g., the net force on the back wall), and the value of p derived from our fitting procedure is used primarily as a check of the efficacy of the model. Given that p is available from some alternative source, that value can be used with Eq. 12 along with the value of F ex to derive 2 by a second means.
As shall be seen below, some of the data has a noticeable skewness, and this requires the addition of a cubic term for analysis. One form for adding this is defined in Eq. 15 below.
The assumption for using this form is that the cubic term should be a small correction to the quadratic fit such that κ (r − 1 ) 1 over the fitting range of the data. This particular form was chosen for the analysis so that the definition (and interpretation) of the two roots (intercepts) remains as defined in Eq. 13, and that one can simply write F( 1 ) = − π p (b 0 + 2 ) 2 . Physically, the addition of this multiplicative factor that linearly depends on r scales the original quadratic force as a linear function of r . In the absence of a particular mechanism, there are no means to define where this term has no effect, so the formulation now has an effective pressure p that cannot be related to the true pressure without further information. However, this correction to the pressure should be small. Fitting this cubic equation to the data will then generate unique values for 1 , 2 , p and κ, though the absolute interpretation becomes less clear.
A particular caveat on this fitting form is that the most extreme value of the force, F ex , no longer occurs when the center of the particle is at the wall surface, r = 1 . Instead, the position of the extremum, r ex , has now been shifted slightly and is given in Eq. 16 to first order in κ.
Comparing the location of the minimum relative to the positions of the two roots is then a simple graphical means of evaluating κ. Finally, we note that evaluating the second derivative of the function at the minimum generates an analogous relationship to Eq. 14 shown below.
At r = 1 , the second derivative is then simply 2π p (similar to Eq. 14), while at r = r ex , there will be some correction terms that are second order in κ.
The skewness correction can be interpreted as the consequence of the effective deformation of the original spherical particle. The simplest physical model for such a deformation is an ovaloid particle. Parallel to the wall surface, the particle retains its original radius (r = b 0 + 2 ), but normal to the wall surface the particle now has one long semi-axis (b l = b 0 + 2 + 3 ) and one short semi-axis (b s = b 0 + 2 − 3 ). The asymmetry of those axes can be derived from the cubic correction terms as shown in Eq. 18, where this supplies a physical metric for interpreting the skewness.
Working through that particular model, one could show that the extremum in the force goes back to being the physically interpretable F ex = − π p(b 0 + 2 ) 2 (note p rather than p ), and that the position of that extremum, r ex = 1 + 3 , occurs when the maximum particle cross section (which is no longer at the center of the particle) coincides with the wall surface.

Results
Simulations were run for five different sphere sizes with radii b 0 /σ = 0, 1, 2, 3, and 4. The case where b 0 /σ = 0 corresponds to the "suspended" particle being exactly equivalent to a fluid particle. The remaining cases correspond to an increasingly larger sphere, though the range and energetics of the interactions between the sphere and solvent were maintained fixed. Since the simulation box size is fixed, the number of fluid particles in the box decreases as the sphere size increases by the algorithm we employed in the setup. Similarly, as the sphere passed into the wall, the number of fluid particles was increased to compensate for the sphere's partial presence in the fluid region. Overall, the setup algorithm we employed kept the solvent pressure maintained at a fixed value ± 5% for the radii of b 0 /σ ≤ 3, with only a slight correlation to the sphere's position. For the largest spheres, b 0 /σ = 4, larger deviations were observed (approaching 10%) which were correlated with the sphere going in and out of the wall. Our analyses below are based on assuming a constant pressure for the system. For the b 0 /σ = 4 system we performed an additional analysis that explicitly included the pressure dependence. This changed our estimate of the derived parameters by only a few percent, which is smaller than their error bounds.

Repulsive Weeks-Chandler-Andersen potential
The first set of simulations was carried out using the repulsive WCA potential defined in Eqs. 4 and 6. The z component of the solvation force for b 0 /σ = 0, 1, 2, 3, and 4 at temperatures kT /ε = 1.6 and kT /ε = 1.2 as a function of the position of the particle center was determined. The results for the kT /ε = 1.6 simulations are shown in Fig. 2. The results for kT/ε = 1.2 are very similar except that the values for the forces are reduced by ∼ 20%. The average pressures for the kT /ε = 1.6 simulations are ∼ 8.2, while those for the kT /ε = 1.2 runs were only ∼ 6.7. This difference in pressure is reflected then in the magnitudes of the forces.
This figure reports the solvation force, F (ε/σ ), on the sphere as a function of its position with respect to the wall as outlined in Fig. 1. At left, the particle is completely shielded from the solvent by the wall and experiences no force from the solvent. As the particle emerges through the wall, it experiences a negative force  Table 1. The open symbols are data points not used in the fits, and the gray lines are added as a guide to the eye, especially to highlight the solvent "ringing" at the right Table 1 Summary of fitting parameters for different particle types. The definitions of the parameters are given in the footnotes and described in the analysis section from the solvent that would push it back into the wall. The magnitude of this force increases on a parabolic trajectory until it reaches a maximum size when the sphere should be halfway through the wall. The magnitude then decreases parabolically and would be expected to return to zero as the particle emerges from the wall at right according to the continuum models. Here, however, we observe that the force weakly oscillates between positive (repulsive) and negative (attractive) values with a period that is comparable with the size of the fluid particles. This phenomenon has been noted previously and is associated with discrete layering of the solvent between the colloid and the surface [15].
With the exception of the solvent ringing, the data appear to be very well described by the parabolic form predicted with the continuum model. This is verified by a least squares fitting of the data points with negative forces. The curves resulting from those fits are shown in the figure, showing a very good quality of fit. Overall, we note three major features that are consistent with this model: (1) the maximum force increases approximately as the square of the radius of the particle, (2) the curvature of the parabolas are all close to the same value, and (3) the minima of the parabolas occur at nearly the same location. The parameters derived from these parabolic fits are given in Table 1 and will be discussed in greater detail below.

LJ-spline potential
We used the same simulation approach as for the WCA potential to run these simulations. The values of the sphere radius are again b 0 = 0, 1, 2, 3, and 4, and the temperatures are kT /ε = 1.2 and kT /ε = 1.6. The pressures for these two sets of runs are ∼ 3.2 and ∼ 4.9, respectively. The pressures here are significantly lower than for the WCA potentials at the same density and temperatures because of the addition of the attractive LJ forces. The results for kT /ε = 1.6 are plotted in Fig. 3, where one can see that they are qualitatively similar  Fig. 2, the curves are offset vertically by 400, 300, 200, 100, and 0 units, respectively. The closed symbols were used for the parabolic fits, which are shown as the black lines, and the parameters for those fits are listed in Table 1. The open symbols are data points not used in the fits, and the gray lines are added as a guide to the eye, especially to highlight the solvent "ringing" at the right. Compared to Fig. 2, there is now a noticeable attractive bump on the left-hand side, the solvent "ringing" is somewhat stronger, and the quality of the parabolic fit is slightly poorer to those found for the WCA fluid. The most striking difference is that the forces are significantly smaller than for the WCA potential, where these simply scale with the system pressure. The solvent layer oscillations on the right-hand side are stronger and more persistent, which reflects the higher level of ordering induced by the attractive LJ forces for a fixed temperature. There is also now a noticeable positive force on the left-hand side of the figures when the particle first emerges through the wall. This is because the attractive forces are longer-ranged than the repulsive forces (r −6 vs r −12 ), so an initial attraction of the particle to the solvent is sensed before the more dominating repulsive forces become apparent. The data here appear to be reasonably well approximated by a parabola, though some skewing is apparent, especially for the smaller particles. As before: (1) the maximum force increases approximately as the square of the radius of the particle, (2) the curvatures of the parabolas are close to the same, and (3) the minima occur at nearly the same location. The parameters derived from these parabolic fits are given in Table 1 and will be discussed in greater detail below.

Analysis
The results of the various ways of fitting the data are summarized in Table 1. The quadratic model overall fit the data quite well with r 2 > 0.95 (and typically r 2 > 0.98), except for the LJ case with b 0 = 0 where r 2 = 0.84. The results were generally weakest for the small particles where we collected the fewest number of data points. The addition of a cubic term increased the quality of the fit generally as would be expected, but it is most significant for the LJ particles. The quality of the fits was now all r 2 > 0.98 except for the LJ case with b 0 = 0 where r 2 = 0.94. A comparison of the quadratic vs. cubic fits is illustrated in Fig. 4 for the b 0 = 1 case for the WCA and LJ particles. A very slight skewing to the left is found for the WCA particles (negative value for κ), while a more significant skewing to the right is found for the LJ particles (positive value for κ). We now discuss each of the parameters in turn.
First we consider the value of 2 , which measures the effective radius, b, of the suspended particle, where b = b 0 + 2 . A plot of select 2 values versus the values for b 0 is shown in Fig. 5. The major conclusion here is that the values from the quadratic and cubic polynomial fits, 2 (q) and 2 (c), respectively, are all There is a very slight skewing to the left for the WCA cubic fit and a more marked skewing to the right for the LJ cubic fit remarkably consistent with 2 = 0.75 ± 0.10σ . The data do suggest a minor subdivision in that the results for the WCA particles tend to decrease slightly as a function of b 0 , while those for the LJ particles tend to increase slightly as a function of b 0 . However, there is an obvious lack of dependence on pressure over the range of p = 3-8, and a similar lack of dependence on T over the range kT /ε = 1.2-1.6. There are also no significant differences in the values determined in the quadratic fit versus the cubic fit. At one level, this is not at all surprising since the separation between the two roots is one means to estimate 2 (see discussion associated with Eq. 13), and that attribute does not depend strongly upon the choice of the fitting function. However, it is still noteworthy that we observe remarkable consistency in these values.
We emphasize that this first method of determining these values relies solely on the shape of the forcedisplacement curve and does not rely upon the measured values of the pressure. An alternative approach is to compare the maximum observed force with the bulk pressure, which has been measured independently from the force on the back wall (see Eq. 12). These values determined for the maximum force from the quadratic and cubic fits, 2 (fq) and 2 (fc), respectively, are given in Table 1. These also show good self-consistency, and these are also shown as the gray symbols in Fig. 5 for kT /ε = 1.6.
Next, we consider the values for 1 , which measure the effective position of the wall. A plot of select 1 values versus the values for b 0 is shown in Fig. 6. For the WCA particles, we again observe a remarkably consistent value of 1 = 0.65 ± 0.05σ from the quadratic and cubic polynomial fits, which is independent of pressure/temperature, particle size, and the order of the polynomial fit. For the LJ particles however, there does seem to be a dependence on the particle radius. This is particularly noticeable for kT /ε = 1.2, but less so for kT /ε = 1.6. We note that the values derived from the cubic fits show less of a dependence on the radius, suggesting that the results from the quadratic fits are being perturbed by the source of the anharmonic interaction. For the smallest LJ particles, the values of 1 are consistent with the results from the WCA. We note here that one expects the effective position of the wall to depend on the spacing pattern of the wall atoms (which defines an inherent corrugation or roughness), and the relative size of the wall, solvent, and particle, so these results are not surprising. The attractive forces of the LJ interaction will draw the solvent and particles into closer contact with the wall, effectively moving the wall surface back with respect to the position determined for the WCA interaction.  An interesting outcome of our analysis is that the pressure of the system could be evaluated directly from the curvature of the force-displacement curve (Eqs. 14 and 17). The values determined from the curvature of the quadratic and cubic curves, p(q) and p(c), respectively, and the bulk value for the pressure, p, determined from the net average force on the back wall are listed in Table 1. A comparison of select values as a function of b 0 is shown in Fig. 7. There are obvious discrepancies for the smallest particles, but those data sets also have the fewest number of sample points and the largest uncertainties in the derived values for p. A likely origin of Values of kT /ε are shown in the legend; "c" refers to cubic fit, and "w" to wall value this comes from the fact that the solvent is always allowed to sense forces from both the wall and the particle, especially where the wall and particle strongly overlap.
Why allow such overlap at all? There are two reasons. First, if the particle is not forced into overlap with the wall, one cannot be certain that contact has been made and one cannot then define precisely where the wall is. Again, in previous work, others have assumed what the point of contact is based on their ad hoc assumptions for the atomic radii, leading to unsatisfactory comparisons with the Brenner model and disappointing results for the derived drag coefficients. We need to define the position of the wall from the rapidly increasing repulsive forces that arise during the forced overlap, and the changes in the solvation forces. Second, by allowing the colloidal particle to continue to completely penetrate the wall, the dependence of the solvent forces on the particle position could then be used to define the size of the particle and the position of the wall self-consistently and under the conditions of the thermal simulations we will subsequently use for the dynamic analysis. While the overlap simulation is non-physical, it is mathematically well justified to determine the size of the particle: we are effectively constraining what fraction of the surface of the particle is exposed to the thermalized solvent and using that to determine the size of the particle according to the solvent response.
In a more physically realistic simulation, the conflicting wall particles would be removed (or their interaction forces with the solvent removed), but this then becomes a more complicated and ill-defined procedure in defining which those should be removed for a particular geometry. The net result is that the solvent experiences an extra repulsive force from the WCA wall particles and an extra attractive force from the LJ wall particles and this somewhat biases the average solvent positions. This then diminishes/enhances the net force that the particle experiences from the solvent, so that this pressure measurement is not directly valid. This is obviously a major correction for the b 0 = 0 particles, but significantly diminishes for increasing particle sizes. Overall, the general level of consistency in these values supports the general efficacy of this means of analyzing the system.
For completeness, we also consider the anharmonic terms that are subtly apparent in our analysis. The values for the skewness, κ, versus b 0 are shown in Fig. 8. For the WCA particles, the skewness is of only minor significance, having a slightly negative value for the smallest particles, b 0 = 0, and is otherwise negligible. Consequently, one expects the continuum model to work quite well for those systems generally. For the LJ particles, the skewness is more apparent and positive, with it being greatest for the smallest particles. As discussed above, one method for interpreting the skewness is to define it in terms of how different the effective radii of the particles are on the right-hand side vs. the left-hand side. That is, one can interpret this as a spherical particle becoming deformed into an ovaloid. These values are defined as 3 and are presented in Table 1. The  Fig. 8 Comparisons of values for κ evaluated by the cubic fits for both the WCA and LJ particles of different radii for kT /ε = 1.2 and 1.6. There is a slight negative correction for the WCA potential and a stronger positive correction for the LJ potential calculations of these values are subject to large uncertainties, but they seem to be consistently 0.2-0.3σ , with the right-hand radius (projecting toward the solvent) being shorter than the left-hand radius (projecting toward the wall). This can be attributed to the attractive portion of the LJ potential pulling the colloidal particle into the solvent, shortening the apparent core radius on that side. Then, as the particle emerges from the wall, the mutual attractive forces between the solvent, wall, and particle enhance the effective attraction between the particle and wall, appearing as a net repulsive solvent force on the particle, effectively lengthening the apparent core radius on that side. The scale of these effective radii changes measures the difference in range between the attractive and repulsive parts of the potential, consistent with ∼ 0.2σ .
The other method of analysis suggested by our approach is to consider the dependence of the maximum force on the radius of the particle (Eq. 12). It is most straightforward to linearize this and plot the G = √ −F ex versus b 0 as shown in Fig. 9. Excellent linear correlations are observed for the four sets of data. In the most direct interpretation, the intercepts for G = 0 then define the values for 2 . For this figure, those results are 2 = 0.75σ for the WCA systems and 2 = 0.90σ for the LJ systems. The results are markedly consistent with our other analyses. We note that this interpretation is only rigorously correct if the pressures are constant for the five component points, which is not completely accurate. Our simple algorithm for determining how many solvent molecules should be placed in the box for each particular position of the colloidal particle did not maintain a perfectly constant value of p or solvent density. We can correct for this however, by considering the correlation between H = √ − F ex /π p versus b 0 . These should have slopes of 1 and intercepts of 2 . Those plots are shown in the inset of Fig. 9, and the observed values of the slopes are in the range of 1.00 ± 0.02 and values of 2 in the range 0.75 ± 0.10σ for that correction.

Discussion
The force-displacement curves were all very well fit by parabolic curves predicted by the continuum analysis. This was most surprising for the smallest colloids, especially where the spectator particle was the same size as the solvent. Though anharmonic/cubic terms were apparent, the adjustments for them are relatively minor and easily rationalized in terms of the physics of the problem. The simple continuum analysis model should then be a highly appropriate and consistent manner by which to treat the data.
The values for 2 were determined from the parabolic ( 2 (q)) or cubic ( 2 (c)) fits of the force versus position curves (data columns 5 and 6 in Table 1). These fits consistently produced values within a 0.75±0.10σ range. This range is essentially our measurement error due to fluctuations in the simulations. These values can be essentially determined by the roots of the fit: the points where the solvent force on the particle goes to zero. Our initial expectations were that this value would be 2 1/6 /2σ (∼ 0.56σ ) for the WCA potential, and either 0.56 or 0.5σ for the LJ potential. This is based on the logic that the size of the particle should be defined by where the force of interaction with other particles becomes zero. However, we are clearly seeing something significantly larger than this. One possible flaw in this analysis, especially for the smallest particles, would arise from the discrete size of the solvent molecules. Considering Fig. 1a, when the solvent molecule is of Fig. 9 Dependence of the solvent force extremum, F ex (determined from cubic fits), on b 0 . The y axis is G = √ − F ex and the data labels are: squares: WCA particles at kT /ε = 1.6; circles: WCA particles at kT /ε = 1.2; diamonds: LJ particles at kT /ε = 1.6; triangles: LJ particles at kT /ε = 1.2. The G = 0 intercept provides values for 2 , provided that p remains constant for each series. The inset shows corrections for the non-constant p, where the y-axis is now H = √ − F ex /π p. The data points are now nearly indistinguishable with slopes ≈ 1.0 and H = 0 intercepts consistent with 2 ≈ 0.75 comparable size to the probe particle, one would expect that the probe particle must lift away from the surface by some amount in order that the solvent could apply some force in the outward direction to counterbalance the solvent force in the inward direction. Therefore, we considered an alternative measure for the radius of the particle: the ratio between the maximum solvent force observed and the independently evaluated pressure. Because of the lack of steric hindrance around the half-embedded colloid, we expected this should be the cleanest measure of the particle radius. However, this method also consistently produced values of 2 ∼ 0.75σ , for both the quadratic ( 2 (fc)) and cubic ( 2 (fc)) fits. We consider it especially impressive that very consistent results were obtained for the smallest particle by both methods since that particle then has the exact same potential as the solvent particles (e.g., we are evaluating the solvent properties). We were expecting substantial deviations here because we are essentially using solvent molecules as the probes, and did not expect to be able to evaluate dimensional properties that would be smaller than the probe. However, the statistical sampling we achieve by time-averaging over many distinct geometries enables this finer scale measure, so that we achieve an effective resolution several times smaller than the solvent diameter.
Previous workers had also noted that the expected radius of 0.5 or 0.56σ for b 0 = 0 type particles did not yield results consistent with their expectations from hydrodynamic analyses [4,5,[12][13][14]16,17]. This is often based on analysis of Eq. 1 where the values of the viscosity μ and the slip factor c are presumed to be well defined and known. To achieve consistency with those expectations, those workers used an effective "hydrodynamic" radius, a concept apparently first introduced by Bocquet [5]. Various definitions of this quantity have been proposed, though not rigorously justified. Typical values for b 0 = 0 type particle have been 1.0σ , a value which is too large by our analysis.
A more justifiable evaluation of an "effective" particle radius was developed by Gaskell et al. in their study of the viscosity of simple fluids [24,25]. For their development, they required a mapping between the discrete velocity of their microscopic atomic simulations and a continuous velocity field for the macroscopic continuum analog. To achieve this, they define an effective form factor for the solvent atoms based on the notion that all of the volume in the simulation box must be assignable as belonging to some atom, or unitarily partitioned between nearby atoms, if one is going to make a connection to continuum analysis. All other physical quantities (e.g., velocity) can then be similarly mapped between discrete microscopic atomic quantities and the macroscopic continuum fields using these form factors. A similar approach is used in the analysis of general thermomechanical quantities from atomic scale simulations [26][27][28]. An effective size radius, a, can then be derived from this form factor. In its simplest instantiation, the form factor is approximated as a step function with radius a. The effective volume of the solvent atom, V a , is simply V /N where V is the volume of the simulation box and N is the number of solvent atoms in that simulation box [24]. The radius is then simply determined from (4/3)πa 3 = V a . Balucani et al. also explored the use of smooth form factors, specifically f (r) = exp[− (r/a) n ], but found that these increased the value of a by only a few percent.
This analysis can be directly applied to our systems for b 0 = 0. While we know the precise x, y dimensions of our system from the periodic boundary conditions (13.68σ × 13.68σ ), the z dimension is not as precisely known. From our system definition, the left-hand boundary should be z left > 0, and the right-hand boundary z right < 32.20σ . For estimation purposes (which are justified by our results below), we assume that the actual positions of each of these walls are shifted ∼ 0.5σ to the interior, so that the z dimension of the box is ∼ 31.2σ . The effective spherical radius can be then defined from this volume, where the result is not very sensitive to our uncertainty in z (± 3%). For our simulations, with ρ = 0.8 relative to a contact fcc structure, the b 0 = 0 systems contain 4238 solvent atoms. This defines an effective radius of ∼ 0.70σ , which is clearly consistent with the results of our analysis above.
For the systems with b 0 > 0, we cannot directly apply this logic since we do not a priori know the relative size of the colloidal and solvent particles. However, it seems valid to assume that the solvent molecules have not changed effective size for our different simulations, where we have tried to maintain a constant solvent density, as well as pressure and temperature, for the related simulations. We can then subtract the appropriate solvent volume determined by the number of solvent atoms, from the simulation box size to determine the volume for the suspended particle. The radii determined by this means are all consistent with our determinations above, a ≈ b 0 + 0.75σ . We note that our imperfect algorithm for removing solvent molecules to compensate for the size of the colloidal particle and maintain a constant solvent density could well be the cause of slight drifts we observed for 2 in Fig. 5.
The wall positions defined by 1 require a bit of additional discussion. The overall values, especially for the WCA potential, are comparable to, but always slightly less than, the values for 2 . The primary reason for this is that the walls were deliberately set up as a loosely packed fcc lattice (ρ = 0.8 relative to the WCA spheres being in contact), so that they would have a significant amount of roughness for our subsequent fluid flow simulations. Although the walls remain impervious to the solvent, the surface is crenulated at the atomic level, and the positions of the outermost layer of atoms (from which 1 is defined) are more representative of the extreme positions of a crenulated wall rather than its average value. The values of 1 = 0.65σ being slightly less for 2 for the WCA potential are rationalized by this means.
The 1 values for the LJ potentials were always found to be smaller than those for the WCA potential and show a noticeable dependence on the particle radii. The attractive part of the LJ potential will pull the solvent closer into the wall particles, effectively pushing the wall back. The additional attraction between the solvent particles and the colloidal particle would further pull the solvent particles into the wall as the colloidal particle moves through it. Consequently, 1 becomes effectively smaller for the larger radii particles for the way in which we have set up our analysis. These values are also the most sensitive parameters with respect to the quadratic versus cubic fits. The cubic corrections are a means for accounting for these additional forces, and they decrease the sensitivity of 1 with respect to the particle radius. Therefore, we feel that the cubic fits are a necessary addition to evaluate the wall position for LJ-type potentials. This aspect in particular highlights some complexities that must be considered if a quantitative analysis of hydrodynamic flow near a surface is desired.

Conclusions
We have used molecular dynamics simulations to characterize the effective particle radii and wall positions for MD simulations of WCA and LJ particles immersed in a solvent and interacting with an atomically defined wall. By analyzing the force-response curves with a continuum model, we found that this simple analysis gave very consistent results for a range of particle sizes, extending down to where the particles are the same size as the solvent. In that limiting case, we found a radius of ∼ 0.75σ for our systems with ρ = 0.8. These results do not correspond to either of the two commonly assumed alternatives for such systems: the bare radius of ∼ 0.5σ or the "hydrodynamic" radius of ∼ 1.0σ . Instead, the results correspond best to the "effective" radius defined by Balucani et al. in their analysis of the viscosity of solvent systems [24]. This is essentially based on partitioning all the volume of the simulation box between the particles as a means of defining their "effective" volume, a required step in mapping MD-derived quantities onto a continuum framework. These values will depend on the exact conditions prescribed in the simulation and should not be taken as independent fixed quantities. Similarly, the position of the "wall" in the simulations was also a variable quantity that depends on the conditions of the simulation, though the effective radii of the wall atoms are similar to that of the solvent atoms. The radius of the particle (and solvent) is a critical feature in the analysis of the hydrodynamics properties of atomistic systems (such as the viscosity and slip coefficient c). Since the "effective" radius is intermediate between the bare radius and hydrodynamic radius that have been assumed in many previous studies, this result will affect their conclusions. This will be discussed further in subsequent publications.
Another issue to consider in future studies is the effect of flow (possibly time-dependent) near the wall that would lead to changes in the fluid structure, such as density oscillations observed within a few fluid atomic diameters of the wall. The approach discussed here is appropriate for quasi-static flow, and higher-rate dynamic effects deserve a separate study. The analysis of these more complex phenomena will still depend on the ability to define the sizes of the particles that are determined in this study. Finally, while the simulations described here have been performed for a single type of fluid and fluid density for each of WCA and LJ potentials, a study of multiple cases of different fluid density and interatomic potential parameters remains a task for the future.