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()