Comparative molecular dynamics simulations of charged solid–liquid interfaces with different water models

Abstract

Aqueous solid–liquid interfaces (SLIs) are ubiquitous in nature and technology, often hosting molecular-level processes with macroscopic consequences. Molecular dynamics (MD) simulations offer a tool of choice to investigate interfacial phenomena with atomistic precision, but there exists a large number of water models, each optimised for a different purpose. Here we compare the ability of common water models to accurately simulate the interface between a charged silica surface and an aqueous solution containing NaCl. We first compare the bulk dielectric constant of water and its dependence on salt concentration for SPC/Fw, SPC/e, TIPS3p, H₂O/DC, TIP3P-Fw, OPC3, TIP3P, TIP3P-FB, TIP3P-ST, FBA/e, and TIPS3p-PPPM, revealing large variations between models. Simulating the interface with silica for the most suitable water models (SPC/Fw, H₂O/DC, TIP3-ST and TIPS3p-PPPM) shows some intrinsic consistency with continuum predictions (Poisson–Boltzmann) whereby the free energy minima obtained from MD simulations and the analytical model are in agreement, provided the latter includes the MD-determined total charge of ions in the Stern layer and dielectric constant. This consistency stands even for water models with a dielectric constant off by 100%. For salt concentrations higher than 0.21 M NaCl, the formation of random ion–ion pairs limits the reproducibility of the MD results and the applicability of the analytical method. The results highlight the applicability of the analytical model down to the nanoscale, provided a priori knowledge of the Stern layer charge is available. The findings could have significant implications for MD simulations of SLIs, especially at charged or electrified interfaces.

---

Introduction

Aqueous solid–liquid interfaces (SLIs) are ubiquitous in nature, industry and technology, for example, in mineral growth,^1,2 the survival and bioenergetics of living organisms,^3,4 electrochemical processes and energy storage,⁵ sensing,⁶ water purification⁷ and in controlling the stability of colloidal suspension.⁸

In all cases, ions are present in the water, even if only from the natural dissociation of water molecules. Most systems contain additional metal ions, which play a significant role, especially in processes taking place at the interface with contacting solids.

Several continuum analytical models describing the behaviour of ions at SLIs have been developed and refined in the last century^9,10 and currently serve as a point of reference for most SLI-based systems. The solid usually carries an electrical potential, either naturally due to the presence of surface charges¹¹ or applied externally in technological applications.¹² This causes the ions in solution to accumulate near the SLI to ensure electroneutrality. Models often refer to the electrical double layer⁹ (EDL) for the region near the solid where the concentration of counter-ions exceeds that of the bulk solution, with triple layer models¹³ also taking into account some solvation aspect of the ions in contact with the solid. These models are remarkably successful for providing a quantitative description of ionic densities near SLIs,^9,10,14 but usually cannot predict the lateral arrangement of ions along the interface since discrete molecular details are needed. Phenomena such as charge–charge correlations, charge inversion and overscreening are also difficult to capture under continuum assumptions.^15,16

If molecular details of the water and the ions near the SLI are necessary, studies often rely on advanced experimental techniques such as scattering techniques,^17,18 sum frequency generation¹⁹ and scanning probing methods,^20–25 or computer simulations.²⁶

Computational methods are particularly well suited for investigating the short time- and length-scales of molecular processes at SLIs. Molecular dynamics (MD) simulations, in particular, are extensively used to study ions at SLIs,^27–31 including under externally applied electrical potentials,^32–34 and for comparison with experimental studies.^22,24,35,36 In MD simulations, the forces experienced by each atom are obtained from interatomic force fields,³⁷ with the trajectory of every atom tracked throughout the simulation. For simulations at SLIs, force fields are needed to describe both the solid surface and the liquid. Usually, the choice of force field for the surface is dictated by the known physical behaviour of the solid, with limited options. However, things are more complicated for the liquid with a vast variety of water models and forcefields.

The available water models can be grouped into different categories based on their resolution or their flexibility. The common water models, based on their resolution, can be categorized into 3-site,^38–41 4-site,^40,42–48 5-site^49–51 and 6-site^52,53 water models. In 3-site water models, only the water atoms are considered,^38–41 while adding a dummy atom with a negative charge to the bisector of the HOH angle brings about the 4-site model.^40,42–48 The 5-site water model has two dummy atoms alongside the three atoms of the water molecules^49–51 and finally all the dummy atoms of the lower resolution water models are present in the 6-site model.^52,53 Also, depending on how the motion of water atoms is integrated in MD simulations, the water models can be either rigid or flexible.^40,54–56 The models are usually calibrated to reproduce selected known water properties,^{40,41,54–59} with their bulk density and dielectric constant as the most common points of reference.

Correctly capturing the dielectric constant of water is particularly important for MD simulations of SLIs since the constant modulates ionic interactions in an aqueous medium. Ongoing optimisation efforts^{54,55,57,60,61} have improved the models, but there remain important challenges. For example, the salt-induced reduction of the dielectric constant is well established,^62–65 but the impact on the strength of ion–ion correlations at SLIs^22,34 is not clearly understood, nor is the salt concentration limit for the applicability of each water model. While various water models are routinely applied to describe SLIs,^{22,24,35,66,67} there does not seem to be a general consensus or a set of guidelines to select the most suitable model for a given SLI.

In this study, we systematically compare MD simulations obtained with common water models of the SLI between saline aqueous solutions and silica. We selected silica as a solid for its ubiquity in numerous natural⁶⁸ and technological systems,^69,70 including for SLI investigation.^34,71 We expect analytical predictions^9,10 to correctly describe the ionic density near the SLI for this system, hence providing a point of reference for comparison with MD simulations. The goal is to quantify similarities and differences between the predictions obtained from the different water models and establish the range of applicability in each case. We start by exploring the dielectric constant of pure and saline bulk water derived from each model.

Methods

System setup

Two different systems were utilized to compare water models in the presence of ions and at SLIs: (i) bulk water with various salt concentrations to obtain the variation in dielectric constant (Fig. 1a) and (ii) a saline aqueous solution near a charged silica surface to investigate the ionic distribution at the SLI (Fig. 1b). For bulk water simulations with each of the water models of SPC/Fw,⁵⁴ SPC/e,⁴¹ TIPS3p,⁷² H₂O/DC,⁵⁷ TIP3P-Fw, OPC3,⁵⁸ TIP3P,⁴⁰ TIP3P-FB,⁵⁹ TIP3P-ST,⁵⁵ FBA/e,⁵⁶ and TIPS3p-PPPM,⁷² 3000 molecules were placed inside a box with dimensions of 5 × 5 × 5 nm³ (see Table 1 for the associated model parameters). Subsequently, the effect of added salt on the water dielectric constant was investigated for four of the models, namely SPC/Fw, H₂O/DC, TIP3-ST and TIPS3p-PPPM. In each case, the same simulation box was used with six different NaCl concentrations: 0 M, 0.5 M, 1 M, 1.5 M, 2 M and 2.5 M.

