Dear all,
I am trying to calculate the electrical conductivity in a cubic system in
different temperatures, where chemical potential is set to zero (see below).
I am using some of the tips given in the mailing list, I have a couple of
questions though.
1- When I solve V=RI to get the voltages, the voltage values are very
large. Does it make sense or am I making a mistake somewhere?
2- In order to calculate the conductance, I need to integrate over an
energy interval. should I integrate over the entire energy spectrum of the
system? The code is very slow, how do I choose the energy interval and the
step (least requirements) to ensure something physically meaningful yet
speed up the code, any idea?
Bests
Chung
--------------------------------------------------------
import matplotlib.pyplot as pyplot
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 time
import itertools as IT
from scipy.spatial import cKDTree
start = time.time()
sigma_x = array([[0, 1], [1, 0]])
sigma_y = array([[0, -1j], [1j, 0]])
sigma_0 = array([[1, 0], [0, 1]])
sigma_z = array([[1, 0], [0, -1]])
a=1.0
t=1.0
L = 4
cutoff_dist = 8
theta=0.24434609527920614
lam=1.0E-3
v0=0.1
hbns= 6.582119514E-16*9.0E22
alpha = 4.
mu = 0.
e0 = 0.6
nu=3.0E11
K=8.6173303E-5
nuoft = 1
#tmp = 40.
se = zeros((int(L**3.), nuoft), float)
mo = zeros((int(L**3.), nuoft), float)
for ii in range(nuoft):
se[:,ii] = linspace(-e0, e0, int(L**3.))
random.shuffle(se[:,ii])
mo[:,ii] = random.uniform(0., theta, int(L**3.))
TL=arange(20, 322, 20)
pltl = TL**(-0.25)
sigma_xx=[]
sigma_xy=[]
tosim = 3.8740458655E-5
fintemp=True
if fintemp:
def dF(energy, Ef, beta):
return 0.25*beta/(cosh(0.5*beta*(energy-Ef))**2)
def get_trans(beta, lin, lout, Ef=0.):
energies=linspace(-e0, e0, 10)
#energies=[i/(15*beta) for i in range (-150,150)]
trans=[]
for energy in energies:
smatrix = kwant.smatrix(sys, energy)
T=smatrix.transmission(lin, lout)
trans.append(T)
return trans,energies
def Use_temp(beta, lin, lout, Ef=0.):
trans,energies=get_trans(beta, lin, lout, Ef=0.)
spd=abs(energies[0]-energies[1])
couple_T_E=list(zip(trans,energies))
Conductance=[]
Tr=[T*dF(energy,Ef=Ef,beta=beta) for T,energy in couple_T_E]
g=trapz(Tr, energies, dx=spd)
return g
def clcte_sigmas(sys):
G = [[Use_temp(beta, i, j, Ef=0.) for i in range(6)] for j in
range(6)]
G -= diag(sum(G, axis=0))
temp = [0, 1, 2, 3, 4]
G = G[temp, :]
G = G[:, temp]
r = linalg.inv(G)
V = r @ array([1, -1, 0, 0, 0])
print(V)
E_y = V[3] - V[2]
E_x = V[1] - V[0]
s_xx = (E_x / (E_x**2 + E_y**2))#*(tosim/(L*1E-7))
s_xy = (E_y / (E_x**2 + E_y**2))#*(tosim/(L*1E-7))
sigma_xx.append(s_xx)
sigma_xy.append(s_xy)
return s_xx, s_xy
def spin_conductance(G, lead_out, lead_in, sigma):
"""Calculate the spin conductance between two leads.
Parameters
----------
G : an instance of `kwant.solvers.common.GreensFunction`
The Greens function of the system as returned by
`kwant.greens_function`.
lead_out : integer
The lead where spin current is collected
lead_in : integer
The lead where spin current is injected
sigma : `numpy.ndarray` of shape (2, 2)
The Pauli matrix of the quantization axis along
which to measure the spin current
Notes
-----
Calculates the spin conductance between two leads
p and q according to:
G_{pq} = Tr[ σ_{α} Γ_{q} G_{qp} Γ_{p} G^†_{qp} ]
Where Γ_{q} is the coupling matrix to lead q ( = i[Σ - Σ^†] )
and G_{qp} is the submatrix of the system's Greens function
between sites interfacing with lead p and q.
(see also the function below (spin_conductance2))"""
ttdagger = G._a_ttdagger_a_inv(lead_out, lead_in)
sigma_z_matrix = kron(eye(ttdagger.shape[0]//2), sigma)
return trace(sigma_z_matrix.dot(ttdagger)).real
def family_color(sites):
return 'black'
def hopping_lw(site1, site2):
return 0.08
def onsite(site):
oe=2.0*e0* kwant.digest.uniform(repr(site))
if oe <= e0:
return oe*sigma_0
else:
return (e0-oe)*sigma_0
def hop1(site1, site2, tempt):
ni = int(site1[0]*L**2. + site1[1]*L + site1[2])
nj = int(site2[0]*L**2. + site2[1]*L + site2[2])
dij = norm(array(site1 - site2))
rs = dij**2.0
eij = (abs(se[ni,0]-mu)+abs(se[nj,0]-mu)+abs(se[ni,0]-se[nj,0]))
zij = nu*exp(-2.0*alpha*(dij))*exp(-(0.5/(K*tempt))*eij)
return zij*sigma_0
def cuboid_shape(pos):
x, y, z = pos
return 0 <= x < L and 0 <= y < L and 0 <= z < L
def xlead_shape(pos):
x, y, z = pos
return x==0 and 0 <= y < L and 0 <= z < L
def ylead1_shape(pos):
x, y, z = pos
return 0 <= x < L/2 and y==0 and 0 <= z < L
def ylead2_shape(pos):
x, y, z = pos
return round(L/2) <= x < L and y==0 and 0 <= z < L
prim_vecs=array([(a,0.,0.),(0.,a,0.),(0.,0.,a)])
offset1=array((0.0, 0.0, 0.0))
lat=kwant.lattice.Monatomic(prim_vecs, offset1, norbs=2)
for tmp in TL:
print('T=',tmp,'K')
print('The conductivity is being calculated')
sys = kwant.Builder()
beta = 1.0/(K*tmp)
sys[lat.shape(cuboid_shape, (0, 0, 0))] = 0*sigma_0
sites = set(sys.sites())
for i in sites:
sys[i] = se[int(i.tag[0]*L**2. + i.tag[1]*L +
i.tag[2]),0]*sigma_0
for jump in IT.combinations(sites, 2):
sys[jump[0], jump[1]] = hop1(jump[0].tag, jump[1].tag, tmp)
lead0 = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0)))
lead0[lat.shape(xlead_shape, (0, 0, 0))] = 0*sigma_0
lead0[lat.neighbors()] = -t*sigma_0
sys.attach_lead(lead0)
sys.attach_lead(lead0.reversed())
lead2 = kwant.Builder(kwant.TranslationalSymmetry((0, -a, 0)))
lead2[lat.shape(ylead1_shape, (0, 0, 0))] = 0*sigma_0
lead2[lat.neighbors()] = -t*sigma_0
sys.attach_lead(lead2)
sys.attach_lead(lead2.reversed())
lead3 = kwant.Builder(kwant.TranslationalSymmetry((0, -a, 0)))
lead3[lat.shape(ylead2_shape, (L-1, 0, 0))] = 0*sigma_0
lead3[lat.neighbors()] = -t*sigma_0
sys.attach_lead(lead3)
sys.attach_lead(lead3.reversed())
system=kwant.plot(sys, site_lw=0.0001, site_color=family_color,
hop_lw=hopping_lw)
sys = sys.finalized()
s_xx, s_xy = clcte_sigmas(sys)
print('σ_xx=', s_xx, 'σ_xy=', s_xy)
print()
end = time.time()
print()
print ("The execution time:")
if (end - start < 60):
print (end - start, "seconds")
if (end - start > 60) and (end - start < 3600):
print ((end - start)/60, "minutes")
if (end - start > 3600) :
print ((end - start)/3600, "hours")
print()