Coverage for /builds/kinetik161/ase/ase/transport/stm.py: 10.00%

1# flake8: noqa 

2import time 

3 

4import numpy as np 

5 

6from ase.parallel import world 

7from ase.transport.greenfunction import GreenFunction 

8from ase.transport.selfenergy import LeadSelfEnergy 

9from ase.transport.tools import dagger 

10 

11 

12class STM: 

13 def __init__(self, h1, s1, h2, s2, h10, s10, h20, s20, 

14 eta1, eta2, w=0.5, pdos=[], logfile=None): 

15 """XXX 

16 

17 1. Tip 

18 2. Surface 

19 

20 h1: ndarray 

21 Hamiltonian and overlap matrix for the isolated tip 

22 calculation. Note, h1 should contain (at least) one 

23 principal layer. 

24 

25 h2: ndarray 

26 Same as h1 but for the surface. 

27 

28 h10: ndarray 

29 periodic part of the tip. must include two and only 

30 two principal layers. 

31 

32 h20: ndarray 

33 same as h10, but for the surface 

34 

35 The s* are the corresponding overlap matrices. eta1, and eta 

36 2 are (finite) infinitesimals. """ 

37 

38 self.pl1 = len(h10) // 2 # principal layer size for the tip 

39 self.pl2 = len(h20) // 2 # principal layer size for the surface 

40 self.h1 = h1 

41 self.s1 = s1 

42 self.h2 = h2 

43 self.s2 = s2 

44 self.h10 = h10 

45 self.s10 = s10 

46 self.h20 = h20 

47 self.s20 = s20 

48 self.eta1 = eta1 

49 self.eta2 = eta2 

50 self.w = w # asymmetry of the applied bias (0.5=>symmetric) 

51 self.pdos = [] 

52 self.log = logfile 

53 

54 def initialize(self, energies, bias=0): 

55 """ 

56 energies: list of energies 

57 for which the transmission function should be evaluated. 

58 bias. 

59 Will precalculate the surface greenfunctions of the tip and 

60 surface. 

61 """ 

62 self.bias = bias 

63 self.energies = energies 

64 nenergies = len(energies) 

65 pl1, pl2 = self.pl1, self.pl2 

66 nbf1, nbf2 = len(self.h1), len(self.h2) 

67 

68 # periodic part of the tip 

69 hs1_dii = self.h10[:pl1, :pl1], self.s10[:pl1, :pl1] 

70 hs1_dij = self.h10[:pl1, pl1:2 * pl1], self.s10[:pl1, pl1:2 * pl1] 

71 # coupling between per. and non. per part of the tip 

72 h1_im = np.zeros((pl1, nbf1), complex) 

73 s1_im = np.zeros((pl1, nbf1), complex) 

74 h1_im[:pl1, :pl1], s1_im[:pl1, :pl1] = hs1_dij 

75 hs1_dim = [h1_im, s1_im] 

76 

77 # periodic part the surface 

78 hs2_dii = self.h20[:pl2, :pl2], self.s20[:pl2, :pl2] 

79 hs2_dij = self.h20[pl2:2 * pl2, :pl2], self.s20[pl2:2 * pl2, :pl2] 

80 # coupling between per. and non. per part of the surface 

81 h2_im = np.zeros((pl2, nbf2), complex) 

82 s2_im = np.zeros((pl2, nbf2), complex) 

83 h2_im[-pl2:, -pl2:], s2_im[-pl2:, -pl2:] = hs2_dij 

84 hs2_dim = [h2_im, s2_im] 

85 

86 # tip and surface greenfunction 

87 self.selfenergy1 = LeadSelfEnergy(hs1_dii, hs1_dij, hs1_dim, self.eta1) 

88 self.selfenergy2 = LeadSelfEnergy(hs2_dii, hs2_dij, hs2_dim, self.eta2) 

89 self.greenfunction1 = GreenFunction(self.h1 - self.bias * self.w * self.s1, self.s1, 

90 [self.selfenergy1], self.eta1) 

91 self.greenfunction2 = GreenFunction(self.h2 - self.bias * (self.w - 1) * self.s2, self.s2, 

92 [self.selfenergy2], self.eta2) 

93 

94 # Shift the bands due to the bias. 

95 bias_shift1 = -bias * self.w 

96 bias_shift2 = -bias * (self.w - 1) 

97 self.selfenergy1.set_bias(bias_shift1) 

98 self.selfenergy2.set_bias(bias_shift2) 

99 

100 # tip and surface greenfunction matrices. 

101 nbf1_small = nbf1 # XXX Change this for efficiency in the future 

102 nbf2_small = nbf2 # XXX -||- 

