Dear all,

I am trying a graphene bulk with magnetic field, and I use

spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
                                     vector_factory=s_factory,
                                     rng=0)
energies, densities = spectrum()

to calculate the dos and compare with some previous studies. My definition of the hopping is:

    def Hop_magnetic(site1,site2):
        B_magneic=200. # units Tesla
        a0=0.142       # nm
        x1,y1=site1.pos
        x2,y2=site2.pos
        xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0    # the units is nm**2
        phb=1j*2.*pi*B_magneic/4.135667    # phi0=h/e=4.135667*e-15 V*s
        return t*exp(xy*phb*0.001)

I found that the effect of the magnetic field is too small and I can not obtain the correct results.

 Could you give me some suggestions?

My code is pasted:

import kwant
from matplotlib import pyplot
from numpy import sqrt,pi,exp
import numpy as np

def make_system(r=8, t=1, tp=-0.1):
    lat = kwant.lattice.honeycomb(norbs=1)
    a, b = lat.sublattices
   
    def circle(pos):
        x, y = pos
        return x**2 + y**2 < 100**2  #-100<x<100 and -100<y<100

    def Hop_magnetic(site1,site2):
        B_magneic=200. #  the unit of  B_magneic is Tesla=V*s/m**2
        a0=0.142       # nm
       
        x1,y1=site1.pos
        x2,y2=site2.pos

        xy=0.5*(x1+x2)*sqrt(3)*a0*(y1-y2)*sqrt(3)*a0    # the units change it to nm**2 
        phb=1j*2.*pi*B_magneic/4.135667    # phi0=h/e=4.135667*e-15 V*s, the unit of  B_magneic is Tesla=V*s/m**2
        return t*exp(xy*phb*0.001)
   
    def Nearest_Hop(site1,site2):
        return Hop_magnetic(site1,site2)

    syst = kwant.Builder()
    syst[lat.shape(circle, (0, 0))] = 0
    syst[lat.neighbors()] = Nearest_Hop
    syst.eradicate_dangling()

    return lat, syst.finalized()

lat, fsyst = make_system()
#kwant.plot(fsyst)
where = lambda s: np.linalg.norm(s.pos) < 1
s_factory = kwant.kpm.LocalVectors(fsyst, where)
spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None,
                                     vector_factory=s_factory,
                                     rng=0)
energies, densities = spectrum()

pyplot.figure()
pyplot.plot(energies, densities)
pyplot.xlabel("energy [t]")
pyplot.ylabel("DOS")
pyplot.show()