Calculation of surface potentials at the silica–water interface using molecular dynamics: Challenges and opportunities

Continuum-based methods are important in calculating electrostatic properties of interfacial systems such as the electric field and surface potential but are incapable of providing sufficient insight into a range of fundamentally and technologically important phenomena which occur at atomistic length-scales. In this work a molecular dynamics methodology is presented for interfacial electric field and potential calculations. The silica–water interface was chosen as an example system, which is highly relevant for understanding the response of field-effect transistors sensors (FET sensors). Detailed validation work is presented, followed by the simulated surface charge/surface potential relationship. This showed good agreement with experiment at low surface charge density but at high surface charge density the results highlighted challenges presented by an atomistic definition of the surface potential. This methodology will be used to investigate the effect of surface morphology and biomolecule addition; both factors which are challenging using conventional continuum models.


Introduction
Accurate calculation of the interfacial potential is vital to a range of fields of study such as membrane biophysics, 1,2) semiconductor-oxide 3) or oxide-oxide 4) interfacial dipole layer physics. Changes in this electrostatic potential are also crucial for interpreting, predicting and engineering the response of potentiometric biosensing platforms such as ion-sensitive field-effect transistors (IS-FETs) 5) and fieldeffect biosensors (BioFETs). 6,7) The vast majority of BioFET models incorporate a model of the electrical double layer primarily based on either the Poisson-Boltzmann equation [8][9][10][11][12][13][14][15][16][17][18] or the linearised Poisson-Boltzmann equation, [19][20][21][22] which were first developed in the early 20th century. 23) For example, the Grahame equation is often used to model DNA response, where the DNA is essentially assumed to be part of the surface charge. [24][25][26] Such predictions cannot provide quantitative accuracy for biosensor response, 21,27) and sometimes provide qualitatively incorrect results. 20) Some more recent attempts at modelling BioFET response treat the charges on the biomolecule discretely in a continuum solvent; 20,21,28,29) such approaches provide a good compromise between atomistic and continuum techniques, but neglect effects such as steric hindrance of counter-ions, biomolecule dynamics, surface morphology and water polarisation which are important for some systems.
Molecular dynamics calculation of electrostatic potential differences has the ability to provide a significant advantage over these continuum models, in that it can in principle describe biomolecular systems without system-specific empirical parametrisation (e.g., choice of a permittivity constant for a specific system of biomolecules 30) ). This is achieved via instead using standard transferable molecular dynamics simulation parameters ("forcefield") for the atomic interactions. Using molecular dynamics thus reduces the risk of model overfitting when using the simulation to predict response due to biomolecular analyte. A further advantage of molecular dynamics simulation in this field is that steric effects are modelled, which means that it can provide insight into electrostatic potential changes due to, for example, ion displacement by analyte molecules, or due to changes in the silica surface morphology.
In previous work we used ab initio dynamics to provide new insight into the energetics, dynamics and mechanisms of proton-transfer reactions involved in the charging of silica surfaces. 31) Subsequently, molecular dynamics simulations were used to investigate the silica-water-(bio) interface 7,[32][33][34][35] using the COMPASS II forcefield. 36) This forcefield has the advantage of accurate simulation of water dynamics, including H-O bond stretching. The study highlighted the importance of water polarisation and salt concentration in controlling electrostatic properties at these interfaces, 7,[32][33][34][35] and demonstrated the effects of divalent electrolyte, such as charge inversion. 7) One limitation of that work was that the computational expense of the COMPASS II forcefield permitted simulation times potentially too short to fully explore the range of conformational space involved in counter-ion adsorption=desorption processes.
In this work, the less computationally expensive INTERFACE forcefield is used 37) which permits significantly longer simulation duration (hundreds of nanoseconds) than COMPASS II, whilst also offering the benefit of being wellexperimentally validated explicitly for silica-water systems. A molecular dynamics methodology for accurate electrostatic property evaluation is validated using a silica-water system as a model system, with careful attention paid to ensuring reliable results such as the choice of electrostatics evaluation algorithm and model system setup. Proceeding from our previous conference paper, 38) which included experimental data collated for the surface potential=pH relationship at the silica-water interface, we herein present simulation data for comparison. This will provide a foundation for future work in which the model is used to investigate the effect of biomolecule addition (biosensing) and changes in surface morphology (IS-FET characteristics).
The following background is split into two subsections: firstly, background relating to accurate evaluation of electrostatics, followed by a brief overview of relevant background for charging at the silica-water interface.
1.1 Accurate electrostatics evaluation for surface systems: EW3DC A method for the rapid evaluation of Coulombic interactions between periodic images is required for the accurate calculation of the electrostatics properties of the system in molecular dynamics simulations. Molecular dynamics simulation commonly operate using a three-dimensional Ewald summation 39) to model a system that is periodic in all three dimensions ðx; y; zÞ. Whilst the method originated as far back as the 1920s, the optimal choice of electrostatic boundary conditions remains an active research area and can have a large effect on the accuracy of the result. 40) The boundary conditions associated with the 3D Ewald summation have resulted in some conceptual ambiguity because the total electrostatic energy per unit cell ("coulombic lattice sum", U) is not well defined; it is mathematically dependent upon the order in which the summation is performed, [40][41][42] i.e., it is a conditionally convergent series. A more detailed description of this topic can be found in the online supplementary data 1 at http://stacks.iop.org/JJAP/57/04FM02/mmedia, but in brief, the boundary term of the Ewald summation includes a boundary term, U IB , in the form of a series expansion in reciprocal space (k) taken at the limit k → 0. At this limit in reciprocal space, the behavior is dominated by the behaviour at infinite distance in real space. 41) When the term is neglected (U IB = 0) we arrive at what is called the "tinfoil" or "conducting" boundary condition as it can be interpreted as embedding the lattice in a conducting medium with an infinitely large dielectric constant. 40,43) The 3D Ewald summation with these tinfoil boundary conditions has become the de facto standard for handling electrostatics in molecular simulations, 40) referred to as simply EW3D in this work. While EW3D (i.e., tinfoil boundary conditions) is accurate for systems with no net dipole moment, systems which have a dipole moment over the cell (i.e., net polarization) demonstrate unphysical coupling between replicas of the interfacial dipole. As a result, different boundary conditions are required. 1) When modelling polarised interfacial systems (such as lipid bilayers or oxide-liquid interfaces), the long-range electrostatic interactions can be most accurately evaluated using a 2D Ewald (EW2D) summation where only interactions between x; y periodic images are considered, 41,42,44) however, despite much work in optimising the implementation, this technique is often prohibitively computationally expensive.
Aside from EW2D, other approaches such as local molecular field (LMF) 45) and symmetry-preserving mean field (SPMF) 46) theory have shown promise for both accurate and efficient evaluation for interfacial systems, 40) but their implementation has not yet become widespread amongst available simulation software packages. An attractive alternative that is commonly practised is to instead simulate the interfacial system as a 2D slab which is periodically repeated in the z-axis with either a vacuum layer or a solvent layer separating each slab in order to minimise unphysical periodic interactions. 41,44) Spohr increased the empty space between periodic images by increasing the size of the cell in the z-dimension (L z ), and showed that the EW3D results converge to the EW2D results with increasing L z , but that extremely large cell sizes are required. Even when L z was five times the lateral dimensions, convergence was not satisfactory. 44) To highlight this issue, Yeh and Wallqvist showed that using EW3D for a quartz=water interface led to unphysical water polarisation and inaccurate electrostatic potentials, even when L z was greater than eight times the lateral dimensions. 1) To remedy this issue, Yeh and Berkowitz proposed utilising the planar vacuum boundary condition, which they considered a corrected form of EW3D and termed EW3DC. 1,47) This term is derived by performing the summation of the infinite boundary term (U IB ) with the xy plane built up first (|x| → ∞ and |y| → ∞) before taking the limit in the z-direction (|z| → ∞). 41,47) This can implemented by simply applying an effective electric field at each frame of the simulation corresponding to the cell polarisation: where 0 , V, M z , and P z are the vacuum permittivity, the volume of the unit simulation cell, the z-component of the total dipole moment of the unit simulation cell, and the z-component of polarisation, respectively. 1,47) Importantly for the present work, without the EW3DC correction, unphysical water polarisation results in an electric field in the bulk of the water (i.e., electrostatic potential gradient in the bulk) making accurate determination of the surface potential impossible. This EW3DC methodology was further validated by Gurtovenko and Vattulainen who found excellent agreement between the electrostatic potential profile of a control system and an EW3DC corrected system. 2) Given the importance of Ewald methods for molecular simulation, other authors have extended upon the aforementioned methods. A more computationally efficient 3D Ewald summation than EW3DC has been proposed for the case of simulating electrolyte confined between infinite charged walls. 48) A efficient method for incorporating the polarisation contribution for systems confined between planar polarisable walls has also been developed. 49) When simulating systems with a net charge on the simulation cell, a further background charge correction is required because M z in Eq. (1) then becomes dependent on the choice of origin of the unit cell. 50) The Ewald 3D summation has also been recast in terms of a sum of pairwise additive interactions, which does not provide an efficient algorithm for computation but instead can be used to help physical interpretation of more efficient Ewald 3D summations. 41,51) 1.2 Surface charging at the silica-water interface Chemical reactions at the surface of silica between reactive silanol groups (Si-OH) and H + =OH − in water are thought to be the primary surface charging mechanism for hydrated silica, 52) with electrolyte effects having a measurable but less significant effect. 53) Silica is negatively charged at physiological pH (pH 6-9), with a point-of-zero charge (i.e., the pH for which the net charge of the surface is zero) of approximately pH 2-4 [54][55][56][57][58][59][60] and only gathers a positive charge at extremely acidic conditions (<pH 2). 61,62) The density of silanol groups is important in determining the surface electrodynamic properties. In this work, ∼5 OH nm −2 was used, as it has been used by other authors to successfully reproduce experimental IS-FET data, 63) is often used in IS-FET modelling 64) and has been reported as the density of fully hydrated silica. 65,66)

