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:
The system is the Fermi-Hubbard modelon a square lattice:
where \(j=(x,y)\) denotes the lattice site. We choose a mean-field initial state,
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¶
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)
|