Shell model extension to the valence force field: application to singlelayer black phosphorus
Received
20th September 2018
, Accepted 21st November 2018
First published on 21st November 2018
We propose an extension of the traditional valence force field model to allow for the effect of electronic polarization to be included in the interatomic potential. Using density functional theory as a reference, this model is parameterized for the specific case of singlelayer black phosphorus by fitting the phonon dispersion relation over the entire Brillouin zone. The model is designed to account for the effect of induced dipole interaction on the longwavelength ( → 0) modes for the case of homopolar covalent crystals. We demonstrate that the near Γpoint frequencies of the IRactive modes are substantially damped by the inclusion of the induced dipole interaction, in agreement with experiment. The fitting procedure outlined here allows for this model to be adapted to other materials, including but not limited to twodimensional crystals.
1 Introduction
Despite the broad availability and popularity of more descriptive and predictive modern theoretical and computational methods such as density functional theory (DFT), force fields for molecular dynamics remain of considerable relevance as they can be used to perform largescale calculations and simulations that would be too computationally expensive (in terms of processing time and memory requirement) for DFT Born–Oppenheimer dynamics.^{1} Force fields have proven remarkably useful and have therefore become nearly ubiquitous in computational condensed matter physics,^{2–4} computational chemistry,^{5} and computational biology.^{6} The valence force field (VFF) is one example and has been used to study the thermomechanical properties of graphene and the electronic properties of GaAs,^{7,8} for instance. In addition, the importance of dynamical properties of lowdimensional materials is ubiquitous in many areas of research and applications relying on light–matter interactions. In this context, models capable of faithfully accounting for optical phonon modes are projected to play a key role in advancing nanoscience. However, because this potential only includes information about interactions related to the formation of atomic bonds, it fails to account for the effects of electronic polarization.^{9} These effects become important in the study of the lattice vibrations of ionic and homopolar covalent crystals (e.g., black phosphorus, wurtzite, diamond, scheelite, etc.), where the induced dipole interactions arising from the lattice vibrations themselves have an important effect on the longwavelength ( → 0) vibrational modes.^{10} In this paper, we develop an extension to the VFF model to allow for the effect of electronic polarization to be accounted for. This extension is based off of the shell model proposed by Dick and Overhauser.^{11} It is designed to capture both the success of the VFF model at reproducing the mechanical properties of materials^{12,13} and the ability of the shell model at describing the effects of electronic polarization,^{10} and to combine them into one cohesive model. Here, we introduce the theoretical model and employ the nonlinear Nelder–Mead method^{14} to determine a good set of numerical parameters by fitting the prediction of the model with those made using DFT. We focus on singlelayer black phosphorus (SLBP), as it is a prototypical example of a homopolar crystal for which VFF fails at reproducing some specific vibrations.^{9} In particular, we demonstrate the resulting model's ability to improve the description of the longwavelength phonon behavior in comparison to the pure VFF case. These improvements are indicative of the model's successful ability to describe both the atomic interactions mediated by covalent bonds and the interactions that arise from the induced dipoles that occur in vibrating crystals.
2 Model
The shell model, which is analogous to the Drude oscillator, is a phenomenological theory for describing polarization effects in molecular systems. It does so by placing charged massless shells around each atomic core in the system; these shells are then assumed to interact with each other and with the atomic cores as if they were attached via springs.^{10,11} While the notion that atoms are surrounded by virtual charged shells interacting via springs may seem a crude approximation, the approach has proven to be an effective method of heuristically describing the effects of electronic polarization while keeping the number of free parameters reasonably small.^{10} Note that since we intend merely to adapt the shell model for use alongside the more common VFF model, we will not provide details for the atomic–atomic interactions between the cores present in the original formulation of the shell model and assume that we can write our new total potential function in the form: 
Φ = ϕ^{(VFF)} + ϕ^{(shell)},  (1) 
here ϕ^{(VFF)} corresponds to the usual valence force field potential and ϕ^{(shell)} corresponds to the part of the potential arising from the shells and their interactions. Before moving on to our discussion of the shell model we will provide a quick review of valence force models and, more specifically, the SLBPspecific VFF that we use here.^{12}
VFFs are a general class of simple empirical models for describing the interatomic forces of a given structure or molecule. They do so by describing the atomic potential in terms of simple quadratic functions of bond length and bond angle. From these potentials, a complete description of the molecular forces can be obtained and used to carry out molecular dynamic simulations. An example of such a VFF model is provided here for the specific case of SLBP:^{12}

 (2) 
