school
Custom Lattices and Models
LM-01 Honeycomb Antiferromagnet

LM-01 Honeycomb Antiferromagnet

Defining a Custom Lattice: the Honeycomb Antiferromagnet

In this tutorial you will learn how to define a lattice geometry that is not part of the standard ALPS library, wire it to a quantum spin model, and run a Monte Carlo simulation on the result. As the concrete example we use the honeycomb (hexagonal) lattice, which underlies graphene, the Kitaev honeycomb model, and many frustrated magnets.

By the end of this tutorial you will be able to:

  1. Write an ALPS lattice library XML file from scratch.
  2. Reference that file in a simulation parameter file.
  3. Run the classical Heisenberg model on your lattice with spinmc.
  4. Verify your definition by comparing against the built-in ALPS honeycomb lattice.

Background: the honeycomb lattice

The honeycomb lattice is a two-dimensional lattice with a triangular Bravais lattice and a two-site basis. Each unit cell contains two inequivalent sites — conventionally called the A and B sublattices — and each site has exactly three nearest neighbors. Because A-sites bond only to B-sites and vice versa, the honeycomb is bipartite, making it a natural host for antiferromagnetic order.

The primitive lattice vectors are

$$ \mathbf{a}_1 = a(1, 0), \qquad \mathbf{a}_2 = a\left(\tfrac{1}{2}, \tfrac{\sqrt{3}}{2}\right), $$

where $a$ is the lattice constant (set to 1 below). The two basis sites sit at

$$ \text{A}: \tfrac{2}{3}\mathbf{a}_1 - \tfrac{1}{3}\mathbf{a}_2 = \left(\tfrac{2}{3}, -\tfrac{1}{3}\right), \qquad \text{B}: \tfrac{1}{3}\mathbf{a}_1 + \tfrac{1}{3}\mathbf{a}_2 = \left(\tfrac{1}{3}, \tfrac{1}{3}\right) $$

when expressed in fractional coordinates of the unit cell. These coordinates place the center of the hexagon at the origin of the unit cell.

Step 1: Write the lattice library file

ALPS reads lattice definitions from an XML file with a <LATTICES> root element. The file must be bare XML — no <?xml ...?> declaration and no <!DOCTYPE> preamble — so that the ALPS parser can read it directly.

Create a file called my_honeycomb.xml with the following contents:

<LATTICES>

  <!-- 1. Bravais lattice: triangular, basis vectors a1 and a2 -->
  <LATTICE name="triangular" dimension="2">
    <BASIS>
      <VECTOR>1 0</VECTOR>
      <VECTOR>0.5 0.8660254037844387</VECTOR>
    </BASIS>
  </LATTICE>

  <!-- 2. Unit cell: 2 sites (A and B sublattices) and 3 nearest-neighbor bonds -->
  <UNITCELL name="honeycomb" dimension="2" vertices="2">
    <VERTEX id="1" type="0"><COORDINATE>0.6666666666666666 -0.3333333333333333</COORDINATE></VERTEX>
    <VERTEX id="2" type="1"><COORDINATE>0.3333333333333333 0.3333333333333333</COORDINATE></VERTEX>
    <EDGE type="0"><SOURCE vertex="1"/><TARGET vertex="2"/></EDGE>
    <EDGE type="1"><SOURCE vertex="2"/><TARGET vertex="1" offset="0 1"/></EDGE>
    <EDGE type="2"><SOURCE vertex="1"/><TARGET vertex="2" offset="1 -1"/></EDGE>
  </UNITCELL>

  <!-- 3. Finite L x W honeycomb lattice with periodic boundary conditions -->
  <LATTICEGRAPH name="honeycomb L x W">
    <FINITELATTICE>
      <LATTICE ref="triangular"/>
      <PARAMETER name="W" default="L"/>
      <EXTENT dimension="1" size="L"/>
      <EXTENT dimension="2" size="W"/>
      <BOUNDARY type="periodic"/>
    </FINITELATTICE>
    <UNITCELL ref="honeycomb"/>
  </LATTICEGRAPH>

</LATTICES>

What each element does

<LATTICE name="triangular"> defines the infinite Bravais lattice by giving its basis vectors. The second basis vector $(0.5,, \sqrt{3}/2)$ tilts at 60° relative to the first, producing the triangular tiling.

<UNITCELL name="honeycomb"> places the two basis sites and declares the three edges that connect them to nearest neighbors:

Edge typeSourceTargetMeaning
0A (cell 0,0)B (cell 0,0)intra-cell bond
1B (cell 0,0)A (cell 0,1)bond along $\mathbf{a}_2$
2A (cell 0,0)B (cell 1,−1)bond along $\mathbf{a}_1 - \mathbf{a}_2$

Together, these three edges per unit cell give every A-site three B-neighbors and every B-site three A-neighbors.

The vertex coordinates are given in fractional coordinates of the unit cell (i.e., as coefficients of the Bravais lattice vectors).

<LATTICEGRAPH name="honeycomb L x W"> assembles the finite simulation cell. Setting W to default to L means that specifying only L=4 in the parameter file produces a square $4\times 4$ arrangement of unit cells, giving $4 \times 4 \times 2 = 32$ sites total. <BOUNDARY type="periodic"/> wraps both dimensions into a torus.

Step 2: Write the parameter file

Create a parameter file parm_honeycomb in the same directory as my_honeycomb.xml:

LATTICE_LIBRARY="my_honeycomb.xml"
LATTICE="honeycomb L x W"
MODEL="Heisenberg"
J=-1
L=4
W=4
UPDATE="cluster"
THERMALIZATION=5000
SWEEPS=50000
{T=0.1;}
{T=0.3;}
{T=0.5;}
{T=0.8;}
{T=1.0;}
{T=1.5;}
{T=2.0;}
{T=3.0;}

Key points:

  • LATTICE_LIBRARY tells ALPS where to find your custom lattice definitions.
  • LATTICE must match exactly the name you gave in <LATTICEGRAPH name="...">.
  • MODEL="Heisenberg" uses the built-in classical Heisenberg model. J=-1 makes the coupling ferromagnetic in the spinmc convention; the susceptibility peak behavior is symmetric, so you can use $J = \pm 1$.
  • L=4, W=4 selects a $4 \times 4$ unit-cell tile (32 sites).
  • UPDATE="cluster" uses Wolff cluster updates, which are much more efficient than local updates near the ordering temperature.

Step 3: Run the simulation

Convert the parameter file to XML and run spinmc:

parameter2xml parm_honeycomb
spinmc --Tmin 5 --write-xml parm_honeycomb.in.xml

The simulation runs eight tasks in sequence (one per temperature), writing results to parm_honeycomb.task1.out.xml through parm_honeycomb.task8.out.xml.

Step 4: Extract results

The Python snippet below reads the susceptibility and energy density from the output files:

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

data = pyalps.loadMeasurements(
    pyalps.getResultFiles(prefix='parm_honeycomb'),
    ['Susceptibility', 'Energy'])

chi = pyalps.collectXY(data, x='T', y='Susceptibility')
energy = pyalps.collectXY(data, x='T', y='Energy')

plt.figure()
pyalps.plot.plot(chi)
plt.xlabel('Temperature $T/|J|$')
plt.ylabel('Susceptibility $\\chi$')
plt.title('Classical Heisenberg on honeycomb (L=4)')
plt.show()

Expected output (values will vary slightly due to Monte Carlo noise):

  T |          chi |       E/site
------------------------------------
0.1 |     0.111780 |    -1.401380
0.3 |     0.122365 |    -1.191363
0.5 |     0.134909 |    -0.941824
0.8 |     0.145397 |    -0.595319
1.0 |     0.140822 |    -0.478350
1.5 |     0.120825 |    -0.326241
2.0 |     0.105280 |    -0.247001
3.0 |     0.080478 |    -0.165197

The susceptibility $\chi(T)$ peaks near $T \approx 0.8,|J|$, consistent with the onset of short-range antiferromagnetic correlations in the classical Heisenberg honeycomb model. At high $T$, $\chi$ decreases towards the Curie law $\chi \propto 1/T$, while the energy density approaches zero.

Step 5: Verify against the built-in lattice

ALPS already ships a "honeycomb lattice" definition in its standard library. You can use it as a quick check that your custom XML file is correct:

LATTICE="honeycomb lattice"
MODEL="Heisenberg"
J=-1
L=4
W=4
UPDATE="cluster"
THERMALIZATION=5000
SWEEPS=50000
{T=0.5;}
{T=1.0;}
{T=2.0;}

Note the absence of LATTICE_LIBRARY — the built-in lattice needs no extra file. Running spinmc on this parameter file should produce susceptibilities that agree with your custom-lattice results within statistical uncertainty (a few percent).

Questions and extensions

  1. Coordination number: The honeycomb lattice has coordination number $z = 3$. How does the susceptibility peak temperature compare with the square lattice ($z = 4$)? Try replacing LATTICE="honeycomb L x W" with LATTICE="square lattice" (no custom file needed) and repeating the calculation.

  2. Finite-size effects: Increase L from 4 to 6 and 8. How do the susceptibility and energy density change? For what range of temperatures is a $4 \times 4$ cell a good approximation to the thermodynamic limit?

  3. Bond anisotropy: The unit cell defines three edge types (0, 1, 2). Modify the Hamiltonian to assign different coupling strengths $J_0$, $J_1$, $J_2$ to each bond type. In ALPS the spin model accepts per-bond-type couplings via the J0, J1, … parameters — see the model definition documentation. What happens to the susceptibility when $J_2 = 0$ (reducing to a chain of dimers)?

  4. Quantum model: Replace MODEL="Heisenberg" with MODEL="spin", add local_S=0.5 and ALGORITHM="loop", and switch the executable to loop. This runs the quantum $S = 1/2$ Heisenberg antiferromagnet on your custom honeycomb lattice. Compare the susceptibility curves for the classical and quantum models.