Interpretable molecular models for molybdenum disulfide and insight into selective peptide recognition

Abstract

Molybdenum disulfide (MoS₂) is a layered material with outstanding electrical and optical properties. Numerous studies evaluate the performance in sensors, catalysts, batteries, and composites that can benefit from guidance by simulations in all-atom resolution. However, molecular simulations remain difficult due to lack of reliable models. We introduce an interpretable force field for MoS₂ with record performance that reproduces structural, interfacial, and mechanical properties in 0.1% to 5% agreement with experiments. The model overcomes structural instability, deviations in interfacial and mechanical properties by several 100%, and empirical fitting protocols in earlier models. It is compatible with several force fields for molecular dynamics simulation, including the interface force field (IFF), CVFF, DREIDING, PCFF, COMPASS, CHARMM, AMBER, and OPLS-AA. The parameters capture polar covalent bonding, X-ray structure, cleavage energy, infrared spectra, bending stability, bulk modulus, Young's modulus, and contact angles with polar and nonpolar solvents. We utilized the models to uncover the binding mechanism of peptides to the MoS₂ basal plane. The binding strength of several 7mer and 8mer peptides scales linearly with surface contact and replacement of surface-bound water molecules, and is tunable in a wide range from −86 to −6 kcal mol⁻¹. The binding selectivity is multifactorial, including major contributions by van-der-Waals coordination and charge matching of certain side groups, orientation of hydrophilic side chains towards water, and conformation flexibility. We explain the relative attraction and role of the 20 amino acids using computational and experimental data. The force field can be used to screen and interpret the assembly of MoS₂-based nanomaterials and electrolyte interfaces up to a billion atoms with high accuracy, including multiscale simulations from the quantum scale to the microscale.

---

Introduction

Two-dimensional (2D) materials, including graphene, hexagonal boron nitride, and transition metal dichalcogenides (TMDs or TMDCs) such as molybdenum disulfide (MoS₂) have received widespread attention in recent years due to their unique structural, electronic, and conductive properties.^1–3 The compounds have strong polar covalent bonds within the layers and weaker inter-layer interaction. Specifically, MoS₂ is widely used in electrochemical catalysts,^4,5 sensors,⁶ composites,^7,8 and photovoltaics.^9,10 The surface, interfacial, and mechanical properties are important for the design of functional materials and increased control over performance. For example, the surface properties of MoS₂ influence the performance as a catalyst in the hydrogen evolution reaction,¹¹ the conductivity upon specific binding of analytes,¹² as well as the performance in electrode materials and polymer/MoS₂ composites.^13,14

To-date, however, it has remained challenging to understand and predict interfacial and mechanical properties of MoS₂ and related 2D materials to support the rational synthesis and performance in applications. Electronic structure calculations such as density functional theory (DFT), coupled with experimental data, have yielded critical information on the electron density and band structures.¹⁵ However, DFT calculations are limited to small scales of hundreds of atoms, insignificant dynamics, and the uncertainties in property predictions for transition metals and their compounds such as MoS₂ are high.^16–18 For example, computed surface energies^19–21 and elastic constants deviate more than 50% relative to experiment (Table 1).²² Binding energies of small molecules on heavy metal surfaces can have several 100% error depending on the density functional and available corrections.^16–18

Table 1 Comparison of key properties of 2H-MoS₂ according to experiments and computational methods. From left to right: method, lattice parameters, Mo–S bond distance (r₀), cleavage energy of the (0001) basal plane (E_cleav), major IR/Raman peak, in-plane Young’ modulus (Y), bulk modulus (B), contact angle with liquids at 298 K, atom mobility, compatibility and transferability of the potential

Method	Lattice parameter		r 0 (Å)	E cleav (mJ m^−2)	IR/Raman^a (cm^−1)	Y (GPa)	B (GPa)	Contact angle (°)		All atoms mobile	Applicable potentials and compatibility^b
	a (Å)	c (Å)						H2O	CH2I2
a Major IR and Raman peak (Fig. 3c). b Details of the potential function are given in Table S1 in the ESI. c We estimate the cleavage energy (surface energy) as 150 ± 10 mJ m^−2, close to the experimental value 121 mJ m^−2 for natural bulk MoS2 samples (ref. 54) and to the lower end of values from DFT calculations of 160 mJ m^−2 (ref. 19). d Reactive empirical bond-order (REBO) potential. e The energy function comprises multiple exponential terms.
Expt	3.16 ^55	12.3 ^55	2.37 ^55	45–121 ^57	384 ^58	170–300 ^59–62	43 ± 3 ^63	69 ± 4 ^52	15 ± 2 ^52	Yes	NA
This work	3.16	12.25	2.38	150 ± 2^c	360 ± 5	175 ± 2	41.7 ± 1	69 ± 2	15 ± 2	Yes	IFF, CVFF, PCFF, OPLS-AA, Dreiding, CHARMM, AMBER
DFT	3.19	12.45 ^64	2.43 ^64	160–284 ^19–21	387 ^65	220 ^66	68 ^22	NA	NA	Yes	NA
SR(2017)^49	3.16	12.16	NA	241	160	226	43.7	69.6	0	Yes	OPLS
LU(2016)^50	NA	NA	NA	NA	NA	NA	NA	86–90	NA	Yes	OPLS
JI(2013)^47	3.21	11.25 ^67	2.39 ^67	NA	NA	229	NA	NA	NA	Yes	Stillinger–Weber (SW)
DA(2012)^46	3.20	11.59 ^67	2.42 ^67	NA	NA	NA	NA	NA	NA	Yes	QEq
VA9(2010)^45	3.25	12.23 ^67	2.45 ^67	NA	404	NA	NA	0 ^49	NA	No	CVFF
VA8(2010)^45	3.21	12.27 ^67	2.41 ^67	NA	584	NA	NA	0 ^49	NA	No	CVFF
LI(2009)^48	3.16	12.31 ^67	2.45 ^67	NA	NA	NA	NA	97.9 ^49	NA	Yes	REBO^d
ON(2009)^44	3.16	12.45 ^67	2.40 ^67	NA	NA	NA	NA	NA	NA	No	Buckingham
MO(2008)^43	3.05	13.89 ^67	2.34 ^67	NA	NA	NA	NA	69.4 ^49	NA	No	Buckingham
BE(2003)^42	3.16	12.31 ^67	2.43 ^67	NA	NA	NA	NA	0 ^49	NA	No	Embedded^e
FA(1996)^41	3.16	12.05 ^67	2.42 ^67	NA	NA	NA	NA	NA	NA	Yes	Embedded empirical^e
BR(1992)^40	3.17	12.28 ^67	2.42 ^67	NA	NA	NA	NA	NA	NA	No	Empirical without Coulomb^e
DR(1988)^39	3.19	11.74 ^67	2.42 ^67	NA	NA	NA	NA	NA	NA	Yes	Empirical without Coulomb^e

---

Molecular dynamics (MD) simulation has become a reliable tool to gain insight into the details of molecular motion, nanoscale assembly, crystal growth, and catalysis. Examples include peptide–metal binding,^23–26 silica–electrolyte interactions,^27,28 mineralization of apatite,²⁹ nanoscale forces in building materials³⁰ and polymer–carbon nanomaterials composites.^31–34 Interpretability of the force field parameters is hereby essential to explain and predict the physical and chemical behavior. Parameters that are consistently derived and explained, rather than numerically fitted, achieve multiple times increased reliability and compatibility as shown by the interface force field (IFF) and its surface model database.^35–38

In this study, we introduce parameters for MoS₂ with record accuracy (Fig. 1) including validation and application to explain specific binding of peptides and amino acid residues. Interestingly, a number of modeling studies since the 1980s have attempted to characterize the interfacial and surface properties of 2H-MoS₂ (Table 1).^39–50 However, the results are poor since earlier force fields consist of empirical fitting parameters that are missing a chemical interpretation. Available potentials neglect the polarity of Mo–S bonds and use energy expressions other than harmonic, Coulomb, and Lennard-Jones, that are difficult to use with common models for water, minerals, organic compounds, and biomolecules.⁵¹ Due to the lack of interpretability and rationale, the parameters by DR (1988),³⁹ BR (1992),⁴⁰ FA (1996),⁴¹ BE (2003),⁴² MO (2008),⁴³ LI (2009),⁴⁸ VA9 (2010),⁴⁵ DA (2012),⁴⁶ and JI (2013)⁴⁷ have considerable deviations (>2%) in Mo–S bond length and lattice parameters in comparison to X-ray diffraction data for 2H-MoS₂ (Table 1). VA8 (2010)⁴⁵ and VA9 (2010)⁴⁵ parameters were tested to calculate vibration spectra, whereby VA9 (2010)⁴⁵ achieves good agreement in vibration frequencies. Nevertheless, a discussion of the polarity of chemical bonding and polarizability was not included, leading to a computed water contact angle of 0° that deviates 100% from measurements of 69° on freshly cleaved MoS₂ surfaces.⁵² Similar difficulties to reproduce surface and interfacial energies, or contact angles, are seen for LI (2009),⁴⁸ BE (2003),⁴² and LU (2016).⁵⁰ The SR (2017)⁴⁹ parameters for 2H-MoS₂ report better performance in lattice parameters, surface wetting, and some mechanical properties. However, bond angles are oversimplified to be octahedral and deviate more than 30° from X-ray data. Vibration constants lack a justification and are over 50% (200 cm⁻¹) redshifted relative to the experimental IR spectrum, allowing lipid-like buckling of the 2D MoS₂ layers, which are known to be stiff structural reinforcements in experiments (Fig. 2). Lennard-Jones parameters have likewise no documentation and include well-depths for sulfur that are about twice as high as typical values in similar chemical environments such as sulfides or sulfoxides. As a result, the parameters offset the cleavage energy more than +60% from a realistic value and yield a contact angle of diiodomethane of 0° in simulations, which differs 100% from the value of 15° obtained in experiments (Table 1). In summary, lack of chemical understanding and interpretation of prior potentials cause multiple discrepancies, disables compatibility and meaningful property predictions.

