Hamiltonian simulation

Requirement: `isqdeployer`

isqdeployer is an all-in-one Python package. It serves as an Object-Oriented Interface for isq. You can easily install it using pip without the need for any additional installations or configurations.

pip install isqdeployer

Spin model

Wavefunction evolving with time

For a quantum system with Hamiltonian $H$ , its wavefunction is evolving as $|\psi(t)\rangle = e^{-i\int dt H(t)} |\psi(0)\rangle$ , where we have assume that $\hbar=1$ . If we further have $H\ne H(t)$ , then we have $|\psi(t)\rangle = e^{-i H t} |\psi(0)\rangle = \sum_ie^{-iE_i t}\langle i|\psi(0)\rangle|i\rangle$ , where we have used the spectral decomposition of $H=\sum_i E_i |i\rangle\langle i|$ . This usually needs us to obtain a diagonalization of $H$ . Practically, a more efficient way [1] to simulate the time evolution is used if we have an efficient way to implement operators directly. Here we recall its major ideal.

For any Hamiltonian in the form of $H=\sum_k h_k$ , we have $e^{-iHt}=( e^{-iH\tau} )^N$ , where $\tau=t/N$ for some large $N$ and therefore $\tau\to0$ . Then we have the following approximation:
$e^{-i\sum_k h_k \tau}\approx \prod_k e^{-i h_k \tau}.$
So the basic element for simulating the whole Hamiltonian is to simulate each term $e^{-i h_k \tau}$ .

Basic gates in spin model

In spin model, the Hamiltonian can be represented as $H= \sum _{ij;mn}\lambda^{mn}_{ij}\sigma_i^m \sigma_j^n,$ where ${\{\sigma_i^m}\}=\{I,X,Y,Z\}$ at each spin site $m$ . Correspondingly the basic gate we need to implement in the circuit is in the form of $e^{i \theta \sigma_i \sigma_j}$ for the pair of any two sites.

$e^{i a \sigma}$ gate implementation

We first consider the term $e^{i a Z}$ on one site. The corresponding implementation is given by

Rz(2a)

Similarly, we can directly implement $e^{i a X}$ and $e^{i a Y}$ by Rx and Ry. These terms represent the effect of the magnetic field in each direction.

$e^{i a \sigma_i^m \sigma_j^n}$ gate implementation

This is a two-spin interaction with action $\sigma_i$ on one site $m$ and $\sigma_j$ on the other site $n$ . We first consider the term $e^{i a ZZ}$ . This term represents an Ising type of interaction between two spins. Its implementation is given by:

iaZZ

For this gate, we can calculate the output state of each input state in $\{ |00\rangle,|01\rangle,|10\rangle,|11\rangle\}$ and see that it is equivalent to $e^{iaZZ}$ . Based on this gate, we can easily implement other types. For implementing $e^{iaXZ}$ , for instance, we can use the transformation $XZ =(U^+ZU)Z=U^+ZZU$ , where $U=Ry(\frac{\pi}{2})$ acting on the first qubit (spin). Its implementation is shown by

iaXZ

By the same spirit, we can implement any $\sigma_i^m\sigma_j^n$ terms.

Example: Single spin in the magnetic field

Consider a system with Hamiltonian $H=X$ , for input state $|\psi\rangle=|0\rangle$ , the evolution of the wavefunction is give by

$|\psi(t)\rangle= e^{-it}\frac{1}{\sqrt{2}}|+\rangle +e^{it}\frac{1}{\sqrt{2}}|-\rangle = \cos t |0\rangle - i \sin t |0\rangle .$

We then simulate this process by quantum circuit.

from isqdeployer import Circuit 
from isqdeployer.circuitDeployer import PauliHamDeployer 
from isqdeployer.utils.pauliHamiltonian import PauliHamiltonian
import numpy as np 
import matplotlib.pyplot as plt 

ham = PauliHamiltonian(nq=1) # define a Hamiltonian
ham.setOneTerm(xi=1.0,p=[1]) # Hamiltonian = Pauli X, (0,1,2,3) represents (I,X,Y,Z)

cir = Circuit(
    nq=1, # only 1 qubit
    workDir="/home/user/Desktop/test", # see how it works by inspecting internal process 
    isInputArg=True, # use additional parameter (t) to increase the speed of the calculation
    ) # circuit, init state is |0>

t = cir.getInputArg(FI='F', id=0) # set the circuit to contains a parameter t

dep = PauliHamDeployer(circuit=cir) # use a deployer to implement gates
dep.expHt(
    h=ham,
    t=t,
    N=15, # Suzuki decomposition
    ) 
cir.setMeasurement([0])

Tlist = np.arange(0,3,0.1)# the time range we study
results = cir.runJob(paramList=[{'F':[t]} for t in Tlist ]) # batch submit all jobs

#---- plot 
prob0 = [res.get('0',0) for res in results]
prob1 = [res.get('1',0) for res in results]
plt.plot(Tlist,prob0,label="Prob(0)")
plt.plot(Tlist,prob1,label="Prob(1)")
plt.xlabel("Time (t)")
plt.ylabel('Probability')
plt.legend()
plt.show()

In the example above, we have set workDir="/home/user/Desktop/test". You can observe intermediate resources in that directory. The original isq file is named resource.isq. The results are calculated using an internal default simulator and are displayed below:

fig1

Example: Spin chain model

We now study a less trivial model whose Hamiltonian is given by
$H = h\sum_{n}Z_n+\sum_{n}[J_{x}^{(n)}X_nX_{n+1}+J_{y}^{(n)}Y_nY_{n+1}+J_{z}^{(n)}Z_nZ_{n+1}].$
The first term represents a uniform magnetic field on the system. The second term represents complex XYZ interactions. This model has very important applications in condensed matter physics.

In this demo, we set up a 3-site system. Parameters are set as $h=-1$ , $J_x=-0.2$ , $J_y=-0.3$ , and $J_z=-0.1$ . Here we do not consider the connection between the first and third sites. One interested in periodic boundary conditions can set additional interaction. In addition, in the beginning, we set the first spin in the state $|1\rangle$ while others are in $|0\rangle$ .

from isqdeployer import Circuit 
from isqdeployer.circuitDeployer import PauliHamDeployer 
from isqdeployer.utils.pauliHamiltonian import PauliHamiltonian
import numpy as np 
import matplotlib.pyplot as plt 

# --- define the Hamiltonian ---
ham = PauliHamiltonian(nq=3) # define a Hamiltonian
# --- hz term
ham.setOneTerm(xi=-1.,p=[3,0,0])
ham.setOneTerm(xi=-1.,p=[0,3,0])
ham.setOneTerm(xi=-1.,p=[0,0,3])
# ---jx term 
ham.setOneTerm(xi=-0.2,p=[1,1,0])
ham.setOneTerm(xi=-0.2,p=[0,1,1])
# ---jy term 
ham.setOneTerm(xi=-0.3,p=[2,2,0])
ham.setOneTerm(xi=-0.3,p=[0,2,2])
# ---jz term 
ham.setOneTerm(xi=-0.1,p=[3,3,0])
ham.setOneTerm(xi=-0.1,p=[0,3,3])


cir = Circuit(nq=3,isInputArg=True,)  

t = cir.getInputArg(FI='F', id=0) # set the circuit to contains a parameter t

cir.X(0) # set the first spin up

dep = PauliHamDeployer(circuit=cir) # use a deployer to implement gates
dep.expHt( 
    h=ham, 
    t=t, 
    N=50, # Suzuki decomposition
    ) 
cir.setMeasurement([0,1,2])

Tlist = np.arange(0,10,0.1)# the time range we study
results = cir.runJob(paramList=[{'F':[t]} for t in Tlist ]) # batch submit all jobs


#---- plot 
def countProb(i,res):
    r = 0.0
    for k, v in res.items():r += (k[i] == "0")*v
    return r*2-1

Z0 = [ countProb(0,res) for res in results]
Z1 = [ countProb(1,res) for res in results]
Z2 = [ countProb(2,res) for res in results]

plt.plot(Tlist,Z0)
plt.plot(Tlist,Z1)
plt.plot(Tlist,Z2)

The result of the time evolution of each spin state is shown as follows:

As a comparison, an exact calculation by QuTip is shown in the figure. The precision can further increase by increasing $N$ . To be more clear, we check $\langle\sigma_z^1(t)\rangle$ more carefully by the following figure.

