Dear Enrique,
Below is a bit of code (basically, it's not mine, it was also taken
from Kwant discussions), that plots the current density.
Feel free to improve it.
Best wishes,
Sergey
from scipy.sparse import spdiags
def current2d(fsys,lead=0, params=None, h=None, psi=None):
x, y = numpy.array([s.pos for s in fsys.sites]).T
n = len(x)
xd = spdiags(x, 0, n, n)
yd = spdiags(y, 0, n, n)
if h is None:
h = fsys.hamiltonian_submatrix(sparse=True, args=params)
Jx = h * xd - xd * h
Jy = h * yd - yd * h
if psi is None:
wf = kwant.wave_function(fsys, 0, args=params)
# psi = numpy.vstack([wf(l) for l in leads])
psi = wf(lead)
jx = numpy.zeros_like(psi, dtype=float)
jy = numpy.zeros_like(psi, dtype=float)
j2 = numpy.zeros(len(fsys.sites), dtype=float)
for f in xrange(len(psi)):
psi_d_conj = spdiags(psi[f].conj(), 0, n, n)
jx[f] = -2*numpy.imag( psi_d_conj * (Jx * psi[f]) )
jy[f] = -2*numpy.imag( psi_d_conj * (Jy * psi[f]) )
## plot the wave-function for the first mode in the given lead
# kwant.plotter.map(fsys, numpy.absolute(psi[f]),
num_lead_cells=10,a=1,vmin=0,show=False)
jxtot=numpy.sum(jx,axis=0)
jytot=numpy.sum(jy,axis=0)
j2=numpy.sqrt(jxtot**2+jytot**2) ###add currents from all the
modes and then compute the square
kwant.plotter.map(fsys, j2, num_lead_cells=10,a=1,vmin=0,
vmax=0.1,show=False)
return jxtot , jytot
# If we want to plot the current density a bit of coarse graining helps.
from scipy.interpolate import griddata
from scipy.ndimage.filters import gaussian_filter
def smooth(x, y, jx, jy, L, W, sigma=0, downsample=1, dpu=3):
x_grid = numpy.linspace(-L-1., L, 2*L*dpu+1)
y_grid = numpy.linspace(-1, W+1., W*dpu+1)
grid = numpy.array([_.flatten() for _ in
numpy.meshgrid(x_grid,y_grid)]).T
jx_grid = griddata((x,y), jx, grid, method='nearest')
jx_grid = jx_grid.reshape(len(y_grid), len(x_grid))
dd = dpu*downsample
#x_gr = x_grid[::dd]
jx_gr = gaussian_filter(jx_grid, sigma*dpu)[::dd,::dd]
jy_grid = griddata((x,y), jy, grid, method='nearest')
jy_grid = jy_grid.reshape(len(y_grid), len(x_grid))
#y_gr = y_grid[::dd]
jy_gr = gaussian_filter(jy_grid, sigma*dpu)[::dd,::dd] ### takes
each third element
return jx_gr, jy_gr
# the (coarse-grained) current density looks good in a stream-plot:
def plot_stream(jx, jy, ttl='J', cmap='rainbow', vmax=0.03, Wl=40):
fig, ax = pyplot.subplots(figsize=(11, 7))
jval = numpy.sqrt(jx**2 + jy**2)
sy, sx = jx.shape
# xext =
# yext = (sy-1)
print (sx,sy)
xgr = numpy.linspace(xmin, xmax, sx)
ygr = numpy.linspace(0, Ly, sy)
# jval_th = numpy.tanh(jval/vmax)
ax.streamplot(xgr, ygr, jx, jy, density=2, arrowsize=1,
color = numpy.tanh(jx/vmax), cmap=cmap,
norm = colors.Normalize(-1,1),
linewidth=numpy.tanh(jval/vmax)*2)
ax.set_title(ttl)
ax.margins(0.05); ax.set_aspect(1);
# I use bars to indicate the position of the leads around the
scattering
# region, you need to adjust this to your geometry
## ax.broken_barh([(0, -5)] , (0, xext), edgecolor='grey',
facecolor='grey')
## ax.broken_barh([(yext, 5)] , (0, xext), edgecolor='grey',
facecolor='grey')
## ax.broken_barh([(0, Wl)] , (0, -5), edgecolor='grey',
facecolor='grey')
## ax.broken_barh([(yext, -Wl)] , (0, -5), edgecolor='grey',
facecolor='grey')
## ax.broken_barh([(0, yext)] , (0, xext), edgecolor='grey',
facecolor='None')
On 07/03/17 16:20, enrique colomes wrote:
Hello,
I was successful in introducing disorder in my Haldane (second nearest
neighbours) system.
I am using kwant1.3, I have read the new two threads [1] and [2], and I
have read [3]. However, I am still unable to compute and plot the local
density current for my system in the scattering region.
I would really appreciate any help. I attach the code.
Once I have it I will attach it, so people can actually use it and
solve future problems.
Regards,
Enrique
[1]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html>
[2]
https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html>
[3]
<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html>
https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current
from scipy.sparse import spdiags
from math import pi, sqrt, tanh, e
import numpy as np
import kwant
import random
import math
import tinyarray
import scipy.sparse.linalg as sla
from matplotlib import pyplot
# Define the graphene lattice
graphene = kwant.lattice.honeycomb()
a, b = graphene.sublattices
L=4
W=4
t=1
tt= 0.03 * t
phi=1*pi/2
def make_system(pot=0.000):
nnn_hoppinga=tt*e**(1j*phi)
nnn_hoppingb=tt*e**(1j*phi)
#### Define the scattering region. ####
def rec(pos):
x, y = pos
return 0 < x < L and 0 < y < W
sys = kwant.Builder()
# w: width and pot: potential maximum of the p-n junction
def potential(site):
(x, y) = site.pos
d = y * 2 + x * 3
return 0
sys[graphene.shape(rec, (1, 1))] = potential
# specify the hoppings of the graphene lattice
nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a))
nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] =
-t
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
sys[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb
#### Define the leads. ####
# left lead
sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
def lead0_shape(pos):
x, y = pos
return (0 < y < W)
lead0 = kwant.Builder(sym0)
lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb
# right lead
sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))
sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)])
sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
def lead1_shape(pos):
x, y = pos
return (0 < y < W)
lead1 = kwant.Builder(sym1)
lead1[graphene.shape(lead1_shape, (0, 0))] = pot
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]]
= -t
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_a]] = -nnn_hoppinga
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in
nnn_hoppings_b]] = -nnn_hoppingb
sys.attach_lead(lead1)
sys.attach_lead(lead0)
for site in sys.leads[0].interface+sys.leads[1].interface:
sys[site]=0
return sys
def jcurrent(sys):
J=kwant.operator.Current(sys)
wfs = kwant.wave_function(sys)
wf= wfs[0]
j=J(wf)
kwant.plotter.map(sys, j)
def main():
sys = make_system()
# Finalize the system.
sys = sys.finalized()
# Plot local current
jcurrent(sys)
# Call the main function if the script gets executed (as opposed to
imported).
# See<http://docs.python.org/library/__main__.html>.
if __name__ == '__main__':
main()
El 09/02/2017 a las 2:04, Joseph Weston escribió:
Hi again,
b) When I asked about computing it directly from Kwant it was because I read in another
thread that it was going to be in Kwant1.3, but I was not sure if it was already
avalaible. I will "fight" then with the code and will try to use your
guidelines Joseph. Definetively doing it with neighbors will be much faster than itering
the whole system. Will this new version be ready soon? :-)
Well it is already available in the 'development' version of Kwant. If
you're using 'conda' to install kwant you can easily install the
development version with:
conda install -c kwant kwant
otherwise you can follow the instructions on the website to install the
development version of Kwant "from source".
There are still a couple of features that we want to implement before
making the 1.3 release, but I would imagine that this release would
be before the end of the month.
Joe
