Commit 62bc9210 authored by Kloss's avatar Kloss
Browse files

add default W(t) and checks

parent b84bb0aa
......@@ -1457,7 +1457,8 @@ class State:
except TypeError:
len_int = 0
if len_int == 0:
logger.warning('interval list has length zero')
logger.warning('no occupied states found, the chemical potential'
'is probably wrong.')
for interval in intervals:
logger.debug(interval)
if boundaries is None:
......
......@@ -158,7 +158,30 @@ class WaveFunction:
Schrödinger equation.
solver_type : `tkwant.onebody.solvers.default`, optional
The solver used to evolve the wavefunction forward in time.
"""
"""
# The size of the central scattering region. The central scattering
# region can be smaller as the total size (kernel.size) of the system
# in case of boundary conditions.
syst_size = psi_init.size
hamiltonian_size = H0.shape[0]
if syst_size > hamiltonian_size:
raise ValueError('initial condition size={} is larger than the '
'Hamiltonian matrix H0 size={}'
.format(syst_size, hamiltonian_size))
# The perturbation W(t) must, if present, always match the size of the
# central scattering region. The true leads are never time dependent.
if W is not None:
if not W.size == syst_size:
raise ValueError('initial condition size={} must be equal '
'to the perturbation W size={}'
.format(syst_size, W.size))
else:
def W(*args, **kwargs):
pass
W.size = syst_size
# add initial time to the params dict
tparams = add_time_to_params(params, time_name=time_name,
......@@ -173,17 +196,12 @@ class WaveFunction:
# we are starting from an eigenstate, so we need to
# be solving: (H0 - E) @ psibar + W(t) @ (psibar + psi_st)
# and psi = (psibar + psi_st) * exp(-1j * energy * time)
kernel = kernel_type(H0 - energy * sp.eye(H0.shape[0]),
kernel = kernel_type(H0 - energy * sp.eye(hamiltonian_size),
W, tparams, np.asarray(psi_init))
# get the object that will actually do the time stepping
self.solver = solver_type(kernel)
# The size of the central scattering region. The central scattering
# region can be smaller as the total size (kernel.size) of the system
# in case of boundary conditions.
syst_size = psi_init.size
# transform initial psi to psibar and psi_st
# note that the dgl is always solved in the variable psibar
psibar = np.zeros((kernel.size,), complex)
......
......@@ -12,7 +12,9 @@ from math import cos, sin
import itertools as it
import functools as ft
import numpy as np
from numpy.testing import assert_array_almost_equal
import tinyarray as ta
import scipy
import scipy.integrate._ode as ode
import scipy.sparse.linalg as spla
import kwant
......@@ -22,6 +24,7 @@ from scipy.special import erf
import collections
import copy
from ... import leads, system
from .. import solvers, kernels, WaveFunction, ScatteringStates
......@@ -1167,3 +1170,87 @@ def test_time_name_and_start_default_change_scattering_state():
ScatteringStates(fsyst_zeit, energy, lead, tmax, params=params,
wavefunction_type=WaveFunction_zeit)
assert str(exc.value) == "time must be a finite real number"
def test_wavefunction_from_hamiltonian():
def create_system(L):
# system building
lat = kwant.lattice.square(a=1, norbs=1)
syst = kwant.Builder()
# central scattering region
syst[(lat(x, 0) for x in range(L))] = 1
syst[lat.neighbors()] = -1
return syst
num_sites = 400
syst = create_system(num_sites).finalized()
# lattice sites and time steps
xi = np.arange(num_sites)
times = np.arange(0, 600, 50)
# initial condition
k = np.pi / 6
psi_init = np.exp(- 0.001 * (xi - 100)**2 + 1j * k * xi)
# hamiltonian matrix
diag = np.ones(len(xi))
offdiag = - np.ones(len(xi) - 1)
H0 = scipy.sparse.diags([diag, offdiag, offdiag],
[0, 1, -1])
H0_kwant = syst.hamiltonian_submatrix(params={'time': 0}, sparse=True)
assert_array_almost_equal(H0.todense(), H0_kwant.todense())
wave_func = WaveFunction(H0, W=None, psi_init=psi_init)
wave_func_kwant = WaveFunction.from_kwant(syst, psi_init)
for time in times:
wave_func.evolve(time)
wave_func_kwant.evolve(time)
psi = wave_func.psi()
psi_kwant = wave_func_kwant.psi()
assert_array_almost_equal(psi, psi_kwant)
def test_wavefunction_from_hamiltonian_raises():
num_sites = 400
# lattice sites and time steps
xi = np.arange(num_sites)
# initial condition
k = np.pi / 6
psi_init = np.exp(- 0.001 * (xi - 100)**2 + 1j * k * xi)
# hamiltonian matrix which is too small
diag = np.ones(len(xi) - 1)
offdiag = - np.ones(len(xi) - 2)
H0 = scipy.sparse.diags([diag, offdiag, offdiag], [0, 1, -1])
with pytest.raises(ValueError) as exc:
WaveFunction(H0, W=None, psi_init=psi_init)
assert 'initial condition size=400 is larger than the Hamiltonian matrix H0 size=399' in str(exc.value)
# hamiltonian matrix
diag = np.ones(len(xi) - 1)
offdiag = - np.ones(len(xi) - 2)
H0 = scipy.sparse.diags([diag, offdiag, offdiag], [0, 1, -1])
# perturbation which is too small
diag = np.ones(len(xi))
offdiag = - np.ones(len(xi) - 1)
H0 = scipy.sparse.diags([diag, offdiag, offdiag], [0, 1, -1])
def W(*args, **kwargs):
pass
W.size = num_sites - 1
with pytest.raises(ValueError) as exc:
WaveFunction(H0, W, psi_init=psi_init)
assert 'initial condition size=400 must be equal to the perturbation W size=399' in str(exc.value)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment