Nanomaterials are used in many important applications such as imaging, 1-3 catalysis, 4,5 energy 6 and food production. 7,8 To control their stability in solution and interaction with other materials, nanomaterials are often functionalized with surface ligands. 9 For guiding the design of surface functionalization, it is essential to understand how surface ligands impact the physical and chemical properties of nanomaterials. A simple but important physical property is the effective size of the nanomaterial, which can be measured using different techniques, such as transmission electron microscopy (TEM), 10 dynamic light scattering (DLS), 11,12 small angle X-ray scattering (SAXS) 13 and fluctuation correlation spectroscopy (FCS). 14 Due to the different experimental conditions, the measured effective size can vary rather significantly. For example, in a recent study 15 of MTAB ((16mercaptohexadecyl)trimethylammonium bromide) functionalized gold nanoparticules, the effective radii determined from DLS for a series of particles of different gold cores are approximately three times of the values measured with TEM. Such large differences suggest that the nanoparticles aggregate to a considerable degree in solution; this could occur despite the charged nature of surface ligands because the asymmetric distribution of MTAB ligands [15][16][17][18][19] leaves a considerable degree of exposed hydrophobic areas, which promote aggregation. Taking the larger effective radii of the gold nanoparticules into consideration was important to matching the computed and measured T 2 relaxation times of surface ligands, which were observed to exhibit generally fast dynamics at the sub-nanosecond time scale; 19 whether ligand dynamics are significantly perturbed by aggregation remains to be analyzed systematically.
Another factor that might have contributed, at least in part, to the difference between TEM and DLS measured radii is that the interfacial solvent and counter ions contribute significantly to the hydrodynamic radius (R H ) inferred from the DLS experiment. For proteins in solution, 20 for example, R H differs from the radius of gyration (R g ) by 20-30% as typical for spherical objects. 21 Since (gold) nanoparticles are often significantly charged, with the bare surface charge density exceeding 0.1 C/m 2 , we expect that interfacial solvent and ion contribute significantly to both equilibrium and dynamic properties of nanoparticles, as observed for ζ potential and electrophoretic mobility in recent studies. 19,22,23 At a qualitative level, we anticipate that charged ligands at the surface slow down diffusion and therefore increase R H , as discussed for molecular ions 24,25 and recently, for small metal oxide nanoparticles. 26 These previous analyses focused on relatively low values of charge (< 5e), thus whether highly charged nanoparticles experience substantially stronger electrostatic friction (thus a larger difference between R H and R g ) remains to be clarified. Moreover, the difference between R H and R g for nanoparticles with flexible surface ligands has not been systematically analyzed. A particular question of interest is whether flexible and charged surface ligands are coupled more strongly with the surrounding solvent and counter ions than rigid ligands, due to the soft modes associated with the collective rearrangements of flexible ligands. Along this line, we note that the behaviors of nanoparticles might differ from colloids, 27 for which length scale separation is better justified. 28 In this work, we study the diffusion of three small gold nanoparticles with different surface ligands using explicit solvent molecular dynamics simulations, which allow us to explicitly compare R H and R g with molecular level of insights. Specifically, we compare three types of surface ligands that differ in length and charge states, and the results help make clear the quantitative contributions of surface charge and ligand flexibility to diffusion and therefore the value of R H . By analyzing the force-force correlation function of the nanoparticle, we discuss how compensation between electrostatic and van der Waals interactions leads to a rapid decay of the force autocorrelation function 29 and reduced total friction for the cationic particles despite their high apparent charges, as recently recognized for proteins by Matyushov and co-workers; 30 we also explicitly analyze the impact of including the first layer of solvent as part of the diffusing particle, which significantly reduces the relaxation times of electrostatic and van der Waals forces but only slightly modifies the relaxation time of the total force on the particle. These results have implications to the applicability of the Stokes-Einstein model to nanoparticles with a molecularly rough and charged surface.
In the following, we first summarize the computational models and simulation details.
Next, we present the computational results, which include finite size effect on computed diffusion constants, impacts of ligand properties and counter ion concentration on nanoparticle diffusion, and behaviors of velocity and force autocorrelation functions. Finally, we draw a few conclusions.
We study three small gold nanoparticle (Au-NP) systems, which are functionalized with short hydrophobic (-S(CH 2 ) 5 H), short cationic (-S(CH 1). A very small gold core is chosen in this work so that the molecular nature of the nanoparticle/water interface is prominent (see snapshots in Fig. 1). The structures of the Au-NP systems are set up based on the crystal structures 31,32 of similar gold clusters using the CHARMM-GUI 33 Nanomaterial Modeler. Since the small nanoparticles feature a high surface curvature, all primary amines are assumed to be protonated as they can effectively avoid each other. 34,35 a. Each gold nanoparticle contains 25 Au atoms and is functionalized with 18 ligands; thus the functionalized particles are referred to as Au 25 SC5H 18 , Au 25 SC5N 18 and Au 25 SC15N 18 . The amine groups are considered to be fully protonated, which is likely a good approximation considering the large curvature of the small Au-NPs studied here. 34,35 Explicit solvent molecular dynamics simulations are conducted for the nanoparticles in salt solution. Most production runs are conducted with a simulation box of 100 Å in length. and SC15N (see Table 1 for their chemical structures). Among those, Au 25 SC5H 18 can be considered to be largely spherical with a uniform surface, while the other two Au-NPs exhibit a significant degree of surface roughness due to the anisotropic packing of charged ligands.
For the cationic particles, simulations >100 ns using smaller boxes are first carried out to more thoroughly sample the ligand configurations; snapshots of the nanoparticles are then solvated in larger boxes for the computation of diffusion constants. For Au 25 SC15N 18 , we explicitly study the finite size effect on the computed diffusion constant by comparing results with three box sizes (80, 100 and 120 Å). The simulations are conducted with 150 mM NaCl, which is the physiological salt concentration; this choice is motivated by the consideration that gold nanoparticles are extensively used in biological studies. [36][37][38] For Au 25 SC15N 18 , a set of simulations with minimal salt (i.e., with only Cl -ions to neutralize the bare charge of the nanoparticle, which is +18 e) is also conducted for comparison. For a summary of system sizes and lengths of simulations, see Table 2.
The nanoparticles are treated with the INTERFACE force field, 39 which has been shown
to describe gold surface and nanoparticles well. [39][40][41] The salt ions, water molecules and ligands are treated using the standard models available as part of the CHARMM36/CgenFF model; [42][43][44][45] since the chemical groups in the ligands are highly similar to those in the protein force field in CHARMM36, the penalty scores are small. The long simulations for examining the solvation structure of nanoparticles and equilibrating ligand configurations have been conducted with OpenMM 7.2 46 using the NPT ensemble, while shorter simulations for the computation of velocity and force autocorrelation functions for the analysis of diffusion and friction, respectively, have been carried out with the CHARMM package 47 interfaced 47,48 with OpenMM 7.2 in the NVE ensemble to minimize the impact of thermostat on the computed correlation functions; the impact of adjusting the solvent density is also explicitly tested and shown to be a few percent on the computed diffusion constant. For the NPT simulations, the Langevin thermostat 49 with a friction constant of 1 ps -1 is used for temperature control (T =300 K), and the Monte Carlo barostat is used for pressure control (p=1 bar). For the computation of non-bonded interactions, the Particle-mesh-Ewald method 50 with a grid size ∼ 1 Å is applied to electrostatic interactions. Switching with a switch distance of 10.0 Å and a cutoff of 12.0 Å is applied to the van der Waals interactions. Shake 51 is applied to constrain bonds involving hydrogen atoms, enabling the use of a 2 fs integration time step.
While the diffusion constant can be computed using mean square displacement and the Einstein relation, we use the integration of velocity autocorrelation function since the object of interest is a single nanoparticle; previous studies 52 suggested that it might take a long time for a large Brownian solute to reach its asymptotic behavior, an effect closely related to the long-time tail in the velocity autocorrelation function. 29,53 The apparent diffusion constant (D P BC ) can be computed in the following fashion, 29,54
in which we explicitly indicate the long-time tail of C vv (t) and the box size (L) dependence of the computed diffusion constant; the choice of t * and a 0 is discussed below. Finite size corrections due to hydrodynamic considerations have been reported; 55,56 for a cube geometry in 3D, we have
in which η is the dynamic solvent viscosity. As discussed below, we observe that the correction seems to work well for the present systems and leads to a consistent estimate of D 0 largely independent of the simulation box length L.
3 Results and Discussion
Before discussing nanoparticle diffusion, we first briefly summarize the structural and dynamical features of surface ligands in the three nanoparticle systems. As illustrated by the snapshots in Fig. 1, the particle with short, hydrophobic ligands (-SC5H) features a rather homogeneous surface; by comparison, the short but cationic ligands (-SC5N) feature less homogeneous packing, while the long cationic ligands (-SC15N) feature even more heterogeneous distributions, which lead to a rather rough surface at the molecular scale. These trends are further illustrated by a spherical Voronoi tessellation 57 of the terminal heavy atoms in the ligands (Fig. 2a-c), which indicate a notable degree of ligand clustering in Au 25 SC15N 18 .
Nevertheless, the difference of the maximal and minimal eigenvalues of the moment of inertia tensor averages to about 10%, 11% and 30% of the largest eigenvalue for Au 25 SC5H 18 , Au 25 SC5N 18 and Au 25 SC15N 18 , respectively, suggesting that the three nanoparticles on av-erage do not differ significantly from a spherical shape.
As discussed in our recent analysis of similar gold nanoparticles, 19 surface ligands rearrange their conformations at the nanosecond time scale, even for the ligands that exhibit transient clustering or bundling behavior [15][16][17][18][19] due to the long hydrophobic segments. Here
The VAC functions for the various systems are shown in Fig. 3; the top panels are semi-log plots for normalized VAC functions, while the bottom panels are log-log plots for unnormalized VAC functions, which help clarify the short-time and long-time behaviors, respectively.
For a simple Brownian particle of mass m that experiences instantaneous (thus constant) friction ζ, the VAC function is expected to be a single exponential decay 29 with the short time expansion in the following form, in which the "Einstein frequency" Ω 0 is given by Ω
B T ; F is the total force on the particle due to the environment. As shown in Fig. 3a-b, the VAC functions are clearly not single exponentials for all the systems studied here, highlighting the complexity of the nanoparticle/environment interactions (see Sect.3.4). We note that even with the shortest ligand, -SC5H, the mass ratio of the nanoparticle and solvent is ∼ 6780/18 ∼ 377, thus the nanoparticle is expected to be well in the Brownian limit. 54 Plotting the unnormalized VAC functions against t (Fig. 4) indicates that Ω 0 is the highest for Au 25 SC5N 18 (also see Sect. 3.4). This can be explained as the following. For the hydrophobic particle, Au 25 SC5H 18 , the interaction with the environ-ment is expected to be weak; a small δF 2 thus leads to a lower Ω 0 . For the two cationic particles, both feature positively charged amines and therefore the magnitude of force fluctuation due to the environment is expected to be similar; however, the particle with the longer ligand is more massive by almost 25% (9595 vs. 7070 amu), leading to a smaller Ω 0 for Au 25 SC15N 18 . In fact, the initial slopes of the logarithm of VAC functions (expected to be -ζ/m for a Langevin model) for the two cationic particles have a ratio ∼0.75, which is close to their mass ratio, further supporting that the friction experienced by the cationic particles are similar in magnitude. For additional discussions of friction, see Sect.3.4, which explicitly analyzes the time dependence of particle-environment interactions.
After t reaches the characteristic "hydrodynamic time",
which is the time taken by a viscous shear wave in the solvent with dynamic viscosity η and mass density ρ m to propagate across the Brownian particle of radius R, the VAC function is expected reach the asymptotic form that decays algebraically, 29
in which a theoretical prediction 58 for a 0 is
, for the limit of heavy particle at infinite dilute concentration and when the particle's diffusion constant is negligible compared to the solvent's kinematic viscosity. For the nanoparticles studied here, the size is in the range of 5-10 Å (see below), the value of t h is in the range of 0.8 to 3.1 ps; the estimate is made with the viscosity of the TIP3P water model used here, 59 which is lower than the experimental value by almost a factor of 3. 60 Accordingly, the predicted value of
Inspection of Fig. 3c-d Å is fairly distinct from the larger ones for t ∼ 2 -5 ps (see Fig. 3c), highlighting the finite size effect on VAC function. For more thorough understanding of the long-time tail of C vv (t)
for charged particles, more extensive sampling of simulation conditions is likely required.
For the computation of diffusion constant, for the few cases where C vv (t) fits well to the form of t -3/2 , we fit the prefactor based on ln[C vv (t)] for t in the range of 2.5 and 6.0 ps, and integrate the contribution from the long-time tail analytically starting from t * =6.0 ps in Eq. 1. As shown in Table 3, the difference between this approach with truncating the integration at 10 ps is a few percent, which is adequate for the current purpose of comparing different systems. Therefore, for the rest cases and discussions, we simply conduct numerical integration of C vv (t) up to t * =10 ps.
Using Au 25 SC15N 18 as an example, we explicitly examine the finite size effect, as often observed for protein diffusion 55 and other subtle structural features that are sensitive to interfacial water structure and density 61,62 . As shown in Fig. 4a and summarized in Table 3, there is significant finite size effect; with box sizes of 80 and 100 Å, for example, the computed apparent diffusion constants, D P BC , differ by more than 10%. Considering that R g is only about 8 Å, it is remarkable that with L on the order of 100 Å, the finite size correction derived by Hummer and co-workers 55,56 is almost of the same magnitude as the computed apparent diffusion constant. Encouragingly, once the correction is included, the estimated large-box limit, D 0 , is almost identical using simulations with the three box sizes (Table 3). Although the instantaneous shape of the cationic particles is clearly not spherical (see Fig. 1), the ligand rearrangement is fast as discussed in the above subsection, thus the average shape does not deviate significantly from being spherical, explaining the apparent lack of requiring a shape correction factor 63 for the finite size correction. Comparing the diffusion constants for the three nanoparticles, it is evident that both ligand charge and length make notable contributions. While Au 25 SC5H 18 and Au 25 SC5N 18 have rather similar R g values of ∼5 Å, the cationic particle has almost 25% reduction in the diffusion constant, highlighting the limited contribution from electrostatic friction (see discussion in Sect.3.4). Increasing the ligand length by ten CH 2 units to Au 25 SC15N 18 further decreases the diffusion by another ∼15%. The diffusion constants can be converted to R H values without parentheses are obtained by truncating the integration of the velocity autocorrelation function at t=10 ps; values with parentheses are obtained by fitting C vv (t) to the form of a 0 t -3/2 in the range of 2-6 ps, and then integrating the contributions from t > 6 ps analytically using the fitted a 0 value. b. The system density is equilibrated with NPT simulations, which reduces the box size by ∼1-2 Å. c. With minimal salt rather than 150 mM NaCl.
through the Stokes-Einstein relation with the stick boundary condition (see discussion below for its applicability) as summarized in Table 3. The results indicate that even with long and flexible ligands as -SC15N, the difference between R H and R g is limited. The difference of around 4 Å corresponds to one layer of solvent surrounding the cationic groups in the ligand, as supported by the radial distribution function of water oxygen around the ligand nitrogen (Fig. 5a); the strong solvent-nitrogen peak suggests that the residence time for the first layer of solvent is expected to be longer than t h (Eq. 4) and thus should be considered to be part of the diffusing particle. For the short cationic ligand, -SC5N, the solvent distribution is almost identical (Fig. 5a), thus R H and R g differ also by ∼4 Å. For the hydrophobic particle, the solvent doesn't exhibit any significant feature around the ligand (Fig. 5a), thus R H and R g differ merely by approximately the excluded solvent size of ∼3 Å, which corresponds to the onset of the solvent-ligand radial distribution function. We note that due to the significant flexibility of surface ligands, water molecules are able to penetrate close to the surface of the nanoparticle, as was illustrated in our previous studies of similar nanoparticles 23,35,64 as well as prior simulation analyses of chromatographic systems. 65 Since both cationic particles bear a considerable degree of charge (+18 e), we also evaluate the contribution of counter ion concentration by running a set of simulations with minimal salt for Au 25 SC15N 18 . It turns out that the physiological salt concentration only perturbs the diffusion constant by a few percent (Table 3), similar to a previous analysis of metal oxide (phosphotungstate) nanoparticle. 26 The relatively small effect is likely because despite significant enhancement of chloride distribution near the nanoparticle surface, the number of chloride ions tightly associated with the nanoparticle is modest; for Au 25 SC15N 18 , it is ∼0.15 ion per ligand (Fig. 5b), which corresponds approximately to three chloride ions for the particle; the somewhat low level of association is likely due to the fact that some cationic groups are shielded due to the wrapping of the flexible ligands on the particle surface (see a snapshot in Fig. 1c), 19 which is further supported by the higher level of chloride association observed for the cationic particle with shorter ligands (Fig. 5b). This level of ion association is still possible even under the minimal salt condition, which explains the modest change of diffusion constant.
The Stokes-Einstein model [66][67][68] relates the diffusion constant of a spherical Brownian particle with its effective size and solvent viscosity,
in which c is 4 and 6 for slip and stick boundary conditions, respectively. 69,70 It is derived using a Markovian approximation to the friction kernel (see below), leading to D E = k B T ζ ; continuum analysis with different boundary conditions is then used to relate the friction ζ to liquid viscosity η and the (spherical) size of the Brownian particle. As discussed previously, 54,71 for a particle with surface roughness comparable to the solvent size, the solvent can exert torque on the diffusing particle and therefore the stick boundary condition is likely more appropriate. Therefore, for the current study, we use c = 6; indeed, if c = 4 was used, the R H will be substantially larger than R g (e.g., R H would be 12.3 Å even for the hydrophobic particle, whose R g =5.1 Å), which is not as physical.
To understand the applicability of the Stokes-Einstein model, it is thus important to examine the behavior of the translational friction, which in general is time dependent and enters the description of diffusion via the generalized Langevin equation (GLE); 29 the time dependence can be written in the form of a memory kernel, M ζ (t),
in which ζ(0) = β 3 δF 2 = mΩ 2 0 . Standard derivations using Fourier-Laplace transform of the generalized Langevin equation lead to the result, 29,30
dt is the relaxation time of the memory function.
It is generally not straightforward to compute M ζ (t), which is defined in terms of the auto-correlation function of the random force due to the environment. 29,72 It is much more straightforward to compute the auto-correlation function of the total force on the particle due to the environment from microscopic simulations (note that the total force includes both frictional and random forces in the GLE framework; 72 in the projection operator framework, 29 the random force is orthogonal to the velocity of the particle),
in which M F (t) is the normalized total force auto-correlation function with a correlation time τ F . The normalized VAC function, C vv (t)/C vv (0), M ζ (t) and M F (t) are related to each other with well-establish expressions in terms of their Fourier-Laplace transforms; 29,73 e.g.,
To further understand factors that dictate the diffusion constant, it is thus useful to evaluate τ F and also various contributions to δF 2 ; for example, since δF contains both van der Waals and electrostatic contributions, δF 2 should contain the corresponding contributions and the cross-correlation of van der Waals and electrostatic terms:
Remarkably, in previous work, 25,30 it was observed for proteins that, δF elec • δF vdW ≈ -δF 2 elec , thus,
it is of interest to investigate if this observation applies more broadly to other systems, such as the nanoparticles in the current study.
In Figure 6, we show the normalized force autocorrelation function and its decomposition at both short (< 1 ps) and long (1 ns) time scales for Au 25 SC15N 18 ; the short-time behavior , respectively, which highlight the minimal contribution from electrostatics in the latter case. In (c), the dashed lines are results when the first solvation shell of the ligand nitrogens (see Fig. 5a) is included as part of the diffusing particle.
is also shown for the other two particles for comparison. Similar to previous analysis of protein, 30 for the cationic particles, we also observe that the electrostatic and van der Waals contributions have long correlation times on the nanosecond time scale; moreover, they have strong and persistent anti-correlations, which have an equal (zero) time value of ∼-0.8. As a result, the total force auto-correlation function decays at a much faster time scale of merely ∼ 20 fs (Fig. 6a), which is at least 4-5 orders of magnitude faster than the electrostatic and van der Waals components. The cationic particle with the short ligand exhibits generally similar features, although the electrostatic and van der Waals components relax somewhat faster than Au 25 SC15N 18 ; for example, by ∼ 1 ps, the normalized electrostatic component of M F (t) decays to ∼0.5 and 0.75 for Au 25 SC5N 18 and Au 25 SC15N 18 , respectively. For the hydrophobic gold particle, the electrostatic component is much weaker and decays on the picosecond time scale (Fig. 6d). The total force correlation is dictated by the van der Waals component, which decays at the timescale of ∼50 fs; i.e., without strong electrostatic interaction, the solvent-particle force correlation in fact decays slower by a factor of two than the cationic particles, again highlighting the significance of compensation between electrostatic and van der Waals components for charged particles. The compensation between electrostatic and van der Waals forces can also be appreciated by examining the absolute magnitude of the force fluctuations, which are summarized in Table 4). The dominant contribution of the van der Waals component is evident and expected for the hydrophobic particle, for which δF 2 elec is merely ∼3-4% of δF 2 vdw . For the charged particles, strong electrostatic interaction between the particle and the environment also increase δF 2 vdw by a factor of ∼20 relative to the hydrophobic particle; due to the anti-correlation between electrostatic and van der Waals force, however, the total force fluctuations for the cationic particles are about only one order of magnitude higher than the hydrophobic particle. In fact, although the two cationic particles exhibit different values for δF 2 elec and δF 2 vdw , they have similar total force fluctuations; for both cases, we note that δF 2 tot is approximately 30-40% of δF 2 vdw , which is remarkably similar to previous observation for proteins. 30 One is tempted to estimate the diffusion constant using Eq. 8 but replacing τ ζ by τ F , which is defined as tcut 0
M F (t)dt, where t cut is an intermediate time between the molecular collision time and relaxation time of the diffusing particle; 73,74 we note that the low-frequency Fourier components of M F (t) and M ζ (t) have distinct trends, and MF (0) = ∞ 0 M F (t)dt in fact approaches 0. 73 Considering their much larger total force fluctuations, the cationic particles are expected to experience stronger frictions and thus expected to have lower diffusion constants than the hydrophobic particle. Indeed, while the strong compensation between electrostatics and van der Waals forces leads to shorter τ F , using τ F in Eq. 8 would predict that the D GLE values (see Table 4) for the cationic particles to be ∼5 times lower than the hydrophobic particle, while explicit computations of D P BC in Table 3 show only 30% difference between Au 25 SC5N 18 and Au 25 SC5H 18 . In other words, while D GLE deviates from D P BC by only ∼10% for the hydrophobic particle, the discrepancies are substantially larger for the cationic particles (compare Tables 3 and4); note that since the key quantities in Eq. 8 are evaluated using the 100 Å box without considering finite size effects, the estimated D GLE should be compared to computed apparent diffusion constant, D P BC , and the difference is about a factor of 2-3. For charged proteins, we note that Matyushov and co-workers 30 also observed that diffusion constant estimated using a simple hydrodynamic (Stokes-Einstein) model did not quantitatively agree with MD results, which found rather weak dependence of diffusion constant on protein charge (up to |Q| = 5e).
One possible origin for the significant discrepancies observed for the cationic particles is that, as discussed in Sect.3.3, the first solvation shell should be considered part of the diffusing cationic particle. To examine this possibility, we re-compute the force fluctuations and autocorrelation functions for Au 25 SC5N 18 . As shown in Fig. 6c, including the first solvation shell as part of the diffusing particle leads to substantially faster relaxation of both electrostatic and van der Waals forces, as well as their anti-correlation; accordingly, the relaxation time estimated for the total force also becomes somewhat shorter (see Table 4).
However, the magnitude of the total force fluctuation, in fact, increases (see Table 4), thus the predicted D GLE value remains substantially different from D P BC .
These observations suggest that for the cationic particles, M F (t) and M ζ (t) are sufficiently different. To match Eq. 8 with the computed D P BC values, the memory function would feature an even shorter relaxation time than the total force autocorrelation function (i.e., τ ζ < τ F ) for the cationic particles. Considering that τ F is in the range of 20-30 fs, even shorter τ ζ in the range of 6-10 fs is required; such short relaxation times for the memory function might suggest that a Markovian model (i.e., the Stokes-Einstein model) might be a reasonable approximation, although the VAC functions clearly do not follow the single exponential behavior as expected for a Langevin model (Fig. 3). Therefore, it is worthwhile explicitly investigating the memory function, 29,72 which in principle can be obtained via analysis of the VAC function 75,76 or Fourier transform of the total force auto-correlation function (see above). 29,73 Alternatively, it is also possible that many more independent trajectories with different initial ligand configurations are needed to more accurately compute the short-time behaviors of M F (t); we note that τ F was estimated to be about 5-6 fs for azurin in solution, 30 which featured larger force fluctuations ( δF 2 ∼10 4 kcal 2 •mol -2 • Å-2 ). More systematic analyses along these lines are left for future work.
To help better interpret experimental measurement of nanoparticle size, it is important to understand how their diffusion depends on the physical and chemical features of surface ligands.
In this study, using explicit solvent molecular dynamics simulations, we have examined the diffusion behaviors of three gold nanoparticles functionalized with ligands of different charge and flexibility. The results suggest that charge reduces the diffusion constant by a modest amount (25%) despite the significant value of the bare charge (+18 e), while increasing the ligand length by 10 more CH 2 units further reduces diffusion by 15%. Accordingly, the estimated R H values based on the Stokes-Einstein model do not differ significantly from R g for the cationic particles. These results help confirm that the discrepancies found in the size of similar nanoparticles by TEM and DLS measurements in Ref. 15 were due to a notable degree of aggregation under solution conditions. The limited impact of electrostatic friction on the diffusion of the nanoparticles is due in large part to the significant anti-correlation between electrostatic and van der Waals interactions between the nanoparticle and environment, as was demonstrated recently for proteins by Matyushov and co-workers, 30 and such strong compensation persists if the first shell of solvent molecules is considered to be part of the diffusing particle. As a result of the compensation, the total force autocorrelation function decays at a much faster time scale (tens of femtoseconds) than the electrostatic and van der Waals components, which relax at the nanosecond time scale; including the first shell of solvent molecules as part of the particle reduces the degree of disparity in the relaxation time scales without significantly affecting the total force autocorrelation function.
Nevertheless, using the time scale associated with the total force autocorrelation function as a proxy of the memory function relaxation time underestimates the diffusion constant for charged nanoparticles, while the approximation works fairly well for the hydrophobic nanoparticle. The impact of memory effects of the electrostatic friction on the diffusion of charged nanoparticles deserves to be further analyzed. Along this line, the effect of explicit electronic polarization on the dynamics of highly charged nanoparticles is also of major interest considering recent discussions of the topic for electrolyte solutions. 77