Dear Anton,
To confirm my observation, I remove the non-essential part of the script to
test the mumps solver and default solver. The system is a zigzag graphene
ribbon, with 200X200 unit cells, there are 4 atoms in one unit cell so there
are 1.6x10^5 atoms. The OS is window 10, kwant is 1.2.2 and python is 3.4.4.
the time for calculating 20 smatrix is around 30 minutes for both default
solver and mumps. Please find attached the script and the runtime. The first
runtime is for mumps. The second runtime is for default. The kwant was
installed with the package prepared by Gohlke.
For the sparse solver the time is longer, 40 minutes but there is no umfpack
installed in the window system. So I do not show the runtime as it may not be
accurate.
KS
On 16/5/2017, 5:51 PM, "[email protected] on behalf of Anton Akhmerov"
<[email protected] on behalf of [email protected]> wrote:
Dear KS,
> I have compared the performance of default, mumps solvers in windows and
> Ubuntu machines. Their performance is very similar without significant
gain
> when using mumps. Do I miss anything in the installation? For example, not
> installing the mumps. Similar performance is found for the sparse solver.
> The system I considered has 10^5 to 6X10^5 sites and machine has a quad
core
> i7 cpu with 16G RAM. My installation procedures are those found in the
> website for pre complled packages.
This is contrary to our experience. Can you please provide more details:
* How did you install Kwant in both cases?
* What specific script did you run?
* What were the execution times that you observed?
Thanks,
Anton
> KS
>
>
>
> Disclaimer: This email (including any attachments) is for the use of the
> intended recipient only and may contain confidential information and/or
> copyright material. If you are not the intended recipient, please notify
the
> sender immediately and delete this email and all copies from your system.
> Any unauthorized use, disclosure, reproduction, copying, distribution, or
> other form of unauthorized dissemination of the contents is expressly
> prohibited.
Disclaimer: This email (including any attachments) is for the use of the
intended recipient only and may contain confidential information and/or
copyright material. If you are not the intended recipient, please notify the
sender immediately and delete this email and all copies from your system. Any
unauthorized use, disclosure, reproduction, copying, distribution, or other
form of unauthorized dissemination of the contents is expressly prohibited.
from __future__ import division # so that 1/2 == 0.5, and not 0
# this is needed for my mac version of python
import matplotlib
import os
from time import ctime
matplotlib.use('TKAgg')
# the above two statement is needed for mac os
from math import pi, sqrt, tanh
from cmath import exp
from matplotlib import pyplot
from storefile import store
#HIDDEN_BEGIN_dwhx
#numpy is needed to handle matrix
import kwant
import kwant.solvers.sparse as ksparse
import kwant.solvers.mumps as kmumps
import numpy as np
from matplotlib.colors import LogNorm
#setup pdf for storing images
from matplotlib.backends.backend_pdf import PdfPages
###setup directory
##curdir=os.getcwd()
##os.chdir(curdir+'/test')
print (ksparse.uses_umfpack)
#define function for magnetic field phase shift
def hopphase_x(sitei, sitej, flux):
xi, yi = sitei.pos
xj, yj = sitej.pos
return exp(-0.5*1j*flux*(xi-xj)*(yi+yj))
def hopphase_y(sitei, sitej, flux):
xi, yi = sitei.pos
xj, yj = sitej.pos
return exp(0.5*1j* flux * (xi + xj) * (yi - yj))
#HIDDEN_END_dwhx
print ('start',ctime())
sys = kwant.Builder()
sys2=kwant.Builder()
#HIDDEN_END_goiq
# Here, we are only working with square lattices
#HIDDEN_BEGIN_suwo
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
#a = 1
# lat = kwant.lattice.honeycomb(a)
# this lattice is not a honeycomb lattice. it is like a rectangle
# this lattice gives correct zigzag and armchair edges
graphene= kwant.lattice.general([(sqrt(3), 0), (0, 1)],
[(0,1/2), (1/(2*sqrt(3)), 0),(sqrt(3)/2, 0), (2/sqrt(3),1/2)])
def rect(pos):
x, y= pos
lx=L*sqrt(3);ly=W*1;
return 0<=x<lx and 0<=y<ly
#Lstruct and Lshift are labels for file name indicating the structure and shift
pyplot.close("all")
Lstruct="_W40L40"
Lshift="no shift"
pp = PdfPages('ribbon'+Lstruct+Lshift+'.pdf')
t = 1.0
W = 200 # y dimension scatteribg region and horizontal lead
Vposit=0 #position of the horizontal lead
WY=200 # width of horizontal lead which can be different from W
L = 200 # x dimension scattering region
WL=200 #WL is width of vertical lead; it could be different from L
lpos=0 # position of vertical lead; the L and WL are different
lpos2=0
Lstub=0 #length of vertical stub which can have potential
WL1=2
WL2=7
flux=0.00
shiftcenter=0.00 #shift the center scattering region
shiftlead=-0.0
enpt=0.3 # shift the energy range.
a, b, c, d =graphene.sublattices
# Define the scattering region
sys2[graphene.shape(rect,(0.,0.))]=0.0*t
sys2[graphene.neighbors()]=-t
for i in range(L):
wtmp=0 #wtmp define the width of the central ribbon
if lpos+WL-1>=i>=lpos: wtmp=-Lstub
for j in range(wtmp,W):
# On-site Hamiltonian
# if i<23 :
# sys[a(i, j)] = 0. * t
# sys[b(i,j)]=0.*t
# sys[c(i,j)]=0.*t
# sys[d(i,j)]=0.*t
# else :
sys[a(i, j)] = (0.+i*0.0)* t+shiftcenter
sys[b(i,j)]=(0.+i*0.0)*t+shiftcenter
sys[c(i,j)]=(0.+i*0.0)*t+shiftcenter
sys[d(i,j)]=(0.+i*0.0)*t+shiftcenter
#new hopping
sys[a(i,j),b(i,j)]=-t*hopphase_x(a(i,j),b(i,j),flux)
sys[b(i,j),c(i,j)]=-t*hopphase_x(b(i,j),c(i,j),flux)
sys[c(i,j),d(i,j)]=-t*hopphase_x(c(i,j),d(i,j),flux)
if j<0: icheck=lpos
else : icheck=0
if i>icheck:
# (xa,ya)=a(i,j).pos
# (xd,yd)=d(i-1,j).pos
# sys[a(i,j), d(i-1,j)]=-t*exp(flux*1j*0.5*(xa-xd)*(ya+yd))
sys[a(i,j), d(i-1,j)]=-t*hopphase_x(a(i,j) ,d(i-1,j),flux)
if lpos+WL>=i>=lpos: jcheck=-Lstub
else : jcheck=0
if j>jcheck:
## (xa,ya)=a(i,j-1).pos
## (xb,yb)=b(i,j).pos
## sys[b(i,j), a(i,j-1)]=-t*exp(flux*1j*0.5*(xb-xa)*(ya+yb))
sys[b(i,j), a(i,j-1)]=-t*hopphase_x(b(i,j),a(i,j-1),flux)
## (xc,yc)=c(i,j-1).pos
## (xd,yd)=d(i,j).pos
## sys[c(i,j),d(i,j-1)]=-t*exp(flux*1j*0.5*(xc-xd)*(yc+yd))
sys[c(i,j), d(i,j-1)]=-t*hopphase_x(c(i,j),d(i,j-1),flux)
# vertical lead
sym2 = kwant.TranslationalSymmetry(graphene.vec((0,-1)))
lead2=kwant.Builder(sym2)
for i in range(lpos, WL+lpos):
lead2[a(i,0)]=0.*t
lead2[b(i,0)]=0.*t
lead2[c(i,0)]=0.*t
lead2[d(i,0)]=0.*t
# lead2[graphene.neighbors()]=-t
lead2[a(i,0),b(i,0)]=-t*hopphase_y(a(i,0),b(i,0),flux)
lead2[b(i,0),c(i,0)]=-t*hopphase_y(b(i,0),c(i,0),flux)
lead2[c(i,0),d(i,0)]=-t*hopphase_y(c(i,0),d(i,0),flux)
lead2[b(i,0),a(i,-1)]=-t*hopphase_y(b(i,0),a(i,-1),flux)
lead2[c(i,0),d(i,-1)]=-t*hopphase_y(c(i,0),d(i,-1),flux)
if i>lpos:
lead2[a(i,0),d(i-1,0)]=-t*hopphase_y(a(i,0),d(i-1,0),flux)
sym3=kwant.TranslationalSymmetry(graphene.vec((0,1)))
lead3=kwant.Builder(sym3)
for i in range(lpos2, WL+lpos2):
lead3[a(i,0)]=0.0*t
lead3[b(i,0)]=0*t
lead3[c(i,0)]=0.0*t
lead3[d(i,0)]=0*t
# lead1[graphene.neighbors()]=-t
lead3[a(i,0),b(i,0)]=-t*hopphase_y(a(i,0),b(i,0),flux)
lead3[b(i,0),c(i,0)]=-t*hopphase_y(b(i,0),c(i,0),flux)
lead3[c(i,0),d(i,0)]=-t*hopphase_y(c(i,0),d(i,0),flux)
lead3[b(i,0),a(i,-1)]=-t*hopphase_y(b(i,0),a(i,-1),flux)
lead3[c(i,0),d(i,-1)]=-t*hopphase_y(c(i,0),d(i,-1),flux)
if i>lpos2:
lead3[a(i,0),d(i-1,0)]=-t*hopphase_y(a(i,0),d(i-1,0),flux)
#attach the lead here
sys.attach_lead(lead2)
sys.attach_lead(lead3)
#sys2.attach_lead(lead2)
#sys2.attach_lead(lead3)
#plot the tight-binding system
fig_structure=kwant.plot(sys, show=False)
pp.savefig(fig_structure)
#---------------
# plot band structure of lead2
lead2plot = lead2.finalized()
kwant.plotter.bands(lead2plot, show=False)
pyplot.title("band structure of lead")
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
pp.savefig()
#pyplot.show()
#
# pyplot.show()
#-------------------
#fsys=sys2.finalized()
fsys = sys.finalized()
energies = []
data = [] #data array for lead 2
data_chck=[]
print ("start", ctime())
for ie in range(-650, 650,65):
energy = ie * 0.0005+0.00001+enpt
print (ie,ctime())
# print energy
# compute the scattering matrix at a given energy
smatrix = kwant.smatrix(fsys, energy)
# smatrix=ksparse.smatrix(fsys,energy)
# smatrix=kmumps.smatrix(fsys,energy)
# compute the transmission probability
energies.append(energy)
data_chck.append(smatrix.transmission(1, 0))
print ('end1',ctime())
#end of energy block
# Use matplotlib to write output
# We should see conductance steps
f=open('data','w')
#plot and store 2nd lead
pyplot.figure()
pyplot.plot(energies, data_chck)
pyplot.title("total conductance ")
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pp.savefig()
print("end",ctime())
pp.close()
#pyplot.show()
Python 3.4.4 (v3.4.4:737efcadf5a6, Dec 20 2015, 20:20:57) [MSC v.1600 64 bit
(AMD64)] on win32
Type "copyright", "credits" or "license()" for more information.
>>>
======== RESTART: C:\Users\K S Chan\Desktop\test\ribbon_testsolver.py ========
False
start Tue May 16 23:29:36 2017
start Tue May 16 23:33:47 2017
-650 Tue May 16 23:33:47 2017
-585 Tue May 16 23:34:00 2017
-520 Tue May 16 23:34:13 2017
-455 Tue May 16 23:34:23 2017
-390 Tue May 16 23:34:44 2017
-325 Tue May 16 23:35:39 2017
-260 Tue May 16 23:36:38 2017
-195 Tue May 16 23:37:14 2017
-130 Tue May 16 23:37:51 2017
-65 Tue May 16 23:38:32 2017
0 Tue May 16 23:39:16 2017
65 Tue May 16 23:40:36 2017
130 Tue May 16 23:41:58 2017
195 Tue May 16 23:42:57 2017
260 Tue May 16 23:44:03 2017
325 Tue May 16 23:45:16 2017
390 Tue May 16 23:47:05 2017
455 Tue May 16 23:48:21 2017
520 Tue May 16 23:49:44 2017
585 Tue May 16 23:51:47 2017
end1 Tue May 16 23:53:21 2017
end Tue May 16 23:53:21 2017
>>>
======== RESTART: C:\Users\K S Chan\Desktop\test\ribbon_testsolver.py ========
False
start Wed May 17 00:05:50 2017
start Wed May 17 00:08:08 2017
-650 Wed May 17 00:08:08 2017
-585 Wed May 17 00:08:23 2017
-520 Wed May 17 00:08:35 2017
-455 Wed May 17 00:08:51 2017
-390 Wed May 17 00:09:12 2017
-325 Wed May 17 00:10:13 2017
-260 Wed May 17 00:11:11 2017
-195 Wed May 17 00:11:44 2017
-130 Wed May 17 00:12:23 2017
-65 Wed May 17 00:13:05 2017
0 Wed May 17 00:13:53 2017
65 Wed May 17 00:15:05 2017
130 Wed May 17 00:16:19 2017
195 Wed May 17 00:17:21 2017
260 Wed May 17 00:18:26 2017
325 Wed May 17 00:19:35 2017
390 Wed May 17 00:21:16 2017
455 Wed May 17 00:22:36 2017
520 Wed May 17 00:24:01 2017
585 Wed May 17 00:25:59 2017
end1 Wed May 17 00:27:32 2017
end Wed May 17 00:27:33 2017
>>>