Fig. 1 System setup for the current study. The simulation box is shown for the bulk water simulations (a) and the simulation of the SLI between water and silica (b). Added ions are not shown for clarity.

Table 1 List of the different water models used in the study, together with the relevant simulation parameters and literature reference

Model	σ O	ε O	σ H	ε H	q H	k OH/2	r OH	K θ	θ	Flexible?	LJ/Coul inner cutoff	LJ/Coul outer cutoff	PPPM	Ref.
SPC/Fw	3.16549	0.155425	0	0	0.41	529.581	1.012	37.95	113.24	Yes	9	9	Yes	43
SPC/e	3.166	0.1553	0	0	0.4238	—	1	—	109.47	No	9	9	Yes	46
TIPS3P	3.1507	0.1521	1.7753 (0.4)	0.0836 (0.046)	0.4170	—	0.9584	—	104.5	No	10	12	No	47
H2O/DC	3.1840	0.14173	0	0	0.45495	—	0.958	—	109.471	No	8	10 (12)	Yes	48
TIP3P-Fw	3.188	0.102	0	0	0.415	529.581	0.9572	34.0435	104.52	Yes	10	12	Yes	29
OPC3	3.17427	0.163406	0	0	0.447585	—	0.97888	—	109.47	No	8	8	Yes	49
TIP3P	3.188	0.102	0	0	0.415	—	0.9572	—	104.52	No	13	13	Yes	29
TIP3P-FB	3.178	0.15586	0	0	0.42422	—	1.0118	—	108.15	No	9	7	Yes	50
TIP3P-ST	3.19257	0.143858	0	0	0.42556	—	1.02030	—	108.11	No	9	9	Yes	44
FBA/e	3.1776	0.18937	0	0	0.4225	358.509	1.024	45.7696	114.7	Yes	9	9	Yes	45
TIPS3p-PPPM	3.1507	0.1521	1.7753 (0.4)	0.0836 (0.046)	0.4170	—	0.9584	—	104.5	No	10	12	Yes	47
Na^+	2.117	0.4715												78
Cl^−	5.029	0.00655												78
Na^+ (TIPS3P)	2.9399	0.08												79
Cl^− (TIPS3P)	4.0445	0.15												79

---

The same four water models were utilized to study the ionic distribution near a charged silica surface and compare the derived ionic distribution with predictions derived from the Poisson–Boltzmann theory, specifically the Gouy–Chapman^9,73,74 analytical SLI model. The silica Q3 surface was taken from the interface forcefield database,⁷⁵ with the assumption of a pH of 7 and 0.67 SiO⁻Na⁺ groups per nm². The cross-sectional dimension of the simulation box was chosen equal to the cross-section of silica (3.7 × 3.7 nm²) to assure an infinite silica with periodic boundary conditions, while a lateral size of 16.5 nm for the box was selected to avoid the interaction of the silica slab of 2.3 nm thickness with its periodic image.

All the simulations were conducted in several steps: (i) water and ions were randomly placed inside the simulation box, followed by (ii) an energy minimization and subsequently (iii) several short equilibration stages before the (iv) main simulation, which ran for 50 ns for bulk water and 30 ns for dielectric calculations and SLI simulations. The temperature and pressure were kept constant at 300 K and 1 atm using a Nose–Hoover thermostat and barostat. The electrostatic interactions were separated into short-range and long-range interactions, with the latter being calculated through particle–particle–particle–mesh (PPPM) algorithms⁷⁶ (see Table 1). The non-bonded and short-range electrostatic interactions were cutoff with the CHARMM potential cutoff function, in which the force and energy smoothly reach zero in the region between inner and outer radius values.⁷⁷ The ionic parameters for Na⁺ and Cl⁻ that solvated all the water molecules except for the TIPS3P water model are taken from the work of Yagasaki et al.⁷⁸ and calibrated for the SPC/e water molecule, while for the TIPS3P water model, the CHARMM forcefield parameters are taken from the study of Vanommeslaeghe et al.⁷⁹ For a salt concentration of 0.1 M, the total number of ions in the simulation box is 24. All the simulations were performed with the 23 Jun 2022 version of LAMMPS.⁸⁰ OVITO⁸¹ and the Matplotlib package of the Python⁸² were used for the visualisation and several C++ and Python codes were developed for analysis of the results.

Dielectric constant calculations

The water dielectric constant ε for each water model was calculated from a simulated water box with the desired salt concentration (Fig. 1a) and using eqn (1), where M, ε₀, V, k_B, T and 〈 〉 represent the total electric dipole, vacuum permittivity, volume, Boltzmann constant, temperature and ensemble average. A dedicated C++ code was written to calculate the electric dipole (eqn (1)). The ensemble averages for the M and M² parameters were calculated from their variation over the course of the simulation (eqn (2) and (3)).

The Gouy–Chapman model

The MD simulation results for the silica slab in contact with four different water models (Fig. 1b) were compared with Poisson–Boltzmann predictions through the Gouy–Chapman model^9,73,74 as described by eqn (4)–(6). In this model, the ions balance electrostatic interaction with the silica surface (Poisson equation⁹) with entropic considerations through the Boltzmann distribution. For a planar interface,⁹ this leads to an exponential decrease (respectively increase) in the concentration of counter-ions (co-ions) from the surface of the solid to the bulk ionic concentration. The characteristic length scale of the exponential evolution is characterized by the so-called Debye length λ (eqn (4)), which depends on ε, ε₀, k_B, T, the electron charge e, the Avogadro number N_A and the bulk concentration c₀ of salt in the system. Here, since a monovalent salt was used, c₀ represents the concentration of both anions and cations. The value of ε was taken from the experimentally measured values by Buchner et al.⁶² for different NaCl concentrations unless stated otherwise. The electrical potential ψⁱ₀ at the surface of silica and the potential distribution ψⁱ(z) along the direction normal to the surface are obtained through eqn (5) and (6) where σⁱ is the surface charge density, z − zⁱ₀ is the distance from the surface and the index i refers to either of the two surfaces of the silica slab. Finally, the density ρⁱ(z) of ions in the solution at any distance z from the surface of the silica is obtained through eqn (7) in which the positive and negative signs are used for cations and anions respectively.

An ion condensation parameter φ(z) is used to compare the MD results obtained with the analytical model with φ(z), defined through eqn (8), where a is chosen as the distance where the counter-ion concentration reaches its maximum value near the surface in MD simulations.

Free energy calculations

