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