Hydration structure of diamondoids from reactive force fields

We optimize 21 random initial geometries of a single AD molecule (see Figure 1) surrounded by a shell of 143 water molecules (and also ten geometries of a single DI, TR and HE each surrounded by 160 water molecules) using the ORCA quantum chemistry program package ^{67, 68} version 4.2.1. We employ the revPBE functional ^{69, 70} and the def2-SVP basis set ^{71, 72} with dispersion correction ⁷³ while using loose optimization criteria. The functional is no hybrid functional and the basis set is relatively small, yet the lack of detail is a necessary compromise for optimizing bigger geometries and the method has proven itself before ^{70, 74}. One geometry-optimized structure is depicted in Figure 1 both with the full water environment (panel a) and with only sublayers of the same (panel b). The structures are then analyzed to identify properties that can be used as benchmarks to compare different MD force fields to.

After establishing our DFT benchmarks, we proceed with MD simulations. We perform MD minimization ⁷⁵ at $T=0$ K using the ReaxFF force field with the different parameter sets of Aryanpour ⁴⁷, Wang ⁶⁵ and Zhang ⁶⁶. The algorithm iteratively adjusts the atomic coordinates in such a way that, governed by the force field, all atoms move into the local potential energy minima that are closest to the initial structures. With each parameter set, we minimize all 21 previously geometry-optimized AD$+$water structures. The simulations are computationally inexpensive, however, they do not exhaustively sample the phase space, so the results can not make any statement about the statistical probability of the configurations. Regardless, as long as the DFT structures are representative of stable states in the MD simulations, then the structures should change as little as possible during the minimization. After that, we compare the MD-minimized structures to the DFT-optimized structures in terms of the benchmark properties.

We modify the Zhang parameter set such that the MD-minimized structures achieve the optimal agreement with the DFT-optimized structures with respect to the benchmarks. We choose the Zhang set for its strength in regards to water simulations ⁶⁶. A full optimization of all ReaxFF force field parameters is, however, outside of the scope of this study. For the modifications we tune the parameters associated with the QEq charge equilibration for the carbon atoms of the AD, as we will explain later. The QEq method is based on the electronegativity equalization principle which states that, at equilibrium, the electron density transfers between all atoms such that the electronegativities at all atomic sites are equal ⁴⁸. The combination of partial charges that fulfills this condition is found by minimizing the electrostatic energy of the system after every reasonable number of time steps. There are three QEq parameters: the shielding, the electronegativity and the hardness. We tune the shielding value in the range of 0.1 to 1 Å^-1 with steps of 0.025 Å^-1, the electronegativity value in the range of 1 to 12 eV with steps of 0.25 eV, and the hardness value in the range of 1 to 8 eV with steps of 1 eV. For each parameter combination, we again minimize the previously DFT-optimized structures and calculate the averaged values for the benchmark properties. To find the optimal parameter set, we minimize the cost function

$$\delta x=\sum_{i}\left[\frac{x_{i}^{{\rm DFT}}-x_{i}^{{\rm ReaxFF}}}{\sigma_{i}}\right]^{2},$$

(1)

where $x_{i}^{{\rm DFT}}$ are the values of the benchmark quantities from the DFT calculations, and $x_{i}^{{\rm ReaxFF}}$ are the their counterparts from the MD simulations. The coefficient $\sigma_{i}$ gives each benchmark quantity a weight, necessary to balance the influence that each quantity has on the cost function. This is an upper tolerance limit for the differences between the DFT and MD results. The weights are all about 10% of the maximum value in which the corresponding benchmark quantities from the DFT-optimization range. The parameter set which yields energy-minimized structures that minimize the cost function is next used to predict the water structure around the AD in bulk-water at room temperature ($T=300$ K).

Next, we prepare the simulations of an AD molecule solvated in bulk-water. We use a cubic simulation box with periodic boundary conditions. We first fill the box with 343 water molecules and equilibrate them for 10 ps in the $NVT$-ensemble at $T=300$ K and at a particularly low density of $0.5$ g/cm³ so that we can later insert an AD molecule into the box without its atoms overlapping with water molecules. We seperately pre-equilibrate the AD molecule in vacuum without periodic boundary conditions and then insert it into the water box so that the center-of-mass (COM) of the AD coincides with the center of the box. If the AD is still overlapping with water molecules, these molecules are removed from the simulation. In the subsequent equilibration in the $NPT$-ensemble, the temperature and pressure are controlled by a Nose-Hoover thermostat and barostat. They are set to $T=300$ K and $P=1$ atm, while the damping constants are set to 25 fs and 250 fs, respectively. The timestep is 1 fs. Meanwhile, the AD is retained in the box center by subtracting the COM velocity of the AD from all atoms in the system. The equilibration is finished after 3 ns, the density is constant and slightly above 1 g/cm³, and the system’s total energy is constant as well. The production simulation runs for 1.5 ns with a time step of 0.5 fs.

