Mailman 3 - Kwant-discuss

Kwant fix for osx-arm64 conda
by Anton Akhmerov June 11, 2025

June 11, 2025

Hi all, I was finally able to identify and fix a long-standing issue for Kwant packaged via conda on apple silicon. The fix should work with the latest version of Kwant (1.5) and latest release of python-mumps. Happy Kwanting, Anton

2 1

Potential barrier in Bernal-stacked graphite
by cw14mi＠brocku.ca May 5, 2025

May 5, 2025

Hi all, My name is Connor. I'm a PhD student at Brock University. As an introduction to using Kwant, I want to simulate the transmission though a slab of Bernal-stacked graphene layers with a potential barrier (that is, all 4 basis atoms in a single unit cell, denoted {A1, B1, A2, B2}, have an on-site energy denoted g, while all other atoms have onsite energy zero. The first-neighbour hopping parameter (between {A1,B1} and {A2,B2}) is denoted t1, and the second-neighbour hopping parameter (between {B1,A2} in Bernal graphite) is denoted t2. I treat the barrier as my scattering region, and the rest of slab as the two leads. Periodic boundary conditions are applied along both in-plane primitive directions. My code (below) runs without error. However, when I plot my scattering region after attaching one lead (the one pointing in the positive stacking direction, taken to be the z-axis), the lead does not appear in the plot (i.e. in this Kwant tutorial: https://kwant-project.org/doc/dev/tutorial/first_steps, the lead appears as a faded red region, which is absent in my plot). Could someone tell me how to go about attaching the leads to my scattering region? Thank you for your time. Please let me know if I've left out any information, Connor ----- Code below ----- import numpy as np import kwant from matplotlib import pyplot # unit cell extrusion La = 1 Lb = 1 Lc = 1 # lattice parameters a0 = 1.42 c0 = 3.35 # hopping parameters / defect on-site energy t1 = 1.0 t2 = 0.5 g = 0.7 # lattice vectors a1 = a0 * np.array([1.5, 0.5*np.sqrt(3), 0.0]) a2 = a0 * np.array([1.5, -0.5*np.sqrt(3), 0.0]) a3 = c0 * np.array([0.0, 0.0, 2.0]) # basis atoms b1 = a0 * np.array([0.0, 0.0, 0.0]) + c0 * np.array([0.0, 0.0, 0.0]) #A1 b2 = a0 * np.array([1.0, 0.0, 0.0]) + c0 * np.array([0.0, 0.0, 0.0]) #B1 b3 = a0 * np.array([1.0, 0.0, 0.0]) + c0 * np.array([0.0, 0.0, 1.0]) #A2 b4 = a0 * np.array([2.0, 0.0, 0.0]) + c0 * np.array([0.0, 0.0, 1.0]) #B2 # Bernal graphite lattice lat = kwant.lattice.general([a1, a2, a3], basis=[b1, b2, b3, b4], norbs=1) A1, B1, A2, B2 = lat.sublattices # left lead lead = kwant.Builder(kwant.TranslationalSymmetry(lat.vec([La, 0, 0]), lat.vec([0, Lb, 0]), lat.vec([0, 0, 1]) )) lead[lat.neighbors(1)] = -t1 lead[lat.neighbors(2)] = -t2 for i in range(La): for j in range(Lb): lead[A1(i,j,0)] = 0 lead[B1(i,j,0)] = 0 lead[A2(i,j,0)] = 0 lead[B2(i,j,0)] = 0 lead = kwant.wraparound.wraparound(lead, keep=2) # scattering region sr = kwant.Builder(kwant.TranslationalSymmetry(lat.vec([La, 0, 0]), lat.vec([0, Lb, 0]))) sr[lat.neighbors(1)] = -t1 sr[lat.neighbors(2)] = -t2 for i in range(La): for j in range(Lb): for k in range(Lc): sr[A1(i,j,k)] = g sr[B1(i,j,k)] = g sr[A2(i,j,k)] = g sr[B2(i,j,k)] = g sr = kwant.wraparound.wraparound(sr) # attach leads sr.attach_lead(lead) kwant.plot(sr); exit()

1 0

Tkwant
by Tiago Silva May 1, 2025

May 1, 2025

