Source code for ase.transport.calculators

import numpy as np
from numpy import linalg

from ase.transport.greenfunction import GreenFunction
from ase.transport.selfenergy import BoxProbe, LeadSelfEnergy
from ase.transport.tools import (cutcoupling, dagger, fermidistribution,
                                 rotate_matrix, subdiagonalize)
from ase.units import kB


[docs]class TransportCalculator: """Determine transport properties of a device sandwiched between two semi-infinite leads using a Green function method. """ def __init__(self, **kwargs): """Create the transport calculator. Parameters: h : (N, N) ndarray Hamiltonian matrix for the central region. s : {None, (N, N) ndarray}, optional Overlap matrix for the central region. Use None for an orthonormal basis. h1 : (N1, N1) ndarray Hamiltonian matrix for lead1. h2 : {None, (N2, N2) ndarray}, optional Hamiltonian matrix for lead2. You may use None if lead1 and lead2 are identical. s1 : {None, (N1, N1) ndarray}, optional Overlap matrix for lead1. Use None for an orthonomormal basis. hc1 : {None, (N1, N) ndarray}, optional Hamiltonian coupling matrix between the first principal layer in lead1 and the central region. hc2 : {None, (N2, N} ndarray), optional Hamiltonian coupling matrix between the first principal layer in lead2 and the central region. sc1 : {None, (N1, N) ndarray}, optional Overlap coupling matrix between the first principal layer in lead1 and the central region. sc2 : {None, (N2, N) ndarray}, optional Overlap coupling matrix between the first principal layer in lead2 and the central region. energies : {None, array_like}, optional Energy points for which calculated transport properties are evaluated. eta : {1.0e-5, float}, optional Infinitesimal for the central region Green function. eta1/eta2 : {1.0e-5, float}, optional Infinitesimal for lead1/lead2 Green function. align_bf : {None, int}, optional Use align_bf=m to shift the central region by a constant potential such that the m'th onsite element in the central region is aligned to the m'th onsite element in lead1 principal layer. logfile : {None, str}, optional Write a logfile to file with name `logfile`. Use '-' to write to std out. eigenchannels: {0, int}, optional Number of eigenchannel transmission coefficients to calculate. pdos : {None, (N,) array_like}, optional Specify which basis functions to calculate the projected density of states for. dos : {False, bool}, optional The total density of states of the central region. box: XXX YYY If hc1/hc2 are None, they are assumed to be identical to the coupling matrix elements between neareste neighbor principal layers in lead1/lead2. Examples: >>> import numpy as np >>> h = np.array((0,)).reshape((1,1)) >>> h1 = np.array((0, -1, -1, 0)).reshape(2,2) >>> energies = np.arange(-3, 3, 0.1) >>> calc = TransportCalculator(h=h, h1=h1, energies=energies) >>> T = calc.get_transmission() """ # The default values for all extra keywords self.input_parameters = {'energies': None, 'h': None, 'h1': None, 'h2': None, 's': None, 's1': None, 's2': None, 'hc1': None, 'hc2': None, 'sc1': None, 'sc2': None, 'box': None, 'align_bf': None, 'eta1': 1e-5, 'eta2': 1e-5, 'eta': 1e-5, 'logfile': None, 'eigenchannels': 0, 'dos': False, 'pdos': []} self.initialized = False # Changed Hamiltonians? self.uptodate = False # Changed energy grid? self.set(**kwargs) def set(self, **kwargs): for key in kwargs: if key in ['h', 'h1', 'h2', 'hc1', 'hc2', 's', 's1', 's2', 'sc1', 'sc2', 'eta', 'eta1', 'eta2', 'align_bf', 'box']: self.initialized = False self.uptodate = False break elif key in ['energies', 'eigenchannels', 'dos', 'pdos']: self.uptodate = False elif key not in self.input_parameters: raise KeyError('%r not a vaild keyword' % key) self.input_parameters.update(kwargs) log = self.input_parameters['logfile'] if log is None: class Trash: def write(self, s): pass def flush(self): pass self.log = Trash() elif log == '-': from sys import stdout self.log = stdout elif 'logfile' in kwargs: self.log = open(log, 'w') def initialize(self): if self.initialized: return print('# Initializing calculator...', file=self.log) p = self.input_parameters if p['s'] is None: p['s'] = np.identity(len(p['h'])) identical_leads = False if p['h2'] is None: p['h2'] = p['h1'] # Lead2 is idendical to lead1 identical_leads = True if p['s1'] is None: p['s1'] = np.identity(len(p['h1'])) if identical_leads: p['s2'] = p['s1'] else: if p['s2'] is None: p['s2'] = np.identity(len(p['h2'])) h_mm = p['h'] s_mm = p['s'] pl1 = len(p['h1']) // 2 pl2 = len(p['h2']) // 2 h1_ii = p['h1'][:pl1, :pl1] h1_ij = p['h1'][:pl1, pl1:2 * pl1] s1_ii = p['s1'][:pl1, :pl1] s1_ij = p['s1'][:pl1, pl1:2 * pl1] h2_ii = p['h2'][:pl2, :pl2] h2_ij = p['h2'][pl2: 2 * pl2, :pl2] s2_ii = p['s2'][:pl2, :pl2] s2_ij = p['s2'][pl2: 2 * pl2, :pl2] if p['hc1'] is None: nbf = len(h_mm) h1_im = np.zeros((pl1, nbf), complex) s1_im = np.zeros((pl1, nbf), complex) h1_im[:pl1, :pl1] = h1_ij s1_im[:pl1, :pl1] = s1_ij p['hc1'] = h1_im p['sc1'] = s1_im else: h1_im = p['hc1'] if p['sc1'] is not None: s1_im = p['sc1'] else: s1_im = np.zeros(h1_im.shape, complex) p['sc1'] = s1_im if p['hc2'] is None: h2_im = np.zeros((pl2, nbf), complex) s2_im = np.zeros((pl2, nbf), complex) h2_im[-pl2:, -pl2:] = h2_ij s2_im[-pl2:, -pl2:] = s2_ij p['hc2'] = h2_im p['sc2'] = s2_im else: h2_im = p['hc2'] if p['sc2'] is not None: s2_im = p['sc2'] else: s2_im = np.zeros(h2_im.shape, complex) p['sc2'] = s2_im align_bf = p['align_bf'] if align_bf is not None: diff = ((h_mm[align_bf, align_bf] - h1_ii[align_bf, align_bf]) / s_mm[align_bf, align_bf]) print('# Aligning scat. H to left lead H. diff=', diff, file=self.log) h_mm -= diff * s_mm # Setup lead self-energies # All infinitesimals must be > 0 assert np.all(np.array((p['eta'], p['eta1'], p['eta2'])) > 0.0) self.selfenergies = [LeadSelfEnergy((h1_ii, s1_ii), (h1_ij, s1_ij), (h1_im, s1_im), p['eta1']), LeadSelfEnergy((h2_ii, s2_ii), (h2_ij, s2_ij), (h2_im, s2_im), p['eta2'])] box = p['box'] if box is not None: print('Using box probe!') self.selfenergies.append( BoxProbe(eta=box[0], a=box[1], b=box[2], energies=box[3], S=s_mm, T=0.3)) # setup scattering green function self.greenfunction = GreenFunction(selfenergies=self.selfenergies, H=h_mm, S=s_mm, eta=p['eta']) self.initialized = True def update(self): if self.uptodate: return p = self.input_parameters self.energies = p['energies'] nepts = len(self.energies) nchan = p['eigenchannels'] pdos = p['pdos'] self.T_e = np.empty(nepts) if p['dos']: self.dos_e = np.empty(nepts) if pdos != []: self.pdos_ne = np.empty((len(pdos), nepts)) if nchan > 0: self.eigenchannels_ne = np.empty((nchan, nepts)) for e, energy in enumerate(self.energies): Ginv_mm = self.greenfunction.retarded(energy, inverse=True) lambda1_mm = self.selfenergies[0].get_lambda(energy) lambda2_mm = self.selfenergies[1].get_lambda(energy) a_mm = linalg.solve(Ginv_mm, lambda1_mm) b_mm = linalg.solve(dagger(Ginv_mm), lambda2_mm) T_mm = np.dot(a_mm, b_mm) if nchan > 0: t_n = linalg.eigvals(T_mm).real self.eigenchannels_ne[:, e] = np.sort(t_n)[-nchan:] self.T_e[e] = np.sum(t_n) else: self.T_e[e] = np.trace(T_mm).real print(energy, self.T_e[e], file=self.log) self.log.flush() if p['dos']: self.dos_e[e] = self.greenfunction.dos(energy) if pdos != []: self.pdos_ne[:, e] = np.take(self.greenfunction.pdos(energy), pdos) self.uptodate = True def print_pl_convergence(self): self.initialize() pl1 = len(self.input_parameters['h1']) // 2 h_ii = self.selfenergies[0].h_ii s_ii = self.selfenergies[0].s_ii ha_ii = self.greenfunction.H[:pl1, :pl1] sa_ii = self.greenfunction.S[:pl1, :pl1] c1 = np.abs(h_ii - ha_ii).max() c2 = np.abs(s_ii - sa_ii).max() print(f'Conv (h,s)={c1:.2e}, {c2:2.e}') def plot_pl_convergence(self): self.initialize() pl1 = len(self.input_parameters['h1']) // 2 hlead = self.selfenergies[0].h_ii.real.diagonal() hprincipal = self.greenfunction.H.real.diagonal[:pl1] import pylab as pl pl.plot(hlead, label='lead') pl.plot(hprincipal, label='principal layer') pl.axis('tight') pl.show() def get_current(self, bias, T=0., E=None, T_e=None, spinpol=False): '''Returns the current as a function of the bias voltage. **Parameters:** bias : {float, (M,) ndarray}, units: V Specifies the bias voltage. T : {float}, units: K, optional Specifies the temperature. E : {(N,) ndarray}, units: eV, optional Contains energy grid of the transmission function. T_e {(N,) ndarray}, units: unitless, optional Contains the transmission function. spinpol: {bool}, optional Specifies whether the current should be calculated assuming degenerate spins **Returns:** I : {float, (M,) ndarray}, units: 2e/h*eV Contains the electric current. Examples: >> import numpy as np >> import pylab as plt >> from ase import units >> >> bias = np.arange(0, 2, .1) >> current = calc.get_current(bias, T = 0.) >> plt.plot(bias, 2.*units._e**2/units._hplanck*current) >> plt.xlabel('U [V]') >> plt.ylabel('I [A]') >> plt.show() ''' if E is not None: if T_e is None: self.energies = E self.uptodate = False T_e = self.get_transmission().copy() else: assert self.uptodate, ('Energy grid and transmission function ' 'not defined.') E = self.energies.copy() T_e = self.T_e.copy() if not isinstance(bias, (int, float)): bias = bias[np.newaxis] E = E[:, np.newaxis] T_e = T_e[:, np.newaxis] fl = fermidistribution(E - bias / 2., kB * T) fr = fermidistribution(E + bias / 2., kB * T) if spinpol: return .5 * np.trapz((fl - fr) * T_e, x=E, axis=0) else: return np.trapz((fl - fr) * T_e, x=E, axis=0) def get_transmission(self): self.initialize() self.update() return self.T_e def get_dos(self): self.initialize() self.update() return self.dos_e def get_eigenchannels(self, n=None): """Get ``n`` first eigenchannels.""" self.initialize() self.update() if n is None: n = self.input_parameters['eigenchannels'] return self.eigenchannels_ne[:n] def get_pdos(self): self.initialize() self.update() return self.pdos_ne def subdiagonalize_bfs(self, bfs, apply=False): self.initialize() bfs = np.array(bfs) p = self.input_parameters h_mm = p['h'] s_mm = p['s'] ht_mm, st_mm, c_mm, e_m = subdiagonalize(h_mm, s_mm, bfs) if apply: self.uptodate = False h_mm[:] = ht_mm.real s_mm[:] = st_mm.real # Rotate coupling between lead and central region for alpha, sigma in enumerate(self.selfenergies): sigma.h_im[:] = np.dot(sigma.h_im, c_mm) sigma.s_im[:] = np.dot(sigma.s_im, c_mm) c_mm = np.take(c_mm, bfs, axis=0) c_mm = np.take(c_mm, bfs, axis=1) return ht_mm, st_mm, e_m.real, c_mm def cutcoupling_bfs(self, bfs, apply=False): self.initialize() bfs = np.array(bfs) p = self.input_parameters h_pp = p['h'].copy() s_pp = p['s'].copy() cutcoupling(h_pp, s_pp, bfs) if apply: self.uptodate = False p['h'][:] = h_pp p['s'][:] = s_pp for alpha, sigma in enumerate(self.selfenergies): for m in bfs: sigma.h_im[:, m] = 0.0 sigma.s_im[:, m] = 0.0 return h_pp, s_pp def lowdin_rotation(self, apply=False): p = self.input_parameters h_mm = p['h'] s_mm = p['s'] eig, rot_mm = linalg.eigh(s_mm) eig = np.abs(eig) rot_mm = np.dot(rot_mm / np.sqrt(eig), dagger(rot_mm)) if apply: self.uptodate = False h_mm[:] = rotate_matrix(h_mm, rot_mm) # rotate C region s_mm[:] = rotate_matrix(s_mm, rot_mm) for alpha, sigma in enumerate(self.selfenergies): sigma.h_im[:] = np.dot(sigma.h_im, rot_mm) # rotate L-C coupl. sigma.s_im[:] = np.dot(sigma.s_im, rot_mm) return rot_mm def get_left_channels(self, energy, nchan=1): self.initialize() g_s_ii = self.greenfunction.retarded(energy) lambda_l_ii = self.selfenergies[0].get_lambda(energy) lambda_r_ii = self.selfenergies[1].get_lambda(energy) if self.greenfunction.S is not None: s_mm = self.greenfunction.S s_s_i, s_s_ii = linalg.eig(s_mm) s_s_i = np.abs(s_s_i) s_s_sqrt_i = np.sqrt(s_s_i) # sqrt of eigenvalues s_s_sqrt_ii = np.dot(s_s_ii * s_s_sqrt_i, dagger(s_s_ii)) s_s_isqrt_ii = np.dot(s_s_ii / s_s_sqrt_i, dagger(s_s_ii)) lambdab_r_ii = np.dot(np.dot(s_s_isqrt_ii, lambda_r_ii), s_s_isqrt_ii) a_l_ii = np.dot(np.dot(g_s_ii, lambda_l_ii), dagger(g_s_ii)) ab_l_ii = np.dot(np.dot(s_s_sqrt_ii, a_l_ii), s_s_sqrt_ii) lambda_i, u_ii = linalg.eig(ab_l_ii) ut_ii = np.sqrt(lambda_i / (2.0 * np.pi)) * u_ii m_ii = 2 * np.pi * np.dot(np.dot(dagger(ut_ii), lambdab_r_ii), ut_ii) T_i, c_in = linalg.eig(m_ii) T_i = np.abs(T_i) channels = np.argsort(-T_i)[:nchan] c_in = np.take(c_in, channels, axis=1) T_n = np.take(T_i, channels) v_in = np.dot(np.dot(s_s_isqrt_ii, ut_ii), c_in) return T_n, v_in