We start by summarizing the results of the DFT geometry optimizations of the water droplets around the DDs. As a basis for comparing the structures that we will obtain from MD simulations to the DFT-optimized structures, we define a set of benchmark quantities that characterize the alignment of the water molecules around the AD.

The first quantity is the radial density distribution of the water O and H atoms as function of the distance $r$ to the AD COM (Figure 2 a). In the first hydration shell, which stretches from $4$ to $6.2$ Å (as defined by the first minimum of the O density), the water H atoms tend to rest $0.3$ Å closer to the AD than the O atoms (see also Figure 1 b). Above $6.2$ Å, the region between the two peaks of the water density distribution forms a $1$ Å wide hydrogen-dominated layer (see also Figure 1 b). Note that both densities go to zero at high $r$, because the optimized structures are droplets and not periodic bulk water.

The second quantity is the orientation angle $\alpha$ of the water dipole moments $\vec{\mu}$ as a function of the distance to the AD COM as defined in Figure 1 c. Each dipole moment vector $\vec{\mu}$ is pointing from the O atom to the space halfway between the H atoms. The angle $\alpha$ is defined between $\vec{\mu}$ and the vector $\vec{r}$ pointing from the AD COM to the water O. One might assume then from the relative positions of the H and O density peaks in Figure 2 (a) that most water dipole moments in the first hydration shell are pointing towards the AD. However, Figure 2 (b) reveals that the opposite is true. Close to the AD, $\cos(\alpha)$ is positive, which means that $\vec{\mu}$ is pointing away from the AD. Yet, panels (a) and (b) do not contradict each other. Close to the AD, the average $\alpha$ value is 53°. As shown in the inset of Figure 2 (b), there is a range of geometries where $\alpha=53$° while one H atom is $\sim 0.26$ Å closer to the AD COM than the O atom. This explains the relative shift between the O and H peak in Figure 2 (a) despite $\vec{\mu}$ pointing away from the AD.

For the third benchmark quantity, we additionally introduce the angle $\beta$ between the normal of the water molecular plane and $\vec{r}$ (see Figure 1 c). In combination with the angle $\alpha$, we obtain a detailed description of the alignment of the water molecules in the first hydration shell. Figure 2 (c) shows the relationship between the angles $\alpha$ and $\beta$ found in the DFT-optimized structures. The map can be divided up in sections, with each section corresponding to similarly oriented water molecules. For instance, the molecules corresponding to the top and bottom edges of that map ($\alpha\approx 0.5\pi$, $\beta\approx\pm 1\pi$) are aligned parallel to the AD surface, with only small variations, and the dipole directions are oriented tangential to the AD surface. The frequent occurrence of such oriented molecules is consistent with previous Monte Carlo and MD simulations that used non polarizable force fields ^{31, 76}. According to these studies, the water molecules in the first hydration shell are tangential to the skeleton surface of the AD. However, this is the single only orientational mode reported. Our simulations yield a spectrum of eight distinct modes of orientation, which we each number with a roman numeral. Many molecules occupy a bulge-like area close to the left edge of the map. These molecules are oriented as illustrated in the inset of Figure 2 (b). The other orientational modes are illustrated in Figure 2 (c). We identify modes \Romanbar1, \Romanbar2, \Romanbar5 and \Romanbar6, because they occur with a particularly high rate in the DFT results. Conversely, the other $\alpha$-$\beta$ combinations (modes \Romanbar3, \Romanbar4, \Romanbar7 and \Romanbar8) have a remarkably low occurrence on the map. The orientations of the water molecules are also related to the presence of dangling OH bonds in the first hydration shell, which is further discussed in appendix C. The $\alpha$-$\beta$ map is a fingerprint of the water structure in the first hydration shell and serves as the third benchmark for the upcoming MD simulations.

The fourth quantity is the net charge of the AD. The net charges of the AD in the DFT calculations are determined with the ESP fit method and range from -0.356 e to +0.025 e with an average of -0.183 e. In accordance with previous investigations, there is a small electron transfer from water towards the AD ⁷⁷.

