Vibrational Analysis of Brucite Surfaces and the Development of an Improved Force Field for Molecular Simulation of Interfaces.

We introduce a nonbonded three-body harmonic potential energy term for Mg-O-H interactions for improved edge surface stability in molecular simulations. The new potential term is compatible with the Clayff force field and is applied here to brucite, a layered magnesium hydroxide mineral. Comparisons of normal mode frequencies from classical and density functional theory calculations are used to verify a suitable spring constant (k parameter) for the Mg-O-H bending motion. Vibrational analysis of hydroxyl librations at two brucite surfaces indicates that surface Mg-O-H modes are shifted to frequencies lower than the corresponding bulk modes. A comparison of DFT and classical normal modes validates this new potential term. The methodology for parameter development can be applied to other clay mineral components (e.g., Al, Si) to improve the modeling of edge surface stability, resulting in expanded applicability to clay mineral applications.

Introduction

Clays and clay minerals are of noted interest in a variety of applied fields, including environmental remediation,[1] high-level waste disposal,[2,3] and molecular sieve catalysis,[4] because of their unique chemical and physical properties. However, interactions of clay minerals with the environment are complex and not fully understood on a molecular level. The ultrafine grain size of clay minerals, typically less than 2 μm, limits the use of many conventional analytical methods for their characterization. Additionally, the layered nature of clay minerals typically creates stacking disorder, and numerous compositional site substitutions and vacancies lead to challenges in structural analysis. In recent years, molecular simulations have provided a supplemental approach to more traditional methods for characterizing clay minerals and evaluating associated geochemical processes.[5−12] A number of atomic-scale models have been proposed that have shown excellent agreement with published experimental results.[7,13,14]

Molecular models of minerals are developed to replicate bulk structural and physical properties, typically ignoring surface properties. Fortunately, the basal (001) surfaces of many clay minerals are accurately described and stable using the force field parameters derived from only bulk behavior.[13,15] However, other clay mineral surfaces (“edge surfaces”) are not typically stable using these classical-based models; often nonsensical frayed mineral edges develop when the edge surface is simulated in contact with an aqueous environment. In some flexible models, where all atoms are allowed full translational motion, surface hydroxyl groups may migrate away from the surface. One solution to the problem of edge surface instability is adding constraints to atoms on the immediate edge, bonding or tethering these atoms in order to establish a stable surface.[16−18] Unfortunately, these unphysical modifications do not accurately capture the dynamical properties of clay edges.

The goal of the current investigation is to develop a new capability for edge hydroxyl parametrization based on the Clayff force field,[13] resulting in improved structural integrity of mineral edge structures while maintaining the flexibility of Clayff. Specifically, we have examined the (11̅0) edge surface of Mg(OH)2 (brucite) to develop a suitable Mg–O–H three-body potential energy term. This term will prevent instances of hydroxyl group distortion and edge surface disruption, particularly for molecular simulations involving clay edges in contact with an aqueous solution. Additionally, the vibrational behavior of surface (and bulk) hydroxyl groups is more accurately described. Brucite was chosen because of its simple layer structure in which all surface hydroxyl groups have identical local environments and orientations with respect to the basal plane.

Originally, Clayff was developed for classical simulations of hydrated mineral systems and is based on a parametrization involving structural, spectroscopic, and density functional theory (DFT) data.[13] Clayff is considered a nonbonded force field as most of the system energy is accounted for by Lennard-Jones and electrostatic terms, while only a small part of the potential energy is included in bonded terms; the O–H bond is the only explicit bond used in Clayff applications. The lack of metal–oxygen bonded terms in Clayff reduces the number of adjustable parameters, yet the increased number of electrostatic and Lennard-Jones interactions in a nonbonded force field is still manageable for simulation of large model systems. Indeed, the capability of modeling large-scale and more complex clay systems, often involving over 100,000 atoms, has been demonstrated.[10,19−22]

In order to preserve the flexibility of Clayff, while introducing a new term that can tether hydroxyl groups to metal atoms without introducing unnecessary bonded terms, we have introduced a nonbonded three-body term for Mg–O–H bending interactions. There are no explicit bonds between Mg and O atoms, such that the Lennard-Jones and electrostatic energies originally parametrized for Mg–O are preserved. We have chosen a simple harmonic functional form Ebend = k(θ – θ0)2, where k is a spring constant and θ0 is the equilibrium angle, with parameters fit to structural and vibrational properties of brucite models derived by DFT calculations. We note that this three-body term does not absolutely prevent hydroxyl groups from leaving the surface during a molecular dynamics (MD) simulation, but significant nonbonded attraction would be required to overcome the harmonic Mg–O–H energy term for such an event to occur. Our simulations of brucite surfaces in an aqueous environment indicate that no hydroxyl groups leave the surface. Surface hydroxyl groups unrealistically detach from the surface when the three-body term is excluded.

This new term is compatible with the previous Clayff energy terms and is applicable to all Mg–O–H interactions (within a prescribed cutoff distance) in the brucite structure. In the original parametrization of Clayff,[13] a nonstandard nonbonded three-body term was used; however, it was found to be impractical in many large-scale simulations because of difficulties in its implementation. A nonbonded three-body term was recently developed for Mg–O–H interactions on the (11̅0) hydrotalcite (Mg6Al2(OH)16CO3·4H2O) surface; however, this implementation required parametrization of several other Clayff terms for application to the specific hydrotalcite structure.[23] In this study, we also replace the O–H bonded interaction with a recently developed Morse potential.[24] The Morse potential provides a description of the O–H bond stretch interaction that is more accurate than a harmonic function that describes an equivalent potential energy response with either bond compression or bond extension.

Methods

The brucite model is based on the published crystal structure for Mg(OH)2.[25] Brucite is comprised of two-dimensional layers of single sheets of hydroxylated edge-sharing Mg–O octahedra (Figure 1). Models of brucite surfaces were derived from the bulk structure using the Materials Studio software (Accelrys Inc.). The basal surface was created by cleaving between layers along the (001) plane. The edge surface was created by cleaving along the (11̅0) plane such that two equivalent surfaces were formed, where each surface magnesium atom was identically coordinated to five oxygen atoms. Bulk models were constructed from 2 × 2 × 2 (6.43 × 6.43 × 9.63 Å3, 40 atoms) and 7 × 8 × 5 (22.59 × 22.36 × 23.47 Å3, 1400 atoms) supercells for DFT calculations and classical simulations, respectively. Corresponding models for the (001) basal surface consisted of 2 × 2 × 2 (6.43 × 6.43 × 29.63 Å3, 40 atoms, 35.81 Å2 surface area for each surface) and 7 × 8 × 2 (22.59 × 22.36 × 50.00 Å3, 560 atoms, 505.11 Å2 surface area) supercells, while those for the edge surface consisted of 3 × 4 × 2 (9.65 × 30.232 × 9.63 Å3, 120 atoms, 92.93 Å2 surface area) and 7 × 5 × 7 (22.55 × 23.01 × 60.00 Å3, 1435 atoms, 518.88 Å2 surface area) supercells. Structures for MD simulations were orthogonalized to facilitate analysis. Each surface model consisted of a single slab of brucite with a vacuum gap of 20 Å and two surfaces that were allowed to completely relax. Larger vacuum gaps were tested for the DFT structures, but only small changes in frequencies were seen.

Figure 1

Brucite (Mg(OH)2) structures: (A) bulk (ac plane) and (B) (11̅0) edge surface (Mg, O, and H atoms are represented by green, red, and white spheres, respectively).

