Kwant-discuss

kwant-discuss@python.org

4 participants
933 discussions

Resolution of Tkwant
by Clarke, Isobel 05 Dec '22

05 Dec '22

Dear Kwant developers, I was wondering if there is a limit on the resolution of Kwant and Tkwant when calculating the density of an evolved wavefunction. I am currently working on a code to simulate the response of a carbon nanotube spin qubit to an electric dipole spin resonance drive. I am using Kwant to define the carbon nanotube as a scattering region and then apply a Hamiltonian to the system. The lowest eigenvector of the Hamiltonian at t=0 is determined and associated with a tkwant.onebody.WaveFunction. This wavefunction is evolved over time and the density of the wavefunction calculated using kwant.operator.Density. The relevant part of my code is given below: class System: def __init__(self, hamiltonian, pertubation_type, number_of_lattices, potential_type): self.potential_type = potential_type self.hamiltonian = hamiltonian self.number_of_lattices = number_of_lattices self.pertubation_type = pertubation_type self.evolve_state = False def potential(self, z, time): if self.potential_type == 0: #infinite square well total_potential = 0 elif self.potential_type == 1: #parabolic potential total_potential = .5 * (z * self.omega_0_au) ** 2 if self.pertubation_type == "cos": self.pertubation = self.cosine_v_ac else: self.pertubation = self.sine_v_ac total_potential += self.pertubation(time, z, self.eV_0_au, self.pulse_frequency_au, self.total_length_au) return total_potential def kwant_shape(self, site): (z,) = site.pos return (-self.total_length_au / 2 <= z < self.total_length_au / 2) def make_system(self): self.template = kwant.continuum.discretize(self.hamiltonian, grid=self.lattice_size_au) self.syst = kwant.Builder() # add the nanotube to the system self.syst.fill(self.template, self.kwant_shape, (0,)) self.syst = self.syst.finalized() def B_constant(z): return -self.g * self.mu_B_au * self.B_z_au * self.hbar_au / 2 def C_constant(z): return -self.g * self.mu_B_au * self.b_sl_au * z * self.hbar_au / 2 def D_constant(z): return -self.g * self.mu_B_au * self.B_y_au * self.hbar_au / 2 # coefficient for the kinetic energy term A_constant = self.hbar_au ** 2 / (2 * self.m_au) # parameters of hamiltonian self.params = dict(A=A_constant, V=self.potential, B=B_constant, C=C_constant, D=D_constant) self.tparams = self.params.copy() # copy the params array self.tparams['time'] = 0 # initial time = 0 # compute the Hamiltonian matrix hamiltonian = self.syst.hamiltonian_submatrix(params=self.tparams) # compute the eigenvalues (energies) and eigenvectors (wavefunctions). eigenValues, eigenVectors = np.linalg.eig(hamiltonian) # sort the eigenvectors and eigenvalues according the ascending eigenvalues. idx = eigenValues.argsort() self.initial_eigenvalues = eigenValues[idx] eigenVectors = eigenVectors[:, idx] # initial wave functions unperturbed by time dependent part of hamiltonian self.psi_1_init = eigenVectors[:, 0] self.psi_2_init = eigenVectors[:, 1] # an object representing the spin-up state self.spin_up_state = tkwant.onebody.WaveFunction.from_kwant(syst=self.syst, psi_init=self.psi_1_init, energy=eigenValues[0], params=self.params) return self.syst def evolve(self, time_steps): self.evolve_state = True total_osc_time = self.second_to_au(0.5e-8) times_au = np.linspace(0, total_osc_time, num=time_steps) times_SI = np.linspace(0, self.au_to_second(total_osc_time), num=time_steps) density_operator = kwant.operator.Density(self.syst) psi = self.spin_up_state data = {'states': []} for time in times_au: psi.evolve(time) # evolve the wavefunction density = psi.evaluate(density_operator) # compute density data['states'].append({'pdf': density.tolist()}) json_string = json.dumps(data) with open(r'\{}.json'.format(timestr), 'w') as outfile: outfile.write(json_string) self.evolve_state = False return True def main(): lattices = 100 potential = 0 system = System("((A * k_z**2) + V(z, time)) * identity(2) + B(z) * sigma_z + C(z) * sigma_x + D(z) * sigma_y", pertubation_type="sin", number_of_lattices=lattices, potential_type=potential) system.make_system() system.evolve(10) if __name__ == '__main__': main() Please note that I haven't included some conversion functions in the above code. Using this code, I obtain the results shown in the animation attached. The first subplot shows the density of the wavefunction across the carbon nanotube evolving with time. The animation shows an expected smooth Gaussian response for the first few time steps, however the graph then becomes very noisy. I was wondering why this occurs and if it could be due to a limit on the resolution of Tkwant? The second subplot shows the energy of the time-dependent perturbation over time, displaying a sin wave perturbation as a simple 'see-sawing' potential. In addition, I was wondering if you had any advice for making the code run faster? I ran the simulation for 10 time steps over a very short time period to obtain the attached animation and it took around 4 hours on my desktop. Thank you for any help with this. Best wishes, Isobel

