Molybdenum disulfide and water interaction parameters

Water in nanopores and nanochannels has unique properties that differ from those of bulk water.^1,2 These unique properties of confined water can be used to exploit novel functionalities and applications. For example, water in carbon-based nano confinements has been found to possess fast translational motion,^3–7 rotation-translation coupling,⁸ and fast rotational motion^9,10 leading to potential applications such as seawater desalination^11,12 and energy storage.¹³ Compared to nanoporous graphene and carbon nanotubes, nanopores in molybdenum disulfide (MoS₂) membranes have shown enhanced biological compatibility and higher signal-to-noise ratio for biological sensing applications.^14–16 In addition, nanoporous MoS₂ membranes have been found to be excellent candidates for water desalination¹⁷ and electric power generation.¹⁸ As these potential applications involve MoS₂ in water environments, it is important to develop accurate MoS₂ force fields for nanofluidic applications.

Molecular dynamics (MD) simulations have been used extensively to study the performance of MoS₂ membranes in DNA sequencing,¹⁵ water desalination,¹⁷ osmotic power generation¹⁸ and for a variety of other applications. However, the accuracy of MoS₂/water interfacial properties predicted by MD simulations depends on the interatomic potential employed in the simulations. Therefore, developing accurate force field parameters can result in a more quantitative estimation of MoS₂/water interfacial properties. Currently, high-accuracy force fields describing the non-bonded interaction between MoS₂ and water molecules are lacking. The force field parameters are generally developed using two approaches. First, the van der Waals (vdW) parameters are developed using the combination rule, known as the Lorentz-Berthelot¹⁹ rule, but the accuracy of this approach is uncertain and it can lead to erroneous results.^20,21 The force field parameters (both vdW and Coulombic) can also be developed from the total interaction energy obtained from first-principles calculations where the force field interaction potential functional form is fitted to the calculated interaction energy curve.^22,23 In this approach, the accuracy of the force field parameters directly follows the fidelity of the first-principles calculations.

High-accuracy ab initio calculations, which can describe electron-electron correlations explicitly, are needed as the vdW interactions are affected by the correlation energy. The random phase approximation method (RPA)^24,25 treats electron-electron correlation explicitly. Prior research has shown that the parameters developed entirely from the first-principles calculations with no adjustable parameters are able to predict properties that are comparable to those of experiments.^22,23,26,27

In this article, we use the RPA method to compute the potential energy surface between MoS₂ and water. The force field parameters are obtained from the potential energy surface and no fitting to experimental data is used in this step. Finally, the solid-liquid (MoS₂ sheet and water) work of adhesion, which is calculated by the dry-surface method^28,29 using the developed parameters, is compared with the work of adhesion values obtained from experimentally measured contact angles to validate the proposed force field parameters.

The RPA method, which is used to treat electron-electron correlations, was first introduced by Bohm and Pines.^30–32 Later, Langreth and Perdew showed that the RPA method could be formulated using the adiabatic-connection method (AC)^33–35 along with the fluctuation-dissipation theorem (FD)^36,37 within the density functional theory (DFT).³⁸ This approach, known as ACFD RPA method, has been one of the most efficient RPA computation methods during the last decade.^39–41 The RPA correlation energy is written as

where χ⁰ is the response function of the Kohn-Sham non-interacting system, V is the interacting potential operator, and ω is the frequency (see the  Appendix for more information about the RPA method).

The VASP package^42,43 was used to carry out the RPA calculations. The projector-augmented wave method was applied for more efficient DFT calculations. The energy cutoff was 408 eV. An energy cutoff of up to 272 eV was used to expand the density response function in plane waves. Different MoS₂ supercell sizes and Brillouin zone sampling Monkhorst-Pack k-meshes were used for the exchange and correlation energy calculations. An 8 × 8 supercell with a 2 × 2 k-mesh and a 2 × 2 supercell with a 4 × 4 k-mesh were used for the exchange and correlation energy, respectively. The lattice constant along the direction normal to the MoS₂ surface (z-axis) was 15 Å. The RPA method is used to compute the potential energy surface between a MoS₂ cluster and a water molecule for different separation distances between the two as well as different orientations of the water molecule. The geometry of the MoS₂ and water molecule was kept rigid for simplicity. The structure of the MoS₂ cluster follows the geometry obtained from molecular statics calculations,⁴⁴ with S–S distance of 3.241 Å between two layers of sulfur atoms and a bond length of 3.17 Å, 2.42 Å, and 3.17 Å for Mo–Mo, Mo–S, and S–S, respectively.

