Parameterization of Molybdenum Disulfide Interacting with Water

Subscriber access provided by BUFFALO STATE

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

Parameterization of Molybdenum Disulfide Interacting with Water Using Free Energy Perturbation Method Leili Zhang, Binquan Luan, and Ruhong Zhou J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b02797 • Publication Date (Web): 01 Aug 2019 Downloaded from pubs.acs.org on August 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Parameterization of Molybdenum Disulfide Interacting with Water Using Free Energy Perturbation Method Leili Zhang, Binquan Luan and Ruhong Zhou* Computational Biology Center, IBM Thomas J. Watson Research Center, Yorktown Heights, NY 10598, USA

*All Correspondence should be addressed to: [email protected]

Abstract Water contact angles (WCA) are often used to parametrize force field parameters of novel 2D nanomaterials, such as molybdenum disulfide (MoS2), which has emerged as a promising nanomaterial in many biomedical applications due to its unique and impressive properties. However, there is a wide range of water-MoS2 contact angles in the literature depending on the aging process on the surface of a MoS2 nanosheet and/or substrate materials. In this study, we revisit and optimize existing parameters for the basal plane of MoS2 with two popular water models, TIP3P and SPC/E, using the wide range of WCAs from various experiments. We develop and deploy the free energy perturbation method for parameterizing MoS2 with experimentally determined WCAs for both fresh and aged surfaces. Energy decomposition analysis on the simulation trajectories reveals that MoS2-water interaction is dominated by van der Waals interaction, which mainly comes from the top layer of MoS2. We conclude that to describe both fresh and aged MoS2 surfaces it is convenient to only adjust the Lennard-Jones parameter εS (the depth of the potential well of a sulfur atom), which displays a surprisingly linear correlation with WCAs.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Molybdenum disulfide (MoS2) is a naturally occurring solid inorganic compound classified as a transition metal dichalcogenide (TMD). Like graphene, monolayer MoS2 belongs to the 2D materials family, but its three-atom-thick monolayer is composed of a honeycomb sheet of molybdenum atoms sandwiched between two honeycomb sheets of sulfur atoms with strong covalent bonding between neighboring sheets. The structure of bulk MoS2 comprises multiple monolayers stacked vertically and held together by rather weak van der Waals (VDW) forces. Figure 1 depicts the atomic structure of an example trilayer MoS2 used in this study. Traditionally, due to its abundance in the Earth’s crust, bulk MoS2 has been used for many years in conventional industries as dry lubricants1 in mechanics or cathode materials for lithium ion batteries. However, recent discovery of the capability of realizing 2D nanomaterials from layered crystalline materials by exfoliating into individual layers which exhibit unique and fascinating properties quite distinct from the bulk counterparts has attracted great interest from the scientific community to further explore TMD materials such as MoS2. As a result, over the last decades, a spectrum of applications2-3 envisaged for MoS2 with different morphologies, particle sizes and porous features has been significantly extended, including electronic devices, optical devices, mechanical devices, thermoelectric devices, photocatalysis, catalyst synthesis, organic synthesis, hydrodesulfurization, hydrodeoxygenation, hydrogen evolution reaction (HER), energy storage and conversion, gas sensors and electrolysis of water.4 Moreover, MoS2 has been reported to show biocompatibility5-6 and biodegradability7, making it suitable for biomedical applications. For examples, it has been demonstrated that chemically exfoliated MoS2 sheets could induce both membrane and oxidation stress, acting as antibacterial agents8. In the field of photothermal therapy and drug delivery, multifunctional theranostic nanocomposite based on MoS2 for multimodal imaging-guided cancer therapy has achieved effective tumor ablation in an animal tumor model9-10. Since MoS2 nanomaterial is a semiconductor, it is also ideal for applications as biosensors for bio-detection systems involving protein and DNA11-12. Furthermore, MoS2 can be applied as a valid contrast agent in X-ray computed tomography imaging due to the considerable X-ray adsorption properties of the Mo element in the compound10. Despite all the advances made on the MoS2 nanomaterial, many knowledge gaps remain for a thorough understanding how it 2 ACS Paragon Plus Environment

Page 2 of 26

