normalize the total density of states to the system size

 The total density of states now integrates to `N`, the size of the system.

doc/source/images/kernel_polynomial_method.py.diff

@@ -46,7 +46,7 @@
  def simple_dos_example():
      fsyst = make_syst().finalized()
  
-@@ -74,19 +85,25 @@
+@@ -74,18 +85,24 @@
      plot_dos([
          ('densities', (energies, densities)),
          ('density subset', (energy_subset, density_subset)),
@@ -66,16 +66,15 @@
  
      # Fermi energy 0.1 and temperature 0.2
      fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
-     n_states = len(fsyst.sites)  # 1 degree of freedom per site
  
--    print('number of filled states:', n_states * spectrum.integrate(fermi))
+-    print('number of filled states:', spectrum.integrate(fermi))
 +    with open('kpm_total_states.txt', 'w') as f:
 +        with redirect_stdout(f):
-+            print('number of filled states:', n_states * spectrum.integrate(fermi))
++            print('number of filled states:', spectrum.integrate(fermi))
  
  
  def increasing_accuracy_example(fsyst):
-@@ -100,7 +117,9 @@
+@@ -99,7 +116,9 @@
      plot_dos([
          ('density', original_dos),
          ('higher energy resolution', increased_resolution_dos),
@@ -86,7 +85,7 @@
  
      spectrum.add_moments(100)
      spectrum.add_vectors(5)
-@@ -110,7 +129,9 @@
+@@ -109,7 +128,9 @@
      plot_dos([
          ('density', original_dos),
          ('higher number of moments', increased_moments_dos),
@@ -97,7 +96,7 @@
  
  
  def operator_example(fsyst):
-@@ -140,7 +161,9 @@
+@@ -139,7 +160,9 @@
      plot_ldos(fsyst, axes,[
          ('energy = 0', zero_energy_ldos),
          ('energy = 1', finite_energy_ldos),

doc/source/tutorial/kernel_polynomial_method.py

@@ -95,9 +95,8 @@ def dos_integrating_example(fsyst):
 #HIDDEN_BEGIN_int2
     # Fermi energy 0.1 and temperature 0.2
     fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
-    n_states = len(fsyst.sites)  # 1 degree of freedom per site
 
-    print('number of filled states:', n_states * spectrum.integrate(fermi))
+    print('number of filled states:', spectrum.integrate(fermi))
 #HIDDEN_END_int2
 
 

doc/source/tutorial/kpm.rst

@@ -66,7 +66,7 @@ states is
 
 .. math::
 
-  ρ(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k),
+  ρ(E) = \sum_{k=0}^{D-1} δ(E-E_k),
 
 :math:`D` being the Hilbert space dimension, and :math:`E_k` the eigenvalues.
 
@@ -124,11 +124,10 @@ distribution is used:
 
 .. literalinclude:: ../images/kpm_normalization.txt
 
-We see that the integral of the density of states is normalized to 1. If
-we wish to calculate, say, the number of states populated in
-equilibrium, then we should integrate with respect to a Fermi-Dirac
-distribution and multiply by the total number of available states in
-the system:
+We see that the integral of the density of states is normalized to the total
+number of available states in the system. If we wish to calculate, say, the
+number of states populated in equilibrium, then we should integrate with
+respect to a Fermi-Dirac distribution:
 
 .. literalinclude:: kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_int2

kwant/kpm.py

@@ -32,9 +32,9 @@ class SpectralDensity:
     .. math::
        ρ_A(e) = ρ(e) A(e),
 
-    where :math:`ρ(e)` is the density of states, and :math:`A(e)` is
-    the expectation value of :math:`A` for all the eigenstates with
-    energy :math:`e`.
+    where :math:`ρ(e) = \sum_{k=0}^{D-1} δ(E-E_k)` is the density of
+    states, and :math:`A(e)` is the expectation value of :math:`A` for
+    all the eigenstates with energy :math:`e`.
 
     Parameters
     ----------
@@ -352,9 +352,7 @@ class SpectralDensity:
         # sum moments of all random vectors
         moments = np.sum(np.asarray(self._moments_list).real, axis=0)
         # divide by the number of random vectors
-        moments = moments / self.num_vectors
-        # divide by the length of the vectors to normalize
-        moments = moments / self.hamiltonian.shape[0]
+        moments /= self.num_vectors
         # divide by scale factor to reflect the integral rescaling
         return moments / self._a
 

kwant/tests/test_kpm.py

@@ -502,8 +502,9 @@ def test_integrate():
     ham = kwant.rmt.gaussian(dim)
     spectrum = make_spectrum(ham, p)
     ones = lambda x: np.ones_like(x)
-    assert np.abs(spectrum.integrate() - simps(spectrum.densities,
-                                             x=spectrum.energies)) < TOL_SP
+    assert np.abs(
+        (spectrum.integrate() - simps(spectrum.densities, x=spectrum.energies))
+        / spectrum.integrate()) < TOL_SP
     assert np.abs(spectrum.integrate() - spectrum.integrate(
         distribution_function=ones)) < TOL