r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import numpy as np

4class GreenFunction:

5 """Equilibrium retarded Green function."""

7 def __init__(self, H, S=None, selfenergies=[], eta=1e-4):

8 self.H = H

9 self.S = S

10 self.selfenergies = selfenergies

11 self.eta = eta

12 self.energy = None

13 self.Ginv = np.empty(H.shape, complex)

15 def retarded(self, energy, inverse=False):

16 """Get retarded Green function at specified energy.

18 If 'inverse' is True, the inverse Green function is returned (faster).

19 """

20 if energy != self.energy:

21 self.energy = energy

22 z = energy + self.eta * 1.j

24 if self.S is None:

25 self.Ginv[:] = 0.0

26 self.Ginv.flat[:: len(self.S) + 1] = z

27 else:

28 self.Ginv[:] = z

29 self.Ginv *= self.S

30 self.Ginv -= self.H

32 for selfenergy in self.selfenergies:

33 self.Ginv -= selfenergy.retarded(energy)

35 if inverse:

36 return self.Ginv

37 else:

38 return np.linalg.inv(self.Ginv)

40 def calculate(self, energy, sigma):

41 """XXX is this really needed"""

42 ginv = energy * self.S - self.H - sigma

43 return np.linalg.inv(ginv)

45 def apply_retarded(self, energy, X):

46 """Apply retarded Green function to X.

48 Returns the matrix product G^r(e) . X

49 """

50 return np.linalg.solve(self.retarded(energy, inverse=True), X)

52 def dos(self, energy):

53 """Total density of states -1/pi Im(Tr(GS))"""

54 if self.S is None:

55 return -self.retarded(energy).imag.trace() / np.pi

56 else:

57 GS = self.apply_retarded(energy, self.S)

58 return -GS.imag.trace() / np.pi

60 def pdos(self, energy):

61 """Projected density of states -1/pi Im(SGS/S)"""

62 if self.S is None:

63 return -self.retarded(energy).imag.diagonal() / np.pi

64 else:

65 S = self.S

66 SGS = np.dot(S, self.apply_retarded(energy, S))

67 return -(SGS.diagonal() / S.diagonal()).imag / np.pi