The fifth quantity is the number of C-H$\cdots$O-H₂ dipole bonds (quasi hydrogen bonds) between the AD’s H atoms and the water O atoms, where the C atoms are the donors and the water O atoms are the acceptors. We define a non-covalent bond as a dipole bond if the distance between the acceptor and the donor is lower than 3.5 Å  and the C-H$\cdots$O bond angle is smaller than 30°. The time-averaged total number of AD-water dipole bonds in the DFT results is 1.3, i.e., only approx. $8$% of the AD hydrogens are engaged in a dipole bond.

After the benchmarks have been established, we minimize the previously DFT-optimized structures using the ReaxFF force fields with the Aryanpour, Wang and Zhang parameter sets. The deviations from the DFT-optimized structures are analyzed and presented in the appendix section A.

In short, we find that the commonly used ReaxFF parameters do not yield a good agreement with the DFT structures. The density distributions and orientations of the water molecules in respect to the AD COM are not correctly reproduced, the net charge of the AD is too high, and the number of AD-water dipole bonds is overestimated by a factor of $\sim 4$. Thus, to fit the ReaxFF simulations to the DFT-optimized structures, it is necessary to modify one of the existing parameter sets. The procedure of fitting the QEq parameters, that govern the charge equilibration of the C atoms, is discussed in the appendix section B.

We find that the cost function (equation 1) is minimized by QEq parameters that compromise on all benchmark quantities.

shielding	$\gamma=0.225$ Å
electronegativity	$\chi=7.00$ eV
hardness	$\eta=7.0601$ eV

The electronegativity corresponds to 2.788 on the Pauling scale and is close to the value for an isolated C atom (2.55) and much closer than the default value from the original Zhang parameter set (1.819). The net charges and number of AD-water dipole bonds are summarized in Table 2. The water O and H atom density distributions and the orientations of the water molecules are compared in Figure 3. The modified Zhang parameter set generally overestimates the distances between the AD and the water molecules by 0.2 to 0.4 Å,   but the relative distance between the O and the H distribution is consistent with the benchmarks (panel a). The angles of the water dipole moments are slightly overestimated on average but reproduce the orientations of the water molecules quite well, including the slope of $\cos(\alpha)$ in the first hydration shell (panel b). We sometimes find molecules with a relatively high $\cos(\alpha)$ value of approx. 0.7 at 4.2 Å, but their number is not statistically significant. The relationship between the $\alpha$ and $\beta$ angles in the first hydration shell is shown in panel (c). In the form of a map, the angles are difficult to compare between MD and DFT, so we group the data into the orientational modes to help us interpret the hydration structure. Panel (d) shows the probability of a water molecule to be in one of the eight orientational modes. We see that in the MD-minimized structures it is 25% less probable to find a water molecule in modes \Romanbar1 to \Romanbar3, compared to the DFT-optimized structures, while modes \Romanbar5 to \Romanbar8 have a 25% higher probability in MD than in DFT. However, the qualitative trends from the DFT benchmarks are all maintained upon MD-minimization.

To further validate the new parameter set, we analyze whether the force field extrapolates to different DDs without loss of accuracy. We use the new parameter set to minimize DFT-optimized structures of a DI, a TR, and a HE molecule, each surrounded by 160 water molecules. For this, we apply the same DFT and MD methods as described in sections 2.1 and 2.2 (but with 10 structures per molecule). In Figure 4, we compare the MD-minimized structures to the DFT-optimized structures in terms of the water O and H density distributions. As with AD, we see a 0.2 to 0.4 Å  shift between the DFT and the MD results, but the relative distributions between O and H are the same. Furthermore, in Table 1 we compare the MD-minimized structures to the DFT-optimized structures in terms of the number of weak DD-water dipole bonds and the DD net charges. In the case of DI and TR, the results are consistent with the AD results. As for HE, however, our modified Zhang parameter set overestimates the net charge. Thus, the analysis shows that the force field trained on AD successfully extrapolates to the larger DDs DI and TR, but further modifications may be needed when attempting to simulate larger structures.

The ReaxFF force field with the modified Zhang parameters for the C atoms is now used to investigate the water structure around the DDs which are each hydrated in bulk water at $T=300$ K.

Figure 5 shows the normalized density distributions of the water O and H atoms as function of the distance to the DD COM. The results reflect typical signatures of hydrophobic solvation of smooth charged spheres, but with quantitative corrections ⁷⁸. From AD to TR, the heights of the first peaks decrease monotonically due to increasing attractive forces from the bulk. Also, each DDs O and H peaks are located at nearly the same distance from the DD’s COM. However, the density distributions around the DDs show a weaker structuring compared to the smooth spheres at the same pressure and temperature. In the case of HE, the height of the first peaks increases again due to a compensating effect of the positive HE net charge (see Table 1). The water molecules respond by a cooperative reorientation, as indicated by the H peak shifting slightly away from the position of the O peak. Further corrections have their origin in the anisotropy of the DD’s structures and small-scale interface effects.