Here the sum includes every atom in the structure, all of the same pucker nearest neighbors, the single different pucker nearest neighbor and both same pucker nearest neighbors to atom
i, respectively. The equilibrium bond length is given as
d. δ
r_{ij} is the relative change in bond length between atoms
i and
j. δ
θ_{jik} is the change in angle between atoms
j,
i and
k with
i at the apex. Lastly, the
r subscript on the force constant
K denotes bond stretching terms, the
θ subscript denotes bond bending terms, and the
rr and
rθ subscripts refer to bond–bond and bond–angle interactions, respectively. The primed terms indicate that the bond or angle in question involves atoms in different puckers (as opposed to unprimed terms which only involve atoms in the same pucker). Values for these force constants are given in
Table 1. A more detailed discussion of this potential, or of VFF models in general, can be found in the literature,
^{12,13,15} and we will now focus our attention on
ϕ^{(shell)} for the remainder of the discussion.
Table 1 Table of original VFF force constants (eV Å^{−2}) for singlelayer black phosphorus given by Midtvedt and Croy^{12} and force constants for VFF found via fitting. The r subscript denotes bond stretching terms, the θ terms denote bond bending terms, and the rr and rθ terms refer to bond–bond and bond–angle interactions, respectively. The primed terms indicate that the bond or angle in question involves atoms in different puckers (as opposed to unprimed terms which only involve atoms in the same pucker)

K
_{
r
}

K
_{
r
}′ 
K
_{
θ
}

K
_{
θ
}′ 
K
_{
rr′}

K
_{
rr′}′ 
K
_{
rθ
}

K
_{
rθ
}′ 
K
_{
rθ
}′′ 
Ref. 12

11.17 
10.3064 
1.18 
0.9259 
−0.6763 
1.2449 
0.58 
1.932 
0.797 
This work 
10.60 
10.45 
0.86 
0.63 
−0.77 
1.36 
0.61 
2.02 
0.94 
Because all of the shell–shell and the shell–core terms are presumed to vary according to a harmonic potential, it is easy to write down ϕ^{(shell)} in a general form as follows:

 (3) 
Here the interatomic bond vectors are denoted by and the displacement vectors of the shells from their equilibrium position around their cores are denoted by . The nindex represents groups of nth nearest neighbors. The force constant of the springs are represented by the ɀ coefficients, where the “s” and “a” subscripts denote that they belong to shell–shell terms and shell–other atom core terms, respectively. The iindex subscript denotes that the particular atom belongs to a shell–selfatom core term. The sum over i covers all atoms in the supercell while the sum over j∈_{n}i runs over all elements j in the set of nth nearest neighbors of i. Lastly, the sum over n runs over groups of nth nearest neighbors of atom i.
From this expression for the potential due to the shell terms, we obtain the force on an atom i from the shell around atom j by:

 (4) 
Similarly, we can take the force on the shell around atom i due to the shell around atom j to be:

 (5) 