The MD simulations have been performed using the LAMMPS package⁴⁵ to obtain the water–MoS₂ work of adhesion using the force field parameters developed in this work. Initially, a cubic box of water molecules is placed on top of a single-layer MoS₂ sheet. The water box contains 2400 water molecules. The dimension of MoS₂ sheet is about 10 nm × 10 nm to ensure that the water molecules do not interact with their periodic images. The entire simulation box has a size of 10 nm × 10 nm × 20 nm in x, y, and z, respectively (the MoS₂ sheet lies in the xy plane). Both the water box and MoS₂ sheet are created using the visual molecular dynamics.⁴⁶ First, the energy of the system is minimized for 10 000 steps. Next, the system is equilibrated in the NVT ensemble for 4 ns at a pressure of 1 atm and a temperature of 300 K. The constant temperature is enforced using the Nosè-Hoover thermostat.^47,48 During the simulations, the MoS₂ atoms are kept frozen to investigate the water–MoS₂ interface properties without the effect of wall motion. A time step of 1.0 fs is used for MD integration. The short-range vdW interactions are calculated using a cutoff scheme with a cutoff distance of 1.4 nm. The long-range electrostatic interactions are calculated using the Particle-Particle-Particle-Mesh (PPPM) method.⁴⁹ During the equilibration simulation, the cubic water box spreads over the MoS₂ sheet where the energy and temperature of the system stabilize and reach constant values. The resulting configuration is used as the starting point for production simulations that are run for further 4 ns to compute the work of adhesion.

Before computing the potential energy surface between MoS₂ and water, the accuracy of the calculations needs to be examined. In addition to the approximations made in the RPA mathematical formulation, other sources of error arise from the computational setup. The effect of all the sources (the incomplete basis sets, the finite lattice constant, the finite supercell, the finite k-points, and the non-zero water coverage) on the interaction energy is considered to ensure the accuracy of the calculations. First, the effect of the basis set is investigated by increasing the plane wave energy cutoff. This results in a change of about 1 meV in the energy when the energy cutoff is increased to 600 eV and 400 eV for the plane wave and the response function, respectively. Next, the lattice constant along the z-axis is increased from 15 to 20 Å in which a difference of less than 1 meV is observed in the interaction energy. The effect of k-points is studied by changing the k-grids for the exchange and correlation energies separately (see Fig. 1). As shown, the correlation energy changes with increasing the k-grid for a fixed 2 × 2 supercell while the exchange energy is weakly dependent on the k-grid. The correlation energy starts to converge for a 3 × 3 k-grid and a 2 × 2 supercell where increasing the k-grid to 4 × 4 results in a change of about 2 meV in the energy. Finally, to test the effects of the supercell size and water coverage on the energy, the supercell size is changed for the exchange energy which is shown to have a strong dependence on the size of the supercell. The exchange energy starts to converge for a 4 × 4 supercell and a 2 × 2 k-grid, while increasing the supercell size to 8 × 8 only improves the energy by about 2 meV. Based on the error analyses, the maximum accumulated error in the interaction energy by the RPA method is estimated to be less than 4 meV which is insignificant compared to the exchange and correlation energies.

With the errors controlled, using the RPA method, the potential energy surface between MoS₂ and water is computed to implement a force field for use in MD simulations. The total interaction energy between a water molecule and a MoS₂ cluster is computed by

where ΔE is the total interaction energy, $EMoS2−water$ is the energy of the MoS₂–water system together, and $EMoS2$ and E_water are the energies of the MoS₂ and water molecule, respectively. These three energies are calculated separately using RPA. The partial charges for atoms of the water molecule are based on the SPC/E water model.⁵⁰ To calculate the partial charges of Mo and S atoms, the electrostatic potential of the MoS₂ (no water molecules present) is first obtained from the 2nd-order Møller-Plesset perturbation (MP2) theory calculations.^23,40 Then the partial charges are obtained by fitting to the calculated electrostatic potential.