As shown in the figure, for $N=2$ case, in a small time regime the result of the calculation matches the exact result very well. However, at large t regime, to obtain a better result we need a larger $N$ .

Any other spin mode

One can easily extend the code to make a circuit to simulate an arbitrary spin model with more complicated interaction types. For any (2-local) spin system, we can represent it as $H=\sum_k \xi_k \sigma_{i(k)}^{m(k)} \sigma_{j(k)}^{n(k)}$ . An it can be calcualted by isqdeployer.

Fermion-Hubbard model

The Hubbard model is one of the most important systems in condensed matter physics. Studying this model will reveal much exotic physics in strongly correlated electron systems such as superconductivity and quantum Hall effect. In this tutorial, we study the time evolution of a wave function in the system with such a Hamiltonian.

Hamiltonian

We start from a most general Hamiltonian given by $H=H_0+H_I$ where $H_0=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}$ and
$H_I=\sum_{ijkl}u_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l}.$ For a system with spins, the indices $i$ and $j$ should also be implicit. $H_0$ represents the single particle terms, such as onsite energy and hopping terms. $H_I$ contains many-body terms. For example, if we set $i=k$ and $j=l$ , we obtain the Coulomb interaction $H_U = U\sum_i n_{i\uparrow} n_{i\downarrow}.$ In the most general case, $H_I$ can consider other interactions, such as Hund's interaction, and is written as $H_I = H_U + H_J + ...$ In the standard Hubbard, we consider $H_0$ and $H_U$ . Nevertheless, the following discussion is also suited for other two-particle interactions.

Our task is to find an encoding method to transform this Hamiltonian into one with Pauli tensors only. By Jordan–Wigner transformation $c_{i}=\sigma_{i}^{-}\prod_{\alpha<i}\left(-Z_{\alpha}\right)$ , the diagonal term of $H_0$ is given by $t_{ii}c^\dagger_i c_i = \frac{t_{ii}}{2}\left(1+Z_{i}\right),$ the hopping term (off-diagonal term) is given by

$\begin{aligned} t_{ij}c^\dagger_{i}c_j+t_{ji}c^\dagger_jc_i= -\chi_{ij}X_{2-}^{i}Y_{2+}^{j}\mathcal{Z}_{ij}X_{2+}^{i}Y_{2-}^{j} -\xi_{ij}X_{2-}^{i}X_{2-}^{j}\mathcal{Z}_{ij}X_{2+}^{i}X_{2+}^{j}\\ -\xi_{ij}Y_{2+}^{i}Y_{2+}^{j}\mathcal{Z}_{ij}Y_{2-}^{i}Y_{2-}^{j} +\chi_{ij}Y_{2+}^{i}X_{2-}^{j}\mathcal{Z}_{ij}Y_{2-}^{i}X_{2+}^{j},\\ \end{aligned}$
where $\chi_{ij}= i\frac{t_{ij}-t_{ji}}{2}\frac{\left(-1\right)^{j-i}}{2}$ , $\xi_{ij}=\frac{t_{ij}+t_{ji}}{2}\frac{\left(-1\right)^{j-i}}{2}$ , and $\mathcal{Z}_{ij}= Z_i\otimes Z_{i+1}...\otimes Z_{j},$ and $X_{2\pm} = R_x (\pm\frac{\pi}{2} )$ , and $Y_{2\pm}=R_y(\pm\frac{\pi}{2})$ . For real $t_{ij}$ case, we have $\chi_{ij}=0$ . Finally the Coulomb interaction is given by $H_U=\frac{U}{4}\left(1+Z_{j}+Z_{i}+Z_{i}Z_{j}\right).$

To streamline configuration, we employ a specialized package, qsci, for deploying circuits in scientific calculations. One can install it by:

pip install qsci

With qsci, it's straightforward to deploy various quantum algorithms such as VQE, Hubbard model, and many others into circuits with different types of backends.

Example

In the following, we will provide a step-by-step example demonstrating coding with isqdeployer and qsci. Let's get started!

