Dear Joe, concerning to your answer https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01808.html , I checked that I do not have any sign problem, as those you described, in my system/code. Therefore I am wandering if some fuzzy problem of matching of wavefunctions is known, or expected, to occur when metals (leads) and semimetals (scattering region) are put in contact (perhaps due to symmetries of whatever), not depending on the sign problem. I also attach below the code that I used, and that looks me correct, after various checks on the spectrum of the scattering region and of the leads. I assumed a cubic lattice with periodic boundary conditions along x and z, while the leads are along y. However, I checked that the problem persists also assuming open boundary conditions. Thank you very much for the help and best regards L. ----------------------------------------------------------------------------------- import kwant import math import cmath import scipy.sparse.linalg as sla import numpy # For plotting from matplotlib import pyplot # For matrix support import tinyarray #----------------------------------------------------- # We define Pauli-matrices sigma_0 = tinyarray.array([[1, 0], [0, 1]]) sigma_x = tinyarray.array([[0, 1], [1, 0]]) sigma_y = tinyarray.array([[0, -1j], [1j, 0]]) sigma_z = tinyarray.array([[1, 0], [0, -1]]) # We define the variables and the constants t = 1.0 v = 1.0 L = 13 # length y W = 13 # length x and z k0 = math.pi/10. #k0 = -0.4636476090008061 m0 = math.cos(k0) vtilde = v/math.sin(k0) h = vtilde*m0 + 2*v a=1 #--------------------------------------------------------------------------------------------------------- syst = kwant.Builder() lat = kwant.lattice.cubic(a) for i in range(W): for j in range(L): for k in range(W): # On-site Hamiltonian syst[lat(i, j, k)] = h*sigma_x # Hopping in x-direction if i > 0: syst[lat(i, j, k), lat(i - 1, j, k)] = -(vtilde/2)*sigma_x # Hopping in y-direction if j > 0: syst[lat(i, j, k), lat(i, j - 1, k)] = -(v/2)*((1j)*sigma_y + sigma_x) # Hopping in z-direction if k > 0: syst[lat(i, j, k), lat(i, j, k - 1)] = -(v/2)*((1j)*sigma_z + sigma_x) # I close hopping in x-direction, connecting W-1 and 0 for j in range(L): for k in range(W): syst[lat(0, j, k), lat(W-1, j, k)] = -(vtilde/2)*sigma_x # I close the hopping in z-direction, connecting W-1 and 0 for i in range(W): for j in range(L): syst[lat(i, j, 0), lat(i, j, W-1)] = -(v/2)*((1j)*sigma_z + sigma_x) #--------------------------------------------------------------------------------------------------------- # We define the right lead, along y sym_right_lead = kwant.TranslationalSymmetry((0, a, 0)) right_lead = kwant.Builder(sym_right_lead) for i in range(W): for k in range(W): right_lead[lat(i, 0, k)] = mu*sigma_0 if i > 0: right_lead[lat(i, 0, k), lat(i-1, 0, k)] = -t*sigma_0 if k > 0: right_lead[lat(i, 0, k), lat(i, 0, k-1)] = -t*sigma_0 right_lead[lat(i, 1, k), lat(i, 0, k)] = -t*sigma_0 # I close hopping in x-direction, connecting W-1 and 0 for k in range(W): right_lead[lat(0, 0, k), lat(W-1, 0, k)] = -t*sigma_0 # I close hopping in z-direction, connecting W-1 for i in range(W): right_lead[lat(i, 0, 0), lat(i, 0, W-1)] = -t*sigma_0 syst.attach_lead(right_lead) left_lead = right_lead.reversed() syst.attach_lead(left_lead)