The model used in MD simulations to describe the non-bonded interactions between the atoms of MoS₂ and water consists of two components, namely, the vdW interactions (ΔE_vdW) and the electrostatic interactions (ΔE_elec). The vdW interactions are characterized by the Lennard-Jones (LJ) potential parameters. The 12-6 LJ potential form is defined by

where ε_ij is the depth of the potential well between any two atoms (i and j atoms of MoS₂ and water, respectively), σ_ij is the distance at which the potential is zero and r_ij is the distance between atoms i and j. The electrostatic interactions are modeled by Coulomb’s law

where q_i and q_j are the partial point charges of atom i and atom j. The total interaction between the MoS₂ and water is then obtained from

For the vdW interactions, ε_ij and σ_ij are obtained by fitting to the RPA energies where the electrostatic component is excluded from the total energy. The water molecule has two types of atoms (O and H atoms) and MoS₂ has Mo and S atoms; therefore, four pairs of ε_ij and σ_ij parameters are required for fitting. However, for simplicity, the H–MoS₂ vdW interactions can be excluded from the parameterization fitting given the O–MoS₂ parameters can be corrected to match the potential energy from the RPA data. This simplification is valid if the vdW interaction due to O dominates compared to that of H. This is tested by plotting the vdW component energies from the RPA method for five different water orientations. As shown in Fig. 2, the energies for the five orientations are close enough to simplify the parametrization and rewrite the vdW interactions as

The Boltzmann averaged vdW interaction energies for the five different water orientations (Fig. 2) are fitted, using the least-squares method, to the potential function in Eq. (6) to obtain the vdW parameters. The vdW parameters from the fitting are tabulated in Table I. Also included in the table are the parameters obtained from the combination rule used in the literature. In addition, the water models and partial charges on Mo and S atoms used in each work are summarized in the table. In order to validate the accuracy of the proposed parameters, MD simulations can be used to predict properties of water–MoS₂ interface which are then compared with those of experimental measurements.

The contact angle of water on MoS₂ surfaces is one of the water–MoS₂ properties reported by several experimental measurements.^51–56 Most of the experimental contact angles lie in the range of 65°–75°. The variations in these experimental values are possibly due to the susceptibility of MoS₂ to surface contaminations.^51,52 However, the droplet sizes, simulated in MD, are of nanometers, which are much smaller than the macroscopic droplets used in experiments. The contact angle values, θ, from simulations, which are extrapolated using the modified Young’s equation, are a function of the liquid-vapor surface tension of the water droplet, γ_LV. But this may not be an accurate extrapolation since γ_LV itself changes with the droplet size at the nanoscale. In addition, water models, in MD simulations, result in γ_LV values that are different from the experimental value. It has been shown that these differences can significantly affect the θ estimation in MD leading to possibly unreliable comparison with the experimental contact angles.²⁹ Alternatively, the solid-liquid work of adhesion, W_SL, is a recommended interface property to validate force field parameters against the experimental data as it is shown to be independent of water models and can be calculated explicitly in MD with no dependence on γ_LV and θ.^28,29

The work of adhesion is obtained explicitly using the dry-surface in MD simulations.²⁸ The work of adhesion is the work per unit area of the solid-liquid interface required to isolate the liquid droplet from the solid surface such that there is no interaction between the liquid and solid. In the dry-surface method, this is achieved by varying the interaction parameters (ε and q) from the actual interactive state (A) to a state (B) in which the interaction energy per unit area between the surface and liquid is insignificant relative to the work of adhesion. The work of adhesion for the system of a water droplet and MoS₂ surface is obtained from the following integrations:

where a is the surface area of the MoS₂ membrane, and $ULJn−m$ and $Uelecn−m$ are the interaction energies between atoms of type n and m due to the LJ and Coulombic interactions, respectively. For the force fields with no partial charges, the last two integrals are ignored. The differentiation and integration are performed numerically using the second-order accurate finite difference method and trapezoidal rule. For experiments, W_SL is obtained from the Young-Dupré equation

