Coarse-Grained Simulation of the Adsorption of Water on Au(111) Surfaces Using a Modified Stillinger-Weber Potential.

For reproducing the behavior of water molecules adsorbed on gold surfaces in terms of density of both bulk and interfacial water and in terms of structuring of water on top of gold atoms, the implementation of a multibody potential is necessary, thus the Stillinger-Weber potential was tested. The goal is using a single nonbonded potential for coarse-grained models, without the usage of explicit charges. In order to modify the angular part of the Stillinger-Weber potential from a single cosine to a piecewise function accounting for multiple equilibrium angles, employed for Au-Au-Au and Au-Au-water triplets, it is necessary to create a version of the simulation package LAMMPS that supports the assignment of multiple favored angles. This novel approach is able to reproduce the data obtained using quantum mechanical calculations and density profiles of both bulk and adsorbed water molecules obtained using classical polarizable force fields.

Introduction

In the last few years, in nanoscience, the use of functionalized nanoparticles with biomolecules has gained a lot of interest, with applications in drug delivery,[1] theranostics,[2] heterogeneous catalysis,[3] and clinical[4] and environmental[5] analyses. Functionalization is possible by using covalently bonded molecules on Au surfaces by means of thiol groups. Noncovalent passivation is also used in water media to promote ligand exchange and to control the nucleation-growth of nanoparticles.[6] With the aim of studying these systems with classical molecular dynamics simulations, traditional force fields using exclusively two-body potentials for the nonbonded interactions are no longer effective. An important example is the simple adsorption of water on the Au(111) surface. Density functional theory (DFT) calculations state that a single water molecule has to be on top, forming an angle of 90° with two Au surface atoms and that hydrogen atoms do not have to point directly to the Au surface.[7] These results cannot be reproduced by using only two body terms, even if we are talking about water models at full atomistic resolution (including electrostatics). In fact, calculations performed using a nonbonded attractive-repulsive two-body term for the interaction of a single water molecule with a face centered cubic (FCC) Au(111) surface will give rise to a large free-energy minimum on the “hollow” site (because of the contribution from interactions with three surface atoms).

To overcome this problem, a multibody potential becomes necessary and the three-body Stillinger–Weber Potential (SWP) can be successfully used. SWP parameters for water molecules in bulk,[8] polyethylene,[9] and organic solutes[10] in water are already present in the literature. In the latter case, the authors were able to reproduce the potential of mean force (PMF) of a methane dimer in water, establishing a basis for a CG force-field for aqueous soft matter systems based only on nonbonded interaction without electrostatic contributions. Other examples of efficient usage of the SWP are the original works on silicon,[11] gallium,[12] tetrahedral carbon,[13] and other elements in a cubic lattice.[14] Furthermore, other CG applications with SWP exist for calcite/water nucleation-growth mechanisms.[15] An implementation of three-body interactions, coupled with a polarizability term in the two-body function, has been successfully used in classical molecular dynamics simulations, in full atomistic representation, for the adsorption of water on Au surfaces.[16]

Full atomistic models using electrostatic interactions alongside a fully parametrized three body term for the physisorption or chemisorption of water on noble metal 111 surfaces are reported by Steinmann et al.[17] and Clabaut et al.[18]

As far as we know, a simulation model for the adsorption of water on an Au surface using a CG approach and without any explicit polarizability and/or electrostatic term does not exist so far.

In this article, we searched for an accurate description of the interaction between water molecules represented by the Molinero water (mW) model[8] and an FCC Au(111) surface, parameterizing Au and water two- and three-body interactions in a modified version of the SWP for CG molecular dynamics simulations.

The SWP functional form (eq ) is explicitly separable in two parts: an attractive-repulsive term for the two-body interaction (eq ) and a purely repulsive three-body angular dependence (eq )Here, A and λ are the strengths of the two-body and three-body interactions and alongside B, σ, a, cos θ0, γ, and ε represent the full set of parameters to completely describe the SWP:[11] ε = ε = ε according to the LAMMPS definition.