The DMol3 code[26] was used for localized numerical orbital DFT calculations of periodic cells. Cell parameters and atomic positions were optimized using the GGA/PW91 functional[27] with polarization functions on heavy atoms using the double numerical plus d-functions (DND) basis set (comparable in size to Gaussian 6-31G*). Evaluation of the mass-weighted Hessian matrix, followed by normal mode analysis, was performed for the final optimized structures at 0 K to yield harmonic frequencies and the vibrational density of states. Structures were optimized with a 2 × 10–5 Ha energy tolerance for convergence. Table 1 shows that the cell parameters obtained from DFT optimization are in reasonably good agreement with experiment[25] and similar DFT methods.[28] As expected, our DFT results show a small systematic overestimation of the experimental cell parameters (2.3% and 1.0% for a and c, respectively), as would be anticipated for a generalized gradient approximation (GGA) exchange-correlation functional.[29] The good comparison between DFT and experiment suggests that our DFT method is sufficient for our main purpose: providing accurate structural and vibrational results for both Mg–O bond lengths and Mg–O–H angles. Although empirical dispersion corrections to our DFT method would affect nonbonded interactions (e.g., layer–layer), such a correction is unlikely to strongly affect Mg–O–H bending modes or O–H stretch modes relevant to this work. A recent comparison of vibrational frequencies of the benzene monomer using the B3LYP and dispersion-corrected B3LYP-D methods showed almost no difference in monomer bend or stretch modes when dispersion corrections were included.[30] Similar DFT methods have been used to study structure and vibrations of similar surfaces and clusters with edge or surface hydroxyl groups, including silica,[31] iron oxyhydroxides, and magnesium oxide.[32,33]

Table 1

Comparison of Brucite Cell Lengths (Angstroms), Angles (Degrees), and Volumes (Cubic Angstroms) between Calculated and Experimental Resultsa

	exptl	DFT	CLAYFF	NB3Bb
a	3.142	3.215	3.229	3.237
b	3.142	3.215	3.229	3.237
c	4.766	4.815	4.696	4.680
α	90.0	90.0	90.0	90.0
β	90.0	90.0	90.0	90.0
γ	120.0	120.0	120.0	120.0
V	40.75	43.10	42.40	42.46

Cell length parameters were calculated for 2 × 2 × 2 supercells then scaled to a single unit cell for direct comparison with experimental values.

NB3B refers to Clayff with the additional three-body Mg–O–H term.

DFT normal modes were visualized within the Materials Studio graphical interface. Modes were divided into three categories (bulk, surface, and mixed) based on the observed motion of the atoms.

In order to accurately determine parameters for the nonbonded three-body potential, the GULP software[34] was used to optimize parameter fits to the results of the DFT calculations. While GULP has been previously used to perform fits to potential energy surfaces,[35] this has usually been carried out by fitting to a series of energies and/or forces for different structures. Although the inclusion of vibrational frequencies in least-squares fits was also in principle possible, this required the mode number to be specified at the outset and made no allowance for the possibility that the order of the frequencies might change during the fit. In the present work we have extended the capability of the fitting algorithm by including the vibrational frequencies as observables with normal mode projection. Here the eigenvectors from the DFT phonon calculations were used to determine the corresponding vibration from the force field that has the maximum overlap with the mode of interest. In this manner it was possible to ensure that the correct vibrational modes that correspond to the Mg–O–H angle bend were properly fitted. During the least-squares refinement of the three-body force constant, both the structure and the vibrational modes were used as observables to determine the value that would yield the best compromise between reproducing the DFT structures and vibrational properties. Specifically, relaxed fitting[36] was used in which the displacements of the structure on relaxation are used to define the sum of squares during refinement. This approach has the advantage that the vibrational modes were always computed for a valid harmonic minimum, as opposed to computing the Hessian for the initial unoptimized structure.

GULP was also used to optimize bulk cell parameters (Table 1) and atomic positions using the original Clayff parametrization (including the trioctahedral version of the Morse O–H term[24]), as well as the updated parameter set, which adds the new three-body term that will hereafter be referred to as NB3B. Not surprisingly, cell parameters are largely unaffected by the short-range three-body term. The classically derived a-cell parameters are slightly greater than DFT values, while c-cell parameters are somewhat greater for the DFT case because of lack of dispersion correction in the DMol3 calculations. Normal mode analysis was performed on the resulting structures for comparison with DFT-derived frequencies with the GDIS code.[37] GDIS allows for visualization of eigenvectors for each normal mode. These vibrational modes were classified as bulk (Figure 2A), surface (Figure 2B), or mixed (significant presence of bulk and surface eigenvectors).

Figure 2

Top view of (11̅0) edge surface simulation; edge surfaces are at top and bottom of figure. The normal mode with frequency 718 cm–1 (A) shows primarily “bulk” eigenvectors and the mode with frequency 720 cm–1 (B) shows primarily surface eigenvectors. Black arrows indicate the relative magnitude and directions of the eigenvectors. Mg, O, and H atoms are represented by green, red, and white spheres, respectively.

The LAMMPS code[38] was used to perform all classical molecular dynamics (MD) simulations. Simulations were performed for 1.2 ns in the NPT ensemble at 298 K followed by 540 ps in the NVT ensemble (N, V, T, and P refer to the number of atoms, volume, temperature, and pressure, respectively). While a shorter simulation time would have been sufficient to investigate dry surfaces such as these, the simulation time was extended to monitor the stability of surface hydroxyl groups. In order to incorporate a three-body nonbonded term, a new pair_style interaction (NB3B, or nonbond three-body) was derived for implementation within LAMMPS. The three-body term was used for all Mg–O–H interactions, provided that the Mg–O separation was no greater than 2.8 Å and the O–H separation was no greater than 1.2 Å. Separation cutoffs were determined from the minimum in the corresponding radial distribution functions calculated from an MD simulation of bulk brucite. It should be noted that while the three-body potential is classified as being of nonbonded type, the functional form is discontinuous at the cutoff boundary; therefore, this model is formally not applicable to dissociation of surface hydroxyl groups. Simulations with and without the three-body Mg–O–H term are labeled NB3B and Clayff, respectively.

Classical power spectra were calculated from the velocity autocorrelation function[39] using the TINKER code[40] over the final 40 ps of an MD run (10 000 configurations). A simulation time of 40 ps has been shown to be sufficiently long to obtain vibrational data, particularly given that molecular vibrations have periods much less than 1 ps.[39] Power spectra were used to compare Mg–O–H vibrational modes with and without the new angle-bending potential. All force field parameters are given in Table 2.

Table 2

Force Field Parameters Used in the Classical Simulationsa

element	atom type	q (e)	σ (Å)	ε (kcal mol–1)
Mg	mgh	1.0500	5.2643	9.030 × 10–7
O	oh	–0.9500	3.1655	0.1554
H	ho	0.4250	0.0000	0.0000

σ = (σσ)1/2 and ε = (εε)1/2 represent the van der Waals radius and energy well depth for the atomic pair. r is the interatomic distance. D0 is the depth of the bond potential well. α is related to the width of the bond potential energy curve. r0 is the equilibrium bond distance. k is the angle bend force constant. θ0 is the equilibrium angle.

Ref (13). A potential cutoff of 10 Å was used.

Ref (23).

The DMol3 code was also used to perform ab initio MD (AIMD) calculations in the NVT ensemble for the same three structures used in the DFT optimizations above. The same GGA/PW91 functional and DND basis set were used. The systems were initially equilibrated for 7.5 ps (0.5 fs time step) at 298 K, and data were collected over an additional 2.5 ps production run. The Mg–O–H bond angle distributions from the final 2.5 ps of the equilibration runs were nearly identical to the distributions from the respective productions runs, indicating that equilibration had been reached. Bond angle information for Mg–O–H was compared to the MD-derived distribution of angles for validation of the optimized θ0-parameter.

Results and Discussion

Equilibrium Mg–O–H Angle (θ0)

For bulk brucite, the experimental Mg–O–H angle is 120.2° as determined from neutron diffraction experiments.[25] The bond angle distribution for the DFT energy-optimized (0 K) structure of a 2 × 2 × 2 simulation cell of bulk brucite (Figure 3A, top) is relatively narrow (as expected for a 0 K minimized structure) with a mean angle of approximately 119°. The small deviations from perfect symmetry are a consequence of numerical factors during the DFT optimization, rather than being due to any physically significant distortion mechanism. Yu and Schmidt[23] derived an Mg–O–H term for a specific hydrotalcite edge structure that uses a θ0 value of 132°. This value was determined solely for (11̅0) edge surfaces. The new three-body term presented here is designed to be applied to all Mg–O–H interactions (within specified cutoffs), including those in the bulk, as well as basal and edge surfaces. The bulk Mg–O–H angles are therefore assumed to dominate the overall distribution, even for surface models.

Figure 3

Comparison of Mg–O–H bond angle distributions for brucite from static DFT and classical optimization (panels A and B) and MD and AIMD simulation at 298 K (panels C and D). Angle distributions are shown with (NB3B) and without (Clayff) a nonbonded three-body term for the Mg–O–H interaction. Surface angles (panels B and D) are distinguished from bulk angles (panels A and C).

Figure 3(A,B) shows the equilibrium Mg–O–H angles calculated for the bulk and surface structures from DFT and classical geometry optimizations. Distributions are shown for simulations applying (1) only previous Clayff parameters, including a Morse potential for O–H bond stretch, and (2) the new three-body term in addition to previous Clayff parameters (noted as NB3B). To implement the three-body term, a value of 120° was assumed for θ0 (from the DFT-derived distribution) and a value of 6.35 kcal mol–1 rad–2 was taken for k (from the Mg–O–H term used by Yu and Schmidt[23]). The choice of the force constant value will be discussed below. Surface angles (Figure 3B) are distinguished from internal (bulk) angles (Figure 3A). For the bulk structures, Clayff and NB3B force fields yield an equilibrium angle of 123°, slightly greater than the DFT angle of 119°. For both surfaces, angles determined using NB3B are in better agreement with DFT than those using Clayff, which exhibits a greater range of values.

For the (001) basal surface model (Figure 3A, middle), the bulk angles calculated from NB3B and DFT are nearly at their respective bulk values, while Clayff shows widely varying angles (94°, 116°, 150°). For the (11̅0) edge surface model (Figure 3A, bottom), NB3B and DFT have a similar narrow distribution of angles (121°–124° for NB3B, 115°–121° for DFT) while Clayff values extend well beyond this range (99°–160°). For this model, there is a clear difference between angles for OH groups at the center of the cell and those just below the surface, which are influenced by the relaxation of the surface, so only those angles associated with OH groups at the center of the slab are included in the edge surface bulk angles. Trends in calculated Mg–O–H angles at the surfaces (Figure 3B) are similar to trends observed with bulk angles. NB3B and DFT show narrow angle distributions, while the Clayff distribution is much broader (101°–159°). The new three-body term greatly improves the treatment of Mg–O–H angles at both basal surfaces and edges.

Angle distributions from MD simulations (298 K) of larger simulation cells are plotted in Figure 3(C,D). While thermal effects and larger system sizes result in a more continuous distribution of angles, the same trends seen in the static calculations are reproduced. The bulk angle distributions from Clayff are relatively wide, extending from about 70° to 170°. The three-body term implemented in NB3B has the effect of narrowing the distributions to about 95°–150° and makes them more symmetric. The NB3B results have angle distributions with maxima at about 120°, which reflects the choice of 120° for the θ0-parameter. The AIMD results are in good agreement with the NB3B results, with ranges of nearly the same width and distributions shifted ∼5–10° lower in the AIMD cases.

Mg–O distances were calculated for nearest neighbors in each optimized structure, comparing surface and bulk hydroxyl groups. Clayff and NB3B consistently overestimate the Mg–O distance (2.22–2.23 Å) relative to DFT (2.12–2.13 Å) for bulk Mg–O pairs or at the basal surface. This same trend was also seen for Mg–O distances at inner edges of sepiolite and palygorskite.[41] However, NB3B correctly predicts the decrease in Mg–O distance at the (11̅0) edge surface (2.17 Å), as observed in DFT (2.04 Å).

