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
