Hi,

Thanks for the reply. I am running a simple QPC simulation in B field like the QHE example, but with a realistic potential (with compressible and incompressible strips). I will try using the current operator from the bleeding edge version.

I should say that there is not yet tutorial documentation for using these operators (it will be written before the release), and the API documentation is currently a bit broken [1]. If you have any questions once you get it working, post back here. Joe [1]: https://gitlab.kwant-project.org/kwant/kwant/issues/75

Hi Adel, Thanks a lot for the code snippet! I now understand why my code takes so long. I was wondering what kind of map Python uses for storing sites? The lookup is constant time for unordered or hash maps in C++, and the size of the map shouldn't matter too much. I use these kinds of maps in C++ all the time to store data mapped to points. The maps sometimes have 10-20 million entries and the code still runs fast. Thanks, Harshad On Tue, Jan 10, 2017 at 3:00 PM, Abbout Adel <abbout.adel@gmail.com> wrote:

Dear Harshad,

To complement Joseph's answer, I would like to come back to your code:

You are using a double loop in which you are calling 'sys.sites' many times. This takes too much time especially for large systems like yours (>200 000 sites).

Instead of doing this, you can use the product of numpy.arrays and use the fact that the elements of the wavefunction and the Hamiltonian are organized in the same way. (which is the same as sys.sites)

you can do something like :

#m is the mode number def Current(m,lead_nbr=0): current=2* array([Wf(lead_nbr)[m]]).T * tsys.hamiltonian_submatrix(args=[phi])* (Wf(lead_nbr)[m].conj()) return current.imag

A toy example is provided below.

I would like to recommend for you, since you are studying a Quantum Point Contact (QPC), to cut your potential and delete all the sites whose potential is larger than some value (3 times the Fermi energy for example). This makes your program faster and may prevent you from facing some memory problems and at the same time does not really change your results.

I Hope that this helps Adel

import kwant from numpy import * from matplotlib import pyplot

def make_system(a=1, t=1.0, W=150, L=150): lat = kwant.lattice.square(a)

sys = kwant.Builder() def hopping(sitei, sitej, phi): xi, yi = sitei.pos xj, yj = sitej.pos return -exp(-0.5j * phi * (xi - xj) * (yi + yj)) #### Define the scattering region. #### sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t sys[lat.neighbors()] = hopping

sys[(lat(90,i) for i in range(W))]=9 sys[(lat(90+i,W/2) for i in range(30))]=9

lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0))) lead[(lat(0, j) for j in range(W))] = 4 * t lead[lat.neighbors()] = hopping

sys.attach_lead(lead) sys.attach_lead(lead.reversed())

return sys

def plot_conductance(sys, energies,phi): # Compute conductance data = [] for energy in energies: smatrix = kwant.smatrix(sys, energy,args=[phi]) data.append(smatrix.transmission(1, 0))

pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel("conductance [e^2/h]") pyplot.show()

sys = make_system()

def color(site): if sys[site]>4: return 'g' else: return 'k' def size(site): if sys[site]>4: return 0.6 else: return 0.3

kwant.plot(sys, site_color=color,site_size=size) tsys = sys.finalized() sites=[site for site in tsys.sites] hoppings=[hop for hop in sys.hoppings()] E=2 phi=0.1 Wf=kwant.wave_function(tsys,E,args=[phi])

def Current(m,lead_nbr=0): result=2* array([Wf(lead_nbr)[m]]).T * tsys.hamiltonian_submatrix(args=[phi])* (Wf(lead_nbr)[m].conj()) return result.imag

#Calculating the current for mode "mode_number", for the wave coming from lead "lead_nbr" mode_number=1 I=Current(mode_number,lead_nbr=0)

#Plotting the result takes much more time than obtaining the results itself.

#the blue color is for the currents in the positive directions for x and y # the red color is for the opisit directions def Bond_current(site1,site2,phi=phi): i,j = sites.index(site1),sites.index(site2) return 6*abs(I[i,j])