4 3

Re: Resolution of Tkwant
by Clarke, Isobel 05 Dec '22

05 Dec '22

Dear Christoph, Thank you so much for your email and for your advice. I have attached the full code for my carbon nanotube system (simplified slightly) which, when run, should output the animation I have also attached. The system is evolved for 5 time steps over a short period of time (0 to 0.5e-8s). It seems that when I increase the number of lattice points, the probability density function becomes 'spiky' over time and the resolution of the graph deteriorates. However, this might just be a problem with my code. Thank you again for any help with this. Best wishes, Isobel ________________________________ From: Christoph Groth Sent: Thursday, December 01, 2022 15:34 To: Clarke, Isobel Cc: kwant-discuss(a)python.org; Buitelaar, Mark Subject: Re: [Kwant] Resolution of Tkwant Clarke, Isobel wrote: > I was wondering if there is a limit on the resolution of Kwant and > Tkwant when calculating the density of an evolved wavefunction. Hi, I’m not aware of any particular issue where a Tkwant simulation “explodes” after some time. Obviously, the numerical calculations have a limited precision. It would be most helpful if you could provide a complete example that demonstrates the problem that you experience. It doesn’t have to be the same physical system that you use in research. It’s actually better if it’s simpler, but one should be able to observe the concrete problems that trouble you. > In addition, I was wondering if you had any advice for making the code > run faster? I ran the simulation for 10 time steps over a very short > time period to obtain the attached animation and it took around > 4 hours on my desktop. As a vague rule of thumb, one needs to step up to a a computing cluster for tkwant when a laptop/desktop is sufficient with kwant alone. (This assumes that one wants to treat similar systems with comparable numbers of orbitals.) Or one has to wait longer... Tkwant calculations need to be integrated over many energies, while with Kwant at zero temperature that’s not necessary. So there’s literally more to do. There are bottlenecks and problems that could be improved in both Kwant and Tkwant, and there may be also possible workarounds. It can help to profile a run to see where time is actually spent. Having a small but representative running example (that ideally does not take four hours!) would help here as well. Cheers Christoph

1 0

Resolution of Tkwant
by zcapicc＠ucl.ac.uk 30 Nov '22

30 Nov '22

