Liouville-von Neumann Equation using the MKL-enhanced Sparse Matrix Product [courtesy of J. Verlage]

This example shows how one can combine QuSpin with the external library sparse_dot which supports an MKL-parallelized sparse matrix product.

To this end, we consider the numerical solution of the Liouville-von Neumann (LvN) equation for the density matrix:

\[i\partial_t \rho(t) = [H,\rho(t)].\]

The system is the Fermi-Hubbard modelon a square lattice:

\[H = -J\sum_{j,\sigma} \left( c^\dagger_{j+1\sigma}c_{j\sigma} + \mathrm{h.c.} \right) + U\sum_j n_{j\uparrow}n_{j\downarrow},\]

where \(j=(x,y)\) denotes the lattice site. We choose a mean-field initial state,

\[\rho(0)=\bigotimes_j \rho_j, \qquad \mathrm{where} \qquad \rho_j = \frac{1}{2}\left( |\uparrow_j\rangle\langle \uparrow_j|+ |\downarrow_j\rangle\langle \downarrow_j| \right),\]

that cannot be written as a pure state [hence the necessity to solve the LvN equation rather than Schroedinger’s equation].

Note that the initial state \(\rho(0)\) is diagonal in the particle number basis; therefore, since the Hamiltonian \(H\) is also sparse, we expect that the time-evolved density operator will remain sparse at least for small times [compared to \(U^{-1}, J^{-1}\)]. Since we are limited to small system sizes by the exponentially growing Hilbert space dimension, we need a memory-efficient way to store the quantum state, e.g., using a sparse matrix. In turn, this requires:

  • an efficient, ideally parallelized, sparse-spase matrix product;

  • a solver for differential equations that allows us to keep the variable [here \(\rho\)] in sparse format at all times.

