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

DMRG runs

In this notebook, we use the SimpleDMRGEngine class from tenpy_toycodes/d_dmrg.py to run DMRG and find MPS ground states.

# 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

In previous notebooks, we learned how to initialize SimpleMPS…

from tenpy_toycodes.a_mps import SimpleMPS, init_FM_MPS, init_Neel_MPS

L = 12

psi_FM = init_FM_MPS(L=L, d=2, bc='finite')
print(psi_FM)
SigmaZ = np.diag([1., -1.])
print(psi_FM.site_expectation_value(SigmaZ))

def init_PM_MPS(L, bc='finite'):
    """Return a paramagnetic MPS (= product state with all spins pointing in +x direction)"""
    d = 2
    B = np.zeros([1, d, 1], dtype=float)
    B[0, 0, 0] = 1./np.sqrt(2)
    B[0, 1, 0] = 1./np.sqrt(2)
    S = np.ones([1], dtype=float)
    Bs = [B.copy() for i in range(L)]
    Ss = [S.copy() for i in range(L)]
    return SimpleMPS(Bs, Ss, bc=bc)

… and how to initialize the TFIModel.

from tenpy_toycodes.b_model import TFIModel

g = 1.2
model = TFIModel(L=L, J=1., 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!)

For small enough system size \(L \lesssim 16\), you can compare the energies to exact diagonalization in the full Hilbert space (which is exponentially expensive!):

from tenpy_toycodes.tfi_exact import finite_gs_energy

if L <= 16:
    energy_exact = finite_gs_energy(L=L, J=1., g=g)

The DMRG algorithm

The file tenpy_toycodes/d_dmrg.py implements the DMRG algorithm. It can be called like this:

from tenpy_toycodes.d_dmrg import SimpleDMRGEngine, SimpleHeff2

chi_max = 15

psi = init_FM_MPS(model.L, model.d, model.bc)
eng = SimpleDMRGEngine(psi, model, chi_max=chi_max, eps=1.e-10)
for i in range(10):
    E_dmrg = eng.sweep()
    E = np.sum(psi.bond_expectation_value(model.H_bonds))
    print("sweep {i:2d}: E = {E:.13f}".format(i=i + 1, E=E))
print("final bond dimensions: ", psi.get_chi())
mag_x = np.mean(psi.site_expectation_value(model.sigmax))
mag_z = np.mean(psi.site_expectation_value(model.sigmaz))
print("magnetization in X = {mag_x:.5f}".format(mag_x=mag_x))
print("magnetization in Z = {mag_z:.5f}".format(mag_z=mag_z))
if model.L <= 16:
    E_exact = finite_gs_energy(L=model.L, J=model.J, g=model.g)
    print("err in energy = {err:.3e}".format(err=E - E_exact))

Exercise: read d_dmrg.py

Read the code of tenpy_toycodes/d_dmrg.py and try to undertstand the general structure of how it works.

Exercise: measure correlation functions

Just looking at expectation values of local operators is not enough. Also measure correlation functions \(\langle Z_{L/4} Z_{3L/4} \rangle\).

Exercise: DMRG runs

Try running DMRG for various different parameters:

• Change the bond dimension chi and truncation threashold eps.

• Change the system size L

• Change the model parameter g (at fixed \(J=1\)) to both the ferromagnetic phase \(g<J\) and the paramagnetic phase \(g>1\).

• Change the initial state.

Exercise: Phase diagram

To map out the phase diagram, it can be convenient to define a function that just runs DMRG for a given model. Fill in the below template. Use it obtain and plot the energy, correlations and magnetizations for different \(L=16, 32, 64\) against \(g\).

def run_DMRG(model, chi_max=50):
    print(f"runnning DMRG for L={model.L:d}, g={model.g:.2f}, bc={model.bc}, chi_max={chi_max:d}")

    raise NotImplementedError("TODO: this is an Exercise!")

    return psi

