A first principles phase diagram calculation, that included van der Waals interactions, was performed for the 3D bulk system (1X)·MoS2(X)·MoTe2. Surprisingly, the predicted phase diagram has at least two ordered phases, at X0.46, even though all calculated formation energies are positive; in a ground-state analysis that examined all configurations with 16 or fewer anion sites. The lower-temperature I-phase is predicted to transform to a higher-temperature I-phase at T500K, and I disorders at T730K. Both these transitions are predicted to be first-order, and there are broad two-phase fields on both sides of the ordered regions. Both the I- and I-phases are predicted to be incommensurate, i.e., aperiodic: I-phase in three dimensions; and I-phase in two dimensions.

Recently there has been great interest in two-dimensional (2D) transition metal dichalcogenide (TMD) materials such as MoS2, MoSe2, and MoTe2, their solid solutions, and related 2D materials.1,2 Traditionally, MoS2 has been used as a dry lubricant3 that is stable up to 623 K. Currently, interest is focused on applications as: band-gap engineering materials;4,5 nano-electronic devices;2,6–8 photovoltaic devices;9,10 valleytronics applications;11,12 2D building blocks for electronic heterostructures;13 and as sensor materials.14 

The individual, three-atom-thick, 2D-layers of the bulk system are bonded by van der Waals forces, hence these forces influence the bulk and multilayer synthesis and therefore anion order-disorder and/or phase separation in solid solutions. The results presented below, for 3D bulk MoS2MoTe2, imply that van der Waals interactions may strongly affect phase stabilities, either between adjacent layers in bulk or few-layer samples, or between monolayers and heterogeneous substrates.

Of the three quasibinary solid solutions (MoS2MoSe2, MoSe2MoTe2,MoS2MoTe2) (1X)·MoS2(X)·MoTe2 has the greatest difference in anionic radii (RS = 1.84 Å; RTe = 2.21 Å),15 which suggests that it is the most likely to exhibit interesting solution behavior. One expects a simple miscibility gap as reported by Kang et al.4 for monolayer MoS2MoTe2, hence the prediction of two configurational entropy (Scon) stabilized incommensurate, i.e., aperiodic, phases is extraordinary (stable phases that have positive formation energies must be entropy stabilized).

Total structure energies, ΔEStr were calculated for fully relaxed MoS2, MoTe2 (2H-structure, space group P63/mmc, AB-stacking of three-atom-thick layers), and for 233 Mom+n(SmTen)2 supercells. The Vienna abinitio simulation program (VASP, version 5.3.316,17) was used for all density-functional theory (DFT) calculations, with projector augmented waves (PAW) and a generalized gradient approximation (GGA) for exchange energies. Electronic degrees of freedom were optimized with a conjugate gradient algorithm. Valence electron configurations were: Mo_pv 4p55s4d; S_s2p4; Te_s2p4. Van der Waals interactions that bond the three-atom thick 2D X-Mo-X layers (X = S, Te) together were modeled with the non-local correlation functional of Klimes et al.18 Total energies were also calculated without van der Waals interactions, but, up to a basis of 140 structures, ground-state analyses always predicted false ground-states (supplementary material). Convergence with respect to k-point meshes was achieved by increasing the number of k-points until the total energy converged. A 500 eV cutoff-energy was used in the “high precision” option, which converges absolute energies to within a few meV/mol (a few tenths of a kJ/mol of exchangeable S- and Te-anions). Precision is at least an order of magnitude better. Residual forces of order 0.02 eV or less were typical. Often, convergence with respect to hexagonal c-axis length was not automatic, and it was necessary to choose an initial c-axis value that is close to the converged value. Calculated interlayer spacings in MoS2 and MoTe2 are 2.992 Å and 3.513 Å, respectively, corresponding experimental values are: 2.977 Å (Ref. 19) and 3.382 Å.20 

Formation energies (ΔEf) for 233 Mom+n(SmTen)2 supercells are plotted in Fig. 1, in which values for ΔEf are normalized per mol of exchangeable anions, S and Te

ΔEf=(EStrmEMoS2nEMoTe2)/(2(m+n)).
(1)

Here EStr is the total energy of the Mom+n(SmTen)2 supercell; EMoS2 is the energy/mol of MoS2; EMoTe2 is the energy/mol of MoTe2.

FIG. 1.

Comparison of formation energies, ΔEf, for the 235 DFT calculations (solid circles, green online) to Cluster Expansion (CE) formation energies: ΔEFit (large open squares, red online) is the CE-fit to the DFT set; ΔEGS (smaller open squares, blue online) are the CE-based ground-state analysis; ΔEI (solid black diamond at X = 0.46, ΔEI0.03 eV) is the I-phase formation energy. All ΔEf > 0 implies that there are no ordered ground-states, and suggests that the phase diagram will have a miscibility gap.

FIG. 1.

Comparison of formation energies, ΔEf, for the 235 DFT calculations (solid circles, green online) to Cluster Expansion (CE) formation energies: ΔEFit (large open squares, red online) is the CE-fit to the DFT set; ΔEGS (smaller open squares, blue online) are the CE-based ground-state analysis; ΔEI (solid black diamond at X = 0.46, ΔEI0.03 eV) is the I-phase formation energy. All ΔEf > 0 implies that there are no ordered ground-states, and suggests that the phase diagram will have a miscibility gap.

Close modal

All supercell energies are positive which suggests a miscibility gap system, unless one or more entropy stabilized phases are stable.

A cluster expansion Hamiltonian (CEH),21 for the (1-X)· MoS2–(X)· MoTe2 quasibinary system was fit to the set of 235 formation energies, ΔEVASP, solid dots (green online) in Fig. 1 with a cross validation score of (CV)2 = 0.00723896. Fitting of the CEH was performed with the Alloy Theoretic Automated Toolkit (ATAT)17,22–25 which automates most of the tasks associated with CEH construction. A complete description of the algorithms underlying the code can be found in Ref. 23. Large open squares in Fig. 1 (red online) indicate values of the 235 ΔEFit that were calculated with the CEH. Smaller open squares (ΔEGS, blue online) indicate the results of a ground-state analysis in which the CE was used to calculate the formation energies for all ordered configurations with 16 or fewer anion sites, 151 023 structures.

A first principles phase diagram (FPPD) calculation was performed with grand-canonical, and canonical, Monte Carlo (MC) simulations using the emc2 and phb codes which are part of the ATAT package.22–24 Most phase boundaries were calculated with the phb program which uses equilibration tests to set the numbers of equilibration- and MC-passes.24 To draw high-T extensions of the two-phase fields, and to locate the II transition, a 48 × 48 × 12 unit cell simulation box was used, with 2000 equilibration passes and 2000 MC-passes (see supplementary material for comparisons of various equilibration- and MC-pass settings in calculations of the II phase transition). The predicted phase diagram is shown in Fig. 2; where (MoS2) denotes the MoS2-rich solution phase, and similarly for (MoTe2).

FIG. 2.

Calculated phase diagram. The isothermal line at 82.5 K indicates that the I-phase, at X0.46, is entropy stabilized (not a ground-state); Here: (MoS2) and (MoTe2) indicate MoS2-rich and MoTe2-rich solid solutions, respectively. Dotted lines in the region of the Idisordered transition indicate inferred extensions of calculated phase boundaries.

FIG. 2.

Calculated phase diagram. The isothermal line at 82.5 K indicates that the I-phase, at X0.46, is entropy stabilized (not a ground-state); Here: (MoS2) and (MoTe2) indicate MoS2-rich and MoTe2-rich solid solutions, respectively. Dotted lines in the region of the Idisordered transition indicate inferred extensions of calculated phase boundaries.

Close modal

Kang et al.4 performed first principles phase diagram calculations for four dichalcogenide monolayer systems: MoSe2(1x)Te2x,WSe2(1x)Te2x, MoS2(1x)Te2x, and WS2(1x)Te2x. Their CEs were fit to about 40 structures per system, and van der Waals interactions (with substrate) omitted. All systems were predicted to have miscibility gaps, and surprisingly, all consolute points are on the Te-rich sides. One expects the consolute point to be on the S-rich side, because it typically requires less energy to substitute a smaller S-ion into a larger Te-ion site, than vice versa. Figure 2 also has reduced solubility on the Te-rich side, but this is related to immiscibility between the I- (I)-phase, and the Te-rich phase.

It is not clear that the Kang et al. MoS2MoTe2 phase diagram would still be a simple miscibility gap had they included hundreds of structures in their CE-fit rather than about 40. Hence the monolayer vs. bulk comparison is uncertain.

