ARTICLE

pubs.acs.org/JPCC

A Comparative Study of Lattice Dynamics of Three- and

Two-Dimensional MoS2
C. Ataca,†,‡ M. Topsakal,‡ E. Akt€urk,‡ and S. Ciraci*,†,‡
†
Department of Physics and ‡UNAM—Institute of Materials Science and Nanotechnology, Bilkent University, Ankara 06800, Turkey

ABSTRACT: This paper presents a comparative study of the

lattice dynamics of three-dimensional layered MoS2 and two-
dimensional single layer MoS2 based on the density functional
theory. A comprehensive analysis of energetics and optimized
structure parameters is performed using different methods. It
is found that the van der Waals attraction between layers of
three-dimensional (3D) layered MoS2 is weak but is essential
to hold the layers together with the equilibrium interlayer
spacing. Cohesive energy, phonon dispersion curves, and
corresponding density of states and related properties, such
as Born-effective charges, dielectric constants, Raman and infrared active modes are calculated for 3D layered as well as 2D single
layer MoS2 using their optimized structures. These calculated values are compared with the experimental data to reveal
interesting dimensionality effects. The absence of a weak interlayer interaction in 2D single layer MoS2 results in the softening of
some of Raman active modes.

’ INTRODUCTION Despite the fact that 2H-MoS2 is a layered material, where

Three-dimensional (3D) MoS2, a well-known transition- MoS2 layers were bound by weak interlayer interaction, sig-
metal dichalcogenide, has two different stable structures: 3R- nificant dimensionality effects have been observed. For exam-
MoS2 polytype1 and layered 2H-MoS2. The latter consists of ple, while 3D MoS2 is an indirect band gap semiconductor, the
the stacking of MoS2 layers and is the subject matter of the band gap increases and becomes direct in 2D single layer
present study. Various properties of 2H-MoS2 (refs 2
17), in MoS2.20 This dimensionality effect may lead to photolumines-
particular lattice dynamics and electronic energy band struc- cence applications in nanoelectronics.47 While the lattice
ture, have been studied extensively. Recently 2D suspended dynamics of 2H-MoS2 have been studied actively in the past
single layer MoS2 sheets, i.e., 1H-MoS2 having hexagonal lattice by using inelastic neutron scattering and Raman
infrared
have been produced.18
21 Single layer MoS2 nanocrystals of spectroscopy6,24 and its phonon dispersion curves, phonon
∼30 Å width were also synthesized on the Au(111) surface22 density of states, infrared and Raman active modes are calcu-
and the first direct real space STM images of single layer MoS2 lated in terms of force constants derived from experimental
nanosheets have been reported. In the meantime, theoretical data, yet an ab initio treatment including van der Waals
studies (refs 16 and 23
30) on 1H-MoS2 have appeared. interaction (vdW) is absent. Recent papers21,46 investigating
Three-dimensional 2H-MoS2 and 2D 1H-MoS2,30 quasi-1D the Raman spectra of 3D and 2D MoS2 came up with conflicting
nanotubes,25 and nanoribbons16,26,29 of MoS2 share the hon- conclusions. In this paper, we present our theoretical investiga-
eycomb structure and display interesting dimensionality effects. tion of the lattice dynamics and related properties of 2H- and
Properties of MoS2 nanocrystals are explored in diverse 1H-MoS2. Our study is carried out from the first principles
fields, such as nanotribology,31 hydrogen production,32,33 within van der Waals (vdW) density functional theory (DFT),
hydrodesulfurization catalyst used for removing sulfur com- where atomic structure, lattice constants, and relevant ener-
pounds from oil,34
40 solar cells,41 and photocatalysis.42 getics are obtained from structure optimization including vdW
Triangular MoS2 nanocrystals of diverse sizes were investi- correction. This method provides a proper treatment of the
gated using atom-resolved scanning tunneling microscopy.43 interaction between the layers of 2H-MoS2, as well as their
A superlow coefficient of sliding friction between surfaces coated spacings. Finally, calculated properties of 2H- and 1H-MoS2,
with 1H-MoS2 has been measured recently.44 A transistor such as Raman and infrared active modes, bulk modulus,
fabricated from the single layer MoS2 has heralded the features dielectric constants, and effective charges are compared to
of 1H-MoS2, which is superior to graphene.45 Studies to date reveal dimensionality effects. Even though the calculated
suggest that MoS2 sheets can be promising for optoelectronic
devices, solar cells, and LEDs. Most recently, the Raman spectra Received: May 31, 2011
of MoS2 sheets have been measured as a function of their Revised: July 18, 2011
thickness.21,46 Published: July 22, 2011

r 2011 American Chemical Society 16354 dx.doi.org/10.1021/jp205116x | J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

consists of a monatomic Mo plane between two monatomic S

