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