Hello everyone, I'm klarc. I recently discovered Tkwant, and the part about obtaining the green lesser function of a system seemed interesting to me. Using the tutorial available in section 2.7 I increased the onsite energies of a potential using a Heaviside step. But for some strange reason I got an error. Here's the code: def onsite_C(site, time): (x,) = site.pos return (V/2) if time >= 0 else 0 def make_double_impurity_system(gamma_L, gamma_C, gamma_R, eps_L, eps_R): # system building lat = kwant.lattice.chain(a=1, norbs=1) syst = kwant.Builder() # central scattering region syst[(lat(x) for x in [-1, 2])] = 0 syst[lat(0)] = onsite_C syst[lat(1)] = onsite_C syst[(lat(0), lat(-1))] = -gamma_L syst[(lat(1), lat(0))] = -gamma_C syst[(lat(2), lat(1))] = -gamma_R # add leads sym1 = kwant.TranslationalSymmetry((-1,)) lead_left = kwant.Builder(sym1) lead_left[lat(0)] = onsite_L lead_left[lat.neighbors()] = -1 syst.attach_lead(lead_left) sym2 = kwant.TranslationalSymmetry((1,)) lead_right = kwant.Builder(sym2) lead_right[lat(2)] = onsite_R lead_right[lat.neighbors()] = -1 syst.attach_lead(lead_right) return syst, lat occupations = [tkwant.manybody.lead_occupation(chemical_potential=mu_L), tkwant.manybody.lead_occupation(chemical_potential=mu_R)] green = tkwant.manybody.GreenFunction(syst, max(times), occupations) idx = tkwant.system.siteId(syst, lat) site_0 = idx(0) site_1 = idx(2) green_lesser = [] for time in tqdm(times, desc='Calculando'): green.evolve(time, 0) green.refine_intervals() green_lesser.append(green.lesser(site_0, site_1)) and error: --------------------------------------------------------------------------- TypeError Traceback (most recent call last) Cell In[299], line 9 6 green.evolve(time, 0) 7 green.refine_intervals() ----> 9 green_lesser.append(green.lesser(site_0, site_1)) File ~/anaconda3/envs/tkwant_env/lib/python3.10/site-packages/tkwant/greenfunctions.py:378, in GreenFunction.lesser(self, i, j, root) 376 if self._green_lesser is None: 377 self._init_lesser() --> 378 self._do_jobs(self._jobs_lesser) 379 return self._green_lesser.evaluate(i, j, root) File ~/anaconda3/envs/tkwant_env/lib/python3.10/site-packages/tkwant/greenfunctions.py:477, in GreenFunction._do_jobs(self, joblist) 475 logger.debug('jobname={}, args={}, kwargs={}'.format(jobname, args, kwargs)) 476 job = getattr(self, jobname) --> 477 job(*args, **kwargs) 478 joblist.clear() File ~/anaconda3/envs/tkwant_env/lib/python3.10/site-packages/tkwant/greenfunctions.py:490, in GreenFunction._evolve_lesser(self, time0, time1, **kwargs) 488 def _evolve_lesser(self, time0, time1, **kwargs): 489 logger.debug('evolve G^<') --> 490 self._green_lesser.evolve(time0, time1) File ~/anaconda3/envs/tkwant_env/lib/python3.10/site-packages/tkwant/greenfunctions.py:53, in _GreenGeneric.evolve(self, time0, time1) ... File tkwant/onebody/kernels.pyx:361, in tkwant.onebody.kernels.Interpolation.eval() File tkwant/onebody/kernels.pyx:416, in tkwant.onebody.kernels.Interpolation._estimate_stepsize() TypeError: Cannot convert 'complex' with non-zero imaginary component to 'double' (this most likely comes from the '**' operator; use 'cython.cpow(True)' to return 'nan' instead of a complex number). Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

1 0

Current visualization in 3D nanowire
by Thomas Bredewoud April 28, 2025

April 28, 2025

Dear all, Visualizing current flow in 3D systems is challenging. Kwant provides built-in features for 2D systems, but does not support 3D current visualization. To still gain some insight into the current flow in my 3D nanowire system, I analyze the current on a single 2D facet of the nanowire. However, the resulting current distribution appears inhomogeneous, which is surprising given that the physical system is homogeneous. I have attached a minimal example (current_visualization.py) that demonstrates the issue. Has anyone encountered a similar situation? Could this be a visualization artifact, or might it indicate an issue with how the current is computed in 3D systems? Any suggestions or insights would be greatly appreciated. Best regards, Thomas

2 2

Issue with importing kwant in python
by Li Jingyi April 28, 2025

April 28, 2025

Dear Kwant dev/maintaining group: I was trying to use kwant for setting up a small simulation with many body interaction in tunneling junctions. However the first step I tried importing kwant with ‘import kwant’ I encountered the kernel crushing with the following message : And when i tried to test ```python -c "import kwant; print('Success’)``` it returns ```zsh: abort python -c "import kwant; print('Success’)”````` which suggestes mismatched libraries. I tried other ways such as build from sources but all kind of did not work very well and the kernel keep crushing, and trying to run I am currently using MacBook Apple M1 MacOS 15.3.2 (24D81), I need some further instructions and help on how to proceed to install kwant and use it on my system. Best Jingyi