Fig. 1 Model of a 2H-MoS₂ unit cell (from ref. 55). The lattice parameters are a = b = 3.16 Å and c = 12.295 Å under standard conditions. The atomic charges for Mo and S are +0.5e and −0.25e, respectively. We describe the structure using two atom types for Mo and 4 atom types for S, which are necessary to identify unique bonds and angles. We distinguish one Mo–S bond of 2.37 Å length and five bond angles θ₁, θ₂, θ₃, θ₄, θ₅, of which the θ₂ angle (∠S1–Mo1–S3) can assume two values θ_2a = 78.38° and θ_2b = 134.4° and has no constraints in the model.

Fig. 2 Comparison of the performance of prior parameters and new parameters for the structure of bulk and single-layer MoS₂. 3D structure of MoS₂ after 50 ns NPT simulation at 298.15 K and 1.013 MPa with (a) prior parameters (SR 2017, ref. 49) and (b) new parameters. The internal structure of the prior model collapses while it is maintained very close to X-ray data in the new model. (c) An AFM image of single-layer MoS₂ with a thickness of ∼0.8 nm, showing a planar stiff geometry (reproduced with permission from ref. 56). (d and e) Snapshots of an MoS₂ layer after 0.5 ns NVT simulation at 298.15 K in vacuum with prior parameters (SR 2017, ref. 49) and new parameters, respectively. The old parameters lead to unphysical lipid-like buckling, whereas the new model preserves the shape and allows minor undulations.

Step-by-step, we resolved these shortcomings by the logical development and explanation of a Hamiltonian, applicable to other 2D materials, and demonstrate improvements in reliability by more than 10 times compared to earlier models for multiple properties.^37,53,54 Our focus is on reproducing chemical bonding, lattice parameters, surface energies for 2H-MoS₂, as well as mechanical properties with a minimum of random parameter search and fitting. We explain compatibility with at least 7 force fields including the interface force field (IFF), Consistent Valence Force Field (CVFF), Dreiding, the Polymer Consistent Force Field (PCFF), COMPASS, Chemistry of Harvard Macromolecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), All-Atom Optimized Potentials for Liquid Simulations (OPLS-AA), and others.³⁷ As an example of an application, we explain the affinity of various peptides to the basal MoS₂ surface in agreement with experimental observations and elucidate molecular controls for tuning the binding strength.

The outline of this paper is as follows. Following a brief introduction of the interatomic potential in different formats (IFF, CVFF, PCFF, Dreiding, CHARMM, AMBER and OPLS-AA), we describe the chemistry of MoS₂ and the transcription into interpretable force field parameters. We discuss, step-by-step, critical information about bonding, structure, vibrations, and interatomic interactions to derive a consistent classical Hamiltonian that reproduces structures and energies. We identify the origin of errors in prior models and explain how limitations are overcome, extensible to any other layered materials. Then, we discuss the validation of structural, vibration, interfacial, and mechanical properties. As an application, we elucidate selective peptide binding mechanisms to 2H-MoS₂. The manuscript ends with conclusions. The ESI† provides computational details, force field files, and molecular models ready to use.

Results and discussion

Chemical features and their translation into force field parameters

Summary of energy expressions. We utilize the energy expressions of IFF,³⁷ CVFF,⁶⁸ CHARMM,⁶⁹ Dreiding,⁷⁰ AMBER,⁷¹ and OPLS-AA,⁷² which are widely used and suitable for compounds across the periodic table:

The energy expressions employ a 12-6 Lennard-Jones potential, including small differences in combination rules and more notable differences in the scaling of nonbonded interactions between 1,4 bonded atoms.^29,30 In addition, we report force field parameters for the IFF, CFF,⁷³ PCFF,⁷⁴ and COMPASS⁷⁵ energy expressions that employ a 9-6 Lennard-Jones potential, Waldmann-Hagler combination rules, and include nonbonded interactions between all 1,4 bonded atoms:

The force field parameters in eqn (1) and (2) include equilibrium bond lengths r_0,ij, harmonic bond stretching constants K_r,ij, equilibrium bond angles θ_0,ijk, angle bending constants K_θ,ijk, atomic charges q_i, as well as the equilibrium nonbond diameters σ_0,ii and the nonbond well depth ε₀,_ii. The latter two parameters are specific to the respective LJ potentials. The bonded parameters r_0,ij, θ_0,ijk, K_r,ij, and K_θ,ijk are used for atoms that are part of predominantly covalent bonds (bonds with less than half ionic character). The nonbonded parameters q_i, σ_0,ii and ε₀,_ii apply to all pairs of atoms, except bonded atoms with 1,2 and 1,3 covalent connections. Further additive terms for torsions, out-of-plane, higher order (cubic, quartic) and cross-terms are not required for MoS₂ and have no impact on compatibility with existing parameters for biomolecules and inorganic compounds (see Section S2 in the ESI† for more details).³⁷

Novel representation of structure and chemical bonding. A critical and central feature is the polarity of the chemical bonds. Atomic charges represent internal dipoles and multipoles in the all-atom model and are considered to be the same for all energy expressions.³⁶ We assign the atomic charges using the extended Born model, which captures the nature of chemical bonding and reproduces interfacial properties with solvents and other inorganic compounds and organic molecules.^36,38 Prior studies include no rationale for the atomic charges chosen and rely on ad-hoc assumptions from DFT calculations, which vary from one density functional or from partition method to another by several 100%. These approaches cannot consistently describe internal polarity using point-based charges.^36,38,76,77

We further utilized equilibrium bond lengths and bond angles from reproducible X-ray data. Minor adjustments within few percent were permitted to account for small superimposed contributions by the nonbond interactions (last two terms in eqn (1) and (2)). Prior methods have similarly relied on X-ray data for the bond length. However, bond angles were oversimplified or included as empirical fit parameters, disregarding existing experimental knowledge. The parameters for bonds and angles, for physical reasons, are furthermore about the same for all energy expressions regardless of 12-6 or 9-6 Lennard-Jones potentials. A small adjustment is recommended for OPLS-AA, related to different scaling of nonbond interactions between 1, 4 bonded atoms in OPLS-AA (0.5 vdW, 0.5 Coulomb) relative to AMBER (0.5 vdW, 5/6th Coulomb) and other force fields (1.0 for vdW and Coulomb).³⁷ Differences in combination rules for 12-6 LJ parameters have negligible influence for MoS₂ (arithmetic mean of ε₀,_ii in CHARMM and AMBER versus geometric mean in CVFF and OPLS-AA). Our force field considers details of compatibility and portability of chemical features and parameters for 2D materials for the first time.

We assigned vibration constants for bonds and angles to qualitatively reproduce Infrared and Raman frequencies and obtain reliable mechanical properties. Most earlier models do not consider realistic vibration frequencies and encounter stability problems, such as unphysical buckling, detrimental for the function of the model (Fig. 3). Another critical contribution are the Lennard-Jones parameters that represent relative atomic radii, short-range repulsion between atoms, and van-der-Waals attraction. Earlier studies include no interpretation and validation for LJ parameters while we discuss and assign the values consistent with nearby elements in the periodic table and with cleavage energies. The nonbond diameters σ_0,ii correlate with van-der-Waals diameters of Mo and S. The well depths ε₀,_ii play a critical role to reproduce atomic polarizability, cleavage energies and, as a result, interactions with solvents, organic molecules, and other compounds, as verified by contact angles. The downside of neglecting such balances can be seen in the work of (SR2017).⁴⁹ Random (and inconsistent) fitting to several properties including water contact angles was performed, arriving at a well depth of sulfur that is about 100% too large relative to similar sulfides. For another solvent, diiodomethane, a 100% mismatch of the contact angle relative to experiment was then reported, as well as other coarse deviations in surface energy by +60% and in energy derivatives.

