Hi Jannis,

You mentioned that you system is infinite and has translational symmetry, yet you are using 'hamiltonian_submatrix'.

'hamiltonian_submatrix' probably does not return what you want; what should it return when the system is infinite, as it is in your case? Probably you want to construct the k-space Hamiltonian and work with that. You could do that with the following function:

    def H_k(k, params):

        H_0 = syst.cell_hamiltonian(params=params)

        V = np.zeros_like(H_0)

        V += syst.inter_cell_hopping(params=params)

        return H_0 + V * exp(1j * k) + V.conjugate().transpose() * exp(-1j * k)

If you are just interested in the eigenvalues/vectors you can also use 'kwant.physics.Bands', which will allow you to evaluate the eigenvalues/vectors given a k (check out the documentation to find out how to get the eigenvectors in addition to the eigenvalues).

I would first rectify this before looking any further.

Happy Kwanting,

Joe

I am trying to calculate the total energy of the electrons in an s-wave superconducting lead with spins. But I am not sure, if my method is correct.

Currently I separate the particle and hole Hamiltonians with a projector

    H = syst.hamiltonian_submatrix(params=params)
    P_electron = np.kron(np.eye(L+1), np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0,0,0], [0, 0,0,0]]))
    H_electron = P_electron.conjugate().transpose() @ H @ P_electron

where L is the number of sites. L+1 because I get a dimension mismatch otherwise. My guess here is, that kwant generates an extra site because of the translational symmetry. For the conservation laws I used

    particle_hole=pauli.sxs0, conservation_law=-pauli.szs0

where pauli.sisj is the kroneker product of the i-th and the j-th pauli matrix. For my Hamiltonian I used the base (e_up, e_down, h_down, h_up).

Then I calculate the eigenvalues of H_electron and multiply them with the fermi distribution to get the total energy.

    spectrum=  np.linalg.eigvalsh(H_electron)
    fermi_dist=kwant.kpm.fermi_distribution(spectrum, params["mu"], T)
    total_energy=np.dot(fermi_dist,spectrum)

Now to the part, why I suspect a mistake in my method. I use a spin-dependent hopping term and scan over a parameter b, that changes the spin dependency. When doing so the qualitative behavior of the total energy E_total(b) changes with the size of my unit cell. But all sites in my unit cell are identical. So I would suspect, that the results should only differ by a factor. For large systems the curvature of E_total(b) converges, but for small systems the results vary quit a lot.

I'd really appreciate some help on this.
Best regards,
Jannis