Coverage for kwant/physics/leads.py : 97%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# Copyright 2011-2016 Kwant authors. # # This file is part of Kwant. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution and at # http://kwant-project.org/license. A list of Kwant authors can be found in # the file AUTHORS.rst at the top-level directory of this distribution and at # http://kwant-project.org/authors.
csr_matrix)
# TODO: Once Kwant depends on numpy >= 1.11, remove the fix if split_fixed: # skip coverage split = np.split else: # skip coverage def split(array, n, axis=0): if array.shape[axis] != 0: return np.split(array, n, axis) else: return n * [array]
# TODO: Use scipy block_diag once we depend on scipy>=0.19 # Throws ValueError, but if fixed ensure that works as intended except ValueError: # skip coverage bdiag_broken = True if bdiag_broken: # skip coverage def block_diag(*matrices): """Construct a block diagonal matrix out of the input matrices.
Like scipy.linalg.block_diag, but works for zero size matrices.""" rows, cols = np.sum([mat.shape for mat in matrices], axis=0) b_mat = np.zeros((rows,cols), dtype='complex') rows, cols = 0, 0 for mat in matrices: new_rows = rows + mat.shape[0] new_cols = cols + mat.shape[1] b_mat[rows:new_rows, cols:new_cols] = mat rows, cols = new_rows, new_cols return b_mat
"""Least squares version that works also with 0-shaped matrices."""
"""Check whether a discrete symmetry relation between two blocks of the Hamiltonian vanishes or not.
For a discrete symmetry S, the projection that maps from the block i to block j of the Hamiltonian with projectors P_i and P_j is S_ji = P_j^+ S P_i. This function determines whether S_ji vanishes or not, i. e. whether a symmetry relation exists between blocks i and j.
Here, discrete symmetries and projectors are in canonical form, such that S maps each block at most either to itself or to a single other block. In other words, for each j, the block S_ji is nonzero at most for one i, while all other blocks vanish. If a nonzero block exists, it is unitary and hence has norm 1. To identify nonzero blocks without worrying about numerical errors, we thus check that its norm is larger than 0.5. """
"""Split and rearrange a list of arrays.
Combine a list of arrays into a single array where the first halves of each array appear first: `[a b], [c d], [e f] -> [a c e b d f]` """
# Container classes
"""The calculated propagating modes of a lead.
Attributes ---------- wave_functions : numpy array The wave functions of the propagating modes. momenta : numpy array Momenta of the modes. velocities : numpy array Velocities of the modes. block_nmodes: list of integers Number of left or right moving propagating modes per conservation law block of the Hamiltonian.
Notes ===== The sort order of all the three arrays is identical. The first half of the modes have negative velocity, the second half have positive velocity. Within these halves the modes are ordered by the eigenvalues of any declared conservation laws. Within blocks with the same conservation law eigenvalue the modes with negative velocity are ordered by increasing momentum, and the modes with positive velocity are ordered by decreasing momentum. Finally, modes are ordered by the magnitude of their velocity. To summarize, the modes are ordered according to the key `(sign(v), conserved_quantity, sign(v) * k , abs(v))` where `v` is velocity, `k` is momentum and `conserved_quantity` is the conservation law eigenvalue.
In the above, the positive velocity and momentum directions are defined with respect to the translational symmetry direction of the system.
The first dimension of `wave_functions` corresponds to the orbitals of all the sites in a unit cell, the second one to the number of the mode. Each mode is normalized to carry unit current. If several modes have the same momentum and velocity, an arbitrary orthonormal basis in the subspace of these modes is chosen.
If a conservation law is specified to block diagonalize the Hamiltonian to N blocks, then `block_nmodes[i]` is the number of left or right moving propagating modes in conservation law block `i`. The ordering of blocks is the same as the ordering of the projectors used to block diagonalize the Hamiltonian. """
"""Stabilized eigendecomposition of the translation operator.
Due to the lack of Hermiticity of the translation operator, its eigendecomposition is frequently poorly conditioned. Solvers in Kwant use this stabilized decomposition of the propagating and evanescent modes in the leads. If the hopping between the unit cells of an infinite system is invertible, the translation eigenproblem is written in the basis `psi_n, h_hop^+ * psi_(n+1)`, with ``h_hop`` the hopping between unit cells. If `h_hop` is not invertible, and has the singular value decomposition `u s v^+`, then the eigenproblem is written in the basis `sqrt(s) v^+ psi_n, sqrt(s) u^+ psi_(n+1)`. In this basis we calculate the eigenvectors of the propagating modes, and the Schur vectors (an orthogonal basis) in the space of evanescent modes.
`vecs` and `vecslmbdainv` are the first and the second halves of the wave functions. The first `nmodes` are eigenmodes moving in the negative direction (hence they are incoming into the system in Kwant convention), the second `nmodes` are eigenmodes moving in the positive direction. The remaining modes are the Schur vectors of the modes evanescent in the positive direction. Propagating modes with the same eigenvalue are orthogonalized, and all the propagating modes are normalized to carry unit current. Finally the `sqrt_hop` attribute is `v sqrt(s)`.
Attributes ---------- vecs : numpy array Translation eigenvectors. vecslmbdainv : numpy array Translation eigenvectors divided by the corresponding eigenvalues. nmodes : int Number of left-moving (or right-moving) modes. sqrt_hop : numpy array or None Part of the SVD of `h_hop`, or None if the latter is invertible. """
""" Compute the self-energy generated by lead modes.
Returns ------- Sigma : numpy array, real or complex, shape (M,M) The computed self-energy. Note that even if `h_cell` and `h_hop` are both real, `Sigma` will typically be complex. (More precisely, if there is a propagating mode, `Sigma` will definitely be complex.) """
# Auxiliary functions that perform different parts of the calculation. """Make an eigenvalue problem for eigenvectors of translation operator.
Parameters ---------- h_cell : numpy array with shape (n, n) Hamiltonian of a single lead unit cell. h_hop : numpy array with shape (n, m), m <= n Hopping Hamiltonian from a cell to the next one. tol : float Numbers are considered zero when they are smaller than `tol` times the machine precision. stabilization : sequence of 2 booleans or None Which steps of the eigenvalue problem stabilization to perform. If the value is `None`, then Kwant chooses the fastest (and least stable) algorithm that is expected to be sufficient. For any other value, Kwant forms the eigenvalue problem in the basis of the hopping singular values. The first element set to `True` forces Kwant to add an anti-Hermitian term to the cell Hamiltonian before inverting. If it is set to `False`, the extra term will only be added if the cell Hamiltonian isn't invertible. The second element set to `True` forces Kwant to solve a generalized eigenvalue problem, and not to reduce it to the regular one. If it is `False`, reduction to a regular problem is performed if possible.
Returns ------- linsys : namedtuple A named tuple containing `matrices` a matrix pencil defining the eigenproblem, `v` a hermitian conjugate of the last matrix in the hopping singular value decomposition, and functions for extracting the wave function in the unit cell from the wave function in the basis of the nonzero singular exponents of the hopping.
Notes ----- The lead problem with degenerate hopping is rather complicated, and the details of the algorithm will be published elsewhere.
"""
if not np.any(h_hop): # skip coverage # Inter-cell hopping is zero. The current algorithm is not suited to # treat this extremely singular case. raise ValueError("Inter-cell hopping is exactly zero.")
# If both h and t are real, it may be possible to use the real eigenproblem.
# First check if the hopping matrix has singular values close to 0. # (Close to zero is defined here as |x| < eps * tol * s[0] , where # s[0] is the largest singular value.)
# The hopping matrix is well-conditioned and can be safely inverted. # Hence the regular transfer matrix may be used.
else:
# The hopping matrix has eigenvalues close to 0 - those # need to be eliminated.
# Recast the svd of h_hop = u s v^dagger such that # u, v are matrices with shape n x n_nonsing. # pad v with zeros if necessary
# Eliminating the zero eigenvalues requires inverting the on-site # Hamiltonian, possibly including a self-energy-like term. The # self-energy-like term stabilizes the inversion, but the most stable # choice is inherently complex. This can be disadvantageous if the # Hamiltonian is real, since staying in real arithmetics can be # significantly faster. The strategy here is to add a complex # self-energy-like term always if the original Hamiltonian is complex, # and check for invertibility first if it is real
not matrices_real) # Check if there is a chance we will not need to add an imaginary term.
else:
# Matrices are complex or need self-energy-like term to be # stabilized.
# If the condition number of the stabilized h is # still bad, there is nothing we can do. raise RuntimeError("Flat band encountered at the requested " "energy, result is badly defined.")
# The function that can extract the full wave function psi from # the projected one (v^dagger psi lambda^-1, s u^dagger psi).
dot(u, psi[n_nonsing:] * lmbdainv))
# Setup the generalized eigenvalue problem.
# Solving a generalized eigenproblem is about twice as expensive # as solving a regular eigenvalue problem. # Computing the LU factorization is negligible compared to both # (approximately 1/30th of a regular eigenvalue problem). # Because of this, it makes sense to try to reduce # the generalized eigenvalue problem to a regular one, provided # the matrix B can be safely inverted.
# A more stringent condition is used here since errors can # accumulate from here to the eigenvalue calculation later.
else:
"""A helper routine for modes(), that wraps eigenproblems.
This routine wraps the regular and general eigenproblems that can arise in a unfied way.
Parameters ---------- a : numpy array The matrix on the left hand side of a regular or generalized eigenvalue problem. b : numpy array or None The matrix on the right hand side of the generalized eigenvalue problem. tol : float The tolerance for separating eigenvalues with absolute value 1 from the rest.
Returns ------- ev : numpy array An array of eigenvalues (can contain NaNs and Infs, but those are not accessed in `modes()`) The number of eigenvalues equals twice the number of nonzero singular values of `h_hop` (so `2*h_cell.shape[0]` if `h_hop` is invertible). evanselect : numpy integer array Index array of right-decaying modes. propselect : numpy integer array Index array of propagating modes (both left and right). vec_gen(select) : function A function that computes the eigenvectors chosen by the array select. ord_schur(select) : function A function that computes the unitary matrix (corresponding to the right eigenvector space) of the (general) Schur decomposition reordered such that the eigenvalues chosen by the array select are in the top left block. """
# Right-decaying modes. # Propagating modes.
else:
# Right-decaying modes. # Propagating modes.
# Note: the division is OK here, since we later only access # eigenvalues close to the unit circle
calc_ev=False)[2]
"""Makes the wave functions that have the same velocity at a time-reversal invariant momentum (TRIM) particle-hole symmetric.
If P is the particle-hole operator and P^2 = 1, then a particle-hole symmetric wave function at a TRIM is an eigenstate of P with eigenvalue 1. If P^2 = -1, wave functions with the same velocity at a TRIM come in pairs. Such a pair is particle-hole symmetric if the wave functions are related by P, i. e. the pair can be expressed as [psi_n, P psi_n] where psi_n is a wave function.
To ensure proper ordering of modes, this function also returns an array of indices which ensures that particle-hole partners are properly ordered in this subspace of modes. These are later used with np.lexsort to ensure proper ordering.
Parameters ---------- wfs : numpy array A matrix of propagating wave functions at a TRIM that all have the same velocity. The orthonormal wave functions form the columns of this matrix. particle_hole : numpy array The matrix representation of the unitary part of the particle-hole operator, expressed in the tight binding basis.
Returns ------- new_wfs : numpy array The matrix of particle-hole symmetric wave functions. TRIM_sort: numpy integer array Index array that stores the proper sort order of particle-hole symmetric wave functions in this subspace. """
"""Apply the particle-hole operator to an array. """
# Take P in the subspace of W = wfs: U = W^+ @ P @ W^*. # Check that wfs are orthonormal and the space spanned # by them is closed under ph, meaning U is unitary. raise ValueError('wfs are not orthonormal or not closed under particle_hole.') else: raise ValueError('particle_hole must square to +1 or -1. P_squared = {}'.format(P_squared))
# Use the matrix square root method from # Applied Mathematics and Computation 234 (2014) 380-384. # Schur decomposition guarantees that vecs are orthonormal. # U should be normal, so vals is diagonal. # Need to take safe square root of U, the branch cut should not go # through any eigenvalues. Otherwise the square root may not be symmetric. # Find largest gap between eigenvalues # Take matrix square root with branch cut in largest gap # For symmetric U sqrt(U) is also symmetric. # We want a new basis W_new such that W_new^+ @ P @ W_new^* = 1. # This is achieved by W_new = W @ sqrt(U). # If P^2 = 1, there is no need to sort the modes further. else: # P^2 = -1. # Iterate over wave functions to construct # particle-hole partners. # The number of modes. This is always an even number >=2. # If there are only two modes in this subspace, they are orthogonal # so we replace the second one with the P applied to the first one. # Store psi_n and P psi_n. # If there are more than two modes, iterate over wave functions # and construct their particle-hole partners one by one. else: # We construct pairs of modes that are particle-hole partners. # Need to iterate over all pairs except the final one. # Take a mode psi_n from the basis - the first column # of the matrix of remaining modes. # Store psi_n and P psi_n. # Remove psi_n and P psi_n from the basis matrix of modes. # First remove psi_n. # Now we project the remaining modes onto the orthogonal # complement of P psi_n. projector: np.outer(P_wf, P_wf.T.conj()) # After the projection, the mode matrix is rank deficient - # the span of the column space has dimension one less than # the number of columns. # Remove the redundant column. # If this is the final iteration, we only have two modes # left and can construct particle-hole partners without # the projection. # Store psi_n and P psi_n. np.eye(wfs.shape[1])) col in new_wfs]) # Store sort ordering in this subspace of modes
time_reversal, chiral): """ Find, normalize and sort the propagating eigenmodes.
Special care is taken of the case of degenerate k-values, where the numerically computed modes are typically a superposition of the real modes. In this case, also the proper (orthogonal) modes are computed. """
# Array for the velocities.
# Array of indices to sort modes at a TRIM by PHS.
# Calculate the full wave function in real space.
# Cast the types if any of the symmetry operators is complex
# Find clusters of nearby eigenvalues. Since the eigenvalues occupy the # unit circle, special care has to be taken to not introduce a cut at # lambda = -1. > eps).flatten() + 1
# Detect the singular case of all eigenvalues equal.
# If there is a degenerate eigenvalue with several different # eigenvectors, the numerical routines return some arbitrary # overlap of the real, physical solutions. In order # to figure out the correct wave function, we need to # have the full, not the projected wave functions # (at least to our current knowledge).
# Finding the true modes is done in two steps:
# 1. The true transversal modes should be orthogonal to each other, as # they share the same Bloch momentum (note that transversal modes with # different Bloch momenta k1 and k2 need not be orthogonal, the full # modes are orthogonal because of the longitudinal dependence e^{i k1 # x} and e^{i k2 x}). The modes with the same k are therefore # orthogonalized. Moreover for the velocity to have a proper value the # modes should also be normalized.
# If the eigenvectors were purely real up to this stage, # they will typically become complex after the rotation. psi = psi.astype(np.common_type(psi, r)) full_psi = full_psi.astype(np.common_type(psi, q))
# 2. Moving infinitesimally away from the degeneracy # point, the modes should diagonalize the velocity # operator (i.e. when they are non-degenerate any more) # The modes are therefore rotated properly such that they # diagonalize the velocity operator. # Note that step 2. does not give a unique result if there are # two modes with the same velocity, or if the modes stay # degenerate even for a range of Bloch momenta (and hence # must have the same velocity). However, this does not matter, # since we are happy with any superposition in this case.
# If the eigenvectors were purely real up to this stage, # they will typically become complex after the rotation.
# With particle-hole symmetry, treat TRIMs individually. # Particle-hole conserves velocity. # If P^2 = 1, we can pick modes at a TRIM as particle-hole eigenstates. # If P^2 = -1, a mode at a TRIM and its particle-hole partner are # orthogonal, and we pick modes such that they are related by # particle-hole symmetry.
# At a TRIM, propagating translation eigenvalues are +1 or -1. (np.abs(np.abs(lmbdainv[indx].real) - 1) < eps).all()): # Set the eigenvalues to the exact TRIM values. else: # Momenta are the negative arguments of the translation # eigenvalues, as computed below using np.angle. np.angle of -1 # is pi, so this assigns k = -pi to modes with translation # eigenvalue -1.
# Original wave functions
# Modes are already sorted by velocity in ascending order, as # returned by la.eigh above. The first half is thus incident, # and the second half outgoing. # Here we work within a subspace of modes with a fixed velocity. # Mostly, this is done to ensure that modes of different velocities # are not mixed when particle-hole partners are constructed for # P^2 = -1. First, we identify which modes have the same velocity. # In each such subspace of modes, we construct wave functions that # are particle-hole partners. # Velocities are sorted in ascending order. Find the indices of the # last instance of each unique velocity. if np.abs(vel-vels[ind+1])>vel_eps] # Now possess an iterator over tuples, where each tuple (i,j) # contains the starting and final indices i and j of a submatrix # of the modes matrix, such that all modes in the submatrix # have the same velocity.
# Iterate over all submatrices of modes with the same velocity. # Pick out wave functions that have a given velocity # Make particle-hole symmetric modes # Store sorting indices of the TRIM modes with the given # velocity. # Gather into a matrix of modes # Store the sort order of all modes at the TRIM. # Used later with np.lexsort when the ordering # of modes is done. # Replace the old modes. # For both cases P^2 = +-1, must rotate wave functions in the # singular value basis. Find the rotation from new basis to old. # Rotate the wave functions in the singular value basis
# Ensure proper usage of chiral symmetry. # No least squares below because the modes should be orthogonal.
raise RuntimeError("Found a mode with zero or close to zero velocity.") raise RuntimeError("Numbers of left- and right-propagating " "modes differ, possibly due to a numerical " "instability.")
# Sort the modes. The modes are sorted first by velocity and momentum, # and finally TRIM modes are properly ordered. -np.sign(velocities) * momenta, np.sign(velocities)])
# TODO: Remove the check once we depend on numpy>=1.8.
# Use particle-hole symmetry to relate modes that propagate in the # same direction but at opposite momenta. # Modes are sorted by velocity (first incident, then outgoing). # Incident modes are then sorted by momentum in ascending order, # and outgoing modes in descending order. # Adopted convention is to get modes with negative k (both in and out) # by applying particle-hole operator to modes with positive k. # With particle-hole symmetry, N must be an even number. # Incident modes # Original wave functions with negative values of k # For incident modes, ordering of wfs by momentum as returned by kwant # is [-k2, -k1, k1, k2], if k2, k1 > 0 and k2 > k1. # To maintain this ordering with ki and -ki as particle-hole partners, # reverse the order of the product at the end. (full_psi[:, :N][:, positive_k]).conj())[:, ::-1] psi[:, :N][:, positive_k[::-1]].dot(rot)
# Outgoing modes # Original wave functions with negative values of k # For outgoing modes, ordering of wfs by momentum as returned by kwant # is like [k2, k1, -k1, -k2], if k2, k1 > 0 and k2 > k1.
# Reverse order at the end to match momenta of opposite sign. full_psi[:, N:][:, positive_k].conj())[:, ::-1] psi[:, N:][:, positive_k[::-1]].dot(rot)
# Modes are ordered by velocity. # Use time-reversal symmetry to relate modes of opposite velocity. # Note: within this function, nmodes refers to the total number # of propagating modes, not either left or right movers.
time_reversal, particle_hole, chiral): """Calculate modes corresponding to a single projector. """
# Defer most of the calculation to helper routines. *(matrices + (tol,)))
# v is never None. # h_hop.shape[0] and v.shape[1] not always the same. # Adding this makes the difference for some tests
# Compute the propagating eigenvectors. # Compute their velocity, and, if necessary, rotate them.
# prop_vecs here is 'psi' in make_proper_modes, i.e. the wf in the SVD # basis. It is in turn used to construct vecs and vecslmbdainv (the # propagating parts). The evanescent parts of vecs and vecslmbdainv come # from evan_vecs above. extract, tol, particle_hole, time_reversal, chiral)
# Prepare output for a single block
particle_hole=None, chiral=None): """Transform the modes data for a given block of the Hamiltonian using a discrete symmetry (see arguments). The symmetry operator can also be specified as can also be identity, for the case when blocks are identical.
Assume that modes_data has the form returned by compute_block_modes, i.e. a tuple (wave_functions, momenta, velocities, nmodes, vecs, vecslmbdainv, v) containing the block modes data.
Assume that the symmetry operator is written in the proper basis (block basis, not full TB).
"""
# Copy to not overwrite modes from previous blocks
nmodes * (np.arange(2*nmodes) // nmodes)) nmodes * (np.arange(2*nmodes) < nmodes)) else: # skip coverage raise ValueError("No relation between blocks was provided.")
discrete_symmetry=None, projectors=None, time_reversal=None, particle_hole=None, chiral=None): """Compute the eigendecomposition of a translation operator of a lead.
Parameters ---------- h_cell : numpy array, real or complex, shape (N,N) The unit cell Hamiltonian of the lead unit cell. h_hop : numpy array, real or complex, shape (N,M) The hopping matrix from a lead cell to the one on which self-energy has to be calculated (and any other hopping in the same direction). tol : float Numbers and differences are considered zero when they are smaller than `tol` times the machine precision. stabilization : sequence of 2 booleans or None Which steps of the eigenvalue prolem stabilization to perform. If the value is `None`, then Kwant chooses the fastest (and least stable) algorithm that is expected to be sufficient. For any other value, Kwant forms the eigenvalue problem in the basis of the hopping singular values. The first element set to `True` forces Kwant to add an anti-Hermitian term to the cell Hamiltonian before inverting. If it is set to `False`, the extra term will only be added if the cell Hamiltonian isn't invertible. The second element set to `True` forces Kwant to solve a generalized eigenvalue problem, and not to reduce it to the regular one. If it is `False`, reduction to a regular problem is performed if possible. Selecting the stabilization manually is mostly necessary for testing purposes. particle_hole : sparse or dense square matrix The unitary part of the particle-hole symmetry operator. time_reversal : sparse or dense square matrix The unitary part of the time-reversal symmetry operator. chiral : sparse or dense square matrix The chiral symmetry operator. projectors : an iterable of sparse or dense matrices Projectors that block diagonalize the Hamiltonian in accordance with a conservation law.
Returns ------- propagating : `~kwant.physics.PropagatingModes` Contains the array of the wave functions of propagating modes, their momenta, and their velocities. It can be used to identify the gauge in which the scattering problem is solved. stabilized : `~kwant.physics.StabilizedModes` A basis of propagating and evanescent modes used by the solvers.
Notes ----- The sorting of the propagating modes is fully described in the documentation for `~kwant.physics.PropagatingModes`. In simple cases where bands do not cross, this ordering corresponds to "lowest modes first". In general, however, it is necessary to examine the band structure -- something this function is not doing by design.
Propagating modes with the same momentum are orthogonalized. All the propagating modes are normalized by current.
This function uses the most stable and efficient algorithm for calculating the mode decomposition that the Kwant authors are aware about. Its details are to be published. """
raise ValueError("Incompatible matrix sizes for h_cell and h_hop.")
# Avoid the trouble of dealing with non-square hopping matrices. # TODO: How to avoid this while not doing a lot of book-keeping?
# Provide default values to not deal with special cases.
(time_reversal, particle_hole, chiral))
# We need the extra transposes to ensure that sparse dot is used. # b.T.conj() @ a @ b.conj() else: # b.T.conj() @ a @ b
# Conservation law basis
# Check that the Hamiltonian has the conservation law # Symmetries that project from block x to block y symm in symmetries] # We did not compute this block yet. *symmetries) else: # Modes in the block already computed. continue
np.allclose(hop_cons[x], hop_cons[y])): else: *symmetries) vecs, vecslmbdainv, sqrt_hops) = zip(*block_modes)
# Reorder by direction of propagation zip(wave_functions, projectors)]).T # Propagating modes object to return group_halves(momenta)) # Store the number of modes per block as an attribute. # nmodes[i] is the number of left or right moving modes in block i. # In the module that makes leads with conservation laws, this is necessary # to keep track of the block structure of the scattering matrix.
for n, v in zip(nmodes, vecs)))
for n, v in zip(nmodes, vecslmbdainv)))
zip(projectors, sqrt_hops)])[:h_hop.shape[1]]
""" Compute the self-energy generated by the lead.
Parameters ---------- h_cell : numpy array, real or complex, shape (N,N) The unit cell Hamiltonian of the lead unit cell. h_hop : numpy array, real or complex, shape (N,M) The hopping matrix from a lead cell to the one on which self-energy has to be calculated (and any other hopping in the same direction). tol : float Numbers are considered zero when they are smaller than `tol` times the machine precision.
Returns ------- Sigma : numpy array, real or complex, shape (M,M) The computed self-energy. Note that even if `h_cell` and `h_hop` are both real, `Sigma` will typically be complex. (More precisely, if there is a propagating mode, `Sigma` will definitely be complex.)
Notes ----- For simplicity this function internally calculates the modes first. This may cause a small slowdown, and can be improved if necessary. """
""" Calculate analytically the self energy for a square lattice.
The lattice is assumed to have a single orbital per site and nearest-neighbor hopping.
Parameters ---------- width : integer width of the lattice """
# Following appendix C of M. Wimmer's diploma thesis: # http://www.physik.uni-regensburg.de/forschung/\ # richter/richter/media/research/publications2004/wimmer-Diplomarbeit.pdf
# p labels transversal modes. i and j label the sites of a cell.
# Precalculate the transverse wave function.
# Precalculate the integrals of the longitudinal wave functions. else:
# Put everything together into the self energy and return it. psi_p_i[:, i] * psi_p_i[:, j].conj() * f_p).sum() |