Quantitative Modeling of MoS2–Solvent Interfaces: Predicting Contact

Article pubs.acs.org/JPCC

Quantitative Modeling of MoS2−Solvent Interfaces: Predicting Contact Angles and Exfoliation Performance using Molecular Dynamics Vishnu Sresht,†,¶ Ananth Govind Rajan,†,¶ Emilie Bordes,†,‡ Michael S. Strano,† Agílio A.H. Pádua,*,†,‡ and Daniel Blankschtein*,† †

Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States Institut de Chimie de Clermont-Ferrand, Université Blaise Pascal and CNRS, 63178 Aubiére, France

‡

S Supporting Information *

ABSTRACT: The large-scale synthesis of molybdenum disulfide (MoS2) using liquid-phase exfoliation, as well as several of its intended applications, including desalination membranes and biosensors, involve liquids coming into intimate contact with MoS2 surfaces. Molecular dynamics (MD) simulations offer a robust methodology to investigate nanomaterial/liquid interactions involving weak van der Waals forces. However, MD force fields for MoS2 currently available in the literature incorrectly describe not only the cohesive interactions between layers of MoS2 but also the adhesive interactions of MoS2 with liquids such as water. Here, we develop a set of force-field parameters that reproduce the properties of bulk 2H-MoS2, with special attention to the distinction between the covalent, intralayer terms and the noncovalent, interlayer Coulombic and van der Waals interactions. The resulting force field is compatible with MD force fields for liquids and can correctly describe interactions at MoS2−liquid interfaces, yielding excellent agreement with experimental contact angles for water (a polar solvent) and diiodomethane (a nonpolar solvent). In light of these results, previously published simulations studies on the desalination potential and biocompatibility of MoS2 devices need to be reevaluated. Potential of mean force (PMF) calculations demonstrate that use of our new force field can explain the current selection of solvents used in experimental studies of the liquid-phase exfoliation of MoS2 flakes, including the colloidal stability of the resulting dispersion. Our new force field enables an accurate description of MoS2 interfaces and will hopefully pave the way for simulation-aided design in applications including membranes, microfluidic devices, and sensors.

S

applications is hindered by the lack of adequate models for the interfacial interactions between MoS2 and liquids. Molecular dynamics (MD) simulations offer a convenient computational framework to study these interactions. Nevertheless, as highlighted in a recent, comprehensive review,19 the atomistic models for bulk and single-layer MoS2 currently available in the literature are not suitable for the calculation of Coulombic and van der Waals interactions with other materials, liquids, or gases. The four most cited articles20 with respect to MoS2 force fields applicable to MoS2 fluid simulations are by (1) Becker et al.21 (BE), (2) Morita et al.22 (MO), (3) Varshney et al.23 (VA), and (4) Liang et al.24 (LI). It is noteworthy that the Coulombic and van der Waals interaction parameters in these four force fields have been combined with various intralayer bonded interaction potentials to yield other force-field models.25,26 The BE force field includes fixed partial charges

ingle-layered and few-layered nanosheets of molybdenum disulfide (MoS2) have exceptional optoelectronic,1,2 structural,3 and barrier4 properties that promise revolutionary advances in membrane5 and sensor6 technologies. The rational design of MoS2-based devices is predicated on a fundamental understanding of the interactions between MoS2 surfaces and the surrounding substrates and liquids. For example, compared to chemical vapor deposition,7 which is only possible on certain substrates, liquid-phase exfoliation enables the large-scale synthesis of the 2D forms of layered materials while also permitting chemical functionalization, deposition on arbitrary substrates, and circuit printing.8−11 However, the absence of an accurate model for 2D material−exfoliation media interactions has precluded a systematic molecular-level screening of exfoliation media, forcing experimentalists to try a variety of organic solvents,12 surfactants,13 and even protein14 solutions, to exfoliate and disperse MoS2 layers. Once synthesized, MoS2 monolayers can serve as excellent microfluidic materials,15 electrodes,16 membranes,17 and ionic-liquid-gated field-effect transistors.18 However, here too, the optimization of these © 2017 American Chemical Society

Received: January 16, 2017 Revised: March 5, 2017 Published: March 6, 2017 9022

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

flux through MoS2 nanopores could be underpredicted or overpredicted, depending on whether the MoS2 force field is overly hydrophilic (VA) or overly hydrophobic (LI), respectively. With the above in mind, in this work, we develop a new force field for MoS2 which overcomes the limitations of the currently available force fields. To this end, we parametrize, from scratch, a transferrable atomistic force field for MoS2 that incorporates the decoupling between the intralayer and the interlayer interactions. Our new force field for MoS2 consists of partial charges and 6−12 LJ parameters for S and Mo, as well as of harmonic potentials for bonds and angles (Figure 1a). In contrast to the efforts by Zhou, Leroy, and co-workers,29,34 we use contact-angle simulations (with water and diiodomethane) not to parametrize our force field, but instead to validate it− thereby demonstrating (a) the ability of the force field to simulate both mechanical and interfacial systems, and (b) its transferability to nonaqueous systems.

in addition to a van der Waals potential and was developed to reproduce the growth of metal islands on molybdenite surfaces.21 The MO force field was parametrized by Morita et al. to model the tribology of MoS2 films.22 The VA force field borrows partial charges from the MO force field and slightly modifies the S−S nonbonded interaction in order to reproduce the crystal structure and thermal properties of 2H MoS2 sheets.23 However, none of these three force fields correctly accounts for the van der Waals interactions. While the MO force field does not include Mo−Mo van der Waals interactions, in the BE and VA force fields, the Mo−Mo nonbonded interactions are extremely strongmore than an order of magnitude larger than the corresponding S−S interactions, in spite of the similar atomic sizes of Mo and S. In the BE, MO, and VA force fields, the artificially larger nonbonded interactions compensate for the neglected bonded interactions, resulting in accurate predictions of the experimentally measured intralayer mechanical properties of MoS2 sheets but in inaccurate predictions of the van der Waals forces between MoS2 sheets. Consequently, the nonbonded interactions between MoS2 sheets modeled using the BE, MO, and VA force fields cannot be used with combining rules to derive the interactions of MoS2 with other materials. Moreover, as demonstrated below, the BE and VA force fields significantly overestimate the cohesive intersheet interactions and also yield predicted water contact angles on MoS2 that are not consistent with those determined experimentally.27,28 Luan and Zhou29 attempted to refine the parameters of the VA force field by fitting Lennard-Jones (LJ) parameters to reproduce the water contact angle on MoS2. However, they used the experimental value of the water contact angle on an MoS2 surface that was later shown to have been contaminated by exposure to air.27,28,30 Liang et al. developed the LI force field by incorporating a reactive bond-order (REBO) potential for intralayer interactions in addition to a LJ potential for the interactions between MoS2 layers. However, the LI force field does not have Coulombic charges and, therefore, cannot describe Coulombic interactions between MoS2 and polar molecules. Furthermore, the strength of the Mo−Mo nonbonded interactions in the LI force field is 2 orders of magnitude smaller than that of the corresponding S−S interactions, which is once again inconsistent with the sizes of the two atoms. The discussion above indicates that the available force fields for MoS2 are not reliable to simulate the interactions of MoS2 with other materials or fluids, because they do not distinguish between the roles of intrasheet bonded terms and intersheet van der Waals and Coulombic terms. The mechanical rigidity of the MoS2 layers should result from the intrasheet bonded terms, whereas the interlayer interactions should result from the intersheet van der Waals and Coulombic terms. In spite of these limitations, the parameters in the BE, MO, VA, and LI force fields have been utilized to make quantitative predictions of the following: (i) the friction behavior and wettability of MoS2 surfaces,29 (ii) the interactions of proteins with MoS2 nanotubes and sheets,31,32 and (iii) the desalination of water using nanoporous MoS2.5,33 In the case of (i), an underestimation of the interactions between water molecules and the MoS2 substrate will lead to an inaccurate estimate of the friction coefficient. In the case of (ii), the high hydrophobicity of MoS2 responsible for the spontaneous denaturation of proteins is likely exaggeratedsuggesting that MoS2 devices could, in fact, be used in biological applications. In the case of (iii), the water

