ED-06 Full Diagonalization

In this tutorial we will use the fulldiag application to calculate the complete spectrum of small quantum spin systems, and use it to compute exact finite-temperature thermodynamic quantities — energy, entropy, specific heat, and susceptibility — as a weighted sum over every eigenstate, entirely free of the statistical noise of quantum Monte Carlo or the truncation errors of DMRG. This complements the previous tutorials, which used sparsediag to target only the lowest few eigenstates of a Hamiltonian: here we need the entire spectrum, because every state, however high in energy, carries Boltzmann weight $e^{-E_n/T}$ at any finite temperature.

Thermodynamics for one-dimensional spin models

The spin-1 Heisenberg chain

Our first simulation will be for a spin-1 Heisenberg quantum antiferromagnetic chain,

on a periodic chain, whose finite-chain thermodynamics were first studied systematically by J.C. Bonner and M.E. Fisher, Phys. Rev. 135, A640 (1964).

Parameters

Lattice

The same periodic chain lattice used throughout this tutorial series (see the ALPS lattice library):

    J     J     J     J     J     J     J
o-------o-------o-------o-------o-------o-------o-------o
0       1       2       3       4       5       6       7
|_______________________________________________________|
                      J   (bond 7-0, periodic)

Method

Thermodynamic quantities require summing Boltzmann weights over the entire spectrum, not just the ground state, so fulldiag is used instead of sparsediag. The full Hilbert space of this S=1, L=8 chain has dimension $3^8=6561$; using CONSERVED_QUANTUMNUMBERS=Sz splits it into 17 much smaller blocks (of at most $1107$ states, the $S_z=0$ sector) that are each diagonalized directly, which is both faster and enough to reconstruct the exact partition function at any temperature.

Using the command line

The parameter file parm6a sets up a full diagonalization for a spin-1 Heisenberg quantum antiferromagnet on a one-dimensional chain with 8 sites.

LATTICE="chain lattice"
MODEL="spin"
local_S = 1
J       = 1
CONSERVED_QUANTUMNUMBERS="Sz"
{L = 8}

The CONSERVED_QUANTUMNUMBERS parameter helps fulldiag split the Hilbert space into invariant subspaces and diagonalize them separately. Using the standard sequence of commands you can first convert the input parameters to XML and then compute the full spectrum of this quantum Hamiltonian using fulldiag:

parameter2xml parm6a
fulldiag --write-xml parm6a.in.xml

The output file now contains results for all the eigenvectors, and you can produce XML plot files for the thermodynamic and magnetic observables using fulldiag_evaluate:

fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T  0.1  parm6a.task1.out.xml

This will generate the following XML plot files:

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_susceptibility.xml
parm6a.task1.plot.magnetization.xml

To extract the calculated results from the XML plot files generated by fulldiag_evaluate, you can use the plot2text tool, and then view this data with your favorite plotting tool. For example, to extact the data of the energy density vs. temperature, use

plot2text parm6a.task1.plot.energy.xml

In a similar way, you can extract the data from the other XML plot files. When Grace is your favorite plotting tool, you can also directly generate a Grace project file from the XML plot file using the plot2xmgr tool. For example, to generate a Grace project file of the energy vs. temperature, use

plot2xmgr parm6a.task1.plot.energy.xml > energy.agr

Similarly the tool plot2gp produces Gnuplot scripts and plot2text converts the file to plain text. However, the preferred method for data evaluation and plotting is using Python.

Using Python

To set up and run the simulation in Python we use the script tutorial6a.py. The first parts of this script imports the required modules and then prepares the input files as a list of Python dictionaries, and then runs the simulation

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
parms = [{ 
        'LATTICE'                   : "chain lattice", 
        'MODEL'                     : "spin",
        'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
        'local_S'                   : 1,
        'J'                         : 1,
        'L'                         : 8
    }]

input_file = pyalps.writeInputFiles('parm6a',parms)
res = pyalps.runApplication('fulldiag',input_file)

We next run the evaluation program on all output files

data = pyalps.evaluateFulldiagVersusT(pyalps.getResultFiles(prefix='parm6a'),DELTA_T=0.1, T_MIN=0.1, T_MAX=10.0)

and finally show all the plots:

