# How to realize Su-Schrieffer-Heeger model in Qiskit

This is a very specific question, which I try to implement a simple dimerized tight-binding Hamiltonian on qiskit. The model is one dimensional, and defined as

$$H = \sum_{\langle i,j\rangle} t_{ij} c^\dagger_i c_j$$

where $$\langle i,j\rangle$$ denotes the nearest neighbor coupling with strength $$t_{ij}$$. The important aspect of the model is that the coupling is dimerized in the following sense. If we label the sites as 1,2,3,4...., then the nearest neighbor coupling $$t_{ij}$$ reads, for example,

$$t_{1,2} = 1, t_{2,3} = 2 , t_{3,4} = 1 , t_{4,5} = 2, etc$$

which is alternating.

I try to pretend this is a molecular Hamiltonian and realize it in the same way with "FermionicOperator". However, I am not sure how to properly index the "one-body" integral here, which is the nearest neighbor coupling $$t_{ij}$$ here.

I don't think it really matter how you index your $$t_{ij}$$. As you mentioned you can use FermionicOperator all you need to do is define the one body integral. This of course can be done in many ways. Here is a convenient way.

def ssh_ham(gamma, lamda, n):
sigmax = np.array([[0,1],[1,0]], dtype=np.complex_)
sigmay = np.array([[0,-1j],[1j,0]], dtype=np.complex_)

op_eye_x = np.eye(n)
op_cos_x = 1/2*(np.eye(n, k=1) + np.eye(n,k=-1))
op_sin_x = 1j/2*(np.eye(n, k=1) - np.eye(n,k=-1))

h = np.kron(gamma*op_eye_x + lamda*op_cos_x, sigmax) + np.kron(lamda*op_sin_x, sigmay)

return h


This function would return the one-body integral for you. In this function, the $$\sigma$$ degree of freedom represents the two atoms in the unit cell, and $$n$$ represents how many unit cells you have in the SSH chain, and for the dimmerized limit you asked for you can take $$\gamma = 1,\ \lambda = 2$$. You can use the output of this function fer_op = FermionicOperator(h1 = ssh_ham(1,2,10))

That being said, I'm not sure what the benefit would be. Trying to solve this Hamiltonian using VQE for example is an over-kill since this Hamiltonian can be readily diagonalized from its one-body Hamiltonian. Putting this on a quantum computer would give you $$2n$$ qubits, this is a Hilbert space of $$2^n$$ dimensions, whereas you can diagonalize the the Hamiltonian by diagonalizing the $$2n \times 2n$$ matrix $$t_{ij}$$.

• Hi @A. Jahin, thanks for the help. I do think that the indexing is important. Please see my answer below. I have also checked that the "fer_op" you obtained does not give the same ground state energy, after map it to qubits. Could you please check that? Thanks!
– fagd
Oct 2, 2020 at 18:15
• I had my constants set for the fully dimmerized limit. Sorry you didn't ask about that. I fixed the coefficients . I tried my code for the fully dimmerized limit, (that is $\gamma = 0$ and $\lambda = 1$) and it gives the correct answer. For example for $n = 4$ I get the energy to be $-3$ at half filling. Oct 2, 2020 at 18:40
• Ok, that is strange. With $\gamma=0,\lambda=1, n=4$ and fer_op defined in your method, I did qubitOp = ferOp.mapping(map_type='parity' , threshold=0.00000001) print( NumPyEigensolver( qubitOp ).run()['eigenvalues'] ) . The answer is -1.9576174.
– fagd
Oct 2, 2020 at 18:48
• Sorry, not sure how to format code in the comment...
– fagd
Oct 2, 2020 at 18:50
• Try using 'NumPyMinimumEigensolver' instead? In this case the true many-body ground state for $4$ electrons is the ground state of the full Hamiltonian. Oct 2, 2020 at 18:55

This is the code I wrote, and I do feel that the index of the single particle integral is important. Below, h1 is the single particle integral in the basis ABABAB, where A,B are the sublattices (similar to spin), and I use "unitary_matrix" to change it AAA...BBB... which is the basis for h1_qubit. The construction of "unitary_matrix" can be found in the qiskit function "FermionicOperator".

Now, one important checking is that the ground state energy of two Hamiltonians, one is the fermionic single particle Hamiltonian (which is h1), and one is the Hamiltonian in the qubit basis (which is qubitOp), are identical. For the latter, it is a many-body Hamiltonian, and the ground state is the many-body ground state. If we naively compare this to the "lowest energy" of h1, then we will find that they are different, for h1 is the single-particle Hamiltonian. To properly obtain the many-body ground state from h1, we need to fill all the states below the Fermi-level of the spectrum of h1. Once we have done that, we found that the two Hamiltonians give the same many-body ground state energy.

import numpy as np
from scipy.linalg import eig
from qiskit.aqua.algorithms import VQE, NumPyEigensolver
from qiskit.chemistry import FermionicOperator
from qiskit.aqua.operators import Z2Symmetries

t1 = np.random.rand() # intra
t2 = np.random.rand() # inter
t0 = np.random.rand() # onsite energy

N = 4 # number of unit cell

BC = 1 # BC=1 is PBC, BC=0 is OBC

h1 = np.zeros( (2*N , 2*N ) )

for i in range(0,(2*N-1)):
if i % 2==1:
h1[ i , i ] = t0
h1[ i , i+1 ] = t2
h1[ i+1 , i ] = t2
else:
h1[ i , i ] = -t0
h1[ i , i+1 ] = t1
h1[ i+1 , i ] = t1

h1[2*N-1,2*N-1] = t0

if BC==1:
h1[ 0 , 2*N-1 ] = t2
h1[ 2*N-1 , 0 ] = t2

print( 'The fermionic Hamiltonian is\n' , h1 )

eigh1 = eig(h1)[0]

idx = np.argsort(eigh1)
eigh1 = eigh1[idx]

print( 'The *single* particle spectrum is \n' , eigh1 )

unitary_matrix = np.zeros((h1.shape), h1.dtype)
n = unitary_matrix.shape[0]
j = np.arange(n // 2)
unitary_matrix[j, 2 * j] = 1.0
unitary_matrix[j + n // 2, 2 * j + 1] = 1.0

h1_qubit = unitary_matrix.dot(h1.dot(unitary_matrix.T.conj()))

ferOp = FermionicOperator(h1=h1_qubit )

qubitOp = ferOp.mapping(map_type='parity' , threshold=0.00000001)

numEigenvalue = 1
result = NumPyEigensolver( qubitOp ,  k = numEigenvalue ).run()
eigvals = result['eigenvalues']
eigvecs = result['eigenstates']
print( 'The *many-body* grond state energy is \n ' , eigvals[0] )

print( 'We check the *many-body* ground state energy is the same as the sum of the energy below the fermi level of the *single* particle spectrum')

print( eigvals[0] - sum( eigh1[0:int(len(eigh1)/2)]))