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