To this end, we can use the open-source python library sparse_dot, which provides the MKL-paralellized function dot_product_mkl. We use it to write our own fourth-order Runge-Kutta (RK) solver for the LvN equation. Note that, unlike the RK solver provided in Scipy where the step size is chosen adaptively, our RK implementation has a fixed step size; however, scipy’s solver does not allow us to keep the state as a sparse matrix at all times.

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"] = "4"  # set number of MKL threads to run in parallel
  9#
 10
 11###########################################################################
 12#                             example 27                                  #
 13#  In this script we demonstrate how to use QuSpin to generate            #
 14#  a Hamiltonian and solve the Louiville-von Neumann equation starting	  #
 15#  from a mixed initial state in the Fermi Hubbard model. We also         #
 16#  show how to write a simple fixed time-step Runge-Kutta solver          #
 17#  that makes use of an MKL-parllelized dot function for sparse matrices. #
 18###########################################################################
 19# import sparse_dot library, see https://github.com/flatironinstitute/sparse_dot.git
 20from sparse_dot_mkl import dot_product_mkl
 21from quspin.tools.misc import get_matvec_function
 22from quspin.operators import hamiltonian
 23from scipy.sparse import csr_matrix
 24from quspin.basis import spinful_fermion_basis_general
 25import numpy as np
 26import time
 27
 28#
 29##### define model parameters #####
 30#
 31Lx, Ly = 2, 2  # expect to see an MKL speedup from Lx,Ly = 2,3 onward
 32N_2d = Lx * Ly  # total number of lattice sites
 33# model params
 34J = 1.0  # hopping amplitude
 35U = np.sqrt(2.0)  # interaction strength
 36# time parameters
 37t_max = 40.0  # total time
 38dt = 0.1  # time step size
 39N_T = int(t_max // dt)
 40time_vec = np.linspace(0.0, t_max, int(2 * t_max + 1))
 41#
 42##### create Hamiltonian to evolve unitarily
 43#
 44# basis
 45basis = spinful_fermion_basis_general(N_2d)
 46# define translation operators for 2D lattice
 47s = np.arange(N_2d)  # sites [0,1,2,...,N_2d-1] in simple notation
 48x = s % Lx  # x positions for sites
 49y = s // Lx  # y positions for sites
 50T_x = (x + 1) % Lx + Lx * y  # translation along x-direction
 51T_y = x + Lx * ((y + 1) % Ly)  # translation along y-direction
 52# site-coupling lists
 53hop_left = [[-J, i, T_x[i]] for i in range(N_2d)] + [
 54    [-J, i, T_y[i]] for i in range(N_2d)
 55]
 56hop_right = [[+J, i, T_x[i]] for i in range(N_2d)] + [
 57    [+J, i, T_y[i]] for i in range(N_2d)
 58]
 59int_list = [[U, i, i] for i in range(N_2d)]
 60# static opstr list
 61static = [
 62    ["+-|", hop_left],  # up hop left
 63    ["-+|", hop_right],  # up hop right
 64    ["|+-", hop_left],  # down hop left
 65    ["|-+", hop_right],  # down hop right
 66    ["n|n", int_list],
 67]
 68# construct Hamiltonian
 69Hcsc = hamiltonian(
 70    static, [], dtype=np.complex128, basis=basis, check_symm=False
 71).tocsr()
 72#
 73##### create the mean-field groundstate we start from
 74#
 75# compute basis with single occupancies only
 76basis_reduced = spinful_fermion_basis_general(
 77    N_2d, Nf=([(j, N_2d - j) for j in range(N_2d + 1)]), double_occupancy=False
 78)
 79# create empty list to store indices of nonzero elements for initial DM
 80rho_inds = []
 81for s in basis_reduced.states:  # loop over singly-occupied states
 82    rho_inds.append(np.argwhere(basis.states == s)[0][0])
 83# create initial state in csr format
 84rho_0 = csr_matrix(
 85    (np.ones(basis_reduced.Ns) / basis_reduced.Ns, (rho_inds, rho_inds)),
 86    shape=(basis.Ns, basis.Ns),
 87    dtype=np.complex128,
 88)
 89
 90
 91#
 92##### define Runge-Kutta solver for sparse matrix
 93#
 94# MKL-parallel function using the sparse_dot library
 95def LvN_mkl(rho):
 96    # define right-hand side of Liouville-von Neumann equation
 97    # see https://github.com/flatironinstitute/sparse_dot.git, needs v0.8 or higher
 98    return -1j * (
 99        dot_product_mkl(Hcsc, rho, cast=False) - dot_product_mkl(rho, Hcsc, cast=False)
100    )
101
102
103#
104# scipy function
105def LvN_scipy(rho):
106    return -1j * (Hcsc @ rho - rho @ Hcsc)
107
108
109#
110# define fixed step-size Runge-Kutta 4th order method
111def RK_solver(rho, dt, LvN):
112    k1 = LvN(rho)
113    k2 = LvN(rho + (0.5 * dt) * k1)
114    k3 = LvN(rho + (0.5 * dt) * k2)
115    k4 = LvN(rho + dt * k3)
116    return rho + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0
117
118
119#
120##### evolve DM by solving the LvN equation
121#
122# empty list to store the solution in
123rho_t = []
124# initial state
125rho_mkl = rho_0.copy()
126# time evolution loop
127starttime = time.time()
128for i in range(N_T):
129    rho_mkl = RK_solver(rho_mkl, dt, LvN_mkl)
130    rho_t.append(rho_mkl)
131    # print("finished step {0:d}/{1:d}.".format(i+1,int(t_max/dt)-1),flush=True)
132#
133print(
134    "\nMKL time evo done in {0:0.4f} secs.".format(time.time() - starttime), flush=True
135)
136#
137# empty list to store the solution in
138rho_t = []
139# initial state
140rho_scipy = rho_0.copy()
141# time evolution loop
142starttime = time.time()
143for i in range(N_T):
144    rho_scipy = RK_solver(rho_scipy, dt, LvN_scipy)
145    rho_t.append(rho_scipy)
146    # print("finished step {0:d}/{1:d}.".format(i+1,int(t_max/dt)-1),flush=True)
147#
148print(
149    "\nScipy time evo done in {0:0.4f} secs.".format(time.time() - starttime),
150    flush=True,
151)