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