This page was generated from exercise_2_mps.ipynb (download).

MPS and model basics

In this notebook, we introduce the SimpleMPS class from tenpy_toycodes/a_mps.py and the model class from tenpy_toycodes/b_model.py.

# standard imports and cosmetics

import numpy as np

import matplotlib.pyplot as plt

np.set_printoptions(precision=5, suppress=True, linewidth=100, threshold=50)
plt.rcParams['figure.dpi'] = 150

SimpleMPS class from tenpy_toycodes/a_mps.py

The file tenpy_toycodes/a_mps.py defines a SimpleMPS class, that provides methods for expectation values and the entanglement entropy.

You can initialize an inital product state MPS with the provided functions - init_FM_MPS to initialize the state \(\lvert \uparrow \uparrow \cdots \uparrow \uparrow \rangle\), and - init_Neel_MPS to initialize the Neel state \(\lvert \uparrow \downarrow \cdots \uparrow \downarrow \rangle\)

import tenpy_toycodes.a_mps
from tenpy_toycodes.a_mps import SimpleMPS, init_FM_MPS, init_Neel_MPS

psi_FM = init_FM_MPS(L=10)
print(psi_FM)

Z = np.diag([1., -1.])

print(psi_FM.site_expectation_value(Z))

Exercise: expectation values and entropy

• Initialize a Neel state MPS. Print the expectation values of Z

• Print the entanglement entropy. What do you expect? Why do you get so many numbers, and not just one? Tip: read the code ;-)

• Extract the half-chain entanglement entropy, i.e., the entropy when cutting the chain into two equal-length halves.

• Read the code of a_mps.py to find out how to get the correlation \(\langle \psi| Z_1 Z_6 |\psi \rangle\). Try it out!

Exercise: init_PM_MPS()

Write a function init_PM_MPS to initialize the state \(\lvert \rightarrow \rightarrow \cdots \rightarrow \rightarrow \rangle\), where \(\lvert \rightarrow \rangle = \frac{1}{\sqrt{2}} \big( \lvert\uparrow \rangle + \lvert\downarrow\rangle \big)\) is the spin-1/2 state pointing in plus x direction.

Tip: the code should be similar to init_FM_MPS.

Model class from tenpy_toycodes/b_model.py

The file tenpy_toycodes/b_model.py defines a TFIModel class representing the transverse field Ising model

\[H = - J \sum_{i} Z_i Z_{i+1} - g \sum_{i} X_i\]

It provides the Hamiltonian both in the form of bond-terms H_bonds (as required for TEBD) and in the form of an MPO H_mpo (as required for DMRG). You can use H_bonds with SimpleMPS.bond_expectation_values to evalue the energy:

from tenpy_toycodes.b_model import TFIModel

L = 10
J = 1.
g = 1.2
model = TFIModel(L=L, J=J, g=g, bc='finite')

print("<H_bonds> = ", psi_FM.bond_expectation_value(model.H_bonds))
print("energy:", np.sum(psi_FM.bond_expectation_value(model.H_bonds)))
# (make sure the model and state have the same length and boundary conditions!)

print("should be", (L-1)* (-J) * (1. * 1.) + L * (-g) * (0.) )

Exercise

• Find the code where the MPO W for this model is defined to be \(W = \begin{pmatrix} \mathbb{1} & \sigma^Z & -g \sigma^X \\ & & -J \sigma^Z \\ & & \mathbb{1} \end{pmatrix}\)

• Check the energies for the other initial states and make sure it matches what you expect.

Exercises (optional, if time left - for the experts, and those who want to become them)