Page 3 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interacts with biological systems, which requires powerful tool such as molecular dynamics (MD) simulations for large scale computer modellings to complement state-of-the-art experimental techniques. Therefore, an accurate computational model of MoS2 compatible with biological systems is highly desirable in order to garner mechanistic insights into the organization of molecules and commonly employed solvents Early parameterizations on MoS2 mainly focused on its metal adsorption properties13, tri-layer formation dynamics14, thermal transport properties15 or mechanical properties such as sliding friction16. However, these parameters generally cannot describe the interactions between MoS2 and water when water models are combined in MD simulations. To overcome this, it is necessary to make modifications based upon previous works to match the experimental work of adhesion (Wad) of a water model. The experimental way of obtaining Wad of water on MoS2 surface is to measure the WCA on the material surface and apply the Young-Dupré equation (Equation 1): 𝑊𝑎𝑑 = 𝛾𝐿𝑉(cos 𝜃 + 1)

(1)

where γLV is the water-vapor surface tension, θ is the WCA. An early study found that the WCA on MoS2 basal plane (consisted of purely sulfur atoms) was ~98° for a monolayer MoS2 on silica17, indicating the MoS2 surface is quite hydrophobic. However, a more recent study revealed that an aging (or polluting) process can change the WCA of MoS2 surface dramatically, increasing it from 69.0° ± 3.8° on the fresh basal plane of a bulk MoS2, to about 90° after oneday exposure to the ambient air due to the hydrocarbon contamination.18 Such an aging process was independently discovered by another work19 where WCA on the fresh bulk MoS2 was determined to be 67.2°. Therefore, a fresh MoS2 basal surface is drastically different from an aged one in terms of hydrophobicity. Such discrepancy in experiments inevitably led to a wide range of available parameters for MoS2 in computational studies. For example, our past parameterization study20 focused on the WCA on the aged surface of MoS2, while some other studies21-23 focused on the fresh surface of MoS2. In order to address the wide selection of parameters, in this study we aim to account for both fresh and aged MoS2 surfaces with a coherent set of parameters. The methodology of parameterizing MoS2 against WCAs includes directly fitting the high level quantum predictions23 and fitting the experimental WCAs by measuring simulated water droplet20-22, 24. However, parameters fitted against quantum methods are not readily compatible 3 ACS Paragon Plus Environment


with the commonly used water models such as TIP3P and SPC/E. As for the WCA measured from simulated nanoscale water droplet (in contrast to much larger droplets in experiment), the size-dependency can significantly affect the accuracy. To directly benchmark against water models and to avoid the size-dependency in the parametrization process, we developed an approach to estimate Wad using the free energy perturbation (FEP) method. Furthermore, in previous parameterization studies21-23 on fresh MoS2 surface, the SPC/E water model25 was used for material simulations. Because of the popularity of the TIP3P water model26 used in biological simulations along with many modern force fields27-29 and potential usages9, 3039

of MoS2 in biological systems, a compatible force field with TIP3P water model is highly

demanded. Therefore, in addition to the SPC/E water model, we provide a separate set of parameters suited for the TIP3P water model, which also demonstrates the transferability of this method.

Method Molecular dynamics simulations Our simulated systems consisted of a trilayer MoS2 and a water box with 2866 water molecules (Figure 1). A trilayer was used in our current simulations because once going beyond three layers, the interaction from the bottom layers to the water molecules absorbed on the top layer is minimal. The MoS2 thickness varies depending on the conditions and synthesizing methods.40-42 In this study we used the distances from our previous study20: the distance between any Mo and S honeycomb sheets in a monolayer was 1.557 Å and the distance between any two neighboring layers of Mo atoms was 6.15 Å. All atoms in the MoS2 nanomaterial were restrained (frozen) in the simulations. Periodic conditions were used in all three dimensions. The simulation box was set to 150×43.9×41.7 Å3 which provided sufficient room (~75 Å in distance) between the periodic image of water film and MoS2 surface. All simulations were performed using NAMD2.11 MD package43. The partial charges of MoS2 were adopted from a previous study15 where Mo has a partial charge of +0.76e and S -0.38e (e is the elementary charge). For Lennard-Jones (LJ) interactions, we kept Rmin (= 21/6σ) and εMo the 4 ACS Paragon Plus Environment

Page 4 of 26

Page 5 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