Table 1 also shows that for all DDs the number of DD-water dipole bonds further decreases at room temperature compared to the results from MD-minimization. This seems to indicate that C-H$\cdots$O-H₂ dipole bonds are not disturbing the water hydrogen bond network surrounding the DDs. The few existing dipole bonds may even be a byproduct of the stabilization of the water hydrogen bond network very close to the interface, as we discuss in appendix C.

In Figure 6 we analyze the water structure of the first hydration shells around the DDs in more detail. The figure shows how the increase in molecular size affects the specific orientational modes of the water molecules in the first hydration shell. The blue lines show the results for the ReaxFF with the modified Zhang parameters. Every column corresponds to an orientational mode and the dots in each line stand for the DD, i.e., AD to HE from left to right. Every dot indicates the probability to find a water molecule oriented in the particular mode around the respective DD. With increasing DD size we find a strong decrease of the probabilities of modes \Romanbar1 and \Romanbar2, only small changes in modes \Romanbar3 to \Romanbar6 and a strong increase of modes \Romanbar7 and \Romanbar8. Generally, it becomes less likely for water molecules to point with their H atoms to the DD as the DD size increases. We suspect three reasons for this trend. First, the change of the DD net charge (see Table 1). Second, the decrease of the DD curvature with increasing DD size. Third, effects on small length scales at the DD-water interfaces due to the discrete partial charges.

We can isolate the first two contributions from the third one by repeating the simulations with Lennard-Jones smooth spheres (SMS) instead of the DDs. For this, we replace the DD in our simulations by a Lennard-Jones particle with size $\sigma_{\rm SMS}$ and charge $q_{\rm SMS}$. The mass $m_{\rm SMS}$ is the sum of the DD’s constituent atomic masses (e.g. $m_{\rm SMS}=386.529$ amu for AD), and $\varepsilon_{\rm SMS}=0.11$ kcal/mol, which is the same value that is used to simulate aliphatic carbon atoms in the General Amber Force Field (GAFF) ⁷⁹. We compare the results from simulations with different $\sigma_{\rm SMS}$ and $q_{\rm SMS}$, to find out which of those properties causes the strongest change of the water orientations in the first hydration shell. For the simulations with different $\sigma_{\rm SMS}$, the charge of the SMS is set to zero, but the size is varied from $\sigma_{\rm SMS}=10$ Å to $15$ Å in steps of $1$ Å. To determine which $\sigma_{\rm SMS}$ corresponds to which DD, we match the density distributions of the water O atoms around the DDs to those around the SMSs. We particularly match the position where the density of the rising first peak reaches the bulk value of 0.03344 Å^-3. For the simulations with different $q_{\rm SMS}$, the size of the hydrated SMSs is the AD size ($\sigma_{\rm SMS}=10$ Å) and the charges range between $q_{\rm SMS}=-0.3$ e and $+0.2$ e in steps of $0.1$ e, since the DD charges obtained with the modified Zhang parameters range from $-0.2979$ e (AD) to $+0.1874$ e (HE). To ensure charge neutrality, we compensate the solute charge by adding a counter charge to one randomly selected water molecule.

In Figure 6 (red and green lines) we review the eight specific orientational modes of the water molecules in the SMS simulations. In general, the probability to find a water molecule in each mode around a SMS is qualitatively similar to that of the corresponding DD from the atomistic simulation. Most of the upward and downward trends with increasing solute size also coincide, albeit with much smaller slopes, so that the water molecules are much stronger affected by the atomistic DDs than by the SMSs. We suspect that this is due to the distribution and strength of the partial charges of the DDs. To confirm this, we re-simulate the atomistic DDs but this time using GAFF and the TIP3P model for water, with all DD partial charges permanently set to zero. The response of the water molecules to the uncharged atomistic DDs (purple line in Figure 6) is just as weak as it is to the uncharged SMSs (green line). Apparently, the water molecules do not distinguish between the atomistic representation of a DD and an equivalently sized SMS, as long as both carry no charge. Neither do the water molecules respond particularly strongly to the change in the solute size as we see from the uncharged SMSs and from the uncharged GAFF model. Consequently, it must be the discrete distribution of partial charges in the modified Zhang model which strongly interact with the water molecules on small length scales.