Forcefield and molecular dynamics setup
The INTERFACE forcefield (v1.5) was used for this molecular dynamics simulation work. 37,[67][68][69][70] One of many available forcefields, the INTERFACE forcefield is one of the only forcefields parameterised to treat charged silicawater interfaces and likely represents the state-of-the-art for accurate molecular dynamics simulation of the electrified silica-water interface. Under physiologically relevant conditions (i.e., pH ∼6-9), the surface is highly charged and most silica-water forcefields do not have parameters to describe the negatively charged silanolate groups at the surface. The INTERFACE forcefield shows excellent agreement with experiment for a large range of experimental data such as water contact angle, adsorption energy of peptides (at 20% surface ionization), water adsorption isotherms, immersion energy in water and the cell parameters of quartz. 37) The forcefield has been used to simulate protein binding, showing good agreement to measured binding isotherms. 69) The NAMD simulation software 71) was used for this work, using 2 fs timesteps, NVT Langevin dynamics at 298.15 K with the SETTLE algorithm for all hydrogen atoms. 72) The particle mesh Ewald (PME) method is used for efficient electrostatics evaluation with periodic boundary conditions, using a 1 Å grid spacing and 1 × 10 −6 tolerance and a local interaction distance of 12 Å.
Over the duration of the simulation (320 ns), several water molecules were found to evaporate into the vacuum layer and cross the z-periodic boundary and so a simple repulsive potential was employed to prevent this [Eq. (2)]. More specifically, a list of atoms which had a z coordinate greater than that of the wall (z 0 = 120 Å) was generated and updated every 100 steps (200 fs), and every step a force, F, was applied to atoms in the list: where k is a harmonic force constant of −1 kcal mol −1 Å −2 .

Molecular model
The model of the silica surface used in this work was based on the crystalline model supplied with the INTERFACE forcefield. 67) The model, shown in Fig. 1, is derived from the ð10 1Þ cleavage plane of α-cristobalite and has Q 3 isolated silanols groups of density 4.807 OH per nm −2 , with a surface area of 33.4076 × 34.8705 Å 2 and a slab thickness of ∼25 Å. Silica atoms within ∼11 Å of the bottom were rigidly constrained. The model was fully hydroxylated on both top and bottom, but in order to generate systems with surface charge, randomly selected silanol groups on the top surface were replaced with negatively charged silanolate groups as discussed in the following.
As the pH becomes increasingly alkaline, the negative surface charge density increases resulting in a more negative surface potential. The relationship between surface charge density and surface potential can be obtained experimentally by potentiometric titration of colloidal silica particles. This relationship is non-linear, but to a first order approximation can be treated as a linear relationship. In this work the linear empirical relationship presented by Emami et al. 67,68) is used in order to compare simulated systems of a given surface charge density to the corresponding effective pH. Specifi-cally, approximately 0.024 C m −2 pH −1 was used, corresponding to that measured for silica surfaces terminated by Q 3 isolated silanols groups, with a density of ∼4.7 OH per nm −2 , at an ionic strength of 0.1 to 0.3 M.
In all systems, a TIP3P water box was initialised on the surface with a height of 73 Å, and NaCl was added to produce a 300 mM solution (15 Na: 15 Cl ions, for this cell size). Six systems were considered with surface charge densities of 0, 0.0413, 0.0825, 0.1238, 0.1650, and 0.3851 C m −2 , corresponding to effective pH values of approximately 3, 4.7, 6.4, 8.1, 9.9, and 19, respectively. In all systems, electroneutrality was maintained by introduction of additional Na + ions. A summary schematic of the simulated systems in shown in Fig. 2. As described in the validation section, two repeats were performed for each system with different initial configurations of sodium ions, one in which these additional sodium ions were initialised in proximity to  the surface charge, and in the other in which they were initialised randomly throughout the bulk.

Simulation analysis methodology
In order to calculate the electric potential as a function of z, the Poisson Equation was solved by first averaging the charge density into xy slabs of z thickness δ using the VMD density profile tool, 73) and then integrating the resulting charge density using the trapezium rule, with the electric potential and field set to zero at z = 0 Å. 1,2) For the grid-spacing convergence testing figure, a 300 mM NaCl electrolyte system with silica surface charge of 0.1238 C m −2 was used. The charge density was calculated over the last 90 ns with 100 ps between each frame used for analysis (i.e., 900 frames in total). For subsequent surface potential calculation figure, the mean charge density was calculated for frames every picosecond over the last 180 ns. This was subdivided into three 60 ns windows from which the standard error across all windows (n = 3) was calculated, and multiplied by 1.96 to obtain the 95% confidence intervals shown in the figure.
Unless otherwise stated, the surface potential was calculated as the potential difference between the surface and bulk water (mean value between z = 80 to 85 Å). The bulk potential after the 180 ns equilibration period was stable over time in all simulations, showing a standard deviation of 0.6 ± 0.1 mV (SD) across all simulations therefore not a significant source of error. The surface was defined as the deepest position at which silanolate groups (SiO − ) were located (z = 26.9 to 27.1 Å). To determine the specific reference point, the minima in the potential within this interval was used.

