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 <[email protected]> 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 <[email protected] > > 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 >