The atom–atom forces are all taken care of by the valence force field part of the potential and we do not need to elaborate on them here, as they have been described in detail elsewhere.^{12}
Now that we have formulated the potential and specified how to find the forces, it is possible to use this new model to perform molecular dynamics simulations provided that one has a good set of VFF and shell force constants. While VFF force constants are widely available in the literature for a large range of materials,^{12,13} the shell model parameters will need to be found via fitting. In the rest of this work, we will focus on SLBP and use results from DFT as reference for the fitting procedure. The DFT phonon dispersion for SLBP used as a reference to parameterize the VFF is in qualitative agreement with the partial phonon dispersion of bulk BP obtained from inelastic neutron scattering experiments.^{16,17} In any case, rigorous parameterizations of force fields using only firstprinciples data constitutes an active area of research.^{18–22}
For the sake of comparison, before we begin fitting the VFF + shell model to the DFT data, we optimize the VFF model by itself to ensure that we find a set of parameters for SLBP that is similar to the literature values.^{12} The resulting set of parameters are given, along with the literature values, in Table 1. We can see that this optimization does indeed yield a similar set of parameters to those reported by Midtvedt and Croy. Small differences can be explained by the use of a different reference for the fitting procedure and the fact that fitting is performed over the entire Brillouin zone (BZ).
3 Singlelayer black phosphorous and the failure of VFF
BP consists of weakly interacting ABstacked puckered layers, each with inequivalent armchair and zigzag directions. It has an orthorhombic lattice and, in the fewlayer configuration, it belongs to the Pmna space group for an odd number of layers and to Pbcm for an even number of layers. In both cases, as well as in bulk, BP belongs to the D_{2h} point group. With knowledge of the space group symmetry it is possible to assign symmetry to each vibrational mode, as shown in Fig. 1. BP's electronic band gap is direct irrespective of the number of layers and increases with decreasing number of layers due to quantum confinement and a change in screening.^{23} In particular, the band gap of the pristine material is predicted to vary from 0.4 eV in the bulk up to 2.0 eV in the singlelayer,^{24} and can be strongly modified by encapsulation due to a strong screening effect.^{25} Furthermore, the ambipolar mobility along the armchair direction can reach a value^{26} of 55000 cm^{2} V^{−1} s^{−1}, and is only 3.6 times smaller than the largest value ever measured (graphene).^{27} Understanding and optimizing such a desirable property requires a description of the crystal vibrations as scattering by phonons is known to contribute significantly to the reduction of carrier mobility. Accurate calculations of the phonon band structure and related properties from firstprinciples DFT can be time consuming, whereas ultrafast simulations can be performed with classical potentials. As mentioned earlier, a VFF potential has already been optimized for SLBP; however, it is based on experimental Raman and infrared frequencies (which only includes the Γpoint) and, because it neglects the description of polarization, the infrared modes are not described very well.^{12} Therefore we suggest to optimize the potential using a completely firstprinciples DFT phonon band structure including all highsymmetry directions, and to add nonlocal terms to the potential which include electronic polarization effects with the goal of improving the description of the infrared active modes, thus illustrating the usefulness of the model presented.

 Fig. 1 Irreducible representations of acoustic and optical modes for singlelayer black phosphorus, at the Γ point. Atomic displacements are indicated by red arrows. The circle and cross indicate displacements going out of the plane of the page and going into it, respectively.  
4 Computational details
Before moving to the fitting procedure, we first provide technical details on how the DFT calculations were performed. In order to obtain the firstprinciples vibrational data for freestanding SLBP, we used the Vienna Ab initio Simulation Package (VASP).^{28–31} The ion cores were modeled using the projector augmented wave (PAW) pseudopotentials.^{32} For these calculations a planewave energy cutoff of 500 eV and a Gaussian smearing of 0.005 eV were used to yield well converged total energies and forces. In order to ensure that there were no unwanted interactions with the periodic images of our structure in the ẑdirection, we added a 12 Å vacuum region. To include van der Waals (vdW) corrections in our calculations, we used the optB86bvdW scheme.^{33,34} In fact, this functional improves the match with experiment as SLBP has a puckered structure, which suggests that noncovalent interactions take place between nonbonded atoms. For the sake of comparison, we have compared local, semilocal and nonlocal exchange–correlation (xc) functionals. The particular choice of (xc) functional is based of the study of thermodynamic and vibrational properties of freestanding and supported (on Au(111)) monolayer BP.^{35} In contrast to predictions made using other (xc) functionals, calculations performed with the optB86bvdW scheme are in agreement with existing experimental evidence.^{36} During relaxation, all of the atoms were relaxed to a force cutoff of 10^{−4} eV Å^{−1} and the kpoint sampling that we used was based off of a Γpoint centered grid. The primitive cells were optimized on a (10 × 15 × 1) kpoint grid.
The Phonopy code^{37} was used for the harmonic phonon calculations. For the DFT phonon dispersion relation we chose an (8 × 8 × 1) supercell to ensure that the forces, which were calculated using VASP, were well converged. For the VFF and VFF + shell calculations we used a (7 × 7 × 1) supercell. Prior to the VFF and VFF + shell force calculation, the atomic positions of the initial structure were relaxed with the lattice vectors (a = 4.51 Å and b = 3.30 Å)^{35} held constant to ensure that the kpath remained the same length. This is important as the residual cannot be properly calculated during the parameter optimization if the DFT and VFF or VFF + shell frequency data have a different kspace domain.
5 Fitting procedure
In this section we set out to demonstrate that the proposed shell extension to the VFF model resolves issues noted for SLBP. Initially, we randomly select a starting set of parameters, β, for the first through eighth nearest neighbor groups. As there are two terms per nth nearest neighbor group (a shell–shell term and a shell–core term) plus an additional term for the atom to “selfshell” interaction, the procedure requires, in principle, a total of 17 parameters to describe these eight terms. Terms beyond the eighth nearest neighbor group were left out as their contributions are increasingly negligible and the results are not significantly changed by their inclusion. Once a random set of parameters is chosen, the phonon dispersion is calculated using the proposed model. After this is done, a residual function is defined as follows: 
 (6) 
