[Kwant] Hubbard model and self-consistent procedure

[郭宇飞]

Dear Kwant developers and users,

Thank you for your contributions to the great software. I am trying to 
reproduce the result demonstrated in Figure 7(a) 
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.79.205430

 

The model includes hopping term t and Hubbard potential U in a hexagonal 
lattice. I did the self-consistent process by Kwant to calculate the spin 
density distribution of the system. But the results seem to be incorrect on the 
distribution of spin density.

 

I attached my code:

sys = my_sys().finalized()

def Fermi(E, mu, T):

    return 1/(exp((E - mu)/T) + 1)

def sorted_eigs(ev):

    evals,evecs=ev

    evals,evecs=map(np.array,zip(*sorted(zip(evals,evecs.transpose()))))

    return evals,evecs.transpose()

def Hartree_Fock(Du,Dd,ncc=200):   

    m=50

    H=sys.hamiltonian_submatrix(sparse=True)

    for nc in range(ncc):

        H0=H

        D=list(itertools.chain(*zip(Du,Dd)))

        H1=np.diag(D)

        H2=H0+H1

        evals,evecs=sorted_eigs(sla.eigsh(H2,k=m,sigma=0))

        up=evecs[::2]

        down=evecs[1::2]

        fermi_up= Fermi(evals,mu, T).reshape(1,m)

        fermi_down= Fermi(evals,mu, T).reshape(1,m)

        Du=(up*up.conj()*fermi_up).sum(axis=1,dtype='double')

        Dd=(down*down.conj()*fermi_down).sum(axis=1,dtype='double')

        if nc%5==0:

            print(f"N_step {nc} Loop")

    return Du,Dd

Du=0.5*np.ones(605)

Dd=np.zeros(605)

Du,Dd=Hartree_Fock(Du,Dd)

Dz=Du-Dd

Dv=Du+Dd   

kwant.plotter.density(sys,Dd,relwidth=0.02,colorbar=True)

kwant.plotter.density(sys,Du,relwidth=0.02,colorbar=True)

kwant.plotter.density(sys,Dz,relwidth=0.02,colorbar=True)

kwant.plotter.density(sys,Dv,relwidth=0.02,colorbar=True)

 

1.     In the above code I use sys.hamiltonian_submatrix in Kwant to calculate 
the Hamitonian of the system. This Hamiltonian is not only a scattering region 
but it might have lead part in it, so how do I distinguish the matrix elements 
in wires and scattering center? Is it would induce the error in my results?

2.     Since I did not find the self-consistent process of Hubbard model in the 
kwant documentation, I set a loop to calculate the spin density distribution of 
the system and draw the spin density distribution with kwant.plotter.density in 
kwant, I am not sure whether this part is correct. Anyone can give some 
suggestions or comments to have a different method on this part?

 

Thank you in advance for your kindly help,

 

Best wishes,

Yufei Guo