Commit af53b154 authored by Kloss's avatar Kloss

cosmethic changes

parent 0f0c2d28
Pipeline #20333 passed with stages
in 9 minutes and 6 seconds
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
import sys
import cmath
import functools
import warnings
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
import cmath
from mpi4py import MPI
import functools
import warnings
import kwant
import kwant_spectrum
import tkwant
from tkwant import manybody
import kwant_spectrum
warnings.simplefilter('ignore') # cache possibly performance warning
......@@ -73,7 +70,8 @@ def main():
if time <= 200: # refine the intervals only at the beginning for performance
solver.refine_intervals()
error_after = solver.estimate_error()
print_master("time= {}, error before/after refinement= {:10.4e}, {:10.4e}".
print_master("time= {}, error before/after refinement= "
"{:10.4e}, {:10.4e}".
format(time, error_before, error_after))
else:
print_master("time= {}, error estimate= {:10.4e}".
......@@ -111,9 +109,9 @@ def main():
# helper routines for wave function / fixed quadrature solver
spectra = kwant_spectrum.spectra(syst.leads)
boundaries = tkwant.leads.automatic_boundary(spectra, tmax)
ModInterval = functools.partial(manybody.Interval, order=15,
quadrature='gausslegendre')
intervals = manybody.calc_intervals(spectra, occupation, Interval=ModInterval)
Interval = functools.partial(manybody.Interval, order=15,
quadrature='gausslegendre')
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)
intervals = manybody.split_intervals(intervals, 10)
# manybody wave function (non-adaptive)
......
#! /usr/bin/env python3
# -*- coding: utf-8 -*-
import warnings
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
import warnings
import kwant
import tkwant
from tkwant import manybody
warnings.simplefilter('ignore') # cache possibly performance warning
def v(time, t0=20, tau=4):
"""Time dependent perturbation V(t)"""
return (1 + erf((time - t0) / tau)) / 2
def create_system(L):
def create_system(length):
def onsite_potential(site, time):
"""Time dependent onsite potential (static part + V(t))"""
return 1 + v(time)
# system building
......@@ -28,9 +25,9 @@ def create_system(L):
syst = kwant.Builder()
# central scattering region
syst[(lat(x, 0) for x in range(L))] = 1
syst[(lat(x, 0) for x in range(length))] = 1
syst[lat.neighbors()] = -1
# time dependent onsite-potential V(t) at the leftmost site
# time dependent onsite-potential at the leftmost site
syst[lat(0, 0)] = onsite_potential
# add leads
......@@ -71,7 +68,7 @@ def main():
density_operator = kwant.operator.Density(syst)
# do the actual tkwant simulation
solver = manybody.State(syst, tmax=tmax)
solver = tkwant.manybody.State(syst, tmax=tmax)
charge_density = []
for time in times:
......@@ -81,7 +78,8 @@ def main():
# plot the result
charge_density = np.array(charge_density)
charge_density -= charge_density[0] # substract offset at initial time
[plt.plot(times, charge_density[:, i], label='site ' + str(i)) for i in range(length)]
for site in range(length):
plt.plot(times, charge_density[:, site], label='site {}'.format(site))
plt.xlabel(r'time $t$')
plt.ylabel(r'charge density $n$')
plt.legend()
......
import time as timer
import functools as ft
from math import sin, pi
import warnings
import numpy as np
from mpi4py import MPI
import time as timer
import matplotlib.pyplot as plt
from mpi4py import MPI
with warnings.catch_warnings():
warnings.simplefilter('ignore')
......@@ -16,8 +14,8 @@ def am_master():
return MPI.COMM_WORLD.rank == 0
def plot_currents(times, r):
plt.plot(times, r, lw=3, label='tkwant')
def plot_currents(times, current):
plt.plot(times, current, lw=3, label='tkwant')
plt.legend(loc=4)
plt.xlabel(r'time $t$')
plt.ylabel(r'current $I$')
......@@ -33,19 +31,19 @@ def faraday_flux(time, V_b, V_bias, tau):
omega = pi / tau
if time <= 0:
return 0
elif 0 < time < tau:
if 0 < time < tau:
return V_bias * (time - sin(omega * time) / omega) / 2
else:
return V_bias * (time - tau) + V_bias * tau / 2
return V_bias * (time - tau) + V_bias * tau / 2
def make_system(L):
def make_system(length):
lat = kwant.lattice.chain(norbs=1)
syst = kwant.Builder()
syst[map(lat, range(L + 2))] = 2
syst[map(lat, range(length + 2))] = 2
syst[lat.neighbors()] = -1
syst[lat(0)] = syst[lat(L)] = barrier
syst[lat(0)] = syst[lat(length)] = barrier
lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
lead[lat(0)] = 2
......@@ -57,11 +55,19 @@ def make_system(L):
def main():
L = 400
Ef = 1.
length = 400
chemical_potential = 1.
tmax = 500
lat, syst = make_system(L)
V_b = 0.5
V_bias = 0.2
tau = 10
params = {'V_b': V_b, 'V_bias': V_bias, 'tau': tau}
lat, syst = make_system(length)
# add the "dynamic" part of the voltage
tkwant.leads.add_voltage(syst, 0, faraday_flux)
fsyst = syst.finalized()
......@@ -70,34 +76,30 @@ def main():
kwant.plot(syst)
# equal occupation for all lead
occupation = tkwant.manybody.make_occupation(Ef)
V_b = 0.5
V_bias = 0.2
tau = 10
occupation = tkwant.manybody.make_occupation(chemical_potential)
# Create the solver
solver = tkwant.manybody.State(fsyst, tmax, occupation,
params={'V_b':V_b, 'V_bias':V_bias, 'tau':tau})
solver = tkwant.manybody.State(fsyst, tmax, occupation, params=params)
# observables
left_lead_interface = [(lat(L + 1), lat(L))]
current_operator = kwant.operator.Current(fsyst, where=left_lead_interface, sum=True)
left_lead_interface = [(lat(length + 1), lat(length))]
current_operator = kwant.operator.Current(fsyst, where=left_lead_interface,
sum=True)
# loop over time, calculating the current as we go
r = []
current = []
start = timer.process_time()
times = range(tmax)
for time in times:
solver.evolve(time)
result = solver.evaluate(current_operator)
r.append(result)
current.append(result)
# if am_master(): # remove comments to have more verbose output
# print('[{:.2f}s] t={}, current is {}'
# .format(timer.process_time() - start, time, r[-1]))
if am_master():
plot_currents(times, r)
plot_currents(times, current)
if __name__ == '__main__':
......
......@@ -7,9 +7,10 @@ import kwant
import tkwant
def make_system(a=1, t=1.0, r=10):
# Start with an empty tight-binding system and a single square lattice.
# `a` is the lattice constant (by default set to 1 for simplicity).
def make_system(a=1, t=1.0, radius=10):
"""Make a tight binding system on a single square lattice"""
# `a` is the lattice constant and `t` the hopping integral
# both set by default to 1 for simplicity.
lat = kwant.lattice.square(a, norbs=1)
......@@ -18,14 +19,13 @@ def make_system(a=1, t=1.0, r=10):
# Define the quantum dot
def circle(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return rsq < r ** 2
return x ** 2 + y ** 2 < radius ** 2
syst[lat.shape(circle, (0, 0))] = 4 * t
syst[lat.neighbors()] = -t
lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lead[(lat(0, j) for j in range(-r//2 + 2, r//2 - 1))] = 4 * t
lead[(lat(0, j) for j in range(-radius//2 + 2, radius//2 - 1))] = 4 * t
lead[lat.neighbors()] = -t
syst.attach_lead(lead)
......@@ -39,7 +39,7 @@ def faraday_flux(time):
return 0.1 * (time - 10 * sin(0.1 * time)) / 2
def evolve(times, psi, operator, title):
def evolve(times, psi, operator):
ops = []
for time in times:
psi.evolve(time)
......@@ -58,13 +58,14 @@ def main():
# multiplying the hoppings to these sites by exp(-1j * faraday_flux(time))
extra_sites = tkwant.leads.add_voltage(syst, 0, faraday_flux)
lead_syst_hoppings = [(s, site) for site in extra_sites
for s in syst.neighbors(site)
if s not in extra_sites]
for s in syst.neighbors(site)
if s not in extra_sites]
syst = syst.finalized()
# create an observable for calculating the current flowing from the left lead
J = kwant.operator.Current(syst, where=lead_syst_hoppings, sum=True)
current_operator = kwant.operator.Current(syst, where=lead_syst_hoppings,
sum=True)
energy = 1.
tmax = 200 * pi
......@@ -73,7 +74,6 @@ def main():
scattering_states = kwant.wave_function(syst, energy=1., params={'time': 0})
psi_st = scattering_states(0)[0]
# boundary conditions typically used by `tkwant.solve`
simple_boundaries = [tkwant.leads.SimpleBoundary(tmax=tmax)
for l in syst.leads]
......@@ -81,32 +81,30 @@ def main():
# boundary conditions with an absorbing potential that increases
# according to x**n where `x` is the distance into the lead
absorbing_boundaries = [tkwant.leads.MonomialAbsorbingBoundary(
num_cells=100, strength=10, degree=6)
num_cells=100, strength=10, degree=6)
for l in syst.leads]
# create time-dependent wavefunctions that starts in a scattering state
# originating from the left lead, using two different types of boundary
# conditions
solver = tkwant.onebody.solvers.default
psi_simple = tkwant.onebody.WaveFunction(syst=syst, boundaries=simple_boundaries,
psi_init=psi_st, energy=energy,
solver_type=solver)
psi_alt = tkwant.onebody.WaveFunction(syst=syst, boundaries=absorbing_boundaries,
psi_init=psi_st, energy=energy,
solver_type=solver)
psi_init=psi_st, energy=energy)
psi_alt = tkwant.onebody.WaveFunction(syst=syst, boundaries=absorbing_boundaries,
psi_init=psi_st, energy=energy)
# evolve forward in time, calculating the current
times = np.arange(0, tmax)
# Simple boundary conditions
start = timer.process_time()
current_simple = evolve(times, psi_simple, J, 'Simple boundary conditions')
current_simple = evolve(times, psi_simple, current_operator)
stop = timer.process_time()
print('Simple Boundary conditions elapsed time: {:.2f}s'.format(stop - start))
# Absorbing boundary conditions
start = timer.process_time()
current_alt = evolve(times, psi_alt, J, 'Absorbing boundary conditions')
current_alt = evolve(times, psi_alt, current_operator)
stop = timer.process_time()
print('Absorbing Boundary conditions elapsed time: {:.2f}s'.format(stop - start))
......
......@@ -6,33 +6,32 @@ import matplotlib.pyplot as plt
import kwant
import tkwant
def make_system(a=1, t=1.0, r=10, r_time_dep=3):
# Start with an empty tight-binding system and a single square lattice.
# `a` is the lattice constant (by default set to 1 for simplicity).
lat = kwant.lattice.square(a, norbs=1)
def make_system(a=1, t=1.0, radius=10, radius_time_dep=3):
"""Make a tight binding system on a single square lattice"""
# `a` is the lattice constant and `t` the hopping integral
# both set by default to 1 for simplicity.
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
# Define the quantum dot
def circle(r):
def _(pos):
def inner(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return rsq < r ** 2
return _
return x ** 2 + y ** 2 < r ** 2
return inner
def hopx(site1, site2, V, B):
# The magnetic field is controlled by the parameter B
def hopx(site1, site2, magnetic_field_strength):
y = site1.pos[1]
return -t * exp(-1j * B * y)
return -t * exp(-1j * magnetic_field_strength * y)
def potential(site, time, V, B):
def potential(site, time, voltage):
x, y = site.pos
return 4 * t + sqrt(x**2 + y**2) * V * (1 - cos(time))
return 4 * t + sqrt(x**2 + y**2) * voltage * (1 - cos(time))
syst[lat.shape(circle(r), (0, 0))] = 4 * t
syst[lat.shape(circle(r_time_dep), (0, 0))] = potential
syst[lat.shape(circle(radius), (0, 0))] = 4 * t
syst[lat.shape(circle(radius_time_dep), (0, 0))] = potential
# hoppings in x-direction
syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
# hoppings in y-directions
......@@ -41,37 +40,50 @@ def make_system(a=1, t=1.0, r=10, r_time_dep=3):
# It's a closed system, so no leads
return syst
syst = make_system().finalized()
kwant.plot(syst, show=False)
plt.show()
def main():
# construct a tight binding system and plot it
syst = make_system().finalized()
kwant.plot(syst, show=False)
plt.show()
# set parameters
params = {'voltage': 1.0, 'magnetic_field_strength': 1.0}
tparams = params.copy()
tparams['time'] = 0 # the initial time
# create an observable for calculating the average radius of a wavefunction
def radius(site):
x, y = site.pos
return sqrt(x**2 + y**2)
V = 1.0
B = 1.0
density_operator = kwant.operator.Density(syst, radius, sum=True)
H = syst.hamiltonian_submatrix(params={'time':0, 'V':V, 'B':B})
evals, evecs = np.linalg.eigh(H)
# create a time-dependent wavefunction that starts in the ground state
hamiltonian = syst.hamiltonian_submatrix(params=tparams)
eigenvalues, eigenvectors = np.linalg.eigh(hamiltonian)
# create an observable for calculating the average radius of a wavefunction
def radius(site):
x, y = site.pos
return sqrt(x**2 + y**2)
ground_state = tkwant.onebody.WaveFunction(syst=syst, params=params,
energy=eigenvalues[0],
psi_init=eigenvectors[:, 0])
R = kwant.operator.Density(syst, radius, sum=True)
# evolve forward in time, calculating the average radius of the wavefunction
times = np.arange(0, 10*pi, 0.1)
average_radius = []
for time in times:
ground_state.evolve(time)
density = ground_state.evaluate(density_operator)
average_radius.append(density)
# create a time-dependent wavefunction that starts in the ground state
ground_state = tkwant.onebody.WaveFunction(syst=syst, params={'V':V, 'B':B},
energy=evals[0], psi_init=evecs[:, 0])
# plot the result
plt.plot(times, average_radius)
plt.xlabel(r'time $t$')
plt.ylabel(r'average radius $r$')
plt.show()
# evolve forward in time, calculating the average radius of the wavefunction
times = np.arange(0, 10*pi, 0.1)
rad = []
for time in times:
ground_state.evolve(time)
result = ground_state.evaluate(R)
rad.append(result)
plt.plot(times, rad)
plt.xlabel(r'time $t$')
plt.ylabel(r'average radius $r$')
plt.show()
if __name__ == '__main__':
main()
......@@ -157,8 +157,8 @@ occupation = [None] * len(syst.leads)
occupation[0] = tkwant.manybody.make_occupation(chemical_potential, bands=0)
spectra = kwant_spectrum.spectra(syst.leads)
ModInterval = ft.partial(manybody.Interval, order=15)
intervals = manybody.calc_intervals(spectra, occupation, Interval=ModInterval)
Interval = ft.partial(manybody.Interval, order=15)
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)
intervals = manybody.split_intervals(intervals, 28)
boundaries = tkwant.leads.automatic_boundary(spectra, tmax)
tasks = manybody.calc_tasks(intervals, spectra, occupation)
......
......@@ -332,7 +332,7 @@ The quadrature order, as an example, can be changed via:
.. jupyter-execute::
Interval = functools.partial(manybody.Interval, order=10)
intervals = manybody.calc_intervals(spectra, occupation, Interval=Interval)
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)
Finally, the intervals are passed to the state with the keyword ``intervals``:
......@@ -369,7 +369,7 @@ values of the interval. The actual values change depending on the
specific problem. The other parameters have default values: ``order`` is
the quadrature order applied to the interval (generally the number of
point at which the interval is discretized) and ``quadrature`` is the
concrete quadrature rule that will be applied. ``integral_type`` is an
concrete quadrature rule that will be applied. ``integration_variable`` is an
additional information to switch between energy and momentum
integration.
......@@ -387,7 +387,7 @@ Intervals can then be calculated with
.. jupyter-execute::
intervals = manybody.calc_intervals(spectra, occupation, Interval=Interval)
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)
Alternatively, the interval order used by the manybody state is changed via
......@@ -472,24 +472,24 @@ function that selects only the positive velocities. We refer to Ref.
`[1] <#references>`__ for details.
One can switch between the two ways to perform the manybody average with the keyword ``integral_type``.
By default, ``integral_type`` = *momentum*, and the integration is performed over the momentum,
One can switch between the two ways to perform the manybody average with the keyword ``integration_variable``.
By default, ``integration_variable`` = *momentum*, and the integration is performed over the momentum,
which corresponds the second equation.
With ``integral_type`` = *energy* one switches to the first equation and integrates explicitly over
With ``integration_variable`` = *energy* one switches to the first equation and integrates explicitly over
the energy:
.. jupyter-execute::
Interval = functools.partial(manybody.Interval, integral_type='energy')
Interval = functools.partial(manybody.Interval, integration_variable='energy')
Intervals can then be calculated with
.. jupyter-execute::
intervals = manybody.calc_intervals(spectra, occupation, Interval=Interval)
intervals = manybody.calc_intervals(spectra, occupation, interval_type=Interval)
Note that ``manybody.Intervals`` returned from ``manybody.intervals()`` always store
momentum intervals, independent of ``integral_type``.
momentum intervals, independent of ``integration_variable``.
Alternatively, the integral type used by the manybody state is changed via
......
from math import sin, pi
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from math import sin, cos, pi
import kwant
import tkwant
lat = kwant.lattice.square(norbs=1)
def make_circle_system(radius):
def make_system(a=1, t=1.0, radius=10, width=7):
"""Make a tight binding system on a single square lattice"""
# `a` is the lattice constant and `t` the hopping integral
# both set by default to 1 for simplicity.
lat = kwant.lattice.square(a, norbs=1)
syst = kwant.Builder()
syst[lat.shape(lambda pos: np.linalg.norm(pos) < radius, (0, 0))] = 4
syst[lat.neighbors()] = -1
return syst
def make_lead(width):
# Define the quantum dot
def circle(pos):
(x, y) = pos
return x ** 2 + y ** 2 < radius ** 2
syst[lat.shape(circle, (0, 0))] = 4 * t
syst[lat.neighbors()] = -t
lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lead[(lat(0, j) for j in range(-(width-1)//2, width//2 + 1))] = 4
lead[lat.neighbors()] = -1
return lead
lead[(lat(0, j) for j in range(-(width-1)//2, width//2 + 1))] = 4 * t
lead[lat.neighbors()] = -t
syst.attach_lead(lead)
syst.attach_lead(lead.reversed())
return syst
syst = make_circle_system(radius=10)
lead = make_lead(width=7)
syst.attach_lead(lead)
syst.attach_lead(lead.reversed()) ;
def faraday_flux(time):
return 0.1 * (time - 10 * sin(0.1 * time)) / 2
added_sites = tkwant.leads.add_voltage(syst, 0, faraday_flux)
interface_hoppings = [(a, b)
for b in added_sites
for a in syst.neighbors(b) if a not in added_sites]
fsyst = syst.finalized()
def main():
# create a system andadd a time-dependent voltage to lead 0 -- this is
# implemented by adding sites to the system at the interface with the lead
# and multiplying the hoppings to these sites by exp(-1j * faraday_flux(time))
syst = make_system()
added_sites = tkwant.leads.add_voltage(syst, 0, faraday_flux)
interface_hoppings = [(a, b)
for b in added_sites
for a in syst.neighbors(b) if a not in added_sites]
fsyst = syst.finalized()
Q = kwant.operator.Density(fsyst)
J = kwant.operator.Current(fsyst, where=interface_hoppings, sum=True)
current_operator = kwant.operator.Current(fsyst, where=interface_hoppings,
sum=True)
# create a time-dependent wavefunction that starts in a scattering state
# originating from the left lead
scattering_states = kwant.wave_function(fsyst, energy=1., params={'time':0})
psi_st = scattering_states(0)[0]
# create a time-dependent wavefunction that starts in a scattering state
# originating from the left lead
scattering_states = kwant.wave_function(fsyst, energy=1., params={'time': 0})
psi_st = scattering_states(0)[0]
tmax = 200 * pi
times = np.arange(0, tmax)
boundaries = [tkwant.leads.SimpleBoundary(tmax=tmax)] * len(syst.leads)
psi = tkwant.onebody.WaveFunction(psi_st, fsyst, boundaries=boundaries,
energy=1.)
tmax = 200 * pi
times = np.arange(0, tmax)
boundaries = [tkwant.leads.SimpleBoundary(tmax=tmax)] * len(syst.leads)
psi = tkwant.onebody.WaveFunction(psi_st, fsyst, boundaries=boundaries, energy=1.)
# evolve forward in time, calculating the current
current = []
for time in times:
psi.evolve(time)
current.append(psi.evaluate(current_operator))
rho = []
current = []
for time in times:
psi.evolve(time)
rho.append(psi.evaluate(Q))
current.append(psi.evaluate(J))
# plot results
plt.figure(figsize=(18, 7))
plt.figure(figsize=(18, 7))
gs = matplotlib.gridspec.GridSpec(2, 2)
gs = matplotlib.gridspec.GridSpec(2, 2)
ax1 = plt.subplot(gs[0, 0])
ax1.plot(times, current)
ax1.set_ylabel('Current')
ax1 = plt.subplot(gs[0, 0])
ax1.plot(times, current)
ax1.set_ylabel('Current')
ax2 = plt.subplot(gs[1, 0])
ax2.plot(times, 0.5 * np.cos(0.1 * times))
ax2.set_xlabel('Time')
ax2.set_ylabel('Voltage')
ax2 = plt.subplot(gs[1, 0])
ax2.plot(times, 0.5 * np.cos(0.1 * times))
ax2.set_xlabel('Time')
ax2.set_ylabel('Voltage')
ax3 = plt.subplot(gs[:, 1])
kwant.plot(syst,
ax=ax3,
site_size=lambda s: 0.4 if s in added_sites else 0.25,
hop_lw=lambda a, b: 0.3 if (a, b) in interface_hoppings else 0.1,
site_color=lambda s: 'r' if s in added_sites else 'b',
hop_color=lambda a, b: 'r' if (a, b) in interface_hoppings else 'k',
lead_color='grey', num_lead_cells=3)
plt.show()
ax3 = plt.subplot(gs[:, 1])
kwant.plot(syst,
ax=ax3,
site_size=lambda s: 0.4 if s in added_sites else 0.25,
hop_lw=lambda a, b: 0.3 if (a, b) in interface_hoppings else 0.1,
site_color=lambda s: 'r' if s in added_sites else 'b',
hop_color=lambda a, b: 'r' if (a, b) in interface_hoppings else 'k', lead_color='grey', num_lead_cells=3)
plt.show()
if __name__ == '__main__':
main()
import time
from math import sin, pi
import numpy as np
from matplotlib import pyplot as plt
......@@ -6,9 +5,10 @@ import kwant
import tkwant
def make_system(a=1, t=1.0, r=10):
# Start with an empty tight-binding system and a single square lattice.
# `a` is the lattice constant (by default set to 1 for simplicity).
def make_system(a=1, t=1.0, radius=10):
"""Make a tight binding system on a single square lattice"""
# `a` is the lattice constant and `t` the hopping integral
# both set by default to 1 for simplicity.
lat = kwant.lattice.square(a, norbs=1)
......@@ -17,14 +17,13 @@ def make_system(a=1, t=1.0, r=10):
# Define the quantum dot
def circle(pos):
(x, y) = pos
rsq = x ** 2 + y ** 2
return rsq < r ** 2
return x ** 2 + y ** 2 < radius ** 2