results_all = {}
# TIP: comment this out after the first run to avoid overriding your data!

L = 8

mag_X = []
mag_Z = []
E = []
corr = []
max_chi = []

gs = [0.1, 0.5, 0.8, 0.9, 1.0, 1.1, 1.2, 1.5]
for g in gs:
    model = TFIModel(L=L, J=1., g=g, bc='finite')
    psi = run_DMRG(model)
    mag_X.append(np.mean(psi.site_expectation_value(model.sigmax)))
    mag_Z.append(np.mean(psi.site_expectation_value(model.sigmaz)))
    E.append(np.sum(psi.bond_expectation_value(model.H_bonds)))
    corr.append(psi.correlation_function(model.sigmaz, model.L//4, model.sigmaz, model.L * 3 // 4))
    max_chi.append(max(psi.get_chi()))

results_all[L] = {
    'g': gs,
    'mag_X': mag_X,
    'mag_Z': mag_Z,
    'E': E,
    'corr': corr,
    'max_chi': max_chi
}

key = 'corr'

plt.figure()
for L in results_all:
    res_L = results_all[L]
    plt.plot(res_L['g'], res_L[key], marker='o', label="L={L:d}".format(L=L))
plt.xlabel('g')
plt.ylabel(key)
plt.legend(loc='best')

Infinite DMRG

The given DMRG code also works with bc='infinite' boundary conditions of the model and state. The given SimpleDMRG code also allows to run infinite DMRG, simply by replacing the bc='finite' for both the model and the MPS.

• Look at the implementation of d_dmrg.py (and a_mps.py) to see where the differences are.

Again, we can compare to analytic calculations possible for the Transverse Field Ising model.

from tenpy_toycodes.tfi_exact import infinite_gs_energy

The L parameter now just indices the number of tensors insite the unit cell of the infinite MPS. It has to be at least 2, since we optimize 2 tensors at once in our DMRG code. Note that we now use the mean to calculate densities of observables instead of extensive quantities:

model = TFIModel(L=2, J=1., g=0.8, bc='infinite')  # just change bc='infinite' here

chi_max = 10

psi = init_FM_MPS(model.L, model.d, model.bc)
eng = SimpleDMRGEngine(psi, model, chi_max=chi_max, eps=1.e-7)
for i in range(10):
    E_dmrg = eng.sweep()
    E = np.mean(psi.bond_expectation_value(model.H_bonds))
    #if i % 10 == 9:
    print("sweep {i:2d}: E/L = {E:.13f}".format(i=i + 1, E=E))
print("final bond dimensions: ", psi.get_chi())
mag_x = np.mean(psi.site_expectation_value(model.sigmax))
mag_z = np.mean(psi.site_expectation_value(model.sigmaz))
print("magnetization density in X = {mag_x:.5f}".format(mag_x=mag_x))
print("magnetization density in Z = {mag_z:.5f}".format(mag_z=mag_z))
E_exact = infinite_gs_energy(model.J, model.g)
print("err in energy = {err:.3e}".format(err=E - E_exact))

Exercise: Infinite DMRG

Try running the infinite DMRG code. How many sweeps do you now need to perform to converge now? Does it depend on g and chi?

From the exercise, you should see that you need significantly more sweeps when the correlation length gets larger, in particular at the critical point! At the critical point, the correlation length (of infinite MPS) scales as \(\xi \propto \chi^\kappa\), so you again need more sweeps to converge at larger bond dimensions.

Advanced exercises for the experts (and those who want to become them ;-) )

• Obtain the ground state of the transverse field ising model at the critical point with DMRG for large L. Try to plot the corrlation function as a function of j-i. What form does it have? Is an MPS a good ansatz for that?

• Compare running DMRG and imaginary time evolution with TEBD from tenpy_toycodes/c_tebd.py for various parameters of L, J, g, and bc. Which one is faster? Do they always agree?