Hi Adel, Running the code you provided also produces the following warning.

UserWarning: Hamiltonian breaks Conservation law, ignoring the symmetry in the computation.

This should explain what is happening. Best, Anton On Sun, 15 Aug 2021 at 14:19, Adel Belayadi wrote:

Dear Anton;

Thank you for your reply

I have used the conservation_law as you have suggested. In fact, a dictionary or function definition of the conservation law is working fine now. I am not getting an error while defining the conservation law. However, when calling smatrix.transmission((1, 1), (0, 1)) which stands for spin up-up transmittance it raises: index 2 is out of bounds for axis 0 with size 2.

Please see for example the following script which have two site families, one with 2 orbitals and the second with 6 orbitals

Cheers, Adel

import kwant from matplotlib import pyplot import numpy as np from math import sqrt lat = kwant.lattice.honeycomb(a=1, norbs=[2, 6] , name=["a", "b"]) a, b =lat.sublattices def make_sys(l, w, cons_law): def rect(pos): x, y = pos return abs(x) < l and abs(y) < w sys=kwant.Builder() sys[a.shape(rect, (0, 0))] = np.eye(2) sys[b.shape(rect, (0, 0))] = np.eye(6) ab_hopps= (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in ab_hopps]] = np.eye(2,6) sys.eradicate_dangling() sys[a.neighbors(1)] = np.eye(2) sys[b.neighbors(1)] = np.eye(6) ### creating leads def lead_shape(pos): x, y = pos return abs(y)< w sym = kwant.TranslationalSymmetry(lat.vec((-1, 0))) sym.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym, conservation_law=cons_law) lead[a.shape(rect, (0, 0))] = np.eye(2) lead[b.shape(rect, (0, 0))] = np.eye(6) lead[[kwant.builder.HoppingKind(*hopping) for hopping in ab_hopps]] = np.eye(2,6) lead.eradicate_dangling() lead[a.neighbors(1)] = np.eye(2) lead[b.neighbors(1)] = np.eye(6) sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys

syst = make_sys(l=3, w=2, cons_law= None) kwant.plot(syst ,fig_size=(8, 8)) pyplot.show()

print('************ calling the conservation law **************') sigma_z = np.array([[1, 0], [0, -1]]) cons_law = {a: -sigma_z , b: -np.kron(sigma_z, np.eye(3))} syst = make_sys(l=3, w=2, cons_law=cons_law) smatrix = kwant.smatrix(syst.finalized(), 0.01) total_trans = smatrix.transmission(1, 0) up_up_trans = smatrix.transmission((1, 0), (0, 0)) dn_dn_trans = smatrix.transmission((1, 1), (0, 1)) up_dn_trans = smatrix.transmission((1, 1), (0, 0)) dn_up_trans = smatrix.transmission((1, 0), (0, 1)) print(total_trans) print(up_up_trans) print(dn_dn_trans) print(up_dn_trans) print(dn_up_trans) print(round(total_trans, 3)==round(up_up_trans+dn_dn_trans+up_dn_trans+dn_up_trans, 3))

Le sam. 14 août 2021 à 17:06, Anton Akhmerov a écrit :

...

Hi Adel,

Please check the documentation [1]. Since you have multiple site families, you need to specify the conservation law as a dictionary, so {metal_like_sub: -np.kron(sigma_z, np.eye(3)), C_like_sub1: -sigma_z, C_like_sub2: -sigma_z}.

Best, Anton

https://kwant-project.org/doc/1/reference/generated/kwant.builder.Builder#kw...

On Sat, 14 Aug 2021 at 16:08, Adel Belayadi wrote:

...

Dear Kwant user, Good day and I hope this email finds you well. Using conservation law was fully explained, in the Kwant mailing list, while we are having sublattices with the same degree of freedom (same number of orbitals). For instance, if we have two sites with 3 orbitals in the case of a spin full Hamiltonian, the onsite matrix is (onsite_matrix(6*6) for each site ). So in this case, the conservation law will be straightforward: conservation_law=-np.kron(sigma_z, np.eye(3)).But what would be the case if we are dealing with sublattices which have different numbers of orbitals?

In my case, I am interested in spin conductance. I am using 3 sublattices with different orbitals. Let us say for instance sublattice0 has (6-orbitals), sublattice1 with (2-orbitals) and sublattice2 with (2-orbitals). In fact I am working on hetero-structure with two sites of C-like(pz with spin up and down) and one site of metal-like(3d-orbitals with spin up and down). I can resume my lattice as follow

lat = kwant.lattice.general(prim_vecs, basis=[ [...], [...], [...]] , norbs=[6, 2, 2]) metal-like_sub, C-like_sub1, C-like_sub2 = lat.sublattices

According to my understanding the conservation law would be something like conservation_law=-np.kron(sigma_z, np.eye(5)). However Kwant raises an error with my choice of the conservation law.

Is there something missing with my understanding of the conservation law!. Or in the case of different orbitals per site, we need to define two lattices (lattice1 for spin up and lattice2 for spindown). Thanks in advance, Best

Dear Anton, I have adopted my conservation law to my system (graphene/TMDs like heterostructure). Band calculation is with correct agreement compared to dft calculation. Everything is fine unless we add hopping terms between spin up and spin down. In this case, we cannot call spin transmission for spin_dn_dn; up_dn and dn_up. To answer my doubts, I have used graphene with a rashba spin-orbit coupling (we are hopping elements between spin up and down like in my case) and the same thing happens to me when calling *dn_dn_trans = smatrix.transmission((1, 1), (0, 1)). * It works only if we set params.rashba to zero. *It is not correct to use the conservation law when we have coupling between spins, like in rashba soc!* I have attached below the script code of graphene with intrinsic and rashba spin orbit terms Thanks in advance import kwant import numpy as np from math import sqrt import matplotlib.pyplot as plt from types import SimpleNamespace import warnings warnings.filterwarnings('ignore') ## ========================================================================================================================== ### ============ Pauli matrices =============== ### -------------------------------------------- sigma_0 = np.array([[1, 0], [0, 1]]) # identity sigma_x = np.array([[0, 1], [1, 0]]) sigma_y = np.array([[0, -1j], [1j, 0]]) sigma_z = np.array([[1, 0], [0, -1]]) ### --------- Onsite------- def onsite(site, params): return (params.Ep+params.mass)*sigma_0 if site.family ==a else (params.Ep-params.mass)*sigma_0 ### ============ Hopping terms for the diffent parts of the system =============== ### --------------------------------------------------------------------------- def kinetic(site_1, site_2, params): """ use it for nearest interaction""" return params.hop*sigma_0 ### ------------------------------------------------------------------------- def Rashba_SOC(site_1, site_2, params): """ use it for nearest interaction""" d_12 = site_2.pos - site_1.pos so = (1j*2/3)*params.lamda_rash*(sigma_x * d_12[1] - sigma_y * d_12[0]) return so ### ------------------------------------------------------------------------- def Intrinsic_SOC(site_1, site_2, params): """ use it for second nearest interaction""" so = (+1j/(3*sqrt(3)))*params.lamda_Int*sigma_z return so ### ------------------------------------------------------------------------- def hopping(site_1, site_2, params): return kinetic(site_1, site_2, params)+Rashba_SOC(site_1, site_2, params) ### ------------------------------------------------------------------------- def hopping_lead(site_1, site_2, params): return kinetic(site_1, site_2, params)+Rashba_SOC(site_1, site_2, params) ## ========================================================================================================================== ### ============ Make lattice ============== params = SimpleNamespace() params.lat_constant=0.24595 lat_constant = params.lat_constant lat = kwant.lattice.honeycomb(a=lat_constant, norbs=2, name=["a", "b"]) a, b =lat.sublattices def make_sys(l, w): def rect(pos): x, y = pos return abs(x) < l and abs(y) < w sys=kwant.Builder() sys[lat.shape(rect, (1, 1))] = onsite sys[lat.neighbors(1)] = hopping sys.eradicate_dangling() ## clockwise sites for SOC (+1) clocksites_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) clocksites_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in clocksites_a]] = Intrinsic_SOC sys[[kwant.builder.HoppingKind(*hopping) for hopping in clocksites_b]] = Intrinsic_SOC ### creating leads def lead_shape(pos): x, y = pos return abs(y)< w sym = kwant.TranslationalSymmetry(lat.vec((-1, 0))) sym.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym, conservation_law=-sigma_z) lead[lat.shape(lead_shape, (1, 1))] = onsite lead[lat.neighbors(1)] = hopping_lead lead.eradicate_dangling() # clockwise sites for SOC (+1) clocksites_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) clocksites_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) lead[[kwant.builder.HoppingKind(*hopping) for hopping in clocksites_a]] = Intrinsic_SOC lead[[kwant.builder.HoppingKind(*hopping) for hopping in clocksites_b]] = Intrinsic_SOC sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys, lead ## ----------------- DIFINING TB AND SYSTEM PARAMETERS ---------------------- params = SimpleNamespace(Ep = 0, hop = -1, lat_constant= 0.24595, center=(0, 0) ) ## ----------------- PLOTTING THE SYSTEM ---------------------- syst, lead = make_sys(l=1, w=1) fsyst= syst.finalized() flead= lead.finalized() # from ipywidgets import interact # @interact(Int=(0, 1.5, .2), rash= (0, 2, .2), mass=(0, 0.51, .05),) # def band_vs_Tranv(Int=0, rash =0, mass=0.,): # fig = plt.figure(figsize=(7, 7)) # ax = fig.add_subplot(111) # ax.set_xlabel("momentum [(lattice constant)^-1]", fontsize=15) # ax.set_ylabel("Energy (eV) of the lead",fontsize=20) # params.mass=mass # params.lamda_Int = Int # params.lamda_rash = rash # kwant.plotter.bands(flead, show=False, ax =ax, args=[params]) params.mass=0.5 params.lamda_Int = 0.2 params.lamda_rash = 0.5 energy = params.mass+0.01 smatrix = kwant.smatrix(syst.finalized(),energy, args=[params]) total_trans = smatrix.transmission(1, 0) up_up_trans = smatrix.transmission((1, 0), (0, 0)) print(total_trans) print(up_up_trans) dn_dn_trans = smatrix.transmission((1, 1), (0, 1)) up_dn_trans = smatrix.transmission((1, 1), (0, 0)) dn_up_trans = smatrix.transmission((1, 0), (0, 1)) print(total_trans) print(up_up_trans) print(dn_dn_trans) print(up_dn_trans) Le dim. 15 août 2021 à 15:30, Adel Belayadi a écrit :