103 coupling_list1 = list(range(nbf1_small)) # XXX -||- 

104 coupling_list2 = list(range(nbf2_small)) # XXX -||- 

105 self.gft1_emm = np.zeros((nenergies, nbf1_small, nbf1_small), complex) 

106 self.gft2_emm = np.zeros((nenergies, nbf2_small, nbf2_small), complex) 

107 

108 for e, energy in enumerate(self.energies): 

109 if self.log is not None: # and world.rank == 0: 

110 T = time.localtime() 

111 self.log.write(' %d:%02d:%02d, ' % (T[3], T[4], T[5]) + 

112 '%d, %d, %02f\n' % (world.rank, e, energy)) 

113 gft1_mm = self.greenfunction1.retarded(energy)[coupling_list1] 

114 gft1_mm = np.take(gft1_mm, coupling_list1, axis=1) 

115 

116 gft2_mm = self.greenfunction2.retarded(energy)[coupling_list2] 

117 gft2_mm = np.take(gft2_mm, coupling_list2, axis=1) 

118 

119 self.gft1_emm[e] = gft1_mm 

120 self.gft2_emm[e] = gft2_mm 

121 

122 if self.log is not None and world.rank == 0: 

123 self.log.flush() 

124 

125 def get_transmission(self, v_12, v_11_2=None, v_22_1=None): 

126 """XXX 

127 

128 v_12: 

129 coupling between tip and surface 

130 v_11_2: 

131 correction to "on-site" tip elements due to the 

132 surface (eq.16). Is only included to first order. 

133 v_22_1: 

134 corretion to "on-site" surface elements due to he 

135 tip (eq.17). Is only included to first order. 

136 """ 

137 

138 dim0 = v_12.shape[0] 

139 dim1 = v_12.shape[1] 

140 

141 nenergies = len(self.energies) 

142 T_e = np.empty(nenergies, float) 

143 v_21 = dagger(v_12) 

144 for e, energy in enumerate(self.energies): 

145 gft1 = self.gft1_emm[e] 

146 if v_11_2 is not None: 

147 gf1 = np.dot(v_11_2, np.dot(gft1, v_11_2)) 

148 gf1 += gft1 # eq. 16 

149 else: 

150 gf1 = gft1 

151 

152 gft2 = self.gft2_emm[e] 

153 if v_22_1 is not None: 

154 gf2 = np.dot(v_22_1, np.dot(gft2, v_22_1)) 

155 gf2 += gft2 # eq. 17 

156 else: 

157 gf2 = gft2 

158 

159 a1 = (gf1 - dagger(gf1)) 

160 a2 = (gf2 - dagger(gf2)) 

161 self.v_12 = v_12 

162 self.a2 = a2 

163 self.v_21 = v_21 

164 self.a1 = a1 

165 v12_a2 = np.dot(v_12, a2[:dim1]) 

166 v21_a1 = np.dot(v_21, a1[-dim0:]) 

167 self.v12_a2 = v12_a2 

168 self.v21_a1 = v21_a1 

169 T = -np.trace(np.dot(v12_a2[:, :dim1], v21_a1[:, -dim0:])) # eq. 11 

170 assert abs(T.imag).max() < 1e-14 

171 T_e[e] = T.real 

172 self.T_e = T_e 

173 return T_e 

174 

175 def get_current(self, bias, v_12, v_11_2=None, v_22_1=None): 

176 """Very simple function to calculate the current. 

177 

178 Asummes zero temperature. 

179 

180 bias: type? XXX 

181 bias voltage (V) 

182 

183 v_12: XXX 

184 coupling between tip and surface. 

185 

186 v_11_2: 

187 correction to onsite elements of the tip 

188 due to the potential of the surface. 

189 v_22_1: 

190 correction to onsite elements of the surface 

191 due to the potential of the tip. 

192 """ 

193 energies = self.energies 

194 T_e = self.get_transmission(v_12, v_11_2, v_22_1) 

195 bias_window = sorted(-np.array([bias * self.w, bias * (self.w - 1)])) 

196 self.bias_window = bias_window 

197 # print 'bias window', np.around(bias_window,3) 

198 # print 'Shift of tip lead do to the bias:', self.selfenergy1.bias 

199 # print 'Shift of surface lead do to the bias:', self.selfenergy2.bias 

200 i1 = sum(energies < bias_window[0]) 

201 i2 = sum(energies < bias_window[1]) 

202 step = 1 

203 if i2 < i1: 

204 step = -1 

205 

206 return np.sign(bias) * \ 

207 np.trapz(x=energies[i1:i2:step], y=T_e[i1:i2:step])