Validation
In order to accurately provide electrostatic properties from molecular dynamics simulations, careful attention must be paid to ensuring that the electrostatics are evaluated appropriately and that the simulation parameters are high enough resolution to provide the desired level of numerical precision. In this section, validation of the EW3DC method, numerical integration grid density, vacuum spacing and simulation timescale are presented. For this work, the desired precision of the calculated surface potential is within the 1-100 mV range, as detection of pH changes in IS-FET sensors involves a change in surface potential of approximately 30 mV=pH, 56) while BioFET biosensing is typically in the 1-100 mV range. 27) 2.4.1 EW3DC validation. For this work, the EW3DC correction was implemented based on the code of Yeh et al. 1,47) and according to Eq. (1). To validate the EW3DC implementation, the force acting on a pair of test charges separated in the z dimension was calculated and shown in Fig. 3, showing agreement with the literature. 1,44,47) The correct long range limit for a 2D periodic system is such that the force approaches a constant value of F z ¼ q i q j =ð2 0 L x L y Þ 44) (Fig. 3, dotted black line), which shows excellent agreement with the EW3DC result (Fig. 3, green  diamonds). This convergence to a single value is explained by the fact that, as the distance between the charges increases, the point-charge nature (described by Coulombs law) becomes decreasingly important and the system becomes increasingly akin to the interaction between charges on two plates of a parallel-plate capacitor.
In contrast to the EW3DC and EW2D methods which emulate and utilise 2D periodic boundary conditions, Coulombs Law describes the behaviour of point charges separated in 3D space (non-periodic), and EW3D models a 3D periodic system. Figure 3 demonstrates that for periodic systems in which the z interactions are undesirable, such as interfacial systems with 2D xy periodicity, the EW3D approach for simulating long-range forces becomes highly inaccurate even with large vacuum spacing between periodic images. This result highlights that researchers should carefully consider the method of evaluating electrostatics depending on the system of interest. The EW3DC methodology is a computationally efficient choice that can simply be implemented and is presently available in many popular molecular dynamics packages (e.g., GROMACS: "ewaldgeometry:3dc" option, LAMMPS: "kspace_modify:slab" option and CASTEP: "dipole_correction").

Numerical integration convergence.
Calculation of the electrostatic potential requires double numerical integration of the charge density and thus a fine mesh is required for precise integration. The thickness of the grid spacing upon which the charge density is discretised, δ, should be suitably thin to provide sufficiently small integration error. In order to investigate convergence, the electric potential was calculated for δ = 0.001, 0.005, 0.01, 0.1, and 1 Å as shown in Fig. 4. The finest grid of 0.01 Å was sufficiently fine to accurately model the surface potential, as demonstrated by its convergence to within 2 mV of the calculated value for the 0.001 Å grid. Coarse 1 Å bins provided extremely inaccurate surface potentials, highlighting the importance of such convergence studies in accurate potential calculations.

Thermodynamic equilibrium convergence.
When calculating equilibrium properties such as the surface potential, it is important to run the simulation for sufficient time such that thermodynamic equilibrium can be reached. This was exemplified in a study of the surface potential Forces simulated using NAMD with the standard EW3D method (blue circles) and EW3DC method (green diamonds). Coulombs law is also shown for comparison (i.e., two point charges in vacuum, non-periodic system) (red line). The correct long range limit for a 2D periodic system is such that the force approaches a constant value of F z ¼ q i q j =ð2 0 L x L y Þ 44) (dotted black line), which shows excellent agreement with the EW3DC result.
across lipid membranes which found that at least 100 ns was required. 2) Values between 3 ns 7,32-34) and 20 ns 1) have been used in the literature to calculate the surface potential at the silica-water interface. In this work, we significantly extend upon these timeframes to investigate the long-term (to 320 ns) convergence of the silica-water interface surface potential as a function of time.
In order to investigate convergence as a function of time, the surface potential was calculated at time, t, based on the charge density average from t = 0 to t. The time taken to reach a stable average potential is indicative of convergence and the results, shown in the online supplementary data 2 at http://stacks.iop.org/JJAP/57/04FM02/mmedia, suggested that at least 100 ns are required for an average potential which is stable within the millivolt range.
Based on this result, simulations were performed for 320 ns, with the final 180 ns used for analysis. In order to quantify error from temporal fluctuations, the mean surface potential was calculated over three 60 ns windows and the standard error across these (n = 3) windows was calculated. In order to further ensure thermodynamic equilibrium has been achieved, two separate simulations were run for each system with different initial conditions: both simulations should reach the same surface potential if converged to thermodynamic equilibrium.

