Coverage for ase / transport / calculators.py: 60.30%
267 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3import numpy as np
4from numpy import linalg
5from scipy.integrate import trapezoid
7from ase.transport.greenfunction import GreenFunction
8from ase.transport.selfenergy import BoxProbe, LeadSelfEnergy
9from ase.transport.tools import (
10 cutcoupling,
11 dagger,
12 fermidistribution,
13 rotate_matrix,
14 subdiagonalize,
15)
16from ase.units import kB
19class TransportCalculator:
20 """Determine transport properties of a device sandwiched between
21 two semi-infinite leads using a Green function method.
22 """
24 def __init__(self, **kwargs):
25 """Create the transport calculator.
27 Parameters
28 ----------
30 h : (N, N) ndarray
31 Hamiltonian matrix for the central region.
32 s : {None, (N, N) ndarray}, optional
33 Overlap matrix for the central region.
34 Use None for an orthonormal basis.
35 h1 : (N1, N1) ndarray
36 Hamiltonian matrix for lead1.
37 h2 : {None, (N2, N2) ndarray}, optional
38 Hamiltonian matrix for lead2. You may use None if lead1 and lead2
39 are identical.
40 s1 : {None, (N1, N1) ndarray}, optional
41 Overlap matrix for lead1. Use None for an orthonormal basis.
42 hc1 : {None, (N1, N) ndarray}, optional
43 Hamiltonian coupling matrix between the first principal
44 layer in lead1 and the central region.
45 hc2 : {None, (N2, N} ndarray), optional
46 Hamiltonian coupling matrix between the first principal
47 layer in lead2 and the central region.
48 sc1 : {None, (N1, N) ndarray}, optional
49 Overlap coupling matrix between the first principal
50 layer in lead1 and the central region.
51 sc2 : {None, (N2, N) ndarray}, optional
52 Overlap coupling matrix between the first principal
53 layer in lead2 and the central region.
54 energies : {None, array_like}, optional
55 Energy points for which calculated transport properties are
56 evaluated.
57 eta : {1.0e-5, float}, optional
58 Infinitesimal for the central region Green function.
59 eta1/eta2 : {1.0e-5, float}, optional
60 Infinitesimal for lead1/lead2 Green function.
61 align_bf : {None, int}, optional
62 Use align_bf=m to shift the central region
63 by a constant potential such that the m'th onsite element
64 in the central region is aligned to the m'th onsite element
65 in lead1 principal layer.
66 logfile : {None, str}, optional
67 Write a logfile to file with name `logfile`.
68 Use '-' to write to std out.
69 eigenchannels: {0, int}, optional
70 Number of eigenchannel transmission coefficients to
71 calculate.
72 pdos : {None, (N,) array_like}, optional
73 Specify which basis functions to calculate the
74 projected density of states for.
75 dos : {False, bool}, optional
76 The total density of states of the central region.
77 box: XXX
78 YYY
80 If hc1/hc2 are None, they are assumed to be identical to
81 the coupling matrix elements between nearest neighbor
82 principal layers in lead1/lead2.
84 Examples
85 --------
87 >>> import numpy as np
88 >>> h = np.array((0,)).reshape((1,1))
89 >>> h1 = np.array((0, -1, -1, 0)).reshape(2,2)
90 >>> energies = np.arange(-3, 3, 0.1)
91 >>> calc = TransportCalculator(h=h, h1=h1, energies=energies)
92 >>> T = calc.get_transmission()
94 """
96 # The default values for all extra keywords
97 self.input_parameters = {'energies': None,
98 'h': None,
99 'h1': None,
100 'h2': None,
101 's': None,
102 's1': None,
103 's2': None,
104 'hc1': None,
105 'hc2': None,
106 'sc1': None,
107 'sc2': None,
108 'box': None,
109 'align_bf': None,
110 'eta1': 1e-5,
111 'eta2': 1e-5,
112 'eta': 1e-5,
113 'logfile': None,
114 'eigenchannels': 0,
115 'dos': False,
116 'pdos': []}
118 self.initialized = False # Changed Hamiltonians?
119 self.uptodate = False # Changed energy grid?
120 self.set(**kwargs)
122 def set(self, **kwargs):
123 for key in kwargs:
124 if key in ['h', 'h1', 'h2', 'hc1', 'hc2',
125 's', 's1', 's2', 'sc1', 'sc2',
126 'eta', 'eta1', 'eta2', 'align_bf', 'box']:
127 self.initialized = False
128 self.uptodate = False
129 break
130 elif key in ['energies', 'eigenchannels', 'dos', 'pdos']:
131 self.uptodate = False
132 elif key not in self.input_parameters:
133 raise KeyError('%r not a vaild keyword' % key)
135 self.input_parameters.update(kwargs)
136 log = self.input_parameters['logfile']
137 if log is None:
138 class Trash:
139 def write(self, s):
140 pass
142 def flush(self):
143 pass
145 self.log = Trash()
146 elif log == '-':
147 from sys import stdout
148 self.log = stdout
149 elif 'logfile' in kwargs:
150 self.log = open(log, 'w')
152 def initialize(self):
153 if self.initialized:
154 return
156 print('# Initializing calculator...', file=self.log)
158 p = self.input_parameters
159 if p['s'] is None:
160 p['s'] = np.identity(len(p['h']))
162 identical_leads = False
163 if p['h2'] is None:
164 p['h2'] = p['h1'] # Lead2 is idendical to lead1
165 identical_leads = True
167 if p['s1'] is None:
168 p['s1'] = np.identity(len(p['h1']))
170 if identical_leads:
171 p['s2'] = p['s1']
172 else:
173 if p['s2'] is None:
174 p['s2'] = np.identity(len(p['h2']))
176 h_mm = p['h']
177 s_mm = p['s']
178 pl1 = len(p['h1']) // 2
179 pl2 = len(p['h2']) // 2
180 h1_ii = p['h1'][:pl1, :pl1]
181 h1_ij = p['h1'][:pl1, pl1:2 * pl1]
182 s1_ii = p['s1'][:pl1, :pl1]
183 s1_ij = p['s1'][:pl1, pl1:2 * pl1]
184 h2_ii = p['h2'][:pl2, :pl2]
185 h2_ij = p['h2'][pl2: 2 * pl2, :pl2]
186 s2_ii = p['s2'][:pl2, :pl2]
187 s2_ij = p['s2'][pl2: 2 * pl2, :pl2]
189 if p['hc1'] is None:
190 nbf = len(h_mm)
191 h1_im = np.zeros((pl1, nbf), complex)
192 s1_im = np.zeros((pl1, nbf), complex)
193 h1_im[:pl1, :pl1] = h1_ij
194 s1_im[:pl1, :pl1] = s1_ij
195 p['hc1'] = h1_im
196 p['sc1'] = s1_im
197 else:
198 h1_im = p['hc1']
199 if p['sc1'] is not None:
200 s1_im = p['sc1']
201 else:
202 s1_im = np.zeros(h1_im.shape, complex)
203 p['sc1'] = s1_im
205 if p['hc2'] is None:
206 h2_im = np.zeros((pl2, nbf), complex)
207 s2_im = np.zeros((pl2, nbf), complex)
208 h2_im[-pl2:, -pl2:] = h2_ij
209 s2_im[-pl2:, -pl2:] = s2_ij
210 p['hc2'] = h2_im
211 p['sc2'] = s2_im
212 else:
213 h2_im = p['hc2']
214 if p['sc2'] is not None:
215 s2_im = p['sc2']
216 else:
217 s2_im = np.zeros(h2_im.shape, complex)
218 p['sc2'] = s2_im
220 align_bf = p['align_bf']
221 if align_bf is not None:
222 diff = ((h_mm[align_bf, align_bf] - h1_ii[align_bf, align_bf]) /
223 s_mm[align_bf, align_bf])
224 print('# Aligning scat. H to left lead H. diff=', diff,
225 file=self.log)
226 h_mm -= diff * s_mm
228 # Setup lead self-energies
229 # All infinitesimals must be > 0
230 assert np.all(np.array((p['eta'], p['eta1'], p['eta2'])) > 0.0)
231 self.selfenergies = [LeadSelfEnergy((h1_ii, s1_ii),
232 (h1_ij, s1_ij),
233 (h1_im, s1_im),
234 p['eta1']),
235 LeadSelfEnergy((h2_ii, s2_ii),
236 (h2_ij, s2_ij),
237 (h2_im, s2_im),
238 p['eta2'])]
239 box = p['box']
240 if box is not None:
241 print('Using box probe!')
242 self.selfenergies.append(
243 BoxProbe(eta=box[0], a=box[1], b=box[2], energies=box[3],
244 S=s_mm, T=0.3))
246 # setup scattering green function
247 self.greenfunction = GreenFunction(selfenergies=self.selfenergies,
248 H=h_mm,
249 S=s_mm,
250 eta=p['eta'])
252 self.initialized = True
254 def update(self):
255 if self.uptodate:
256 return
258 p = self.input_parameters
259 self.energies = p['energies']
260 nepts = len(self.energies)
261 nchan = p['eigenchannels']
262 pdos = p['pdos']
263 self.T_e = np.empty(nepts)
264 if p['dos']:
265 self.dos_e = np.empty(nepts)
266 if pdos != []:
267 self.pdos_ne = np.empty((len(pdos), nepts))
268 if nchan > 0:
269 self.eigenchannels_ne = np.empty((nchan, nepts))
271 for e, energy in enumerate(self.energies):
272 Ginv_mm = self.greenfunction.retarded(energy, inverse=True)
273 lambda1_mm = self.selfenergies[0].get_lambda(energy)
274 lambda2_mm = self.selfenergies[1].get_lambda(energy)
275 a_mm = linalg.solve(Ginv_mm, lambda1_mm)
276 b_mm = linalg.solve(dagger(Ginv_mm), lambda2_mm)
277 T_mm = np.dot(a_mm, b_mm)
278 if nchan > 0:
279 t_n = linalg.eigvals(T_mm).real
280 self.eigenchannels_ne[:, e] = np.sort(t_n)[-nchan:]
281 self.T_e[e] = np.sum(t_n)
282 else:
283 self.T_e[e] = np.trace(T_mm).real
285 print(energy, self.T_e[e], file=self.log)
286 self.log.flush()
288 if p['dos']:
289 self.dos_e[e] = self.greenfunction.dos(energy)
291 if pdos != []:
292 self.pdos_ne[:, e] = np.take(self.greenfunction.pdos(energy),
293 pdos)
295 self.uptodate = True
297 def print_pl_convergence(self):
298 self.initialize()
299 pl1 = len(self.input_parameters['h1']) // 2
301 h_ii = self.selfenergies[0].h_ii
302 s_ii = self.selfenergies[0].s_ii
303 ha_ii = self.greenfunction.H[:pl1, :pl1]
304 sa_ii = self.greenfunction.S[:pl1, :pl1]
305 c1 = np.abs(h_ii - ha_ii).max()
306 c2 = np.abs(s_ii - sa_ii).max()
307 print(f'Conv (h,s)={c1:.2e}, {c2:2.e}')
309 def plot_pl_convergence(self):
310 self.initialize()
311 pl1 = len(self.input_parameters['h1']) // 2
312 hlead = self.selfenergies[0].h_ii.real.diagonal()
313 hprincipal = self.greenfunction.H.real.diagonal[:pl1]
315 import pylab as pl
316 pl.plot(hlead, label='lead')
317 pl.plot(hprincipal, label='principal layer')
318 pl.axis('tight')
319 pl.show()
321 def get_current(self, bias, T=0., E=None, T_e=None, spinpol=False):
322 '''Returns the current as a function of the
323 bias voltage.
325 **Parameters:**
326 bias : {float, (M,) ndarray}, units: V
327 Specifies the bias voltage.
328 T : {float}, units: K, optional
329 Specifies the temperature.
330 E : {(N,) ndarray}, units: eV, optional
331 Contains energy grid of the transmission function.
332 T_e {(N,) ndarray}, units: unitless, optional
333 Contains the transmission function.
334 spinpol: {bool}, optional
335 Specifies whether the current should be
336 calculated assuming degenerate spins
338 **Returns:**
339 I : {float, (M,) ndarray}, units: 2e/h*eV
340 Contains the electric current.
342 Examples
343 --------
345 >> import numpy as np
346 >> import pylab as plt
347 >> from ase import units
348 >>
349 >> bias = np.arange(0, 2, .1)
350 >> current = calc.get_current(bias, T = 0.)
351 >> plt.plot(bias, 2.*units._e**2/units._hplanck*current)
352 >> plt.xlabel('U [V]')
353 >> plt.ylabel('I [A]')
354 >> plt.show()
356 '''
357 if E is not None:
358 if T_e is None:
359 self.energies = E
360 self.uptodate = False
361 T_e = self.get_transmission().copy()
362 else:
363 assert self.uptodate, ('Energy grid and transmission function '
364 'not defined.')
365 E = self.energies.copy()
366 T_e = self.T_e.copy()
368 if not isinstance(bias, (int, float)):
369 bias = bias[np.newaxis]
370 E = E[:, np.newaxis]
371 T_e = T_e[:, np.newaxis]
373 fl = fermidistribution(E - bias / 2., kB * T)
374 fr = fermidistribution(E + bias / 2., kB * T)
376 if spinpol:
377 return .5 * trapezoid((fl - fr) * T_e, x=E, axis=0)
378 else:
379 return trapezoid((fl - fr) * T_e, x=E, axis=0)
381 def get_transmission(self):
382 self.initialize()
383 self.update()
384 return self.T_e
386 def get_dos(self):
387 self.initialize()
388 self.update()
389 return self.dos_e
391 def get_eigenchannels(self, n=None):
392 """Get ``n`` first eigenchannels."""
393 self.initialize()
394 self.update()
395 if n is None:
396 n = self.input_parameters['eigenchannels']
397 return self.eigenchannels_ne[:n]
399 def get_pdos(self):
400 self.initialize()
401 self.update()
402 return self.pdos_ne
404 def subdiagonalize_bfs(self, bfs, apply=False):
405 self.initialize()
406 bfs = np.array(bfs)
407 p = self.input_parameters
408 h_mm = p['h']
409 s_mm = p['s']
410 ht_mm, st_mm, c_mm, e_m = subdiagonalize(h_mm, s_mm, bfs)
411 if apply:
412 self.uptodate = False
413 h_mm[:] = ht_mm.real
414 s_mm[:] = st_mm.real
415 # Rotate coupling between lead and central region
416 for alpha, sigma in enumerate(self.selfenergies):
417 sigma.h_im[:] = np.dot(sigma.h_im, c_mm)
418 sigma.s_im[:] = np.dot(sigma.s_im, c_mm)
420 c_mm = np.take(c_mm, bfs, axis=0)
421 c_mm = np.take(c_mm, bfs, axis=1)
422 return ht_mm, st_mm, e_m.real, c_mm
424 def cutcoupling_bfs(self, bfs, apply=False):
425 self.initialize()
426 bfs = np.array(bfs)
427 p = self.input_parameters
428 h_pp = p['h'].copy()
429 s_pp = p['s'].copy()
430 cutcoupling(h_pp, s_pp, bfs)
431 if apply:
432 self.uptodate = False
433 p['h'][:] = h_pp
434 p['s'][:] = s_pp
435 for alpha, sigma in enumerate(self.selfenergies):
436 for m in bfs:
437 sigma.h_im[:, m] = 0.0
438 sigma.s_im[:, m] = 0.0
439 return h_pp, s_pp
441 def lowdin_rotation(self, apply=False):
442 p = self.input_parameters
443 h_mm = p['h']
444 s_mm = p['s']
445 eig, rot_mm = linalg.eigh(s_mm)
446 eig = np.abs(eig)
447 rot_mm = np.dot(rot_mm / np.sqrt(eig), dagger(rot_mm))
448 if apply:
449 self.uptodate = False
450 h_mm[:] = rotate_matrix(h_mm, rot_mm) # rotate C region
451 s_mm[:] = rotate_matrix(s_mm, rot_mm)
452 for alpha, sigma in enumerate(self.selfenergies):
453 sigma.h_im[:] = np.dot(sigma.h_im, rot_mm) # rotate L-C coupl.
454 sigma.s_im[:] = np.dot(sigma.s_im, rot_mm)
456 return rot_mm
458 def get_left_channels(self, energy, nchan=1):
459 self.initialize()
460 g_s_ii = self.greenfunction.retarded(energy)
461 lambda_l_ii = self.selfenergies[0].get_lambda(energy)
462 lambda_r_ii = self.selfenergies[1].get_lambda(energy)
464 if self.greenfunction.S is not None:
465 s_mm = self.greenfunction.S
466 s_s_i, s_s_ii = linalg.eig(s_mm)
467 s_s_i = np.abs(s_s_i)
468 s_s_sqrt_i = np.sqrt(s_s_i) # sqrt of eigenvalues
469 s_s_sqrt_ii = np.dot(s_s_ii * s_s_sqrt_i, dagger(s_s_ii))
470 s_s_isqrt_ii = np.dot(s_s_ii / s_s_sqrt_i, dagger(s_s_ii))
472 lambdab_r_ii = np.dot(np.dot(s_s_isqrt_ii, lambda_r_ii), s_s_isqrt_ii)
473 a_l_ii = np.dot(np.dot(g_s_ii, lambda_l_ii), dagger(g_s_ii))
474 ab_l_ii = np.dot(np.dot(s_s_sqrt_ii, a_l_ii), s_s_sqrt_ii)
475 lambda_i, u_ii = linalg.eig(ab_l_ii)
476 ut_ii = np.sqrt(lambda_i / (2.0 * np.pi)) * u_ii
477 m_ii = 2 * np.pi * np.dot(np.dot(dagger(ut_ii), lambdab_r_ii), ut_ii)
478 T_i, c_in = linalg.eig(m_ii)
479 T_i = np.abs(T_i)
481 channels = np.argsort(-T_i)[:nchan]
482 c_in = np.take(c_in, channels, axis=1)
483 T_n = np.take(T_i, channels)
484 v_in = np.dot(np.dot(s_s_isqrt_ii, ut_ii), c_in)
486 return T_n, v_in