where ω_{i}() is the frequency of the ith normal mode at the phonon wave vector . It is important to note that the sum over is actually an integral (crystal momentum being a continuous variable spanning the entire first BZ), however for our purposes here it is sufficient to use a discretized sum over the path. G_{ij}() is a weighting function defined to be: 
 (7) 
where _{i}() represents the normalized eigenvector of the ith normal mode at phonon wave vector . This optimization procedure was carried out several times, with different random initial parameters, to ensure that the results were indeed consistent and not highly dependent on the randomly selected initial parameter set.
We will now explain why we introduce the weighting function in the form provided above. Care must be taken to ensure that only bands belonging to the corresponding normal mode are compared to one another. While it is straightforward to do this for bands at highsymmetry points (provided that there is no accidental degeneracy), it becomes considerably more complicated as one moves along the path from one symmetry point to another. This is because the accepted method for labeling bands is to order them by frequency,^{38} and this can introduce singular critical points that are not required by symmetry. A singular critical point is a point along the band where one or more component of changes sign discontinuously. These unnecessary critical points are called accidental degeneracies (as opposed to essential degeneracies, which are required by symmetry) and are not manifestations of underlying physics, but are instead a result of the method of band labeling. As a result, if care is not taken to correct the appearance of accidental degeneracies then bands belonging to different normal modes may get compared during optimization. In addition, bands can be initially ordered incorrectly, especially when the force field parameters are far from their optimal values. These two effects often lead to a numerical frustration of the minimization that prevents the optimization algorithm from converging to a physically significant minimum, as the residual will have been calculated by comparing bands that do not belong to the same normal mode. The weighting term is included in the residual function to alleviate this situation, as the eigenvectors of two distinct normal modes (at a given point and computed by the two different methods) are orthonormal to each other when the procedure is close to convergence. Numerically, the introduction of the weighting function makes it convenient to just take the difference between all bands and then multiply them by the magnitude of the inner product of the eigenvectors to filter out the incorrect pairs.
The residual function S can then be minimized with respect to β in order to obtain a complete set of good force parameters for a given material. We used the nonlinear Nelder–Mead algorithm to carry out this minimization as it does not require derivatives of the residual function to be taken with respect to the parameters, which in our case would be difficult.^{14}
In the specific case of SLBP, we found that the fourth, fifth, and eighth nearest neighbor terms introduced nonphysical instabilities into the structure that were manifested as imaginary frequencies in the phonon dispersions, and they were therefore neglected. This leaves us with only the first, second, third, sixth and seventh nearest neighbor terms in the shell portion of the potential (11 shell parameters). The atomic relations that these terms refer to are depicted in Fig. 2. For the initial optimization of the shell force constants we used the VFF and corresponding force constants outlined by Midtvedt and Croy for SLBP.^{12} The initial force constants are given in Table 1.

 Fig. 2 Neighbors included in the shell extension to the VFF model for singlelayer black phosphorus. Gray and purple denote atoms in the top and bottom plane (puckers) and red denotes the reference atom. (1) to (7) represent the first to seventh nearest neighbors, respectively.  
Once the initial shell parameter optimization is complete, with the VFF parameters held constant, the entire potential is reoptimized with all 20 parameters (9 VFF + 11 shell). The final set of force constants is given in Table 2.
Table 2 Table of optimized parameters for VFF and shell model (eV Å^{−2})
K
_{
r
}

K
_{
r
}′ 
K
_{
θ
}

K
_{
θ
}′ 
K
_{
rr′}

K
_{
rr′}′ 
K
_{
rθ
}

K
_{
rθ
}
^{′}

K
_{
rθ
}′′ 
10.4079 
8.8248 
0.8489 
0.8113 
0.0322 
0.9197 
0.2593 
1.3175 
0.4056 
ɀ
_{i}

ɀ
_{a}

ɀ
_{s}

ɀ
_{a}′ 
ɀ
_{s}′ 
ɀ
^{(2)}_{a}