In particular, for the two-body function (φ2), it is possible to optimize the parameters in order to reproduce the force constant, energy minimum, and equilibrium distance, matching the zeroth, first, and second derivatives of the potential (eq ), calculated at the equilibrium value (r = re).With this approach, the large number of parameters in the SWP can be reduced to only two quantities, namely re/σ and a.

Methods

All the simulations have been performed using the simulation package LAMMPS[19] with a 5 fs time step in the isothermal–isobaric (NPT) ensemble with a semi-isotropic pressure coupling. Temperature and pressure were controlled using a Nosè–Hoover[20] thermostat and a barostat using 100 fs, and 1 ps dumping constants, for temperature and pressure, respectively. The equilibrium external pressures were 0, 0, and 1 atm for X, Y, and Z, respectively, where Z is the axis normal to the Au surface. The simulations were performed by using a box (5 × 5 × 9 nm3) containing a slab of 2880 Au atoms (FCC) and 5785 water molecules and applying the periodic boundary conditions.

A piecewise function for the angular dependence of the SWP has been implemented in C++ as a new pair style (link for LAMMPS downloadable input files and for the modified code: https://github.com/gripani/LAMMPSinput_file and https://github.com/gripani/LAMMPSbuild_SWmod). Analysis of the simulations has been performed by means of homemade programs.

Results and Discussion

The parameters for the water-Au pair have been refined to reproduce the a posteriori feature (mass density) of the GolP-CHARMM polarizable force field.[21]

The parameters for the Au bulk concerning the two-body part only have been obtained starting from the cohesive energy (EC = ε), lattice parameter (), and bulk modulus ()already existing in the literature.[22]

The three-body function (φ3) depends on λ, γ, and cos θ0. For the gamma value, we have used the already parametrized value, γ = 1.2[11] for all the possible triplets in the system. According to Molinero, the tetrahedrality of bulk water is provided by cos θ0 = −1/3 and λ = 23.15.[8] Concerning the Au FCC lattice in the bulk phase, it is not possible to reproduce the three-body energy using a single value of the equilibrium angle.[14] The same issue, with more critical aspects, appears in the case of the mW–Au–Au triplet as well.

Indeed, one adsorbed water and its nearest neighbor Au atom, for the FCC lattice, form an angle of 90° with an Au atom of the first layer and an angle of 144.7° with a nearest neighbor Au of the second layer (Figure ).

Figure 1

Schematic representation of the two different angles formed by the normal to the Au surface and two adjacent Au atoms, connected by red lines (corresponding to the Au–Au equilibrium distance). The 90° angle is shown in light green and the 144.7° angle is shown in violet.

The functional form of the piecewise function (eq ) concerning N different equilibrium angles is given as follows (xmin and xmax are the limiting values of the n-th piece interval)

In Table , the parameters of the piecewise function for an FCC lattice are reported. The same g function with modified parameters (reported in Table ) can also be used to describe the behavior of water molecules interacting with the Au atoms, assuring the preferential top site geometry.

Table 2

Parameters for the Piecewise Function Associated with the Au–Au–Au Triplet (FCC)

n	xnmin	xnmax	cos θne	gn0
1	–1	–7/8	–1	0
2	–7/8	–5/8	–3/4	1/32
3	–5/8	–3/8	–1/2	0
4	–3/8	–1/8	–1/4	1/32
5	–1/8	1/8	0	0
6	1/8	3/8	1/4	1/32
7	3/8	1	1/2	0

Table 3

Parameters for the Piecewise Function Associated with the mW–Au–Au Tripleta

N	xnmin	xnmax	cos θne	gn0
1	–1	–3/2/2	–2/3	0
2	–3/2/2	–1/26	–1/6	1/12
3	–1/26	1	0	0

The value g0 ensures continuity and derivability to the function that has to be evaluated numerically.

In order to implement the modified angular dependence for the multiangle approach, we defined a new pairstyle in the LAMMPS source code (SWPmod).

Specifically, the pairwise function g(cos θ) (eq ) goes in place of (cos θ – cos θ0)2 in the three-body function (eq ).

Because of technical reasons, we have introduced in our code a switch parameter, with the role of choosing which piecewise function is necessary for which triplet (alongside the full set of parameters n, xmin, xmax, cos θe, and g0). In this way, the approach can be easily extended to other molecular systems. In particular, if switch = 0, no piecewise function will be used (single-angle harmonic cosine function); if switch = 1, the two angle piecewise function will be used for the mW, Au, Au triplet, as reported in Table ; if switch = 2 the four angle piecewise function will be used (for bulk Au as reported in Table ).

Table 1

Optimized Two-Body SWP Parameters for all the Possible Pairs in the System

I	j	re (Å)	σ (Å)	a	A	B
mW	mW	2.685	2.393	1.800	7.04955628	0.60222
Au	Au	2.951	2.045	2.053	1.76300000	2.223000
mW	Au	3.000	2.250	2.150	5.10293601	1.053210

It is noteworthy that it is possible to realize a fine tuning of the correct top site configuration by changing the strength of the repulsive three-body term (λ) for the mW–Au–Au triplet. Figure (top) shows the average orientation of water molecules with respect to the Au surface sites. The figure nicely illustrates how it is possible to shift the equilibrium configuration of the adsorbed water molecules from the hollow to the top sites, performing NPT simulations with different λ values (Table ).

Figure 2

(Top) Water densities at different lambda values for the mW, Au, Au triplet (10, 20, 30, 40, and 50, from left to right). The polar plots report the water density (see color code on the right) as a function of polar and azimuthal angles referred to the vector connecting a single Au atom to the first neighbor water molecule, averaged on all the Au atoms of the surface. The water adsorption preferences of the angles on the azimuthal axes refer to about 20° for hollow sites and 0° and about 40° for top sites (central and next Au atoms). (Bottom) Time-averaged Z density of water (g/mL) of NPT simulations (in red) and GolP reference (in green) vs the distance (Å) of water molecules from the Au surface.

Table 4

Three-Body SWPmod Parameters Used in all the Simulations

i	j	k	ε (kcal/mol)	λ	cos θ0	γ	switch
mW	mW	mW	6.189	23.15	–1/3	1.2	0
Au	Au	mW	2.000	10–50		1.2	1
Au	Au	Au	73.94	10.00		1.2	2

From the water density in the z-direction (i.e., normal to the surface), the PMF of the adsorption of the water molecules on the Au surface can be obtained (eq ).

The two minima correspond to the first and second layers of the adsorbed molecules, from which it is possible to estimate both equilibrium (eq ) and kinetic (eq ) constants of the exchange process between the two layers (corresponding to the first and second minima of Figure in the case λ = 50).ΔW0 is the difference between the values of the PMF of the first and the second minimum and ΔW† is the difference between the value of the PMF calculated at the first maximum and the value of the PMF calculated at the first minimum. In the Eyring equation (eq ), the transmission coefficient is set equal to 1, as it is usually assumed for transport processes.

Figure 3

PMW obtained by the density profile (λ = 50 in red) and GolP reference (in green). The two minima correspond to the first (3.0 Å) and second (4.5 Å) layers of adsorbed mW.

Conclusions

In conclusion, with the support of a multibody potential, we were able to reproduce the correct configuration of the adsorbed water molecules on the Au(111) surface and the water density profile as obtained using the GolP-CHARMM force field. The model is also able to describe the first and second layers of the adsorbed water molecules. The advantages of this approach are the absence of explicit partial charges so that it does not need any electrostatic contribution on the energy and the reduced number of simulated particles with the CG approach. In contrast, the GolP-CHARMM force field makes use of explicit dummy charges bound to Au atoms and is fully atomistic. Finally, the versatility and general form of the three-body part of the Stillinger–Weber potential allows a possible extension of the proposed approach to other molecular systems.