Commit faba2761 by Pablo Piskunow Committed by Joseph Weston

 # Tutorial 2.8. Calculating spectral density with the Kernel Polynomial Method # ============================================================================ # # Physics background # ------------------ # - Chebyshev polynomials, random trace approximation, spectral densities. # # Kwant features highlighted # -------------------------- # - kpm module,kwant operators. import scipy import _defs # For plotting from matplotlib import pyplot as plt #HIDDEN_BEGIN_sys1 # necessary imports import kwant import numpy as np # define the system def make_syst(r=30, t=-1, a=1): syst = kwant.Builder() lat = kwant.lattice.honeycomb(a, norbs=1) def circle(pos): x, y = pos return x ** 2 + y ** 2 < r ** 2 syst[lat.shape(circle, (0, 0))] = 0. syst[lat.neighbors()] = t syst.eradicate_dangling() return syst #HIDDEN_END_sys1 # Plot several density of states curves on the same axes. def plot_dos(labels_to_data): for label, (x, y) in labels_to_data: plt.plot(x, y, label=label, linewidth=2) plt.legend(loc=2, framealpha=0.5) plt.xlabel("energy [t]") plt.ylabel("DoS [a.u.]") plt.show() plt.clf() def site_size_conversion(densities): return 3 * np.abs(densities) / max(densities) # Plot several local density of states maps in different subplots def plot_ldos(fsyst, axes, titles_to_data, file_name=None): for ax, (title, ldos) in zip(axes, titles_to_data): site_size = site_size_conversion(ldos) # convert LDoS to sizes kwant.plot(fsyst, site_size=site_size, site_color=(0, 0, 1, 0.3), ax=ax) ax.set_title(title) ax.set(adjustable='box-forced', aspect='equal') plt.show() plt.clf() def simple_dos_example(): #HIDDEN_BEGIN_kpm1 fsyst = make_syst().finalized() spectrum = kwant.kpm.SpectralDensity(fsyst) #HIDDEN_END_kpm1 #HIDDEN_BEGIN_kpm2 energies, densities = spectrum() #HIDDEN_END_kpm2 #HIDDEN_BEGIN_kpm3 energy_subset = np.linspace(0, 2) density_subset = spectrum(energy_subset) #HIDDEN_END_kpm3 plot_dos([ ('densities', (energies, densities)), ('density subset', (energy_subset, density_subset)), ]) def dos_averaging_example(fsyst): spectrum = kwant.kpm.SpectralDensity(fsyst) #HIDDEN_BEGIN_av1 print('identity resolution:', spectrum.average()) #HIDDEN_END_av1 #HIDDEN_BEGIN_av2 # 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.average(fermi)) #HIDDEN_END_av2 def increasing_accuracy_example(fsyst): spectrum = kwant.kpm.SpectralDensity(fsyst) original_dos = spectrum() # get unaltered DoS #HIDDEN_BEGIN_acc1 spectrum.increase_energy_resolution(tol=0.03) #HIDDEN_END_acc1 increased_resolution_dos = spectrum() plot_dos([ ('density', original_dos), ('higher energy resolution', increased_resolution_dos), ]) #HIDDEN_BEGIN_acc2 # we supply the *total* number of moments and sampling points spectrum.increase_accuracy(num_moments=200, num_rand_vecs=5) #HIDDEN_END_acc2 increased_moments_dos = spectrum() plot_dos([ ('density', original_dos), ('higher number of moments', increased_moments_dos), ]) def operator_example(fsyst): #HIDDEN_BEGIN_op1 # identity matrix matrix_op = scipy.sparse.eye(len(fsyst.sites)) matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op) #HIDDEN_END_op1 #HIDDEN_BEGIN_op2 # 'sum=True' means we sum over all the sites kwant_op = kwant.operator.Density(fsyst, sum=True) operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op) #HIDDEN_END_op2 plot_dos([ ('identity matrix', matrix_spectrum()), ('kwant.operator.Density', operator_spectrum()), ]) def ldos_example(fsyst): #HIDDEN_BEGIN_op3 # 'sum=False' is the default, but we include it explicitly here for clarity. kwant_op = kwant.operator.Density(fsyst, sum=False) local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op) #HIDDEN_END_op3 #HIDDEN_BEGIN_op4 zero_energy_ldos = local_dos(energy=0) finite_energy_ldos = local_dos(energy=1) _, axes = plt.subplots(1, 2, figsize=(12, 7)) plot_ldos(fsyst, axes,[ ('energy = 0', zero_energy_ldos), ('energy = 1', finite_energy_ldos), ]) #HIDDEN_END_op4 def vector_factory_example(fsyst): spectrum = kwant.kpm.SpectralDensity(fsyst) #HIDDEN_BEGIN_fact1 # construct vectors with n random elements -1 or +1. def binary_vectors(n): return np.rint(np.random.random_sample(n)) * 2 - 1 custom_factory = kwant.kpm.SpectralDensity(fsyst, vector_factory=binary_vectors) #HIDDEN_END_fact1 plot_dos([ ('default vector factory', spectrum()), ('binary vector factory', custom_factory()), ]) def bilinear_map_operator_example(fsyst): #HIDDEN_BEGIN_blm rho = kwant.operator.Density(fsyst, sum=True) # sesquilinear map that does the same thing as rho def rho_alt(bra, ket): return np.vdot(bra, ket) rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho) rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt) #HIDDEN_END_blm plot_dos([ ('kwant.operator.Density', rho_spectrum()), ('bilinear operator', rho_alt_spectrum()), ]) def main(): simple_dos_example() fsyst = make_syst().finalized() dos_averaging_example(fsyst) increasing_accuracy_example(fsyst) operator_example(fsyst) ldos_example(fsyst) vector_factory_example(fsyst) bilinear_map_operator_example(fsyst) # Call the main function if the script gets executed (as opposed to imported). # See . if __name__ == '__main__': main()
 ... ... @@ -54,6 +54,8 @@ where we observe the cosine-like dispersion of the square lattice. Close to k=0 this agrees well with the quadratic dispersion this tight-binding Hamiltonian is approximating. .. _closed-systems: Closed systems .............. ... ...
 ############################################################## Calculating spectral density with the kernel polynomial method ############################################################## We have already seen in the ":ref:closed-systems" tutorial that we can use Kwant simply to build Hamiltonians, which we can then directly diagonalize using routines from Scipy. This already allows us to treat systems with a few thousand sites without too many problems. For larger systems one is often not so interested in the exact eigenenergies and eigenstates, but more in the *density of states*. The kernel polynomial method (KPM), is an algorithm to obtain a polynomial expansion of the density of states. It can also be used to calculate the spectral density of arbitrary operators. Kwant has an implementation of the KPM method that is based on the algorithms presented in Ref. [1]_. Roughly speaking, KPM approximates the density of states (or any other spectral density) by expanding the action of the Hamiltonian (and operator of interest) on a (small) set of *random vectors* as a sum of Chebyshev polynomials up to some order, and then averaging. The accuracy of the method can be tuned by modifying the order of the Chebyshev expansion and the number of random vectors. See notes on accuracy_ below for details. .. seealso:: The complete source code of this example can be found in :download:tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py> .. _accuracy: .. specialnote:: Performance and accuracy The KPM method is especially well suited for large systems, and in the case when one is not interested in individual eigenvalues, but rather in obtaining an approximate spectral density. The accuracy in the energy resolution is dominated by the number of moments. The lowest accuracy is at the center of the spectrum, while slightly higher accuracy is obtained at the edges of the spectrum. If we use the KPM method (with the Jackson kernel, see Ref. [1]_) to describe a delta peak at the center of the spectrum, we will obtain a function similar to a Gaussian of width :math:σ=πa/N, where :math:N is the number of moments, and :math:a is the width of the spectrum. On the other hand, the random vectors will *explore* the range of the spectrum, and as the system gets bigger, the number of random vectors that are necessary to sample the whole spectrum reduces. Thus, a small number of random vectors is in general enough, and increasing this number will not result in a visible improvement of the approximation. Introduction ************ Our aim is to use the kernel polynomial method to obtain the spectral density :math:ρ_A(E), as a function of the energy :math:E, of some Hilbert space operator :math:A. We define .. math:: ρ_A(E) = ρ(E) A(E), where :math:A(E) is the expectation value of :math:A for all the eigenstates of the Hamiltonian with energy :math:E, and the density of states is .. math:: ρ(E) = \frac{1}{D} \sum_{k=0}^{D-1} δ(E-E_k), :math:D being the Hilbert space dimension, and :math:E_k the eigenvalues. In the special case when :math:A is the identity, then :math:ρ_A(E) is simply :math:ρ(E), the density of states. Calculating the density of states ********************************* In the following example, we will use the KPM implementation in Kwant to obtain the density of states of a graphene disk. We start by importing kwant and defining our system. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_sys1 :end-before: #HIDDEN_END_sys1 After making a system we can then create a ~kwant.kpm.SpectralDensity object that represents the density of states for this system. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_kpm1 :end-before: #HIDDEN_END_kpm1 The ~kwant.kpm.SpectralDensity can then be called like a function to obtain a sequence of energies in the spectrum of the Hamiltonian, and the corresponding density of states at these energies. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_kpm2 :end-before: #HIDDEN_END_kpm2 When called with no arguments, an optimal set of energies is chosen (these are not evenly distributed over the spectrum, see Ref. [1]_ for details), however it is also possible to provide an explicit sequence of energies at which to evaluate the density of states. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_kpm3 :end-before: #HIDDEN_END_kpm3 .. image:: ../images/kpm_dos.* In addition to being called like functions, ~kwant.kpm.SpectralDensity objects also have a method ~kwant.kpm.SpectralDensity.average which can be used to integrate the density of states against some distribution function over the whole spectrum. If no distribution function is specified, then the uniform distribution is used: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_av1 :end-before: #HIDDEN_END_av1 .. 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 average 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: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_av2 :end-before: #HIDDEN_END_av2 .. literalinclude:: ../images/kpm_total_states.txt .. specialnote:: Stability and performance: spectral bounds The KPM method internally rescales the spectrum of the Hamiltonian to the interval (-1, 1) (see Ref [1]_ for details), which requires calculating the boundaries of the spectrum (using scipy.sparse.linalg.eigsh). This can be very costly for large systems, so it is possible to pass this explicitly as via the bounds parameter when instantiating the ~kwant.kpm.SpectralDensity (see the class documentation for details). Additionally, ~kwant.kpm.SpectralDensity accepts a parameter epsilon, which ensures that the rescaled Hamiltonian (used internally), always has a spectrum strictly contained in the interval (-1, 1). If bounds are not provided then the tolerance on the bounds calculated with scipy.sparse.linalg.eigsh is set to epsilon/2. Increasing the accuracy of the approximation ******************************************** ~kwant.kpm.SpectralDensity has two methods for increasing the accuracy of the method, each of which offers different levels of control over what exactly is changed. The simplest way to obtain a more accurate solution is to use the increase_energy_resolution method: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_acc1 :end-before: #HIDDEN_END_acc1 This will update the number of calculated moments and also the default number of sampling points such that the maximum distance between successive energy points is tol (see notes on accuracy_). .. image:: ../images/kpm_dos_acc.* Alternatively, you can directly increase (or decrease) the number of moments, random vectors and sampling points with the method ~kwant.kpm.SpectralDensity.increase_accuracy. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_acc2 :end-before: #HIDDEN_END_acc2 .. image:: ../images/kpm_dos_r.* In the above example we only increased the number of moments and decreased the number of random vectors. The keyword argument num_sampling_points can be also specified by passing the keyword argument to ~kwant.kpm.SpectralDensity.increase_accuracy. .. _operator_spectral_density: Calculating the spectral density of an operator *********************************************** Above, we saw how to calculate the density of states by creating a ~kwant.kpm.SpectralDensity and passing it a finalized Kwant system. When instantiating a ~kwant.kpm.SpectralDensity we may optionally supply an operator in addition to the system. In this case it is the spectral density of the given operator that is calculated. ~kwant.kpm.SpectralDensity accepts the operators in a few formats: * *explicit matrices* (numpy array of scipy sparse matrices will work) * *operators* from kwant.operator If an explicit matrix is provided then it must have the same shape as the system Hamiltonian. .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_op1 :end-before: #HIDDEN_END_op1 Or, to do the same calculation using kwant.operator.Density: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_op2 :end-before: #HIDDEN_END_op2 Using operators from kwant.operator allows us to calculate quantities such as the *local* density of states by telling the operator not to sum over all the sites of the system: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_op3 :end-before: #HIDDEN_END_op3 ~kwant.kpm.SpectralDensity will properly handle this vector output, which allows us to plot the local density of states at different point in the spectrum: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_op4 :end-before: #HIDDEN_END_op4 .. image:: ../images/kpm_ldos.* This nicely illustrates the edge states of the graphene dot at zero energy, and the bulk states at higher energy. Advanced topics *************** Custom distributions for random vectors ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ By default ~kwant.kpm.SpectralDensity will use random vectors whose components are unit complex numbers with phases drawn from a uniform distribution. There are several reasons why you may wish to make a different choice of distribution for your random vectors, for example to enforce certain symmetries or to only use real-valued vectors. To change how the random vectors are generated, you need only specify a function that takes the dimension of the Hilbert space as a single parameter, and which returns a vector in that Hilbert space: .. literalinclude:: kernel_polynomial_method.py :start-after: #HIDDEN_BEGIN_fact1 :end-before: #HIDDEN_END_fact1 Reproducible calculations ^^^^^^^^^^^^^^^^^^^^^^^^^ Because KPM internally uses random vectors, running the same calculation twice will not give bit-for-bit the same result. However, similarly to the funcions in ~kwant.rmt, the random number generator can be directly manipulated by passing a value to the rng parameter of ~kwant.kpm.SpectralDensity. rng can itself be a random number generator, or it may simply be a seed to pass to the numpy random number generator (that is used internally by default). Defining operators as sesquilinear maps ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Above__, we showed how ~kwant.kpm.SpectralDensity can calculate the spectral density of operators, and how we can define operators by using kwant.operator. If you need even more flexibility, ~kwant.kpm.SpectralDensity will also accept a *function* as its operator