Simulating 2D Materials with DFT: Graphene and Beyond

Two-dimensional materials — graphene, hexagonal boron nitride, transition-metal dichalcogenides, and a growing family of monolayers — sit at the frontier of condensed-matter physics and device engineering. Density functional theory (DFT) is the primary computational lens for understanding their electronic structure, but simulating a genuinely two-dimensional system inside a code built for three-dimensional periodicity takes care. This post covers the practical modeling choices and the physics DFT predicts for the most important 2D materials.

The Challenge of Modeling a 2D System in a 3D Code

Plane-wave DFT codes like Quantum ESPRESSO assume periodic boundary conditions in all three dimensions. A monolayer is periodic in the plane but isolated perpendicular to it. The standard trick is the slab model: place the monolayer in a unit cell that repeats normally in-plane but includes a large region of empty space — vacuum — along the out-of-plane direction.

The vacuum decouples each layer from its periodic images above and below, so the calculation approximates a single, freestanding sheet. Getting this right is the first and most common source of error in 2D simulations.

flowchart TB
  subgraph cell [Periodic cell along z]
    V1[Vacuum]
    L[2D monolayer]
    V2[Vacuum]
    V1 --- L
    L --- V2
  end
  L --> C[Converge vacuum 15-20 A]
  V1 -.->|periodic image| V2

Choosing the Vacuum Spacing

Too little vacuum and neighboring image layers interact spuriously through overlapping wavefunctions and, for polar or charged systems, long-range electrostatics. Too much vacuum wastes plane waves and inflates cost. The reliable procedure is a convergence test: increase the vacuum thickness and monitor the total energy, work function, or band structure until they stop changing.

As a rule of thumb, 15-20 Å of vacuum suffices for nonpolar monolayers like graphene and MoS2. Polar slabs, adsorbate-covered surfaces, or charged systems need more, often combined with a dipole correction to cancel the artificial field created by asymmetry across the slab.

Avoiding Periodic Image Interaction

Beyond vacuum thickness, two effects deserve attention.

• Electrostatic coupling. Any net dipole perpendicular to the layer generates a field that repeats with the cell. Quantum ESPRESSO’s dipole correction (dipfield and edir) removes this artifact, which matters greatly for Janus monolayers (e.g. MoSSe) and for adsorption studies where symmetry is broken.
• Coulomb truncation. For charged 2D systems and accurate GW-level calculations, a 2D Coulomb cutoff restricts the interaction to a single layer, eliminating the need for enormous vacuum. This is increasingly standard for defect and doping studies in monolayers.

For neutral, nonpolar monolayers, a well-converged vacuum is usually enough, and the workflow reduces to a careful slab relaxation followed by a band-structure calculation along the high-symmetry path.

Graphene: Dirac Cones and Semimetallic Behavior

Graphene — a single hexagonal sheet of carbon — is the canonical 2D material. Its signature electronic feature is the pair of Dirac cones at the K and K’ corners of the hexagonal Brillouin zone, where the valence and conduction bands meet at a single point. Near these points the bands are linear, so electrons behave as massless Dirac fermions with an effective “speed of light” of about \(10^6\,\mathrm{m/s}\).

DFT reproduces the Dirac cones cleanly at the GGA level. The practical requirements are:

• A fine k-point mesh and a band path that passes exactly through K, since the physics lives at that single point.
• Recognition that graphene is a zero-gap semimetal — there is no band gap, so it cannot be used as a transistor channel without opening a gap by other means (nanoribbons, substrate interaction, or a perpendicular field in bilayers).

Because standard DFT slightly renormalizes the Fermi velocity relative to experiment (many-body effects steepen the cones), quantitative transport work sometimes moves to GW. But for band topology and qualitative structure, GGA-DFT is sufficient and inexpensive.

Transition-Metal Dichalcogenides: The Direct-Gap Story

The TMD family, with formula MX2 (M = Mo, W; X = S, Se, Te), transformed 2D research because these materials are semiconductors — unlike graphene — with technologically useful gaps.

The Indirect-to-Direct Transition in MoS2

Bulk MoS2 is an indirect-gap semiconductor with a gap around 1.2-1.3 eV. Remarkably, when thinned to a single monolayer it becomes a direct-gap semiconductor near 1.8-1.9 eV, with the gap at the K point. This transition is why monolayer MoS2 photoluminesces strongly while the bulk does not — a landmark result that DFT explained by tracking how the conduction-band minimum and valence-band maximum shift with layer count.

Band-Gap Underestimation and Spin-Orbit Coupling

