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

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

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']))

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']),

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

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 '

353 'not defined.')

354 E = self.energies.copy()

355 T_e = self.T_e.copy()

357 if not isinstance(bias, (int, float)):

358 bias = bias[np.newaxis]

359 E = E[:, np.newaxis]

360 T_e = T_e[:, np.newaxis]

362 fl = fermidistribution(E - bias / 2., kB * T)

363 fr = fermidistribution(E + bias / 2., kB * T)

365 if spinpol:

366 return .5 * np.trapz((fl - fr) * T_e, x=E, axis=0)

367 else:

368 return np.trapz((fl - fr) * T_e, x=E, axis=0)

370 def get_transmission(self):

371 self.initialize()

372 self.update()

373 return self.T_e

375 def get_dos(self):

376 self.initialize()

377 self.update()

378 return self.dos_e

380 def get_eigenchannels(self, n=None):

381 """Get ``n`` first eigenchannels."""

382 self.initialize()

383 self.update()

384 if n is None:

385 n = self.input_parameters['eigenchannels']

386 return self.eigenchannels_ne[:n]

388 def get_pdos(self):

389 self.initialize()

390 self.update()

391 return self.pdos_ne

393 def subdiagonalize_bfs(self, bfs, apply=False):

394 self.initialize()

395 bfs = np.array(bfs)

396 p = self.input_parameters

397 h_mm = p['h']

398 s_mm = p['s']

399 ht_mm, st_mm, c_mm, e_m = subdiagonalize(h_mm, s_mm, bfs)

400 if apply:

401 self.uptodate = False

402 h_mm[:] = ht_mm.real

403 s_mm[:] = st_mm.real

404 # Rotate coupling between lead and central region

405 for alpha, sigma in enumerate(self.selfenergies):

406 sigma.h_im[:] = np.dot(sigma.h_im, c_mm)

407 sigma.s_im[:] = np.dot(sigma.s_im, c_mm)

409 c_mm = np.take(c_mm, bfs, axis=0)

410 c_mm = np.take(c_mm, bfs, axis=1)

411 return ht_mm, st_mm, e_m.real, c_mm

413 def cutcoupling_bfs(self, bfs, apply=False):

414 self.initialize()

415 bfs = np.array(bfs)

416 p = self.input_parameters

417 h_pp = p['h'].copy()

418 s_pp = p['s'].copy()

419 cutcoupling(h_pp, s_pp, bfs)

420 if apply:

421 self.uptodate = False

422 p['h'][:] = h_pp

423 p['s'][:] = s_pp

424 for alpha, sigma in enumerate(self.selfenergies):

425 for m in bfs:

426 sigma.h_im[:, m] = 0.0

427 sigma.s_im[:, m] = 0.0

428 return h_pp, s_pp

430 def lowdin_rotation(self, apply=False):

431 p = self.input_parameters

432 h_mm = p['h']

433 s_mm = p['s']

434 eig, rot_mm = linalg.eigh(s_mm)

435 eig = np.abs(eig)

436 rot_mm = np.dot(rot_mm / np.sqrt(eig), dagger(rot_mm))

437 if apply:

438 self.uptodate = False

439 h_mm[:] = rotate_matrix(h_mm, rot_mm) # rotate C region

440 s_mm[:] = rotate_matrix(s_mm, rot_mm)

441 for alpha, sigma in enumerate(self.selfenergies):

442 sigma.h_im[:] = np.dot(sigma.h_im, rot_mm) # rotate L-C coupl.

443 sigma.s_im[:] = np.dot(sigma.s_im, rot_mm)

445 return rot_mm

447 def get_left_channels(self, energy, nchan=1):

448 self.initialize()

449 g_s_ii = self.greenfunction.retarded(energy)

450 lambda_l_ii = self.selfenergies[0].get_lambda(energy)

451 lambda_r_ii = self.selfenergies[1].get_lambda(energy)

453 if self.greenfunction.S is not None:

454 s_mm = self.greenfunction.S

455 s_s_i, s_s_ii = linalg.eig(s_mm)

456 s_s_i = np.abs(s_s_i)

457 s_s_sqrt_i = np.sqrt(s_s_i) # sqrt of eigenvalues

458 s_s_sqrt_ii = np.dot(s_s_ii * s_s_sqrt_i, dagger(s_s_ii))

459 s_s_isqrt_ii = np.dot(s_s_ii / s_s_sqrt_i, dagger(s_s_ii))

461 lambdab_r_ii = np.dot(np.dot(s_s_isqrt_ii, lambda_r_ii), s_s_isqrt_ii)

462 a_l_ii = np.dot(np.dot(g_s_ii, lambda_l_ii), dagger(g_s_ii))

463 ab_l_ii = np.dot(np.dot(s_s_sqrt_ii, a_l_ii), s_s_sqrt_ii)

464 lambda_i, u_ii = linalg.eig(ab_l_ii)

465 ut_ii = np.sqrt(lambda_i / (2.0 * np.pi)) * u_ii

466 m_ii = 2 * np.pi * np.dot(np.dot(dagger(ut_ii), lambdab_r_ii), ut_ii)

467 T_i, c_in = linalg.eig(m_ii)

468 T_i = np.abs(T_i)

470 channels = np.argsort(-T_i)[:nchan]

471 c_in = np.take(c_in, channels, axis=1)

472 T_n = np.take(T_i, channels)

473 v_in = np.dot(np.dot(s_s_isqrt_ii, ut_ii), c_in)

475 return T_n, v_in