Mg–O–H Angle Bend Force Constant (k)

The originally published Mg–O–H three-body term in Clayff used a k-parameter of 30.0 kcal mol rad–2 for all metal–O-H surface hydroxyl interactions.[13] However, it was noted that the implementation of the three-body term was problematic, and it was suggested as an optional feature for surface simulations.[13] On the basis of the Mg–O–H parameters developed by Yu and Schmidt[23] for the (11̅0) hydrotalcite edge, we initially implemented a k-parameter of 6.35 kcal mol rad–2. This value was determined by comparing the classical and DFT phonon spectra of the brucite surface. One distinct difference between the Yu and Schmidt model and that of this study is that modified atomic charges were required for all species in the Yu and Schmidt model for hydrotalcite. In contrast, we have kept the charge scheme defined in the original Clayff.

To validate the three-body potential in combination with the original Clayff parameters, the magnitude of the force constant was refined by fitting against the DFT structural and vibrational results, as previously described. Depending on whether bulk brucite or the (001) surface was used as the basis of the fit, the force constants obtained were 5.08 kcal mol–1 rad–2 and 5.81 kcal mol–1 rad–2, respectively. For the case of bulk brucite, a comparison of the high-frequency Mg–O–H angle-bending mode is as follows: the value obtained directly from DFT is 968.6 cm–1, while the original and fitted force field values are 1023 and 1001 cm–1, respectively. Although this suggests that the force constant should be lowered further, it is necessary to recall that the fit is a compromise between multiple vibrational modes and the need to describe the structure correctly. Overall, we can conclude that the fitted values are largely consistent with the force constant from Yu and Schmidt of 6.35 kcal mol–1 rad–2, suggesting that this is a reasonable choice given that the DFT results will also have systematic errors relative to experiment.

Normal modes obtained following DFT optimization were used as a baseline for comparison with GULP-derived classical normal modes using Clayff and NB3B (Figure 4). Modes below a frequency of ≈400 cm–1 represent motion of Mg and O atoms and were not categorized or plotted. Mode frequencies between 400 and 1200 cm–1 correspond to Mg–O–H librational modes, while O–H stretch modes have frequencies of approximately 3600 cm–1. Our DFT frequencies are in good agreement with previous DFT normal mode calculations of bulk brucite using a single unit cell (5 atoms).[42]

Figure 4

Comparison of DFT-derived normal mode frequencies for bulk (top), basal surface (001) (middle), and edge (11̅0) surface (bottom) of brucite. DFT modes are shown in green, Clayff modes in red, and NB3B modes in blue.

Both Clayff and NB3B produce similar bulk bending modes in good agreement with DFT (Figure 4, top). However, DFT surface models are much better represented by NB3B than by Clayff. The highest-frequency Mg–O–H mode at the basal surface (Figure 4, middle) is 592 cm–1 for Clayff, compared to 704 cm–1 for NB3B and 687 cm–1 for DFT. Likewise, the highest frequency at the (11̅0) surface (Figure 4, bottom) in the lower range is 489 cm–1 for Clayff, compared to 735 cm–1 for NB3B and 752 cm–1 for DFT. For reference, frequencies from vibrational spectroscopy of brucite samples show hydroxyl bend modes up to 725 cm–1.[43,44] The improvement of surface mode frequency agreement for NB3B is another indication that inclusion of the new three-body term is an improvement in the overall model.