Let study a 4-site system with Hamiltonian given by $H=t\sum_{i=1}^{3}(c^\dagger_{i\sigma}c_{i+1\sigma}+h.c.)-\mu\sum_{i\sigma} c^\dagger_{i\sigma } c_{i\sigma}+U\sum_{i}c^\dagger_{i\uparrow}c_{i\uparrow}c^\dagger_{i\downarrow}c_{i\downarrow},$

example

where we have defined $i=4$ means $i=0$ .We define the Hamiltonian of a Fermion system as follows:

from qsci.secondQuantization.hamiltonian.fermion import SpinfulHamiltonian
hamHubbard = SpinfulHamiltonian()

The hopping terms are acting like a ring. We set the parameters as $t=-1.0$ , $\mu=1.0$ , $U=2.0$ by:

t=-1.0
U=2.0
mu=1.0
for i in range(4): hamHubbard.setHopping(i=i,j=(i+1)%4,tij=t) # set hopping
for i in range(4): hamHubbard.setCoulombU(i=i,U=U) # set U
for i in range(4): hamHubbard.setOnSiteEnergy(i=i,m=-mu) # set mu

Since there is the spin degree of freedom, we need to use 2 qubits to represent a single site in the Hubbard model. The initial state, as shown in figure, is set by:

cir.X(0).X(1)

Next, we need to choose the Jordan-Wigner (JW) transformation to encode the Hamiltonian into a circuit model:

from qsci.secondQuantization.qubitEncoding import JordanWignerEncoding
jw = JordanWignerEncoding() # use Jordan-Wigner encoding
ham = hamHubbard.exportPauliOperator(jw).exportPauliHamiltonian() # Pauli Hamiltonian

The following section is similar to previous examples, where the simulation of time evolution (controlled by the exp factor) is managed by isqdeployer. A complete runnable code is provided below:

from qsci.secondQuantization.hamiltonian.fermion import SpinfulHamiltonian
from qsci.secondQuantization.qubitEncoding import JordanWignerEncoding
from isqdeployer.circuit import Circuit
from isqdeployer.circuitDeployer.pauliHamiltonian import Deployer
import matplotlib.pyplot as plt
import numpy as np 
#---- parameters of Hubbard model ----
t=-1.0
U=2.0
mu=1.0
#-------------- set Hubbard model Hamiltonian -------------------
hamHubbard = SpinfulHamiltonian()
for i in range(4): hamHubbard.setHopping(i=i,j=(i+1)%4,tij=t) # set hopping
for i in range(4): hamHubbard.setCoulombU(i=i,U=U) # set U
for i in range(4): hamHubbard.setOnSiteEnergy(i=i,m=-mu) # set mu

#-------------- encode Hamiltonian into Pauli type
jw = JordanWignerEncoding() # use Jordan-Wigner encoding
ham = hamHubbard.exportPauliOperator(jw).exportPauliHamiltonian() # Pauli Hamiltonian
#-------------- run on circuit
cir = Circuit(nq=8,isInputArg=True) # 4 site Hubbard model map into 8 qubits
cir.X(0).X(1) #-- set initial state 
dpl = Deployer(circuit=cir) 
dpl.expHt(h=ham,t=cir.getInputArg('F',0),N=50) 
cir.setMeasurement(range(8))

#---- run job
tList = np.arange(0,5,0.1)
results = cir.runJob(paramList=[{'F':[t]} for t in tList])

#--- plot 
def N(res,i):
    r = 0 
    for k,v in res.items():
        r += (k[2*i] == "1")*v
    return r 
niList = [ [N(res,i) for res in results] for i in range(4) ]
for i in range(4):plt.plot( tList, niList[i], '-o', label=f"$N_{{{i}\\uparrow}}$" )
plt.legend()
plt.title("simulation with N=50")
plt.xlabel("t")
plt.ylabel("particle number")
plt.show()

The results with $N=50$ and $N=500$ are shown below. We also calculate the result by exact method (classical simulation). As we can see, for the $N=50$ case, the result at a small time region agrees with the exact results very well. At the large time region, the time slide $\tau=T/N$ becomes large, therefore the error increase. As one can see from the geometry of the system, sites 1 and 3 are symmetric, therefore we know that $\langle N_{1\sigma}\rangle=\langle N_{3\sigma}\rangle$ . To obtain a better result, we can set $N=500$ as shown in the figure.