Surprisingly, the phase diagram predicted here has multiple two-phase fields, separated by two ordered incommensurate phases, neither of which is a ground state. To seven digits, the calculated bulk composition of the I-phase, just above its 82.5 K minimum temperature of stability, is X = 0.4642857 = 13/28; i.e., Mo14S15Te13. Stability of the I-phases is a robust result, in Monte-Carlo simulations: (1) CEH fits to 128, 153, 162, 182, 225, and 235 formation energies all predict I-type ordering; (2) I-phase forms spontaneously on cooling of the I-phase at T500K, and on heating of a low-T equilibrium (MoS2)+(MoTe2) assemblage to T82K; (3) I-phase forms spontaneously on heating of the I-phase, or cooling of a disordered solid solution with X0.46. This calculation considers only Scon, and ignores excess vibrational entropy, Svib, which could conceivably destabilize the I-phases. In light of (1) above, however, this seems highly unlikely. Also, there is no fully satisfactory way to model Svib for an aperiodic phase, and a reasonable approximate structure (with I-phase like ordering) would require at least a low symmetry 84-atom cell; which is beyond the scope of this study.

Figure 3 shows how the CE-calculated I-phase formation energy ΔEf(Iphase)ΔEI, varies as a function of MC-supercell size and shape. The flat minimum at ΔEI0.03 eV/anion is the calculated I-phase formation energy which is plotted as the solid black diamond in Fig. 1. Supercell dimensions were chosen to accommodate Mo14S15Te13 stoichiometry and I-phase ordering. Note that many of the ΔEf plotted in Fig. 1, are lower in energy than ΔEI, but that they are for periodic structures with 16 or fewer anion sites in which Scon0 as T0K, and clearly (Fig. 4) the I-phase is incommensurate, i.e., aperiodic with Scon> 0. Figure 4 exhibits the S:Te (yellow:brown online, respectively) ordering at: (a) 200 K; and (b) 575 K, i.e., below and above the II phase transition. Mo-atoms are omitted for clarity, and labels (001)D, (100)D, and (010)D refer to corresponding crystallographic planes in the high-T P63/mmcdisordered-phase.

FIG. 3.

Minimum energies for various Monte-Carlo supercells; Nz is the length of the supercell, in c-axis units of the P63/mmc disordered-phase cell constants. The flat minimum at ΔEfΔEI0.03 eV/anion is interpreted as the I-phase formation energy.

FIG. 3.

Minimum energies for various Monte-Carlo supercells; Nz is the length of the supercell, in c-axis units of the P63/mmc disordered-phase cell constants. The flat minimum at ΔEfΔEI0.03 eV/anion is interpreted as the I-phase formation energy.

Close modal
FIG. 4.

Monte-Carlo-snapshots of S:Te-ordering in the I- and I-phases at: (a) X0.46 and T = 200 K; and (b) X0.475 and T = 575 K, respectively (online S = yellow, Te = brown, Mo omitted for clarity). Labels (001)D, (100)D, and (010)D, refer to corresponding crystallographic planes in the high-T P63/mmc disordered phase.

FIG. 4.

Monte-Carlo-snapshots of S:Te-ordering in the I- and I-phases at: (a) X0.46 and T = 200 K; and (b) X0.475 and T = 575 K, respectively (online S = yellow, Te = brown, Mo omitted for clarity). Labels (001)D, (100)D, and (010)D, refer to corresponding crystallographic planes in the high-T P63/mmc disordered phase.

Close modal

In the (001)D- and (100)D-planes,. SmTen. chains in the ⟨010D direction, most often have m=5or6 and n=4or5. Also, in (100)D, the SmTen chains exhibit irregular alignments relative to one another. Note however, that SmTen-chains in (001)D-planes order along the ⟨001D direction, such that Sm-units alternate with Ten-units in adjacent 3-atom thick 2D-layers; i.e., …SmTenSmTen… chains are stacked on top of …TenSmTenSm… chains with inescapable misfits, owing to the different and variable values of m and n. Thus I-phase ordering is inevitably imperfect, aperiodic, and incommensurate, which suggests that the II phase transition is first-order.

The difference between I- and I-phases appears to be a distinction between 3D-ordering in the low-T I-phase and 2D-ordering in the high-T I-phase. Clearly, ordering is stronger in the basal (001)D-plane than in the (100)D- or (010)D-planes, and striped order within (001)D persists above the II transition.

