1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
import numpy as np
class GreenFunction:
"""Equilibrium retarded Green function."""
def __init__(self, H, S=None, selfenergies=[], eta=1e-4):
self.H = H
self.S = S
self.selfenergies = selfenergies
self.eta = eta
self.energy = None
self.Ginv = np.empty(H.shape, complex)
def retarded(self, energy, inverse=False):
"""Get retarded Green function at specified energy.
If 'inverse' is True, the inverse Green function is returned (faster).
"""
if energy != self.energy:
self.energy = energy
z = energy + self.eta * 1.j
if self.S is None:
self.Ginv[:] = 0.0
self.Ginv.flat[:: len(self.S) + 1] = z
else:
self.Ginv[:] = z
self.Ginv *= self.S
self.Ginv -= self.H
for selfenergy in self.selfenergies:
self.Ginv -= selfenergy.retarded(energy)
if inverse:
return self.Ginv
else:
return np.linalg.inv(self.Ginv)
def calculate(self, energy, sigma):
"""XXX is this really needed"""
ginv = energy * self.S - self.H - sigma
return np.linalg.inv(ginv)
def apply_retarded(self, energy, X):
"""Apply retarded Green function to X.
Returns the matrix product G^r(e) . X
"""
return np.linalg.solve(self.retarded(energy, inverse=True), X)
def dos(self, energy):
"""Total density of states -1/pi Im(Tr(GS))"""
if self.S is None:
return -self.retarded(energy).imag.trace() / np.pi
else:
GS = self.apply_retarded(energy, self.S)
return -GS.imag.trace() / np.pi
def pdos(self, energy):
"""Projected density of states -1/pi Im(SGS/S)"""
if self.S is None:
return -self.retarded(energy).imag.diagonal() / np.pi
else:
S = self.S
SGS = np.dot(S, self.apply_retarded(energy, S))
return -(SGS.diagonal() / S.diagonal()).imag / np.pi
|