ɀ
^{(2)}_{s}

ɀ
^{(6)}_{a}

ɀ
^{(6)}_{s}

ɀ
^{(5)}_{a}

ɀ
^{(5)}_{s}

−1.108 
−0.852 
0.1957 
−10.2613 
−0.014 
0.0624 
11.6562 
0.1005 
0.1075 
0.0584 
0.1451 
Once this final set of optimized parameters is found, we can use it with the proposed model to calculate the phonon dispersion relation. A plot of the resulting spectrum along with the DFT result is shown in Fig. 3.

 Fig. 3 Phonon dispersion of singlelayer black phosphorus calculated with (a) VFF + shell model (b) DFT. Comparison between the phonon dispersion calculated with DFT and with (c) VFF + shell model and (d) VFF only.  
In order to get a better understanding of how the inclusion of the electronic polarization via the shell model extension affects the normal modes of lattice vibrations for SLBP, we have superimposed the phonon dispersion plots obtained with our VFF + shell model over the DFT result and displayed it alongside a similar comparison between the plain VFF and DFT dispersions. These juxtaposed plots are displayed in Fig. 3 and the Γpoint frequencies of each normal mode, as calculated with each of these three different models, are given in Table 3.
Table 3
ΓPoint frequency (cm^{−1}) of each mode as calculated using DFT, VFF^{12} and VFF + shell, respectively
Mode 
DFT 
VFF 
VFF + shell 
B_{1u} 
126.1 
246.1 
162.6 
B_{1g} 
188.0 
194.7 
209.1 
B^{1}_{3g} 
220.4 
233.2 
230.1 
A^{1}_{g} 
349.1 
352.0 
348.9 
A_{u} 
413.7 
439.5 
407.5 
B_{2g} 
419.7 
437.0 
410.5 
B^{2}_{3g} 
420.6 
465.1 
412.0 
A^{2}_{g} 
445.1 
473.3 
443.2 
B_{2u} 
460.0 
485.9 
452.1 
From these comparisons it is clear that, aside from a general improvement of fit, the apparent “twisting” of the B_{1u} and B^{1}_{3g} modes in the VFFonly result was resolved and the overall behavior of the highfrequency modes are greatly improved by the shell model extension. Note that the B_{1u} mode, which suffered a longwavelength ( → 0) damping of ∼83 cm^{−1} from the addition of the shell extension to the VFF model, is an IRactive mode.^{9,39} This observation is important as it is understood^{9,10,12} that IRactive modes are highly sensitive to the effects of electronic polarization in the near Γpoint region of kspace and their behavior cannot be accurately described without taking these interactions into account.
6 Conclusion
We proposed an extension to the traditional VFF model, based on the historical shell model, to allow for the effects of electronic polarization to be accounted for in the study of lattice vibrations using classical methods. We showed how the developed model can be fitted to reproduce information for a specific material (in our case we chose to fit it to SLBP). This allowed us to calculate the phonon dispersion relation using the new model and compare it and the VFFonly result to the phonon dispersion found using DFT. We highlighted the importance of the effects of electronic polarization in accurately describing the longwavelength behavior of the IRactive modes, as well as the highfrequency modes of homopolar covalent crystals. This enhancement is important as it allows for an improved use of classical methods in the study of the lattice vibrations of lowdimensional covalent crystals.
Conflicts of interest
There are no conflicts to declare.
Acknowledgements
Part of this work was performed using supercomputing resources provided by the Center for Computational Innovations (CCI) at Rensselaer Polytechnic Institute. DT acknowledges support from the Office of Naval Research. AC was supported by NSF Grant EFRI 2DARE (EFRI1542707).
References
 T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. EngelHerbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama and A. C. T. van Duin, npj Comput. Mater., 2016, 2, 15011 CrossRef CAS.
 D. W. Brenner and B. J. Garrison, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 34, 1304 CrossRef CAS.
 C. Pryor, J. Kim, L. W. Wang, A. J. Williamson and A. Zunger, J. Appl. Phys., 1998, 83, 2548–2554 CrossRef CAS.

