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