Hide keyboard shortcuts

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 

2 

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 

9 

10 

11class TransportCalculator: 

12 """Determine transport properties of a device sandwiched between 

13 two semi-infinite leads using a Green function method. 

14 """ 

15 

16 def __init__(self, **kwargs): 

17 """Create the transport calculator. 

18 

19 Parameters: 

20 

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 

70 

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. 

74 

75 Examples: 

76 

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() 

83 

84 """ 

85 

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': []} 

107 

108 self.initialized = False # Changed Hamiltonians? 

109 self.uptodate = False # Changed energy grid? 

110 self.set(**kwargs) 

111 

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) 

124 

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 

131 

132 def flush(self): 

133 pass 

134 

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

141 

142 def initialize(self): 

143 if self.initialized: 

144 return 

145 

146 print('# Initializing calculator...', file=self.log) 

147 

148 p = self.input_parameters 

149 if p['s'] is None: 

150 p['s'] = np.identity(len(p['h'])) 

151 

152 identical_leads = False 

153 if p['h2'] is None: 

154 p['h2'] = p['h1'] # Lead2 is idendical to lead1 

155 identical_leads = True 

156 

157 if p['s1'] is None: 

158 p['s1'] = np.identity(len(p['h1'])) 

159 

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

165 

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] 

178 

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 

194 

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 

209 

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 

217 

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

235 

236 # setup scattering green function 

237 self.greenfunction = GreenFunction(selfenergies=self.selfenergies, 

238 H=h_mm, 

239 S=s_mm, 

240 eta=p['eta']) 

241 

242 self.initialized = True 

243 

244 def update(self): 

245 if self.uptodate: 

246 return 

247 

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

260 

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 

274 

275 print(energy, self.T_e[e], file=self.log) 

276 self.log.flush() 

277 

278 if p['dos']: 

279 self.dos_e[e] = self.greenfunction.dos(energy) 

280 

281 if pdos != []: 

282 self.pdos_ne[:, e] = np.take(self.greenfunction.pdos(energy), 

283 pdos) 

284 

285 self.uptodate = True 

286 

287 def print_pl_convergence(self): 

288 self.initialize() 

289 pl1 = len(self.input_parameters['h1']) // 2 

290 

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

298 

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] 

304 

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() 

310 

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. 

314 

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 

327 

328 **Returns:** 

329 I : {float, (M,) ndarray}, units: 2e/h*eV 

330 Contains the electric current. 

331 

332 Examples: 

333 

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() 

344 

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() 

356 

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

358 bias = bias[np.newaxis] 

359 E = E[:, np.newaxis] 

360 T_e = T_e[:, np.newaxis] 

361 

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

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

364 

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) 

369 

370 def get_transmission(self): 

371 self.initialize() 

372 self.update() 

373 return self.T_e 

374 

375 def get_dos(self): 

376 self.initialize() 

377 self.update() 

378 return self.dos_e 

379 

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] 

387 

388 def get_pdos(self): 

389 self.initialize() 

390 self.update() 

391 return self.pdos_ne 

392 

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) 

408 

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 

412 

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 

429 

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) 

444 

445 return rot_mm 

446 

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) 

452 

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

460 

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) 

469 

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) 

474 

475 return T_n, v_in