Modeling Zeta Potential for Nanoparticles in Solution: Water Flexibility Matters

Nonequilibrium molecular dynamics simulations were performed to study the electrokinetic properties of five mainstream TIPxP water models (namely, TIP3P-FB, TIP3Pm, TIP4P-FB, TIP4P-Ew, and TIP4P/2005) in NaCl aqueous solutions in the presence of a negatively charged TiO2 surface. The impact of solvent flexibility and system geometry on the electro-osmotic (EO) mobility and flow direction was systematically assessed and compared. We found that lack of water flexibility decelerates the forward EO flow of aqueous solutions at moderate (0.15 M) or high (0.30 M) NaCl concentrations, in some special cases to such an extent that EO flow reversal occurs. Zeta potential (ZP) values were then determined from the bulk EO mobilities using the Helmholtz–Smoluchowski formula. The straight comparison against available experimental data strongly suggests that water flexibility improves the ZP determination of NaCl solutions adjacent to a realistic TiO2 surface under neutral pH conditions.


S2
Section S1. The choice of electric field strength Upon analysis of Figure S1, we notice that the electric field strength of 0.02 V/Å impacts mostly the first layer of Na + ions adsorbed on the negatively charged TiO2 surface while the distribution of Na + ions in the diffuse layer is less affected by it. Noteworthy, the difference in the first peak intensity between EMD and NEMD simulations corresponds to less than 2 Na + ions being removed from the TiO2 surface by the action of the external electric field. Therefore, we conclude that for this specific system setup (highly negatively charged TiO2 surface at moderate salt concentration) the chosen electric field strength provides fair performance and keeps an acceptable response on the stream velocity of water. Further, we observe little disturbance compared to the distribution of Na + ions from the EMD simulations. One should be aware that other systems (e.g., neutral or barely charged surface interfacing with an aqueous solution) may behave differently and the impact of a particular electric field strength has to be verified beforehand. Figure S1 shows the ion density profile of Na + ions normal to the TiO2 surface in the presence (NEMD simulation) or absence (EMD simulation) of an external electric field. Figure S1. Number density profiles of Na + ions normal to the TiO2 slab in the presence (orange) or absence (black) of an external electric field.

Section S2. The influence of reflective walls on the water and ion distribution
Upon analysis of the water density profiles in Figure S2, we identify that the water density is kept at its bulk value up to 52.5 Å from the TiO2 slab center, evidenced by a clear plateau in the water density profiles, whether the reflective walls are present or not. In both cases, we identify that the bulk-water density starts consistently dropping from 55 Å outwards, suggesting the formation of a liquid-vapor interface. In the presence of reflective walls, we identify that this liquid-vacuum interface takes place at closer distances from the TiO2 surface due to the presence of reflective walls. While for the simulation without reflective walls at the z-boundaries, we notice that the liquid-vapor region becomes more extensive outwards compared to the simulation setup where the particle confinement is imposed at the z-boundaries.
Moreover, we notice that the ion density starts deviating from its bulk-like value at shorter distances from the TiO2 surface than the water density does. Figure S2 shows that the bulk-like behavior in the ion density profiles, either in the presence or absence of reflective walls, is lost at distances farther than ~50 Å from the TiO2 slab center. Importantly, we do not identify any significant deviation in the ion concentration due to the placement of reflective walls, but a slightly lower ion density in their absence between 40 and 50 Å from the TiO2 slab center. Hence, we argue that the lowering in the water and ion density near the edge of the simulation box is caused by the formation of a liquid-vapor interface at the z-boundaries rather than simulation artifacts arising from spurious interaction between ions and the reflective walls. S4 Figure S2. 1D number density profile for water and ions along the z-direction to the TiO2 slab in the presence and absence of reflective walls. Color code: water density with reflective walls (black curve); water density without reflective wall (grey curve); ion density with reflective walls (dark blue); ion density without reflective walls (light blue); region used for zeta potential estimation (region within the two dashed red lines); reflective walls (solid red line). Figure S3. Na + Cl -RDF profiles for the different water/ion force fields adopted in this work.   Table S1. Alternative ZP estimation using the corresponding dielectric constant (Table 2)

Fluid viscosity
Fluid viscosity was calculated using the periodic perturbation method 1 . Herein, an external force a is applied to produce a velocity field v in the liquid such that ay and az are zero and ax is a function of z only. To satisfy the system periodicity and obtain a smooth and periodic velocity profile with small local shear rates, an acceleration given by a periodic cosine function is applied, which can be expressed as: in which is the acceleration amplitude (herein 0.02e-5 Å/fs 2 ) and is the z-length of the simulation box. At steady state, the velocity profile generated by this acceleration is where is the generated cosine-shaped velocity amplitude and its ensemble average is related to the shear viscosity accordingly to where , , and are the atomic mass, x-component of velocity, and z-coordinate of the particles, respectively. S16

Dielectric constant
The static dielectric constant was estimated by analyzing the fluctuation in the total dipole moment of pure water systems composed of TIPxP models in their rigid or flexible version, using the following relation: where stands for the volume of the simulation box, is the Boltzmann constant, and ⟨ 2 ⟩ − ⟨ ⟩ 2 the directional dipole moment fluctuation in (eÅ) 2 .