Fig. 3 Comparison of experimental data (red) and results from molecular dynamics simulation (black) for structural, mechanical, and vibrational properties of 2H-MoS₂. (a) XRD patterns show a good match (experimental data from ref. 80). (b) The compressibility of bulk MoS₂ is nearly identical in experiment and simulation (experimental data from ref. 63). (c) Experimental infrared data (from ref. 58 and 82) are reasonably approximated by MD simulation, although the intensities are difficult to reproduce in MD. MD simulation also captures additional Raman peaks (see experimental details in ref. 65). A simulated IR spectrum from DFT calculations with the rPBE functional (green) approximately reproduces the wavenumbers and relative intensities of the experimental data.

Rationale and implementation. Therefore, our rationale drastically differs from earlier blind and incomplete fitting methods. We present (1) the first stable force field of MoS₂ and (2) at least 10 times higher accuracy in key interfacial properties (see also ref. 28, 30, 37 and 53). 2H-MoS₂ consists of S–Mo–S sandwich layers with the space group P6/mmc (Fig. 1).⁵⁵ It is the most stable polymorph compared to two other polymorphs 1T-MoS₂ and 3R-MoS₂ under standard conditions (298 K and 101.3 kPa). The two S–Mo–S sandwich layers in the unit cell are related by a 180° rotation along the c-axis and, accordingly, atom types in the 1^st layer can be used to describe the 2^nd layer. The atomic positions from X-ray data suggest strong covalent Mo–S bonding within the layers and weaker interlayer forces. Atomic charges were identified as +0.5e for Mo and −0.25e for S, respectively, with about 10% uncertainty. These values are consistent with a high atomization energy of Mo, similar to C, which indicates strong covalent bonding to S according to the extended Born model.³⁶ The ionization energy of Mo is lower than that of carbon and the coordination number higher than in carbon sulfides and thiols, indicating more ionic character (up to +0.8e).³⁶ However, the Mo–S bond is also longer at 239 pm compared to C–S bonds in CS₂ and in thiocarbonyl compounds at ∼150 pm, and therefore comparable internal dipole moments result from a lower Mo charge of 0.50 ± 0.05e. The specific value of +0.50e was confirmed by a good match to computed contact angles with water and a less polar liquid (CH₂I₂), which are both sensitive to internal polarity.

2H-MoS₂ has only one Mo–S bond type and more than five types of bond angles. For simplicity, we can work with two atom types of Mo and four atom types of S. The model can then be designed with five types of bond angles θ₁, θ₂, θ₃, θ₄, θ₅,^39,40 whereby angle θ₂ can assume two possible values (78.38° and 134.4°) and an angle bending constant K_θ,ijk of zero (Table 2). When one θ₂ angle (S1–Mo–S3 angle) would change, the other θ₂ angle on the same Mo atom would change by the same amount (Fig. 1). Our chosen assignment avoids the definition of further atom types, which is possible but would increase the complexity of the model. The five angles help reproduce the lattice parameters of MoS₂, and the force constants were chosen to reproduce approximate IR/Raman vibration frequencies from experiments. Bond rotation and out-of-plane force constants are zero due to the absence of such deformations among bonded atoms in MoS₂. In some prior models,^43,49 only two types of angles, Mo–S–Mo and S–Mo–S, were included, then resulting in failure to capture the geometry and bonding environment (Fig. 2a and b), as well as in bending artifacts and unphysical crumpling of larger layers (Fig. 2c–e). The new models replicate bending stability consistent with the flat geometry and uniform height of MoS₂ layers seen in AFM images, analogous to earlier IFF models for clay minerals (Fig. 2c–e).^53,78

Table 2 Force field parameters for 2H-MoS₂ (IFF)

I. Nonbond	Atomic charge (e)	σ (pm)		ε (kcal mol^−1)
		CVFF, CHARMM, AMBER, DREIDING, OPLS-AA (12-6 LJ)	PCFF, COMPASS (9-6 LJ)	CVFF, CHARMM, AMBER, DREIDING, OPLS-AA (12-6 LJ)	PCFF, COMPASS (9-6 LJ)
a The equilibrium bond length (r0,ij) and bond stretching constant (Kr) are the same for all force fields. b The equilibrium bond angles (θ0,ijk) and angle bending constants (Kθ) are the same for all the force fields. The values for θ0,ijk are equal to data from X-ray diffraction (Fig. 1). We recommend a minor modification of angles θ1 and θ4 for OPLS-AA, which uses significantly different scaling rules for nonbond interactions between 1, 4 bonded atoms. A complete list of individual angles and angle bending constants is given in Table S2 in the ESI. c The value in brackets is recommended for OPLS-AA to best reproduce bond length and lattice parameters. The default angle of 84.32° still facilitates good performance. d Two values are possible for θ2: θ2a and θ2b. The values are represented using an angle bending constant of zero without loss of accuracy (see Table S2 in the ESI for details).
Mo	+0.5	480	501	0.07	0.054
S	−0.25	384	398	0.30	0.26

---

II. Bond^a	r 0,ij (pm)	K r (kcal mol^−1 Å^−2)
Mo–S	239	118

---

III. Angles^b	θ 0,ijk (°)	K θ (kcal mol^−1 rad^−2)
θ 1 and θ4	84.32 (OPLS only: 83.5)^c	205
θ 2	78.38 (θ2a), 134.4 (θ2b)^d	0
θ 3	134.4	3.6
θ 5	78.38	3.6

---

The assignment of Lennard-Jones equilibrium nonbond diameters σ_0,ii was guided by known crystallographic radii,⁷⁹ according to which Mo is somewhat larger than S, followed by numerical refinement (Table 1). σ_0,ii assumes a 4% larger value in the 9-6 LJ potential versus the 12-6 LJ potential to compensate decreased repulsion.^28,30,53 The exact magnitude of the equilibrium nonbond diameters σ_0,ii was determined so that lattice parameters closely agree with experimental data, and refined in combination with ε₀,_ii to reproduce the cleavage energy and contact angles. The nonbond well depths ε₀,_ii represent atomic polarizabilities and assume a minor repulsive role to counterbalance internal Coulomb attraction. The ε₀,_ii value for Mo is lower than for S since the positively charged metal atom is less polarizable than the negatively charged S atom. The well depth of S atoms in MoS₂ (−0.30 or −0.26 kcal mol⁻¹) is about the same as that of S atoms in thiols and in similar chemical environments.³⁷ In contrast, well depths in SR (2017) were fitted to about twice this value without explanation.⁴⁹ We obtain mechanical properties in the correct range without further parameter adjustments.

Validation of structural, vibrational, and mechanical properties

Lattice parameters were tested on a supercell constructed from multiples of the unit cell (Table 3). The deviation from experiment is consistently <0.5% and the density agrees better than 1% with available measurements. The use of different energy expressions with minor changes in combination rules causes only minor changes (Table 3). The accuracy is better than with DFT and indicates that the force field mimics the real lattice configuration with minor discrepancies (Table 1).

Table 3 Lattice parameters for a (18 × 18 × 2) super cell of 2H-MoS₂ in experiment (ref. 55) and calculated from NPT molecular dynamics with the new parameters (IFF) using CVFF, Dreiding, PCFF, COMPASS, CHARMM, AMBER and OPLS-AA energy expressions at 298.15 K and 1.013 MPa pressure

Method	a (Å)	b (Å)	c (Å)	α = β (°)	γ (°)	r 0 (Å)	ρ (g cm^−3)
Expt (XRD)	56.89	56.89	24.59	90	120	2.38	5.00
CVFF, Dreiding	56.89	56.89	24.50	90	120	2.36	5.03
PCFF, COMPASS	56.92	56.92	24.58	90	120	2.36	5.00
CHARMM	56.92	56.92	24.50	90	120	2.36	4.99
AMBER	57.11	57.20	24.51	90	120	2.37	4.97
OPLS-AA	57.23	57.19	24.57	90	120	2.37	4.96

---

