Hi Michael, Thanks for getting back to me. Sorry for the delay in getting back to you too. 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).
Neat; that clarifies a lot of things. 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?
FiPy's solvers perform well, and it's easy to set up on windows and linux
(I haven't tried OSX), but the documentation isn't great and it has trouble
with complex boundary conditions (particularly internal
boundary conditions; look elsewhere if you want to fix internal values
easily and reliably). Some other technical and specific things: I had to
download the dev version and tinker with the source a bit to get certain
kinds of periodic 3d boundary conditions working, but those changes should
all be in the dev branch now; sampling from FiPy grids (even uniform grids)
can be very slow unless you pre-compute the nearest cells to your sampling
points (the default behaviour does an O(n) lookup to find the nearest cell
to every point you want to sample from).
I've got some decent code for interfacing between kwant and FiPy in 2D
(with matching square lattices), and some uglier code for arbitrary 3D
lattices; I'll put it up on github and send you the link in a week or so.
Regards,
Daniel R-P
On 25 October 2014 09:27, Michael Wimmer
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