for s in pyalps.flatten(data):
    plt.figure()
    plt.title("Antiferromagnetic Heisenberg chain")
    pyalps.plot.plot(s)
    plt.show()

Output data

Diagonalizing the same Hamiltonian independently gives the following per-site thermodynamic quantities (energy $E/L$, specific heat $C/L$, uniform susceptibility $\chi/L$):

The specific heat has a broad maximum around $T\approx J$, while the susceptibility peaks near $T\approx1$–$1.5\,J$ and then drops sharply towards $T=0$ — the finite-size remnant of the exponentially activated behaviour $\chi\sim e^{-\Delta/T}$ expected from the Haldane gap $\Delta\approx0.41\,J$ computed in ED-02.

Exercises

• Repeat the computation of the susceptibility, but for L=9 (this will take a few minutes, use L=7 if you are impatient), and compare the result with the one for L=8 ! Recall that all steps starting with parameter2xml … must be repeated.

• Give a rough estimate of the temperature range in which one may use the finite-size result for an infinite chain.

The spin-1/2 Heisenberg ladder

We now change the model to a two-leg ladder,

of length L=6 (i.e. 6 rungs, 12 sites), setting the parameter LATTICE to “ladder” and the couplings J0 and J1 to 1 — the same lattice used in ED-03:

      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 0
      |        |        |        |        |        |
      J1       J1       J1       J1       J1       J1
      |        |        |        |        |        |
      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 1

The parameters are in parm6b, the Python script in tutorial6b.py. Run it exactly as parm6a above:

parameter2xml parm6b
fulldiag --write-xml parm6b.in.xml
fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6b.task1.out.xml

parm6b is otherwise identical to parm6a, but with LATTICE="ladder", J0=1, J1=1, and {L = 6} in place of the chain parameters.

Output data

Diagonalizing this L=6 ladder independently (with open leg boundaries) gives a finite-size singlet-triplet gap of $\Delta\approx0.68\,J$, somewhat above the accepted infinite-ladder value of $\approx 0.5\,J$ quoted in the hint below — a reminder that a 6-rung system is still fairly far from the thermodynamic limit for a gapped state. The corresponding per-site thermodynamics:

Exercises

• Discuss the position of the maximum of the specific heat (hint: the infinite ladder has a gap of approximately J/2)!
• Repeat this computation for 7 rungs (or 5 rungs, if you are impatient).
• What can you infer from the comparison about the temperature range where the finite-size results are a good approximation to the infinite system?

Remarks: If you are familiar with Quantum Monte Carlo (QMC) simulations, you will know that larger systems can be treated and hence better approximations to the thermodynamic limit obtained. (Try it ! You may find that it is in fact not so easy to do better with QMC than with full diagonalization for the above examples.)

Exact diagonalization definitely remains the method of choice under certain conditions. Firstly, if the system is intrinsically finite (and possibly small), full diagonalization yields the exact result at once. Secondly, exact diagonalization does not suffer from the sign problem such that it can also be used for fermionic or frustrated models without any additional constraint. Both conditions are simultaneously satisfied in the models for magnetic molecules which will be discussed in the second part.

Thermodynamics of magnetic molecules

We will now use fulldiag to simulate the exact thermodynamics of small spin clusters, by calculating their complete spectrum.

Two coupled dimers

Setting up a custom lattice/graph and model

Consider the following system of two coupled dimers: an “upper” dimer of two spins $S_0$ coupled by $J_0$, a “lower” dimer of two spins $S_1$ also coupled by $J_0$, and every spin of the upper dimer coupled to every spin of the lower dimer by $J_1$,

        J0
   1 ---------- 2        (upper dimer, spin S0)
   |  \        /  |
   |   \  J1  /   |
   J1    \  /    J1
   |     /  \     |
   |   /      \   |
   3 ---------- 4        (lower dimer, spin S1)
        J0

This is precisely the graph encoded below: sites 1,2 (type 0) form the upper dimer, sites 3,4 (type 1) the lower dimer, edges of type 0 are the intra-dimer bonds ($J_0$) and edges of type 1 are the four inter-dimer bonds ($J_1$).

First, we need a graph to represent this problem. This is defined by the following entry in the file dd-graph.xml:

