# 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
.. 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


[ ]: