Oops, forgot to attach dump_system.py. Sorry about that.
On Fri, Oct 29, 2021 at 2:37 PM Mani Chandra
Hi Christoph,
Attached are the relevant files. The total code is ~100 lines.
Copying the instructions from my previous email:
time python dump_system.py time python solve.py time python densities_and_currents.py
The parameters for the run are set in params.py. The numbers currently set will produce a 5000x5000 site model. The first script will dump the system to disk, the second will compute the scattering states, and the third will load the system and the scattering states and compute the densities and the currents.
------------------------------ To load the data, do:
data = np.load(f'densities_and_currents_d_c_{params.d_c:.4f}.npz') rho_left = data['rho_left'] rho_right = data['rho_right']
---------------------------------- To plot the raw density:
rho = (rho_right - rho_left).reshape([params.length, params.height]) norm = matplotlib.colors.TwoSlopeNorm(vmin=rho.min(), vcenter=0, vmax=rho.max())
plt.contourf(rho.T, 1000, cmap='bwr', norm=norm) plt.gca().set_aspect('equal') plt.colorbar()
------------------------------------- To plot the interpolated density:
rho_interpolated = rho_right_interpolated - rho_left_interpolated rho_interpolated = rho_interpolated/rho_interpolated.max()
norm = matplotlib.colors.TwoSlopeNorm(vmin=rho_interpolated.min(), vcenter=0, vmax=rho_interpolated.max())
plt.contourf(x_rho, y_rho, rho_interpolated[:, :, 0], 1000, cmap='bwr', norm=norm) plt.gca().set_aspect('equal') plt.colorbar()
----------------------------------- Finally, check this out:
from scipy.ndimage import gaussian_filter
rho = (rho_right - rho_left).reshape([params.length, params.height]) rho_filtered = gaussian_filter(rho, sigma=50) norm = matplotlib.colors.TwoSlopeNorm(vmin=rho_filtered.min(), vcenter=0, vmax=rho_filtered.max())
plt.contourf(rho_filtered.T, 1000, cmap='bwr', norm=norm) plt.gca().set_aspect('equal')
The above code gives a plot (attached) quite similar to the interpolated kwant density, but scipy's gaussian filter takes only 10 secs to blur out the 5000x5000 raw density data. Granted, it will probably only work with rectangular lattices. I wish I could do the same with currents, but I'm not sure how to get the currents into a format like length x height (as I did from the density above) because of their definition on the hoppings. Is there a function to obtain the current in a convenient matrix-like format, like the above?
Cheers, Mani
On Fri, Oct 29, 2021 at 1:27 AM Christoph Groth
wrote: Mani Chandra wrote:
- I added you to the code repo on github:
It would be preferable if you could prepare a small self-contained example Python script (rather not a notebook) that you could post to the list. Notably, this would allow others to follow along (now and in the future).
I imagine that that script could be just 20 lines or so: Make up a system with a similar graph (the values of the Hamiltonian do not matter), then make up some makeshift current/density (just a random array for example), and use this to run the interpolator.
Christoph