where the experimental values of γ_LV and θ are used to calculate the work of adhesion. The work of adhesion values for different force fields in MD simulations as well as the experimental values as a function of the binding energy are shown in Fig. 3. As shown, the parameters developed in this work, using the RPA energies, lead to a work of adhesion value which is closest to the range of work of adhesion values obtained from the experimental contact angle measurements. Note that no fitting from experimental data has been applied to develop the force field. The agreement with experimental data shows that using ab initio calculations directly is a reasonably rigorous method to develop accurate force field parameters without any fit to experiments. This has been shown consistently in our other force field developments for graphitic-carbon-water interaction²² and hexagonal boron-nitride-water interactions.^23,40

In summary, we use the ACDF RPA method to calculate the interaction energy between MoS₂ and water. The RPA results, without any fitting to the experimental data, are used to develop force field parameters that are subsequently used in MD simulations. The corresponding water–MoS₂ work of adhesion in MD simulations is found to be comparable to the work of adhesion obtained from the experimental contact angle measurements. The good agreement with experiments shows that developing force field parameters directly from first-principle calculations is an accurate method to predict properties of emerging systems at molecular scales.

This work is supported by AFOSR under Grant No. FA9550-12-1-0464 and by NSF under Grant Nos. 1264282, 1420882, 1506619, and 1545907. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) which is supported by National Science Foundation Grant No. OCI-1053575. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Award Nos. OCI0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.

For a given electronic system, the total energy (E) of the ground state, in the Kohn-Sham (KS) DFT, is written as

where T_KS is the kinetic energy of the KS non-interacting system, E_ion-e is the ion-electron interaction energy, E_H is the Hartree energy, and E_xc is the electron-electron interaction energy known as the exchange-correlation energy. The GGA-PBE functional⁵⁷ is used to obtain the electron density based on the single-particle KS orbitals. All the terms in Eq. (A1) are functionals of the electron density, ρ(r), which, in the KS DFT, is obtained by

where N is the number of particles and $ϕir$ is the KS orbital for a single particle i. In this approach, the difficulty arises in computing the electron-electron interaction term (E_xc). To do that, the AC formulation is first used to write the ground state total energy of the system. The Hamiltonian of the system (⁠$Ĥ$⁠) can be written as a summation of non-perturbative (λ = 0) and perturbative (λ) terms as follows:

where λ is a dimensionless parameter (0 ≤ λ ≤ 1). Considering the Hamiltonian for an electronic system and that $HΨλ=EΨλ$ (⁠$Ψλ$ is the wave function of the ground state), the total energy in the ground state is obtained by

In KS DFT, the electron density is kept fixed and the zeroth order ground state wave function, $Ψ0$⁠, is simply the Slater determinant of the single-electron orbitals. Given that the electron-electron interaction operator is defined as $V=∑i,j≠iNλri−rj$ (r_i is position vector of electron i), by equating equations (A1) and (A4), the exchange-correlation energy is written as

where $ρ^r=∑i=1Nδr−ri$ is the electron-density operator, $δρ^r=ρ^r−ρr$ is the density fluctuation operator, and ρ(r) is the expectation value of the electron-density operator.

The zero-temperature fluctuation-dissipation theorem is used to express the expectation value of the two-particle density fluctuation operator in Eq. (A5) as a frequency integration of the dissipation of the system

where $χλ(ω,r,r′)$ is the linear density-density response function of the interacting system which can be expressed in terms of the response function of the KS non-interacting system, $χ0(ω,r,r′)$⁠, as follows:³⁹ 

where V is the interacting potential operator and $fxcλ$ is the exchange-correlation kernel. $χ0(ω,r,r′)$ is simply obtained by the knowledge of the single-particle KS orbitals. In order to perform the integration in Eq. (A6) practically, the response function, given in Eq. (A7), needs to be approximated by neglecting $fxcλ$⁠. This approximation is known as the random phase approximation [we write the approximated response function as $χRPAλ(ω,r,r′)$]. Using Eqs. (A5)–(A7), the exchange-correlation energy (E_xc) is given by

where the first and second terms in {} are the exchange energy (E_x) and RPA correlation energy (⁠$EcRPA$⁠), respectively. The RPA correlation energy can be written in the following reduced form: