diff --git a/TODO b/TODO
deleted file mode 100644
index d1689a87594d9a608b461dcc6361fb5a876bff9d..0000000000000000000000000000000000000000
--- a/TODO
+++ /dev/null
@@ -1,154 +0,0 @@
-Roughly in order of importance.                                     -*-org-*-
-
-* Document the order of sites/orbitals in finalized builders
-
-* Add calculation of current density
-
-* Re-design the interface of low level systems
-  considering the following
-  - We want support for multiple symmetry directions
-  - Consider making System.hamiltonian a vectorized function.  Value functions
-    should be vectorizable as well (when marked with a decorator?)
-  - Consider using two-way-compressed uni-directional graphs in systems
-    advantages:
-    - The number of edges gets reduced by 1/2.
-    - The Hermicity of the Hamiltonian is enforced in a more obvious way.
-    - It's clear at system level which Hamiltonian elements are defined
-      directly.
-  - Verify that the graph slicing algorithm continues to work.
-
-* Improve experience when working with sites in a finalized builder
-  (It should be possible to map sites to nodes efficiently.  Probably, the best
-  way to do it is to sort the sites in finalized builders.)
-  One could then also RLE compress sites in finalized builders.
-
-* Make Kwant objects pickleable
-
-* Improve lattice.neighbors()
-  to support the following cases:
-  - The lattice has lower dimensionality than the space
-  - Hoppings between sites from two separate lattices
-
-  The function should be rewritten in a more robust manner using Delaunay
-  tessellation from Scipy/QHull and a proper choice of the set of sites on
-  which to perform the tessellation.
-
-* Implement "inhomogeneous arrays" for ldos and wavefunc
-  The point is to have ldos or wavefunc return objects that behave like 2d
-  arrays, where the site index and orbital index are separate.  However, this
-  must work in the general case where the number of orbitals per site varies.
-
-* Go through the documentation and apply consistent linking/naming practice:
-  We could adopt the scheme that is used in Python's stdlib.  See for example
-  http://docs.python.org/2/library/threading.html
-  This would mean:
-  - When referring to functions, we always append "()" to their name.
-  - Names of objects (types, functions, ...) that are a documented part of
-    Kwant are always hyperlinks _when_ the object itself is meant.
-  - Concepts that are meant in a more general sense are not hyperlinked, even if
-    there exists an object of the same name in Kwant.  (Consider the usage of
-    linking for the term "Thread" in the document linked above.)
-
-* Harmonize numerical tolerances throughout Kwant
-  By default, a channel should open at the same energy in any solver.
-  If this is possible, we may want to have one central tolerance setting.
-
-* Plotter: show site upon click
-
-* fix physics.noise (docstring, __all__, etc.)
-
-* Verify that selective plotting works.
-
-* Improve support for imperfect lattices / amorphous systems.
-  - Provide a "list of atoms" site group
-  - Allow attaching leads in real space: attach_lead(realspace=True)
-  - Add example / tutorial
-
-* Write an RGF solver that uses graph/slicer.
-
-* Add support for easily adding magnetic field to a system.
-
-* Allow plotting of infinite systems
-
-* Use sparse linear algebra to calculate bands
-  However, SciPy's sparse eigenvalues don't seem to work well.
-
-* Allow attaching leads with further than nearest lead unit cell hoppings.
-  The most easy way to do this is increasing the period of the lead.
-  Alternatively, generalize modes and InfiniteSystem format.
-
-* In finalized leads, only keep the sites of a single lead unit cell.
-
-* Add support for optimization of lead fundamental domains.
-
-* Incorprorate efficient correlated disorder using scipy.spatial.
-
-* Implement the C solver interface.
-
-* Wrap TB_SIM as a solver.
-
-* Benchmark MUMPS and check whether nested dissection would be useful.
-  If yes, implement it.
-
-* Consider implementing "models" module, incorporate units into it.
-  - It is potentially related to vectorized value functions (since modesl are
-    spatially homogeneous)
-  - It may also be related to units, which we may borrow from somewhere,
-    or implement on our own.
-
-* Improve handling of s-matrices with important additional degrees of freedom.
-  i.e. spin, or electron-hole.  Currently, one needs to manually process lead
-  modes and to perform a basis transformation of the scattering matrix.
-
-* Implement automatic search for discrete symmetries of a tight binding model
-  and of conservation laws.
-
-* Revisit Symmetry, make it more general. I.e. allow for description of other
-  discrete symmetries, possibly also symmetries which change site and hopping
-  values (=gauge symmetries).
-
-* Make an option for plotting hoppings as arcs, in order for overlapping
-  hoppings to be distinguishable. Also potentially add arrows to the hoppings.
-
-* Adopt mincut/maxflow algorithm from networkx or python-graph to find the best
-  representation of a lead unit cell.  Fixing unit cells could be postponed
-  until a lead is attached.
-
-* Speed up kwant.digest
-  by using SipHash and Ziggurat algorithm in C/Cython
-
-* Improve finding of modes with zero velocity or infinite decay length
-  Right now a mode with truly zero velocity may be incorrectly identified as
-  propagating or evanescent.
-
-* Implement a pretty printer for lists/tuples of sites
-
-* Ensure 100% test coverage
-
-* Ideas for kwant 2
-** Consider whether to use *args or args in System.hamiltonian
-
-** lead_info
-Currently, lead_info is either of two things, a sequence of propagating modes or a sequence of self-energies.
-- There is no need to have lead_info as an attribute of BlockResult.
-- In the case of SMatrix and GreensFunction this attribute should be renamed to what it is, i.e. propagating_modes and selfenergies.
-- In cases where a common name is needed (i.e. return value of _make_linear_sys), a better name might be lead_data or solved_leads.
-- BlockResult itself could be renamed to something more specific, perhaps GeneralSMatrix.
-
-** Unify the order of indices
-For example in wave_function and in PropagatingModes.
-
-** Take advantage of when a system is known to be Hermitian
-Either
-- omit the lead with most open channels by default when calling smatrix and greens_function or
-- verify that the scattering matrix is indeed current-conserving (it might be not due to numerical error)
-
-** lattices
-- In add_site_family, the first argument cannot be a general site family, it must be a lattice from lattice.py.  We should rename site_family -> lattice in this context.  In fact, that argument could be even a Polyatomic lattice (that is not a site family at all).
-- The automatically selected other_vectors could be longer but more orthogonal.
-
-** modes
-When on a band edge, the calculation of modes should fail in a clear way, for example by raising kwant.UndefinedResultError (itself derived from ValueError).  We could also provide a function (tentative name: kwant.tolerate) that turns such exceptions into None values (useful for plotting).
-
-** hamiltonian_submatrix and hamiltonian
-Systems have both methods and its confusing which does what.