Re: [Kwant] Current density
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
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
>
Re: [Kwant] Current density
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
#!/usr/bin/env python3
import cmath
import numpy as np
import kwant
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}
def hopping(site_i, site_j):
# whatever -- just call a few functions
return cmath.exp(-1j * np.linalg.norm(site_i.pos - site_j.pos))
lat = kwant.lattice.square()
syst = kwant.Builder()
syst[(lat(i, j) for i in range(1000) for j in range(200))] = 4
syst[lat.neighbors()] = hopping
fsyst = syst.finalized()
psi = np.random.rand(len(fsyst.sites))
J = current(fsyst, psi)
print('current between (0, 0) and (1, 0):', J[lat(0, 0), lat(1, 0)])
signature.asc
Description: PGP signature
Re: [Kwant] Current density
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 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 > 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
>>
>>
Re: [Kwant] Current density
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 signature.asc Description: PGP signature
Re: [Kwant] Current density
Hi Joe, 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. Thanks, Harshad On Tue, Jan 10, 2017 at 4:18 AM, Joseph Weston wrote: > Hi Harshad, > > > 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? > > There are a number of things that could be done, however in the next > release of Kwant (1.3) there will be direct support for calculating > currents from wavefunctions. We hope to have this release before the end > of the > month. > > In the meantime you can use the bleeding-edge version > of Kwant available on the Kwant Gitlab [1] (although this will require > building the package as described in the "contribute" documentation > [2]). > > If you decide to go this route, could you post a complete example > script for your problem, so that I can test the performance of our > implementation of current calculation? In any case you can post > back in this thread with additional questions/updates. > > Thanks, > > Joe > > Links > - > [1]: https://gitlab.kwant-project.org/kwant/kwant > [2]: https://kwant-project.org/contribute >
Re: [Kwant] Current density
Hi Harshad, > 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? There are a number of things that could be done, however in the next release of Kwant (1.3) there will be direct support for calculating currents from wavefunctions. We hope to have this release before the end of the month. In the meantime you can use the bleeding-edge version of Kwant available on the Kwant Gitlab [1] (although this will require building the package as described in the "contribute" documentation [2]). If you decide to go this route, could you post a complete example script for your problem, so that I can test the performance of our implementation of current calculation? In any case you can post back in this thread with additional questions/updates. Thanks, Joe Links - [1]: https://gitlab.kwant-project.org/kwant/kwant [2]: https://kwant-project.org/contribute signature.asc Description: PGP signature
Re: [Kwant] Current density
Hi Joe, Thanks a lot. I did it on purpose. :)
Re: [Kwant] Current density
Hi again, I noticed in the example code you posted that you have some magnetic field in the scattering region of your system but none in the leads. I don't know anything about what precisely you may be wanting to simulate, but in general when there are abrupt changes in the Hamiltonian you'll generate some unphysical effects in your results. If you want to be able to put magnetic field in the leads as well then read on, else ignore me :). For your case of uniform magnetic field you can choose a gauge such that the vector potential, and hence Peierls phase, is invariant by translation in the symmetry direction of your leads (if you have leads which are not all parallel to each other things become more complicated). You can then set the hoppings in the leads in exactly the same way as in the central scattering region. It seems that the Landau gauge chosen in your example code is the correct one to use when you have leads in the x-direction (as you do). Regards, Joe signature.asc Description: OpenPGP digital signature
Re: [Kwant] Current density
Hi, You can calculate the particle current from site "j" to site "i" due to a state "psi": I_ij = -2 Im(psi_i* H_ij psi_j) Where "Im" is the imaginary part, "H_ij" is the Hamiltonian matrix element between sites "i" and "j" and "*" is complex conjugation. H_ij can be obtained from your finalized system, "sys", with `sys.hamiltonian(i, j, args=args)`. The wavefunctions "psi" can be obtained from `kwant.wave_function(sys, args)`. The following Python snippet should give you what you want:: lat = kwant.lattice.square() B = 1.0 # your magnetic field E = 1.0 # energy @ which you want the current sys = make_system() wf = kwant.wave_function(sys, energy=E, arg=(B,)) site_i = lat(4, 5) # whatever sites you want site_j = lat(4, 6) # get the indices of the sites in the wavefunction i = sys.sites.index(site_i) j = sys.sites.index(site_j) def current(psi): H_ij = sys.hamiltonian(i, j, args=(B,)) return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag I_ij = 0 # add the contributions from all open modes in all leads for l, lead in enumerate(sys.leads): for psi in wf(l): I_ij += current(psi) There are several caveats to the above: 1. It will not work if we have sites with more than 1 orbital as we have assumed that there is a 1-1 mapping between sites and elements in "psi". 2. No statistical physics has been included in the above -- the quantity calculated is *not* strictly the particle current, but the "particle current at energy E". To get the actual particle current you need to do the statistical physics right. This will involve integrating over energy and weighting the contributions from each lead by some distribution function (e.g. Fermi-Dirac) to get the full particle current. Depending on what kind of system you are studying you may be able to justify foregoing the energy integration on physical grounds, but this is a question of the physics you are looking at and is independent of kwant. Hope that helps, Joe P.S. I have not tested the above code snippet, and it's entirely possible that I've got a sign wrong, missed a few syntax errors etc. The point is just to give you the gist of how you would go about implementing such a calculation. signature.asc Description: OpenPGP digital signature
