Graphene: Building a Structure from Scratch

Graphene is a single layer of carbon atoms arranged in a honeycomb, with a C–C bond length of $l = 1.42$ Å. LAMMPS has no built-in “graphene” lattice, so we describe it ourselves. That is the most useful skill in this whole guide, because the same recipe builds almost any crystal.

The idea: a unit cell that tiles the plane

Any crystal can be described by a unit cell (a box defined by lattice vectors $\vec a_1, \vec a_2, \vec a_3$) and the basis atoms inside it. Copy the cell over and over and you get the whole crystal.

For graphene, I use an orthogonal (rectangular) cell. It’s not the smallest possible cell, but its sides line up with x and y, so it fits neatly into a rectangular simulation box.

In a textbook you would usually write basis positions in Cartesian coordinates. LAMMPS instead wants fractional coordinates: each position written as a fraction of the lattice vectors.

To convert, divide each Cartesian coordinate by the length of the matching lattice vector. For atom 2: $\tfrac{\sqrt3/2}{\sqrt3} = \tfrac12$ and $\tfrac{1/2}{3} = \tfrac16$.

Step 1: size the box

variable l equal 1.42 # bond length in angstroms
variable m  equal 20 # number of unit cells along length
variable n  equal 20 # number of unit cells along width
variable L equal $(sqrt(3)*v_l*v_m) # length of the graphene sheet
variable W equal $(3*v_l*v_n) # width of the graphene sheet

region box block 0 ${L} 0 ${W} -10 10 units box
create_box 1 box
mass 1 12.01 # mass of carbon atom

The box length is the number of cells times the cell length: $L = m\sqrt3\,l$ and $W = 3n\,l$. This matters more than it looks. Because the box is periodic, the right edge of the sheet must connect perfectly to the left edge. If the box doesn’t hold a whole number of unit cells, atoms at the edges end up too close together or too far apart, and the sheet has a defect line running along the boundary.

region box block xlo xhi ylo yhi zlo zhi defines a rectangular block, and create_box 1 box turns it into the simulation box with 1 atom type. units box means the numbers are in Å. Without it, LAMMPS would multiply them by the lattice spacing.

The z range (−10 to 10) doesn’t matter here, because this is a 2D simulation (dimension 2).

Step 2: the three building commands

lattice custom ${l} &
  a1 $(sqrt(3)) 0.0 0.0 &
  a2 0.0 3.0 0.0 &
  basis 0.0 0.0 0.0 &
  basis $(1/2) $(1/6) 0.0 &
  basis $(1/2) $(1/2) 0.0 &
  basis 0.0 $(2/3) 0.0

region mybox block 0 ${L} 0 ${W} -0.1 0.1 units box
create_atoms 1 region mybox

lattice style scale keyword values... defines the unit cell.

• custom because graphene isn’t one of the built-in styles (sc, bcc, fcc, hcp, diamond…).
• scale = 1.42: every vector below is multiplied by this, so a1 $(sqrt(3)) 0 0 becomes $(\sqrt3\,l, 0, 0)$.
• a1, a2 (and a3) are the lattice vectors; each basis line adds one atom in fractional coordinates.
• A long command can be split over several lines by ending each line with &.

region ID style args marks out the space to fill. Regions can be blocks, cylinders, spheres, cones and more. Here, a thin slab around z = 0.

create_atoms type style args fills that region with the lattice, using atom type 1 (carbon) and style region.

The order matters: lattice first (what to tile), then region (where), then create_atoms (do it).

Step 3: everything else

The rest of the script is the standard skeleton from the previous chapter, using the Tersoff potential for carbon:

pair_style tersoff/table
pair_coeff * * SiCGe.tersoff C

SiCGe.tersoff holds Tersoff parameters for Si, C and Ge (Tersoff, Phys. Rev. B 39, 5566, 1989). The C at the end maps atom type 1 to the carbon entry. tersoff/table is a faster, tabulated version of the standard tersoff style.

Run it

lmp -in graphene.lmp

The sheet has 1600 atoms ($20 \times 20$ cells $\times$ 4 atoms per cell); if you see a different number, check your L and W. The full run is 200,000 steps (50 ps) and takes about a minute on a few cores. The temperature should settle around 300 K.

Open graphene_data.txt in OVITO, add Create bonds, and load graphene_dump.txt as a trajectory to watch the sheet vibrate.

################### Initialization ###################
variable project string graphene
variable t equal 300.0 # temperature (K)

log           ${project}_log.txt
units           metal # grams/mole, angstroms, ps, eV, K, bars, gm/cm^3
dimension      2
boundary       p p p
atom_style     atomic

# Building the graphene structure

variable l equal 1.42 # bond length in angstroms
variable m  equal 20 # number of unit cells along length
variable n  equal 20 # number of unit cells along width
variable L equal $(sqrt(3)*v_l*v_m) # length of the graphene sheet
variable W equal $(3*v_l*v_n) # width of the graphene sheet

region box block 0 ${L} 0 ${W} -10 10 units box
create_box 1 box
mass 1 12.01 # mass of carbon atom

lattice custom ${l} &
  a1 $(sqrt(3)) 0.0 0.0 &
  a2 0.0 3.0 0.0 &
  basis 0.0 0.0 0.0 &
  basis $(1/2) $(1/6) 0.0 &
  basis $(1/2) $(1/2) 0.0 &
  basis 0.0 $(2/3) 0.0
region mybox block 0 ${L} 0 ${W} -0.1 0.1 units box
create_atoms 1 region mybox

velocity all create $t 12345 dist gaussian

# Definition of interatomic potential
pair_style tersoff/table
pair_coeff * * SiCGe.tersoff C

# Simulation settings
timestep 0.00025
neighbor 3.0 bin
neigh_modify delay 0 every 1 check yes

min_style cg
minimize 1e-10 1e-12 10000 100000


write_data ${project}_data.txt

# Equilibration
fix 1 all nvt temp $t $t $(100.0*dt)
thermo 100
thermo_style custom step temp etotal press
dump 1 all custom 100 ${project}_dump.txt id type x y z vx vy vz
run 200000

Download the graphene folder (.zip): the script and SiCGe.tersoff.

Try this

1. Change m and n to make a long, narrow ribbon. Check that the atom count is $4mn$.
2. Set m to a non-integer value, such as 20.5. Look at the edge of the sheet in OVITO, and compare the starting energy in the log.
3. Change dimension 2 to dimension 3 and the box z range to −20 to 20. In 3D, the sheet is free to ripple out of plane.

Going further

This sheet is the starting point for real graphene studies. My LAMMPS code samples include a script for thermal rectification in graphene, where heat flows more easily in one direction than the other.