Dear Patrik,
I have some remarks for you:
1) In the definition of your variable 'sea' you are adding a scalar to a
matrix. This will give you a result different from what you expect. Start
by correcting this: gd ----> gd *eye(dim)
2) Call your function mount_vlead after you define your system and leads
vli = syst.sites()
mount_vlead(syst, vli, norbs=no)
3) If you mount a virtual lead, you will get an error when calling the
scattering matrix:
AttributeError: 'SelfEnergyLead' object has no attribute 'modes'
because the leads you added have no modes.
you can ask for a submatrix of the scattering matrix by specify only
your real leads:
smatrix(sys,energy, in_lead=[0,1],out_lead=[0,1])
This also will not help and you will get the same error.
In you are case, you are adding the same lead (so the same self energy)
for all your sites. You can avoid this just by
not adding any virtual lead and instead add directly your self energy to
your potential Ex:
sys[site1]=pot+self_energy
By doing this, your Hamiltonian becomes not Hermitian but this does not
prevent you from calculating the scattering matrix:
smatrix(sys,..........., check_hermiticity=False)
Now, this does not help you in specifying the current flowing through your
imaginary leads and put it to zero.
The only way I see to achieve that, is to use 1D real leads and use a
potential on each one and find the combination which makes the current =0.
For this purpose, you need to use the conductance matrix
conductance_matrix() #check documentation.
I hope this helps.
Adel
On Thu, Jul 27, 2017 at 2:49 PM, Patrik Arvoy <arv...@gmail.com> wrote:
>
> Dear users and developers,
>
> I would like to attach Büttiker’s virtual lead (the leads with zero net
> currents) to all sites in the scattering region to mimic dephasing
> processes.
>
> The attached code works fine as long as the switch 'vl=False' (without
> attaching any virtual lead). I get equal conductance for up and down spin.
>
> As I understood, attaching virtual leads can be done via the following
> builder class in kwant:
> https://kwant-project.org/doc/1/reference/generated/kwant.bu
> ilder.SelfEnergyLead#kwant.builder.SelfEnergyLead
>
> I have the retarded self energy of each virtual lead (- 1j*0.0025) that I
> can replace with '*selfenergy_func*' in this builder class (am I right?).
> Now if I change 'vl=True' :
>
> 1- I get an error message, related to symmetry but there is no argument to
> define symmetry...
> 2- how can I ensure zero net currents through these virtual leads?
>
> Thanks in advance for any help
> Patrik
>
>
> --------------------------------------------------
> from mpl_toolkits.mplot3d import Axes3D
> from scipy.spatial import *
> from matplotlib import rcParams
> from numpy import *
> from numpy.linalg import *
> import pickle
> import sys
> import os
> import string
> import heapq
> import kwant
> import tinyarray
> from matplotlib import pyplot
>
>
>
>
>
> chiral=True
> if chiral:
> p = pi/5 #phi
> t = 0.66 #theta
> a = 0.34
> x = 1.4
> e1 = 0
> e2 = 0.3
> t2=0.1
> t1=-x*t2
> t0 = 2
> lam=-0.08
> t_so1 = 0.01 #spin-orbit coupling param
> t_so2 = x*t_so1 #spin-orbit coupling param
> tl=tr=0.3
> gd = 0.005 #the gamma_d, dephasing parameter
> N = 30
> sigma_0 = tinyarray.array([[1, 0], [0, 1]])
> sigma_x = tinyarray.array([[0, 1], [1, 0]])
> sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
> sigma_z = tinyarray.array([[1, 0], [0, -1]])
> no=2 #number of orbitals
> def sigma_v1(ap): # pauli metrix along the vertical axis
> value=sigma_z*cos(t)+sin(t)*(sigma_x*sin(ap)-sigma_y*cos(ap))
> return value
>
> def sigma_v2(ap): # pauli metrix along the vertical axis
> value=sigma_z*cos(t)-sin(t)*(sigma_x*sin(ap)-sigma_y*cos(ap))
> return value
>
> def family_color(sites):
> return 'black' #if site.family == sites
>
> def hopping_lw(site1, site2):
> return 0.08
>
>
> class Amorphous(kwant.builder.SiteFamily):
> def __init__(self, coords):
> self.coords = coords
> super(Amorphous, self).__init__("amorphous", "",no)
>
> def normalize_tag(self, tag):
> try:
> tag = int(tag[0])
> except:
> raise KeyError
>
> if 0 <= tag < len(coords):
> return tag
> else:
> raise KeyError
>
> def pos(self, tag):
> return self.coords[tag]
>
> coords=[(0.0000000000, 0.0000000000, 0.0000000000), (-0.1336881039,
> 0.4114496766, 0.3400000000), (-0.4836881039, 0.6657395614, 0.6800000000),
> (-0.9163118961, 0.6657395614, 1.0200000000), (-1.2663118961, 0.4114496766,
> 1.3600000000), (-1.4000000000, 0.0000000000, 1.7000000000), (-1.2663118961,
> -0.4114496766, 2.0400000000), (-0.9163118961, -0.6657395614, 2.3800000000),
> (-0.4836881039, -0.6657395614, 2.7200000000), (-0.1336881039,
> -0.4114496766, 3.0600000000), (0.0000000000, -0.0000000000, 3.4000000000),
> (-0.1336881039, 0.4114496766, 3.7400000000), (-0.4836881039, 0.6657395614,
> 4.0800000000), (-0.9163118961, 0.6657395614, 4.4200000000), (-1.2663118961,
> 0.4114496766, 4.7600000000), (-1.4000000000, 0.0000000000, 5.1000000000),
> (-1.2663118961, -0.4114496766, 5.4400000000), (-0.9163118961,
> -0.6657395614, 5.7800000000), (-0.4836881039, -0.6657395614, 6.1200000000),
> (-0.1336881039, -0.4114496766, 6.4600000000), (0.0000000000, -0.0000000000,
> 6.8000000000), (-0.1336881039, 0.4114496766, 7.1400000000), (-0.4836881039,
> 0.6657395614, 7.4800000000), (-0.9163118961, 0.6657395614, 7.8200000000),
> (-1.2663118961, 0.4114496766, 8.1600000000), (-1.4000000000, 0.0000000000,
> 8.5000000000), (-1.2663118961, -0.4114496766, 8.8400000000),
> (-0.9163118961, -0.6657395614, 9.1800000000), (-0.4836881039,
> -0.6657395614, 9.5200000000), (-0.1336881039, -0.4114496766, 9.8600000000),
> (-1.4000000000, 0.0000000000, 0.0000000000), (-1.2663118961, -0.4114496766,
> 0.3400000000), (-0.9163118961, -0.6657395614, 0.6800000000),
> (-0.4836881039, -0.6657395614, 1.0200000000), (-0.1336881039,
> -0.4114496766, 1.3600000000), (0.0000000000, -0.0000000000, 1.7000000000),
> (-0.1336881039, 0.4114496766, 2.0400000000), (-0.4836881039, 0.6657395614,
> 2.3800000000), (-0.9163118961, 0.6657395614, 2.7200000000), (-1.2663118961,
> 0.4114496766, 3.0600000000), (-1.4000000000, 0.0000000000, 3.4000000000),
> (-1.2663118961, -0.4114496766, 3.7400000000), (-0.9163118961,
> -0.6657395614, 4.0800000000), (-0.4836881039, -0.6657395614, 4.4200000000),
> (-0.1336881039, -0.4114496766, 4.7600000000), (0.0000000000, -0.0000000000,
> 5.1000000000), (-0.1336881039, 0.4114496766, 5.4400000000), (-0.4836881039,
> 0.6657395614, 5.7800000000), (-0.9163118961, 0.6657395614, 6.1200000000),
> (-1.2663118961, 0.4114496766, 6.4600000000), (-1.4000000000, 0.0000000000,
> 6.8000000000), (-1.2663118961, -0.4114496766, 7.1400000000),
> (-0.9163118961, -0.6657395614, 7.4800000000), (-0.4836881039,
> -0.6657395614, 7.8200000000), (-0.1336881039, -0.4114496766, 8.1600000000),
> (0.0000000000, -0.0000000000, 8.5000000000), (-0.1336881039, 0.4114496766,
> 8.8400000000), (-0.4836881039, 0.6657395614, 9.1800000000), (-0.9163118961,
> 0.6657395614, 9.5200000000), (-1.2663118961, 0.4114496766, 9.8600000000)]
> amorphous_lat = Amorphous(coords)
>
> syst = kwant.Builder()
>
>
> for i in range(N):
> syst[amorphous_lat(i)] = e1*sigma_0
> syst[amorphous_lat(N+i)] = e2*sigma_0
> syst[amorphous_lat(i), amorphous_lat(N+i)] = lam*sigma_0
> if i > 0:
> syst[amorphous_lat(i), amorphous_lat(i-1)] = t1*sigma_0 +
> 1j*t_so1*(sigma_v1(i*p)+sigma_v1((i-1)*p))
> syst[amorphous_lat(N+i),amorphous_lat(N+i-1)] = t2*sigma_0 +
> 1j*t_so2*(sigma_v2(i*p)+sigma_v2((i-1)*p))
>
>
>
>
> #the section that cause the problem
> vl=False
> if vl:
> def mount_vlead(syst, vlead_interface, norbs):
> """Mounts virtual lead to interfaces provided.
>
> :sys: kwant.builder.Builder
> An unfinalized system to mount leads
> :vlead_interface: sequence of kwant.builder.Site
> Interface of lead
> :norb: integer
> Number of orbitals in system hamiltonian.
> """
> dim = len(vlead_interface)*norbs
> sea = zeros((dim, dim), dtype=float) - 1j*0.5*gd
>
> def selfenergy_func(energy, args=()):
> return sea
>
> vlead = kwant.builder.SelfEnergyLead(selfenergy_func,
> vlead_interface)
> syst.leads.append(vlead)
>
> #attaching virtual lead to all sites
> vli = syst.sites()
> mount_vlead(syst, vli, norbs=no)
>
>
>
>
>
>
>
>
> lat=kwant.lattice.cubic(a, norbs=no)
>
> syst[lat(0, 0, -1)] = e1*sigma_0
> syst[amorphous_lat(0), lat(0, 0, -1)] = tl*sigma_0
>
> sym = kwant.TranslationalSymmetry([0, 0, -a])
> dn_lead = kwant.Builder(sym, conservation_law=-sigma_z)
>
> dn_lead[lat(0, 0, -2)] = e1*sigma_0
> dn_lead[lat.neighbors()] = t0*sigma_0
> syst.attach_lead(dn_lead)
>
> prim_vecs=tinyarray.array([(a,0.,0.),(0.,a,0.),(0.,0.,a)])
> offset=tinyarray.array((-1.2663118961, 0.4114496766,0.0))
> lat2=kwant.lattice.Monatomic(prim_vecs, offset, norbs=no)
>
> syst[lat2(0, 0, N)] = e1*sigma_0
> syst[amorphous_lat(2*N-1), lat2(0, 0, N)] = tr*sigma_0
>
> sym1 = kwant.TranslationalSymmetry([0, 0, a])
> up_lead = kwant.Builder(sym1, conservation_law=-sigma_z)
> up_lead[lat2(0, 0, N+1)] = e1*sigma_0
> up_lead[lat2.neighbors()] = t0*sigma_0
> syst.attach_lead(up_lead)
>
> system=kwant.plot(syst, site_lw=0.1, site_color=family_color,
> hop_lw=hopping_lw)
>
>
> trans=True
> if trans:
> syst = syst.finalized()
> energies = []
> datau = []
> datad = []
>
> for ie in range(-320,520):
> energy = ie * 0.001
> smatrix = kwant.smatrix(syst, energy=energy)
> energies.append(energy)
> Gu=smatrix.transmission((1, 0), 0)
> Gd=smatrix.transmission((1, 1), 0)
> datau.append(Gu)
> datad.append(Gd)
>
> fig = pyplot.figure()
> pyplot.plot(energies, datau, 'r--')
> pyplot.plot(energies, datad, 'b:')
> pyplot.legend(['Gu', 'Gd'], loc='upper left')
> pyplot.xlim([-0.4,0.65])
> pyplot.ylim([-0.03,1.0])
> pyplot.show()
>
>
>
--
Abbout Adel
