Coverage for /builds/ase/ase/ase/transport/stm.py: 10.81%
111 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3# flake8: noqa
4import time
6import numpy as np
7from scipy.integrate import trapezoid
9from ase.parallel import world
10from ase.transport.greenfunction import GreenFunction
11from ase.transport.selfenergy import LeadSelfEnergy
12from ase.transport.tools import dagger
15class STM:
16 def __init__(self, h1, s1, h2, s2, h10, s10, h20, s20,
17 eta1, eta2, w=0.5, pdos=[], logfile=None):
18 """XXX
20 1. Tip
21 2. Surface
23 h1: ndarray
24 Hamiltonian and overlap matrix for the isolated tip
25 calculation. Note, h1 should contain (at least) one
26 principal layer.
28 h2: ndarray
29 Same as h1 but for the surface.
31 h10: ndarray
32 periodic part of the tip. must include two and only
33 two principal layers.
35 h20: ndarray
36 same as h10, but for the surface
38 The s* are the corresponding overlap matrices. eta1, and eta
39 2 are (finite) infinitesimals. """
41 self.pl1 = len(h10) // 2 # principal layer size for the tip
42 self.pl2 = len(h20) // 2 # principal layer size for the surface
43 self.h1 = h1
44 self.s1 = s1
45 self.h2 = h2
46 self.s2 = s2
47 self.h10 = h10
48 self.s10 = s10
49 self.h20 = h20
50 self.s20 = s20
51 self.eta1 = eta1
52 self.eta2 = eta2
53 self.w = w # asymmetry of the applied bias (0.5=>symmetric)
54 self.pdos = []
55 self.log = logfile
57 def initialize(self, energies, bias=0):
58 """
59 energies: list of energies
60 for which the transmission function should be evaluated.
61 bias.
62 Will precalculate the surface greenfunctions of the tip and
63 surface.
64 """
65 self.bias = bias
66 self.energies = energies
67 nenergies = len(energies)
68 pl1, pl2 = self.pl1, self.pl2
69 nbf1, nbf2 = len(self.h1), len(self.h2)
71 # periodic part of the tip
72 hs1_dii = self.h10[:pl1, :pl1], self.s10[:pl1, :pl1]
73 hs1_dij = self.h10[:pl1, pl1:2 * pl1], self.s10[:pl1, pl1:2 * pl1]
74 # coupling between per. and non. per part of the tip
75 h1_im = np.zeros((pl1, nbf1), complex)
76 s1_im = np.zeros((pl1, nbf1), complex)
77 h1_im[:pl1, :pl1], s1_im[:pl1, :pl1] = hs1_dij
78 hs1_dim = [h1_im, s1_im]
80 # periodic part the surface
81 hs2_dii = self.h20[:pl2, :pl2], self.s20[:pl2, :pl2]
82 hs2_dij = self.h20[pl2:2 * pl2, :pl2], self.s20[pl2:2 * pl2, :pl2]
83 # coupling between per. and non. per part of the surface
84 h2_im = np.zeros((pl2, nbf2), complex)
85 s2_im = np.zeros((pl2, nbf2), complex)
86 h2_im[-pl2:, -pl2:], s2_im[-pl2:, -pl2:] = hs2_dij
87 hs2_dim = [h2_im, s2_im]
89 # tip and surface greenfunction
90 self.selfenergy1 = LeadSelfEnergy(hs1_dii, hs1_dij, hs1_dim, self.eta1)
91 self.selfenergy2 = LeadSelfEnergy(hs2_dii, hs2_dij, hs2_dim, self.eta2)
92 self.greenfunction1 = GreenFunction(self.h1 - self.bias * self.w * self.s1, self.s1,
93 [self.selfenergy1], self.eta1)
94 self.greenfunction2 = GreenFunction(self.h2 - self.bias * (self.w - 1) * self.s2, self.s2,
95 [self.selfenergy2], self.eta2)
97 # Shift the bands due to the bias.
98 bias_shift1 = -bias * self.w
99 bias_shift2 = -bias * (self.w - 1)
100 self.selfenergy1.set_bias(bias_shift1)
101 self.selfenergy2.set_bias(bias_shift2)
103 # tip and surface greenfunction matrices.
104 nbf1_small = nbf1 # XXX Change this for efficiency in the future
105 nbf2_small = nbf2 # XXX -||-
106 coupling_list1 = list(range(nbf1_small)) # XXX -||-
107 coupling_list2 = list(range(nbf2_small)) # XXX -||-
108 self.gft1_emm = np.zeros((nenergies, nbf1_small, nbf1_small), complex)
109 self.gft2_emm = np.zeros((nenergies, nbf2_small, nbf2_small), complex)
111 for e, energy in enumerate(self.energies):
112 if self.log is not None: # and world.rank == 0:
113 T = time.localtime()
114 self.log.write(' %d:%02d:%02d, ' % (T[3], T[4], T[5]) +
115 '%d, %d, %02f\n' % (world.rank, e, energy))
116 gft1_mm = self.greenfunction1.retarded(energy)[coupling_list1]
117 gft1_mm = np.take(gft1_mm, coupling_list1, axis=1)
119 gft2_mm = self.greenfunction2.retarded(energy)[coupling_list2]
120 gft2_mm = np.take(gft2_mm, coupling_list2, axis=1)
122 self.gft1_emm[e] = gft1_mm
123 self.gft2_emm[e] = gft2_mm
125 if self.log is not None and world.rank == 0:
126 self.log.flush()
128 def get_transmission(self, v_12, v_11_2=None, v_22_1=None):
129 """XXX
131 v_12:
132 coupling between tip and surface
133 v_11_2:
134 correction to "on-site" tip elements due to the
135 surface (eq.16). Is only included to first order.
136 v_22_1:
137 corretion to "on-site" surface elements due to he
138 tip (eq.17). Is only included to first order.
139 """
141 dim0 = v_12.shape[0]
142 dim1 = v_12.shape[1]
144 nenergies = len(self.energies)
145 T_e = np.empty(nenergies, float)
146 v_21 = dagger(v_12)
147 for e, energy in enumerate(self.energies):
148 gft1 = self.gft1_emm[e]
149 if v_11_2 is not None:
150 gf1 = np.dot(v_11_2, np.dot(gft1, v_11_2))
151 gf1 += gft1 # eq. 16
152 else:
153 gf1 = gft1
155 gft2 = self.gft2_emm[e]
156 if v_22_1 is not None:
157 gf2 = np.dot(v_22_1, np.dot(gft2, v_22_1))
158 gf2 += gft2 # eq. 17
159 else:
160 gf2 = gft2
162 a1 = (gf1 - dagger(gf1))
163 a2 = (gf2 - dagger(gf2))
164 self.v_12 = v_12
165 self.a2 = a2
166 self.v_21 = v_21
167 self.a1 = a1
168 v12_a2 = np.dot(v_12, a2[:dim1])
169 v21_a1 = np.dot(v_21, a1[-dim0:])
170 self.v12_a2 = v12_a2
171 self.v21_a1 = v21_a1
172 T = -np.trace(np.dot(v12_a2[:, :dim1], v21_a1[:, -dim0:])) # eq. 11
173 assert abs(T.imag).max() < 1e-14
174 T_e[e] = T.real
175 self.T_e = T_e
176 return T_e
178 def get_current(self, bias, v_12, v_11_2=None, v_22_1=None):
179 """Very simple function to calculate the current.
181 Asummes zero temperature.
183 bias: type? XXX
184 bias voltage (V)
186 v_12: XXX
187 coupling between tip and surface.
189 v_11_2:
190 correction to onsite elements of the tip
191 due to the potential of the surface.
192 v_22_1:
193 correction to onsite elements of the surface
194 due to the potential of the tip.
195 """
196 energies = self.energies
197 T_e = self.get_transmission(v_12, v_11_2, v_22_1)
198 bias_window = sorted(-np.array([bias * self.w, bias * (self.w - 1)]))
199 self.bias_window = bias_window
200 # print 'bias window', np.around(bias_window,3)
201 # print 'Shift of tip lead do to the bias:', self.selfenergy1.bias
202 # print 'Shift of surface lead do to the bias:', self.selfenergy2.bias
203 i1 = sum(energies < bias_window[0])
204 i2 = sum(energies < bias_window[1])
205 step = 1
206 if i2 < i1:
207 step = -1
209 return np.sign(bias) * trapezoid(x=energies[i1:i2:step],
210 y=T_e[i1:i2:step])