planes having the same 2D hexagonal lattice. Mo and S2 occupy
alternating corners of a hexagon to form a honeycomb struc-
ture. Each Mo has six nearest S atoms, and each S atom has three
nearest Mo atoms. Two layers in the unit cell are displaced
relative to each other, as such that Mo atoms of one layer are
situated on top of S atoms in two adjacent layers. In this respect,
the arrangements of layers are different from that of graphite,
where three carbon atoms of one layer are located above the
hollow sites (center of the hexagon) of the adjacent layers. The
Figure 1. (a) Side and top views of atomic structure of 2H-MoS2 with atomic configurations and relevant structural parameters of 3D
hexagonal lattice. The unit cell is delineated, lattice constants |a| = |b|, c 2H-MoS2 are illustrated in Figure 1a. First-principles calcula-
and internal structure parameters are indicated. Honeycomb structure tions of lattice dynamics and related properties of 3D and 2D
consisting of Mo (red ball) and S2 (gray balls) located at the corners of MoS2 are obtained by using the optimized atomic structure.
hexagons is seen in the top view. (b) Corresponding Brillouin zone with Therefore, we start to investigate the energetics, minimum
symmetry directions. energy atomic structure, and structure parameters of these
crystals.
interlayer interaction is a weak vdW, its effects on lattice
dynamics are significant. We showed that the calculated shifts
of the frequencies of Raman active modes as a result of phonon ’ OPTIMIZED STRUCTURE
softening upon lowering the dimensionality from 3D to a 2D The interaction between 1H-MoS2 layers in 2H-MoS2 has
single layer MoS2 are sensitive to the deviation of lattice predominantly vdW character. Therefore, DFT calculations
constants from the experimental values. within GGA but without vdW are known to overestimate the
interlayer distance and the lattice parameter c in Figure 1a. To
’ METHOD present a correct estimation of lattice constants, we included
We carried out first-principles plane wave calculations with- vdW correction to GGA calculations using two different
in DFT and used ultrasoft pseudopotential (UP).48 The methods. The first one (GGA+D) is used mainly for molecules
exchange correlation potential is approximated by a general- and corrects GGA by adding interatomic C6/R6 interaction.
ized gradient approximation (GGA) using the PW9149 func- The C6 coefficient and cutoff distance are deduced from
tional for both spin-polarized and spin-unpolarized cases. The relevant molecules.50 The second method (GGA+DF) aims
vdW corrections50,51 are included to determine the interlayer at solutions from the first principles without empiricism and
spacing in 2H-MoS2. While all discussions in the paper are uses nonlocal exchange-correlation functional to treat vdW
based on the results obtained within GGA combined with vdW interaction in GGA.51 The latter method is tested for molecules
interaction, the calculations within local density approxi- and solids. In order to find the most appropriate method, we
mation52 (LDA) using UP are also performed for the purpose carried out GGA calculations without vdW correction, as well as
of comparison. All structures are treated using the periodic GGA+D and GGA+DF calculations. For the sake of compar-
boundary conditions. Kinetic energy cutoff, Brillouin zone ison, we also performed LDA calculations, which is known to
(BZ) sampling is determined after extensive convergence include vdW interaction partially.56,57 In our analysis we con-
analysis. A large spacing of ∼10 Å between 2D single layers sider two layered 3D crystals both having honeycomb struc-
of MoS2 is taken to prevent interactions between them. A ture, namely, graphite and 2H-MoS2, in which the cohesions of
plane-wave basis set with kinetic energy cutoff of 952 eV is used layers are known to be achieved mainly by the weak vdW
for high accuracy of the results.53 In the self-consistent field interaction. The optimized lattice constants |a| = |b| and c, and
potential and total energy calculations BZ is sampled by special interlayer interaction energy Ei calculated by using these
k-points.54 The numbers of these k-points are (15  15  7) methods are presented in Table 1. For comparison the experi-
for 2H-MoS2 and (25  25  1) for 1H-MoS2, respectively. All mental values are also given. The interlayer interaction energy
atomic positions and lattice constants are optimized by using or cohesion of 2H-MoS2 relative to individual MoS2 layers can
the conjugate gradient method, where the total energy and be deduced by calculating the total energy of 2H-MoS2 as a
atomic forces are minimized. The convergence for energy is function of interlayer spacing z in the direction perpendicular to
chosen as 10
7 eV between two consecutive steps, and the the MoS2 layer, i.e., ET(z) and by setting ET(zf∞) = 0. Then
maximum Hellmann
Feynman force acting on each atom is the absolute value of half of the minimum of ET(z = c) is taken
less than 10
5 eV/Å upon ionic relaxation. The pressure in the as Ei (interlayer interaction energy per layer).
unit cell is kept below 1 kbar. Numerical calculations have been For graphite, GGA without vdW attains the experimental a,
performed by using the PWSCF package.55 The phonon but overestimates c by 25% relative to the experimental
dispersion curves are calculated within density functional value.58
62 Expectantly, the calculated interlayer interaction
perturbation theory (DFPT) using plane wave methods as energy Ei = 5 meV per layer is much smaller than the
implemented in the PWSCF package. experimental value.58
62 On the other hand, the values of c
and Ei calculated for graphite are improved significantly when
vdW correction is included. While GGA+DF overestimates c
’ ENERGETICS AND OPTIMIZED STRUCTURES OF 3D by 6.6%, GGA+D underestimates it by 4.2%. However, both
AND 2D MOS2 methods result overestimate Ei by 14% relative to experimen-
Three dimensional bulk 2H-MoS2 is a stable layered struc- tal value. Interestingly, LDA yield almost the experimental
ture and includes two MoS2 units in the unit cell. Each layer value of c, even if it underestimates the experimental Ei by 8%.
16355 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

Table 1. Lattice Constants |a| = |b| and c, and Interlayer Interaction Energy Ei (per graphene (C2) or MoS2 unit) of Graphite and
2H-MoS2 Calculated Using GGA, GGA+D and GGA+DF, and LDA Methods.a

graphite (graphene) 2H-MoS2 (1H-MoS2)

a (Å) c (Å) Ei (eV/C2, kcal/mol) a (Å) c (Å) Ei (meV/MoS2, kcal/mol)

GGA 2.461 (2.463) 8.407 5, 0.115 3.215 (3.214) 15.540 6, 0.138

