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