Vacuum layer validation and convergence.
A vacuum layer is used to minimise undesirable periodic interactions between periodic z images. For the EW3DC correction, elongation of the simulation cell to include empty space to the extent of 3-5 times the xy dimensions of the box has been recommended in the literature 41,44,47) which would correspond to a cell z dimension of 200 to 285 Å in this work. A z cell size of 150 Å was found to be unstable so larger cell sizes were used. Convergence testing (online supplementary data 3 at http://stacks.iop.org/JJAP/57/04FM02/mmedia) showed that z = 265, 300, and 335 Å were all stable, and surface potentials were converged to within 1 mV of each other, so 265 Å was used in this work.
Over the duration of the simulation (320 ns), several water molecules were found to evaporate into the vacuum layer and cross the z-periodic boundary. This "water noise" prohibits accurate calculation of the electrostatic potential of the system 2) and therefore, similar to other literature, 74) a simple repulsive potential was employed to prevent this [Eq. (2)].

Results and discussion
In this section, the results of the simulated surface potential as a function of pH are presented and discussed in the context of experimental data for this system. Figure 5 shows various literature experimental surface potential measurements [54][55][56][57][58][59][60] as a function of solution pH for silica surfaces from pH 2-10, shown as dotted lines. Experimental potentiometric data does not generally extend beyond pH 10 as silica is known to significantly dissolve under such alkaline conditions 75) and high pH is less technologically relevant. From this data, a non-linear response is seen which is approximately linear at physio-logical=biosensing conditions (pH 5-9) and showed a pH sensitivity of approximately 30-40 mV=pH. While most of these samples were generated using thermal oxidation or native oxide, the recent work of Sakata et al. was generated by sputtering silica and showed an anomalously high pH sensitivity around 60 mV=pH, which is the theoretical "Nernst limit" for systems in which simple proton exchange reactions occur. 76) In the popular site-binding model of pH sensitivity, the variation between different experiments is typically attributed primarily to either (a) surface hydroxyl density 64,77) or (b) variation in the acid-base equilibrium constant of the hydroxyl groups (pK b -pK a ). Changes in the pre-treatment of the silica are expected to significantly vary (a), 66) however the extent to which (b) is relevant when comparing different silica samples is unclear as determination of pK b =pK a values of silica remains a challenging area of active research. 31,62,78,79) Interference from cation exchange from the electrolyte (due to different electrolytes and ionic strengths used between experiments) can also affect the measured surface potential shift. 33,55) Figure 5 also presents the simulated surface potential as a function of surface charge density, shown as solid black lines. In agreement with experiment, with increasing surface charge density (or effective pH) up to 21% surface ionisation (0.1650 C m −2 ) the surface potential decreased, with quantitative agreement for low surface charge density. Error was quantified where possible: deviation in the calculated surface potential between repeat simulations was approximately 1 to 2 mV when the surface charge density was below 0.1650 C m −2 . Error due to the variability of the potential in the bulk water was approximately 1 mV. Such errors were therefore small compared to the changes in potential with surface charge (e.g., ∼65 mV per 0.05 C m −2 ).
The first layer of molecules at the interface is expected to play an important role in the electrostatic properties at polarised interfaces due to accumulation of counter-ions and water polarisation, therefore the first peak in the number density was used to calculate the density of the first layer of sodium ions, chloride ions and water molecules respectively, as a function of surface charge. The density of sodium ions within the first layer was found to increase linearly with surface charge density, the density of chloride was found to be within 1.96 standard error of each other for all surface charge densities, the exception being for the most charged surface (50% surface ionisation, 0.3851 C m −2 ) which showed higher chloride density due to strong electrostatic attraction to the dense sodium ion layer on the surface. Similarly, the water density within the first layer showed no clear correlation with surface charge density. Reorientational water polarisation was found to increase linearly with surface charge density up to a maximum value at 21% surface ionisation (0.1650 C m −2 ) after which no further reorientational polarisation occurred. These findings suggest sodium ion density and water polarisation dominated the electrostatic properties of the system.
In contrast to the systems with surface charge below 0.1650 C m −2 , for higher surface charge density, the surface potential began to increase with increased surface charge, particularly for the systems at effective pH 10 and above. At such high surface charge densities, there was also increased deviation in the calculated surface potential between repeat simulations (approximately 8-15 mV) suggesting poorer convergence to equilibrium. This reversal in the direction of surface potential originates from the definition used to define the "surface", and highlights a particular challenge in the calculation of surface potentials from atomistic systems: In this work the "surface" is defined as the z position corresponding to the potential minima associated with the position of negatively charged silanolate groups within the simulations (z = 26.9 to 27.1 Å). This definition, however, neglects the contribution from water molecules that penetrate beyond this layer and polarise. With increasing surface charge, two effects occurred simultaneously with led to the increasingly positive surface potential: water penetration and polarisation (a) and sodium ion accumulation (b). As for (a), water penetration and polarisation increased which meant that at high surface charge there was a large dipolar contribution from orientationally polarised water below the surface which was neglected under this definition of the surface. As for (b), increased sodium ion concentration near the surface 80) and increased electrostatic attraction between silanolate and sodium ions resulted in increased penetration of sodium ions between negatively charged silanolates which had the effect of reducing the mean charge density at the interfacial layer, and thus making the surface potential less negative. While the effect of (a) is purely methodological due to how the surface is defined, the effect of (b) is a process that can occur in physical systems but which is not described using conventional 1D=2D Poisson-Boltzmann-based or capacitor based models, which neglect such atomic-scale details and prohibit the overlap of negative surface charges and positive counterion charge.
The difficulty in determining the relevant reference point for the surface led Lee et al. to state such surface potential calculations can only be approximate. 80) In principle, however, a semi-empirical model could be designed in which the surface definition was calibrated to experiment. Gaussian smoothing can be applied to the atomistic potential distribution to obtain one which is more directly comparable to continuum (mean-field) models, 4) however this technique does not resolve the issue of identifying the experimentally relevant "surface" reference point for surface potential evaluation.
Continuum models such as the Grahame equation 81) or the Debye-Hückel equation 81) describe the surface as a charged The simulated surface potentials are presented relative to the system with no surface charge by subtraction of the surface potential for the zero surface-charge system. Literature data is displayed with the electrolyte used in the legend, and were measured by methods indicated by marker: X-ray photoelectron spectroscopy (XPS), 57) electrolyte-on-insulator (EOS), 58) impedance, 60) and ion-sensitive field-effect transistor (ISFET). [54][55][56]59) The bottom x-axis corresponds to the pH for the experimentally measured surface potential, whereas the top axis corresponds to the surface charge density used in the simulated system. These axes are aligned according to the empirical relationship between surface charge density and pH for this silica system as described in the methods section. For the simulated system, two repeats with different initial conditions were performed for each data point and showed good agreement, with differences only discernible at this scale at high surface charge (⩾ 0.1650 C m −2 ). For the simulated result, error bars shown are for each simulation, and calculated as 1.96 multiplied by standard error between three 60 ns window averages (95% confidence interval). The graph is limited to the range available for experimental data (up to pH 10). The simulated data point at extremely high surface charge density (0.1650 C m −2 ) showed a showed potential value of 396 ± 20 mV. Lines are drawn to guide the eye. layer which is separated from counter-ion charge by an infinitely thin, well-defined, surface layer. As the negative surface charge increases, the counter-ion charge near the surface increases, and the surface potential becomes more negative, as shown in the online supplementary data 4 at http://stacks.iop.org/JJAP/57/04FM02/mmedia. The differences and similarities between atomistic models of the surface potential versus continuum models of the surface potential are highlighted the schematic Fig. 6. As discussed, if the definition of the surface is the z-position of the negatively charged silanolate groups at the surface, atomistic models show the opposite trends in surface potential at high surface charge when compared to continuum models and experiment. If the definition of the surface was taken to be inside of the silica (i.e., below closest approach of water molecules), then potential differences were found to be significantly larger in magnitude than those measured in experiment.
While concept of surface potential is useful for direct comparison to experimental measurements, the electric field (i.e., the electrostatic force) is the fundamental property of the system and this property has the advantage that it can be calculated accurately from molecular dynamics simulations without resorting to semi-empirical calibration or issues defining the surface layer. For example, FETs operate by changes in the interfacial electric field modulating the electron density in a semiconducting layer below the oxide. In such a system, modelling the "field-effect" would be most directly performed by measuring the electric field (e.g., via calculating the force of a test charge below the oxide). In future work, this method will be used to investigate the field-effect in silica-water-biomolecule systems, with the goal of improving rational design of the surface chemistry of BioFETs.
To the best of the knowledge of the authors, this study represents the first application of molecular dynamics to simulate the surface charge=surface potential relationship for silica-water systems and suggests that such a model can be used to predict surface potential changes for low surface charge systems. This model is not intended to supersede continuum models for this application, indeed quantitative accuracy of the surface charge to surface potential is highly sensitive to the accuracy of the empirical relationship used for pH to surface charge density (as is also often the case when parameterising site-binding models for predicting pHpotential response), and to the definition of the surface plane. This model, however, presents three key advantages over conventional continuum models such as site-binding models 63,82,83) and the Grahame equation for biosensing. [24][25][26] Firstly, complex biomolecular systems can be incorporated with relative ease and without requirement for system specific parametrisation. For example, continuum-models of biomolecule response often require tuning of an empirical parameter specific to the experimental system such as the permittivity of the biomolecule layer. 30) Secondly, experiments show potentiometric response due to binding of neutral molecules such as pentane; 84) this cannot be explained by continuum methods but could be explained by, for example, water polarisation effects of electrolyte displacement, both of which can be modelled using this methodology. Finally, the atomic-scale effects of surface morphology (e.g., roughness, porosity, etc.) can be explicitly investigated. The results of these applications will be presented in a future publication focusing on the application of this technique.

