LM-02 Custom Model Hamiltonians

Defining Custom Model Hamiltonians: Spin, Fermion, and Boson

In tutorial LM-01 you learned how to define a custom lattice. Here you will go one level deeper and write a model library XML file that defines a quantum Hamiltonian from its elementary building blocks — local Hilbert spaces and interaction terms. This lets you simulate any lattice model, not just the handful that ALPS ships by default.

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

1. Describe the four XML building blocks: SITEBASIS, QUANTUMNUMBER, OPERATOR, and HAMILTONIAN.
2. Construct the spin-1 Heisenberg chain with single-ion anisotropy from scratch and extract its spin gap with sparsediag.
3. Write analogous model files for spinless fermions and soft-core bosons.
4. Understand how the matrix-element formulas differ between particle types.

The ALPS model library XML format

A model library file is a bare XML document (no <?xml ...?> preamble) with a <MODELS> root element. A complete model is built in two steps:

Step 1 — SITEBASIS: the local Hilbert space

<SITEBASIS name="my basis">
  <PARAMETER name="P" default="1"/>              <!-- tunable parameter, e.g. spin S -->
  <QUANTUMNUMBER name="Q" min="-P" max="P"
                 type="half_integer"/>           <!-- conserved quantum number       -->
  <OPERATOR name="Raise"
            matrixelement="sqrt(P*(P+1)-Q*(Q+1))">  <!-- off-diagonal raising operator -->
    <CHANGE quantumnumber="Q" change="1"/>
  </OPERATOR>
  <OPERATOR name="Lower"
            matrixelement="sqrt(P*(P+1)-Q*(Q-1))">
    <CHANGE quantumnumber="Q" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Diag" matrixelement="Q"/>      <!-- diagonal operator              -->
</SITEBASIS>

Key rules:

• matrixelement is an algebraic expression evaluated at the source quantum-number values (before the CHANGE is applied).
• <CHANGE quantumnumber="Q" change="±1"/> shifts Q by the given amount; ALPS automatically enforces the min/max bounds and sets the matrix element to zero outside that range.
• type="half_integer" allows half-integer quantum-number values (needed for spin $S=\tfrac{1}{2}$, $\tfrac{3}{2}$, …).

Step 2 — BASIS: quantum-number constraints

BASIS is the connector between a SITEBASIS and a HAMILTONIAN. It declares which quantum numbers are conserved and what simulation parameter controls each sector.

<BASIS name="my basis">
  <SITEBASIS ref="my basis"/>
  <CONSTRAINT quantumnumber="Q" value="Q_total"/>  <!-- fix total Q to parm Q_total -->
</BASIS>

Step 3 — HAMILTONIAN: on-site and bond terms

<HAMILTONIAN name="my model">
  <PARAMETER name="J" default="1"/>
  <BASIS ref="my basis"/>             <!-- reference the BASIS, not SITEBASIS -->
  <SITETERM site="i">
    D*Diag(i)*Diag(i)                <!-- on-site expression                 -->
  </SITETERM>
  <BONDTERM source="i" target="j">
    J*Diag(i)*Diag(j)                <!-- bond expression                    -->
  </BONDTERM>
</HAMILTONIAN>

Common mistake: writing <SITEBASIS ref="..."/> inside <HAMILTONIAN> causes ALPS to abort with “unexpected element”. Always use <BASIS ref="..."/> inside <HAMILTONIAN>, and define a separate <BASIS> element that wraps the <SITEBASIS>.

The complete file has the structure:

<MODELS>
  <!-- SITEBASIS elements (local Hilbert spaces)   -->
  <!-- BASIS elements (constraint declarations)    -->
  <!-- HAMILTONIAN elements (interaction terms)    -->
</MODELS>

Reference it in the parameter file with:

MODEL_LIBRARY = "my_models.xml"
MODEL         = "my model"

Example 1 — Spin: S=1 chain with single-ion anisotropy

Physics background

The Hamiltonian is

where $J > 0$ is the antiferromagnetic exchange and $D$ is the single-ion anisotropy. Writing $\mathbf{S}_i \cdot \mathbf{S}_j$ in raising/lowering operators:

At $D = 0$ the chain is in the Haldane phase: a gapped, symmetry-protected topological phase with a bulk spin gap $\Delta \approx 0.41\,J$ (Haldane 1983; White & Huse 1993). Increasing $D$ closes the gap at a quantum critical point $D_c \approx 0.97\,J$ and drives the system into the large-$D$ phase, whose ground state is close to the product state $|0, 0, \ldots, 0\rangle$ (all sites in the $S^z = 0$ state).

Hamiltonian skeleton (open chain, 4 sites shown):

  D       D       D       D
  o --J-- o --J-- o --J-- o

D  = single-ion anisotropy  D*(Sz)^2      [SITETERM]
J  = Heisenberg exchange                   [BONDTERM]
h  = uniform magnetic field -h*Sz         [SITETERM, default 0]

Lattice choice

We use the chain lattice (ALPS built-in, periodic boundary conditions) with $L=8$ sites. Periodic BC lets sparsediag exploit translational symmetry and block-diagonalize by total momentum $k$ as well as $S^z$, giving a cleaner finite-size view of the bulk gap. See the ALPS lattice library for the full list of available geometries.

Method choice

sparsediag (sparse Lanczos exact diagonalization) is exact and well-suited to systems with $\lesssim 10^6$ basis states. For $L=8$, $S=1$, $S^z_{\text{tot}}=0$ the Hilbert space has $\binom{16}{8} = 12\,870$ states — tiny for Lanczos. Exact diagonalization is the right choice here because we want precise eigenvalues to measure the gap to four decimal places as a function of $D$.

Step 1: Write the model library file

Create my_models.xml:

<MODELS>

  <!-- ===========================================================
       Site basis: spin-S
         * local_S sets the spin magnitude (default 1 for this tutorial)
         * S is the per-site quantum number fixed to local_S
         * Sz runs from -S to S in integer steps
         * Operators: Splus, Sminus (raising/lowering), Sz (diagonal)
  ============================================================ -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>

    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S"      max="S"       type="half_integer"/>

    <!-- S^+ raises Sz by 1; matrix element from Wigner-Eckart theorem -->
    <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
      <CHANGE quantumnumber="Sz" change="1"/>
    </OPERATOR>

    <!-- S^- lowers Sz by 1 -->
    <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
      <CHANGE quantumnumber="Sz" change="-1"/>
    </OPERATOR>

    <!-- Sz is diagonal: matrix element equals the current value of Sz -->
    <OPERATOR name="Sz" matrixelement="Sz"/>
  </SITEBASIS>

  <!-- BASIS links the site basis to a total-Sz constraint -->
  <BASIS name="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <!-- ===========================================================
       Hamiltonian: Heisenberg exchange + single-ion anisotropy
         H = J * sum_{<ij>} (Sz(i)*Sz(j) + 0.5*(S+(i)*S-(j) + S-(i)*S+(j)))
           + D * sum_i  Sz(i)^2
           - h * sum_i  Sz(i)      [optional uniform field]
  ============================================================ -->
  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>

    <BASIS ref="anisotropic spin chain"/>

    <!-- on-site: single-ion anisotropy + field -->
    <SITETERM site="i">
      D*Sz(i)*Sz(i) - h*Sz(i)
    </SITETERM>

    <!-- bond: Heisenberg exchange written in S+/S-/Sz -->
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

What each element does

<SITEBASIS name="my spin"> defines the local Hilbert space. A spin-$S$ site has $2S+1$ states labeled by $S^z \in \{-S,\, -S+1,\, \ldots,\, S\}$. With the default local_S=1 each site has three states: $|-1\rangle$, $|0\rangle$, $|{+}1\rangle$.

<PARAMETER name="local_S"> sets the spin magnitude. ALPS requires the spin quantum number (S) and the tunable parameter to have different names to avoid a name clash at runtime; following the ALPS convention, the tunable parameter is local_S while the quantum number is S.

<QUANTUMNUMBER name="S" min="local_S" max="local_S" ...> pins $S$ to the value of the local_S parameter, so every site carries the same total spin. <QUANTUMNUMBER name="Sz" ...> registers $S^z$ as a running quantum number that ALPS uses to block-diagonalize the Hamiltonian.

<BASIS name="..."> wraps the site basis and adds a <CONSTRAINT> that fixes the total $S^z$ to the simulation parameter Sz_total. Inside <HAMILTONIAN>, write <BASIS ref="..."/> (not <SITEBASIS ref="..."/>); ALPS will report an error if you use the wrong tag.