The Poisson–Boltzmann hypothesis underpinning the Gouy–Chapman model assumes that electrostatic interactions and entropy determine the free energy of the system (eqn (5) and (7)). Here, the free energy was also calculated through MD simulations for comparison. The chemical potential of an ion in bulk and at a distance z from the surface is calculated through eqn (9a) and (9b). When the system reaches equilibrium, the chemical potential is the same for all the z values (eqn (9c)), leading to an expression for the standard free energy difference ΔG⁰(z) (eqn (9d) and (10)). ΔG⁰(z) represents the amount of free energy required to move a mole of ions from the bulk liquid to distance z from the surface. Eqn (11) and (7) both rest on the Poisson–Boltzmann assumption and are hence similar. Here, comparing the free energy G(z) obtained from MD simulation with the eψⁱ₀ calculated from the analytical solution is used to test the Poisson–Boltzmann hypothesis, implicitly ignoring potential energies other than the electrostatic interaction in the Boltzmann distribution.

μ_bulk = μ₀ + kT ln c₀ (9a)

μ(z) = μ₀(z) + kT ln ρⁱ(z) (9b)

μ₀ = kT ln c₀ = μ₀(z) + kT ln ρⁱ(z) (9c)

ΔG⁰(z) = μ₀(z) − μ₀ (9d)

Results

Salt dependence of the dielectric constant for different water models

Results from the simulations show that the derived dielectric constant of bulk water strongly depends on the water model considered (Fig. 2a), with the flexible TIP3p model yielding the highest value (ε = 183.05 ± 27.66), more than twice the experimental value of 78 for pure water at the same temperature.⁸³ The TIPS3p water model which is mainly used with the CHARMM forcefield for biological systems⁶⁹ also overestimates the dielectric constant (ε = 104.68 ± 3.98), but the application of the PPPM algorithm for the calculation of long-range electrostatic interaction slightly reduces the value to ε = 97.66 ± 4.80.

Fig. 2 Dielectric constants of bulk aqueous solutions calculated from MD simulations with different water models. In pure water (a), the 11 different water models exhibit significant differences. The three models best reproducing the experimental dielectric constant (SPC/Fw, H₂O/DC and TIP3P-ST), together with the widely used TIPS3p-PPPM, are selected to further calculate the effect of added NaCl (b). The models correctly capture the decrease in dielectric constant with increasing salt concentration, but they all tend to overestimate the value, with SPC/Fw-PPPM performing best at lower salt concentrations (<0.5 M). For comparison, the experimental values⁶² are also shown.

Among the different water models considered here, SPC/Fw, H₂O/DC, TIP3P-ST, OPC3 and FBA/e yield values closest to the known water dielectric constant. We therefore selected the SPC/Fw, H₂O/DC and TIP3P-ST water models to investigate the effect of salt concentrations on their dielectric constant (Fig. 2b). We also investigated the TIPS3p-PPPM water model due to its widespread use in the field. As expected, the dielectric constant reduces upon increasing the salt concentration of the solution for all the models investigated here. The SPC/Fw model correctly replicates the experimental values up to a concentration of 0.5 M but deviates for larger concentrations. The TIP3P-ST and H₂O/DC models overestimate ε for concentrations below 0.5 M and predict similar values to those of the SPC/Fw water model at higher salt concentrations. The TIPS3p-PPPM water model consistently overestimates ε by the largest amount at all concentrations, with its prediction at 2.5 M NaCl almost twice the experimental value of ε = 44.

The significant deviations between simulated and experimental dielectric constants raise obvious questions regarding their impact on the simulated ionic distribution at interfaces, here with silica. These questions are addressed in the subsequent section.

The solid–liquid interface for different water models

The distribution of cations and anions near a negatively charged silica surface confirms that all the water models qualitatively capture the expected behaviour of the EDL near the surface (SI Fig. S1). All show a gradual decrease in the concentration of co-ions concomitant with an increase in the concentration of counter-ions near the charged silica. Decreasing the bulk salt concentration increases the distance needed for the ions’ concentration to reach the bulk concentration (Debye length), also in line with analytical predictions.

Using the counter-ion/co-ion density profiles in conjunction with eqn (11), the free energy required to move a counter- or co-ion from bulk to the position (z) can be calculated (Fig. 3). The free energy profiles for the counter-ions exhibit a minimum whose value decreases with increasing salt concentration for all water models (see SI Fig. S2 for the detailed profiles). However, the position of the minimum is consistently positioned 1.65–1.70 nm from the surface, regardless of the water model used or the salt concentration. This suggests the presence of a dense counterion layer (Stern layer) the properties of which are dominated by the counterions rather than the water molecules. Here, the Stern layer is defined as the region located between the free energy minima and the silica surface. The total charge of counter-ions and co-ions in the Stern layer was calculated for comparison with silica's surface charge. The charge of the ions in each Stern layer increases with the salt concentration for all the water models. At the highest salt concentration (0.84 M), all the models converge to the same Stern layer total charge except for TIPS3P-PPPM (Fig. 3a). The charge, around +4.5e, is slightly higher than half the surface charge of silica (−8e for each side). This was taken into account for comparison with the Gouy–Chapman model, with σⁱ taken as the effective total charge of the silica surface plus the Stern layer.

Fig. 3 Free energy of counter-ions near the silica–water interfaces obtained for different water models as a function of the solution's salt concentration. The total charge of counter- and co-ions obtained from MD simulations (a) shows some disagreement between the models at lower salt concentrations, with TIPS3P-PPPM consistently underestimating the value. The free energy of a counter-ion (cation) (b)–(e) is obtained through MD simulation using the electrostatic energy of the analytical model for SPC/Fw (b), H₂O/DC (c), TIPS3P-PPPM (d), and TIP3P-ST (e). In each case, simulation results are compared with the analytical solution calculated using the dielectric constant derived from the MD simulation (ε_MD) and that experimentally measured (ε_exp.).⁶² The errors are shown as shaded around the reported values.

To compare the MD simulation results with the analytical models, the free energy minimum obtained through the simulation (minimum of G(z) from eqn (10)) was compared with the minimum electrostatic potential energy (eψⁱ₀ from eqn (5)) of the analytical solution. Two sets of analytical equations were considered; the first with the reported water dielectric constant obtained from experimental measurements⁶² and the second using the dielectric values obtained from each water model at the specific salt concentration considered (Fig. 3b–e). Since two EDLs are present in each simulation, the value reported here is the average of the free energies from both the layers. Despite the variation in the dielectric constant (Fig. 2b) and charge of the Stern layer (Fig. 3a), the MD results exhibit a good agreement with the analytical predictions for all the water models when using the experimental dielectric constant (Fig. 3b–e) except for TIPS3P-PPPM. The agreement is better at higher salt concentrations. Predictions from the TIP3P-PPPM model (Fig. 3d) show a deviation in the range of 0.25 k_BT–0.33 k_BT from the analytical value when using the experimental dielectric constant. Interestingly, if the dielectric constant obtained from simulations is used, no noticeable difference can be seen between the analytical solution and the MD predictions. This demonstrates some intrinsic consistency of the model, despite the large error in the dielectric constant obtained.