Hey Anton, I should have thought of the physics behind my hopping values when calling the conservation law I have used realistic hopping values instead of np.eye(2, 6) and it is working now. Regards, Adel

Le dim. 15 août 2021 à 14:43, Anton Akhmerov a écrit :

...

Hi Adel,

Running the code you provided also produces the following warning.

...

UserWarning: Hamiltonian breaks Conservation law, ignoring the symmetry in the computation.

This should explain what is happening.

Best, Anton

On Sun, 15 Aug 2021 at 14:19, Adel Belayadi wrote:

...

Dear Anton;

Thank you for your reply

I have used the conservation_law as you have suggested. In fact, a

dictionary or function definition of the conservation law is working fine now. I am not getting an error while defining the conservation law. However, when calling smatrix.transmission((1, 1), (0, 1)) which stands for spin up-up transmittance it raises: index 2 is out of bounds for axis 0 with size 2.

...

Please see for example the following script which have two site

families, one with 2 orbitals and the second with 6 orbitals

...

Cheers, Adel

import kwant from matplotlib import pyplot import numpy as np from math import sqrt lat = kwant.lattice.honeycomb(a=1, norbs=[2, 6] , name=["a", "b"]) a, b =lat.sublattices def make_sys(l, w, cons_law): def rect(pos): x, y = pos return abs(x) < l and abs(y) < w sys=kwant.Builder() sys[a.shape(rect, (0, 0))] = np.eye(2) sys[b.shape(rect, (0, 0))] = np.eye(6) ab_hopps= (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in ab_hopps]]

= np.eye(2,6)

...

sys.eradicate_dangling() sys[a.neighbors(1)] = np.eye(2) sys[b.neighbors(1)] = np.eye(6) ### creating leads def lead_shape(pos): x, y = pos return abs(y)< w sym = kwant.TranslationalSymmetry(lat.vec((-1, 0))) sym.add_site_family(lat.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(lat.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym, conservation_law=cons_law) lead[a.shape(rect, (0, 0))] = np.eye(2) lead[b.shape(rect, (0, 0))] = np.eye(6) lead[[kwant.builder.HoppingKind(*hopping) for hopping in ab_hopps]]

= np.eye(2,6)

...

lead.eradicate_dangling() lead[a.neighbors(1)] = np.eye(2) lead[b.neighbors(1)] = np.eye(6) sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys

syst = make_sys(l=3, w=2, cons_law= None) kwant.plot(syst ,fig_size=(8, 8)) pyplot.show()

print('************ calling the conservation law **************') sigma_z = np.array([[1, 0], [0, -1]]) cons_law = {a: -sigma_z , b: -np.kron(sigma_z, np.eye(3))} syst = make_sys(l=3, w=2, cons_law=cons_law) smatrix = kwant.smatrix(syst.finalized(), 0.01) total_trans = smatrix.transmission(1, 0) up_up_trans = smatrix.transmission((1, 0), (0, 0)) dn_dn_trans = smatrix.transmission((1, 1), (0, 1)) up_dn_trans = smatrix.transmission((1, 1), (0, 0)) dn_up_trans = smatrix.transmission((1, 0), (0, 1)) print(total_trans) print(up_up_trans) print(dn_dn_trans) print(up_dn_trans) print(dn_up_trans) print(round(total_trans,

3)==round(up_up_trans+dn_dn_trans+up_dn_trans+dn_up_trans, 3))

...

Le sam. 14 août 2021 à 17:06, Anton Akhmerov <

...

...

Hi Adel,

Please check the documentation [1]. Since you have multiple site families, you need to specify the conservation law as a dictionary, so {metal_like_sub: -np.kron(sigma_z, np.eye(3)), C_like_sub1: -sigma_z, C_like_sub2: -sigma_z}.

Best, Anton

https://kwant-project.org/doc/1/reference/generated/kwant.builder.Builder#kw...

...

On Sat, 14 Aug 2021 at 16:08, Adel Belayadi

wrote:

...

...

Dear Kwant user, Good day and I hope this email finds you well. Using conservation law was fully explained, in the Kwant mailing

...

...

...

In my case, I am interested in spin conductance. I am using 3

sublattices with different orbitals. Let us say for instance

...

sublattice0 has (6-orbitals), sublattice1 with (2-orbitals) and sublattice2 with (2-orbitals). In fact I am working on hetero-structure with two sites of C-like(pz with spin up and down) and one site of

anton.akhmerov+kd@gmail.com> a écrit : list, while we are having sublattices with the same degree of freedom (same number of orbitals). For instance, if we have two sites with 3 orbitals in the case of a spin full Hamiltonian, the onsite matrix is (onsite_matrix(6*6) for each site ). So in this case, the conservation law will be straightforward: conservation_law=-np.kron(sigma_z, np.eye(3)).But what would be the case if we are dealing with sublattices which have different numbers of orbitals? metal-like(3d-orbitals with spin up and down). I can resume my lattice as follow

...

...

...

lat = kwant.lattice.general(prim_vecs, basis=[ [...], [...], [...]]

, norbs=[6, 2, 2])

...

metal-like_sub, C-like_sub1, C-like_sub2 = lat.sublattices

According to my understanding the conservation law would be something like conservation_law=-np.kron(sigma_z, np.eye(5)). However Kwant raises an error with my choice of the conservation law.

Is there something missing with my understanding of the conservation law!. Or in the case of different orbitals per site, we need to define two lattices (lattice1 for spin up and lattice2 for spindown). Thanks in advance, Best