same from our previous study (Table 1) and only optimized εS to best fit the experimental WCAs. Thus, in contrast to the all-atom approach for an aged (or contaminated) surface44-45, we adopted a coarse-grained approach by adjusting εS to include changes to WCA due to the adsorbed hydrocarbons. Instead of obtaining force-field parameters for each system with hydrocarbons and validating them using contact angle simulations21-23, we directly fit force-field parameters to reproduce the experimental contact angles. The long-range electrostatic interactions were calculated with the Particle Mesh Ewald (PME) method with a grid size of 100×40×40, while the VDW interactions were treated with the smooth cutoff with a switching distance of 10 Å and a final cutoff of 12 Å. Following similar protocols in our previous studies4652,

we first performed 10,000 steps of steepest descendant energy minimization, followed by a 20

ns equilibration in the canonical ensemble (NVT) at 300 K before the FEP simulations for any parameter optimization (see below). All bonds were treated as rigid ones thereafter to allow a 2 fs timestep.

Free Energy Perturbation method The FEP method was first introduced by Zwanzig to investigate the high-temperature equation of state for noble gases53. It was then used for methane hydration free energy54, water cavity formation free energy55, methanol hydration free energy56, and more recently for many complex biological systems57-69 . It can be formulated as: Δ𝐴𝜆 = 𝐴𝜆 + Δ𝜆 ― 𝐴𝜆 = ― 𝛽 ―1ln 〈exp [ ―𝛽(𝐸𝜆 + Δ𝜆 ― 𝐸𝜆)]〉𝜆

(2)

Δ𝐴 = ∑𝜆Δ𝐴𝜆

(3)

where ΔA is the Helmholtz free energy difference due to the perturbational change, λ is the FEP step parameter of choice that tunes the non-bonded interactions between MoS2 and the water in our simulations, β is the Boltzmann constant, and Eλ is the total energy of the simulated system at step λ. Specifically, Eλ can be defined as: 𝐸𝜆 = (1 ― 𝜆)𝐸0 + 𝜆𝐸1

(4)

where E0 is the energy when λ = 0, referring to the energy of water film alone (no contact with MoS2) in our simulations, and E1 is the energy when λ = 1, referring to the total energy of the water film in contact with MoS2 (Figure 1). The step parameters (λ) for FEP simulations were set 5 ACS Paragon Plus Environment


Page 6 of 26

to be (0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00). Soft-core VDW interactions were used following the previous work64. Electrostatic interactions were turned on after λ ≥ 0.1. Each FEP simulation lasted at least 16 ns. To prevent the water film from diffusing through the periodic boundary and interacting with MoS2’s periodic image, we put restraints on the water film (see Supporting Information for more detailed description on the restraints). The Helmholtz free energy change ΔA was then used to calculate the work of adhesion: Δ𝐴𝑎𝑑 = Δ𝐴𝑏𝑜𝑢𝑛𝑑 ― Δ𝐴𝑓𝑟𝑒𝑒 = 𝑊𝑎𝑑 × 𝑆

(5)

where ΔAbound (bound state) and ΔAfree (free state) are from two sets of FEP simulations defined in the thermodynamic cycle of Figure 2A, and S is the area of water-MoS2 interface (or the crosssection area of the simulation box, Figure 1). Similarly, we applied the FEP method to calculate the water-vapor surface tension γLV using another thermodynamic cycle as illustrated in Figure 2B. Finally, we estimated the WCAs (θ) using the relationship expressed in Equation 1 and compared with experimental values.

Results and Discussions Water-vapor surface tension As a first step of validation, we calculated the water-vapor surface tension γLV with our FEP based method. The water-water contact free energy ΔAw can be computed with the following equation: Δ𝐴𝑤 = 2𝛾𝐿𝑉 = Δ𝐴𝑏𝑜𝑢𝑛𝑑 ― Δ𝐴𝑓𝑟𝑒𝑒

(6)

where ΔAbound (bound state) and ΔAfree (free state) are from two sets of FEP simulations defined in the thermodynamic cycle of Figure 2B. ΔAbound and ΔAfree are calculated with free energy perturbation (FEP) simulations here. We use TIP3P water model as a test case to find the optimized simulation setup (shown in Figure S1A). The previously calculated γLV of TIP3P water was 52.3 ± 1.5 mJ/m2.70


Page 7 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To avoid the “end-point catastrophe”, we introduced wall-potential restraints on the water films, shown in Figure S1A. We examined how tuning restraints on the water films affected the estimation of water-vapor surface tension. The details are shown in the Supporting Information (SI text, Figure S1, Table S1 and Table S2). The optimized condition yielded 53.8 ± 1.0 mJ/m2 as the γLV value at 300 K for the TIP3P water model, in agreement with the literature.70 With the same approach, γLV of the SPC/E water at 300 K was computed to be 64.5 ± 0.2 mJ/m2 (which also agreed well with the 63.6 ± 1.5 mJ/m2 from the literature70). In addition, we compared γLV of TIP3P and SPC/E water models estimated from our current FEP method to the ones with Kirkwood-Buff method71-72: 𝐿 /2

