7. MPS Evolution with TEBD

This example demonstrates use of the TEBD object to simulate the time evolution of matrix product states.

%config InlineBackend.figure_formats = ['svg']
import quimb as qu
import quimb.tensor as qtn
import numpy as np

To address an interesting and practical case (entanglement doesn’t grow too much) we’ll use as an initial state the all zero state apart from two flipped spins:

L = 44
zeros = '0' * ((L - 2) // 3)
binary = zeros + '1' + zeros + '1' + zeros
print('psi0:', f"|{binary}>")

psi0: |00000000000000100000000000000100000000000000>

We can turn this into a MatrixProductState using MPS_computational_state():

psi0 = qtn.MPS_computational_state(binary)
psi0.show()  # prints ascii representation of state

     1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1    
... >─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─>─ ...
    │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │     
                                 ...                                  
     1 1 1 1 1 1 1 1 
    >─>─>─>─>─>─>─>─<
    │ │ │ │ │ │ │ │ │

To set the hamiltonian, we need to define a nearest neighbour interaction (LocalHam1D) hamiltonian. These can built directly using arrays representing the one and two site terms, or from a SpinHam1D. But, for the Heisenberg hamiltonian with no magnetic field,

\[ \hat{H} = J \sum_{i=1}^{L - 1} \mathbf{\sigma}_i \cdot \mathbf{\sigma}_{i + 1} ~ , \]

we’ll use a built-in:

H = qtn.ham_1d_heis(L)

Now we are ready to set-up the actual TEBD object, which trys to mimic the api of the normal Evolution object.

We could set a default timestep dt here, or instead, a default tolerance tol, which calculates dt automatically based on the total evolution length, trotter decomposition order, and hamiltonian norm.

tebd = qtn.TEBD(psi0, H)

# Since entanglement will not grow too much, we can set quite
#     a small cutoff for splitting after each gate application
tebd.split_opts['cutoff'] = 1e-12

We’ll also set up some parameters and quantities to compute during the evolution:

# times we are interested in
ts = np.linspace(0, 80, 101)

mz_t_j = []  # z-magnetization
be_t_b = []  # block entropy
sg_t_b = []  # schmidt gap

# range of bonds, and sites
js = np.arange(0, L)
bs = np.arange(1, L)

Now we are ready to being the evolution, we use the generator at_times to yield the state at each target time, and set a tol here which will calculate a timestep to use. At each time, we’ll compute the desired quantities and add them to our results.

# generate the state at each time in ts
#     and target error 1e-3 for whole evolution
for psit in tebd.at_times(ts, tol=1e-3):
    mz_j = []
    be_b = []
    sg_b = []
    
    # there is one more site than bond, so start with mag
    #     this also sets the orthog center to 0
    mz_j += [psit.magnetization(0)]
    
    for j in range(1, L):
        # after which we only need to move it from previous site
        mz_j += [psit.magnetization(j, cur_orthog=j - 1)]
        be_b += [psit.entropy(j, cur_orthog=j)]
        sg_b += [psit.schmidt_gap(j, cur_orthog=j)]
        
    mz_t_j += [mz_j]
    be_t_b += [be_b]
    sg_t_b += [sg_b]

t=80, max-bond=15: 100%|###################################################################################################################| 101/101 [02:05<00:00,  1.24s/it]

tebd.pt.show()

     2 4 5 6 7 8 9 10 11 12 13 14 14 14 15 14 14 15 14 14 14 15 15 15 14 1    
... >─>─>─>─>─>─>─>──>──>──>──>──>──>──>──>──>──>──>──>──>──>──>──>──>──>─ ...
    │ │ │ │ │ │ │ │  │  │  │  │  │  │  │  │  │  │  │  │  │  │  │  │  │  │     
                                 ...                                  
    5 14 14 14 15 14 14 13 12 11 10 9 8 7 6 5 4 2 
    ─>──>──>──>──>──>──>──>──>──>──>─>─>─>─>─>─>─●
     │  │  │  │  │  │  │  │  │  │  │ │ │ │ │ │ │ │

TEBD.err contains a (rough) upper estimate for the error incurred by the trotter decomposition so far:

tebd.err  #  should be < tol=1e-3

0.0009938931328611381

We can also check the normalization of final state:

tebd.pt.H @ tebd.pt

(0.9999999994906472+5.981119624341705e-17j)

And that energy has been conserved:

H = qtn.MPO_ham_heis(L)
print("Initial energy:", qtn.expec_TN_1D(psi0.H, H, psi0))
print("Final energy:", qtn.expec_TN_1D(tebd.pt.H , H, tebd.pt))

Initial energy: 8.75
Final energy: (8.749999804810404-8.169106342572505e-16j)

Finally let’s plot the quantities we calculated:

import matplotlib.pyplot as plt

with plt.style.context(qu.NEUTRAL_STYLE):
    plt.figure(figsize=(12, 7))

    # plot the magnetization
    ax1 = plt.subplot(131)
    plt.pcolormesh(js, ts, np.real(mz_t_j), vmin=-0.5, vmax=0.5)
    plt.set_cmap('RdYlBu')
    plt.colorbar()
    plt.title('Z-Magnetization')
    plt.xlabel('Site')
    plt.ylabel('time [ $Jt$ ]')

    # plot the entropy
    ax2 = plt.subplot(132, sharey=ax1)
    plt.pcolormesh(bs, ts, be_t_b)
    plt.setp(ax2.get_yticklabels(), visible=False)
    plt.set_cmap('viridis'), plt.colorbar()
    plt.title('Block Entropy')
    plt.xlabel('Bond')

    # plot the schmidt gap
    ax3 = plt.subplot(133, sharey=ax1)
    plt.pcolormesh(bs, ts, sg_t_b, vmin=0, vmax=1)
    plt.setp(ax3.get_yticklabels(), visible=False)
    plt.set_cmap('magma_r')
    plt.colorbar()
    plt.title('Schmidt Gap')
    plt.xlabel('Bond')

    plt.show()

/tmp/ipykernel_2042863/1525416688.py:6: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(js, ts, np.real(mz_t_j), vmin=-0.5, vmax=0.5)
/tmp/ipykernel_2042863/1525416688.py:8: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.colorbar()
/tmp/ipykernel_2042863/1525416688.py:15: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(bs, ts, be_t_b)
/tmp/ipykernel_2042863/1525416688.py:17: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.set_cmap('viridis'), plt.colorbar()
/tmp/ipykernel_2042863/1525416688.py:23: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(bs, ts, sg_t_b, vmin=0, vmax=1)
/tmp/ipykernel_2042863/1525416688.py:26: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.colorbar()

Where we can see ballistic propagation (and reflection) of the ‘light cone’.

7.1. Non-translationally invariant Hamiltonians

The LocalHam1D class can also handle Hamiltonians with site specific interactions and fields. A good example of this is the Many-Body-Localized spin hamiltonian:

\[ \hat{H} = J \sum_{i=1}^{L - 1} \mathbf{\sigma}_i \cdot \mathbf{\sigma}_{i + 1} + \sum_i^L h_i \sigma^Z_i \]

Where \(h_i\) is a random variable. Here we construct it manually:

builder = qtn.SpinHam1D(S=1/2)

# specify the interaction term (defaults to all sites)
builder += 0.5, '+', '-'
builder += 0.5, '-', '+'
builder += 1.0, 'Z', 'Z'

# add random z-fields to each site
np.random.seed(2)
for i in range(L):
    builder[i] += 2 * np.random.rand() - 1, 'Z'
    
H = builder.build_local_ham(L)

But there is also the built-in ham1d_mbl() which has options for various random noise types etc.

We can perform the same TEBD evolution as above with this new Hamiltonian:

tebd = qtn.TEBD(psi0, H)
tebd.split_opts['cutoff'] = 1e-10

# times we are interested in
ts = np.linspace(0, 80, 101)

mz_t_j = []  # z-magnetization
be_t_b = []  # block entropy
sg_t_b = []  # schmidt gap

# range of bonds, and sites
js = np.arange(0, L)
bs = np.arange(1, L)

# generate the state at each time in ts
#     and target error 1e-3 for whole evolution
for psit in tebd.at_times(ts, tol=1e-3):
    mz_j = []
    be_b = []
    sg_b = []
    
    # there is one more site than bond, so start with mag
    #     this also sets the orthog center to 0
    mz_j += [psit.magnetization(0)]
    
    for j in range(1, L):
        # after which we only need to move it from previous site
        mz_j += [psit.magnetization(j, cur_orthog=j - 1)]
        be_b += [psit.entropy(j, cur_orthog=j)]
        sg_b += [psit.schmidt_gap(j, cur_orthog=j)]
        
    mz_t_j += [mz_j]
    be_t_b += [be_b]
    sg_t_b += [sg_b]

t=80, max-bond=19: 100%|###################################################################################################################| 101/101 [02:21<00:00,  1.40s/it]

And finally, plot the quantities again:

with plt.style.context(qu.NEUTRAL_STYLE):
    plt.figure(figsize=(12, 7))

    # plot the magnetization
    ax1 = plt.subplot(131)
    plt.pcolormesh(js, ts, np.real(mz_t_j), vmin=-0.5, vmax=0.5)
    plt.set_cmap('RdYlBu')
    plt.colorbar()
    plt.title('Z-Magnetization')
    plt.xlabel('Site')
    plt.ylabel('time [ $Jt$ ]')

    # plot the entropy
    ax2 = plt.subplot(132, sharey=ax1)
    plt.pcolormesh(bs, ts, be_t_b)
    plt.setp(ax2.get_yticklabels(), visible=False)
    plt.set_cmap('viridis'), plt.colorbar()
    plt.title('Block Entropy')
    plt.xlabel('Bond')

    # plot the schmidt gap
    ax3 = plt.subplot(133, sharey=ax1)
    plt.pcolormesh(bs, ts, sg_t_b, vmin=0, vmax=1)
    plt.setp(ax3.get_yticklabels(), visible=False)
    plt.set_cmap('magma_r')
    plt.colorbar()
    plt.title('Schmidt Gap')
    plt.xlabel('Bond')

    plt.show()

/tmp/ipykernel_2042863/1525416688.py:6: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(js, ts, np.real(mz_t_j), vmin=-0.5, vmax=0.5)
/tmp/ipykernel_2042863/1525416688.py:8: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.colorbar()
/tmp/ipykernel_2042863/1525416688.py:15: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(bs, ts, be_t_b)
/tmp/ipykernel_2042863/1525416688.py:17: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.set_cmap('viridis'), plt.colorbar()
/tmp/ipykernel_2042863/1525416688.py:23: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.pcolormesh(bs, ts, sg_t_b, vmin=0, vmax=1)
/tmp/ipykernel_2042863/1525416688.py:26: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  plt.colorbar()

Here we can see much more confinement.