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
