Commit ab98335e by Joseph Weston

### use `syst` for Kwant systems instead of `sys`

`this fixes #35`
parent 10025e14
 ... ... @@ -9,18 +9,18 @@ from math import pi @@ -81,6 +82,50 @@ return sys return syst +def make_system_note1(a=1, t=1.0, W=10, r1=10, r2=20): + lat = kwant.lattice.square(a) + sys = kwant.Builder() + syst = kwant.Builder() + def ring(pos): + (x, y) = pos + rsq = x**2 + y**2 + return ( r1**2 < rsq < r2**2) + sys[lat.shape(ring, (0, 11))] = 4 * t + sys[lat.neighbors()] = -t + syst[lat.shape(ring, (0, 11))] = 4 * t + syst[lat.neighbors()] = -t + sym_lead0 = kwant.TranslationalSymmetry((-a, 0)) + lead0 = kwant.Builder(sym_lead0) + def lead_shape(pos): ... ... @@ -29,20 +29,20 @@ + lead0[lat.shape(lead_shape, (0, W))] = 4 * t + lead0[lat.neighbors()] = -t + lead1 = lead0.reversed() + sys.attach_lead(lead0) + sys.attach_lead(lead1) + return sys + syst.attach_lead(lead0) + syst.attach_lead(lead1) + return syst + + +def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20): + lat = kwant.lattice.square(a) + sys = kwant.Builder() + syst = kwant.Builder() + def ring(pos): + (x, y) = pos + rsq = x**2 + y**2 + return ( r1**2 < rsq < r2**2) + sys[lat.shape(ring, (0, 11))] = 4 * t + sys[lat.neighbors()] = -t + syst[lat.shape(ring, (0, 11))] = 4 * t + syst[lat.neighbors()] = -t + sym_lead0 = kwant.TranslationalSymmetry((-a, 0)) + lead0 = kwant.Builder(sym_lead0) + def lead_shape(pos): ... ... @@ -51,16 +51,16 @@ + lead0[lat.shape(lead_shape, (0, 0))] = 4 * t + lead0[lat.neighbors()] = -t + lead1 = lead0.reversed() + sys.attach_lead(lead0) + sys.attach_lead(lead1, lat(0, 0)) + return sys + syst.attach_lead(lead0) + syst.attach_lead(lead1, lat(0, 0)) + return syst + + def plot_conductance(sys, energy, fluxes): def plot_conductance(syst, energy, fluxes): # compute conductance @@ -90,18 +135,31 @@ smatrix = kwant.smatrix(sys, energy, args=[flux]) smatrix = kwant.smatrix(syst, energy, args=[flux]) data.append(smatrix.transmission(1, 0)) - pyplot.figure() ... ... @@ -84,30 +84,30 @@ def main(): sys = make_system() syst = make_system() # Check that the system looks as intended. - kwant.plot(sys) - kwant.plot(syst) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, file="ab_ring_sys." + extension, + kwant.plot(syst, file="ab_ring_syst." + extension, + fig_size=size, dpi=_defs.dpi) + # Finalize the system. sys = sys.finalized() syst = syst.finalized() @@ -111,6 +169,17 @@ for i in range(100)]) + # Finally, some plots needed for the notes + sys = make_system_note1() + syst = make_system_note1() + for extension in ('pdf', 'png'): + kwant.plot(sys, file="ab_ring_note1." + extension, + kwant.plot(syst, file="ab_ring_note1." + extension, + fig_size=size, dpi=_defs.dpi) + sys = make_system_note2() + syst = make_system_note2() + for extension in ('pdf', 'png'): + kwant.plot(sys, file="ab_ring_note2." + extension, + kwant.plot(syst, file="ab_ring_note2." + extension, + fig_size=size, dpi=_defs.dpi) + + ... ...
 ... ... @@ -31,29 +31,29 @@ + fig.savefig("closed_system_result." + extension, dpi=_defs.dpi) def plot_wave_function(sys): def plot_wave_function(syst): + size = (_defs.figwidth_in, _defs.figwidth_in) + # Calculate the wave functions in the system. ham_mat = sys.hamiltonian_submatrix(sparse=True) ham_mat = syst.hamiltonian_submatrix(sparse=True) evecs = sla.eigsh(ham_mat, k=20, which='SM')[1] # Plot the probability density of the 10th eigenmode. - kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, - kwant.plotter.map(syst, np.abs(evecs[:, 9])**2, - colorbar=False, oversampling=1) + for extension in ('pdf', 'png'): + kwant.plotter.map( + sys, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1, + syst, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1, + file="closed_system_eigenvector." + extension, + fig_size=size, dpi=_defs.dpi) def main(): sys = make_system() syst = make_system() - # Check that the system looks as intended. - kwant.plot(sys) - kwant.plot(syst) - # Finalize the system. sys = sys.finalized() syst = syst.finalized()
 ... ... @@ -9,7 +9,7 @@ import kwant @@ -96,22 +97,40 @@ smatrix = kwant.smatrix(sys, energy) smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(0, 1)) - pyplot.figure() ... ... @@ -62,27 +62,27 @@ return 0 if site.family == a else 1 - # Plot the closed system without leads. - kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False) - kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_color=family_colors, site_lw=0.1, colorbar=False, + file="graphene_sys1." + extension, + kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False, + file="graphene_syst1." + extension, + fig_size=size, dpi=_defs.dpi) # Compute some eigenvalues. compute_evs(sys.finalized()) compute_evs(syst.finalized()) @@ -133,9 +155,11 @@ for lead in leads: sys.attach_lead(lead) syst.attach_lead(lead) - # Then, plot the system with leads. - kwant.plot(sys, site_color=family_colors, site_lw=0.1, - kwant.plot(syst, site_color=family_colors, site_lw=0.1, - lead_site_lw=0, colorbar=False) + size = (_defs.figwidth_in, 0.9 * _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_color=family_colors, colorbar=False, site_lw=0.1, + file="graphene_sys2." + extension, + kwant.plot(syst, site_color=family_colors, colorbar=False, site_lw=0.1, + file="graphene_syst2." + extension, + fig_size=size, dpi=_defs.dpi, lead_site_lw=0) # Finalize the system. sys = sys.finalized() syst = syst.finalized()
 ... ... @@ -11,23 +11,23 @@ @@ -22,7 +23,7 @@ return x**2 + y**2 < r**2 sys = kwant.Builder() - sys[lat.shape(circle, (0, 0))] = 0 + sys[lat.shape(circle, (0,0))] = 0 sys[lat.neighbors()] = t sys.eradicate_dangling() syst = kwant.Builder() - syst[lat.shape(circle, (0, 0))] = 0 + syst[lat.shape(circle, (0,0))] = 0 syst[lat.neighbors()] = t syst.eradicate_dangling() if tp: @@ -32,9 +33,11 @@ def plot_system(sys): - kwant.plot(sys) def plot_system(syst): - kwant.plot(syst) - # the standard plot is ok, but not very intelligible. One can do - # better by playing wioth colors and linewidths + # standard plot - not very intelligible for this particular situation + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, file="plot_graphene_sys1." + extension, + kwant.plot(syst, file="plot_graphene_syst1." + extension, + fig_size=size, dpi=_defs.dpi) # use color and linewidths to get a better plot ... ... @@ -36,23 +36,23 @@ def hopping_lw(site1, site2): return 0.04 if site1.family == site2.family else 0.1 - kwant.plot(sys, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw) - kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_lw=0.1, site_color=family_color, + hop_lw=hopping_lw, file="plot_graphene_sys2." + extension, + kwant.plot(syst, site_lw=0.1, site_color=family_color, + hop_lw=hopping_lw, file="plot_graphene_syst2." + extension, + fig_size=size, dpi=_defs.dpi) def plot_data(sys, n): def plot_data(syst, n): @@ -58,7 +65,11 @@ # the usual - works great in general, looks just a bit crufty for # small systems - kwant.plotter.map(sys, wf, oversampling=10, cmap='gist_heat_r') - kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r') + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plotter.map(sys, wf, oversampling=10, cmap='gist_heat_r', + kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r', + file="plot_graphene_data1." + extension, + fig_size=size, dpi=_defs.dpi) ... ... @@ -60,13 +60,13 @@ def family_shape(i): @@ -68,15 +79,22 @@ def family_color(i): return 'black' if sys.site(i).family == a else 'white' return 'black' if syst.site(i).family == a else 'white' - kwant.plot(sys, site_color=wf, site_symbol=family_shape, - kwant.plot(syst, site_color=wf, site_symbol=family_shape, - site_size=0.5, hop_lw=0, cmap='gist_heat_r') + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_color=wf, site_symbol=family_shape, + kwant.plot(syst, site_color=wf, site_symbol=family_shape, + site_size=0.5, hop_lw=0, cmap='gist_heat_r', + file="plot_graphene_data2." + extension, + fig_size=size, dpi=_defs.dpi) ... ... @@ -75,11 +75,11 @@ def site_size(i): return 3 * wf[i] / wf.max() - kwant.plot(sys, site_size=site_size, site_color=(0, 0, 1, 0.3), - kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3), - hop_lw=0.1) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_size=site_size, site_color=(0,0,1,0.3), + kwant.plot(syst, site_size=site_size, site_color=(0,0,1,0.3), + hop_lw=0.1, file="plot_graphene_data3." + extension, + fig_size=size, dpi=_defs.dpi) ... ...
 ... ... @@ -10,27 +10,27 @@ @@ -33,7 +34,10 @@ # checking shapes: sys = make_cuboid() syst = make_cuboid() - kwant.plot(sys) - kwant.plot(syst) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, file="plot_zincblende_sys1." + extension, + kwant.plot(syst, file="plot_zincblende_syst1." + extension, + fig_size=size, dpi=_defs.dpi) # visualize the crystal structure better for a very small system sys = make_cuboid(a=1.5, b=1.5, c=1.5) syst = make_cuboid(a=1.5, b=1.5, c=1.5) @@ -41,8 +45,12 @@ def family_colors(site): return 'r' if site.family == a else 'g' - kwant.plot(sys, site_size=0.18, site_lw=0.01, hop_lw=0.05, - kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05, - site_color=family_colors) + size = (_defs.figwidth_in, _defs.figwidth_in) + for extension in ('pdf', 'png'): + kwant.plot(sys, site_size=0.18, site_lw=0.01, hop_lw=0.05, + kwant.plot(syst, site_size=0.18, site_lw=0.01, hop_lw=0.05, + site_color=family_colors, + file="plot_zincblende_sys2." + extension, + file="plot_zincblende_syst2." + extension, + fig_size=size, dpi=_defs.dpi) ... ...
 ... ... @@ -9,7 +9,7 @@ # For plotting @@ -55,19 +56,25 @@ smatrix = kwant.smatrix(sys, energy, args=[-welldepth]) smatrix = kwant.smatrix(syst, energy, args=[-welldepth]) data.append(smatrix.transmission(1, 0)) - pyplot.figure() ... ... @@ -33,11 +33,11 @@ def main(): sys = make_system() syst = make_system() - # Check that the system looks as intended. - kwant.plot(sys) - kwant.plot(syst) - # Finalize the system. sys = sys.finalized() syst = syst.finalized()
 ... ... @@ -9,17 +9,17 @@ import kwant @@ -69,7 +70,10 @@ sys.attach_lead(right_lead) syst.attach_lead(right_lead) # Plot it, to make sure it's OK -kwant.plot(sys) -kwant.plot(syst) +size = (_defs.figwidth_in, 0.3 * _defs.figwidth_in) +for extension in ('pdf', 'png'): + kwant.plot(sys, file="quantum_wire_sys." + extension, + kwant.plot(syst, file="quantum_wire_syst." + extension, + fig_size=size, dpi=_defs.dpi) # Finalize the system sys = sys.finalized() syst = syst.finalized() @@ -90,8 +94,13 @@ # Use matplotlib to write output ... ...
 ... ... @@ -9,7 +9,7 @@ # For plotting @@ -70,19 +71,24 @@ smatrix = kwant.smatrix(sys, energy) smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission(1, 0)) - pyplot.figure() ... ... @@ -32,11 +32,11 @@ def main(): sys = make_system() syst = make_system() - # Check that the system looks as intended. - kwant.plot(sys) - kwant.plot(syst) - # Finalize the system. sys = sys.finalized() syst = syst.finalized()
 ... ... @@ -30,11 +30,11 @@ def main(): sys = make_system() syst = make_system() - # Check that the system looks as intended. - kwant.plot(sys) - kwant.plot(syst) - # Finalize the system. sys = sys.finalized() syst = syst.finalized()
 ... ... @@ -31,11 +31,11 @@ The `~kwant.builder.Builder` method ``possible_hoppings`` has been rendered obsolete. Where previously one would have had :: for kind in lat.nearest: sys[sys.possible_hoppings(*kind)] = t syst[syst.possible_hoppings(*kind)] = t now it suffices to write :: sys[lat.neighbors()] = t syst[lat.neighbors()] = t This is possible because `~kwant.builder.Builder` now accepts *functions* as keys in addition to `~kwant.builder.Site` objects and tuples of them ... ... @@ -126,7 +126,7 @@ the following prototype:: then the values of ``t``, ``B`` and ``pot`` for which to solve the system can be passed to `~kwant.solvers.default.smatrix` like this:: kwant.smatrix(sys, energy, kwant.smatrix(syst, energy, args=(2., 3., 4.)) With many parameters it can be less error-prone to collect all of them into a ... ... @@ -147,7 +147,7 @@ collection could be a dictionary, or a class instance, for example:: params = SimpleNamespace(t=1, mu=2) for params.B in B_values: kwant.smatrix(sys, energy, args=[params]) kwant.smatrix(syst, energy, args=[params]) Arguments can be passed in an equivalent way to `~kwant.solvers.default.wave_function`, ... ...
 ... ... @@ -35,7 +35,7 @@ To support these applications, ``attach_lead`` now returns a list of all the sites that have been added to the system. Creating a buffer for disorder can be thus done as follows:: sys[sys.attach_lead(lead, add_cells=10)] = onsite syst[syst.attach_lead(lead, add_cells=10)] = onsite Note how we set the onsite Hamiltonians of the sites that have been added to the value used in the system. ... ...
 ... ... @@ -28,7 +28,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20): lat = kwant.lattice.square(a) sys = kwant.Builder() syst = kwant.Builder() #### Define the scattering region. #### # Now, we aim for a more complex shape, namely a ring (or annulus) ... ... @@ -40,8 +40,8 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20): # and add the corresponding lattice points using the `shape`-function #HIDDEN_BEGIN_lcak sys[lat.shape(ring, (0, r1 + 1))] = 4 * t sys[lat.neighbors()] = -t syst[lat.shape(ring, (0, r1 + 1))] = 4 * t syst[lat.neighbors()] = -t #HIDDEN_END_lcak # In order to introduce a flux through the ring, we introduce a phase on ... ... @@ -61,11 +61,11 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20): return iy0 < 0 and ix0 == 1 # ix1 == 0 then implied # Modify only those hopings in x-direction that cross the branch cut def hops_across_cut(sys): for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(sys): def hops_across_cut(syst): for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(syst): if crosses_branchcut(hop): yield hop sys[hops_across_cut] = hopping_phase syst[hops_across_cut] = hopping_phase #HIDDEN_END_lvkt #### Define the leads. #### ... ... @@ -84,20 +84,20 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20): #### Attach the leads and return the system. #### #HIDDEN_BEGIN_skbz sys.attach_lead(lead) sys.attach_lead(lead.reversed()) syst.attach_lead(lead) syst.attach_lead(lead.reversed()) #HIDDEN_END_skbz return sys return syst def plot_conductance(sys, energy, fluxes): def plot_conductance(syst, energy, fluxes): # compute conductance normalized_fluxes = [flux / (2 * pi) for flux in fluxes] data = [] for flux in fluxes: smatrix = kwant.smatrix(sys, energy, args=[flux]) smatrix = kwant.smatrix(syst, energy, args=[flux]) data.append(smatrix.transmission(1, 0)) pyplot.figure() ... ... @@ -108,16 +108,16 @@ def plot_conductance(sys, energy, fluxes): def main(): sys = make_system() syst = make_system() # Check that the system looks as intended. kwant.plot(sys) kwant.plot(syst) # Finalize the system. sys = sys.finalized() syst = syst.finalized() # We should see a conductance that is periodic with the flux quantum plot_conductance(sys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi plot_conductance(syst, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi for i in range(100)]) ... ...
 ... ... @@ -31,7 +31,7 @@ def make_system(a=1, t=1.0, r=10): #HIDDEN_BEGIN_qlyd lat = kwant.lattice.square(a) sys = kwant.Builder() syst = kwant.Builder() # Define the quantum dot def circle(pos): ... ... @@ -44,19 +44,19 @@ def make_system(a=1, t=1.0, r=10): y = site1.pos[1] return -t * exp(-1j * B * y) sys[lat.shape(circle, (0, 0))] = 4 * t syst[lat.shape(circle, (0, 0))] = 4 * t # hoppings in x-direction sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx # hoppings in y-directions sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t # It's a closed system for a change, so no leads return sys return syst #HIDDEN_END_qlyd #HIDDEN_BEGIN_yvri def plot_spectrum(sys, Bfields): def plot_spectrum(syst, Bfields): # In the following, we compute the spectrum of the quantum dot # using dense matrix methods. This works in this toy example, as ... ... @@ -66,7 +66,7 @@ def plot_spectrum(sys, Bfields): energies = [] for B in Bfields: # Obtain the Hamiltonian as a dense matrix ham_mat = sys.hamiltonian_submatrix(args=[B], sparse=True) ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True) # we only calculate the 15 lowest eigenvalues ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False) ... ... @@ -82,36 +82,36 @@ def plot_spectrum(sys, Bfields): #HIDDEN_BEGIN_wave def plot_wave_function(sys): def plot_wave_function(syst): # Calculate the wave functions in the system. ham_mat = sys.hamiltonian_submatrix(sparse=True) ham_mat = syst.hamiltonian_submatrix(sparse=True) evecs = sla.eigsh(ham_mat, k=20, which='SM')[1] # Plot the probability density of the 10th eigenmode. kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, kwant.plotter.map(syst, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1) #HIDDEN_END_wave def main(): sys = make_system() syst = make_system() # Check that the system looks as intended. kwant.plot(sys) kwant.plot(syst) # Finalize the system. sys = sys.finalized() syst = syst.finalized() # The following try-clause can be removed once SciPy 0.9 becomes uncommon. try: # We should observe energy levels that flow towards Landau # level energies with increasing magnetic field. plot_spectrum(sys, [iB * 0.002 for iB in range(100)]) plot_spectrum(syst, [iB * 0.002 for iB in range(100)]) # Plot an eigenmode of a circular dot. Here we create a larger system for # better spatial resolution. sys = make_system(r=30).finalized() plot_wave_function(sys) syst = make_system(r=30).finalized() plot_wave_function(syst) except ValueError as e: if e.message == "Input matrix is not real-valued.": print("The calculation of eigenvalues failed because of a bug in SciPy 0.9.") ... ...
 ... ... @@ -39,7 +39,7 @@ def make_system(r=10, w=2.0, pot=0.1): x, y = pos return x ** 2 + y ** 2 < r ** 2 sys = kwant.Builder() syst = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): ... ... @@ -47,7 +47,7 @@ def make_system(r=10, w=2.0, pot=0.1): d = y * cos_30 + x * sin_30 return pot * tanh(d / w) sys[graphene.shape(circle, (0, 0))] = potential