1 1

Peierls phase convention
by Pericles Philippopoulos April 15, 2025

April 15, 2025

Hello everyone, I am new to KWANT but have been trying to implement a tight-binding model over a (2D) square lattice under the influence of a constant perpendicular magnetic field. This example is explored in section IV B of [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)]. My implementation is inspired by the KWANT tutorials. I start with the following function to set up the model: def make_system(a: float=1, eps0: float=4.0, hop: float=1.0, r: float=20): """ Creates a tight-binding model Args: a (float): Lattice constant eps0 (float): On-site energy hop (float): Hopping energy r (float): Radius of the quantum dot Returns: syst (kwant.Builder): The system lat (kwant.lattice): The lattice """ lat = kwant.lattice.square(a, norbs=1) syst = kwant.Builder() # Define the quantum dot def circle(pos): (x, y) = pos rsq = x ** 2 + y ** 2 return rsq <= r ** 2 def hopx(tosite, fromsite, peierls): phi = peierls(tosite, fromsite) return -hop * t * phi def hopy(tosite, fromsite, peierls): phi = peierls(tosite, fromsite) return -hop * t * phi syst[lat.shape(circle, (0, 0))] = eps0 * t # hoppings in x-direction syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx # hoppings in y-directions syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = hopy # Closed system, no leads return syst, lat And then I diagonalize the resulting Hamiltonian for different values of the magnetic field (in units of the flux quantum hc/e in cgs units): def sorted_eigs(ev): evals, evecs = ev # Sort eigenvalues in ascending order and get the corresponding indices sorted_indices = np.argsort(evals) # Use the sorted indices to sort evals sorted_evals = evals[sorted_indices] # Use the sorted indices to sort eigenvectors accordingly sorted_eigenvectors = evecs[:, sorted_indices] return sorted_evals, sorted_eigenvectors for iB in np.linspace(0, 0.03, 30): B = iB def B_syst(pos): return B peierls_syst = gauge(B_syst) params = dict(peierls=peierls_syst) ham_mat = syst.hamiltonian_submatrix(params=params, sparse=True) # Diaglonalize the Hamiltonian evals, evecs = sorted_eigs(sla.eigsh(ham_mat.tocsc(), k=10, sigma=0.05)) If I compare the results to Fig. 2 (a) of [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)], I get results that are those for a magnetic field B that is ~2 times as small. Looking at the source code, specifically line 791 of https://gitlab.kwant-project.org/kwant/kwant/blob/v1.5.0/kwant/physics/gaug… I find that the phases are written with a factor of pi in the exponent. However, looking at [Graf, M., & Vogl, P. Electromagnetic fields and dielectric response in empirical tight-binding theory. Physical Review B, 51(8), 4940 (1995)], if I express the magnetic field in units of the flux quantum / l^2 (as is suggested in the docstring at line 1005) I expect to get a factor of h / hbar = 2*pi in the exponent, consistent with the discrepancy. Is my reasoning correct? Is the KWANT definition for the flux quantum perhaps hc/(2e) (ie the superconducting flux quantum)? Any help would be greatly appreciated. Thanks.

1 0

F(ferromagnet) - S(superconductor) Junction Conductance in Kwant – Need Help with Spin-Resolved Conductance.
by Jabed Umar April 14, 2025

April 14, 2025