Dear Kwant developers, I was wondering if there is a limit on the resolution of Kwant and Tkwant when calculating the density of an evolved wavefunction. I am currently working on a code to simulate the response of a carbon nanotube spin qubit to an electric dipole spin resonance drive. I am using Kwant to define the carbon nanotube as a scattering region and then apply a Hamiltonian to the system. The lowest eigenvector of the Hamiltonian at t=0 is determined and associated with a tkwant.onebody.WaveFunction. This wavefunction is evolved over time and the density of the wavefunction calculated using kwant.operator.Density. The relevant part of my code is given below: class System: def __init__(self, hamiltonian, pertubation_type, number_of_lattices, potential_type): self.potential_type = potential_type self.hamiltonian = hamiltonian self.number_of_lattices = number_of_lattices self.pertubation_type = pertubation_type self.evolve_state = False def potential(self, z, time): if self.potential_type == 0: #infinite square well total_potential = 0 elif self.potential_type == 1: #parabolic potential total_potential = .5 * (z * self.omega_0_au) ** 2 if self.pertubation_type == "cos": self.pertubation = self.cosine_v_ac else: self.pertubation = self.sine_v_ac total_potential += self.pertubation(time, z, self.eV_0_au, self.pulse_frequency_au, self.total_length_au) return total_potential def kwant_shape(self, site): (z,) = site.pos return (-self.total_length_au / 2 <= z < self.total_length_au / 2) def make_system(self): self.template = kwant.continuum.discretize(self.hamiltonian, grid=self.lattice_size_au) self.syst = kwant.Builder() # add the nanotube to the system self.syst.fill(self.template, self.kwant_shape, (0,)) self.syst = self.syst.finalized() def B_constant(z): return -self.g * self.mu_B_au * self.B_z_au * self.hbar_au / 2 def C_constant(z): return -self.g * self.mu_B_au * self.b_sl_au * z * self.hbar_au / 2 def D_constant(z): return -self.g * self.mu_B_au * self.B_y_au * self.hbar_au / 2 # coefficient for the kinetic energy term A_constant = self.hbar_au ** 2 / (2 * self.m_au) # parameters of hamiltonian self.params = dict(A=A_constant, V=self.potential, B=B_constant, C=C_constant, D=D_constant) self.tparams = self.params.copy() # copy the params array self.tparams['time'] = 0 # initial time = 0 # compute the Hamiltonian matrix hamiltonian = self.syst.hamiltonian_submatrix(params=self.tparams) # compute the eigenvalues (energies) and eigenvectors (wavefunctions). eigenValues, eigenVectors = np.linalg.eig(hamiltonian) # sort the eigenvectors and eigenvalues according the ascending eigenvalues. idx = eigenValues.argsort() self.initial_eigenvalues = eigenValues[idx] eigenVectors = eigenVectors[:, idx] # initial wave functions unperturbed by time dependent part of hamiltonian self.psi_1_init = eigenVectors[:, 0] self.psi_2_init = eigenVectors[:, 1] # an object representing the spin-up state self.spin_up_state = tkwant.onebody.WaveFunction.from_kwant(syst=self.syst, psi_init=self.psi_1_init, energy=eigenValues[0], params=self.params) return self.syst def evolve(self, time_steps): self.evolve_state = True total_osc_time = self.second_to_au(0.5e-8) times_au = np.linspace(0, total_osc_time, num=time_steps) times_SI = np.linspace(0, self.au_to_second(total_osc_time), num=time_steps) density_operator = kwant.operator.Density(self.syst) psi = self.spin_up_state data = {'states': []} for time in times_au: psi.evolve(time) # evolve the wavefunction density = psi.evaluate(density_operator) # compute density data['states'].append({'pdf': density.tolist()}) json_string = json.dumps(data) with open(r'\{}.json'.format(timestr), 'w') as outfile: outfile.write(json_string) self.evolve_state = False return True def main(): lattices = 100 potential = 0 system = System("((A * k_z**2) + V(z, time)) * identity(2) + B(z) * sigma_z + C(z) * sigma_x + D(z) * sigma_y", pertubation_type="sin", number_of_lattices=lattices, potential_type=potential) system.make_system() system.evolve(10) if __name__ == '__main__': main() Please note that I haven't included some conversion functions in the above code. Using the results from this code, I can obtain a graph of density of the wavefunction across the carbon nanotube evolving with time. The graph shows an expected smooth Gaussian response for the first few time steps, however the graph then becomes very noisy and spiky. I was wondering why this occurs and if it could be due to a limit on the resolution of Tkwant? In addition, I was wondering if you had any advice for making the code run faster? I ran the simulation for 10 time steps over a very short time period to obtain the attached animation and it took around 4 hours on my desktop. Thank you for any help with this. Best wishes, Isobel

1 0

About calculating local observables using Green function
by Hari Gautam 25 Nov '22

25 Nov '22

Dear kwant users, I find the wavefunction methods with sufficient example for calculation of local observables using operators in kwant. I would also like to calculate the local observables using Green functions. As far as I know, I understood that the green_function method only provides the surface Green function at the interface of the lead and the scattering system which is helpful to calculate the total conductance. How can we define local observables using Green functions in kwant? Is there any methods like operators included in the kwant? Thank you Hari

2 3

Attaching lead to a system
by Shashank Kumar 03 Nov '22

03 Nov '22

Hi all, I have made a one dimensional system with 3 sub-lattices. However, in order to attach leads I have to make a lead and manually attach it to sublattice A. But, apparently this code below is not working. Any suggestion will be help full. def make_system(nlattice): lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (1/3.0, 0), (2/3.0, 0)]) sub_a, sub_b, sub_c = lat.sublattices # central scattering region syst = kwant.Builder() syst[(sub_a(x,0) for x in range(nlattice))] = 4 syst[(sub_b(x,0) for x in range(nlattice))] = 4 syst[(sub_c(x,0) for x in range(nlattice))] = 4 syst[lat.neighbors()] = -1 # add leads lat_lead = kwant.lattice.square() sym_lead1 = kwant.TranslationalSymmetry((-1, 0)) lead1 = kwant.Builder(sym_lead1) lead1[lat_lead(-1, 0)] = 4 lead1[lat_lead.neighbors()] = -1 syst[lat_lead(-1,0)] = 4 syst[lat_lead(-1,0), sub_a(0,0)] = -1 syst.attach_lead(lead1) syst.attach_lead(lead1.reversed()) kwant.plot(lead1) return syst fsyst = make_system(nlattice=10).finalized() kwant.plot(sys, fig_size=[10,8], site_size=0.2) Shashank Kumar

2 1

Problem in Bandstructure of AGNR defect
by sahu.ajit92＠gmail.com 15 Oct '22

15 Oct '22