Conclusions and future work
The surface charge=potential relationship for a silica-electrolyte system was calculated via molecular dynamics simulation for the first time. Extensive validation work was presented to ensure convergence of the integration density, vacuum spacing and equilibrium, and to ensure accurate choice of periodic boundary conditions and electrostatics evaluation methodology. Error in the calculated surface potential was quantified where possible, with the error between repeat simulations of approximately 10 mV and the error due to variation in the bulk potential of approximately 1 mV. Such validation work is essential for reliability and for enabling quantitative accuracy from molecular dynamics simulation results. An empirical relationship between surface charge density and pH was used to provide a comparison to experimental pH=surface potential measurements, and good agreement was found for low surface charge density. For high surface charge density, disagreement was found between the simulated surface potential and that found from continuum models=experiment, the reason for which was a result of the challenge of presenting a rigorous atomistic definition of the surface potential. For many applications, such as understanding the field effect in field-effect transistors, calculation of the electric field (e.g., as the forces on a test charge) would circumvent this issue while providing a more physically relevant property than the surface potential.
Molecular dynamics simulation of interfacial potentials presents several significant advantages over continuum models for modelling interfacial electrodynamics, applications of which will be further investigated in future work. Example applications include modelling the effects of  and continuum surface potential models. Red-blue dipolar rectangles represent the dipole from oriented water molecules. Continuum models have a clear definition for the surface and therefore a clear definition of the surface potential, but in contrast, the definition of a surface potential for atomistic models is more challenging. In continuum models, increased surface charge results in increased charge separation and more negative surface potentials. In contrast, for atomistic models, at high charge densities, charge separation between negative surface and positive counter-ions or polarised water can decrease. If the surface is defined as the layer of negatively charged silanolate groups, then this can lead to disagreement between the two models. complex biomolecular systems, the effect of neutral analyte on water polarisation and ion displacement, and also the effects of surface morphology on the electric field and electric potential. Given that these effects are currently poorly understood, insight into these would significantly improve the capability for rational design of effective surface chemistry for highly sensitive potentiometric biosensors.