S. Bose, W. J. Fan and D. H. Zhang, presented in part at the 2017 Conference on Lasers and ElectroOptics Pacific Rim (CLEOPR), Singapore, 2017 Search PubMed.
 S. C. Kohale, S. Pratihar and W. L. Hase, J. Phys. Chem. Lett., 2018, 9, 1554–1560 CrossRef CAS PubMed.
 E. Boulanger, L. Huang, C. Rupakheti, A. D. MacKerell Jr and B. Roux, J. Chem. Theory Comput., 2018, 14, 3121–3131 CrossRef CAS PubMed.
 A. Lajevardipour, M. NeekAmal and F. Peeters, J. Phys.: Condens. Matter, 2012, 24, 175303 CrossRef CAS PubMed.
 C. Delerue and M. Lannoo, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 3966 CrossRef CAS.
 C. Kaneta, H. KatayamaYoshida and A. Morita, J. Phys. Soc. Jpn., 1986, 55, 1213–1223 CrossRef CAS.

G. Venkataraman, L. Feldkamp and V. Sahni, Dynamics of perfect crystals, MIT Press, Cambridge, 1975 Search PubMed.
 B. G. Dick Jr and A. W. Overhauser, Phys. Rev., 1958, 112, 90 CrossRef.
 D. Midtvedt and A. Croy, Phys. Chem. Chem. Phys., 2016, 18, 23312–23319 RSC.
 V. Perebeinos and J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 241409 CrossRef.
 J. A. Nelder and R. Mead, Comput. J., 1965, 7, 308–313 CrossRef.

J. Wan, Y.W. Tan, J.W. Jiang, T. Chang and X. Guo, arXiv:1808.01714v1, 2018.
 Y. Fujii, Y. Akahama, S. Endo, S. Narita, Y. Yamada and G. Shirane, Solid State Commun., 1982, 44, 579–582 CrossRef CAS.
 Y. Yamada, Y. Fujii, Y. Akahama, S. Endo, S. Narita, J. Axe and D. McWhan, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 30, 2410 CrossRef CAS.
 F. Ercolessi and J. B. Adams, EPL, 1994, 26, 583 CrossRef CAS.
 S. Izvekov, M. Parrinello, C. J. Burnham and G. A. Voth, J. Chem. Phys., 2004, 120, 10896–10913 CrossRef CAS PubMed.
 P. Brommer and F. Gähler, Modell. Simul. Mater. Sci. Eng., 2007, 15, 295 CrossRef CAS.
 H. Sheng, M. Kramer, A. Cadien, T. Fujita and M. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 134118 CrossRef.
 L. Vanduyfhuys, S. Vandenbrande, J. Wieme, M. Waroquier, T. Verstraelen and V. Van Speybroeck, J. Comput. Chem., 2018, 39, 999–1011 CrossRef CAS PubMed.
 H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomanek and P. D. Ye, ACS Nano, 2014, 8, 4033–4041 CrossRef CAS PubMed.
 L. Liang, J. Wang, W. Lin, B. G. Sumpter, V. Meunier and M. Pan, Nano Lett., 2014, 14, 6400–6406 CrossRef CAS PubMed.
 D. Y. Qiu, F. H. da Jornada and S. G. Louie, Nano Lett., 2017, 17, 4706–4712 CrossRef CAS PubMed.
 J. Yang, S. Tran, J. Wu, S. Che, P. Stepanov, T. Taniguchi, K. Watanabe, H. Baek, D. Smirnov, R. Chen and C. N. Lau, Nano Lett., 2017, 18, 229–234 CrossRef PubMed.
 K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim and H. Stormer, Solid State Commun., 2008, 146, 351–355 CrossRef CAS.
 G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561 CrossRef CAS.
 G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269 CrossRef CAS.
 G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186 CrossRef CAS.
 G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50 CrossRef CAS.
 P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979 CrossRef.
 J. Klimes, D. R. Bowler and A. Michaelides, J. Phys.: Condens. Matter, 2009, 22, 022201 CrossRef PubMed.
 J. Klimes, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195131 CrossRef.
 D. Tristant, A. Cupo and V. Meunier, 2D Mater., 2018, 5, 035044 CrossRef.
 J. L. Zhang, S. Zhao, C. Han, Z. Wang, S. Zhong, S. Sun, R. Guo, X. Zhou, C. D. Gu and K. D. Yuan,
et al.
, Nano Lett., 2016, 16, 4903–4908 CrossRef CAS PubMed.
 A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5 CrossRef CAS.
 J. C. Phillips, Phys. Rev., 1956, 104, 1263–1277 CrossRef CAS.
 Y. Aierken, D. Cakir, C. Sevik and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 081408 CrossRef.

This journal is © the Owner Societies 2019 