Two caveats matter for TMDs. First, semi-local DFT underestimates band gaps due to the well-known band-gap problem; quantitative gaps require hybrid functionals (HSE) or the GW approximation. Second, the heavy metal atoms make spin-orbit coupling (SOC) significant — it splits the valence band at K by hundreds of meV in WSe2, underpinning the valleytronics that make these materials interesting. SOC must be switched on (fully relativistic pseudopotentials) whenever spin-valley physics is the target.

Van der Waals Corrections for Stacking and Layering

Layered materials are held together out-of-plane by weak van der Waals (vdW) forces. Standard GGA functionals do not describe dispersion interactions, so they badly predict interlayer spacing and binding — often failing to bind layers at all. Any study of bilayers, multilayers, heterostructures, or adsorption on a 2D sheet must add a dispersion correction.

Common choices include:

• DFT-D2 / DFT-D3 (Grimme): semi-empirical pairwise corrections, cheap and widely used.
• Tkatchenko-Scheffler (TS): environment-dependent C6 coefficients.
• Nonlocal vdW-DF functionals (vdW-DF2, rVV10): dispersion built into the functional itself.

&SYSTEM
  ...
  vdw_corr = 'grimme-d3'   ! add D3 dispersion for interlayer binding
/

With a vdW correction, DFT recovers realistic interlayer distances (about 3.3-3.4 Å in graphite) and can distinguish stacking orders — AB versus AA in bilayer graphene, or the twist-angle physics behind magic-angle systems. See the Quantum ESPRESSO documentation for the supported dispersion schemes.

Strain Engineering in 2D Materials

Because they are atomically thin and mechanically flexible, 2D materials tolerate enormous elastic strain — graphene withstands over 20 percent before fracture — and strain becomes a knob to tune electronic properties. DFT is ideally suited to explore this “straintronics” design space.

The workflow is straightforward:

1. Take the relaxed monolayer cell.
2. Apply biaxial or uniaxial strain: scale lattice vectors by \((1 + \varepsilon)\).
3. Relax internal atomic positions at fixed cell shape.
4. Recompute the band structure.
5. Map gap, effective mass, or band-edge positions vs. strain.

Real predictions from this approach include: tensile strain reduces the MoS2 gap and can drive a direct-to-indirect transition; strain gradients create pseudo-magnetic fields in graphene; and strain shifts band edges enough to tune band alignment in van der Waals heterostructures for photocatalysis and photovoltaics.

A Reproducible 2D DFT Checklist

To summarize the choices that separate a trustworthy 2D calculation from a misleading one:

• Vacuum: converge it explicitly; 15-20 Å for nonpolar monolayers.
• Dipole/Coulomb corrections: add for polar, asymmetric, or charged slabs.
• k-points: dense in-plane mesh, single point out-of-plane; hit K for graphene/TMDs.
• vdW correction: mandatory for any multilayer or adsorption study.
• SOC and functional: use fully relativistic pseudopotentials plus hybrids/GW when quantitative gaps or valley splitting matter.

These practices apply equally to hBN, phosphorene, MXenes, and the ever-growing catalog of predicted monolayers. For more on the underlying methods, see our computational engines page and the rest of the blog.

Beyond Single Layers: Heterostructures and Moire Physics

The frontier of 2D modeling is van der Waals heterostructures — vertically stacked, chemically distinct monolayers held together by dispersion forces. Because the layers need not be lattice-matched, stacking two sheets with a small lattice mismatch or a relative twist produces a long-wavelength moire superlattice. The magic-angle twisted bilayer graphene that hosts correlated insulating and superconducting phases is the most famous example.

Simulating these systems pushes DFT to its practical limits: a moire cell can contain thousands of atoms, and the flat bands that drive the interesting physics are exquisitely sensitive to interlayer distance, so an accurate vdW correction is non-negotiable. Type-II band alignment in TMD heterobilayers (MoS2/WSe2, for example) spatially separates electrons and holes into different layers, producing long-lived interlayer excitons that DFT band-alignment calculations predict and that underpin proposals for photodetectors and excitonic devices. Getting the band offsets right again demands hybrid functionals plus SOC, applied to cells large enough to accommodate the moire periodicity.

Run it on Simatra

Converging vacuum spacing, running hybrid-functional band structures, mapping strain over many cell distortions, and building large van der Waals heterostructure supercells are exactly the kind of repetitive, memory-hungry jobs where hardware acceleration pays off. Simatra runs these DFT workflows on GPU-accelerated clusters powered by the GPU-Opt-V2 instance, achieving up to 5x faster convergence and handling supercells up to roughly 2,000 atoms — enough for moire and heterostructure cells. Choose Quantum ESPRESSO or our native KRONOS engine and scale from a single monolayer to a full strain sweep. Get started with a free trial and $100 in credits at app.simatra.io.