Hi all, --------- I'm simulating an FS (Ferromagnet–Superconductor) junction in Kwant (code attached below) to study conductance as a function of energy, and I'm observing the following: 1. Sub-gap conductance, which depends on Andreev reflection, reduces with increasing exchange field. 2. Supra-gap conductance remains the same as in an NS junction since the exchange field doesn't affect transport above the gap. 3. A divergence peak in conductance at E = Δ comes from the divergence in the density of states of the s-wave superconductor. These results are consistent and expected. ---------However, I would like to extract spin-resolved conductance, especially in the sub-gap regime, where spin-up and spin-down conductance should differ. This is because minority spins can always find a partner for Andreev's reflection, but majority spins may not. This leads to asymmetry in the Andreev reflection probabilities for different spin channels. ------------I'm using the S-matrix from Kwant to compute conductance. From what I understand: - Normal reflection corresponds to r_ee^{σσ'} - Andreev reflection corresponds to r_he^{σσ'} As I understand it, the issue is that Kwant does not explicitly label spin in the mode indices, which makes it difficult to separate the spin components in the conductance calculation. Has anyone worked on spin-resolved conductance in FS or NS junctions using Kwant? Any advice, tips, or example code would be really helpful. Thank you! ---------------------------------------------------------------------Code for FS junction import kwant import numpy as np from matplotlib import pyplot # Pauli matrices sigma_0 = np.identity(2) sigma_x = np.array([[0, 1], [1, 0]]) sigma_y = np.array([[0, -1j], [1j, 0]]) sigma_z = np.array([[1, 0], [0, -1]]) tau_0 = np.identity(2) tau_x = np.array([[0, 1], [1, 0]]) tau_y = np.array([[0, -1j], [1j, 0]]) tau_z = np.array([[1, 0], [0, -1]]) # Lattice and system parameters a = 1 # lattice constant W, L = 1, 100 # width and length barrier = 0 # barrier height barrierpos = (19, 20) # barrier x-position range mu = 0.5 Delta = 0.01 Deltapos = 20 # ferromagnetic region t = 0.25 h = 0.2 # exchange field strength (ferromagnetism) # Kronecker products for combined spin ⊗ particle-hole basis Delta_term = Delta * np.kron(sigma_x, tau_0) t_hop = -t * np.kron(sigma_z, tau_0) sz_tz = np.kron(sigma_z, tau_z) k = np.kron(sigma_z, tau_0) # Define lattice lat = kwant.lattice.square(norbs=4) # for particle-hole basis syst = kwant.Builder() # Ferromagnetic region for x in range(Deltapos): for y in range(W): syst[lat(x, y)] = (2 * t - mu) * k + h * sz_tz # Barrier region (still ferromagnetic) for x in range(barrierpos[0], barrierpos[1]): for y in range(W): syst[lat(x, y)] = (2 * t + barrier - mu) * k + h * sz_tz # Superconducting region for x in range(Deltapos, L): for y in range(W): syst[lat(x, y)] = (2 * t - mu) * k + Delta_term # Hoppings syst[lat.neighbors()] = t_hop # === Leads === # Ferromagnetic lead (left) sym_left = kwant.TranslationalSymmetry((-a, 0)) lead0 = kwant.Builder(sym_left, conservation_law=-np.kron(sigma_z, tau_0)) for y in range(W): lead0[lat(0, y)] = (2 * t - mu) * k + h * sz_tz lead0[lat.neighbors()] = t_hop # Superconducting lead (right) sym_right = kwant.TranslationalSymmetry((a, 0)) lead1 = kwant.Builder(sym_right) for y in range(W): lead1[lat(0, y)] = (2 * t - mu) * k + Delta_term lead1[lat.neighbors()] = t_hop # Attach leads and finalize syst.attach_lead(lead0) syst.attach_lead(lead1) syst = syst.finalized() # === Conductance plotting === def plot_conductance(syst, energies): data = [] for energy in energies: smatrix = kwant.smatrix(syst, energy) # G = N - R_ee + R_he data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] - smatrix.transmission((0, 0), (0, 0)) + smatrix.transmission((0, 1), (0, 0))) pyplot.figure() #pyplot.plot(energies, normalized_data) pyplot.plot(energies, data) pyplot.xlabel("Energy [t]") pyplot.ylabel("Normalized Conductance") pyplot.title("FS Junction Conductance") pyplot.grid(True) pyplot.show() energies = np.linspace(0, 0.1, 100) plot_conductance(syst, energies)

1 0

Quantum transport
by 西科大刘燕 April 2, 2025

April 2, 2025

Dear all, How to calculate the local current in a graphene nanoribbon at a specific energy, where the device region contains Rashba spin-orbit coupling (RSOC), using the Kwant package?Can you provide a simple example? Any help and guidance would be appreciated! Best regards, Rease

1 0

Density of graphene states
by 18682755968＠163.com March 25, 2025

March 25, 2025

Dear all, Using the Kwant package to calculate the spin-resolved density of states in a graphene nanoribbon system, where the central device region contains Rashba spin-orbit coupling (RSOC) in its Hamiltonian while the leads are free of spin-orbit coupling. We aim to obtain the total spin-up density of states distribution in the device region as the electron energy varies from 0 to 1.The attached files contain our program code. We found that when the wire contains a suspended key, we can run the result, but we don't know if the result is correct. When the wire does not have a hovering key, a strange value appears when executing, and an error is reported. We seek advice on the necessary modifications to correct these errors. For reference and convenience, I attach the code I'm using at the end of this message. Any help and guidance would be appreciated! Best regards, Rease ############################################ import numpy as np import tinyarray as ta import kwant import matplotlib.pyplot as plt sigma_0 = ta.array([[1, 0], [0, 1]]) sigma_x = ta.array([[0, 1], [1, 0]]) sigma_y = ta.array([[0, -1j], [1j, 0]]) sigma_z = ta.array([[1, 0], [0, -1]]) lat = kwant.lattice.honeycomb(norbs=2) def make_system(t=1.0, lambda_rsoc=0.1): syst = kwant.Builder() def onsite(site): return sigma_0 def hopping(site1, site2): dx, dy = site2.pos - site1.pos rsoc = 1j * lambda_rsoc * (sigma_x * dy - sigma_y * dx) return -t * sigma_0 + rsoc def rectangle(pos): x, y = pos return 0 <= x < 10 and 0.5 <= y < 10 syst[lat.shape(rectangle, (0, 0))] = onsite syst[lat.neighbors()] = hopping def lead_onsite(site): return sigma_0 def lead_hopping(site1, site2): return -t * sigma_0 sym_left = kwant.TranslationalSymmetry((-1, 0)) lead_left = kwant.Builder(sym_left, conservation_law=-sigma_z) lead_left[lat.shape(rectangle, (0, 0))] = lead_onsite lead_left[lat.neighbors()] = lead_hopping syst.attach_lead(lead_left) sym_right = kwant.TranslationalSymmetry((1, 0)) lead_right = kwant.Builder(sym_right, conservation_law=-sigma_z) lead_right[lat.shape(rectangle, (0, 0))] = lead_onsite lead_right[lat.neighbors()] = lead_hopping syst.attach_lead(lead_right) return syst syst = make_system().finalized() def spin_up_dos(syst): spin_up = 0.5 * (sigma_0 + sigma_z) return kwant.operator.Density(syst, spin_up, sum=True) energies = np.linspace(0, 1, 100) dos = [] for energy in energies: ldos = kwant.ldos(syst, energy) dos.append(np.sum(ldos[::2])) plt.figure(figsize=(8, 6)) plt.plot(energies, dos, label="Spin-Up DOS") plt.xlabel("Energy [t]") plt.ylabel("Density of States") plt.title("Spin-Up Density of States in Graphene Nanoribbon") plt.legend() plt.grid(True) plt.show()

1 0

Quantum transport
by 18682755968＠163.com March 19, 2025

March 19, 2025

Dear all, Using the Kwant package to calculate the spin-resolved conductance of a graphene nanoribbon, where the device region incorporates Rashba spin-orbit coupling (RSOC) in its Hamiltonian while the leads do not. The left lead is numbered 0 and the right lead is numbered 1. The goal is to obtain the conductance from spin-up electrons in the left lead to spin-down electrons in the right lead. How should the program be written? For reference and convenience, I attach the code I'm using at the end of this message, but the result was unsatisfactory. Any help and guidance would be appreciated! Best regards, Rease ############################################################## import kwant import numpy as np from matplotlib import pyplot from tinyarray import array sigma_0 = array([[1, 0], [0, 1]]) sigma_x = array([[0, 1], [1, 0]]) sigma_y = array([[0, -1j], [1j, 0]]) sigma_z = array([[1, 0], [0, -1]]) graphene = kwant.lattice.general([(1, 0), (0.5, np.sqrt(3) / 2)], [(0, 0), (0, 1 / np.sqrt(3))],norbs=2) a, b = graphene.sublattices def make_system(W=10, L=30, t=1.0, alpha=0.1): syst = kwant.Builder() def circle(pos): x, y = pos return abs(x)<=L and abs(y)<=W syst[graphene.shape(circle, (0, 0))] = (4 * t * sigma_0 + alpha * sigma_y) syst[graphene.neighbors()] = -t * sigma_0 def lead_shape(pos): return abs(pos[1]) < W lead = kwant.Builder(kwant.TranslationalSymmetry((1, 0))) lead[graphene.shape(lead_shape, (0, 0))] = 4 * t * sigma_0 lead[graphene.neighbors()] = -t * sigma_0 syst.attach_lead(lead) syst.attach_lead(lead.reversed()) return syst def plot_conductance(syst, energies): data = [ ] for energy in energies: smatrix = kwant.smatrix(syst, energy) data.append(smatrix.transmission((1, 1), (0, 0))) pyplot.figure() pyplot.plot(energies, data) pyplot.xlabel("energy [t]") pyplot.ylabel("conductance [e^2/h]") pyplot.show() def main(): syst = make_system() kwant.plot(syst) syst = syst.finalized() plot_conductance(syst, energies=[0.01 * i for i in range(100)]) if __name__ == "__main__": main()

2 2