def Current_color(site1,site2,phi=phi): i,j = sites.index(site1),sites.index(site2) if (site1.pos[0]>site2.pos[0] or site1.pos[1]>site2.pos[1]) and I[i,j]>0: return 'r' else: return 'b'

kwant.plot(sys,hop_lw=Bond_current,site_color='w',hop_color=Current_color)

pyplot.show()

On Tue, Jan 10, 2017 at 7:11 AM, Harshad Sahasrabudhe <hsahasra@purdue.edu

...

wrote:

...

Hi All,

I am trying to calculate current density from the wavefunctions using the following code:

for i in range(Np): for j in range(Nc): lat_idx_i = i-floor(Np/2) lat_idx_j = j-floor(Nc/2) site_i = lat(lat_idx_i, lat_idx_j) idx_i = sys.sites.index(site_i)

if i < Np-1: site_j = lat(lat_idx_i+1, lat_idx_j) idx_j = sys.sites.index(site_j) H_ij = sys.hamiltonian(idx_i, idx_j, V, peierls_phase_factor) current_density_bond_x[i][j] = -2 * (wf_orb[idx_i].conjugate() \ * H_ij * wf_orb[idx_j]).imag

if j < Nc-1: site_j = lat(lat_idx_i, lat_idx_j+1) idx_j = sys.sites.index(site_j) H_ij = sys.hamiltonian(idx_i, idx_j, V, peierls_phase_factor) current_density_bond_y[i][j] = -2 * (wf_orb[idx_i].conjugate() \ * H_ij * wf_orb[idx_j]).imag

However, this code snippet takes about one and a half hours to run. The total number of sites in the system is about 201000. Is there any other way to write the code so that it runs faster?

THanks, Harshad

-- Abbout Adel

Hi Joe,

You can also just iterate over the system's graph to get all the hoppings. Something like:

def current(syst, psi, args=()): def hopping_current(i, j): H_ij = syst.hamiltonian(i, j, *args) return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag ## returns a dictionary that maps hopping -> current return {(syst.sites[i], syst.sites[j]): hopping_current(i, j) for i, j in syst.graph}

Thanks! This is faster than other methods. I need to sharpen my Python, got confused between list and dict. Thanks, Harshad On Wed, Jan 11, 2017 at 12:44 PM, Joseph Weston <joseph.weston08@gmail.com> wrote:

Hi again!

...

...

#m is the mode number def Current(m,lead_nbr=0): current=2* array([Wf(lead_nbr)[m]]).T * tsys.hamiltonian_submatrix( args=[phi])* (Wf(lead_nbr)[m].conj()) return current.imag

Thanks a lot for the code snippet! I now understand why my code takes so long.

The only thing I would say about this is that it is using dense linear algebra, so the complexity is O(N**2), also trying to call `hamiltonian_submatrix` will most likely blow up your memory, or not far off (10^5 sites means 10^10 matrix entries).

You can also just iterate over the system's graph to get all the hoppings. Something like:

def current(syst, psi, args=()):

def hopping_current(i, j): H_ij = syst.hamiltonian(i, j, *args) return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag

## returns a dictionary that maps hopping -> current return {(syst.sites[i], syst.sites[j]): hopping_current(i, j) for i, j in syst.graph}

This won't work as-is if you have >1 degree of freedom per site (matrix onsites / hoppings), but from your example it seems that you have 1 degree of freedom per site anyway. There is a full example attached, FYI.

...

I was wondering what kind of map Python uses for storing sites? The lookup is constant time for unordered or hash maps in C++, and the size of the map shouldn't matter too much. I use these kinds of maps in C++ all the time to store data mapped to points. The maps sometimes have 10-20 million entries and the code still runs fast.

I did not previously see that you were calling `sys.sites.index` in your loop. `sys.sites` is just a python list, so `sys.sites.index` is O(N). It only recently became apparent that having the inverse mapping (site -> index) would be useful. In bleeding edge kwant this is available as the `id_by_site` attribute of finalized systems. `id_by_site` is a python dictionary, which has the same complexity as a C++ hash map (i.e. O(1) to fetch an element)

Thanks,

Joe