<LATTICES>
<GRAPH name="double dimer" vertices="4">
<VERTEX id="1" type="0"></VERTEX>
<VERTEX id="2" type="0"></VERTEX>
<VERTEX id="3" type="1"></VERTEX>
<VERTEX id="4" type="1"></VERTEX>
<EDGE type="0" source="1" target="2"/>
<EDGE type="0" source="3" target="4"/>
<EDGE type="1" source="1" target="3"/>
<EDGE type="1" source="1" target="4"/>
<EDGE type="1" source="2" target="3"/>
<EDGE type="1" source="2" target="4"/>
</GRAPH>
</LATTICES>

Note: the file must use bare XML with no <?xml?> declaration or <!DOCTYPE> header — ALPS’s parser does not support those.

Then we also want to assign different Heisenberg exchanges J0 and J1 to the edges of type 0 and 1, respectively. This is achieved by the following entry in the file model-dspin.xml:

<MODELS>
<SITEBASIS name="spin">
  <PARAMETER name="local_spin" default="local_S"/>
  <PARAMETER name="local_S" default="1/2"/>
  <QUANTUMNUMBER name="S" min="local_spin" max="local_spin"/>
  <QUANTUMNUMBER name="Sz" min="-S" max="S"/>
  <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
    <CHANGE quantumnumber="Sz" change="1"/>
  </OPERATOR>
  <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
    <CHANGE quantumnumber="Sz" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Sz" matrixelement="Sz"/>
</SITEBASIS>
<BASIS name="spin">
  <SITEBASIS ref="spin">
    <PARAMETER name="local_spin" value="local_S#"/>
    <PARAMETER name="local_S#" value="1/2"/>
  </SITEBASIS>
  <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
</BASIS>
<HAMILTONIAN name="dimerized spin">
<PARAMETER name="J" default="1"/>
<PARAMETER name="h" default="0"/>
<BASIS ref="spin"/>
<SITETERM site="i">
<PARAMETER name="h#" default="h"/>
    -h#*Sz(i)
</SITETERM>
<BONDTERM source="i" target="j">
<PARAMETER name="J#" default="J"/>
    J#*Sz(i)*Sz(j)+J#/2*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j))
</BONDTERM>
</HAMILTONIAN>
</MODELS>

Note: the file must use bare XML with no <?xml?> declaration or <!DOCTYPE> header, and must include the <SITEBASIS> and <BASIS> definitions — without them fulldiag cannot construct the Hilbert space.