<OPERATOR name="Splus" matrixelement="..."> defines the raising operator $S^+$. The matrixelement gives $\langle S,\, S^z{+}1 \mid S^+ \mid S,\, S^z \rangle = \sqrt{S(S+1)-S^z(S^z+1)}$, which is evaluated at the source value of $S^z$. The child <CHANGE quantumnumber="Sz" change="1"/> tells ALPS that $S^+$ increases $S^z$ by 1; ALPS automatically clamps the matrix element to zero when $S^z = S$ (no room to raise further).

<SITETERM site="i"> specifies the on-site Hamiltonian. The expression D*Sz(i)*Sz(i) evaluates the diagonal operator $S^z_i$ twice on site $i$, giving the single-ion anisotropy $D(S^z_i)^2$.

<BONDTERM source="i" target="j"> specifies the two-site bond interaction. The product Splus(i)*Sminus(j) denotes $S^+_i \otimes S^-_j$; ALPS builds the tensor-product matrix automatically. The factor of $\tfrac{1}{2}$ arises from $S^+S^- = 2(S^xS^x + S^yS^y)$.

Step 2: Write the parameter file

To resolve the spin gap $\Delta = E_1(S^z_{\text{tot}}=1) - E_0(S^z_{\text{tot}}=0)$ you run sparsediag in each $S^z_{\text{tot}}$ sector separately. Create parm_spin1:

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "chain lattice"
MODEL                    = "anisotropic spin chain"
L                        = 8
J                        = 1
NUMBER_EIGENVALUES       = 1
CONSERVED_QUANTUMNUMBERS = "Sz"

{D=0.0; Sz_total=0;}
{D=0.0; Sz_total=1;}
{D=1.0; Sz_total=0;}
{D=1.0; Sz_total=1;}
{D=2.0; Sz_total=0;}
{D=2.0; Sz_total=1;}

Key points:

• LATTICE="chain lattice" — standard ALPS periodic chain; no LATTICE_LIBRARY needed. With periodic boundary conditions sparsediag also exploits translational symmetry (diagonalizing each total-momentum sector independently), which gives a cleaner view of the bulk Haldane gap than open boundaries.
• CONSERVED_QUANTUMNUMBERS="Sz" activates block diagonalization; ALPS only diagonalizes the block with Sz_total equal to the requested value.
• Six tasks are created: two sectors ($S^z_{\text{tot}} = 0,\, 1$) at three values of $D$.
• NUMBER_EIGENVALUES=1 requests only the lowest eigenvalue in each momentum block.

Step 3: Run the simulation

parameter2xml parm_spin1
sparsediag --write-xml parm_spin1.in.xml

sparsediag runs six tasks in sequence and writes results to parm_spin1.task1.out.xml through parm_spin1.task6.out.xml.

Step 4: Extract the spin gap

import xml.etree.ElementTree as ET, glob