𝛾𝐿𝑉 = ∫ ―𝑧 𝐿 /2(𝑝𝑁(𝑧) ― 𝑝𝑇(𝑧))d𝑧 𝑧

where Lz is the size of the simulation box on the z direction, pN is the pressure tensor component on the normal direction of the interface and pT is the pressure tensor component on the tangential direction of the interface. Following this Kirkwood-Buff method, 10 ns MD simulations were conducted for a water box (with the size of 20×20×40 Å3) of TIP3P water and SPC/E water, respectively. The pressure tensor was calculated using the pressure profile function of NAMD2 Furthermore, following Vega et al.70, a long range tail correction (LRC) has been applied to the Kirkwood-Buff integration of surface tension using the following equation, derived by Blokhuis et al.73: 1

∞

𝛾𝑡𝑎𝑖𝑙 = 12𝜋𝜀𝜎6(𝜌𝐿 ― 𝜌𝑉)2∫0𝑑𝑠∫𝑟 𝑑𝑟 coth(𝑟𝑠/𝑑0)(3𝑠3 ― 𝑠)/𝑟3 𝑐

where ε and σ are VDW parameters of the water model, rc is the cutoff of the VDW interactions in the simulation (~11 Å in our simulations, which is the average between the 10 Å switching function and the 12 Å cutoff). ρL, ρV and d0 can be obtained by fitting the density of the water box to the following 4-parameter equation (with an extra z0), assuming hyperbolic tangent function of the density at the water-vapor interface: 1

1

𝜌(𝑧) = 2(𝜌𝐿 + 𝜌𝑉) ― 2(𝜌𝐿 ― 𝜌𝑉)tanh [(𝑧 ― 𝑧0)/𝑑0]. Water surface tension was then reported in Table S3. Error estimations are based on a block average over 10 blocks within each of the simulations. The obtained water surface tension values 7 ACS Paragon Plus Environment


for TIP3P and SPC/E water models agree very well with the previously calculated values from Vega et al.70 This in turn agrees with the water surface tension values calculated with our current FEP method. Interestingly, no tail correction seems to be needed for the FEP method to agree with the Kirkwood-Buff method, which was also found in a previous free energy calculation method using Bennett acceptance ratio74.

FEP simulations as a viable method to estimate WCAs As a free energy calculation method, the dry-surface method was successfully adopted to estimate WCAs of materials in previous studies23, 75. Here we show that with certain restraints, the FEP method76 which is commonly applied in biological systems57-69 can be readily applied to calculations of WCAs on materials. More specifically, we used the same configuration of wallpotential restraints shown in the water-vapor surface tension calculations to avoid the “end-point catastrophe” (see Supporting Information for details). We then calculated the WCAs using Luan et al.’s parameters20 for TIP3P water model and Heiranian et al.’s parameters23 for SPC/E water model as a verification. Our estimated WCA with Luan et al.’s parameters and the TIP3P water is 89.7° ± 0.3° (90.1° in the literature20 with the largest water droplet). Meanwhile, our estimated WCA with Heiranian et al.’s parameters and the SPC/E water is 73.6° ± 0.2° (~ 75° in their work23). These results indicate that our current approach based on the FEP method can successfully reproduce the water-vapor surface tension γLV and WCAs from previous literatures.

Water models significantly influence WCAs on MoS2 Next, we turned to the experimental WCAs for the optimization of MoS2 parameters. We performed a total of 90 FEP calculations for both TIP3P water and SPC/E water with a series of εS (the accumulated simulation time was 1.44 μs). The calculated work of adhesion and WCAs are listed in Table 1 (for TIP3P water model) and Table 2 (for SPC/E water model). As shown in Figure 3, the most striking feature is the linear dependency between WCA (θ) and εS for both the TIP3P water model (Figure 3A) and SPC/E water model (Figure 3B), with the coefficient of determination R2 ~ 1, respectively. This is in contrast to the non-linear correlation between θ and εMo in our previous study20 (more discussions below and in Appendix). 8 ACS Paragon Plus Environment

Page 8 of 26

Page 9 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The experimentally measured WCAs on the fresh MoS2 basal surface (bulk) was determined to be 69.0° by Kozbial et al.18 and 67.2° by Gurarslan et al.19 Here we aim for 68°, which results in a value of -0.40 kcal/mol for εS using TIP3P water model. This is within the εS range of -0.38 ~ 0.47 kcal/mol for organic sulfur from CHARMM force field. In contrast, our recommended value of εS is -0.60 kcal/mol using SPC/E water model, compared to -0.46 ~ -0.50 kcal/mol provided in the literatures for the same water model22-23. This drastic difference results from three factors. (i) VDW interactions dominate the MoS2-water binding free energy; (ii) the VDW attraction of SPC/E (756.85 kcal1/2·Å6/mol1/2, defined as

― 𝜀𝑂 ⋅ 𝑅6𝑚𝑖𝑛)77 is weaker than TIP3P

(762.84 kcal1/2·Å6/mol1/2)78 (therefore, when εS is the same, an overall lower work of adhesion is expected from SPC/E water model, in agreement with the data from Table 1 and Table 2). (iii) The water-vapor surface tension of SPC/E water model is larger than TIP3P water model. Combining all these three causes, one needs a higher εS to compensate the loss of VDW interactions when using the SPC/E water model. The aged MoS2 basal plane (a simplified description including the effect of adsorbed airborne hydrocarbons44), however, has a wide range of WCAs (73° ~ 98°) recorded in the literatures17-18. Because of this, we suggest picking the value of εS that fits a specific experiment of interest with our “linear model”: 𝜃(TIP3P) = 145.7( ± 1.0)𝜀𝑆 + 126.3( ± 0.4) 𝜃(SPC/E) = 99.0( ± 1.5)𝜀𝑆 +127.2( ± 0.7)

(7)

where θ is in degrees and εS is in kcal/mol. See Appendix for the proper range of θ and εS where these equations apply. In the next section, we will further investigate the correlation between WCA and εS from our simulations.

VDW interactions dominate MoS2-water interactions To analyze the contributions from each energy terms, we independently performed 20 ns of equilibrated simulations with the recommended εS values for fresh basal plane of MoS2 from our current work. For consistency, all other conditions (number of atoms, temperature, non-bonded interactions, timestep, etc.) are kept the same from FEP simulations. The energy decomposition results are given in Figure 4 and Table S4 for both water models. 9 ACS Paragon Plus Environment


As shown in Table S4, the total VDW energy (normalized by the interfacial area of the simulation box, same hereafter) is -0.149 ± 0.004 kcal/mol/Å2, while the electrostatic energy (calculated with PME) is only -0.0016 ± 0.0007 kcal/mol/Å2 for the TIP3P model. As for the SPC/E model, the VDW energy is -0.177 ± 0.004 kcal/mol/Å2 and the electrostatic energy is 0.0017 ± 0.0005 kcal/mol/Å2. Clearly, the VDW interaction is stronger between MoS2 and SPC/E water than that between MoS2 and TIP3P water model. This is due to the much deeper Lennard Jones potential well (i.e., larger well-depth εS) required for the SPC/E water model (0.60 kcal/mol) than the TIP3P water model (0.40 kcal/mol). These data clearly show that the VDW interaction is the dominating term in both cases, in agreement with the previous experimental18 and theoretical21 study. To further illustrate the VDW contributions from each layer of MoS2, we plotted the VDW energy between the water film and (1) each honeycomb sheet of Mo or S in a monolayer MoS2, and (2) each layer of the trilayer MoS2 in Figure 4A (TIP3P) and Figure 4C (SPC/E). Note that we define the “upper sheet” as the honeycomb sheet of S that is closest to the water film within a monolayer MoS2. The other honeycomb sheet of S is defined as the “lower sheet”. We define the first layer of the trilayer MoS2 as “the top layer” with the other two as “the middle layer” and “the bottom layer”, respectively. Examining Figure 4A and Figure 4C, we see that the major contribution of VDW energy comes from the upper sheet of S in the top layer of MoS2. This is in agreement with the fact that tuning εS affects much more profoundly than tuning εMo. For example, we found that changing εMo by 0.05 kcal/mol leads to about 2° change on the WCA.20 In this work, changing εS by 0.05 kcal/mol leads to about 5° or 7° change on the WCA (Table 1 and Table 2). Aside from the top layer of MoS2, the VDW contributions from the middle layer is a mere 0.0013 ± 0.0001 kcal/mol/Å2 and 0.0016 ± 0.0001 kcal/mol/Å2 for TIP3P and SPC/E water models, respectively. Meanwhile, basically zero contribution arises from the VDW interaction between the bottom layer of MoS2 and the water box, which is not surprising since the thickness of a monolayer MoS2 is 6.15 Å, while the cut-off distance for VDW interactions is justifiable at 12 Å due to the diminishing effect. Thus, after two layers of MoS2 essentially no VDW interactions from the Mo or S honeycomb sheet will be felt by the water molecules.19 On the other hand, the contribution from the long-range electrostatic interactions is surprisingly negligible as compared to that from the VDW interactions. The entire electrostatic energy comes 10 ACS Paragon Plus Environment