Note that we do not really need this definition since the “spin” Hamiltonian in the default models.xml file already contains suitable definitions. Nevertheless, the above example illustrates how to automatically assign exchange constants Jn to a bond of type n, using the hash symbol (#).

In passing we have also assigned different types to vertices 1,2 and 3,4. Therefore, we are able to assign different local spins to the upper and lower dimer by specifying the corresponding values via local_S0 and local_S1, respectively.

We will now compute the magnetization curve for local spins $S_0=1$ (upper dimer) and $S_1=1/2$ (lower dimer), $J_0=1$, $J_1=0.4$ at a temperature $T=0.02$.

Using the command line

The parameter file parm6c sets up the simulation:

LATTICE="double dimer"
MODEL="dimerized spin"
LATTICE_LIBRARY="dd-graph.xml"
MODEL_LIBRARY="model-dspin.xml"
local_S0=1
local_S1=1/2
J0      = 1
J1      = 0.4
h       = 0 
CONSERVED_QUANTUMNUMBERS="Sz"
{T = 0.02}

Note the new parameters LATTICE_LIBRARY and MODEL_LIBRARY that point to the files containing our custom lattice and model. The computation is performed with the following sequence of commands:

parameter2xml parm6c
fulldiag --write-xml parm6c.in.xml
fulldiag_evaluate --H_MIN 0 --H_MAX 4 --DELTA_H 0.025 --versus h parm6c.task1.out.xml

Note that here we have used fulldiag_evaluate with commandline arguments to specify a magnetic field range. In particular the commandline argument –versus h puts the magnetic field on the x-axis rather than temperature, as in the previous examples.

Method

With only $(2\cdot1+1)^2\times(2\cdot\tfrac12+1)^2=36$ basis states, this system is trivial for fulldiag — the point of the example is not computational difficulty but showing how to wire up a fully custom lattice graph and Hamiltonian, of the kind needed for real magnetic-molecule data, entirely from parameter files.

Output data

At the low temperature $T=0.02\,J_0$ used here, the total spin $S_{\text{total}}$ of the ground state increases in unit steps as $h$ is raised, each step occurring where two levels from the analytic spectrum below cross:

The three crossing fields follow directly from the analytic energies given in the hint below: equating the lowest energy of $S_{\text{total}}=n$ (at $S^z=n$) to that of $S_{\text{total}}=n+1$ (at $S^z=n+1$) gives $h_1=1.0$, $h_2=1.4$ and $h_3=2.4$, in units of $J_0$ — a magnetization staircase with plateaus at $S^z_{\text{total}}=0,1,2,3$, rounded only by the small but nonzero temperature.

Question

• Plot and interpret the result!

Hint: The spectrum of two coupled $S_0=1$ and $S_1=1/2$ dimers can be found analytically. The energies are (with some degeneracies):

• $-11J_0/4$, $-3J_0/4-2J_1$ for total spin $S_{\text{total}}=0$;
• $-7J_0/4$, $-3J_0/4-J_1$, $5J_0/4-3J_1$ for $S_{\text{total}}=1$;
• $-3J_0/4+J_1$, $J_0/4$, $5J_0/4-J_1$ for $S_{\text{total}}=2$;
• $5J_0/4+2J_1$ for $S_{\text{total}}=3$.

The molecular complex $V_{15}$

The final example is a model for the molecular nanomagnet $\mathrm{V}_{15}$, in which 15 $\mathrm{V}^{4+}$ ions ($S=1/2$ each) sit on two hexagons sandwiching a central triangle. Its low-energy magnetism is dominated by strong antiferromagnetic exchange around the two hexagons and weaker exchange linking them to the central triangle, so that at low temperature the 15 spins behave effectively like a frustrated spin-1/2 triangle: two low-lying doublets plus a quadruplet, well separated from the rest of the spectrum (G. Chaboussant et al., Europhys. Lett. 59, 291 (2002)).

        1---2                 outer hexagon (6 spins)
       /     \
      6       3
       \     /
        5---4
       /     \
      o       o                inner hexagon (6 spins),
       \     /                 rotated relative to the outer one
        o---o
          |
        triangle (3 spins)    weakly coupled to both hexagons

For the simulations we make the simplifying assumption of equal Heisenberg exchange $J$ along all edges of this graph. Setting up the custom lattice and model files follows exactly the same recipe as the coupled-dimer example above — a <GRAPH> with 15 vertices and the corresponding <EDGE> list for the hexagon and triangle bonds, and the same dimerized spin (or the built-in spin) Hamiltonian. The associated graph is defined in v15-graph.xml. Running the simulations now proceeds as above. The parameters are in parm6d, the Python script in tutorial6d.py.

Method

The full Hilbert space of 15 spin-1/2 sites has dimension $2^{15}=32768$; even resolved into $S_z$ sectors, fulldiag has to diagonalize a block of a few thousand states exactly to get the low-temperature thermodynamics right, which is why — as the note below says — this run takes noticeably longer than the previous examples.

Question

• How do you explain the low-temperature behavior of the magnetic susceptibility? (Hint: think about the effective low-energy spin-1/2 triangle mentioned above, and how its own level crossings in a field would show up in $\chi(T)$.)

Note that this computation is already a bit challenging and will take a while. So, it may be a good idea to have a break before you come back and look at the result.

Summary

Full diagonalization gives numerically exact thermodynamics for any finite quantum spin system, from translationally invariant chains and ladders to fully custom, frustrated magnetic molecules — the price is the exponential growth of the Hilbert space with system size, which is why the method is confined to some tens of spins even when symmetries are exploited, in contrast to the far larger systems reachable with sparsediag for ground-state and low-energy properties alone.

Additional Exercises

Hubbard model on a square lattice

• Set up a parameter file for the Hubbard model on a square lattice.
  • Find the lattice that you would like to use. You could use either a square lattice or a chain. Pay attention to the boundary conditions!
  • Find the Hubbard model in the model library. Make sure you understand its terms.
  • Switch on the symmetries: conserved are momenta (in x and y direction) and particles.
  • Pick trial parameters that you can check (e.g. t=0, or U=0), and run the simulation.
  • How would you introduce a t'?