# With periodic BC, sparsediag diagonalizes each total-momentum sector
# separately and writes one <EIGENVALUES> block per sector.
# Taking the minimum over all sectors gives the ground state energy.
results = {}
for fname in sorted(glob.glob('parm_spin1.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    D  = float(params['D'])
    Sz = int(params['Sz_total'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results.setdefault(D, {})[Sz] = min(means)

print(f"{'D/J':>5}  {'E(Sz=0)':>12}  {'E(Sz=1)':>12}  {'gap Δ/J':>9}  phase")
print('─' * 60)
for D in sorted(results):
    E0 = results[D][0]
    E1 = results[D][1]
    gap = E1 - E0
    phase = 'Haldane' if D < 0.97 else 'large-D'
    print(f"{D:5.1f}  {E0:12.5f}  {E1:12.5f}  {gap:9.5f}  {phase}")

Expected output (L=8, periodic chain, $J=1$):

  D/J      E(Sz=0)      E(Sz=1)    gap Δ/J  phase
────────────────────────────────────────────────────────────
  0.0   -11.33696   -10.74340    0.59356  Haldane
  1.0    -7.00947    -6.62505    0.38442  large-D
  2.0    -4.33528    -3.66645    0.66883  large-D

At $D = 0$ the gap is $\Delta \approx 0.59\,J$ for $L=8$; the bulk value is $0.41\,J$ and the finite-size gap converges from above as $L$ grows (see extension question 1). Near $D = J$ the gap has narrowed to $0.38\,J$ — the finite-size precursor of the quantum phase transition at $D_c \approx 0.97\,J$. At $D = 2J$ the system is in the large-$D$ phase and the gap has grown to $0.67\,J$.

Example 2 — Fermion: Spinless fermion chain

Spinless fermions are the simplest fermionic model: each site is either empty ($|0\rangle$) or occupied ($|1\rangle$), with no spin degree of freedom. The Hamiltonian

describes a free-hopping band (controlled by $t$) with nearest-neighbor density-density repulsion $V$ and chemical potential $\mu$. When $V > 2t$ the ground state is a charge-density wave; when $V < 2t$ it is a Luttinger liquid (metallic).

Hamiltonian skeleton (open chain, 4 sites shown):

  -μ        -μ        -μ        -μ
   o ------- o ------- o ------- o
      -t,V      -t,V      -t,V

-μ = chemical potential   -μ*n           [SITETERM]
-t = hopping amplitude    -t*(c†c+h.c.)  [BONDTERM]
 V = density repulsion     V*n*n          [BONDTERM]

Add the following to my_models.xml (inside the <MODELS> element):

  <!-- ===========================================================
       Site basis: spinless fermion
         * Occupation N ∈ {0, 1}  (Pauli exclusion from max="1")
         * type="fermion" tells ALPS to apply the Jordan-Wigner
           string for 1D systems automatically
  ============================================================ -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>

    <!-- c† creates a fermion: ⟨1|c†|0⟩ = 1 -->
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- c destroys a fermion: ⟨0|c|1⟩ = 1 -->
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- number operator: diagonal, matrix element = current N -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       Hamiltonian: spinless fermion chain
         H = -t*(cdag(i)*c(j)+h.c.) + V*n(i)*n(j) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>

    <BASIS ref="spinless fermion chain"/>

    <SITETERM site="i">
      -mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>

Why type="fermion" matters

Fermionic operators anticommute: $\{c_i, c^\dagger_j\} = \delta_{ij}$. In one dimension this is handled exactly by the Jordan-Wigner transformation, which maps fermions to hard-core bosons with a non-local sign string. Setting type="fermion" on the SITEBASIS instructs ALPS to insert this string automatically whenever a fermionic operator hops across intermediate sites. Without it, the signs would be wrong and the spectrum would be incorrect.

Lattice choice

We use the open chain lattice (open boundary conditions) with $L=10$ sites. Open BC avoids the fermion sign complications that arise in momentum-sector diagonalization for interacting fermions, and is standard for 1D chain ED. See the ALPS lattice library for further chain variants.

Method choice

sparsediag is exact and efficient here: for $L=10$, $N=5$ the Hilbert space has $\binom{10}{5}=252$ states — essentially instantaneous. We request NUMBER_EIGENVALUES=4 to see both the ground state and the lowest neutral excitations within the half-filled sector.

Parameter file

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "spinless fermion chain"
L                        = 10
t                        = 1
mu                       = 0
NUMBER_EIGENVALUES       = 4
CONSERVED_QUANTUMNUMBERS = "N"

{V=0.0; N_total=5;}
{V=1.0; N_total=5;}
{V=2.0; N_total=5;}
{V=3.0; N_total=5;}

N_total=5 fills half the chain (one fermion per two sites), where the CDW instability is strongest.

Run the simulation

parameter2xml parm_fermion
sparsediag --write-xml parm_fermion.in.xml

Extract and display results

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_fermion.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    V  = float(params['V'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[V] = sorted(means)

print(f"{'V/t':>5}  {'E₀':>10}  {'E₁':>10}  {'E₁−E₀':>8}")
print('─' * 42)
for V in sorted(results):
    Es = results[V]
    print(f"{V:5.1f}  {Es[0]:10.5f}  {Es[1]:10.5f}  {Es[1]-Es[0]:8.5f}")

Expected output (L=10, open chain, $t=1$, $N=5$):

  V/t        E₀          E₁       E₁−E₀
──────────────────────────────────────────
  0.0   -6.02667    -5.45742    0.56926
  1.0   -5.00000    -4.30391    0.69609
  2.0   -4.26419    -3.50162    0.76257
  3.0   -3.74967    -2.97936    0.77031

Summary

The excitation gap within the half-filled sector grows monotonically with $V$: repulsion stiffens the charge sector and increases the cost of neutral (particle-hole) excitations. To measure the true charge gap $\Delta_c = E_0(N+1) + E_0(N-1) - 2E_0(N)$ — which vanishes for $V < 2t$ and opens for $V > 2t$ — add tasks with N_total=4 and N_total=6 to the parameter file (see extension question 3).

Example 3 — Boson: Soft-core Bose-Hubbard chain

The Bose-Hubbard model allows multiple bosons per site (up to Nmax) and describes the superfluid–Mott-insulator transition in optical lattice experiments:

The on-site repulsion $U$ penalizes double (and higher) occupancy; for $U \gg t$ and commensurate filling the ground state is a Mott insulator with a particle-hole gap.

Hamiltonian skeleton (open chain, 4 sites shown):

   U,-μ       U,-μ       U,-μ       U,-μ
    o -------- o -------- o -------- o
         -t         -t         -t

 U  = on-site repulsion  (U/2)*n*(n-1)  [SITETERM]
-μ  = chemical potential  -μ*n           [SITETERM]
-t  = hopping amplitude  -t*(b†b+h.c.)   [BONDTERM]

Add to my_models.xml:

  <!-- ===========================================================
       Site basis: soft-core boson
         * Occupation N ∈ {0, 1, …, Nmax}  (default Nmax=3)
         * b† and b carry sqrt factors from the bosonic commutation
           relation [b, b†] = 1
  ============================================================ -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>

    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>

    <!-- b† raises N by 1; matrix element ⟨N+1|b†|N⟩ = sqrt(N+1) -->
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- b lowers N by 1; matrix element ⟨N-1|b|N⟩ = sqrt(N) -->
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- number operator -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       Hamiltonian: soft-core Bose-Hubbard chain
         H = -t*(bdag(i)*b(j)+h.c.) + U/2*n(i)*(n(i)-1) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"   default="1"/>
    <PARAMETER name="U"   default="0"/>
    <PARAMETER name="mu"  default="0"/>

    <BASIS ref="Bose-Hubbard chain"/>

    <SITETERM site="i">
      U/2*n(i)*(n(i)-1) - mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

Matrix elements: bosons vs fermions vs spins

For soft-core bosons the $\sqrt{N+1}$ factor reflects the enhanced stimulated emission into an already-occupied mode — the essence of Bose-Einstein statistics. Hard-core bosons (occupancy capped at 1) and spinless fermions share the matrix element $= 1$; their distinction is fermionic anticommutation, encoded by type="fermion" on the SITEBASIS.

Lattice choice

We use the open chain lattice with $L=6$ sites. Open BC keeps the Hilbert space manageable and avoids unwanted translational degeneracies. With $N_{\max}=3$ and $L=6$ the Hilbert space size for $N=6$ is manageable for Lanczos ED. See the ALPS lattice library for other geometries (square lattice is commonly used to study the 2D superfluid-Mott phase diagram).

Method choice

sparsediag (Lanczos exact diagonalization) is exact at any $U/t$ and directly gives the lowest eigenvalues in the fixed-$N$ sector. For $L=6$, $N=6$, $N_{\max}=3$ the sector dimension is at most $\binom{6+3}{3}^6$ states screened by $N=6$ — tractable for Lanczos. Quantum Monte Carlo (worm) is preferred for larger systems, but ED is ideal for establishing benchmarks and for small $L$ where sign-problem-free comparisons exist.

Parameter file

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "Bose-Hubbard chain"
Nmax                     = 3
L                        = 6
t                        = 1
NUMBER_EIGENVALUES       = 2
CONSERVED_QUANTUMNUMBERS = "N"

{U=0.0;  mu=0; N_total=6;}
{U=2.0;  mu=0; N_total=6;}
{U=5.0;  mu=0; N_total=6;}
{U=10.0; mu=0; N_total=6;}

N_total=6 gives unit filling (one boson per site) where the Mott gap is largest.

Run the simulation

parameter2xml parm_boson
sparsediag --write-xml parm_boson.in.xml

Extract and display results

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_boson.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    U  = float(params['U'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[U] = sorted(means)

print(f"{'U/t':>5}  {'E₀(N=6)':>12}  {'E₁(N=6)':>12}  {'E₁−E₀':>8}")
print('─' * 48)
for U in sorted(results):
    Es = results[U]
    print(f"{U:5.1f}  {Es[0]:12.5f}  {Es[1]:12.5f}  {Es[1]-Es[0]:8.5f}")

Expected output (L=6, open chain, $t=1$, $N_{\max}=3$, $N=6$):

  U/t      E₀(N=6)      E₁(N=6)     E₁−E₀
────────────────────────────────────────────────
  0.0   -10.26758     -9.55861     0.70897
  2.0    -6.64995     -5.60418     1.04577
  5.0    -3.71673     -1.39313     2.32360
 10.0    -1.96696      4.35707     6.32403

Summary

The excitation gap within the unit-filled sector grows rapidly with $U$: for $U = 0$ the bosons are superfluid and the gap is set by kinetic energy alone ($\approx 0.71\,t$); by $U = 10\,t$ on-site repulsion has made any particle-hole excitation extremely costly ($\approx 6.3\,t$), signalling deep Mott-insulator physics. To measure the true single-particle gap $\Delta = E_0(N+1) + E_0(N-1) - 2E_0(N)$ (which vanishes in the superfluid and opens in the Mott phase), add tasks with N_total=5 and N_total=7 (see extension question 4).

The complete my_models.xml

Putting all three sitebases and Hamiltonians together into one file:

<MODELS>

  <!-- ===== SPIN ===== -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>
    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S"      max="S"       type="half_integer"/>
    <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="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>
    <BASIS ref="anisotropic spin chain"/>
    <SITETERM site="i">D*Sz(i)*Sz(i) - h*Sz(i)</SITETERM>
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== SPINLESS FERMION ===== -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="spinless fermion chain"/>
    <SITETERM site="i">-mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== SOFT-CORE BOSON ===== -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>
    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="U"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="Bose-Hubbard chain"/>
    <SITETERM site="i">U/2*n(i)*(n(i)-1) - mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

Questions and extensions

1. Spin gap vs system size: Run the spin-1 example for $L = 6, 8, 10, 12$ at $D=0$. Plot $\Delta(L)$ vs $1/L$. Does the gap extrapolate to the bulk Haldane value $\Delta_\infty \approx 0.41\,J$ as $L \to \infty$? How do the finite-size corrections compare between open and periodic boundary conditions?

2. Spin-$\tfrac{1}{2}$ vs spin-1: Add local_S=0.5 to parm_spin1 (the SITEBASIS uses local_S as its tunable parameter). The Heisenberg model with $S = \tfrac{1}{2}$ is gapless (a Luttinger liquid). Confirm that the gap extrapolates to zero as $L \to \infty$, in contrast to the $S = 1$ case.

3. CDW gap in spinless fermions: Run parm_fermion with $V = 0, 1, 2, 3$. Extract the charge gap $\Delta_c = E_0(N_{\text{tot}}+1) + E_0(N_{\text{tot}}-1) - 2E_0(N_{\text{tot}})$ by adding tasks with N_total=4 and N_total=6. At which value of $V/t$ does the CDW gap open?

4. Mott gap in bosons: Compute the single-particle gap for the Bose-Hubbard chain as a function of $U/t$ at unit filling. Compare with the spinless-fermion CDW gap — both are charge gaps, but their physical origins (repulsion-driven vs kinetic-energy-driven) differ.

5. Custom lattice + custom model: Combine my_models.xml with the honeycomb lattice from LM-01. Run the spin-1 Heisenberg model on the honeycomb lattice using sparsediag for a small cluster. Does single-ion anisotropy $D$ stabilize a different phase on the frustrated honeycomb than on the chain?

6. Adding a next-nearest-neighbor bond: Add a second <BONDTERM> with type="1" to the anisotropic spin chain Hamiltonian:

   <BONDTERM type="1" source="i" target="j">
     J2*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j)))
   </BONDTERM>

   Use the nnn chain lattice and study the effect of next-nearest-neighbor coupling $J_2$ on the Haldane gap. See the lattice how-tos for how edge type 1 is defined in that lattice.