Figure 5(a) is a Monte-Carlo T-scan (heating) of the total energy ETOT(T) which confirms first-order character for the II transition: a critical (continuous) transition26 would not exhibit the sharp change at T503K; also a transition from the lower-T, higher-ETOT phase to the higher-T, lower-ETOT phase requires that the lower-T phase be superheated, i.e., metastable before it transforms. Figure 5(b) is an idealized schematic that compares equilibrium- and metastable-transition paths (solid black line, and blue line with red dots, online respectively). This transition is more subtle in cooling simulations, but still evident in snapshots. In the Fig. 5 inset, no change of slope at the transition is evident in the Helmholtz energy, F(T), which suggests that the transition is weakly first-order.27 Also, when the MC temperature-increment was decreased from 1.0 K/MC-step to 0.1 K/MC-step (not shown), the predicted transition temperature decreased by about 10 K, which indicates more superheating at 1.0 K/MC-step than at 0.1 K/MC-step, hence the first-order character. This transition is shown as a dotted line in Fig. 2 because the two-phase fields that a first-order transition implies are too narrow to resolve in the MC-simulations.

FIG. 5.

(a) Monte-Carlo T-scans (heating) of: total energy, ETOT(T), and Helmholtz energy, F(T), (inset) as functions of temperature; (b) idealized schematic comparing equilibrium- and metastable transition paths (solid black line and red-dotted blue line, online, respectively). ETOT(T) indicates a phase transition at T503K; the absence of a clear change of slope in F(T) indicates that the transition is only weakly first-order; hence the dotted II transition line in Fig. 2.

FIG. 5.

(a) Monte-Carlo T-scans (heating) of: total energy, ETOT(T), and Helmholtz energy, F(T), (inset) as functions of temperature; (b) idealized schematic comparing equilibrium- and metastable transition paths (solid black line and red-dotted blue line, online, respectively). ETOT(T) indicates a phase transition at T503K; the absence of a clear change of slope in F(T) indicates that the transition is only weakly first-order; hence the dotted II transition line in Fig. 2.

Close modal

To summarize, a first principles phase diagram calculation for the 3D bulk system (1X)·MoS2(X)·MoTe2, that includes van der Waals interactions, predicts the formation of two entropy stabilized incommensurate, i.e., aperiodic phases: the I- and I-phases, at X 0.46. Above the minimum temperature for stability of the I-phase, T82K, the calculation predicts broad two-phase fields between the I- or I-phase and disordered S- or Te-rich solution phases, (MoS2) and (MoTe2), respectively. Both the II and Idisordered transitions are predicted to be first-order. Dramatic changes in phase relations can be induced by arbitrarily small differences in energy, hence van der Waals interactions should not be ignored in layered 2D-systems such as MoS2MoTe2.

See supplementary material for the ground state analysis of cluster expansions based on energies estimated without including van der Waals interactions, calculated transition points obtained from MC simulations for the high-T extensions of the two-phase fields and the comparisons of various equilibration- and MC-pass settings in calculations of the II phase transition.

This work was supported by NIST-MGI, and B. Burton wishes to thank R. Selinger, J. Douglas, K. Migler, A. Davydov, B. Campbell, J. Lau, and L. Bendersky for useful discussions.