Perhaps more intuitive, the comparison between the density profiles of counter- and co-ion near the interface obtained through MD and through analytical predictions (Fig. 4) shows a good agreement for all water models and salt concentrations, at least for larger distances from the interface (>∼3 Å) where no Stern layer is present. Closer to the surface, the analytical model anticipates a gradual decline in co-ion density, but MD results show a sudden and complete elimination of the co-ions, especially at higher concentrations (see 0.84 M in Fig. 4a–d). This effect is less pronounced at lower salt concentrations. All the water models except for TIPS3P-PPPM overestimate the counter-ion density profiles near the surface at the lowest salt concentration of 0.042 M (Fig. 4). This tendency is progressively reversed with increasing the salt concentration (see, e.g., Fig. 4a), with TIPS3P-PPPM consistently underestimating the counter-ion density close to the Stern layer (Fig. 4c, g, k, o, s).

Fig. 4 Comparison of the ion density profiles near the silica surface obtained with the four selected water models at various salt concentrations. In all cases, the simulated profiles are shown with solid lines and the equivalent analytical predictions with dashed lines. The profiles associated with the counter-ions (co-ions) appear blue (red). The SPC/Fw profiles at the bulk NaCl concentrations of 0.84 M, 0.42 M, 0.21 M, 0.1 M, and 0.042 are given in (a), (e), (i), (m), and (q), respectively. The H₂O/DC profiles at 0.84 M, 0.42 M, 0.21 M, 0.1 M, and 0.042 are given in (b), (f), (j), (n), and (r), respectively. The TIP3P-PPPM profiles at 0.84 M, 0.42 M, 0.21 M, 0.1 M, and 0.042 are given in (c), (g), (k), (o), and (s), respectively. The TIP3P-ST profiles at 0.84 M, 0.42 M, 0.21 M, 0.1 M, and 0.042 are given in (d), (h), (l), (p), and (t), respectively.

Until here, the experimentally measured value of ε at each salt concentration has been used in the analytical predictions. It is instructive to also consider analytical predictions obtained with the ε value obtained from the MD simulations. Here, just SPC/Fw and TIPS3p-PPPM water models are selected for the comparison; SPC/Fw is chosen for its more accurate presentation of the dielectric constant, while TIPS3p-PPPM offers a good point of comparison for a dielectric constant significantly deviating from the experimentally measured value. The results, obtained for 0.042 M and 0.84 M NaCl (SI Fig. S3), show that SPC/Fw slightly underestimates the counter-ion density near the silica surface at 0.84 M when using the MD calculated dielectric constant, but the TIPS3p-PPPM predictions are in good agreement with the analytical solution despite its poor prediction of ε.

The condensation parameter values, as calculated through eqn (8), are shown in Fig. 5 (see also Fig. S4 for more details). According to the Poisson–Boltzmann predictions, the condensation parameter should decrease (respectively increase) for co-ions (respectively counter-ions) when moving away from the surface of silica, with the condensation parameter converging to a fixed value beyond a specific distance, here between 2 and 3 Å from the surface (Fig. S4). For co-ions, this fixed value is systematically overestimated by MD, with a larger difference at higher salt concentrations. The MD simulations tend to agree better with analytical predictions for the counter-ions up to a concentration of 0.21 M (Fig. 4c and d, see also Fig. S3), except for the TIP3-ST model, where the agreement stops beyond 0.1 M (Fig. S4p and t).

Fig. 5 Comparison of the condensation parameter calculated analytically with the results from MD simulations with two water models: SPC/Fw and TIPS3P-PPPM. In each case, the results are given at 0.84 M NaCl: SPC/Fw (a) and TIPS3P-PPPM (b), and at 0.042 M NacL: SPC/Fw (c) and TIPS3P-PPPM (d). The results are a subset of Fig. S3 and suggest the formation of ion pairs, which are responsible for the differences between analytical and MD predictions.

At higher salt concentrations, the MD results overestimate the condensation parameter, as for co-ions. This similar behaviour for both counter and co-ions at distances larger than the Debye length suggests that it is due to the formation of ion pairs based on the distance criterion: the difference in the condensation profiles of counter- and co-ions remains fixed for both SPC/Fw (Fig. 5a) and TIPS3p-PPPM (Fig. 5b). Counting the number of pairs in simulation snapshots further supports this interpretation with 37.25 ± 0.23 pairs for SPC/Fw and 34.49 ± 0.98 for TIPS3p-PPPM at 0.84 M. This is also consistent with the concentration profiles of ion pairs at different concentrations for SPC/E (Fig. S5). In contrast, when the salt concentration is low enough to discourage the formation of ion pairs (e.g., 0.1 M yielding 0.26 ± 0.04 and 0.53 ± 0.02 pairs respectively), the two models tend to follow analytical predictions better (Fig. 4c and d).

Discussion

Solid–liquid interfaces are present in many natural and technological systems. The solid surface is usually either naturally charged or can carry an applied potential in technological applications, resulting in a reorganisation of the ions in the vicinity of the surface. MD simulations are routinely used to investigate such interfaces,^33,66 often complementing and interpreting experiment results,^{22,24,34–36} but there is no clear consensus on the best choice of water model for a given application. Previous studies have reported efforts to fine-tune the dielectric constant of different water models,^{54,55,57–61} but little has been done to ensure a dielectric constant correctly reproducing the known evolution with salt concentration. Among the water models studied here, SPC/Fw,⁵⁴ H₂O/DC,⁵⁷ TIP3-ST⁵⁵ and TIPS3p-PPPM⁷² can reproduce the experimentally measured dielectric constant of water up to a salt concentration of ∼0.5 M, but they unilaterally overestimate the dielectric constant at higher salt concentrations. The higher accuracy observed for SPC/Fw at lower salt concentrations might be due to the fact that the associated Na⁺ and Cl⁻ parameters were taken from the values calculated in SPC/e water,⁸⁴ a model very similar to SPC/Fw. Additionally, calibration of these parameters in SPC/e was conducted using the hydration free energy of ions,⁸⁴ something particularly appropriate here. However, the TIPS3p-PPPM water model, which, in combination with the CHARMM forcefield family,⁸⁵ is very common for the simulation of biological systems, overestimated the dielectric constant over the whole salt concentration range probed, despite the Na⁺ and Cl⁻ parameters taken from the CHARMM forcefield.⁸⁶ This is noteworthy considering the careful validation of the CHARMM family of forcefields and its widespread use in the literature,^86–88 suggesting that some errors might cancel each other in the calculation of electrostatic interactions: even though the compatibility of ionic parameters with the water model is important, care should be taken to avoid situations in which errors cancel each other in the forcefield calibration.

