[Kwant] Current density

[Harshad Sahasrabudhe]

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