Hello, I'm confused by some problem when doing time-dependent spin related transport calculations by tkwant. Let me briefly introduce what the main problem is. I'm now considering a 1D chain model with some magnetic spin moment on some of the sites. By applying a voltage to one lead, I want to know how the spin accumulate near the other lead. My problem is that in my code the calculated spin density and charge density near the second lead are always same which I believe should not be true because this means that only one spin channel is always occupied and the other is always empty! My code is as follows: import tkwant import kwant from matplotlib import pyplot as plt from math import pi, sqrt, cos, sin from cmath import atan import numpy as np ############ spin matrix ############# s0=np.array([[1,0],[0,1]]) sx=np.array([[0,1],[1,0]]) sy=np.array([[0,-1j],[1j,0]]) sz=np.array([[1,0],[0,-1]]) length=10 l_M=1 M=0.2 t=1 times=range(200) lat = kwant.lattice.general([(5,0,0),(0,5,0),(0,0,1)], [(2.5,2.5,0)],norbs=2) def scattering_region(pos): x,y,z=pos return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z>=0 and z<=length-1 def mag_region0(pos): x,y,z=pos return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z<=-1 and z>=-l_M def mag_region1(pos): x,y,z=pos return abs(x-2.5)<=1 and abs(y-2.5)<=1 and z>=length and z<=length+l_M-1 syst=kwant.Builder() syst[lat.shape(scattering_region, (2.5, 2.5, 0))]=0*s0 syst[lat.shape(mag_region0, (2.5, 2.5, -1))]=M*sz syst[lat.shape(mag_region1, (2.5, 2.5, length))]=M*sz for i in range (length+2*l_M-1): syst[lat.sublattices[0](0,0,i-l_M),lat.sublattices[0](0,0,i-l_M+1)]=t*s0 def lead_region(pos): x,y,z=pos return abs(x-2.5)<=1 and abs(y-2.5)<=1 lead=kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0,0,-1))),conservation_law=-sz) lead[lat.shape(lead_region, (2.5,2.5,0))]=0*s0 lead[lat.neighbors()]=t*s0 def phi(time): Vb, tau= 0.6, 20 if time > tau: return Vb*time - Vb*tau/2 return Vb*time*time/2/tau syst.attach_lead(lead) syst.attach_lead(lead.reversed()) kwant.plot(syst,site_size=0.2,site_lw=0.01,hop_lw=0.05) tkwant.leads.add_voltage(syst,0,phi) syst=syst.finalized() def current_cut(site_to, site_from): return site_from.pos[2]<=length-1 and site_to.pos[2]>length-1 def density_cut(site): return site.pos[2]>=length and site.pos[2]<=length+l_M-1 current_operator=kwant.operator.Current(syst,where=current_cut) charge_density_operator=kwant.operator.Density(syst,s0,where=density_cut) spin_density_operator=kwant.operator.Density(syst,sz,where=density_cut) occupy=tkwant.manybody.lead_occupation(chemical_potential=-0.5) print('make up of state ...') state=tkwant.manybody.State(syst,tmax=max(times),occupations=occupy,refine=False,combine=True) Currents=[] Charge=[] Spin=[] for time in times: print('%.4f percent ...' % (100*times.index(time)/len(times))) state.evolve(time) current=state.evaluate(current_operator) charge=state.evaluate(charge_density_operator) spin=state.evaluate(spin_density_operator) Currents.append(current) Charge.append(charge) Spin.append(spin) np.savetxt('Times_linear_chain.txt',times) np.savetxt('Current_linear_chain.txt',Currents) np.savetxt('Charge_density_linear_chain.txt',Charge) np.savetxt('Spin_density_linear_chain.txt',Spin) plt.figure() plt.plot(times,Currents) plt.xlabel('time $t$') plt.ylabel('Current $I$') plt.show() plt.figure() plt.subplot(211) plt.plot(times,Charge) plt.xlabel('time $t$') plt.ylabel('Charge density') plt.subplot(212) plt.plot(times,Spin) plt.xlabel('time $t$') plt.ylabel('Spin density') plt.show() Does anyone know where the problem is? -Yizhou LIU