When simulating the distribution of ions near the surface of charged silica, we show that the position of the minimum of free energy (considered as the Stern layer) can be used to compare MD and analytical predictions. Using the nominal surface charge of the silica in the analytical model leads to an overestimate of the free energy with respect to the MD simulation for all the water models. However, including the ions of the Stern layer, defined as the region between the surface and the minimum of free energy, generally leads to a good agreement between the analytical and the MD results except for the TIP3P-PPPM water model.

Before discussing the TIP3P-PPPM simulations in more detail, it is worth noting that the comparison between the other models and the analytical prediction is not straightforward. The analytical prediction uses experimentally measured dielectric constants and the charge of the Stern layer (for 0.042 M–0.84 M salt concentrations). In contrast, the MD simulations rely on overestimated dielectric constants by up to 10, except for the TIPS3P-PPPM model. The agreement between the analytical and MD free energy minima leads to three observations. First, despite recent efforts to fine-tune the dielectric constant of some water models with tolerances of less than 1,^{54,55,57–61} the exact value of the dielectric constant in the water model does not noticeably affect the free energy of ionic adsorption. This might be due to the weak dependence of the free energy on the dielectric constant ( according to eqn (4) and (5)) or due to errors in MD simulations that are larger than the error introduced by the overestimation of the dielectric constant. So, at least for the increase in the accuracy of the free energy values, which matters in interface driven phenomena, the fine-tuning of water models should direct towards other properties than the dielectric constant. Second, while the Gouy–Chapman model can be applied here, there is a need for a priori knowledge of the Stern layer, at least for highly charged surfaces such as silica (0.08 C m⁻²). Third, since the main contribution to the free energy is electrostatic interactions, ion–ion correlations do not appear significant for salt concentrations of up to 0.84 M.

The predictions with the TIPS3-PPPM model tend to deviate from other models and analytical predictions. Replacing the experimentally measured dielectric constant with the TIPS3P-PPPM-deduced dielectric constant in the analytical predictions leads to a good agreement of the free energy between the MD and analytical results. This implies an opposite dependence of the Stern layer charge (and the total charge in eqn (4) which is the sum of Stern charge and the surface charge) on the square root of the dielectric constant according to eqn (4) and (5). The higher dielectric constant of TIP3P-PPPM caused an over-screening of the electrostatic interactions, leading to lower counter-ion density near the surface. The observed agreement between MD simulations and PB theory with the water model dielectric constant, even for TIP3P-PPPM, demonstrates that MD simulations reproduces the PB solution not because of a specific dielectric constant value, but because both enthalpic (electrostatic) and entropic (ionic distribution) contributions are intrinsically present in MD simulations. This agreement holds particularly in the low-concentration regime, where ion–ion correlations and excluded volume effects are minimal, and validates the use of PB as a continuum approximation of the statistical behaviour captured atomistically in MD simulations.

The simulated ionic distribution near the surface is generally in good agreement with the analytical prediction for salt concentrations of up to 0.21 M, except for the TIPS3P-PPPM model. For higher salt concentrations, care should be taken due to the formation of the ion–ion pairs not considered in the analytical model. The random nature of the ionic pair formation also increases errors in MD simulations. The fact that no clear difference could be observed between the ion distribution derived with the different water models (except for TIPS3P-PPPM) might seem surprising considering the different dielectric constants deduced here (79.83 ± 2.16, 83.67 ± 0.46 and 82.2 ± 1.6 for SPC/Fw, H₂O/DC, and TIP3-ST). Two possible reasons can be put forward: first, the uncertainty in the calculation of the ion distribution profile may dominate over a small (up to 10) change in the dielectric constant. Second, there might be a lack of convergence for the calculated dielectric constant, since it is derived through eqn (1) which converges slowly.

While insightful for SLI simulations, the findings presented in this study are not devoid of limitations. First, we used Na⁺ and Cl⁻ calibrated with the SPC/e water model for SPC/Fw, H₂O/DC and TIP3-ST water models due to the unavailability of the relevant parameter. This could be improved through a model-specific calibration of ion interaction parameters. The representation of the surface model is another limitation of the current study: the interface forcefield has been calibrated for the TIPS3p-PPPM forcefield and has cutoff values similar to the CHARMM forcefield. Switching the water model and the associated forcefield cutoff values changes the interface energy based on which the forcefield has been calibrated. However, since the focus of the current study is on the liquid part of the SLI, the silica surface can be taken as a generic charged surface, and its exact surface energy does not significantly affect the results. An additional limitation of the current study is the use of a single, non-polarizable ion parameter set which cannot capture potential variations in ion–water interactions or interfacial energetics arising from different ion models. In particular, anion polarizability can change significantly near charged surfaces. Future work using multiple ion models or explicitly polarizable force fields would help clarify these effects. The present study also uses non-polarizable water models, which treat the molecular charge distribution as fixed and therefore cannot capture field-induced polarization near highly charged interfaces. This approximation may affect the local dielectric response, water structuring, and ion hydration in the Stern layer, especially at high salt concentrations. However, while this may influence the quantitative accuracy of the simulations, the comparative trends across water models remain reliable for the purposes of this study. Finally, limited simulation time and length scale inevitably affect the accuracy of the results, something common to all MD studies. Despite these limitations, the comparative nature of the study offers novel insights into the use of different water models for MD simulations of aqueous SLIs, in particular their suitability and robustness for describing different aspects of the interface.

Conclusion

In the present study, we compare the ability of different water models commonly used in MD simulations to correctly describe the interface between a charged solid surface and a saline aqueous solution. The results show that the choice of model can significantly affect the resulting dielectric constant of the saline solution, potentially also affecting the simulated ionic distribution at the SLI. Some water models correctly predict the variation in the dielectric constant with salt concentration in a dilute regime (<0.2 M), but all models tend to fail at higher salt concentrations, likely due to ion–ion correlation effects. Comparison with the well-established Gouy–Chapman continuum formalism shows a good general agreement for dilute salt concentrations for several water models up to 0.42 M salt concentration but not for the widely adopted TIPS3p-PPPM. The current study sheds more light on the practical aspects of common water models in the context of MD investigations of aqueous SLIs. We anticipate the findings to be of interest for conducting simulations at charged or electrified interfaces.

Author contributions

MT and KV designed the problem. MT conducted the simulations and analysed the results with input from KV. MT and KV wrote the manuscript.

Conflicts of interest

The authors declare no competing financial interests.

Data availability

Data for this article, including all the simulations, the results of which are presented in the main manuscript and supplementary information (SI) figures and tables, are available at the Durham University Research Data Repository at https://collections.durham.ac.uk/files/r147429921x. Supplementary information: simulated ionic profiles at the silica interface with different water models and salt concentrations (Fig. S1); associated free energy profiles (Fig. S2); comparison between MD and analytical predictions for ionic density with SPC/Fw and TIPS3p-PPPM (Fig. S3); and ion condensation parameters at different salt concentrations (Fig. S4). See DOI: https://doi.org/10.1039/d5cp03058g.