In the following we show that there is a connection between the changes in the water structure and the electrostatic potential inside the water. We calculate the radial charge density distributions $\rho_{q}(r)$ of the water with respect to the center of mass of the solutes via

$$\rho_{q}(r)=\sum_{j}q_{j}\rho_{j}(r)$$

and from this the electric field

$$E(r)=\frac{1}{\epsilon_{r}\epsilon_{0}r^{2}}\int_{0}^{r}\rho_{q}(r^{\prime})r^{\prime 2}{\rm d}r^{\prime}$$

(2)

and electrostatic potential that the water molecules create around the solute

$$\phi(r)=-\int_{0}^{r}E(r^{\prime}){\rm d}r^{\prime}.$$

(3)

We do this for the ReaxFF simulations with the modified Zhang parameters and the simulations with the charged SMSs.

The radial distribution of the water charge density around the charged atomistic DDs is shown in Figure 7 (a). Notably, the charge density is much more sensitive to the solute’s net charge at the first positive peak than it is anywhere else, due to the dielectric screening properties of the water. While the charge of the solute increases, the first peak of the charge density gradually decreases. This is because the negative net charge of the solute attracts the positively charged H atoms and repels the negatively charged O atoms. Conversely, if the solute becomes steadily more positive, more H atoms are repelled and O atoms attracted so that the positive peak decreases. However, the influence of the solute charge on the water molecules is in competition with the intermolecular forces that stabilize the hydrogen bond network at the water interface. So, even strongly charged positive solutes will not repel all the H atoms from their water interface layers, which is why the first peak never turns negative in the tested charge range. With the decrease of the first charge density peak we also see a decrease of the electrostatic potential $\phi$ of the system at the interface compared to the bulk water. This is because the first hydration layer becomes less and less positive relative to the regions beyond, where the charge density stays almost the same height. This charge relocation causes a steeper electrostatic potential between the first hydration layer and the bulk water.

In panel (b) we see that the electrostatic properties in the simulations with the charged SMSs qualitatively agree with the charged atomistic DDs. However, the first peak of the charge density is less sensitive to the solute’s net charge and so is, consequently, the electrostatic potential: At the interface it has a total range of 1.7 V in the atomistic simulations and only 1.2 V in the SMS case. The SMS has a strong influence on the left flank of the first peak but the range of the SMS Coulomb potent-ial seems to be decayed before it reaches the peak. In contrast, the atomistic DDs exert a longer ranged influence on the water charge distribution, on the one hand because of the DD’s different shapes and sizes, and on the other hand because of the more complex interactions between individual DD and water atoms on small length scales.

Lastly, in appendix C we show that these interactions are also responsible for the particular distribution of dangling OH groups in the first hydration shells.

The goals of this work were to create a ReaxFF force field for MD simulations of DDs and to characterize their hydration structure in more detail than has ever been done before. We reviewed ReaxFF force fields from literature ^{47, 65, 66} and realized that they do not yield simulation results that are consistent with our DFT benchmarks. Thus, we proposed a ReaxFF force field with refitted parameters.

ReaxFF employs polarizable partial charges that describe many-particle interactions on small length scales more accurately than fixed partial charges do. This allowed us to shed new light on the electrostatic properties of the first hydration shell and its structure, which has been described as very unusual in experiments ¹⁶. We found a spectrum of water orientations, greatly expanding on what was known before in literature ^{31, 76}. We showed that there are almost no interfacial C-H$\cdots$O-H₂ bonds, in which we agree with experiments ¹⁶, except at very close distances to the interface, where they contribute to the stabilization of the water-water hydrogen bond network ²³. However, we also showed that there is a substantial number of dangling OH groups within the first hydration shell, which have also been observed around NDs and other hydrocarbons ^{38, 39, 40, 41, 42, 24}. Generally, we found that for understanding hydrophobic hydration ^{81, 82} it is of fundamental importance to know the effects of small surface perturbations on hydrophobic behavior. From the systematic exclusion of atomistic detail we can conclude that the hydration structure is strongly affected by the interactions between the discrete partial charges on small length scales at the DD-water interface. However, we can also conclude that hydrophobic behavior is qualitatively impervious to small surface perturbations.

Further work is underway to extend the model to simulate the interplay between dopants and electrolytes and to study surface transfer doping and the dynamics of solvated electrons generated from DDs and NDs. Our study may also motivate future investigations with a view to a better transferability of ReaxFF force fields to a wider variety of hydrocarbon/water systems.