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