1.
Q. H.
Wang
,
K.
Kalantar-Zadeh
,
A.
Kis
,
J. N.
Coleman
, and
M. S.
Strano
,
Nat. Nanotechnol.
7
,
699
(
2012
).
2.
R.
Ganatra
and
Q.
Zhang
,
ACS Nano
8
(
5
),
4074
(
2014
).
3.
F. L.
Claus
,
Solid Lubricants and Self-Lubricating Solids
(
Academic Press
,
New York
,
1972
);
W.
Mller-Warmuth
and
R.
Schllhorn
,
Progress in Intercalation Research
(
Springer
,
1994
), ISBN 0-7923-2357-2.
4.
J.
Kang
,
S.
Tongay
,
J.
Li
, and
J.
Wu
,
J. Appl. Phys.
113
,
143703
(
2013
).
5.
A.
Kutana
,
E. S.
Penev
, and
B. I.
Yakobson
,
Nanoscale
6
,
5820
(
2014
).
6.
B.
Radisavljevic
,
A.
Radenovic
,
J.
Brivio
,
V.
Giacometti
, and
A.
Kis
,
Nat. Nanotechnol.
6
(
3
),
147
(
2011
).
7.
S.
Das
,
H.-Y.
Chen
,
A. V.
Penumatcha
, and
J.
Appenzeller
,
Nano Lett.
13
(
1
),
100
(
2013
).
8.
H.
Wang
,
L.
Yu
,
Y.-H.
Lee
,
Y.
Shi
,
A.
Hsu
,
M. L.
Chin
,
L.-J.
Li
,
M.
Dubey
,
J.
Kong
, and
T.
Palacios
,
Nano Lett.
12
(
9
),
4674
(
2012
).
9.
D.
Jariwala
,
V. K.
Sangwan
,
L. J.
Lauhon
,
T. J.
Marks
, and
M. C.
Hersam
,
ACS Nano
8
(
2
),
1102
(
2014
).
10.
M.
Fontana
,
T.
Deppe
,
A. K.
Boyd
,
M.
Rinzan
,
A. Y.
Liu
,
M.
Paranjape
, and
P.
Barbara
,
Sci. Rep.
3
,
1634
(
2013
).
11.
H.
Zeng
,
J.
Dai
,
W.
Yao
,
D.
Xiao
, and
X.
Cui
,
Nat. Nanotechnol.
7
(
8
),
490
(
2012
).
12.
K. F.
Mak
,
K.
He
,
J.
Shan
, and
T. F.
Heinz
,
Nat. Nanotechnol.
7
(
8
),
494
(
2012
).
13.
A. K.
Geim
and
I. V.
Grigorieva
,
Nature
499
,
419
(
2013
).
14.
S. X.
Wu
,
Z. Y.
Zeng
,
Q. Y.
He
,
Z. J.
Wang
,
S. J.
Wang
,
Y. P.
Du
 et al.
Small
8
,
2264
(
2012
);
[PubMed]
T. Y.
Wang
,
H. C.
Zhu
,
J. Q.
Zhuo
,
Z. W.
Zhu
,
P.
Papakonstantinou
,
G.
Lubarsky
 et al 
Anal. Chem.
85
,
10289
(
2013
).
[PubMed]
15.
See http://abulafia.mt.ic.ac.uk/shannon/ptable.php for database of elemental data including atomic and ionic radii.
16.
G.
Kresse
and
J.
Hafner
,
Phys. Rev.
47
,
558
(
1993
);
G.
Kresse
, Ab initio Molekular Dynamik für flüssige Metalle thesis,
Technische Universität Wien
,
1993
;
G.
Kresse
,
Phys. Rev.
B49
(
14
),
251
(
1994
);
G.
Kresse
and
J.
Furthmüller
,
Comput. Mat. Sci.
6
,
15
50
(
1996
);
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
);
17.
Reference to specific software packages does not imply a NIST endorsement.
18.
J.
Klimes
,
D. R.
Bowler
, and
A.
Michaelides
,
Phys. Rev. B
83
,
195131
(
2011
);
J.
Klimes
,
D. R.
Bowler
, and
A.
Michaelides
,
J. Phys. Condens. Matter
22
,
022201
(
2010
).
[PubMed]
19.
Y.
Takeuchi
and
W.
Nowacki
,
Schweiz. Miner. Petrogr. Mitt
44
,
105
(
1964
).
20.
P.
Villars
and
L. D.
Calvert
,
Pearson's Handbook of Crystallographic Data for Intermetallic Phases
(
American Society for Metals
,
Metals Park, OH
,
1985
), Vol. 1, p.
44073
.
21.
J. M.
Sanchez
,
F.
Ducastelle
, and
D.
Gratias
,
Physica
128
,
334
(
1984
).
22.
A.
van de Walle
,
M.
Asta
, and
G.
Ceder
,
CALPHAD J.
26
,
539
(
2002
).
23.
A.
van de Walle
and
G.
Ceder
,
J. Phase Equilib.
23
,
348
(
2002
).
24.
A.
van de Walle
and
M.
Asta
,
Modell. Simul. Mater. Sci. Eng.
10
,
521
(
2002
).
25.
A.
van de Walle
and
G.
Ceder
,
Rev. Mod. Phys.
74
,
11
(
2002
).
26.
M. E.
Fisher
,
Rev. Mod. Phys.
70
(
2
),
653
(
1998
).
27.

The phrase “weakly first-order” is used to describe a transition with a sufficiently small discontinuity in order parameter that some first-order characteristics are undetectable: e.g., a kink in F(T), or I + I′ two-phase fields.

Supplementary Material