The effects of the three-body term are also seen by comparing MD power spectra. MD simulations were performed for larger 7 × 5 × 7 supercells of bulk brucite as well as the basal (001) and edge (11̅0) surfaces. Power spectra shown in Figure 5 were decomposed into “bulk” and “surface” contributions based on the disposition of the H atoms. Further refinement of the power spectra by specific hydroxyl groups in the model is not possible because the spectra are obtained from a time average of atomic velocities for specific atom types. As with the normal mode analysis, a distinct shift to lower frequencies is seen for the surface librational modes compared to the bulk modes, whether the three-body term was used or not. For basal surface modes (Figure 5B), dominant peaks in the NB3B power spectrum at ≈475 and ≈700 cm–1 coincide with DFT normal modes (Mg–O–H stretch) between 415–465 cm–1 and 685 cm–1 (Figure 4). The Clayff spectrum shows poor agreement with the DFT surface modes, with only one broad peak at ≈560 cm–1. Within the Mg–O–H librational regime at the (11̅0) edge surface (Figure 5C), the surface power spectrum using Clayff is diminished around 600 cm–1 (red dashed line), while the surface power spectrum for NB3B (blue dashed line) shows well-defined peaks in the 600–750 cm–1 region, in agreement with DFT modes near 750 cm–1 (Figure 4).

Figure 5

Hydrogen vibrational power spectra from classical MD simulation for (A) bulk brucite, (B) the (001) basal surface, and (C) the (11̅0) edge surface.

A similar comparison is made for the librational motion of bulk Mg–O–H groups in the surface models. Slight differences in frequencies for bulk modes in the three model systems are likely due to the influence of the vacuum interface in the surface models. DFT bulk normal modes for Mg–O–H interactions extend to 1013 cm–1. The bulk power spectrum using Clayff (Figure 5A) begins to fall off at about 900 cm–1 (a slight underestimation of the DFT normal modes), while the NB3B spectrum identifies a vibrational response up to ≈1100 cm–1 (a slight overestimation). When a nonbonded three-body interaction potential is included, there is significant improvement in the agreement of the MD power spectrum with the DFT normal modes. Additional tests were performed with other k-parameter values, but these always resulted in worse agreement with DFT; for the (11̅0) edge, the correspondence of MD power spectra peaks at ≈625 and ≈740 cm–1 with observed normal modes (Figure 4) is particularly convincing. The value of 6.35 kcal mol rad–2 derived by Yu and Schmidt[23] was therefore maintained for additional calculations.

The O–H bond stretch region for brucite is known to exist up to ≈3700 cm–1.[45,46] In the original implementation of Clayff, a harmonic bond was used to describe the O–H stretch, but a Morse potential has been recently introduced.[24] Here, the Morse potential was used for the O–H bond stretch in both Clayff and NB3B, and their resulting power spectra are shown in Figure 6. For both surfaces, the surface O–H bond stretch frequency is higher than the bulk (internal) frequency, which may be attributed to hydrogen bonding between layers in the bulk that results in weaker hydroxyl bonds. For O–H vibrations in the bulk or at the surfaces, the NB3B power spectrum shows a shift toward higher frequencies than the Clayff peaks, in better agreement with the corresponding DFT modes (Figure 4).

Figure 6

Hydrogen vibrational power spectra from classical MD simulation for (A) bulk brucite, (B) the (001) basal surface, and (C) the (11̅0) edge surface. Spectra for both the bulklike and surface regions are shown separately in panels B and C.

Conclusion

We have introduced a nonbonded three-body harmonic term in classical molecular simulations for describing Mg–O–H dynamics while avoiding the tendency of bonded hydroxyl groups to “dissociate” from the surface. Although validated for hydroxyl behavior in bulk brucite and at two surfaces, this methodology can be applied to similar bulk or surface models of layered minerals using the Clayff force field. Detailed comparisons of surface vibrational modes using DFT normal modes as a reference reveal that the addition of an Mg–O–H bending term substantially improves the ability of Clayff to describe librational motion at hydroxylated surfaces. The new Mg–O–H term was also applied to bulk brucite. In this case, the Mg–O–H librational modes predicted by Clayff already show decent agreement with DFT, so a substantial improvement in Clayff was not realized. Direct fitting of the vibrational modes obtained from DFT for the new Mg–O–H term supports the k-value proposed previously[23] for modeling hydrotalcite and establishes a “proof-of-concept” procedure that can be applied to other metal–O-H interactions of interest for layered mineral simulations, especially those involving clay edges and other interfacial interactions associated with environmental processes.