■

THEORETICAL METHODS Partial Charge and Intersheet Potential Energy Calculations. Density functional theory (DFT) calculations on a 3D periodic system defined by the crystallographic unit cell of the 2H allotrope of MoS235 (space group P63/mmc, CIF file ID 9007660 from the Crystallography Open Database) were carried out using the Gaussian and Plane Waves (GPW) method of the Quickstep/CP2K package36,37 with a cutoff of 750 Ry. Calculations were carried out using the PBEsol functional38 (GGA), with DZVP-MOLOPT basis sets and the GTH pseudopotentials distributed with CP2K. The potential energy between two MoS2 sheets, at a series of interlayer distances, were calculated using CP2K, with exchange and correlation functionals for rVV10 from the libXC library.39 Contact-Angle Simulations. We utilized the open-source simulation package GROMACS 4.5,40 in the canonical (NVT) ensemble, to simulate the contact angle of water and diiodomethane (DIM) on monolayer, bilayer, and trilayer MoS2 surfaces. The water molecules were represented using the SPC/E water model.41 This three-site water model has been extensively used to parametrize the force fields for a variety of materials, including minerals and polymers, by reproducing experimental contact angles through MD simulations.42−44 We have also reported the contact angles obtained on monolayer MoS2 using the SPCE-F,45 TIP4P/2005,46 and TIP4P/Ice47 water models, along with the simulated surface tensions of these water models, in the Supporting Information. The force field for DIM was obtained from the Automated Topology Builder (ATB) and is based on the GROMOS 53A6 force field.48 As discussed in the Supporting Information, the surface tension of DIM predicted by this force field (49.2 mN/m) agrees very well with the experimentally determined value (50.8 mN/m). By decoupling van der Waals interactions from intralayer bonded interactions, our new force field enabled the use of geometric-mean combining rules for the LJ interactions between MoS2 and the liquids considered. Long-range Coulombic interactions were accounted for using the ParticleMesh Ewald (PME) technique.49 LJ interactions were cut-off at a distance of 1.2 nm. The Bussi−Parinello thermostat50 was used to maintain the system at a temperature of 298.15 K, and the SETTLE51/P-LINCS52 algorithm was used to keep rigid the bonds in the liquid molecules, with an integration time step of 2 fs. The atoms in the solid were kept frozen at their equilibrium positions.53 Water droplets containing N = 250, 9023

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

the PPPM method with an accuracy of 10−4 in energy. The simulations were carried out in a 10.0 nm × 10.0 nm × 5.0 nm simulation box filled with solvent molecules, and maintained at NPT conditions of 0.1 MPa and 300 K using a Nosé−Hoover barostat and thermostat. Forces on the two sheets in the direction normal to the sheets were computed every 20 fs and integrated using the trapezoidal rule to calculate the PMF between the two sheets in each solvent as a function of their separation d.

500, 1000, 2000, and 4000 water molecules were simulated, and the contact angle (θ) was obtained in each case by following the procedure proposed by Ingebrigtsen and Toxvaerd.54 The MDAnalysis Python module55 was utilized to compute the axial-radial density profiles, by averaging 1000 postequilibrium MD simulation snapshots, distributed over a period of 10 ns. The macroscopic contact angle was then obtained from the yintercept of the linear plot between cos θ and 1/R, where R is the radius of the circular contact line of the droplet and θ is the contact angle of the droplet containing N water molecules. In the Supporting Information, we have also reported the macroscopic contact angle obtained from the method proposed by Blake et al.56 for various water models on monolayer MoS2 modeled using the force field proposed in the present work, as well as for the SPC/E water model on monolayer MoS2 modeled using the LI and MO force fields. Liquid-Phase Exfoliation Simulations. For the liquidphase exfoliation simulations, the initial configuration consisted of a periodic simulation box containing one stack of four layers of MoS2 equilibrated in each of the following five solvents: Nmethylpyrrolidone (NMP), dimethyl sulfoxide (DMSO), dimethylformamide (DMF), methanol (MeOH), and water. The solvents considered were modeled using the OPLS-All Atom force field.57 Water was represented by the SPC/E model.41 The approximate size of the MoS2 stack is 5.0 nm by 4.0 nm with a thickness of 2.46 nm. Peeling takes place by pulling on the shorter edge of 4.0 nm (the armchair edge). The initial dimensions of the simulation box were 10.0 nm × 8.0 nm × 6.0 nm, thus allowing for a sufficient volume of solvent surrounding the stack, as shown in Figure 3a. The systems contained around 25 000 atoms. In all simulations, the cutoff for van der Waals interactions was taken as 1.2 nm, and longrange Coulombic interactions were accounted for using the PPPM method with an accuracy of 10−4 in energy. The simulation boxes containing the stack and the solvents were equilibrated for 5 ns at constant NPT conditions of 0.1 MPa and 350 K using a Nosé−Hoover barostat and thermostat. Note that a temperature of 350 K was chosen in order to accelerate equilibration of the system. All the solvents considered here have boiling points above this temperature. The potential of mean force (PMF) for peeling off the top layer was calculated by applying a harmonic biasing potential (of 70 kcal mol−1) distributed over Mo atoms at the peeled edge of the uppermost MoS2 sheet, and over the S atoms at the edge of the second MoS2 sheet underneath (Figure 3a). The biasing coordinate, d, is the distance between the centers of mass of these two sets of atoms. The PMF was calculated using umbrella sampling and the weighed-histogram analysis method (WHAM)58 implemented in the WHAM code developed by Grossfield and co-workers. The range of peeling separations, from 0.4 to 1.8 nm, was sampled every 0.025 nm. At each peeling separation, the system was equilibrated for 50 ps, and subsequently, data was collected for 150 ps. Reaggregation Simulations. The reaggregation simulations consisted of pulling two rigid 7.0 nm × 7.0 nm ABstacked MoS2 monolayers, immersed in a solvent, toward each other, using the LAMMPS MD code, from a separation of 1.5 nm to a separation of 0.275 nm. The sheets were pulled 0.25 Å closer to each other over 25 ps and then held in place at each fixed separation for 1 ns. Note that the sheets were kept parallel to each other throughout the simulation as depicted in Figure 4a. van der Waals interactions were cutoff of at 1.2 nm, and long-range Coulombic interactions were accounted for through

■

RESULTS AND DISCUSSION Force-Field Parameterization. We begin with a reparameterization of the partial charges for S and Mo atoms, using electronic structure density functional theory (DFT) calculations, as descried in the Methods sections. Instead of using Mulliken population analysis or Qeq electronegativity equalization to assign partial charges, as in the formulation of the existing force fields for MoS2 (Table 1), we utilized the DDAP Table 1. Partial Atomic Charges for MoS2 Calculated from DFT and Using the QEq Electronegativity Equalization Method literature Morita et al.22 Becker et al.21 Dallavale et al.61 this work 1 × 1 × 1 (bulk) 1 × 1 × 1 (bulk) 1 × 1 × 1 (bilayer) 3 × 3 × 1 (bulk, 54 atoms) 4 × 4 × 1 (bulk, 96 atoms) 6 × 6 × 2 (bulk, 432 atoms) 8 × 8 × 1 (bilayer)

method

qMo/e

qS/e

unknown Mulliken/B3LYP QEq QEq

+0.76 +0.70 +0.734 +0.74

−0.38 −0.35 −0.367 −0.37

DDAP/PBEsol DDAP/LDA DDAP/PBEsol DDAP/PBEsol DDAP/PBEsol DDAP/PBEsol QEq

+0.36 +0.38 +0.40 +0.48 +0.50 +0.50 +0.30

−0.18 −0.19 −0.20 −0.24 −0.25 −0.25 −0.15

formalism59 which is better suited for the derivation of point charges for dense periodic systems. Atomic charges were obtained for bulk 2H MoS2 and for a bilayer, by adding 10 Å of empty space above and below the unit cell. Note that removing the periodicity along z, using a bilayer, yields identical partial charges for S and Mo atoms (Table 1), as compared to the bulk material. We carried out calculations on an MoS2 bilayer because this is the smallest system that could be studied to determine the differences between the charges on the sulfur atoms on the surface and the sulfur atoms in the bulk of fewlayered MoS2 flakes. Apart from utilizing the PBE generalized gradient approximation (GGA) exchange−correlation functional, a local-density approximation (LDA) functional was also tested for comparison. Calculations were performed up to a 6 × 6 × 2 supercell, periodic in 3D, to study the convergence of the partial charges as a function of system size. We also calculated Qeq charges (Table 1) for a bilayer of MoS2 with an extent of 8 × 8 × 1 unit cells, using the LAMMPS60 molecular dynamics code. The electronegativity and self-Coulombic parameters were taken from the UFF force field62 and are listed in Table 2. In the bilayer systems, the surface and inner S atoms have partial charges that differ by less than 0.01e from those in both the DFT and the Qeq calculations, suggesting that surface effects on atomic charges are negligible. 9024

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

of the elastic constants, bulk modulus, and the Young’s modulus reported in Table 3 are calculated at the experimental lattice parameters, based on which the fit was performed. In the Supporting Information, we also report the corresponding values obtained for the fully relaxed lattice simulated using the force field proposed in this work. The root-mean-square deviation between the elastic properties obtained for the unrelaxed and relaxed lattices and the experimentally measured elastic properties are 0.300 and 0.407, respectively. To maintain compatibility with the OPLS-AA framework, the nonbonded (LJ and Coulombic) interactions between 1−2 and 1−3 neighbors were set to zero, and those between 1−4 neighbors were scaled by a factor of 0.5.57 Force-Field Validation. In order to validate our new force field for MoS2, we (a) compared the interaction potential between two MoS2 sheets, obtained using our force field, with that obtained using dispersion-corrected DFT calculations, and we (b) compared the simulated contact angle of water and diiodomethane droplets on MoS2 surfaces with the experimentally measured contact angles on MoS2.27,28 The potential energy of interaction between two AB-stacked, parallel monolayers of MoS2 (Figure 1b) was calculated using the PBE density functional, including dispersion through Grimme’s D3 semiclassical method,64 as well as using the rVV10 functional,65,66 which includes dispersion through a nonlocal term. The potential energy curve is shown in Figure 1c, where both the total PBE+D3 and PBE+rVV10 energies, and just the dispersion (D) and nonlocal (NL) terms, are plotted. It is noteworthy that although both DFT methods include dispersion interactions using different approaches, they yield very similar interaction energies. This finding, in turn, substantiates the accuracy of both DFT methods when calculating van der Waals forces. Figure 1c shows that the attraction between the MoS2 layers is almost exclusively due to van der Waals interactions because there is little difference between the full DFT energies and the dispersive or nonlocal contributions. This finding should be contrasted with our recent finding in the case of phosphorene,69 where a nonnegligible extent of orbital overlap between layers was observed, leading to a 25%-stronger interlayer cohesion relative to purely van der Waals interactions. The observed difference in interlayer cohesion in MoS2 and phosphorene reflects the fact that, unlike the phosphorus atoms in phosphorene, the sulfur atoms in MoS2 possess fully occupied orbitals that are all directed away from the surface and toward the Mo atoms.70 The potential energy curve obtained using our new force field is also plotted in Figure 1d, and it agrees very well with the energies obtained using both PBE+D3 and PBE+rVV10. This agreement is significant because no fitting of force-field parameters to the DFT energies was carried out. On the other hand, as shown in Figure 1d, the BE, MO, and VA force fields largely overestimate the minimum intersheet potential energy predicted by dispersion-corrected DFT. As discussed earlier, this overestimation is expected, given the very strong LJ terms used in these force fields in order to reproduce intralayer interactions. The derived force-field parameters for Mo and S should be consistent with experimental quantities that probe van der Waals interactions, such as the cleavage energy of MoS2 or its interfacial interactions with water. Our new force field predicts a surface energy of 240.75 mJ/m2 for 2H MoS2 (the depth of the minimum in Figure 1c), which is consistent with the values

Table 2. QEq Parameters for Mo and S Mo S

χ/eV

J/eV

3.465 6.928

7.510 8.972

Table 1 shows that, as the size of the bulk MoS2 system increases, the charges obtained from the DFT calculations converge to the following values: qMo = +0.50 and qS = −0.25. Consequently, in our new force field, we have chosen these values for the partial charges of Mo and S. Note that these partial charges are smaller than those reported in the literature by Morita et al.,22 Becker et al.,21 and Dallavale et al.61 The agreement between the DDAP charges (qMo = +0.50 and qS = −0.25) and the Qeq charges (qMo = +0.30 and qS = −0.15) is reasonable and consistent with the literature.59 Moreover, the values obtained here (qMo=+0.50 and qS = −0.25) are reasonable based on the electronegativity difference between Mo and S. After determining the partial charges, we obtained 6−12 LJ parameters and harmonic bond and angle force constants by fitting to the experimental lattice parameters, elastic constants, bulk modulus, and Young’s modulus of 2H MoS2 (Table 3), Table 3. Lattice Constants and Mechanical Properties of Bulk 2H MoS2 Calculated Using Our New Atomistic Force Field (Å)

force field

exp35

(GPa)

force field

exp67,68

a b c α β γ

3.163 3.163 12.159 90.000 90.000 120.000

3.161 3.161 12.295 90.000 90.000 120.000

C11 C12 C13 C33 C44 bulk modulus Young’s modulus

239.9 −52.3 −9.8 57.4 13.9 43.7 226.0

238 −54 23 52 19 43 240

using the GULP lattice dynamics code.63 Note that the 1 functional form of the bond potential is 2 kr(r − r0)2 , where kr is the harmonic bond constant and r0 is the equilibrium bond length. Similarly, the functional form of the angular potential is 1 k (θ − θ0)2 , where kθ is the harmonic angle constant and θ0 is 2 θ the equilibrium angle. As shown in Figure 1a, we have only considered the acute S−Mo−S angles in our model. The optimized parameters are listed in Table 4. A higher weight was placed on the C11 = C22 and C33 elastic constants because they are directly related to the rigidity of the layers of the material along x and y, and to the rigidity of the material along z, respectively. The list of weights utilized in the parametrization process is available in the Supporting Information. The values Table 4. Force-Field Parameters for Mo and S bonds

r0/Å

kr/(kJ mol−1 Å−2)

Mo−S angles

2.41 θ0/deg

430.3 kθ/(kJ mol−1 rad−2)

S−Mo−S Mo−S−Mo nonbonded

83.8 83.8

1187.7 2050.0

Mo S

σ/Å

ε/(kJ mol−1)

q/e

4.43 3.34

0.485 2.085

+0.50 −0.25 9025

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

in our simulations (1.2 nm), the dispersion interactions between a liquid droplet and a layer of MoS2 located beyond the third layer were found to be negligible. Moreover, Coulombic interactions are negligible for an extended monolayer MoS2 surface, even in the presence of polar solvents such as water. This is due to the exponential decay of the electric potential over the MoS2 basal plane, as shown in our recent work.53 As described in the Methods section, we used two different methods to compute the contact angle of water on MoS2 surfaces: (a) the method proposed by Ingebrigtsen and Toxvaerd,54 wherein the curvature of the first two layers of water above the solid surface determines the contact angle, and (b) the method proposed by Blake et al.,56 wherein a circular profile is fit to the droplet contour after excluding the first two layers of water above the solid surface, and the tangent to the circle at the solid surface determines the contact angle. In this section, we report the results using method (a). Results using method (b) are summarized in the Supporting Information, wherein we show that the force field developed in this work is in better agreement with experiment than the other MoS2 force fields reported in the literature, irrespective of the method used. The macroscopic contact angle of water on trilayer MoS2 (Figure 2d) was calculated to be 69.6° ± 0.8° (Figure 2e), which is very close to the experimental value of 69°.27,28 For DIM, we obtained a contact angle of 0°, whereas the experimental contact angle is 15°.27,28 Note that if we utilize the Young−Dupre equation to compute the work of adhesion W = γ(1 + cos θ), this 15° difference between the simulated and the experimentally measured contact angles of DIM on MoS2 corresponds to a difference of 1.7% between the simulated and the experimentally measured works of adhesion. The close agreement between the predicted and the experimental contact angles, without the need to adjust any parameters, suggests that our new force field for MoS2 can be utilized, along with available force fields for liquids and geometric-mean combining rules, to quantitatively predict solid−liquid interfacial phenomena. We also simulated the contact angle of SPC/E water using the BE, MO, VA, and LI force fields discussed above. Specifically, the BE and VA force fields predicted complete wetting with a contact angle of 0° (Figure 2a,b). This is due to the high ε parameters in these two force fields, reflecting the combination of intermolecular and intramolecular interactions. This suggests that the recent use of the VA force field to simulate the performance of MoS2 nanopores in water desalination applications33 may incorrectly estimate the water flux through MoS2 nanopores. On the other hand, although the MO force field yielded a water contact angle of 69.4° ± 0.9° (Figure 2b, Figure S1 in the Supporting Information), which is close to the experimentally measured value of 69°,27,28 we have already shown that the MO force field does not correctly model the intersheet interactions (Figure 1d) and hence should not be used. Finally, the LI force field yielded a water contact angle of 97.9° ± 0.8° (Figure 2c, Figure S1 in the Supporting Information), indicating that it underestimates the water−MoS2 interactions which promote wetting. This underestimation, combined with the significant difference in the LJ parameters of Mo and S in the LI force field, suggests that the use of the LI force field to simulate water desalination through nanoporous MoS25 could lead to incorrect predictions of membrane performance. Given the excellent performance of the new force field in modeling the interfacial behavior (contact angles) of MoS2/

Figure 1. (a) New force-field parameters for MoS2 include the parameters σ and ε of the LJ potential for sulfur (S, shown in yellow) and molybdenum (Mo, shown in pink), r0 and kr for the harmonic bond rMo−S, and kθ and θ0 for the harmonic potential on the angles θMo−S−Mo and θS−Mo−S. (b) The 2H allotrope of MoS2 considered throughout this study. (c) Potential energy of interaction between two parallel single layers of 2H MoS2 calculated by DFT and using the new atomistic force field presented here. Blue squares correspond to PBE +D3 energies, and green dots represent the dispersion contribution. Red squares correspond to rVV10 energies, and yellow dots depict the contribution from the nonlocal functional. The black line corresponds to the new force field developed here. (d) Comparison of the potential energies between two parallel single layers of 2H MoS2 calculated using the force field developed here and the BE, MO, VA, and LI force fields.

of: (i) 164.46 mJ/m2 obtained by Bjorkman et al.71 using advanced DFT (adiabatic-connection fluctuation−dissipation theorem within the random-phase approximation (RPA)), (ii) 284 mJ/m2 obtained by Fuhr et al.72 using the full potential linearized augmented plane wave (FP-LAPW) method with LDA for exchange and correlation, and (iii) 260 mJ/m2 obtained by Weiss et al.73 by fitting an empirical LJ potential to the bulk modulus and interlayer separation of 2H MoS2. Note that Bjorkman et al.71 also used the VV10 functional to compute the surface energy of 2H MoS2, and their estimate of 246.9 mJ/m2 is in excellent agreement with the surface energy predicted by our force field. The contact angle of liquids on a surface provides a convenient way to quantitatively probe the magnitude of the intermolecular interactions operating between the surface of the material and the liquid. Here, we validated our new force field by comparing the contact angles of water (a polar solvent) and diiodomethane (DIM, a nonpolar solvent) on monolayer, bilayer, and trilayer MoS2, predicted using MD simulations, against their experimentally measured values, as reported by Kozbial et al. and Gurarslan et al.27,28 Note that trilayer MoS2 is representative of bulk MoS2 because, for the LJ cutoff utilized 9026

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

Figure 2. Representative equilibrium MD simulation snapshots of droplets composed of 4000 SPC/E water molecules on 2H MoS2 trilayers simulated using different force-field models. Atoms of molybdenum, sulfur, oxygen, and hydrogen are shown in pink, yellow, red, and white, respectively. Estimates of the macroscopic contact angle θ of SPC/E water on bulk MoS2 simulated using these force fields are reported on the topright of the corresponding snapshot. (a) SPC/E water completely wets (zero contact angle) MoS2 trilayers modeled using the BE and VA force fields. (b) The MO force field yields a contact angle of 69.4° ± 0.9°. (c) The LI force field renders MoS2 hydrophobic and yields a contact angle of 97.9° ± 0.8°. (d) The new force field for MoS2 proposed here yields a contact angle of 69.6° ± 0.8° on trilayer MoS2, in excellent agreement with experiments.27,28 (e) Variations of the water contact angle θ (left-hand y axis) and of the cosine of the water contact angle cos θ (right-hand y axis) with the inverse of the droplet radius R−1 for SPC/E water droplets on monolayer (ML), bilayer (BL), and trilayer (TL) MoS2. Note that the intercept of the line of least-squares fit to cos θ vs R−1, when R → ∞, is used to estimate the water contact angle for a macroscopic droplet (R → ∞).

possibly because the point of application of the biasing potential is too distant from the front where the layers peel away and because the histograms provide less significant information. We have determined that the reproducibility of the PMF curves is within 2.5%, as illustrated by the two PMF curves plotted for NMP in Figure S2 in the Supporting Information, which were obtained starting from different initial configurations. Simulating the Reaggregation of Dispersed MoS2 Sheets. The second in silico experiment corresponds to the reaggregation of the exfoliated MoS2 monolayers in various solvents. Two rigid 7.0 nm × 7.0 nm AB-stacked MoS2 monolayers, immersed in a solvent, were pulled toward each other from a separation of 1.5 nm to a separation of 0.275 nm. The simulated PMF curves for MoS2 sheets in vacuum, isopropanol (IPA), water, DMSO, DMF, and NMP are shown in Figure 4b. Note that aggregation corresponds to bringing the two sheets together from the extreme right of Figure 4b to the global, primary minimum at the left end of each PMF curve. The presence of the solvent molecules has several effects on the PMF curve. First, for all solvents except water, the intercalated solvent molecules “screen” the intersheet interactions. Therefore, the depth of the PMF well at the primary minimum is smaller when the sheets are immersed in a solvent medium (other than water) than when the sheets interact in vacuum. This stabilization reduces the driving force responsible for aggregation in these media. Second, the presence of the solvent induces an energy barrier between the primary and the secondary minima (see Figure 4b). The difference between the PMFs at the primary maxima and at the secondary minima corresponds to the amount of work required to force out the first solvation shell from between the two sheets. This energy barrier must be overcome during the aggregation of the monolayers, and the height of this barrier determines the energy requirements for this process. The differences in the heights of this free-energy barrier account for the variation in the colloidal stability of dispersed nanomaterial sheets in different solvents. Indeed, the larger the

liquid systems, we proceeded to test its ability to obtain molecular-level insights on the exfoliation and aggregation processes involving MoS2. For this purpose, we carried out two types of “in silico experiments,” where we steered model systems to simulate exfoliation and aggregation in order to calculate the relevant energetics. Simulating the Liquid-Phase Exfoliation of MoS2. The first in silico experiment consists of peeling one layer of MoS2 from a stack of several MoS2 layers, in different solvents, and is representative of the liquid-phase exfoliation process. We calculate the potential of mean force (PMF), that is, the reversible work, required to peel this one layer, to different extents, from the stack, while allowing full flexibility of the material. The PMF curves for peeling the MoS2 stack in the different solvents are shown in Figure 3b. At small separations, there is little difference between the curves, as shown in Figure 3c, because solvent molecules have not yet intercalated between the first and the second sheets. At separations beyond 0.7 nm, differences between the solvents become evident. Exfoliation will be easier in solvents with lower PMF values. Water is the least favorable of the solvents tested, followed by methanol, and then NMP. The most favorable solvents for the exfoliation process are DMF and DMSO, which are indistinguishable based on our simulation results. This order (DMF ≈ DMSO > NMP > MeOH > water) agrees reasonably well with the one observed in exfoliation experiments.74 The experimentally determined ranking of MoS2 dispersion media is NMP > DMF > DMSO > MeOH > water. Therefore, NMP is experimentally observed to be more favorable than DMSO or DMF in yielding a high concentration of dispersed few-layered MoS2 sheets. Note that, as explained below, this discrepancy could be explained by NMP being better than DMSO or DMF at preventing the aggregation of MoS2 sheets after liquid-phase exfoliation. We found that after a certain separation (which depends on the system), it becomes difficult to converge the WHAM algorithm used to calculate the PMF (see the Methods section), 9027

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

Figure 4. (a) Setup for the parallel plate dispersion stability simulations with a few representative solvent (NMP) molecules shown. Two AB-stacked MoS2 sheets, 7.0 nm by 7.0 nm each, are immersed in a 10.0 nm by 10.0 nm by 5.0 nm simulation box containing solvent molecules. Atoms of molybdenum, sulfur, oxygen, nitrogen, carbon, and hydrogen are shown in pink, yellow, red, blue, cyan, and white, respectively. The forces experienced by the two sheets in the direction normal to the sheets are computed as a function of the separation d between the sheets. (b) PMF U(d) curves corresponding to the two MoS2 sheets in vacuum (dashed line) and in the five solvents considered here.

in Figure 4b predict (see section S6 in the Supporting Information for details) the following order for the colloidal suspension lifetimes: NMP > DMSO > DMF > IPA > H2O, in excellent agreement with the experimental reports.74 It is also interesting to note that the stabilization afforded to MoS2 sheets in water when they are at their optimum separation (i.e., at the primary minimum in the PMF curve) is greater than that in vacuum. This allows us to rationalize an observation that has puzzled researchers since Morrison and co-workers first investigated dispersions of exfoliated MoS278contrary to intuition, water is a poor solvent for MoS2, in spite of the fact that MoS2 sheets have polar bonds.53 Interaction Energies between MoS2 Sheets and Exfoliation Media. We expect that the interactions between large (micron-sized) MoS2 sheets will be dominated by van der Waals interactions, with the Coulombic interactions playing an important role only at the edges.53 Indeed, we recently showed that the armchair and zigzag planes in MoS2, which are the exposed planes at edges of nanometer-size flakes, are polar. Therefore, for smaller MoS2 sheets, such as the nanometersized sheets considered here, Coulombic interactions could be significant. In order to assess the relative strengths of these two long-range interactions for nanometer-sized sheets, we computed the energy difference between a set of four, approximately 5.0 nm by 4.0 nm each, fully exfoliated MoS2 sheets (final state) and a four-layer stack of 2H MoS2 (initial state) in the different solvents considered here (Figure S3 in the Supporting Information). Note that in the final state considered, the separations are at least 2 times the cutoff radius for van der Waals interactions. This ensures that sheet−sheet interactions are negligible and that the surfaces of all the MoS2 sheets are entirely exposed to solvent. The averages of the energy differences (ΔEtot), including their Coulombic (ΔECoul) and van der Waals (ΔEvdW) components, obtained from 5 ns

Figure 3. (a) Setup of the simulation box for the peeling PMF calculations. Atoms of molybdenum, sulfur, oxygen, nitrogen, carbon, and hydrogen are shown in pink, yellow, red, blue, cyan, and white, respectively. Four layers of MoS2 of approximately 5.0 nm by 4.0 nm each are immersed in 2000 (NMP) to 10000 (water) molecules. The biasing potential is applied with respect to the distance d as defined in the text. (b) PMF associated with peeling one layer of MoS2 from a stack of MoS2 in vacuum and in the five solvents considered here. The curves shown have all converged to a tolerance of at least 10−5 using the WHAM method. (c) Close-up view of the PMF associated with peeling for d < 0.6 nm.

energy barrier, the more stable the colloidal dispersion. It should be noted that in this in silico experiment, we force the two MoS2 sheets to remain rigid and parallel to each other throughout the simulation. In practice, as reported by Shih et al., 75 the flexibility of the sheets and their entropic fluctuations76 could allow them to traverse a more freeenergetically favorable path toward aggregation. Altabet and Debenedetti have also predicted that accounting for flexibility results in an increase in the critical drying distance, suggesting that flexible MoS2 sheets will have a greater tendency to aggregate than the rigid sheets modeled in this work.77 Consequently, our PMF curves provide upper bounds to the free-energy barriers to aggregation. Accordingly, the lifetimes of colloidal suspensions of MoS2 could be shorter than those predicted by applying transition-state theory75 to the curves in Figure 4b. Nevertheless, our simulations reveal significant distinctions between the efficacy of commonly utilized solvents in preventing the aggregation of MoS2. Specifically, along with the transition-state-theory-based framework for the kinetics of colloidal aggregation proposed by Shih et al.,75 the PMF curves 9028

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C MD trajectories at 350 K, are reported in Table 5. Note that these energy differences are different from the free-energy