GGA+D 2.461 (2.463) 6.425 122, 2.816 3.220 (3.220) 12.411 160, 3.693
GGA+DF 2.463 (2.463) 7.157 116, 2.678 3.258 (3.254) 13.152 176, 4.063
LDA+UP 2.441 (2.441) 6.669 96, 2.216 3.125 (3.118) 12.137 110, 2.539
experiment 2.461
2.46358
62 (2.455) 6.708
6.71258
62 104 ( 10,58,62 3.1610,63,64 12.29,10,63 140 ( 22,67
2.401 ( 0.231 65
(3.20, 3.27 ) 66
12.30 64
3.239 ( 0.498
a
The corresponding values calculated for graphene and single layer MoS2 are given in parentheses. Experimental values are given for the sake of
comparison. Experimental values of lattice constant a of 1H-MoS2 given by refs 65 and 66 appear to be too large and not confirmed.

Table 2. Lattice Constants |a| = |b| and c, Internal Structure Parameters, Such as Bond Lengths dMo
S and dS
S, S
Mo
S Bond
Angles θ (S
Mo
S), Bulk Modulus B0, in-Plane Stiffness C for 2D Single Layer, Cohesive Energy per MoS2 Unit EC, Born
Effective Charges of Constituent Atoms ZB*[Mo] and ZB*[S], High Frequency Intralayer and Interlayer Dielectric Constants,
ε and ε^ Calculated for 2H-MoS2 and 1H-MoS2 Using the GGA+D Method
a (Å) c (Å) dMo
S (Å) dS
S (Å) Θ B0 (GPa)/C (N/m) EC (eV) ZB*(Mo) ZB*(S) ε ε^

)
2H-MoS2 3.220 12.411 2.436 3.150 80.564 44 15.316 1.23
0.57 15.60 6.34
1H-MoS2 3.220 2.437 3.153 80.617 145.82 15.156 1.21
0.57 4.58 1.26

An analysis made for 2H-MoS2 reveals similar trends and 2H-MoS2 by fitting the Murnaghan equation68 as 44 GPa.
indicates vdW interaction as the dominant interaction between The experimental value69 is given as 43 GPa. Using van der
its layers. While interlayer interaction calculated with GGA is Waals included DFT Rydberg et al.11 calculated B0 as 39 GPa.
only 6 meV per MoS2 unit, c is badly overestimated to be Our value calculated for the bulk modulus of B0 is in good
15.54 Å. This is 26.4% longer relative to the experimental value agreement with the experimental value.
of the lattice constant c measured10,63,64 to be c = 12.29
12.30 Å. Single layer, 1H-MoS2 has high planar strength but trans-
While GGA optimizes lateral lattice constants at |a| = |b| = versal flexibility. While Young’s modulus normally characterizes
3.215 Å, the measured lateral lattice constants |a| = |b| = the mechanical strength of bulk materials, owing to the ambi-
3.16 Å.10,63,64 In contrast, GGA+D and GGA+DF methods guities in defining the Young’s modulus for the 2D honeycomb
estimate Ei to be 160 and 176 meV, respectively. Accordingly, c structure, one can use in-plane stiffness C = (1/A0)(∂2ES/∂ε2s )
values calculated by these vdW corrections are 12.41 Å (i.e., in terms of the equilibrium area of the 2D cell, A0.70,71 We focused
overestimated by 0.9%) and 13.15 Å (overestimated by 7.0%), our attention on the harmonic range of the elastic deformation,
respectively. The lateral lattice constant a is optimized by GGA where the structure responded to strain ε linearly. Here εs is the
+D to be ∼3.22 Å (3.258 Å by GGA+DF). On the other hand, elongation per unit length. The strain energy is defined as ES =
LDA underestimates both a and c by ∼2% relative to experi- ET(εs)
ET(εs = 0); namely, the total energy at a given strain εs
mental value and predicts Ei = 110 meV. In view of this analysis minus the total energy at zero strain. The calculated in-plane
and comparison with measured10,63,64 lattice constants, the stiffness of 1H-MoS2 is 145.82 N/m. This value can be compared
GGA+D method appears to be suitable to obtain optimized with the experimental value of graphene, i.e., 340 ( 50 N/m.72
structure and related properties of 2H-MoS2 and 1H-MoS2. For Sun et al.17 obtained high-frequency dielectric constants, ε, and
the rest of the paper, we will use results obtained from this Born effective charge, ZB*, of 2H-MoS2 by fitting to the experi-
method unless it is stated otherwise. mental data. They found Born effective charges, ZB*[Mo] = 1.11
electrons (positive charge) for Mo and ZB*[S] =
0.52 electrons
’ COHESIVE ENERGY AND OTHER PROPERTIES (negative charge) for each S atom and dielectric constants ε =
)

The cohesive energy per MoS2 unit relative to the free Mo 15.2 and ε^ = 6.2 in the intralayer and interlayer directions,
and S atoms, EC =
ET[Crystal] + ET[Mo] + 2ET[S], is respectively. Here we calculate Born effective charges and high-
calculated from the structure optimized total energies of the frequency dielectric constant of 2H-MoS2 to be ZB*[Mo] = 1.23
3D crystal, ET[2H-MoS2]/2 or 2D crystal ET[1H-MoS2] and electrons and ZB*[S] =
0.57 electrons.73 High-frequency di-
the free atom total energies of Mo and S, ET[Mo] and ET[S], electric constants are calculated to be, ε = 15.60 and ε^ = 6.34 in
)

respectively. The calculated values of EC for 2H-MoS2 and 1H- the intralayer and interlayer directions, respectively. The values
MoS2 are 15.316 and 15.156 eV per MoS2 unit, respectively. calculated from the first principles are in good agreement with
Their difference is exactly the interlayer interaction energy of those determined by Sun et al.17 from experimental data. We note
3D MoS2, which was calculated to be 160 meV. This indicates that contour plots of the calculated charge density, F(r) indicate
that the cohesion of the layers in 2H-MoS2 is the same as the that MoS2 layers have directional bonds. These bonds are formed
cohesion of 1H-MoS2. In addition to cohesive energy, the zero by Mo-d and S-p orbital hybridization and have a significant ionic
pressure bulk modulus B0, is an important mechanical property component. The charge transfer estimated by Mulliken74 analysis
of 3D crystals. Here we calculated the bulk modulus of indicates an excess electronic charge of 0.215 electrons on each S
16356 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

atom and depletion of electrons on each Mo atom amounting to

0.430 electrons. Negatively charged S atoms form negative charge
layers on both sides of 2D MoS2. We note the direction of
calculated charge transfer is in compliance with the Pauling’s
electronegativity scale,75 as well as with Born effective charges. On
the basis of the charge density analysis, MoS2 layers can be viewed
as a positively charged Mo atomic plane sandwiched between two
negatively charged S-atomic planes. As for 1H-MoS2, Born
effective charges are calculated to be ZB*[Mo] = 1.21 electrons
for Mo and ZB*[S] =
0.57 electrons for S. These values are close
to those of 2H-MoS2. However, high-frequency dielectric con-
stants of 1H-MoS2, which differ dramatically from 3D MoS2 due
to dimensionality effect, are ε = 4.58 and ε^ = 1.26 in the

)
intralayer and interlayer directions, respectively.
In Table 2 we present structure parameters and related
properties calculated for 2H-MoS2 and 1H-MoS2 using GGA
+D. These are lattice constants a, c, internal structural para-
meters, bulk modulus B0 (in-plane stiffness C for 1H-MoS2),
cohesive energy EC, Born effective charges, and high-frequency
dielectric constants.

’ LATTICE DYNAMICS
Previously, the lattice dynamics of 2H-MoS2 has been inves-
tigated both experimentally and theoretically. The phonon
dispersion curves and density of states have been obtained by
fitting the force constants to experimental data.6,24 Meanwhile
2D sheets of MoS2 including the single layer 1H-MoS2 were
synthesized.18
20,22 Our objective is (i) to calculate phonon
dispersion curves and total density of states of both 2H-MoS2
and 1H-MoS2 from the first principles including vdW correction;
(ii) to reveal the dimensionality effects between 3D and 2D
MoS2; (iii) to provide understanding of phonon anomalies
observed in Raman spectra.21,46
While Mo atoms in 2H-MoS2 occupy sites of D3h symmetry, S
atoms obey C3v symmetry. The overall symmetry of the crystal is
D6h, having 24 symmetry elements and 12 irreducible represen-
tations. Four second-order representations involve the lateral
(in-plane) displacements of Mo and S atoms, as indicated in
Figure 2. First-order representations are coupled with the
displacements perpendicular to the layers of atoms or parallel
to z axis. Similar observations are also valid for monolayer 1H-
MoS2 having D3h symmetry, 12 symmetry elements and 6
irreducible representations. In order for an irreducible represen-
tation to be infrared active mode, it must create a dipole moment
in the system. For 2H-MoS2, E1u and A2u modes are infrared
active. Similarly, E0 and A200 modes are infrared active modes of
1H-MoS2. Raman active modes induce polarization or quadruple
moment in the lattice. A1g, E1g and E2g modes are Raman active
modes for 2H-MoS2, so as A10 , E0 , and E00 modes for 1H-MoS2.
Verble and Wieting2,3 related the interlayer interaction in 2H-
MoS2 with the splitting of the frequencies of E2g and E1u modes.
In both modes the first layer atoms have similar displacements,
but the displacements of the second layer atoms are in opposite
direction as indicated in Figure 2. The minute difference between
Figure 2. (a) Calculated phonon dispersion curves of 2H-MoS2, Ω(k)
the frequencies of E2g and E1u shows that the interlayer vdW
versus k along symmetry directions of BZ, and corresponding density of states
(b). (c) and (d) are the same as (a) and (b) for 1H-MoS2. (e) Difference of interaction in 2H-MoS2 is small. Wakabayashi et al.6 first
the densities of states of 2H-MoS2 and 1H-MoS2 (see text). Phonon branches reported the phonon dispersion of 2H-MoS2 by neutron scatter-
derived from neutron scattering data6 and branches calculated by using a local ing. Because of experimental limitations, they only reported 12
basis set46,77 are indicated in (a) and (c) by green (light) squares, respectively. phonon branches along Γ
M and Γ
A directions out of 18
Infrared (IR) and Raman (R) active modes with symmetry representations available ones. The error term in their experiments is (%5.
and frequencies (cm
1) at the Γ-point are indicated. Bertrand23 reported that surface phonons have frequencies lower
16357 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

Table 3. Calculated Frequencies of Raman (R) and Infrared (IR) Active Modes (in cm
1) of 2H- and 1H-MoS2 at the Γ-Point and
Their Symmetry Analysisa

a
The subscripts u and g represent antisymmetric and symmetric vibrations, respectively. The other subscript i (i = 1, 2, 3) indicates the stretching modes.
IR and R frequencies of 2H- and 1H-MoS2 are calculated for the fully optimized lattice constants and internal structural parameters. Entries of IR and R
frequencies of 2H- and 1H-MoS2 indicated by (*) are calculated using the experimental lattice constants a and c of 2H-MoS2, but optimizing other
internal structural parameters. The entry with (**) is calculated with a = 3.14 Å and corresponds to a, which is smaller than the experimental lattice
constant a of 2H-MoS2.

than those of bulk MoS2 phonons. There is a softening of phonon interaction some branches are slightly split. The difference
modes upon going to the edges of nanocrystal 2H-MoS2. Recent between 2H- and 1H-MoS2 is substantiated by the difference
experimental study by Livneh and Sterer76 revealed the effect of of total density of states, i.e., ΔD(Ω) = D3D(Ω)
2D2D(Ω).
pressure and temperature on the Raman active modes of 2H- The plot of ΔD(Ω) in Figure 2e indicates an overall shift of
MoS2. They reported that upon increase of the temperature of critical point frequencies of 3D MoS2 to slightly higher values,
the system, the frequencies of Raman active modes decrease. while some modes show a reverse trend. The out of plane
Whereas the frequencies of Raman active modes increase with acoustical (ZA) mode of 1H-MoS2 has parabolic dispersion,
increasing pressure. since the transverse forces decay exponentially. Also the
Phonon dispersion curves of 3D and 2D MoS2 and their total LO
TO splitting is properly predicted. We also determined
densities of states calculated within DFPT55 using structure the infrared (IR) and Raman (R) active modes at the Γ-point of
optimized GGA+D are presented in Figure 2. Specific experi- BZ. Our results presented in Table 3.
mental data and earlier calculations are also indicated for the sake Earlier, Raman2,3,9,23,24,78 and infrared spectra2,3,17 of 2H-
of comparison. The phonon branches calculated from the first- MoS2 were studied experimentally. Wieting and Verble2,3 re-
principles for 3D MoS2 are in overall agreement with experi- ported three Raman active modes at 287, 383, and 409 cm
1. On
mental data as well as with that calculated using valence force the other hand, Chen and Wang78 have observed four Raman
field method.6 The calculated acoustical and optical branches of active modes in bulk at E22g = 32 cm
1, E1g = 286 cm
1, E12g =
1H-MoS2, Ω(k) are positive for any k in BZ. This indicates that 383 cm
1, and A1g = 408 cm
1. The small E22g mode is not
the suspended, single layer 1H-MoS2 structure is stable. observed by Wieting and Verble2,3 because of the spectral limit of
Phonon dispersion curves and corresponding density of states the Raman measurements between 20 and 1000 cm
1. Also
for 2H-MoS2 and 1H-MoS2 are similar, except that the number experimentally, at the zone center, IR active modes are observed
of branches in the former are doubled. Owing to the weak vdW at 384 cm
12,3,17 and 470 cm
12,3 (468 cm
117). For the case of
16358 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

1H-MoS2, Lee et al.21 investigated the Raman spectra of 2D as reported by Lee et al.,21 namely, that while A0 softens, E0
MoS2 sheets as a function of thickness down to a single layer 1H- becomes harder by going from 3D to 2D. To address the question,
MoS2. They reported that the frequency of the Raman active A1g whether the Raman active modes of the slabs of 2D MoS2 comply
mode of 2H-MoS2 (i.e., thick sheet) decreases gradually from with the above trends, we calculated the frequencies of two-layer and
408 to ∼403 cm
1 corresponding to the frequency of A10 of 1H- three-layer MoS2 using the experimental value of lattice constant, a of
MoS2 indicating the phonon softening. As for the Raman active 2H-MoS2. Since there are no data available for the spacings of layers
mode of 1H-MoS2 E0 displays a reverse behavior and hence its in slabs, we used again the experimental lattice constant c of 2H-
frequency decreases from 383.4 to 382 cm
1 corresponding to MoS2 and set the spacings equal to c/2. We found that A10 increases
the frequency of E2g mode of 2H-MoS2. Even if the lateral with increasing number of layers (namely, A10 f 404.9 cm
1 for
displacements of atoms in E2g mode are not affected significantly, bilayer MoS2 and A10 f 405.9 cm
1 for three layer MoS2, and
one nevertheless expects that all modes of 3D bulk MoS2 slightly approaches to A10 of 2H-MoS2. Nevertheless, the absence of
soften in single layer MoS2 due to the absence of interlayer experimental data on the structure of bilayer and three layer MoS2
interaction. The observation by Lee et al.21 was surprising. In fact, slabs prevents us from drawing more definite conclusions regarding
a recent study by Ramakrishna et al.46 ended up with a different phonon softening or phonon hardening with dimensionality.
trend; they observed that both A1g and E2g0 modes of 2H-MoS2 Finally, in Table 3 we present the frequencies of the IR active
sheets become softer as their thickness decreases. modes calculated using optimized as well as experimental lattice
In this work we attempted to clarify the controversial results constants. For 2H-MoS2, the present GGA+D calculations using
reported for the above crucial dimensionality effect. To this end, experimental lattice constants can give values in good agreement
we calculated the frequencies of Raman active modes of 2H-MoS2 with experimental data, namely, A2u = 467.6 cm
1 (as compared
and 1H-MoS2 and compared our results with available experi- to experimental value of 468 cm
1) and E1u = 381.6 cm
1 (as
mental data.2,3,17,21,46 Present GGA+D calculations predict compared to experimental value of 384 cm
1).2,17
Raman active modes E1g = 286.6 cm
1, E12g = 378.5 cm
1 and A1g =
400.2 cm
1 using the optimized lattice constants of a = 3.22 Å ’ DISCUSSION AND CONCLUSIONS
and c = 12.41 Å. The optimized lattice constant of 1H-MoS2 did In this paper we investigated lattice dynamics of 3D layered
not alter from the optimized value of 3D bulk despite the MoS2 and 2D single layer MoS2. We revealed that weak vdW
absence of interlayer interaction. Using the optimized lateral interaction between layers is critical to hold the layers together
lattice constant a = 3.22 Å, we obtained the frequencies of the and to calculate the interlayer spacing within 0.8% of the
Raman active modes as E0 = 380.2 cm
1 and A10 = 406.1 cm
1. experimental value. Therefore, the inclusion of the vdW interac-
(See Table 3.) According to these results of the ab initio GGA tion between the layers of MoS2 is found to be essential for ab
+D method, the frequencies of both modes should increase as initio calculations of energetics and optimized structures, cohe-
one goes from 3D to 2D single layer, which disagree with sion, and phonon dispersions. We examined two different
experimental results, except the behavior of E2g f E0 reported approaches to include vdW interaction in DFT calculations
by Lee et al.21 within GGA approximation. The spacings between MoS2 layers
The source of this disagreement between ab initio calculations or perpendicular lattice constant, which are critically overesti-
and experimental data is sought in the lattice constants, which are mated by GGA alone, are improved after the inclusion of vdW
overestimated by GGA+D calculations. We repeated the same interaction. As for the lateral lattice constant a, it is almost
GGA+D calculations using experimental lattice constants,10,63,64 independent of vdW interaction. Lattice dynamics and related
namely, a = 3.16 Å and c = 12.30 Å for 2H-MoS2 and a = 3.16 Å properties, such as phonon dispersions, effective charges, di-
for 1H-MoS2 by assuming that the lateral lattice constant, a, did electric constants, etc., calculated within GGA combined with
not change by going from 3D to single layer. We find that A1g = vdW are found to be in overall agreement with various experi-
407.7, E2g = 381.3 cm
1, and E1g = 277.8 cm
1 for 2H-MoS2, mental data and with the empirical values derived therefrom.
while for 1H-MoS2 A10 = 397.8 cm
1, E0 = 381.2 cm
1 and E00 = However, the shifts of the frequencies of Raman active modes by
274.5. Apparently, Raman active modes of 2H-MoS2 calculated going from 3D to 2D single layer are found to be very sensitive to
with experimental lattice constants are in good agrement with the values of lattice constants used in the calculations. For
observed Raman frequencies.2,78 Moreover, we are able to example, the shifts of Raman active modes predicted by using
reproduce the experimental trend that frequency of the Raman the lattice constants optimized through GGA with vdW correc-
active mode A1g softens for A1g f A0 , i.e., as the dimensionality is tion disagree with experimental data. We showed that by using
reduced from 3D to 2D. The change in the frequency is negligibly experimental lattice constants one can reproduce the experimen-
small for E2g f E0 . Noting the fact that the lattice constant of tally observed shifts. In a similar manner, one can obtain the
graphene, a, gets slightly smaller than that of 3D graphite, despite frequencies of infrared active modes, which agree better with
GGA+D optimizes a of 3D and 2D almost at the same value. experimental data, if the experimental lattice constants are used
Considering the possibility that the lattice constant of 1H-MoS2 in the calculations.
a can get smaller than the lateral lattice constant of 3D MoS2 a =
3.16 Å, we repeated our calculations for 1H-MoS2 using a = 3.14 Å ’ AUTHOR INFORMATION
and found that the frequency of E0 increases from 381.2 to
385 cm
1 confirming the anomalous effect reported by Lee Corresponding Author
et al. (See Table 3.) This is, however, a hypothetical situation *E-mail: [email protected].
and will be clarified when experimental data on the lattice
constant a of free-standing 1H-MoS2 will be available. We also
note that Raman active modes of 2H-MoS2 calculated by LDA, ’ ACKNOWLEDGMENT
which underestimates the lattice constants was able to repro- This work is supported by TUBITAK through Grant No.
duce the same dimensionality effect between 3D and 2D MoS2 104T537 and Grant No. 108T234. S.C. acknowledges TUBA for
16359 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361
The Journal of Physical Chemistry C ARTICLE

partial support. Part of the computational resources has been (33) Jaramillo, T. F.; Jorgensen, K. P.; Bonde, J.; Nielsen, J. H.;
provided by TUBITAK ULAKBIM, High Performance and Grid Horch, S.; Chorkendorff, I. Science 2007, 317, 100–102.
Computing Center (TR-Grid e-Infrastructure). We thank the (34) Lauritsen, J.; Helveg, S.; Laegsgaard, E.; Stensgaard, I.; Clausen,
DEISA Consortium (www.deisa.eu), funded through the EU B.; Topsoe, H.; Besenbacher, E. J. Catal. 2001, 197, 1–5.
FP7 project RI-222919, for support within the DEISA Extreme (35) Lauritsen, J.; Bollinger, M.; Laegsgaard, E.; Jacobsen, K.;
Norskov, J.; Clausen, B.; Topsoe, H.; Besenbacher, F. J. Catal. 2004,
Computing Initiative. 221, 510–522.
(36) Lauritsen, J. V.; Kibsgaard, J.; Olesen, G. H.; Moses, P. G.;
Hinnemann, B.; Helveg, S.; Norskov, J. K.; Clausen, B. S.; Topsoe, H.;
’ REFERENCES Laegsgaard, E.; Besenbacher, F. J. Catal. 2007, 249, 220–233.
(1) Title, R. S.; Shafer, M. W. Phys. Rev. B 1973, 8, 615–620. (37) Moses, P. G.; Hinnemann, B.; Topsoe, H.; Norskov, J. K.
(2) Wieting, T. J.; Verble, J. L. Phys. Rev. B 1971, 3, 4286–4292. J. Catal. 2007, 248, 188–203.
(3) Verble, J. L.; Wieting, T. J. Phys. Rev. Lett. 1970, 25, 362–365. (38) Raybaud, P.; Hafner, J.; Kresse, G.; Kasztelan, S.; Toulhoat, H.
(4) Kasowski, R. V. Phys. Rev. Lett. 1973, 30, 1175–1178. J. Catal. 2000, 189, 129–146.
(5) Mattheiss, L. F. Phys. Rev. Lett. 1973, 30, 784–787. (39) Sun, M.; Nelson, A.; Adjaye, J. J. Catal. 2004, 226, 32–40.
(6) Wakabayashi, N.; Smith, H. G.; Nicklow, R. M. Phys. Rev. B 1975, (40) Todorova, T.; Prins, R.; Weber, T. J. Catal. 2007, 246, 109–117.
12, 659–663. (41) Kline, G.; Kam, K.; Ziegler, R.; Parkinson, B. Sol. Energy Mater.
(7) Coehoorn, R.; Haas, C.; Dijkstra, J.; Flipse, C. J. F.; de Groot, 1982, 6, 337–350.
R. A.; Wold, A. Phys. Rev. B 1987, 35, 6195–6202. (42) Wilcoxon, J.; Thurston, T.; Martin, J. Nanostruct. Mater. 1999,
(8) Kobayashi, K.; Yamauchi, J. Phys. Rev. B 1995, 51, 17085–17095. 12, 993–997, 4th International Conference on Nanostructured Materi-
(9) Frey, G. L.; Tenne, R.; Matthews, M. J.; Dresselhaus, M. S.; als (NANO 98), Stockholm, Sweden, June 14
19, 1998.
Dresselhaus, G. Phys. Rev. B 1999, 60, 2883–2892. (43) Lauritsen, J. V.; Kibsgaard, J.; Helveg, S.; Topsoe, H.; Clausen,
(10) B€oker, T.; Severin, R.; M€uller, A.; Janowitz, C.; Manzke, R.; B. S.; Laegsgaard, E.; Besenbacher, F. Nat. Nanotechnol. 2007, 2, 53–58.
Voss, D.; Kr€uger, P.; Mazur, A.; Pollmann, J. Phys. Rev. B 2001, (44) Lee, C.; Li, Q.; Kalb, W.; Liu, X.-Z.; Berger, H.; Carpick, R. W.;
64, 235305. Hone, J. Science 2010, 328, 76–80.
(11) Rydberg, H.; Dion, M.; Jacobson, N.; Schr€oder, E.; Hyldgaard, (45) Radisavljevic, B.; Radenovic, A.; Brivio, J.; Giacometti, V.; Kis,
P.; Simak, S. I.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2003, A. Nat. Nanotechnol. 2011, 6, 147–150.
91, 126402. (46) Matte, H. S. S. R.; Gomathi, A.; Manna, A. K.; Late, D. J.; Datta,
(12) Fuhr, J.; Saul, A.; Sofo, J. Phys. Rev. Lett. 2004, 92, 026802. R.; Pati, S. K.; Rao, C. N. R. Angew. Chem., Int. Ed. 2010, 49, 4059–4062.
(13) Ivanovskaya, V. V.; Zobelli, A.; Gloter, A.; Brun, N.; Serin, V.; (47) Splendiani, A.; Sun, L.; Zhang, Y.; Li, T.; Kim, J.; Chim, C.-Y.;
Colliex, C. Phys. Rev. B 2008, 78, 134104. Galli, G.; Wang, F. Nano Lett. 2010, 10, 1271
1275, PMID: 20229981
(14) Huang, M.; Cho, K. J. Phys. Chem. C 2009, 113, 5238–5243. (48) Vanderbilt, D. Phys. Rev. B 1990, 41, 7892–7895.
(15) Moses, P. G.; Mortensen, J. J.; Lundqvist, B. I.; Norskov, J. K. (49) Perdew, J.; Chevary, J.; Vosko, S.; Jackson, K.; Pederson, M.;
J. Chem. Phys. 2009, 130, 104709. Singh, D.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671–6687.
(16) Botello-Mendez, A. R.; Lopez-Urias, F.; Terrones, M.; Terrones, (50) Grimme, S. J. Comput. Chem. 2006, 27, 1787–1799.
H. Nanotechnology 2009, 20, 325703. (51) Dion, M.; Rydberg, H.; Schr€oder, E.; Langreth, D. C.; Lundqvist,
(17) Sun, Q.-C.; Xu, X. S.; Vergara, L. I.; Rosentsveig, R.; Musfeldt, B. I. Phys. Rev. Lett. 2004, 92, 246401.
J. L. Phys. Rev. B 2009, 79, 205405. (52) Ceperley, D.; Alder, B. Phys. Rev. Lett. 1980, 45, 566–569.
(18) Joensen, P.; Frindt, R.; Morrison, S. Mater. Res. Bull. 1986, (53) Kresse, G.; Furthm€uller, J. Phys. Rev. B 1996, 54, 11169–11186.
21, 457–461. (54) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192.
(19) Novoselov, K.; Jiang, D.; Schedin, F.; Booth, T.; Khotkevich, V.; (55) Quantum-ESPRESSO is a community project for high-quality
Morozov, S.; Geim, A. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, quantum-simulation software, based on density-functional theory, and
10451–10453. coordinated by Paolo Giannozzi. See http://www.quantum-espresso.
(20) Mak, K. F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T. F. Phys. Rev. org and http://www.pwscf.org.
Lett. 2010, 105, 136805. (56) Arellano, J. S.; Molina, L. M.; Rubio, A.; Alonso, J. A. J. Chem.
(21) Lee, C.; Yan, H.; Brus, L. E.; Heinz, T. F.; Hone, J.; Ryu, S. ACS Phys. 2000, 112, 8114–8119.
Nano 2010, 4, 2695–2700. (57) Girifalco, L. A.; Hodak, M. Phys. Rev. B 2002, 65, 125404.
(22) Helveg, S.; Lauritsen, J. V.; Lægsgaard, E.; Stensgaard, I.; (58) Baskin, Y.; Meyer, L. Phys. Rev. 1955, 100, 544.
Nørskov, J. K.; Clausen, B. S.; Topsøe, H.; Besenbacher, F. Phys. Rev. (59) Trucano, P.; Chen, R. Nature (London) 1975, 258, 136.
Lett. 2000, 84, 951–954. (60) Zhao, Y. X.; Spain, I. L. Phys. Rev. B 1989, 40, 993–997.
(23) Bertrand, P. A. Phys. Rev. B 1991, 44, 5745–5749. (61) Bosak, A.; Krisch, M.; Mohr, M.; Maultzsch, J.; Thomsen, C.
(24) Jimenez Sandoval, S.; Yang, D.; Frindt, R. F.; Irwin, J. C. Phys. Phys. Rev. B 2007, 75, 153408.
Rev. B 1991, 44, 3955–3962. (62) Hanfland, M.; Beister, H.; Syassen, K. Phys. Rev. B 1989, 39,
(25) Seifert, G.; Terrones, H.; Terrones, M.; Jungnickel, G.; Frauen- 12598–12603.
heim, T. Phys. Rev. Lett. 2000, 85, 146–149. (63) Kam, K.; Parkinson, B. J. Phys. Chem. 1982, 86, 463–467.
(26) Li, Y.; Zhou, Z.; Zhang, S.; Chen, Z. J. Am. Chem. Soc. 2008, (64) Park, K.; RichardsBabb, M.; Freund, M.; Weiss, J.; Klier, K.
130, 16739–16744. J. Phys. Chem. 1996, 100, 10739–10745.
(27) Lebegue, S.; Eriksson, O. Phys. Rev. B 2009, 79, 115409. (65) Joensen, P.; Crozier, E.; Alberding, N.; Frindt, R. J. Phys. C:
(28) Bollinger, M.; Lauritsen, J.; Jacobsen, K.; Norskov, J.; Helveg, Solid State Phys. 1987, 20, 4043–4053.
S.; Besenbacher, F. Phys. Rev. Lett. 2001, 87, 196803. (66) Yang, D.; Sandoval, S.; Divigalpitiya, W.; Irwin, J.; Frindt, R.
(29) Ataca, C.; Sahin, H.; Akturk, E.; Ciraci, S. J. Phys. Chem. C 2011, Phys. Rev. B 1991, 43, 12053–12056.
115, 3934–3941. (67) Weiss, K.; Phillips, J. M. Phys. Rev. B 1976, 14, 5392–5395.
(30) Ataca, C.; Ciraci, S. J. Phys. Chem. C 2011, 115, 13303–13311. (68) Murnaghan, F. D. The Compressibility of Media under Ex-
(31) Rapoport, L.; Bilik, Y.; Feldman, Y.; Homyonfer, M.; Cohen, S.; treme Pressures. Proc. Natl. Acad. Sci. U.S.A. 1944, 30, 244.
Tenne, R. Nature 1997, 387, 791–793. (69) Batsanov, S. Effects of explosions on materials: modification
(32) Hinnemann, B.; Moses, P.; Bonde, J.; Jorgensen, K.; Nielsen, J.; and synthesis under high-pressure shock compression. High pressure
Horch, S.; Chorkendorff, I.; Norskov, J. J. Am. Chem. Soc. 2005, shock compression of condensed matter; Springer-Verlag: New York,
127, 5308–5309. 1994.

16360 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361

The Journal of Physical Chemistry C ARTICLE

(70) Yakobson, B. I.; Brabec, C. J.; Bernholc, J. Phys. Rev. Lett. 1996,
76, 2511–2514.
(71) Reddy, C.; Rajendran, S.; Liew, K. Nanotechnology 2006, 17,
864–870.
(72) Lee, C.; Wei, X.; Kysar, J. W.; Hone, J. Science 2008, 321,
385–388.
(73) It is known that Born effective charge is a tensor. The reported
experimental as well as calculated values are gvien for each atom in the xy
plane. The calculated values for the z direction are ZB*[Mo] = 0.58 and
ZB*[S] =
0.31 electrons for 2H-MoS2 and very small values for 1H-
MoS2.
(74) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1841–1846.
(75) It should be kept in mind that there are ambiguities in
calculating charge transfer. In fact, different methods result in different
values of charge transfer.
(76) Livneh, T.; Sterer, E. Phys. Rev. B 2010, 81, 195209.
(77) Soler, J.; Artacho, E.; Gale, J.; Garcia, A.; Junquera, J.; Ordejon,
P.; Sanchez-Portal, D. J. Phys.: Condens. Matter 2002, 14, 2745–2779.
(78) Chen, J.; Wang, C. Solid State Commun. 1974, 14, 857–860.

16361 dx.doi.org/10.1021/jp205116x |J. Phys. Chem. C 2011, 115, 16354–16361