• Write an optimized function correlation_function_all_j(psi, op_i, i, op_j, max_j) that returns the same values as the following, naive snippet: results = []   for j in range(i+1, max_j):   results.append(psi.correlation_function(op_i, i, op_j, j)   return results This snippet is \(\mathcal{O}(L^2)\): for each j it calls the correlation_function, which internally also has an \(\mathcal{O}(L)\) loop for the contractions. You can get this down to a single \(\mathcal{O}(L)\) loop, if you identify and reuse the parts that are the same in the diagrams for the contractions.

• For benchmarks, it’s useful to consider the XX chain in a staggered field, given by the Hamiltonian

  \[ \begin{align}\begin{aligned} H = \sum_{i=0}^{N-2} (\sigma^x_i \sigma^x_{i+1} + \sigma^y_i \sigma^y_{i+1}) - h_s \sum_{i=0}^{N-1} (-1)^i \sigma^z_i = 2 \sum_{i=0}^{N-2} (\sigma^+_i \sigma^-_{i+1} + \sigma^+_i \sigma^-_{i+1}) - h_s \sum_{i=0}^{N-1} (-1)^i \sigma^z_i\\ for the usual Pauli matrices :math:`\sigma^x, \sigma^y, \sigma^z`.\end{aligned}\end{align} \]

A Jordan-Wigner transformation maps the XX Chain to free fermions, which we can diagonalize exactly with a few lines of python codes that are given in tenpy_toycodes/free_fermions_exact.py.

from tenpy_toycodes.free_fermions_exact import XX_model_ground_state_energy

print("E_exact = ", XX_model_ground_state_energy(L=10, h_staggered=0.))

The following code implements the model for the XX Chain hamiltonian, but the MPO lacks some terms. Fill them in!

Tip: In Python, (-1)**i represents \((-1)^i\).

Tip: For the Hamiltonian

\[H = \sum_{i=0}^{N-2} (\sigma^x_i \sigma^x_{i+1} + \sigma^y_i \sigma^y_{i+1} + \sigma^z_i \sigma^z_{i+1}) - h \sum_{i=0}^{N-1} \sigma^z_i,\]

a possible MPO matrix W looks like

\[ \begin{align}\begin{aligned}\begin{split} W = \begin{pmatrix} 1 & \sigma^x & \sigma^y & \sigma^z & -h \sigma^z \\ 0 & 0 & 0 & 0 & \sigma^x \\ 0 & 0 & 0 & 0 & \sigma^y \\ 0 & 0 & 0 & 0 & \sigma^z \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} .\end{split}\\Which parts do we need here?\end{aligned}\end{align} \]

Compare the energies of different states to check that you got it correct. To really test it well, you should also check some states with non-trivial entanglement, e.g. ground states as obtained by DMRG (obtained as discussed in the next exercise notebook), for which you can directly compare to the XX_model_ground_state_energy.

class XXChain:
    """Simple class generating the Hamiltonian of the
    The Hamiltonian reads
    .. math ::
        H = - J \\sum_{i} \\sigma^x_i \\sigma^x_{i+1} - g \\sum_{i} \\sigma^z_i
    """
    def __init__(self, L, hs, bc='finite'):
        assert bc in ['finite', 'infinite']
        self.L, self.d, self.bc = L, 2, bc
        self.hs = hs
        self.sigmax = np.array([[0., 1.], [1., 0.]])  # Pauli X
        self.sigmay = np.array([[0., -1j], [1j, 0.]]) # Pauli Y
        self.sigmaz = np.array([[1., 0.], [0., -1.]]) # Pauli Z
        self.id = np.eye(2)
        self.init_H_bonds()
        self.init_H_mpo()

    def init_H_bonds(self):
        """Initialize `H_bonds` hamiltonian."""
        sx, sy, sz, id = self.sigmax, self.sigmay, self.sigmaz, self.id
        d = self.d
        nbonds = self.L - 1 if self.bc == 'finite' else self.L
        H_list = []
        for i in range(nbonds):
            hL = hR = 0.5 * self.hs
            if self.bc == 'finite':
                if i == 0:
                    hL = self.hs
                if i + 1 == self.L - 1:
                    hR = self.hs
            H_bond = np.kron(sx, sx) + np.kron(sy, sy)
            H_bond = H_bond - hL * (-1)**i * np.kron(sz, id) - hR * (-1)**(i+1) * np.kron(id, sz)
            # H_bond has legs ``i, j, i*, j*``
            H_list.append(np.reshape(H_bond, [d, d, d, d]))
        self.H_bonds = H_list

    # (note: not required for TEBD)
    def init_H_mpo(self):
        """Initialize `H_mpo` Hamiltonian."""
        w_list = []
        for i in range(self.L):
            w = np.zeros((4, 4, self.d, self.d), dtype=complex)
            w[0, 0] = w[3, 3] = self.id

            raise NotImplementedError("add further entries here")

            w_list.append(w)
        self.H_mpo = w_list

#model = XXChain(9, 4., bc='finite')