Accordingly, the computed XRD pattern is in good agreement with experimental data (JCPDS # 37-1492) (Fig. 3a).⁸⁰ Likewise, computed mechanical properties show an impressive match to experiments (Fig. 3b). The compressibility was calculated under different pressures (0.001, 1, 1.5, 2, 3, 4 kbar) in NPT molecular dynamics simulations by recording the volume change in comparison to experimental data.⁶³ The two curves are essentially identical up to about 2 kbar. Even when the pressure exceeds 2 kbar, the difference between computation and experimental data remains under 2%. Further mechanical properties include the Young's modulus in-plane, computed as 175 ± 2 GPa and measured in a range from 170 to 400 GPa.^59–62 The agreement is very good and the computed value is similar to the in-plane modulus of 160 ± 10 GPa of clay minerals.⁸¹ Elastic properties help to validate the consistency of the force field parameters, along with the cleavage energy and hydration energy, as they represent the second derivatives of the energy with respect to coordinates. The original energy landscape and the derivatives agree with experimental data within few %. In contrast, earlier models feature over 50% deviation of the original energy (cleavage energy), having no compatibility with other inorganic and organic compounds.⁴⁹

The experimental vibration spectrum shows two characteristic adsorption peaks at 384 and 470 cm⁻¹ that belong to the infrared-active E¹_1u and A¹_2u modes of MoS₂ (Fig. 3c).^58,82 The simple classical model with only three bond stretching and angle bending constants (Table 2) reproduces the frequencies with a 20–30 cm⁻¹ blue shift (360 and 450 cm⁻¹). In experiment, also a major A_1g Raman band is observed at 409 cm⁻¹. This band is not visible in the IR spectrum, however, MD simulations reveal all vibrations regardless of selection rules and the Raman band can be seen near 410 cm⁻¹ (Fig. 3c). DFT calculations of the IR spectrum with the rPBE functional are close to experimental wavenumbers and qualitatively reproduce the intensities.⁶⁵ An excellent review of vibrational and optical properties is given in ref. 65. The performance of the classical model relative to earlier force fields is excellent while, as a limitation, intensities cannot be reproduced without quantum mechanical details (Table 1).

Surface and interfacial properties

Molecular dynamics simulations of the cleavage energy, MoS₂–water, and MoS₂–diiodomethane interactions were carried out using the NAMD program in comparison to experimental data and to cleavage energies computed by DFT calculations (Fig. 4 and Section S1 in the ESI†). The cleavage energy of the basal plane of a layered material is a key material parameter for its application in devices and for the validation of the model Hamiltonian. Unfortunately, the quantitative analysis of the cleavage energy can be challenging in experiments and a wide range of values from 40 to 121 mJ m⁻² has been reported for MoS₂.^54,57,83 The differences are related to the surface quality (well-ordered and clean) and the analysis method. Tang et al. measured the surface tension of MoS₂ basal plane as 110 ± 10 mJ m⁻² using an in situ transmission electron microscopy probing technique.⁸⁴ In 2017, Otyepková et al. used inverse gas chromatography (IGC) to measure and analyze the surface properties of MoS₂.⁵⁴ For a natural MoS₂ sample, the surface energy was determined in the range from 99 mJ m⁻² at a higher surface coverage of 20% gas probe to 121 mJ m⁻² at a low surface coverage of 1% gas probe.⁵⁴ The higher values are likely closer to the pure system. Separately, we considered data from quantum mechanics at the DFT level for further reference. These data tend to be less reliable and show considerable scatter between 160 and 284 mJ m⁻² (Table 1).^19–21 On balance, we consider the data at the lower end for further reference,¹⁹ arriving at a best estimate of 150 ± 10 mJ m⁻² for the cleavage energy as a reference value. The cleavage energy was computed with the force field using molecular dynamics simulations in the NVT ensemble at 298 K as an energy difference using two boxes (Fig. 4). One box contained a surface slab of greater than 2 nm thickness separated with a 40 Å vacuum layer (Fig. 4a), and the other box a combined 3D periodic model of the same number of atoms without vacuum (Fig. 4b). The results agree with the target data of 150 mJ m⁻² for IFF-CHARMM. The use of different energy expressions consistently yields 148 ± 2 mJ m⁻² (Table S4 in the ESI†).

Fig. 4 Computation of surface and interfacial properties of 2H-MoS₂. (a and b) The models used to calculate the cleavage energy comprise a cleaved surface slab of a thickness of 4 MoS₂ layers and a 40 Å vacuum layer, as well as the equivalent 3D periodic bulk system for reference. A is the surface area. The computed cleavage energy agrees with the reference data from experiments and from DFT. (c and d) Representative equilibrium snapshots of water and CH₂I₂ on the 2H-MoS₂ surface at 298 K. Computed contact angles θ agree with laboratory data without further adjustable parameters. A new method to determine the contact angle was employed that is less ambiguous than prior methods. It involves drawing a circle of a size that best approximates the contour of the liquid and reading out the angle between the diameter that ends at the solid–liquid interface and the diameter aligned with the surface normal (see black construction in c).

The wetting behavior of MoS₂ was studied for two solvents, polar and hydrogen-bonded water as well as somewhat less polar diiodomethane (CH₂I₂), similar to previous studies on silica surfaces, graphite, and clays.^28,34,51 An artificial cubic-shaped cluster of solvent molecules was located on an extended MoS₂ surface and subjected to MD simulations for equilibration into a continuous semi-cylindrical liquid continuous in one dimension (Fig. S1 in the ESI†). We used the convergence of the total energy versus time to determine the equilibrium state (Fig. S3 in the ESI†). After several nanoseconds, the change in amplitude was less than 0.01% of the total energy. The contact angle was measured using a new circle-based method as an average over 100 snapshots in equilibrium for each simulation (Fig. 4c and S2 in the ESI†). The circle-based method reduces biases and uncertainties on the order of ±5° associated with individual readings in earlier studies^28,49,50 to a small uncertainty of ±2° or only ±1°. We obtained a contact angle of water of 69 ± 2° (Fig. 4c) and a contact angle of 15 ± 2° for CH₂I₂ (Fig. 4d) using the MoS₂ parameters in IFF-CHARMM format. The two values closely match data from laboratory measurements on freshly cleaved MoS₂ surfaces, which are 69 ± 3.8° and 15 ± 2°, respectively.⁵² The uncertainty in experimental measurements arises from differences in advancing contact angles versus receding contact angles and is included to characterize the overall reproducibility. Quantitative agreement of computed cleavage energies, contact angles, and mechanical properties with experimental data without further adjustable parameters documents an unprecedented level of internal consistency of the force field for MoS₂ and its suitability for property predictions of other mixed-phase systems.

For the computation of the contact angles, we used the TIP3P water model and further tests with the flexible SPC water model yield nearly identical results within ±1°. A similarly small dependence on the water model was found earlier on silica,²⁸ metal,⁸⁵ and graphite surfaces³⁴ using chemically consistent, interpretable parameters in IFF.³⁷ Deviations may be larger if other, chemically inconsistent and non-interpretable force field parameters are used. The CVFF parameters of CH₂I₂ were updated to reproduce dipole moment, density at 298 K, and vaporization energy at the boiling point 480 K in agreement with experimental data (Table S3 in the ESI†).

Earlier computational studies of the contact angle also involved tuned interaction parameters to match contact angles of fresh and aged MoS₂ surfaces,⁸⁶ however, in this approach the chemical changes to the surface such as contamination by airborne hydrocarbons or oxidation are not considered.⁶⁵ In our approach, the interaction parameters have a clear rationale, do not require customized fits to reproduce contact angles, and changes in surface chemistry can be explicitly included.

Peptide recognition on MoS₂ surfaces and identification of binding residues

Functionalization of MoS₂ with peptides and other functional ligands finds promising applications in sensing, catalysis, hierarchical assembly, and 3D-printed biomaterials.^4–10 Several studies on the adsorption of peptide molecules on the MoS₂ surface have been reported and indicate that the hydrophobicity, aromatic structure and electrostatic properties of amino acids are important factors affecting adsorption and surface assembly.^87–92 However, experimental data remain qualitative in nature, including AFM, TEM imaging, binding information from concentration measurements, QCM, XPS, CD, Raman shifts, and relative ranking from washing cycles in phage display. Quantitative binding energies and conformations in solutions are not accessible. As an example application of this model, we determine the binding strength and conformations of eight different peptides on the MoS₂ surface in aqueous solution at pH 7 (Fig. 5 and 6). We discuss the results in the context of experimental findings and develop criteria for selective peptide design. This analysis is not part of our model derivation.

Fig. 5 Adsorption characteristics of three peptides identified by phage display and a control sequence on 2H-MoS₂ surfaces in aqueous solution at 298.15 K and 101.3 kPa from molecular dynamics simulation. Side and top views of adsorbed conformations in equilibrium are shown, including highlights of specific binding features. Adsorption energies, the average number of displaced water molecules on the surface, and the percentage of contact time for each residue (within 3.5 Å from the MoS₂ surface) are identified. The color code in the bar charts (c, f, i, l) represents the hydrophobicity according to the Eisenberg and Weiss scale (ref. 98).

Fig. 6 Adsorption characteristics of three peptide mutants (from ref. 89) and a control sequence (from ref. 93) on 2H-MoS₂ surfaces in aqueous solution at 298.15 K and 101.3 kPa from molecular dynamics simulations. Side and top views of adsorbed conformations in equilibrium are shown, including highlights of specific binding features. Adsorption energies, the average number of displaced water molecules on the surface, and the percentage of contact time for each residue (within 3.5 Å from the MoS₂ surface) are indicated. The color code in the bar charts (c, f, i, l) represents the hydrophobicity according to the Eisenberg and Weiss scale (ref. 98).

Choice of peptides and characterization of binding. The peptides include YSATFTY,⁹² TSHMSNT, YIPHTPN (this study), GGGGGGG (a control sequence), three mutants YGAGAGAY, RGAGAGAR, EGAGAGAE,⁸⁹ and another control peptide TLTTLTN (a Pt binding peptide)^93,94 with sequences given in N → C order. We monitored the conformations, adsorption energies, average number of displaced water molecules, as well as the contact time of individual amino acid residues on the surface (Fig. 5 and 6). We utilized the IFF-CHARMM parameters for molecular dynamics simulations with multiple replicas and 30 ns and 100 ns simulation time. The models included single peptides dissolved in 2000 water molecules (TIP3P model) at a low surface coverage of 8% to 20%. The onset of equilibrium was determined using convergence of the total energy versus time as previously reported (Fig. S3 in the ESI†).^24,26,95 The adsorption energy corresponds to the difference between the average energy of the peptide in the adsorbed state versus the average energy of the peptide desorbed from the MoS₂ surface at >4 nm distance.²⁴ The average number of water molecules displaced by the peptides upon adsorption equals the difference in the number of water molecules within 3.5 Å distance from the MoS₂ surface atomic layer. 3.5 Å is the typical thickness of the first layer of water molecules in contact with the MoS₂ surface, and the structure of this water layer is most affected by the surface and by adsorption of the peptide. The contact time for each residue with the surface was measured in % of time, equal to the number of snapshots in which any atom of the residue was located within 3.5 Å from the MoS₂ surface in relation to the total number of snapshots in equilibrium (at least 1000). Full details of the methods, including a Python script, are described in Section S1 in the ESI.†⁹⁶

Binding affinity and molecular mechanism. The peptide sequence YSATFTY, containing acylated N-terminal and amidated C-terminal groups, was identified in phage display experiments as the strongest known binder to the MoS₂ surface to-date (Fig. 5a–c). The peptide shows the largest negative adsorption energy (ΔE_ads) of −86 ± 6 kcal mol⁻¹ and enables continuous assembly of ordered peptide films on MoS₂ surfaces (see also Section S3 in the ESI†).⁹² A single molecule of YSATFTY replaces about 25 water molecules on the MoS₂ surface and binds via the backbone with all groups, especially the side groups Tyr and Phe. The phenyl rings allow effective epitaxial packing onto superficial sulfur atoms with optimum van-der-Waals contacts (highlight in Fig. 5b). Simultaneously, the OH groups in Tyr and in other residues maintain favorable interactions with water (Fig. 5a). The peptide YSATFTY also adsorbs strongly due to limited solubility in water, which can be seen in a transition from a globule-like equilibrium conformation in solution to a largely flat-on conformation on the surface within less than 10 ns simulation time (Movie S1 in the ESI†). Next, the peptide TSHMSNT was identified as a strong binder in phage display (Fig. 5d–f). However, the adsorption energy was only −19 ± 2 kcal mol⁻¹ and the peptide only replaced 4 water molecules on the MoS₂ surface. The gap between TSHMSNT and YSATFTY is large since residues for strongest binding (Tyr, Phe) are missing. Even though contact times are similar to YSATFTY, the contact area per residue is low. Asn and Met contribute significantly to binding while Ser and Thr provide mainly backbone flexibility. His, in the non-protonated state, prefers to stay in solution as the nitrogen atoms in the imidazole side group carry negative charges of −0.7e and −0.36e, respectively (highlights in Fig. 5d and S4 in the ESI†).

Peptide GGGGGGG was tested as a control sequence to examine the contribution of the backbone to adsorption (Fig. 5g–i). The computed binding energy of −9 ± 2 kcal mol⁻¹ concurs with no specific binding expectations, although it is quite significant given the absence of side groups. GGGGGGG displaced 2.7 water molecules and aliphatic hydrogen atoms of the Gly backbone were found closest to the surface. The linear conformation, related to the absence of intramolecular hydrogen bonds, shows that flexible backbone contact alone does not lead to strong binding. A similarly small binding energy of −8 ± 3 kcal mol⁻¹ was observed for the peptide YIPHTPN (Fig. 5j–l). YIPHTPN was identified as a weak binder in phage display and contains the potentially strongly binding Tyr. However, two Pro residues introduce local helicity in the backbone and drastically reduce effective surface contact (highlights in Fig. 5j).⁹⁷ Steric hindrance also originates from bulky Ile. Hydrophobic side groups, such as –CH₃ in Ala (Fig. 5a–c) and butyl groups in Ile (and Leu) seek surface contact with MoS₂ to avoid contact with water. Similarly, aliphatic hydrogen atoms in Pro and Thr had van-der-Waals interactions with the MoS₂ surface with limited contact area. Thereby, large Pro, Ile, and Leu groups diminish potential binding of neighbor groups through steric demand. In YIPHTPN, His maintained contact with the surface due to conformation constraints from neighboring Pro residues.

The order of computed adsorption energies for these four 7-mer peptides agrees with the relative binding affinity according to experimental data from phage display. Adsorption energies vary widely from −86 to −8 kcal mol⁻¹ and the data clearly show that side groups have a greater influence on the adsorption than the backbone.

Furthermore, we examined three 8-mer peptide mutants with differently charged end groups YGAGAGAY (±0), RGAGAGAR (+2)·2Cl⁻ and EGAGAGAE (−2)·2Na⁺, which have been characterized in experiments by Hayamizu's team (Fig. 6).⁸⁹ The computed adsorption energies are −26 ± 2 kcal mol⁻¹, −17 ± 3 kcal mol⁻¹, and −6 ± 2 kcal mol⁻¹, respectively, and match the order of the experimental binding affinity. The neutral peptide YGAGAGAY was most attracted to the MoS₂ surface, whereby the aromatic ring in Tyr undergoes coordination with polarizable sulfur atoms on the surface, and the phenolic OH group maintains favorable contact with water (Fig. 6a–c). The GAGAGA backbone, however, does not leverage binding. In comparison, the peptide YSATFTY with the strongest adsorption (−86 kcal mol⁻¹) also contains two Tyr residues located at the ends with similar coordination on the surface, plus additional binding residues in a suitable sequence (Fig. 5a–c). We also noted that Tyr in a terminal position, such as in YIPHTPN, can result in marginal adsorption strength (−8 ± 2 kcal mol⁻¹ in Fig. 5j–l). Then, the 2 Pro residues induce 2× local helical twists that disrupt contact with the surface and make binding of Tyr and Asn ineffective (Fig. 5j–l). Therefore, the effect of individual residues is not the only reason for the high binding affinity, but also the sequence along the backbone. Among the charged mutants, the positively charged peptide RGAGAGAR (+2) clearly prefers the MoS₂ surface (−17 kcal mol⁻¹) over the negatively charged peptide EGAGAGAE (−2) (−6 kcal mol⁻¹). Hereby, multiple positively charged H atoms at the N atoms in the guanidinium group (+0.4e) coordinate in a well-matched geometry at least three negatively charged sulfur atoms (−0.25e) exposed on the MoS₂ surface (Fig. 6d–f). Therefore, Arg significantly binds to the MoS₂ surface. In contrast, the peptide EGAGAGAE with negatively charged carboxylate groups does not coordinate well with the MoS₂ surface, which consists of negatively charged sulfur atoms while the positively charged Mo atoms (+0.5e) are tucked away underneath the surface (Fig. 6g–i). Weak backbone contact (1.9 water molecules replaced) plus some repulsion of Glu lower the binding energy to −6 kcal mol⁻¹, consistent with lowest binding strength and instability of EGAGAGAE peptide films on MoS₂ surfaces in experiments.⁸⁹

A small binding energy of −6 ± 2 kcal mol⁻¹ was also observed for a Pt-binding peptide TLTTLTN (T7)⁹³ (Fig. 6j–l). This control sequence, identified by phage display against a different material, was also expected of random affinity to MoS₂. The simulation predicts, in agreement, no specific attraction towards the MoS₂ surface and only 1.6 displaced water molecules. Leu diminishes binding, like Ile in YIPHTPN (Fig. 5j–l), for steric reasons and lack of strong affinity. The high Thr content in TLTTLTN allows conformational flexibility and favorable interaction with supernatant water, which, in the absence of strongly binding residues, lowers attraction to the MoS₂ surface (highlights in Fig. 6j). We therefore recognize that several neutral polar groups in TLTTLTN, without other residues that induce a specific geometry match to the MoS₂ surface (e.g., Tyr) or a match in polarity (Arg), tend to be more attracted to water and result in weak adsorption. In comparison, GGGGGGG without any specific side groups shows slightly better adsorption (−9 kcal mol⁻¹) (Fig. 5g–i).

Surface interactions can also be further characterized by free energy calculations. However, the computational expense is high and we consider the average contact time of individual residues as a comparable measure.

Prediction of binding selectivity, ranking of amino acids, and discussion in the context of prior studies. The results show in unprecedented detail how simulations can be used to obtain quantitative forecasts of peptide affinity to MoS₂ for any amino acid sequence. Overall, the binding energy shows a linear correlation with the number of displaced water molecules on the MoS₂ surface (Fig. 7). Peptides that facilitate more exchange have a larger negative binding energy. The release of more water molecules into the bulk solution allows the formation of more hydrogen bonds with each other compared to water in contact with the partly hydrophobic MoS₂ surface. Simultaneous adsorption of the peptide backbone to the surface and formation of hydrogen bonds between protruding side groups and supernatant water molecules is also favorable (Fig. 5a). The peptide YSATFTY exhibits only few gaps in contact with the MoS₂ surface, followed by significantly less surface contact of peptides YGAGAGAY, TSHMSNT, and RGAGAGAR (Fig. 5 and 6). Weakly binding peptides GGGGGGG, EGAGAGAE, YIPHTPN, and TLTTLTN replace only a minor amount of water molecules.

Fig. 7 Relationship between the adsorption energy of peptides and the number of replaced water molecules in the first molecular layer on the MoS₂ surface. On average, every replaced water molecule contributes −3.3 kcal mol⁻¹ binding energy.

The wide range of possible binding energies thus arises from concerted interactions of multiple amino acids with the surface. Strongly binding residues include Tyr due to epitaxial fit and a polar OH group that enhances contact with water, as well as Phe (and likely Trp). Binding is also enhanced by Arg due to a good epitaxial fit of the guanidinium group to the MoS₂ surface and matching of positive atomic charges to opposite MoS₂ surface charges. Some electrostatic binding may occur for Lys, although without the geometry advantage of Arg. Asn, Gln, and Met tend to somewhat support binding. His, in the nonprotonated state, has reduced surface attraction due to negative charges on the nitrogen atoms. Ser and Thr allow backbone flexibility and interactions with water. Gly and Ala (and likely Cys) also support flexible conformations and are weak binders themselves. Val has unfavorable effects on adsorption due to steric hindrance. Leu and Ile are sterically demanding and tend to reduce adsorption even further. Negatively charged residues such as Asp and Glu diminish adsorption due to unfavorable electrostatic interactions. Pro induces helicity in a peptide and strongly disfavors binding by lowering the number of displaced water molecules.

These trends are largely supported by further experimental observations. A study of two amyloid peptides in contact with MoS₂ (SNNFGAILSS and GLMVGGVVIA) reported stronger binding of Ser, Asn, and Phe to MoS₂ than binding of Met, Val, and Leu.⁸⁷ Investigations of the binding of a dodecapeptide (HLLQPTQNPFRN) to MoS₂ suggested significant binding of Phe and His to MoS₂, as well as some binding of Arg, and no binding of Leu.⁹⁰ These data, however, are at least partly based on simulations and have some uncertainty. Other reported peptides with high binding constants in experiments include GVIHRNDQWTAPGGG and DRWVARDPASIFGGG, as well as a weakly binding/non-binding peptide SVMNTSTKDAIEGGG.⁸⁸ The interpretation of binding residues in ref. 88 relied exclusively on MD simulation with a misleading force field, indicating Val, Pro, and Trp as binding residues, whereas Val and Pro are expected to bind weakly according to other studies. A more likely interpretation of the experimental data in ref. 88 is to consider Trp, Phe, Arg, Asn, Gln, and perhaps His as significantly binding residues in the first two peptides, with Ser and Thr serving as linkers for conformational flexibility. These binding residues, except 1 Asn, are missing in the 3^rd nonbinding peptide, and only 1 intermediately binding Met is present, explaining the difference in relative binding strength and far lower adsorbed amount in experiment. Tyr was clearly the strongest binder in the series of X-(GA)_n-X peptides, followed by Arg, as analyzed in detail here.⁸⁹ Studies of the self-assembly of two other peptides on MoS₂, LEY (LEYLLEYLL) and LKY (LKYLLKYLL) show strong binding of Tyr and a trend towards more adsorption of LKY than LEY at low peptide concentration (comparable to single peptide adsorption considered here).⁹¹ Preferred adsorption of the positively charged, Lys-containing peptide LKY over the negatively charged Glu-containing peptide LEY agrees with our interpretation.⁸⁹ At the same time, very different assembly patterns were observed for multiple LKY and LEY peptides.⁹¹ LKY showed higher affinity to MoS₂ and larger film thickness at low peptide concentration in contrast to higher stability of LEY at high peptide concentration. The major driver for binding was clearly Tyr in both cases.^89,91

In summary, simulations and available experimental data enable a consistent interpretation of peptide binding to MoS₂ and design rules for MoS₂-binding single peptides. Computational screening allows quantitative forecasts of peptide affinity for new amino acid sequences in atomic detail, many times better than with existing force fields. Our analysis further shows that traditional descriptors such as amino acid polarity, hydrophilicity, aromaticity, and charge are not sufficient to identify the multifactorial rules of binding. The reason is that such identifiers only characterize the peptides and neglect the unique surface chemistry and geometry of the substrate. For example, there is no direct correspondence between the contact time of the peptides and the number of hydrophobic groups or the sum of hydrophobicity indices. Design rules are uncovered when considering the chemistry of both the substrate and the peptides, as previously shown for metals and oxides.^{95,97,99–102} The binding strength also further depends on surface coverage, pH, and ionic strength.¹⁰³ The assembly of multiple peptides and films across large areas remains challenging to predict, however, computational screening and experimental studies can take these conditions into account and be used in tandem to accelerate progress.

Potential broader uses and applications

The models can be applied to study protein-, DNA-, lipid-, and polymer-MoS₂ nanomaterials in 4D resolution, MoS₂ electrode–electrolyte interfaces, as well as mixed MoS₂ materials with graphene, metals, and minerals in composites and capacitors. Molecular dynamics simulations allow the quantitative, dynamic analysis at the level of individual atoms and functional groups up to the 1000 nm scale and support the discovery of specific materials formulations and properties. Substitutions by dopant elements can be implemented by changes in atomic charges and LJ parameters (similar to clay minerals, see ref. 53). The reliability of the simulations relative to experimental data is multiple times higher than prior force fields and clearly better than DFT, opening up interesting perspectives for computation-driven discoveries. The consistent representation of chemical bonding, structure, and energies allows coupling with quantum-mechanical methods at the local scale to analyze the complete electronic structure, as well as integration with coarse-grain and continuum models towards the micrometer scale. Stress–strain curves including failure can be computed by replacing the harmonic bonded potential for Mo–S bonds with a Morse potential (IFF-R).

Conclusions

We present interpretable force field parameters for atomistic simulations of MoS₂ and related nanomaterials and biomaterials, compatible with major force fields and programs for molecular simulations (IFF, CVFF, PCFF, DREIDING, CHARMM, AMBER, and OPLS-AA). The models reproduce lattice parameters, infrared frequencies, surface and interfacial properties, as well as mechanical properties relative to experiments in multiple times higher accuracy compared to prior models. Prior force fields for MoS₂ fail to reproduce the physical structure and energies relative to experimental data due to lack of a rationale for parameters and severe internal inconsistencies. The model overcomes these limitations and bypasses questionable assumptions in DFT by reliance on reproducible experimental data and established theory, performing clearly better than current density functionals in interfacial and mechanical properties.

We utilized the models to analyze the binding mechanism of 8 different peptides on the MoS₂ surface in aqueous solution. Specific amino acid sequences lead to a wide range of adsorption energies from −6 to −86 kcal mol⁻¹. More surface contact and release of surface-bound water increase the negative binding energy. We ranked all 20 amino acids in relative binding strength, consistent with the interpretation of extensive experimental data and insight from simulation. Strong binding is mediated by aromatic amino acids via favorable epitaxial coordination to superficial sulfur atoms (e.g., Tyr, Phe), by positively charged Arg through complementary charge distribution and geometry of guanidinium groups to the superficial sulfur lattice. Conformationally flexible, weakly polar residues with side chains can support binding and simultaneously contribute favorable hydrogen bonding to water molecules in the solution phase. Negatively charged residues diminish binding due to mild electrostatic repulsion, and bulky hydrophobic residues (Pro, Ile, Leu) strongly diminish binding and surface contact. The binding selectivity arises from unique properties of the peptides and of the surface, and is not possible to derive using traditional criteria such as hydrophobicity that only focus on the organic molecules and neglect the substrate.

The models explore the performance limit of current energy expressions and can be used for property predictions of MoS₂-containing nanomaterials and biomaterials, as well as in-depth, full atomic-scale characterization of assemblies up to a billion atoms. Full electrolyte conditions, realistic dynamics, and combinations with electronic structure calculations at the local scale may be used to assess band gaps, conductivity, and catalytic activity. Potential uses include electrocatalysts, supercapacitors, polymer composites, biosensors, and mixed layered materials containing MoS₂. Coupling with atomic electron tomography measurements can further correlate crystal defects with the physical and chemical properties of 2D materials.^104,105 Extensions to other 2D chalcogenides are in progress.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Science Foundation (DMREF-1623947, CBET-1530790, OAC 1931587), and the University of Colorado-Boulder. J. M. acknowledges the support by STROBE: A National Science Foundation Science & Technology Center (DMR-1548924). We acknowledge computing resources at the Summit supercomputer, a joint effort of the University of Colorado Boulder and Colorado State University, which is supported by the National Science Foundation (ACI-1532235 and ACI-1532236), as well as resources at the Argonne Leadership Computing Facility, which is a DoE Office of Science User Facility supported under Contract DE-AC02-06CH11357.

References

1. Z. Wang and B. Mi, Environ. Sci. Technol., 2017, 51, 8229–8244.
2. H. Wang, C. Li, P. Fang, Z. Zhang and J. Z. Zhang, Chem. Soc. Rev., 2018, 47, 6101–6127.
3. Z. Li, X. Meng and Z. Zhang, J. Photochem. Photobiol., C, 2018, 35, 39–55.
4. Z. Wu, B. Fang, Z. Wang, C. Wang, Z. Liu, F. Liu, W. Wang, A. Alfantazi, D. Wang and D. P. Wilkinson, ACS Catal., 2013, 3, 2101–2107.
5. L. B. Huang, L. Zhao, Y. Zhang, Y. Y. Chen, Q. H. Zhang, H. Luo, X. Zhang, T. Tang, L. Gu and J. S. Hu, Adv. Energy Mater., 2018, 1800734.
6. S. Barua, H. S. Dutta, S. Gogoi, R. Devi and R. Khan, ACS Appl. Nano Mater., 2018, 1, 2–25.
7. S. Wu, J. Wang, L. Jin, Y. Li and Z. Wang, ACS Appl. Nano Mater., 2018, 1, 337–343.
8. B. Dou, J. Yang, R. Yuan and Y. Xiang, Anal. Chem., 2018, 90, 5945–5950.
9. Z. Guan, C.-S. Lian, S. Hu, S. Ni, J. Li and W. Duan, J. Phys. Chem. C, 2017, 121, 3654–3660.
10. Y. Almadori, N. Bendiab and B. Grévin, ACS Appl. Mater. Interfaces, 2017, 10, 1363–1373.
11. G. R. Bhimanapati, T. Hankins, Y. Lei, R. A. Vilá, I. Fuller, M. Terrones and J. A. Robinson, ACS Appl. Mater. Interfaces, 2016, 8, 22190–22195.
12. T. H. Kim, Y. H. Kim, S. Y. Park, S. Y. Kim and H. W. Jang, Chemosensors, 2017, 5, 15.
13. K.-K. Liu, W. Zhang, Y.-H. Lee, Y.-C. Lin, M.-T. Chang, C.-Y. Su, C.-S. Chang, H. Li, Y. Shi and H. Zhang, Nano Lett., 2012, 12, 1538–1544.
14. K. Liu, Q. Yan, M. Chen, W. Fan, Y. Sun, J. Suh, D. Fu, S. Lee, J. Zhou and S. Tongay, Nano Lett., 2014, 14, 5097–5103.
15. X. Tian, D. S. Kim, S. Yang, C. J. Ciccarino, Y. Gong, Y. Yang, Y. Yang, B. Duschatko, Y. Yuan and P. M. Ajayan, Nat. Mater., 2020, 1–7.
16. V. G. Ruiz, W. Liu and A. Tkatchenko, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 93, 035118.
17. P. L. Silvestrelli and A. Ambrosetti, J. Low Temp. Phys., 2016, 185, 183–197.
18. W. Liu, A. Tkatchenko and M. Scheffler, Acc. Chem. Res., 2014, 47, 3369–3377.
19. T. Björkman, A. Gulans, A. V. Krasheninnikov and R. M. Nieminen, Phys. Rev. Lett., 2012, 108, 235502.
20. J. Fuhr, J. Sofo and A. Saúl, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 8343–8347.
21. K. Weiss and J. M. Phillips, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 14, 5392.
22. V. Alexiev, R. Prins and T. Weber, Phys. Chem. Chem. Phys., 2000, 2, 1815–1827.
23. H. Heinz, K. C. Jha, J. Luettmer-Strathmann, B. L. Farmer and R. R. Naik, J. R. Soc., Interface, 2011, 8, 220–232.
24. I. L. Geada, H. Ramezani-Dakhel, T. Jamil, M. Sulpizi and H. Heinz, Nat. Commun., 2018, 9, 716.
25. H. Ramezani-Dakhel, P. A. Mirau, R. R. Naik, M. R. Knecht and H. Heinz, Phys. Chem. Chem. Phys., 2013, 15, 5488–5492.
26. H. Ramezani-Dakhel, L. Y. Ruan, Y. Huang and H. Heinz, Adv. Funct. Mater., 2015, 25, 1374–1384.
27. T. V. Ndoro, E. Voyiatzis, A. Ghanbari, D. N. Theodorou, M. C. Böhm and F. Müller-Plathe, Macromolecules, 2011, 44, 2316–2327.
28. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Chem. Mater., 2014, 26, 2647–2658.
29. T. Z. Lin and H. Heinz, J. Phys. Chem. C, 2016, 120, 4975–4992.
30. R. K. Mishra, L. Fernández-Carrasco, R. J. Flatt and H. Heinz, Dalton Trans., 2014, 43, 10602–10616.
31. D. Roccatano, E. Sarukhanyan and R. Zangi, J. Chem. Phys., 2017, 146, 074703.
32. S. Frankland, A. Caglar, D. Brenner and M. Griebel, J. Phys. Chem. B, 2002, 106, 3046–3048.
33. J. R. Gissinger, C. Pramanik, B. Newcomb, S. Kumar and H. Heinz, ACS Appl. Mater. Interfaces, 2018, 10, 1017–1027.
34. C. Pramanik, J. R. Gissinger, S. Kumar and H. Heinz, ACS Nano, 2017, 11, 12805–12816.
35. H. Heinz, H. J. Castelijns and U. W. Suter, J. Am. Chem. Soc., 2003, 125, 9500–9510.
36. H. Heinz and U. W. Suter, J. Phys. Chem. B, 2004, 108, 18341–18352.
37. H. Heinz, T.-J. Lin, R. K. Mishra and F. S. Emami, Langmuir, 2013, 29, 1754–1765.
38. J. Liu, E. Tennessen, J. Miao, Y. Huang, J. M. Rondinelli and H. Heinz, J. Phys. Chem. C, 2018, 122, 14996–15009.
39. M. Drew, S. Edmondson, G. Forsyth, R. Hobson and P. Mitchell, Catal. Today, 1988, 2, 633–641.
40. T. M. Brunier, M. G. Drew and P. C. Mitchell, Mol. Simul., 1992, 9, 143–159.
41. P. Faye, E. Payen and D. Bougeard, J. Chem. Soc., Faraday Trans., 1996, 92, 2437–2443.
42. U. Becker, K. M. Rosso, R. Weaver, M. Warren and M. F. Hochella Jr, Geochim. Cosmochim. Acta, 2003, 67, 923–934.
43. Y. Morita, T. Onodera, A. Suzuki, R. Sahnoun, M. Koyama, H. Tsuboi, N. Hatakeyama, A. Endou, H. Takaba and M. Kubo, Appl. Surf. Sci., 2008, 254, 7618–7621.
44. T. Onodera, Y. Morita, A. Suzuki, M. Koyama, H. Tsuboi, N. Hatakeyama, A. Endou, H. Takaba, M. Kubo and F. Dassenoy, J. Phys. Chem. B, 2009, 113, 16526–16536.
45. V. Varshney, S. S. Patnaik, C. Muratore, A. K. Roy, A. A. Voevodin and B. L. Farmer, Comput. Mater. Sci., 2010, 48, 101–108.
46. M. Dallavalle, N. Sändig and F. Zerbetto, Langmuir, 2012, 28, 7393–7400.
47. J.-W. Jiang, H. S. Park and T. Rabczuk, J. Appl. Phys., 2013, 114, 064307.
48. T. Liang, S. R. Phillpot and S. B. Sinnott, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 245110.
49. V. Sresht, A. Govind Rajan, E. Bordes, M. S. Strano, A. A. Padua and D. Blankschtein, J. Phys. Chem. C, 2017, 121, 9022–9031.
50. B. Luan and R. Zhou, Appl. Phys. Lett., 2016, 108, 131601.
51. H. Heinz and H. Ramezani-Dakhel, Chem. Soc. Rev., 2016, 45, 412–448.
52. A. Kozbial, X. Gong, H. Liu and L. Li, Langmuir, 2015, 31, 8429–8435.
53. H. Heinz, H. Koerner, K. L. Anderson, R. A. Vaia and B. L. Farmer, Chem. Mater., 2005, 17, 5658–5669.
54. E. Otyepková, P. Lazar, J. Luxa, K. Berka, K. Čépe, Z. Sofer, M. Pumera and M. Otyepka, Nanoscale, 2017, 9, 19236–19244.
55. B. Schönfeld, J. Huang and S. Moss, Acta Crystallogr., Sect. B: Struct. Sci., 1983, 39, 404–407.
56. Z. Yin, H. Li, H. Li, L. Jiang, Y. Shi, Y. Sun, G. Lu, Q. Zhang, X. Chen and H. Zhang, ACS Nano, 2011, 6, 74–80.
57. G. Cunningham, M. Lotya, C. S. Cucinotta, S. Sanvito, S. D. Bergin, R. Menzel, M. S. Shaffer and J. N. Coleman, ACS Nano, 2012, 6, 3468–3480.
58. T. Wieting and J. Verble, Phys. Rev. B: Solid State, 1971, 3, 4286–4292.
59. R. C. Cooper, C. Lee, C. A. Marianetti, X. Wei, J. Hone and J. W. Kysar, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 035423.
60. S. Bertolazzi, J. Brivio and A. Kis, ACS Nano, 2011, 5, 9703–9709.
61. P. Li, Z. You and T. Cui, J. Micromech. Microeng., 2013, 23, 045026.
62. J. Feldman, J. Phys. Chem. Solids, 1976, 37, 1141–1144.
63. A. Webb, J. Feldman, E. Skelton, L. Towle, C. Liu and I. Spain, J. Phys. Chem. Solids, 1976, 37, 329–335.
64. G. Levita, A. Cavaleiro, E. Molinari, T. Polcar and M. C. Righi, J. Phys. Chem. C, 2014, 118, 13809–13816.
65. A. Molina-Sanchez, K. Hummer and L. Wirtz, Surf. Sci. Rep., 2015, 70, 554–586.
66. J. Li, N. V. Medhekar and V. B. Shenoy, J. Phys. Chem. C, 2013, 117, 15842–15848.
67. P. Nicolini and T. Polcar, Comput. Mater. Sci., 2016, 115, 158–169.
68. P. Dauber-Osguthorpe, V. A. Roberts, D. J. Osguthorpe, J. Wolff, M. Genest and A. T. Hagler, Proteins: Struct., Funct., Genet., 1988, 4, 31–47.
69. J. Huang and A. D. MacKerell, J. Comput. Chem., 2013, 34, 2135–2145.
70. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
71. J. M. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
72. W. L. Jorgensen, D. S. Maxwell and J. TiradoRives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
73. H. Sun, S. J. Mumby, J. R. Maple and A. T. Hagler, J. Am. Chem. Soc., 1994, 116, 2978–2987.
74. H. Sun, Macromolecules, 1995, 28, 701–712.
75. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
76. C. C. Dharmawardhana, K. Kanhaiya, T.-J. Lin, A. Garley, M. R. Knecht, J. Zhou, J. Miao and H. Heinz, Mol. Simul., 2017, 43, 1394–1405.
77. K. C. Gross, P. G. Seybold and C. M. Hadad, Int. J. Quantum Chem., 2002, 90, 445–458.
78. Y. T. Fu, G. D. Zartman, M. Yoonessi, L. F. Drummy and H. Heinz, J. Phys. Chem. C, 2011, 115, 22292–22300.
79. S. S. Batsanov, Inorg. Mater., 2001, 37, 871–885.
80. R. Zhang, W. Wan, L. Qiu and Y. Zhou, Mater. Lett., 2016, 181, 321–324.
81. G. D. Zartman, H. Liu, B. Akdim, R. Pachter and H. Heinz, J. Phys. Chem. C, 2010, 114, 1763–1772.
82. T. Livneh and J. E. Spanier, 2D Mater., 2015, 2, 035003.
83. A. P. Gaur, S. Sahoo, M. Ahmadi, S. P. Dash, M. J.-F. Guinel and R. S. Katiyar, Nano Lett., 2014, 14, 4314–4321.
84. D.-M. Tang, D. G. Kvashnin, S. Najmaei, Y. Bando, K. Kimoto, P. Koskinen, P. M. Ajayan, B. I. Yakobson, P. B. Sorokin and J. Lou, Nat. Commun., 2014, 5, 3631.
85. H. Heinz, R. A. Vaia, B. L. Farmer and R. R. Naik, J. Phys. Chem. C, 2008, 112, 17281–17290.
86. L. L. Zhang, B. Q. Luang and R. H. Zhou, J. Phys. Chem. B, 2019, 123, 7243–7252.
87. J. Wang, L. Liu, D. Ge, H. Zhang, Y. Feng, Y. Zhang, M. Chen and M. Dong, Chem. –Eur. J., 2018, 24, 3397–3402.
88. S. Cetinel, W.-Z. Shen, M. Aminpour, P. Bhomkar, F. Wang, E. R. Borujeny, K. Sharma, N. Nayebi and C. Montemagno, Sci. Rep., 2018, 8, 1–12.
89. P. Li, K. Sakuma, S. Tsuchiya, L. Sun and Y. Hayamizu, ACS Appl. Mater. Interfaces, 2019, 11, 20670–20677.
90. C. Muratore, A. Juhl, A. Stroud, D. Wenbi Lai, A. Jawaid, K. Burzynski, J. Dagher, G. Leuty, C. Harsch and S. Kim, Appl. Phys. Lett., 2018, 112, 233704.
91. L. Sun, T. Narimatsu, S. Tsuchiya, T. Tanaka, P. Li and Y. Hayamizu, RSC Adv., 2016, 6, 96889–96897.
92. J. Chen, E. Zhu, J. Liu, S. Zhang, Z. Lin, X. Duan, H. Heinz, Y. Huang and J. J. De Yoreo, Science, 2018, 362, 1135–1139.
93. C. Y. Chiu, Y. J. Li, L. Y. Ruan, X. C. Ye, C. B. Murray and Y. Huang, Nat. Chem., 2011, 3, 393–399.
94. E. Zhu, S. Wang, X. Yan, M. Sobani, L. Ruan, C. Wang, Y. Liu, X. Duan, H. Heinz and Y. Huang, J. Am. Chem. Soc., 2019, 141, 1498–1505.
95. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Chem. Mater., 2014, 26, 5725–5734.
96. R. T. McGibbon, K. A. Beauchamp, M. P. Harrigan, C. Klein, J. M. Swails, C. X. Hernández, C. R. Schwantes, L.-P. Wang, T. J. Lane and V. S. Pande, Biophys. J., 2015, 109, 1528–1532.
97. H. Heinz, B. L. Farmer, R. B. Pandey, J. M. Slocik, S. S. Patnaik, R. Pachter and R. R. Naik, J. Am. Chem. Soc., 2009, 131, 9704–9714.
98. D. Eisenberg, R. M. Weiss, T. C. Terwilliger and W. Wilcox, Faraday Symp. Chem. Soc., 1982, 17, 109–120.
99. E. Pustovgar, R. K. Mishra, M. Palacios, J.-B. d. E. de Lacaillerie, T. Matschei, A. S. Andreev, H. Heinz, R. Verel and R. J. Flatt, Cem. Concr. Res., 2017, 100, 245–262.
100. F. S. Emami, V. Puddu, R. J. Berry, V. Varshney, S. V. Patwardhan, C. C. Perry and H. Heinz, Chem. Mater., 2014, 26, 2647–2658.
101. S. V. Patwardhan, F. S. Emami, R. J. Berry, S. E. Jones, R. R. Naik, O. Deschaume, H. Heinz and C. C. Perry, J. Am. Chem. Soc., 2012, 134, 6244–6256.
102. T. Jamil, A. Javadi and H. Heinz, Green Chem., 2020, 22, 1577–1593.
103. H. Ramezani-Dakhel, N. M. Bedford, T. J. Woehl, M. R. Knecht, R. R. Naik and H. Heinz, Nanoscale, 2017, 9, 8401–8409.
104. J. Zhou, Y. Yang, Y. Yang, D. S. Kim, A. Yuan, X. Tian, C. Ophus, F. Sun, A. K. Schmid, M. Nathanson, H. Heinz, Q. An, H. Zeng, P. Ercius and J. Miao, Nature, 2019, 570, 500–503.
105. R. Xu, C.-C. Chen, L. Wu, M. C. Scott, W. Theis, C. Ophus, M. Bartels, Y. Yang, H. Ramezani-Dakhel, M. R. Sawaya, H. Heinz, L. D. Marks, P. Ercius and J. Miao, Nat. Mater., 2015, 14, 1099–1103.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sc01443e

---