Page 10 of 26

Page 11 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the top layer of MoS2 (Table S4, Figure 4B and Figure 4D), while the middle and bottom layers of MoS2 have no meaningful impact on the total electrostatic energy of the system. This is understandable from a simple physical perspective if we imagine an “ideal” situation, where a “multipole” (as the simplified water film) is interacting with alternating positively charged and negatively charged, infinitely-large planes (e.g., Mo and S honeycomb sheets through periodic boundary condition). Assuming that all charges are uniformly distributed; all planes are equally separated; and the total net charge is zero, from electromagnetic theories, the electric fields out of a uniformly charged and infinitely large surface are constant anywhere above the surface. Thus, regardless of the distance between the multipole and MoS2 layers, as long as the multipole does not penetrate the MoS2 layer, the electrostatic forces should always cancel each other out from the three honeycomb sheets of Mo and S that form the monolayer MoS2 and the total should be zero. Therefore, the electrostatic energy should be zero too. In reality (from our simulation), the small electrostatic energy (~-0.0016 kcal/mol/Å2) can be simply attributed to the discrete nature of the charged atoms on each MoS2 layer (particularly the closest one). This negligible electrostatic energy contribution also explains the fact that little effect comes from tuning partial charges of MoS2 when fitting experimental WCAs, and why a wide range of partial charges exist in the literatures22. In contrast, if we split the electrostatic contributions of MoS2 layer to each honeycomb sheet of Mo or S, the resulting energy is substantial (assuming zero potential on each honeycomb sheet) as one would expect (Table S4, Figure 4B, Figure 4D). The absolute electrostatic energy of the single Mo/S honeycomb sheet from the bottom MoS2 layer is the highest because the zero-point energy reference point is the charged plane itself. As we can see, the major contributing term to the interaction energy between water and MoS2 is the VDW interaction between water and the top layer of MoS2, particularly, the upper sheet of S. This is simply due to the unique structure of MoS2 where the inter-nanosheet distance is relatively large and the nearby electrostatic field is 0. Based on these, two conclusions can be inferred. First, a linear correlation can be found between WCAs and εS under certain circumstances. In the appendix, we derive such correlation with the assumption of the dominating VDW term between the upper sheet of S and water (denoted as EVDW-STU). We also reveal other possible contributing terms for WCAs, such as the entropy from water 11 ACS Paragon Plus Environment


configurations. Second,

― 𝜀𝑆 differs almost ~20% in order to obtain the same WCAs for TIP3P

and SPC/E water models. This is simply due to the fact that water surface tension of these two models differ by 19.9% based on our calculations. Meanwhile, TIP3P and SPC/E water models share nearly identical VDW parameters (with thus, in order to obtain a similar ~68° WCA,

― 𝜀𝑂 differed by 1.0% but compensating εH), ― 𝜀𝑆 has to differ by 22.5% when using with the

two different water models.

Number of MoS2 layers does not significantly affect WCAs To examine the effect of number of MoS2 layers on the WCAs, we then compared the WCAs from the trilayer MoS2 and a monolayer MoS2. The results are shown in Table 3. The WCA is 68.1° ± 0.2° for the trilayer MoS2, and 68.9° ± 0.3° for the monolayer MoS2. The very small difference (consistent with a previous study79) is not surprising, because as shown above the total contribution of the middle and bottom layers to the interaction energy is minimal, only about 0.9% (Table S4). Finally, it might be worthwhile to point out that in reality the base material could be charged, i.e., the net charge is not zero. As discussed above, the total electrostatic energy is only related to the relative distances between multipoles in the water film, and it is irrelevant to the distance between the base material and the water film (db-w). For example, a dipole with charges (-q and +q) and a dipole moment (q·d) has an electrostatic energy of Eqd in an ideal constant electric field (E), regardless of the distance between the dipole and the source of the field (an infinitely large charged plane). However, the variance of the electrostatic energy (which contributes to the entropy, see Equation 9 in Appendix) depends on the magnitude of the electrostatic potential, which is determined by db-w. By adjusting the number of layers of MoS2, db-w varies. Therefore, number of layers of MoS2 may affect WCAs in the experiment80 when the base material is charged or contaminated.