sheets. The ability to use these exfoliation and aggregation PMF curves to obtain a ranking of solvents that is in excellent agreement with experimental observations, provides yet another validation of our new force field for MoS2. The new force field for MoS2 developed here should enable the computational exploration, characterization, and optimization of MoS2-based devices in a variety of applications where interfacial interactions are salient, including semipermeable membranes, energy storage systems, microfluidic sensors, and electrolyte-gated field effect transistors. We anticipate that the new MoS2 force field will also pave the way for engineering newer applications, including the self-assembly of different layered nanomaterials into van der Waals heterostructures and the study of gas adsorption on MoS2 for sensory applications.

Table 5. Energy Differences between Four Fully-Exfoliated MoS2 Sheets and a Four-Layer Stack Immersed in the Five Solvents Considered Herea solvent

ΔECoul

ΔEvdW

ΔEtot

(kJ mol−1) NMP DMSO DMF EtOH H2O vacuum

−526 −85 −840 222 2798 2066

4262 6205 4373 9170 9938 18133

3658 6022 3664 9753 12983 12987

■

a

The simulation boxes contain between 3000 and 18000 solvent molecules. In the exfoliated state, the sheets are separated by 2.4 nm of solvent. The total energy differences are listed, in addition to the Coulombic and the van der Waals components. Uncertainties in energies are of the order of 150 kJ mol−1.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b00484. Simulated water contact angles for the Morita (MO) and REBO (RE) MoS2 force fields, contact angles of cylindrical water droplets evaluated using two different methods for various water models and MoS2 force fields, surface tension of diiodomethane evaluated using MD simulations, reproducibility of the peeling simulations, representations of the simulation boxes used to calculate energy differences between a four-layer stack of MoS2 sheets and four fully exfoliated MoS2 sheets, application of the kinetic theory for colloidal aggregation, weights used in force-field parametrization, and calculated elastic constants at the relaxed geometry of the MoS2 unit cell (PDF)

differences obtained using the PMF calculations reported in Figure 3 and Figure 4. Indeed the free energies include not only energetic but also entropic effects. For the five solvents considered here, Table 5 shows that the van der Waals component is dominant relative to the Coulombic component. This finding is consistent with our recent work on the wettability of the basal plane of 2H MoS2,53 which, for infinitely large sheets, does not depend on the Coulombic charges in MoS2. In the present calculations of energy differences for nanometer-sized MoS2 sheets, the presence of exposed Mo or S atoms in our model of the MoS2 monolayers results in a large electric field at the edges of the MoS2 sheets. As a result, strongly polar solvents are expected to interact more strongly (i.e., with a more negative ECoul) with the edges of the MoS2 sheets. Our final state exposes the edges of the MoS2 sheets to more solvent relative to their exposure in our initial state. Consequently, the more polar the solvent, the more negative is the value of ΔECoul = initial Efinal Coul − ECoul . The ranking of the solvents in terms of their ΔECoul seen in Table 5 is consistent with the experimentally measured dipole moments of the solvents. These findings are also significant in the context of the intercalation (or expulsion) of solvent molecules into (or from) the gap between MoS2 layers: polar solvents favor the presence of exposed edges, thereby promoting exfoliation and hindering aggregation.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. ORCID

Ananth Govind Rajan: 0000-0003-2462-0506 Daniel Blankschtein: 0000-0002-7836-415X Author Contributions ¶

V.S. and A.G.R. contributed equally to this work.

Notes

The authors declare no competing financial interest.

■

■

ACKNOWLEDGMENTS V.S., A.G.R., M.S.S., and D.B. acknowledge the National Science Foundation (NSF) CBET-1511526 for financial support. A.G.R. and M.S.S. acknowledge the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award Grant DE-FG02-08ER46488 Mod 0008. This work used the Extreme Science and Engineering Discovery Environment (XSEDE),79 which is supported by NSF grant number ACI1053575. A.P. acknowledges support from the Institut Universitaire de France and from the Agence Nationale de la Recherche, project CLINT (ANR-12-IS10-003). E.B.’s stay at MIT was funded by ANR project CLINT.

CONCLUSIONS In this manuscript, we have developed a new atomistic force field for MoS2 using a combination of lattice dynamics and quantum mechanical calculations. The new force field distinguishes between the intralayer bonded interactions, which contribute to the elasticity of each monolayer, and the nonbonded interactions, which determine the interlayer cohesion between MoS2 monolayers as well as the adhesion between MoS2 sheets and other materials and solvents. We have shown that our new force field yields superior performance, when compared to force fields which are available in the literature, in describing interfacial phenomena, including the contact angles of water and diiodomethane on MoS2. Finally, we have utilized our new force field to carry out PMF calculations to shed light on the efficacy of various solvents for both the liquid-phase exfoliation and the dispersion of MoS2

■

REFERENCES

(1) Radisavljevic, B.; Radenovic, A.; Brivio, J.; Giacometti, V.; Kis, A. Single-layer MoS2 Transistors. Nat. Nanotechnol. 2011, 6, 147−50.

9029

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C

