Solving the 1D XY Model with DMRG

Overview

In this tutorial, we'll build a 1D quantum XY model from scratch using the Model class and solve it using the DMRG algorithm to find the ground state.

The XY model is defined by the Hamiltonian:

$H = -\sum_{\langle i,j \rangle} (J_x \sigma_i^x \sigma_j^x + J_y \sigma_i^y \sigma_j^y)$

where $\sigma_i^x$ and $\sigma_i^y$ are Pauli operators at site $i$, and $J_x$, $J_y$ control the coupling strengths.

Step 1: Set up the lattice and model

First, we create a 1D chain lattice with a specified number of sites and open boundary conditions:

// Create a 1D chain with 20 sites
var L = 20
var phys_dim = 2


// Create a linear lattice with open boundary conditions
var lat = lattice("chain",[L])

// Create a model from the lattice
var M = model(lat)

Step 2: Build the XY Hamiltonian

Next, we construct the XY interactions using the Model API. First add coefficients and neighbourhood rules, then define the terms:

// Define coefficients for XY couplings
var Jx = 1.0
var Jy = 1.0

M.add(coefficient("Jx", Jx))
M.add(coefficient("Jy", Jy))

// Define a nearest-neighbor neighbourhood rule
// For a 1D chain, nearest neighbors are displaced by [1]
M.add(neighbourhood_rule("nn", [1]))

// Add X-X interaction term
M.add(term("Jx", [X, X], "nn"))

// Add Y-Y interaction term
M.add(term("Jy", [Y, Y], "nn"))

Each term is defined by:

• A coefficient name (e.g., "Jx")
• A list of operators: first operates at reference site, second at the neighboring site
• A neighbourhood rule name that defines the displacement

Step 3: Convert model to MPO

Convert the Model to an MPO format suitable for DMRG:

// Convert model to complex MPO
var H = mpo(M)

Step 4: Create initial MPS guess

Create an initial random MPS that will be optimized by DMRG:

// Initial state: random MPS with bond dimension 4
var psi0 = make_random_mps(L, phys_dim, 4, as_complex)

Step 5: Configure DMRG parameters

Set up the DMRG optimization parameters:

var dmrg_options = [
    "schedule": "logarithmic",
    "max_bond_dimension": 64,
    "max_iterations": 100,
    "energy_accuracy": 1e-8,
    "truncation_accuracy": 1e-10,
    "verbose": "sweep"
]

Parameter explanation:

• schedule: "logarithmic" grows bond dimension fast enough for most purpose
• max_bond_dimension: 64 is a reasonable upper bound for this system
• max_iterations: 100 sweeps should be sufficient for convergence
• energy_accuracy: Stop when energy changes by less than 1e-8 per sweep
• verbose: "sweep" prints energy after each full sweep

Step 6: Run DMRG

Execute the DMRG algorithm to find the ground state:

var result = dmrg(psi0, H, dmrg_options)

var ground_energy = result[0]
var ground_state = result[1]

print("Ground state energy per site: " + string(ground_energy / L))

Expected output with verbose: "sweep":

sweep: 0 ,Energy: 10000000000.0000, max bond 4
sweep: 1 ,Energy: -21.8996, max bond 8
sweep: 2 ,Energy: -24.4139, max bond 14
sweep: 3 ,Energy: -24.6527, max bond 21
sweep: 4 ,Energy: -24.7203, max bond 27
sweep: 5 ,Energy: -24.7509, max bond 31
sweep: 6 ,Energy: -24.7609, max bond 32
sweep: 7 ,Energy: -24.7628, max bond 33
sweep: 8 ,Energy: -24.7630, max bond 33
sweep: 9 ,Energy: -24.7630, max bond 32
sweep: 10 ,Energy: -24.7630, max bond 37
Ground state energy per site: -1.23814899941336

The algorithm converges quickly, with the energy stabilizing after a few sweeps. Notice that the bond dimension grows adaptively according to the logarithmic schedule.

Step 7: Verify convergence

Check that DMRG converged properly:

// Get the final bond dimension
var final_bond_dim = ground_state.max_bond()
var max_allowed_bond = dmrg_options["max_bond_dimension"]

if (final_bond_dim < max_allowed_bond) {
    print("Converged with bond dimension: " + string(final_bond_dim))
    print("Result is reliable")
} else {
    print("Hit maximum bond dimension - consider increasing it")
}

// Check if MPS is in canonical form
if (ground_state.is_canonical()) {
    print("Ground state is in canonical form")
}

Complete example

Here's the full script:

// Parameters
var L = 20
var phys_dim = 2
var Jx = 1.0
var Jy = 1.0

// Create lattice and model
var lat = lattice("chain", [L])
var M = model(lat)

// Define coefficients and neighbourhood
M.add(coefficient("Jx", Jx))
M.add(coefficient("Jy", Jy))
M.add(neighbourhood_rule("nn", [1]))

// Add XY interaction terms
M.add(term("Jx", [X, X], "nn"))
M.add(term("Jy", [Y, Y], "nn"))

// Convert to MPO
var H = mpo(M)

// Initial state
var psi0 = make_random_mps(L, phys_dim, 4, as_complex)

// DMRG options
var dmrg_options = [
    "schedule": "logarithmic",
    "max_bond_dimension": 64,
    "max_iterations": 100,
    "energy_accuracy": 1e-8,
    "truncation_accuracy": 1e-10,
    "verbose": "sweep"
]

// Run DMRG
var result = dmrg(psi0, H, dmrg_options)
var ground_energy = result[0]
var ground_state = result[1]

print("Ground state energy: " + string(ground_energy))
print("Energy per site: " + string(ground_energy / L))
print("Final bond dimension: " + string(ground_state.max_bond()))

Next steps

• Try different coupling strengths (Jx, Jy) and observe how the ground state changes
• Increase system size L and monitor how computation time scales
• Experiment with different schedules (static, logarithmic, polynomial)