Dear Daniel,

sorry for the long delay - I started to work on a reply before
traveling, and then it took some time to get back to your question.

First, with regards to your units question: When you write down the
Hamiltonian in kwant (in kwant, we just have numbers), you choose a unit
of energy (e.g. if all matrix elements are in meV, this unit is meV).
Also, each of your lattice points will correspond to a finite volume in
real-space (for example, in a cubic lattice with lattice constant a this
volume would be a^3).

All the functions in kwant will reflect these units. For example, the
local density of states (LDOS) has units per energy per volume. Kwant's
ldos thus returns for every site the number of electrons per chosen unit
of energy per lattice point volume.

This units are fixed and do not depend on system size. Hence, when you
integrate ldos over space to get the total DOS D(E), D(E) will indeed
scale with system size, as you expected. I include a little example
below where this is demonstrated. If you do not observe this for your
model, I believe this means your model is not correct.

You correctly assume that ldos is essentially a sum over all incoming
scattering wave functions (and divided by 2 pi). One can show that this
is equivalent to the Green's function expression. This also indicates
how you can compute the ldos in the infinite lead, namely by summing
over all propagating modes and dividing by 2 pi. Alternatively, one can
just make a scattering system that is equivalent to an infinite system.
I compare both approaches in an example below.

I suppose that for your problem you rather need the ldos integrated over
energy rather than integrated over space right? The former will give you
charge density.

It turns out that a student of ours is also working on combining kwant
with the Poisson equation! Out of curiosity: What is your experience
with FiPy? Can you recommend using it?

Best regards,

Michael

================================================
Program codes
================================================
# Example: scaling of DOS with system size

import kwant
import numpy as np

def dos(energy, L):
    lat = kwant.lattice.chain()
    sys = kwant.Builder()

    for x in xrange(L):
        sys[lat(x)] = 2
    sys[lat.neighbors()] = -1

    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
    lead[lat(0)] = 2
    lead[lat.neighbors()] = -1

    sys.attach_lead(lead)
    sys.attach_lead(lead.reversed())

    sys = sys.finalized()

    return np.sum(kwant.ldos(sys, energy))

import matplotlib.pyplot as plt

L_arr = range(1, 20)
dos_arr = []

for L in xrange(1, 20):
    dos_arr.append(dos(0.3, L))

plt.plot(L_arr, dos_arr)
plt.show()

=======================================================
# Example: ldos of infinite system

import kwant

import numpy as np
import matplotlib.pyplot as plt

W = 30
lat = kwant.lattice.square()

sys = kwant.Builder()
for y in xrange(W):
    sys[lat(0, y)] = 4
sys[lat.neighbors()] = -1

lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
for y in xrange(W):
    lead[lat(0, y)] = 4
lead[lat.neighbors()] = -1

sys.attach_lead(lead)
sys.attach_lead(lead.reversed())

sys = sys.finalized()

# Compute ldos in scattering region
# (since system is completely translationally
#  invariant, this is equivalent to an infinite
#  system)
sys_ldos = kwant.ldos(sys, energy=0.6)
y_arr = []
for i in xrange(sys.graph.num_nodes):
    y_arr.append(sys.site(i).pos[1])
y_arr = np.asarray(y_arr)
indx = np.argsort(y_arr)

# Compute ldos of infinite system directly

def lead_ldos(lead, energy):
    prop_modes = lead.modes(energy=energy)[0]

    ldos = np.zeros((prop_modes.wave_functions.shape[0],), float)
    for i in xrange(prop_modes.wave_functions.shape[1]):
        ldos += abs(prop_modes.wave_functions[:, i])**2

    return ldos/2.0/np.pi

l_ldos = lead_ldos(sys.leads[0], energy=0.6)
yl_arr = []
for i in xrange(sys.leads[0].cell_size):
    yl_arr.append(sys.leads[0].site(i).pos[1])
yl_arr = np.asarray(yl_arr)
lindx = np.argsort(yl_arr)

plt.plot(y_arr[indx], sys_ldos[indx])
plt.plot(yl_arr[lindx], l_ldos[lindx], 'o')
plt.show()

On 13.09.2014 17:14, Daniel Rodgers-Pryor wrote:
> Hi,
>
> Firstly, thanks for creating Kwant - it's so nice to use physics code
> written by people who understand software-engineering as well as physics :)
>
> I've got a few questions about units and density-of-states in Kwant,
> please respond if you know anything about any of them; don't feel the
> need to respond to them all at once.
>
> I'm trying to add self-consistent electrostatics to my Kwant system
> (using FiPy as a finite-element/volume Poisson-solver). Obviously, I
> need to calculate the electron-carrier-density of the system by
> integrating over the Fermi-Dirac occupation, then feed that (via
> real-space basis functions) to my Poisson solver. I'm not quite sure how
> to handle the density of states in the context of Kwant:
>
> I assume that just summing over the LDOS (integrating over all space)
> will give the density of states as a function of energy, D(E). Plotting
> it seems to produce reasonable bands, but I'm not quite sure about the
> units, or how it scales with system size. In the system I'm modelling
> (low-temperature p-donors in silicon), every lattice site adds an
> electron to the system (the temperatures are low enough that the silicon
> is frozen-out as a conductor and can just be treated as a background
> dielectric constant). I should be able to integrate over the
> density-of-states until the total equals the (known) number of
> electrons, but the density of states obtained by summing over the LDOS
> calculated by Kwant does not scale properly with the number of sites in
> the system; larger systems always need higher Fermi energies, which
> isn't physical at all.
>
> What am I missing here? Are the units of the LDOS Kwant calculates
> somehow normalised? How can I get a density-of-states which scales
> appropriately with the total number of electrons/sites in my system?
>
> The lead unit-cell of my system will need to be solved self-consistently
> too; how can I calculate the local density of states (and thus, via
> Fermi-Dirac, the electron-density) of a lead?
>
> Am I correct in assuming that the LDOS produced by Kwant is equivalent
> to summing over the state-density-weighted scattering-wavefunctions from
> the modes in all leads (and thus that integrating it over the
> occupied-energies will produce a sensible total electron-density)?
>
> Finally, and slightly unrelated, do my chosen energy-units need to be
> accounted-for anywhere in Kwant's Schroedinger-solutions? I'm writing my
> Hamiltonian terms in meV; will bands, LDOS etc. all naturally scale to
> make this choice transparent? Similarly, does effective electron-mass
> need to be accounted for at all?
>
> Thanks so much for your help,
>
> Daniel R-P