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@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
Animation - 20 lattice points.mp4
Description: Animation - 20 lattice points.mp4
import kwant.continuum
import numpy as np
import matplotlib.pyplot as plt
import kwant
from numpy import vectorize
from matplotlib import animation
from tkwant import onebody
import tkwant
import os.path
import json
import time as timelib
from matplotlib import animation
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
self.a_0_SI = 5.2917721090380e-11
self.total_length_SI = 0.66e-6
self.omega_0_eV = 1e-3 # in eV
self.B_0_au = 2.127659574468085e-08
self.b_sl_au = 2.946499088192983e-12
self.eV_0_au = 0.0003674930360069677
self.pulse_frequency_au = self.ev_to_hartree(self.omega_0_eV)
self.total_length_au = self.m_to_au(self.total_length_SI)
self.lattice_size_au = self.total_length_au / self.number_of_lattices
self.mu_B_au = .5
self.z_au = np.arange(-self.total_length_au / 2, self.total_length_au /
2, self.lattice_size_au,
dtype=np.double)
def ev_to_hartree(self, ev):
return ev / 2.72114e1
def hartree_to_ev(self, hartree):
return hartree * 2.72114e1
def m_to_au(self, m):
return m / self.a_0_SI
def second_to_au(self, time):
return time * 4.1341373336493e16
def au_to_second(self, time):
return time / 4.1341373336493e16
def cosine_v_ac(self, time, z, eV_0_au, pulse_frequency_au,
total_length_au):
return ((eV_0_au * np.cos(2 * np.pi * pulse_frequency_au * time)) * z)
/ total_length_au
def sine_v_ac(self, time, z, eV_0_au, pulse_frequency_au, total_length_au):
return ((eV_0_au * np.sin(2 * np.pi * pulse_frequency_au * time)) * z)
/ total_length_au
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()
self.B_y_au = 0
def B_constant(z):
return -self.B_0_au
def C_constant(z):
return -self.b_sl_au * z
def D_constant(z):
return -self.B_y_au
# coefficient for the kinetic energy term
A_constant = 1
# 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 of unperturbed hamiltonian
self.psi_1_init = eigenVectors[:, 0]
self.psi_2_init = eigenVectors[:, 1]
# tkwant 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, lattices):
self.evolve_state = True
# extract lowest two energies (our qubit states)
E_1 = np.real(self.initial_eigenvalues[0])
E_2 = np.real(self.initial_eigenvalues[1])
# compute the difference in these energies
self.delta_E = np.abs(E_2 - E_1)
# compute the resonant frequency for the rabi oscillations
self.pulse_frequency_au = self.delta_E / (2 * np.pi)
total_osc_time = self.second_to_au(0.5e-8)
times_au = np.linspace(0, total_osc_time, num=time_steps)
density_operator = kwant.operator.Density(self.syst)
psi = self.spin_up_state
if self.pertubation_type == "cos":
y3_func = self.cosine_v_ac
else:
y3_func = self.sine_v_ac
# empty array for density and perturbation function
density = np.zeros((time_steps, lattices))
perturb = np.zeros((time_steps, lattices))
array = np.arange(0, time_steps, 1)
for i in array:
# evolve the state with time
psi.evolve(times_au[i])
density[i] = psi.evaluate(density_operator)
perturb[i] = self.hartree_to_ev(y3_func(times_au[i], self.z_au,
self.eV_0_au, self.pulse_frequency_au, self.total_length_au))
self.evolve_state = False
return density, perturb
def animations(self, density, y3, nsteps):
# create a figure with two subplots
fig_animation, axs = plt.subplots(2,1)
# 1920 x 1080
w_in_inches = 10
h_in_inches = 6
dpi = 100
fig_animation.set_size_inches(w_in_inches, h_in_inches, True)
fig_animation.set_dpi(dpi)
fig_animation.tight_layout(pad=5.0)
# set variables for each of the axes
ax1 = axs[0]
ax2 = axs[1]
# initialise two line objects (one in each axes)
line1, = ax1.plot([], [], lw=2)
line2, = ax2.plot([], [], lw=2, color='r', label='$H_1(t)$')
line = [line1, line2]
z_max_au = self.m_to_au(self.total_length_SI) / 2
z_max_SI = self.total_length_SI / 2
y2_max = self.hartree_to_ev((self.eV_0_au * z_max_au /
self.total_length_au))
y1 = density[0]
y1_max = np.max(y1) * 2
lattice_size_SI = self.total_length_SI / self.number_of_lattices
z_SI = np.arange(-self.total_length_SI / 2, self.total_length_SI / 2,
lattice_size_SI, dtype=np.double)
# PDF
ax1.set_xlim(-z_max_SI, z_max_SI)
ax1.set_ylim(0, np.max(density[-1,:]*1.2))
ax1.set_ylabel("$|\psi(z,t)|^2$")
ax1.set_xlabel("z ($m$)")
ax1.set_title("PDF over time")
# Potential
ax2.set_xlim(-z_max_SI, z_max_SI)
ax2.set_ylim(-y2_max, y2_max)
ax2.set_ylabel("$E$ (eV)")
ax2.set_xlabel("z ($m$)")
# initialization function: plot the background of each frame
def init():
line[0].set_data([], [])
line[1].set_data([], [])
return line
def animate(i):
# update line objects.
line[0].set_data(z_SI, density[i,:])
line[1].set_data(z_SI, y3[i,:])
ax2.legend(loc="upper right")
return line
# call the animator. blit=True means only re-draw the parts that have
changed.
anim = animation.FuncAnimation(fig_animation, animate, init_func=init,
frames=nsteps-1, interval=300, blit=True)
#folder_path = r"\Users\"
#anim.save('{}/animation.mp4'.format(folder_path), writer='ffmpeg')
plt.show()
return True
def main():
lattices = 20
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();
nsteps = 5
density, perturbation = system.evolve(nsteps,lattices)
system.animations(density, perturbation, nsteps)
if __name__ == '__main__':
main()