Acknowledgements

The authors wish to acknowledge the Hamilton supercomputer at Durham University for the provision of the computational facilities and support. This project has received funding from the UK Engineering and Physical Sciences Research Council (EPSRC grant EP/S028234/1).

References

1. M. Ricci, J. J. Segura, B. W. Erickson, G. Fantner, F. Stellacci and K. Voitchovsky, Growth and Dissolution of Calcite in the Presence of Adsorbed Stearic Acid, Langmuir, 2015, 31(27), 7563–7571.
2. J. Weber, J. N. Bracco, K. Yuan, V. Starchenko and A. G. Stack, Studies of Mineral Nucleation and Growth Across Multiple Scales: Review of the Current State of Research Using the Example of Barite (BaSO₄), ACS Earth Space Chem., 2021, 5(12), 3338–3361.
3. A. Springer, V. Hagen, D. A. Cherepanov, Y. N. Antonenko and P. Pohl, Protons Migrate along Interfacial Water without Significant Contributions from Jumps between Ionizable Groups on the Membrane Surface, Proc. Natl. Acad. Sci. U. S. A., 2011, 108(35), 14461–14466.
4. B. Coste, J. Mathur, M. Schmidt, T. J. Earley, S. Ranade, M. J. Petrus, A. E. Dubin and A. Patapoutian, Piezo1 and Piezo2 Are Essential Components of Distinct Mechanically Activated Cation Channels, Science, 2010, 330(6000), 55–60,  DOI:10.1126/science.1193270.
5. I. Källquist, F. Lindgren, M.-T. Lee, A. Shavorskiy, K. Edström, H. Rensmo, L. Nyholm, J. Maibach and M. Hahlin, Probing Electrochemical Potential Differences over the Solid/Liquid Interface in Li-Ion Battery Model Systems, ACS Appl. Mater. Interfaces, 2021, 13(28), 32989–32996.
6. I. Wachta and K. Balasubramanian, Electroanalytical Strategies for Local pH Sensing at Solid–Liquid Interfaces and Biointerfaces, ACS Sens., 2024, 9(9), 4450–4468.
7. D. Deng, W. Aouad, W. A. Braff, S. Schlumpberger, M. E. Suss and M. Z. Bazant, Water Purification by Shock Electrodialysis: Deionization, Filtration, Separation, and Disinfection, Desalination, 2015, 357, 77–83.
8. M. von Smoluchowski, Zur Kinetischen Theorie Der Brownschen Molekularbewegung Und Der Suspensionen, Ann. Phys., 1906, 326(14), 756–780.
9. J. N. Israelachvili, Intermolecular and Surface Forces: Revised Third Edition, Academic press, 2011.
10. W. Trewby, M. Tavakol, Y. M. Jaques and K. Voïtchovsky, Towards Local Tracking of Solvated Metal Ions at Solid-Liquid Interfaces, Mater. Today Phys., 2024, 44, 101441.
11. S. Sjöberg, Silica in Aqueous Environments, J. Non-Cryst. Solids, 1996, 196, 51–57.
12. C. M. A. Brett; A. M. O. Brett; C. M. A. Brett and A. M. O. Brett, Electrochemistry: Principles, Methods, and Applications, Oxford University Press, Oxford, 1993.
13. C. Tournassat, Y. Chapron, P. Leroy, M. Bizi and F. Boulahya, Comparison of Molecular Dynamics Simulations with Triple Layer and Modified Gouy–Chapman Models in a 0.1 M NaCl–Montmorillonite System, J. Colloid Interface Sci., 2009, 339(2), 533–541.
14. C. Zhao, D. Ebeling, I. Siretanu, D. van den Ende and F. Mugele, Extracting Local Surface Charges and Charge Regulation Behavior from Atomic Force Microscopy Measurements at Heterogeneous Solid-Electrolyte Interfaces, Nanoscale, 2015, 7(39), 16298–16311.
15. Q. Tan, G. Zhao, Y. Qiu, Y. Kan, Z. Ni and Y. Chen, Experimental Observation of the Ion–Ion Correlation Effects on Charge Inversion and Strong Adhesion between Mica Surfaces in Aqueous Electrolyte Solutions, Langmuir, 2014, 30(36), 10845–10854.
16. C. Calero, J. Faraudo and D. Bastos-González, Interaction of Monovalent Ions with Hydrophobic and Hydrophilic Colloids: Charge Inversion and Ionic Specificity, J. Am. Chem. Soc., 2011, 133(38), 15025–15035.
17. P. Fenter, C. Park, K. Nagy and N. Sturchio, Resonant Anomalous X-Ray Reflectivity as a Probe of Ion Adsorption at Solid-Liquid Interfaces, Thin Solid Films, 2007, 515(14), 5654–5659.
18. I. C. Bourg, S. S. Lee, P. Fenter and C. Tournassat, Stern Layer Structure and Energetics at Mica-Water Interfaces, J. Phys. Chem. C, 2017, 121(17), 9402–9412.
19. S. Pullanchery, S. Kulik, B. Rehl, A. Hassanali and S. Roke, Charge Transfer across C–H⋯O Hydrogen Bonds Stabilizes Oil Droplets in Water, Science, 2021, 374(6573), 1366–1370.
20. T. Fukuma and R. Garcia, Atomic- and Molecular-Resolution Mapping of Solid-Liquid Interfaces by 3D Atomic Force Microscopy, ACS Nano, 2018, 12(12), 11785–11797.
21. S. A. Contera, K. Voïtchovsky and J. F. Ryan, Controlled Ionic Condensation at the Surface of a Native Extremophilemembrane, Nanoscale, 2010, 2(2), 222–229.
22. M. Ricci, P. Spijker and K. Voitchovsky, Water-Induced Correlation between Single Ions Imaged at the Solid–Liquid Interface, Nat. Commun., 2014, 5, 4400.
23. C. Cafolla and K. Voitchovsky, Real-Time Tracking of Ionic Nano-Domains under Shear Flow, Sci. Rep., 2021, 11(1), 1–9.
24. W. Trewby, J. Faraudo and K. Voïtchovsky, Long-Lived Ionic Nano-Domains Can Modulate the Stiffness of Soft Interfaces, Nanoscale, 2019, 11(10), 4376–4384.
25. I. Siretanu, D. Ebeling, M. P. Andersson, S. L. S. Stipp, A. Philipse, M. C. Stuart, D. van den Ende and F. Mugele, Direct Observation of Ionic Structure at Solid-Liquid Interfaces: A Deep Look into the Stern Layer, Sci. Rep., 2015, 4(1), 4956.
26. E. Bertrand, T. D. Blake and J. D. Coninck, Influence of Solid–Liquid Interactions on Dynamic Wetting: A Molecular Dynamics Study, J. Phys.: Condens. Matter, 2009, 21(46), 464124.
27. Y. Guo, D. Surblys, Y. Kawagoe, H. Matsubara, X. Liu and T. Ohara, A Molecular Dynamics Study on the Effect of Surfactant Adsorption on Heat Transfer at a Solid-Liquid Interface, Int. J. Heat Mass Transfer, 2019, 135, 115–123.
28. P. Robin, N. Kavokine and L. Bocquet, Modeling of Emergent Memory and Voltage Spiking in Ionic Transport through Angstrom-Scale Slits, Science, 2021, 373(6555), 687–691.
29. E. T. Karim, M. He, A. Salhoumi, L. V. Zhigilei and P. K. Galenko, Kinetics of Solid–Liquid Interface Motion in Molecular Dynamics and Phase-Field Models: Crystallization of Chromium and Silicon, Philos. Trans. R. Soc., A, 2021, 379(2205), 20200320.
30. Z. Li, V. G. Ruiz, M. Kanduč and J. Dzubiella, Highly Heterogeneous Polarization and Solvation of Gold Nanoparticles in Aqueous Electrolytes, ACS Nano, 2021, 15(8), 13155–13165.
31. S. K. Natarajan and J. Behler, Neural Network Molecular Dynamics Simulations of Solid–Liquid Interfaces: Water at Low-Index Copper Surfaces, Phys. Chem. Chem. Phys., 2016, 18(41), 28704–28725.
32. M. Tavakol, A. Newbold and K. Voïtchovsky, Electrified Nanogaps under an AC Field: A Molecular Dynamics Study, J. Phys. Chem. C, 2024, 128(49), 21050–21059.
33. M. Tavakol and K. Voïtchovsky, Water and Ions in Electrified Silica Nano-Pores: A Molecular Dynamics Study, Phys. Chem. Chem. Phys., 2024, 26(33), 22062–22072.
34. J. Wang, H. Li, M. Tavakol, A. Serva, B. Nener, G. Parish, M. Salanne, G. G. Warr, K. Voïtchovsky and R. Atkin, Ions Adsorbed at Amorphous Solid/Solution Interfaces Form Wigner Crystal-like Structures, ACS Nano, 2024, 18(1), 1181–1194.
35. M. Ricci, P. Spijker, F. Stellacci, J.-F. Molinari and K. Voitchovsky, Direct Visualization of Single Ions in the Stern Layer of Calcite, Langmuir, 2013, 29(7), 2207–2216.
36. C. Cafolla, T. Bui, T. T. Bao Le, A. Zen, W. J. Tay, A. Striolo, A. Michaelides, H. C. Greenwell and K. Voïtchovsky, Local Probing of the Nanoscale Hydration Landscape of Kaolinite Basal Facets in the Presence of Ions, Mater. Today Phys., 2024, 46, 101504.
37. D. C. Rapaport, The Art of Molecular Dynamics Simulation, 2nd edn, Cambridge University Press, Cambridge, 2004.
38. W. L. Jorgensen, Quantum and Statistical Mechanical Studies of Liquids. 10. Transferable Intermolecular Potential Functions for Water, Alcohols, and Ethers. Application to Liquid Water, J. Am. Chem. Soc., 1981, 103(2), 335–340.
39. H. J. C. Berendsen; J. P. M. Postma; W. F. van Gunsteren and J. Hermans, Interaction Models for Water in Relation to Protein Hydration, in Intermolecular Forces: Proc. Fourteenth Jerusalem Symp. Quantum Chem. Biochem., ed. Pullman, B., Springer Netherlands, Dordrecht, 1981, pp. 331–342.
40. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of Simple Potential Functions for Simulating Liquid Water, J. Chem., 1983, 79(2), 926.
41. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The Missing Term in Effective Pair Potentials, J. Phys. Chem., 1987, 91(24), 6269–6271.
42. J. D. Bernal and R. H. Fowler, A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions, J. Chem. Phys., 1933, 1(8), 515–548.
43. W. L. Jorgensen, Revised TIPS for Simulations of Liquid Water and Aqueous Solutions, J. Chem. Phys., 1982, 77(8), 4156–4163.
44. H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura and T. Head-Gordon, Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew, J. Chem. Phys., 2004, 120(20), 9665–9678.
45. J. L. F. Abascal, E. Sanz, R. García Fernández and C. Vega, A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice, J. Chem. Phys., 2005, 122(23), 234511.
46. J. L. F. Abascal and C. Vega, A General Purpose Model for the Condensed Phases of Water: TIP4P/2005, J. Chem. Phys., 2005, 123(23), 234505.
47. S. Izadi, R. Anandakrishnan and A. V. Onufriev, Building Water Models: A Different Approach, J. Phys. Chem. Lett., 2014, 5(21), 3863–3871.
48. S. Piana, A. G. Donchev, P. Robustelli and D. E. Shaw, Water Dispersion Interactions Strongly Influence Simulated Structural Properties of Disordered Protein States, J. Phys. Chem. B, 2015, 119(16), 5113–5123.
49. F. H. Stillinger and A. Rahman, Improved Simulation of Liquid Water by Molecular Dynamics, J. Chem. Phys., 1974, 60(4), 1545–1557.
50. M. W. Mahoney and W. L. Jorgensen, A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions, J. Chem. Phys., 2000, 112(20), 8910–8922.
51. S. W. Rick, A Reoptimization of the Five-Site Water Potential (TIP5P) for Use with Ewald Sums, J. Chem. Phys., 2004, 120(13), 6085–6093.
52. H. Nada and J. P. J. M. van der Eerden, An Intermolecular Potential Model for the Simulation of Ice and Water near the Melting Point: A Six-Site Model of H₂O, J. Chem. Phys., 2003, 118(16), 7401–7413.
53. H. Nada, Anisotropy in Geometrically Rough Structure of Ice Prismatic Plane Interface during Growth: Development of a Modified Six-Site Model of H₂O and a Molecular Dynamics Simulation, J. Chem. Phys., 2016, 145(24), 244706.
54. Y. Wu, H. L. Tepper and G. A. Voth, Flexible Simple Point-Charge Water Model with Improved Liquid-State Properties, J. Chem. Phys., 2006, 124(2), 024503.
55. Y. Qiu, P. S. Nerenberg, T. Head-Gordon and L.-P. Wang, Systematic Optimization of Water Models Using Liquid/Vapor Surface Tension Data, J. Phys. Chem. B, 2019, 123(32), 7061–7073.
56. R. Fuentes-Azcatl and M. C. Barbosa, Flexible Bond and Angle, FBA/ε Model of Water, J. Mol. Liq., 2020, 303, 112598.
57. C. J. Fennell, L. Li and K. A. Dill, Simple Liquid Models with Corrected Dielectric Constants, J. Phys. Chem. B, 2012, 116(23), 6936–6944.
58. S. Izadi and A. V. Onufriev, Accuracy Limit of Rigid 3-Point Water Models, J. Chem. Phys., 2016, 145(7), 074501.
59. L.-P. Wang, T. J. Martinez and V. S. Pande, Building Force Fields: An Automatic, Systematic, and Reproducible Approach, J. Phys. Chem. Lett., 2014, 5(11), 1885–1891.
60. S. P. Kadaoluwa Pathirannahalage, N. Meftahi, A. Elbourne, A. C. G. Weiss, C. F. McConville, A. Padua, D. A. Winkler, M. Costa Gomes, T. L. Greaves, T. C. Le, Q. A. Besford and A. J. Christofferson, Systematic Comparison of the Structural and Dynamic Properties of Commonly Used Water Models for Molecular Dynamics Simulations, J. Chem. Inf. Model., 2021, 61(9), 4521–4536,  DOI:10.1021/acs.jcim.1c00794.
61. C. Vega and J. L. F. Abascal, Simulating Water with Rigid Non-Polarizable Models: A General Perspective, Phys. Chem. Chem. Phys., 2011, 13(44), 19663–19688.
62. R. Buchner, G. T. Hefter and P. M. May, Dielectric Relaxation of Aqueous NaCl Solutions, J. Phys. Chem. A, 1999, 103(1), 1–9.
63. A. Yu Zasetsky and I. M. Svishchev, Dielectric Response of Concentrated NaCl Aqueous Solutions: Molecular Dynamics Simulations, J. Chem. Phys., 2001, 115(3), 1448–1454.
64. K. Nörtemann, J. Hilland and U. Kaatze, Dielectric Properties of Aqueous NaCl Solutions at Microwave Frequencies, J. Phys. Chem. A, 1997, 101(37), 6864–6869.
65. A. Chandra, Static Dielectric Constant of Aqueous Electrolyte Solutions: Is There Any Dynamic Contribution?, J. Chem. Phys., 2000, 113(3), 903–905.
66. I. C. Bourg and G. Sposito, Molecular Dynamics Simulations of the Electrical Double Layer on Smectite Surfaces Contacting Concentrated Mixed Electrolyte (NaCl-CaCl₂) Solutions, J. Colloid Interface Sci., 2011, 360(2), 701–715.
67. D. A. Welch, B. L. Mehdi, H. J. Hatchell, R. Faller, J. E. Evans and N. D. Browning, Using Molecular Dynamics to Quantify the Electrical Double Layer and Examine the Potential for Its Direct Observation in the In-Situ TEM, Adv. Struct. Chem. Imaging, 2015, 1(1), 1.
68. H. C. Monger and E. F. Kelly, Silica Minerals, in Soil Mineralogy with Environmental Applications, John Wiley & Sons, Ltd, 2002, pp. 611–636.
69. X. Shang, A. V. Benderskii and K. B. Eisenthal, Ultrafast Solvation Dynamics at Silica/Liquid Interfaces Probed by Time-Resolved Second Harmonic Generation, J. Phys. Chem. B, 2001, 105(47), 11578–11585.
70. S. Park, H. Lee, H.-E. Kim, H.-D. Jung and T.-S. Jang, Bifunctional Poly (l-Lactic Acid)/Hydrophobic Silica Nanocomposite Layer Coated on Magnesium Stents for Enhancing Corrosion Resistance and Endothelial Cell Responses, Mater. Sci. Eng., C, 2021, 127, 112239.
71. Y. I. Rabinovich, M. S. Esayanur, K. D. Johanson, J. J. Adler and B. M. Moudgil, Measurement of Oil-Mediated Particle Adhesion to a Silica Substrate by Atomic Force Microscopy, J. Adhes. Sci. Technol., 2002, 16(7), 887–903.
72. W. E. Reiher, Theoretical Studies of Hydrogen Bonding, PhD Thesis, Harvard University, Cambridge, MA, 1985.
73. M. Gouy, Sur La Constitution de La Charge Électrique à La Surface d’un Électrolyte, J. Phys. Theor. Appl., 1910, 9(1), 457–468.
74. D. L. Chapman, A Contribution to the Theory of Electrocapillarity, London Edinburgh Dublin Philos. Mag. J. Sci., 1913, 25(148), 475–481.
75. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Force Field and a Surface Model Database for Silica to Simulate Interfacial Properties in Atomic Resolution, Chem. Mater., 2014, 26(8), 2647–2658.
76. R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles, 1st edn, CRC Press, Boca Raton, 1988.
77. A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins, J. Phys. Chem. B, 1998, 102(18), 3586–3616.
78. T. Yagasaki, M. Matsumoto and H. Tanaka, Lennard-Jones Parameters Determined to Reproduce the Solubility of NaCl and KCl in SPC/E, TIP3P, and TIP4P/2005 Water, J. Chem. Theory Comput., 2020, 16(4), 2460–2473.
79. K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov and A. D. Mackerell Jr., CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields, J. Comput. Chem., 2010, 31(4), 671–690.
80. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS - a Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales, Comput. Phys. Commun., 2022, 271, 108171.
81. A. Stukowski, Visualization and Analysis of Atomistic Simulation Data with OVITO–the Open Visualization Tool, Modell. Simul. Mater. Sci. Eng., 2009, 18(1), 015012.
82. J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng., 2007, 9(03), 90–95.
83. D. P. Fernández, Y. Mulev, A. R. H. Goodwin and J. M. H. L. Sengers, A Database for the Static Dielectric Constant of Water and Steam, J. Phys. Chem. Ref. Data, 1995, 24(1), 33–70.
84. T. Yagasaki, M. Matsumoto and H. Tanaka, Lennard-Jones Parameters Determined to Reproduce the Solubility of NaCl and KCl in SPC/E, TIP3P, and TIP4P/2005 Water, J. Chem. Theory Comput., 2020, 16(4), 2460–2473.
85. J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grubmüller and A. D. MacKerell, CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins, Nat. Methods, 2017, 14(1), 71–73.
86. D. Beglov and B. Roux, Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations, J. Chem. Phys., 1994, 100(12), 9050.
87. T. Lewis-Atwell, P. A. Townsend and M. N. Grayson, Comparisons of Different Force Fields in Conformational Analysis and Searching of Organic Molecules: A Review, Tetrahedron, 2021, 79, 131865.
88. M. Tavakol and T. J. Vaughan, Energy Dissipation of Osteopontin at a HAp Mineral Interface: Implications for Bone Biomechanics, Biophys. J., 2022, 121(2), 228–236.

---