(23) Varshney, V.; Patnaik, S. S.; Muratore, C.; Roy, A. K.; Voevodin, A. A.; Farmer, B. L. MD Simulations of Molybdenum Disulphide (MoS2): Force-Field Parameterization and Thermal Transport Behavior. Comput. Mater. Sci. 2010, 48, 101−108. (24) Liang, T.; Phillpot, S. R.; Sinnott, S. B. Parametrization of a Reactive Many-Body Potential for Mo-S Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 245110. (25) Onodera, T.; Morita, Y.; Suzuki, A.; Koyama, M.; Tsuboi, H.; Hatakeyama, N.; Endou, A.; Takaba, H.; Kubo, M.; Dassenoy, F.; et al. A Computational Chemistry Study on Friction of h-MoS2. Part I. Mechanism of Single Sheet Lubrication. J. Phys. Chem. B 2009, 113, 16526−16536. (26) Jiang, J.-W.; Park, H. S.; Rabczuk, T. Molecular Dynamics Simulations of Single-Layer Molybdenum Disulphide (MoS2): Stillinger-Weber Parametrization, Mechanical Properties, and Thermal Conductivity. J. Appl. Phys. 2013, 114, 064307. (27) Kozbial, A.; Gong, X.; Liu, H.; Li, L. Understanding the Intrinsic Water Wettability of Molybdenum Disulfide (MoS2). Langmuir 2015, 31, 8429−8435. (28) Gurarslan, A.; Jiao, S.; Li, T.-D.; Li, G.; Yu, Y.; Gao, Y.; Riedo, E.; Xu, Z.; Cao, L. Van der Waals Force Isolation of Monolayer MoS2. Adv. Mater. 2016, 28 (45), 10055−10060. (29) Luan, B.; Zhou, R. Wettability and Friction of Water on a MoS2 Nanosheet. Appl. Phys. Lett. 2016, 108, 131601. (30) Gaur, A. P. S.; Sahoo, S.; Ahmadi, M.; Dash, S. P.; Guinel, M. J.f.; Katiyar, R. S. Surface Energy Engineering for Tunable Wettability through Controlled Synthesis of MoS2. Nano Lett. 2014, 14, 4314−21. (31) Ling, Y.; Gu, Z.; Kang, S.-G.; Luo, J.; Zhou, R. Structural Damage of a β-Sheet Protein upon Adsorption onto Molybdenum Disulfide Nanotubes. J. Phys. Chem. C 2016, 120, 6796−6803. (32) Gu, Z.; Li, W.; Hong, L.; Zhou, R. Exploring Biological Effects of MoS2 Nanosheets on Native Structures of α-helical Peptides. J. Chem. Phys. 2016, 144, 175103. (33) Li, W.; Yang, Y.; Weber, J. K.; Zhang, G.; Zhou, R. Tunable, Strain-Controlled Nanoporous MoS2 Filter for Water Desalination. ACS Nano 2016, 10, 1829−1835. (34) Leroy, F. Revisiting the droplet simulation approach to derive force-field parameters for water on molybdenum disulfide from wetting angle measurements. J. Chem. Phys. 2016, 145, 164705. (35) Schonfeld, B.; Huang, J. J.; Moss, S. C. Anisotropic MeanSquare Displacements (MSD) in Single-Crystals of 2H-MoS2 And 3RMoS2. Acta Crystallogr., Sect. B: Struct. Sci. 1983, 39, 404−407. (36) Vandevondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (37) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k: atomistic simulations of condensed matter systems. WIREs: Comput. Mol. Sci. 2014, 4, 15−25. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (39) Marques, M. A. L.; Oliveira, M. J. T.; Burnus, T. Libxc: A Library of Exchange and Correlation Functionals for Density Functional Theory. Comput. Phys. Commun. 2012, 183, 2272−2281. (40) Pronk, S.; Pouya, I.; Lundborg, M.; Rotskoff, G.; Wesén, B.; Kasson, P. M.; Lindahl, E. Molecular Simulation Workflows as Parallel Algorithms: The Execution Engine of Copernicus, a Distributed HighPerformance Computing Platform. J. Chem. Theory Comput. 2015, 11, 2600−2608. (41) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (42) Šolc, R.; Gerzabek, M. H.; Lischka, H.; Tunega, D. Wettability of Kaolinite (001) Surfaces − Molecular Dynamic Study. Geoderma 2011, 169, 47−54. (43) Hirvi, J. T.; Pakkanen, T. A. Molecular Dynamics Simulations of Water Droplets on Polymer Surfaces. J. Chem. Phys. 2006, 125, 144712. (44) Emami, F. S.; Puddu, V.; Berry, R. J.; Varshney, V.; Patwardhan, S. V.; Perry, C. C.; Heinz, H. Force Field and a Surface Model

(2) Eda, G.; Yamaguchi, H.; Voiry, D.; Fujita, T.; Chen, M.; Chhowalla, M. Photoluminescence from Chemically Exfoliated MoS2. Nano Lett. 2011, 11, 5111−5116. (3) Bertolazzi, S.; Brivio, J.; Kis, A. Stretching and Breaking of Ultrathin MoS2. ACS Nano 2011, 5, 9703−9709. (4) Liu, K.; Feng, J.; Kis, A.; Radenovic, A. Atomically Thin Molybdenum Disulfide Nanopores with High Sensitivity for DNA Translocation. ACS Nano 2014, 8, 2504−2511. (5) Heiranian, M.; Farimani, A. B.; Aluru, N. R. Water Desalination with a Single-Layer MoS2 Nanopore. Nat. Commun. 2015, 6, 8616. (6) Li, H. H.; Yin, Z.; He, Q.; Li, H.; Huang, X.; Lu, G.; Fam, D. W. H.; Tok, A. I. Y.; Zhang, Q.; Zhang, H. Fabrication of Single- and Multilayer MoS2 Film-Based Field-Effect Transistors for Sensing NO At Room Temperature. Small 2012, 8, 63−67. (7) Govind Rajan, A.; Warner, J. H.; Blankschtein, D.; Strano, M. S. Generalized Mechanistic Model for the Chemical Vapor Deposition of 2D Transition Metal Dichalcogenide Monolayers. ACS Nano 2016, 10 (4), 4330−4344. (8) Finn, D. J.; Lotya, M.; Cunningham, G.; Smith, R. J.; McCloskey, D.; Donegan, J. F.; Coleman, J. N. Inkjet Deposition of LiquidExfoliated Graphene and MoS2 Nanosheets for Printed Device Applications. J. Mater. Chem. C 2014, 2, 925−932. (9) Torrisi, F.; Hasan, T.; Wu, W.; Sun, Z.; Lombardo, A.; Kulmala, T. S.; Hsieh, G. W.; Jung, S.; Bonaccorso, F.; Paul, P. J.; et al. InkjetPrinted Graphene Electronics. ACS Nano 2012, 6, 2992−3006. (10) Li, J.; Ye, F.; Vaziri, S.; Muhammed, M.; Lemme, M. C.; Ö stling, M. Efficient Inkjet Printing of Graphene. Adv. Mater. 2013, 25, 3985− 3992. (11) Withers, F.; Yang, H.; Britnell, L.; Rooney, A. P.; Lewis, E.; Felten, A.; Woods, C. R.; Sanchez Romaguera, V.; Georgiou, T.; Eckmann, A.; et al. Heterostructures Produced from Nanosheet-Based Inks. Nano Lett. 2014, 14, 3987−3992. (12) Coleman, J. N.; Lotya, M.; O’Neill, A.; Bergin, S. D.; King, P. J.; Khan, U.; Young, K.; Gaucher, A.; De, S.; Smith, R. J.; et al. TwoDimensional Nanosheets Produced by Liquid Exfoliation of Layered Materials. Science 2011, 331, 568−571. (13) Smith, R. J.; King, P. J.; Lotya, M.; Wirtz, C.; Khan, U.; De, S.; O’Neill, A.; Duesberg, G. S.; Grunlan, J. C.; Moriarty, G.; et al. LargeScale Exfoliation of Inorganic Layered Compounds in Aqueous Surfactant Solutions. Adv. Mater. 2011, 23, 3944−3948. (14) Guan, G.; Zhang, S.; Liu, S.; Cai, Y.; Low, M.; Teng, C. P.; Phang, I. Y.; Cheng, Y.; Duei, K. L.; Srinivasan, B. M.; et al. Protein Induces Layer-by-Layer Exfoliation of Transition Metal Dichalcogenides. J. Am. Chem. Soc. 2015, 137, 6152−6155. (15) Huang, Y.; Shi, Y.; Yang, H. Y.; Ai, Y. A Novel Single-Layered MoS2 Nanosheet based Microfluidic Biosensor for Ultrasensitive Detection of DNA. Nanoscale 2015, 7, 2245−9. (16) Feng, C.; Ma, J.; Li, H.; Zeng, R.; Guo, Z.; Liu, H. Synthesis of Molybdenum Disulfide (MoS2) for Lithium Ion Battery Applications. Mater. Res. Bull. 2009, 44, 1811−1815. (17) Cao, T.; Wang, G.; Han, W.; Ye, H.; Zhu, C.; Shi, J.; Niu, Q.; Tan, P.; Wang, E.; Liu, B.; et al. Valley-Selective Circular Dichroism of Monolayer Molybdenum Disulphide. Nat. Commun. 2012, 3, 887. (18) Perera, M. M.; Lin, M. W.; Chuang, H. J.; Chamlagain, B. P.; Wang, C.; Tan, X.; Cheng, M. M. C.; Tománek, D.; Zhou, Z. Improved Carrier Mobility in Few-Layer MoS2 Field-Effect Transistors with Ionic-Liquid Gating. ACS Nano 2013, 7, 4449−4458. (19) Nicolini, P.; Polcar, T. A Comparison of Empirical Potentials for Sliding Simulations of MoS2. Comput. Mater. Sci. 2016, 115, 158−169. (20) Based on the number of citing articles as reported by Google Scholar on April 13,2016. (21) Becker, U.; Rosso, K. M.; Weaver, R.; Warren, M.; Hochella, M. F. Metal Island Growth and Dynamics on Molybdenite Surfaces. Geochim. Cosmochim. Acta 2003, 67, 923−934. (22) Morita, Y.; Onodera, T.; Suzuki, A.; Sahnoun, R.; Koyama, M.; Tsuboi, H.; Hatakeyama, N.; Endou, A.; Takaba, H.; Kubo, M.; et al. Development of a New Molecular Dynamics Method for Tribochemical Reaction and its Application to Formation Dynamics of MoS2 Tribofilm. Appl. Surf. Sci. 2008, 254, 7618−7621. 9030

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031

