Efficient time evolution using the omp-parallelized expm_multiply_parallel

In this example, we demonstrate the usage of the function tools.evolution.expm_multiply_parallel(), designed to compute matrix exponentials for static Hamiltonians.

One particular application of the matrix exponential is unitary dynamics generated by a static Hamiltonian (e.g., quantum quanches, etc.). In some cases, one can also use piecewise-constant functions to define a nontrivial dynamics: examples onclude periodically-driven systems, but also some Trotterization schemes.

The function tools.evolution.expm_multiply_parallel is a modified, omp-parallelized implementation of scipy.sparse.linalg.expm_multiply. See line 5 in the code snippet below to set the number of omp threads [omp version of QuSpin only!].

To showcase the usage of tools.evolution.expm_multiply_parallel(), consider unitary time evolution generated by the \(T\)-periodic spin-1 Heisenberg-like Hamiltonian

\[\begin{split}H(t) = \bigg\{ \!\begin{array}{c}\! &H_0,\qquad 0\leq t\leq T/2 \\ \!&H_1,\qquad T/2<t\leq T \end{array}\end{split}\]
\[H_0 = \sum_{j=0}^{L-1} \frac{1}{2} J_{xy} S^+_{j+1}S^-_j + \mathrm{h.c.}, \qquad H_1 = \sum_{j=0}^{L-1} J_{zz} S^z_{j+1}S^z_j + h_z S^z_j\]

where \(\vec S_j\) is a spin-1 operator of Pauli operators acting on lattice site \(j\). We use periodic bounary conditions and work in the zero magnetization sector.

We choose as the initial state the ground state of the average Hamiltonian \(H_\mathrm{ave} = 1/2(H_0+H_1)\), and evolve it under \(H(t)\). Every evolution cycle, we measure the energy density and the entanglement entropy density of half the system:

\[\begin{split}\mathcal{E}_\mathrm{ave}(\ell T) &= \frac{1}{L}\langle\psi(\ell T)| H_\mathrm{ave} |\psi(\ell T)\rangle, \quad |\psi(\ell T)\rangle = [\exp(-i H_1 T/2)\exp(-i H_0 T/2)]^\ell|\psi(0)\rangle, \\ s_\mathrm{ent}(\ell T) &= \frac{1}{L_A} \mathrm{tr}_A \left[ \rho_A(\ell T)\log \rho_A(\ell T)\right], \quad \rho_A(\ell T) = \mathrm{tr}_{L \backslash L_A} |\psi(\ell T)\rangle\langle\psi(\ell T)|.\end{split}\]

We compare the entanglement entropy density value to the Page value, cf. PRL 71, 1291 (1993), which takes into account finite-size corrections. As expected, the system heats up quickly to an infinite-temperature state.

Script

download script

  1#
  2import sys, os
  3
  4os.environ["KMP_DUPLICATE_LIB_OK"] = (
  5    "True"  # uncomment this line if omp error occurs on OSX for python 3
  6)
  7os.environ["OMP_NUM_THREADS"] = "4"  # set number of OpenMP threads to run in parallel
  8os.environ["MKL_NUM_THREADS"] = "1"  # set number of MKL threads to run in parallel
  9#
 10
 11###########################################################################
 12#                            example 22                                   #
 13#  This example shows the usage of the function `expm_multiply_parallel   #
 14#  to do time evolution for piece-wise constatnt Hamiltonians. For this   #
 15#  purpose, we show a simulation of a periodically-driven Heinseberg-ike  #
 16#  spin-1 system.                                                         #
 17###########################################################################
 18from quspin.basis import spin_basis_1d
 19from quspin.operators import hamiltonian
 20from quspin.tools.evolution import expm_multiply_parallel
 21import numpy as np
 22import time
 23
 24#
 25##### define data type for the simulation
 26dtype_real = np.float64
 27dtype_cmplx = np.result_type(dtype_real, np.complex64)
 28#
 29##### define model parameters #####
 30L = 12  # system size
 31Jxy = np.sqrt(2.0)  # xy interaction
 32Jzz_0 = 1.0  # zz interaction
 33hz = 1.0 / np.sqrt(3.0)  # z external field
 34T = 1.5  # period of switching for periodic drive
 35N_steps = 100  # number of driving cycles to evolve for
 36#
 37##### define Hamiltonians H_0, H_1 and H_ave
 38# build the spin-1 basis in the zero magnetization, positive parity and zero-momentum sector
 39basis = spin_basis_1d(
 40    L,
 41    S="1",
 42    m=0,
 43    kblock=0,
 44    pblock=1,
 45)
 46print("total number of basis states {}.\n".format(basis.Ns))
 47# define operators with OBC using site-coupling lists
 48J_zz = [[Jzz_0, i, (i + 1) % L] for i in range(L)]  # PBC
 49J_xy = [[0.5 * Jxy, i, (i + 1) % L] for i in range(L)]  # PBC
 50h_z = [[hz, i] for i in range(L)]
 51# static and dynamic lists
 52static_0 = [
 53    ["+-", J_xy],
 54    ["-+", J_xy],
 55]
 56static_1 = [
 57    ["zz", J_zz],
 58    ["z", h_z],
 59]
 60dynamic = []
 61# compute the time-dependent Heisenberg Hamiltonian
 62H0 = hamiltonian(static_0, dynamic, basis=basis, dtype=dtype_real)
 63H1 = hamiltonian(static_1, dynamic, basis=basis, dtype=dtype_real)
 64H_ave = 0.5 * (H0 + H1)
 65#
 66##### compute the initial state
 67# calculate ground state of H_ave
 68E, V = H_ave.eigsh(k=1, which="SA")
 69psi_i = V[:, 0]
 70#
 71# preallocate arrays for the observables
 72#
 73E_density = np.zeros(N_steps + 1, dtype=dtype_real)
 74Sent_density = np.zeros(N_steps + 1, dtype=dtype_real)
 75# compute initial values for obsrvables
 76E_density[0] = H_ave.expt_value(psi_i).real / L
 77Sent_density[0] = basis.ent_entropy(psi_i, sub_sys_A=range(L // 2), density=True)[
 78    "Sent_A"
 79]
 80#
 81##### compute the time evolution using expm_multiply_parallel
 82#
 83# construct piece-wise constant unitaries
 84expH0 = expm_multiply_parallel(H0.tocsr(), a=-1j * 0.5 * T, dtype=dtype_cmplx)
 85expH1 = expm_multiply_parallel(H1.tocsr(), a=-1j * 0.5 * T, dtype=dtype_cmplx)
 86#
 87# auxiliary array for memory efficiency
 88psi = psi_i.copy().astype(np.complex128)
 89work_array = np.zeros(
 90    (2 * len(psi),), dtype=psi.dtype
 91)  # twice as long because complex-valued
 92#
 93# loop ober the time steps
 94for j in range(N_steps):
 95    #
 96    # apply to state psi and update psi in-place
 97    expH0.dot(psi, work_array=work_array, overwrite_v=True)
 98    expH1.dot(psi, work_array=work_array, overwrite_v=True)
 99    # measure 'oservables'
100    E_density[j + 1] = H_ave.expt_value(psi).real / L
101    Sent_density[j + 1] = basis.ent_entropy(psi, sub_sys_A=range(L // 2), density=True)[
102        "Sent_A"
103    ]
104    #
105    print("finished evolving {0:d} step".format(j))
106#
107# compute Page-corrected entanglement entropy value
108m = basis.sps ** (L // 2)
109n = basis.sps**L
110s_page = (np.log(m) - m / (2.0 * n)) / (L // 2)
111#
112#
113##### Plot data
114#
115import matplotlib.pyplot as plt  # plotting library
116
117#
118times = T * np.arange(N_steps + 1)
119#
120plt.plot(times, E_density, "-b", label="$\\mathcal{E}_\\mathrm{ave}(\\ell T)$")
121plt.plot(times, Sent_density, "-r", label="$s_\\mathrm{ent}(\\ell T)$")
122plt.plot(times, s_page * np.ones_like(times), "--r", label="$s_\\mathrm{Page}$")
123plt.xlabel("$\\ell T$")
124# plt.xlim(-T,T*(N_steps+1))
125plt.legend()
126plt.grid()
127plt.tight_layout()
128#
129# plt.show()