Hi,

In the discussion thread about calculation of current density, it is mentioned for the function defined for current as follows only work for 1 orbital per site. Does the expression work for system with spin? I guess the Hamiltonian matrix element for sites I and j will be a matrix for a system with spin and the wavefunction for a site is a vector. If this is the case, the function still works for systems with spin degree of freedom (or other internal degree of freedom) Ks chan

def current(psi): H_ij = sys.hamiltonian(i, j, args=(B,)) return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag

This is incorrect because you are conflating the *site index* (the integers passed to 'sys.hamiltonian') and *orbital index* (that you are using to index 'psi'. In addition you are preforming element-wise multiplications when you should be doing matrix multiplications. In this specific example, if you have exactly 2 (spin) degrees of freedom on every site then all of the even-indexed orbitals will be spin up and the odd will be spin down, so you would have to do: def current(psi): H_ij = sys.hamiltonian(i, j, args=(B,)) iorbs = slice(i//2, i//2 + 1) jorbs = slice(j//2, j//2 + 1) return -2 * (psi[iorbs].conjugate().dot(H_ij).dot(psi[jorbs])).imag to get the charge current. For the spin current you would have to put a 'sigma_z' somewhere in the above expression. Because this is often confusing for users we have introduced 'operators' into the (soon to be released) Kwant 1.3 so that you don't have to do this kind of index gymnastics (which gets even more complicated if you have a system where different site families have different numbers of orbitals) Happy Kwanting, Joe