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

Hi Joe, thank you for the fast response. You are probably right, that I need to use the k-space Hamiltonian for the infinite case. Would my method work for a finite system then? Potentially with periodic boundary conditions? The measure, that I am looking for, is the total energy of all occupied particle states. I am not really sure, how to calculate this in an infinite system. First I probably need to isolate the particle Hamiltonian with a projector. But I don't know how to proceed from there in k-space. Best regards, Jannis On 04.12.2019 16:18, Joseph Weston wrote:

Hi again,

I noticed my code example probably won't even run in the generic case; we need to do this wrangling with 'numpy.zeros_like' because 'inter_cell_hopping' returns a rectangular matrix in general (when only a subset of the sites in a unit cell connect via hoppings to the next unit cell). Probably the code should look something like this:

        def H_k(k, params):

            H_0 = syst.cell_hamiltonian(args, params=params)

            hop = syst.inter_cell_hopping(args, params=params)             V = np.zeros_like(H_0)             V[:, : hop.shape[1]] = hop

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

Happy Kwanting,

Joe

Hi Joe

If I understand your question correctly then you need to do do a k-space integration; something like:

∑_n ∫ E_n(k) f(E_n(k)) dk

Where the sum runs over all the bands, the integral runs over the Brillouin zone and 'f' is the occupation at energy E (e.g. Fermi-Dirac). In this way you count the energy contribution for each state, weighted by its occupation. Note that we have elided the density of states, as we are integrating directly in momentum space, as opposed to in energy, and the density of states in 1D is constant in k.

Does that make sense?

It does make sense and helps a lot. Is there a way to calculate E_n(k) of a 2D infinite system in kwant? 2D translational symmetric systems can't be finalized and cell_hamiltionian and inter_cell_hopping don't work with wraparound either. Best regards, Jannis

∫ {\displaystyle \displaystyle \int }