Article

The Journal of Physical Chemistry C Database for Silica to Simulate Interfacial Properties in Atomic Resolution. Chem. Mater. 2014, 26, 2647−2658. (45) Tironi, I. G.; Brunne, R. M.; van Gunsteren, W. F. On the Relative Merits of Flexible versus Rigid Models for Use in Computer Simulations of Molecular Liquids. Chem. Phys. Lett. 1996, 250, 19−24. (46) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (47) Abascal, J. L. F.; Sanz, E.; Garcia Fernandez, R.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/ Ice. J. Chem. Phys. 2005, 122, 234511. (48) Malde, A. K.; Zuo, L.; Breeze, M.; Stroet, M.; Poger, D.; Nair, P. C.; Oostenbrink, C.; Mark, A. E. An Automated Force Field Topology Builder (ATB) and Repository: Version 1.0. J. Chem. Theory Comput. 2011, 7, 4026−4037. (49) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (50) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (51) Miyamoto, S.; Kollman, P. A. An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (52) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (53) Govind Rajan, A.; Sresht, V.; Padua, A. A. H.; Strano, M. S.; Blankschtein, D. Dominance of Dispersion Interactions and Entropy over Electrostatics in Determining the Wettability and Friction of Two-Dimensional MoS2 Surfaces. ACS Nano 2016, 10 (10), 9145− 9155. (54) Ingebrigtsen, T. S.; Dyre, J. C. The Impact Range for Smooth Wall-Liquid Interactions in Nanoconfined Liquids. Soft Matter 2014, 10, 4324−4331. (55) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 2011, 32, 2319−2327. (56) Blake, T. D.; Clarke, A.; De Coninck, J.; de Ruijter, M. J. Contact Angle Relaxation during Droplet Spreading: Comparison between Molecular Kinetic Theory and Molecular Dynamics. Langmuir 1997, 13, 2164−2166. (57) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (58) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (59) Blochl, P. E. Electrostatic Decoupling of Periodic Images of Plane-Wave-Expanded Densities and Derived Atomic Point Charges. J. Chem. Phys. 1995, 103, 7422. (60) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (61) Dallavalle, M.; Sändig, N.; Zerbetto, F. Stability, Dynamics, and Lubrication of MoS2 Platelets and Nanotubes. Langmuir 2012, 28, 7393−7400. (62) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, A Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (63) Gale, J. D. GULP: A Computer Program for the SymmetryAdapted Simulation of Solids. J. Chem. Soc., Faraday Trans. 1997, 93, 629−637. (64) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab-initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (65) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals Density Functional: The Simpler the Better. J. Chem. Phys. 2010, 133, 244103.

(66) Sabatini, R.; Gorni, T.; De Gironcoli, S. Nonlocal van der Waals Density Functional Made Simple and Efficient. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 041108. (67) Webb, A.; Feldman, J.; Skelton, E.; Towle, L.; Liu, C.; Spain, I. High Pressure Investigations of MoS2. J. Phys. Chem. Solids 1976, 37, 329−335. (68) Feldman, J. Elastic Constants of 2H-MoS2 and 2H-NbSe2 Extracted from Measured Dispersion Curves and Linear Compressibilities. J. Phys. Chem. Solids 1976, 37, 1141−1144. (69) Sresht, V.; Pádua, A. A. H.; Blankschtein, D. Liquid-Phase Exfoliation of Phosphorene: Design Rules from Molecular Dynamics Simulations. ACS Nano 2015, 9, 8255−8268. (70) Didziulis, S. V.; Fleischauer, P. D.; Soriano, B. L.; Gardos, M. N. Chemical and Tribological Studies of MoS2 Films on SiC Substrates. Surf. Coat. Technol. 1990, 43−44, 652−662. (71) Björkman, T.; Gulans, A.; Krasheninnikov, A. V.; Nieminen, R. M. van der Waals Bonding in Layered Compounds from Advanced Density-Functional First-Principles Calculations. Phys. Rev. Lett. 2012, 108, 235502. (72) Fuhr, J.; Sofo, J.; Saúl, A. Adsorption of Pd on MoS2(1000): Abinitio Electronic-Structure Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 60, 8343−8347. (73) Weiss, K.; Phillips, J. Calculated Specific Surface Energy of Molybdenite (MoS2). Phys. Rev. B 1976, 14, 5392−5395. (74) Cunningham, G.; Lotya, M.; Cucinotta, C. S.; Sanvito, S.; Bergin, S. D.; Menzel, R.; Shaffer, M. S. P.; Coleman, J. N. Solvent Exfoliation of Transition Metal Dichalcogenides: Dispersibility of Exfoliated Nanosheets Varies Only Weakly Between Compounds. ACS Nano 2012, 6, 3468−80. (75) Shih, C.-J.; Lin, S.; Strano, M. S.; Blankschtein, D. Understanding the Stabilization of Liquid-Phase-Exfoliated Graphene in Polar Solvents: Molecular Dynamics Simulations and Kinetic Theory of Colloid Aggregation. J. Am. Chem. Soc. 2010, 132, 14638−48. (76) Ulissi, Z.; Govind Rajan, A.; Strano, M. S. Persistently Auxetic Materials: Engineering the Poisson Ratio of 2D Self-Avoiding Membranes under Conditions of Non-Zero Anisotropic Strain. ACS Nano 2016, 10 (8), 7542−7549. (77) Altabet, Y. E.; Debenedetti, P. G. The Role of Material Flexibility on the Drying Transition of Water Between Hydrophobic Objects: A Thermodynamic Analysis. J. Chem. Phys. 2014, 141, 18C531. (78) Divigalpitiya, W. M.; Frindt, R. F.; Morrison, S. R. Inclusion Systems of Organic Molecules in Restacked Single-layer Molybdenum Disulfide. Science 1989, 246, 369−71. (79) Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; et al. XSEDE: Accelerating Scientific Discovery. Comput. Sci. Eng. 2014, 16, 62−74.

9031

DOI: 10.1021/acs.jpcc.7b00484 J. Phys. Chem. C 2017, 121, 9022−9031