Oops, forgot to attach dump_system.py. Sorry about that.
On Fri, Oct 29, 2021 at 2:37 PM Mani Chandra <mc0...@gmail.com> wrote:
> 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 <christoph.gr...@cea.fr>
> 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
>>
>
import numpy as np
import kwant
import params
import system
import pickle
print(f'Making system...', end='', flush=True)
sys = system.make_system()
print(f'done.', flush=True)
print(f'Dumping system...', end='', flush=True)
file_sys = open(f'system_{params.length}x{params.height}.obj', 'wb')
pickle.dump(sys, file_sys)
print(f'done.', flush=True)