Dear all, I make AGNR in Kwant and introduce the defect in both the scattering region and lead. When plotting the bandstructure of the lead there is a straight line that comes at E=0. I am confused that there is some error in the code or not. Please help Thanking You code- ________________ import kwant import numpy as np import matplotlib.pyplot as plt pi = np.pi a = 1 d = 1.42; a1 = d*np.sqrt(3); lat = kwant.lattice.general( [[a1, 0], [a1/2, a1*np.sqrt(3)/2]],[[0, 0], [0, a1/np.sqrt(3)]] ) a,b = lat.sublattices def set_pub(): # for plotting the figures plt.rcParams.update({ "font.weight": "bold", # bold fonts "font.size":18, #"tick.labelsize": 15, # large tick labels "lines.linewidth": 2, # thick lines "lines.color": "k", # black lines "grid.color": "0.5", # gray gridlines "grid.linestyle": "-", # solid gridlines "grid.linewidth": 0.5, # thin gridlines "xtick.major.size":10, "xtick.minor.size":5, "ytick.major.size":10, "ytick.minor.size":5, "savefig.dpi": 600, # higher resolution output. }) def make_sys(w,l): def rect(pos): x,y = pos return -w*a1<x<w*a1 and -l*a1<y<l*a1 model1 = kwant.Builder() model1[lat.shape(rect,(1,1))] = 0 model1[lat.neighbors()] = 2.46 return model1 def lead_join(model1,w): def lead_shape(pos): x, y = pos return -w*a1<x<w*a1 Num = 8*(a1*np.sqrt(3)/2) sym = kwant.TranslationalSymmetry([0,Num]) lead = kwant.Builder(sym) lead[lat.shape(lead_shape, (1,1))] = 0 lead[lat.neighbors()] = 2.46 model1.attach_lead(lead) model1.attach_lead(lead.reversed()) return model1,lead def defect(model1,lead): def delet_lead(pos, shift = +a1*np.sqrt(3)): x, y = pos z=x**2+(y-shift)**2 return z<(a1)**2 def delet_str(pos, shift = a1*np.sqrt(3)): x, y = pos z = x**2 + (y-shift)**2 return z<(a1)**2#z<(2*a1)**2 del model1[lat.shape(delet_str, (1,1))] del lead[lat.shape(delet_lead, (2,2))] kwant.plot(model1) set_pub() plt.rc('axes', linewidth = 2.0) plt.show() return model1 def plot_bandstructure(flead, momenta): bands = kwant.physics.Bands(flead) energies = [bands(k) for k in momenta] band_gap = 100 for k in momenta: for band in bands(k): if(band <0): #Lower than fermi band min_band=band else: max_band=band if(max_band-min_band<band_gap): band_gap=max_band-min_band break print("Gap Energy= " + str(band_gap) +" ev ") kwant.plotter.bands(flead) # bandstructure of lead plt.ylim([-3,3]) set_pub() plt.xlabel("K $(\AA^{-1})$",fontweight='bold') plt.ylabel("Energy (eV)",fontweight='bold') plt.show() return band_gap def main(): w = 5 l = 8 sys = make_sys(w,l) syst,lead = lead_join(sys,w) syst = defect(syst,lead) fsys = syst.finalized() momenta = [-pi + 0.02 * pi * i for i in range(101)] Band_Gap = plot_bandstructure(lead.finalized(), momenta) if __name__ == '__main__': main()

2 2

Re: Problem in Bandstructure of AGNR defect
by Ajit Kumar Sahu 15 Oct '22

15 Oct '22

1 0

Does Kwant use Transfer Matrix like Methods or ED?
by Kartikeya Arora 27 Sep '22

27 Sep '22

I tried to read the documentations and everything but got a bit confused. Here's my question and the task I want to complete. I want to simulate a 1d Anderson Model and calculate either Localization Length or Conductance for different system sizes. Now, using normal ED, I can only go up to 5000 system size in 24 hours (100 disorder realizations). Now, I want to calculate quantities like Localization Length/Conductance as a function of system size. At your best estimate, How much system size I'll be able to reach and in what time? And how should I go about coding it? Thanks, Kartikeya

1 0

Re: How the spin is positioned in the Smatrix?
by Diego Bruno 26 Sep '22

26 Sep '22

Can you disponibilize the code?

4 3

Calculation of Andreev Spectrum for NS junction in KWANT
by Subhadeep Chakraborty 26 Sep '22

26 Sep '22

Hi Kwant Users ! Can anyone tell me if one could obtain the andreev spectrum for an NS junction using KWANT. ___________________________________ It would be a great a help if someone shares a basic idea regarding how to do it ? Thanks in advance,

1 0

Jump to page:

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

Kwant-discuss