Conclusion


Page 12 of 26

Page 13 of 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We revisited and optimized the parameterization of MoS2 using a wide range of experimental WCAs in this work. FEP simulations are proved to be suitable for such process. The obtained values of γLV are in reasonable agreement with previous studies for both TIP3P and SPC/E water models. Based on the values of γLV of these two widely used water models, we provided corresponding values of LJ parameter of S (εS) for MoS2, listed in Table 1 and Table 2. For uses with biological systems, TIP3P water model is more prevalent because it was developed along with the popular CHARMM and OPLSAA force fields (also the default water model). Note that in recent studies, water models such as TIP4P/200522, 81, TIP4P/ICE22, 76, 82 and TIP4P-D83 have seen an increased popularity due to their superior properties. Compatible MoS2 parameters with these new water models will be provided in future studies. A surprisingly linear correlation was found between WCA and εS, regardless of the water model, which was shown to be related to the dominant VDW contribution from the nearest honeycomb sheet of S to the total adhesion of work (Wad). The electrostatic interactions, on the other hand, do not affect the total Wad much due to the symmetric nature of the charge distributions on the MoS2 layers, which minimizes the effort of tuning partial charges of MoS2 in fitting experimental WCAs. Additionally, the number of layers of MoS2 do not influence on WCA either, unless the base material has net charges (due to contaminations, etc.). Ideally, an all-atom model of the aged MoS2 surface, including an explicit treatment of hydrocarbon contaminants (as evidenced in recent experiments18 and models44-45), is required to capture kinetic properties (e.g. the diffusion constant) and detailed water structure near the surface. Our simplified (unifying adsorbed contaminants with the MoS2) model, however, provides correct surface energies on average for various aged surfaces. We expect this new FEP based approach that accurately calculates the work of adhesion and WCAs to help parametrize other novel nanomaterials as well.

Appendix The linear correlation between WCA and εS Combining Equation 2, Equation 3 and Equation 5, the adhesion of work (binding free energy of the water film on MoS2) is formulated as:



∑ ln 〈exp [ ―𝛽(𝐸

𝑆𝑊𝑎𝑑 = Δ𝐴𝑎𝑑≅Δ𝐴𝑏𝑜𝑢𝑛𝑑 = ― 𝛽 ―1

𝜆 + Δ𝜆

𝜆

Page 14 of 26

― 𝐸𝜆)]〉𝜆because ΔA

free

~ 0 if

intramolecular energy calculations are decoupled from FEP calculations (alchDecouple = on in NAMD2.11). Throwing in Equation 4, we get:

∑ ln 〈exp [ ―𝛽Δ𝜆(𝐸

𝑆𝑊𝑎𝑑 = ― 𝛽 ―1

1

𝜆

― 𝐸0)]〉𝜆Assume a linear correlation between ΔA

bound

and

ΔAλ, we can simplify the above equation to: 𝑆𝑊𝑎𝑑 = ― 𝛽 ―1ln 〈exp [ ―𝛽(𝐸1 ― 𝐸0)]〉

(8)

where E1-E0 purely comes from non-bonded interactions between MoS2 and water: 𝐸1 ― 𝐸0 = 𝐸VDW(MoS2water) + 𝐸elec(MoS2water) =

∑𝐸

VDW(𝑙𝑎𝑦𝑒𝑟water)

+

𝑙𝑎𝑦𝑒𝑟𝑠

∑𝐸 𝑙𝑎𝑦𝑒𝑟𝑠

elec(𝑙𝑎𝑦𝑒𝑟water)where

Eother is the sum over non-bonded

= 𝐸VDW ― STU + 𝐸other interactions except for the VDW contribution from the S (upper) of the top layer of MoS2. Based on the energy contribution analysis, Eother

Parameterization of Molybdenum Disulfide Interacting with Water

Recommend Documents