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

  \[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 \(\sigma^x, \sigma^y, \sigma^z\).

  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{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?

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')