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
¶
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 eachj
it calls thecorrelation_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
a possible MPO matrix W looks like
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')
[ ]: