1# flake8: noqa 

2import time 

3 

4import numpy as np 

5 

6from ase.transport.tools import dagger 

7from ase.transport.selfenergy import LeadSelfEnergy 

8from ase.transport.greenfunction import GreenFunction 

9from ase.parallel import world 

10 

11 

12class STM: 

13 def __init__(self, h1, s1, h2, s2 ,h10, s10, h20, s20, eta1, eta2, w=0.5, pdos=[], logfile = None): 

14 """XXX 

15  

16 1. Tip 

17 2. Surface 

18  

19 h1: ndarray 

20 Hamiltonian and overlap matrix for the isolated tip 

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

22 principal layer. 

23 

24 h2: ndarray 

25 Same as h1 but for the surface. 

26 

27 h10: ndarray 

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

29 two principal layers. 

30 

31 h20: ndarray 

32 same as h10, but for the surface 

33 

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

35 2 are (finite) infinitesimals. """ 

36 

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

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

39 self.h1 = h1 

40 self.s1 = s1 

41 self.h2 = h2 

42 self.s2 = s2 

43 self.h10 = h10 

44 self.s10 = s10 

45 self.h20 = h20 

46 self.s20 = s20 

47 self.eta1 = eta1 

48 self.eta2 = eta2 

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

50 self.pdos = [] 

51 self.log = logfile 

52 

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

54 """ 

55 energies: list of energies  

56 for which the transmission function should be evaluated. 

57 bias. 

58 Will precalculate the surface greenfunctions of the tip and 

59 surface. 

60 """ 

61 self.bias = bias 

62 self.energies = energies 

63 nenergies = len(energies) 

64 pl1, pl2 = self.pl1, self.pl2 

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

66 

67 #periodic part of the tip 

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

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

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

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

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

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

74 hs1_dim = [h1_im, s1_im] 

75 

76 #periodic part the surface  

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

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

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

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

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

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

83 hs2_dim = [h2_im, s2_im] 

84 

85 #tip and surface greenfunction  

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

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

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

89 [self.selfenergy1], self.eta1) 

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

91 [self.selfenergy2], self.eta2) 

92 

93 #Shift the bands due to the bias. 

94 bias_shift1 = -bias * self.w 

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

96 self.selfenergy1.set_bias(bias_shift1) 

97 self.selfenergy2.set_bias(bias_shift2) 

98 

99 #tip and surface greenfunction matrices. 

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

101 nbf2_small = nbf2 #XXX -||- 

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

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

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

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

106 

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

108 if self.log != None: # and world.rank == 0: 

109 T = time.localtime() 

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

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

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

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

114 

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

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

117 

118 self.gft1_emm[e] = gft1_mm 

119 self.gft2_emm[e] = gft2_mm 

120 

121 if self.log != None and world.rank == 0: 

122 self.log.flush() 

123 

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

125 """XXX 

126 

127 v_12: 

128 coupling between tip and surface  

129 v_11_2: 

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

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

132 v_22_1: 

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

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

135 """ 

136 

137 dim0 = v_12.shape[0] 

138 dim1 = v_12.shape[1] 

139 

140 nenergies = len(self.energies) 

141 T_e = np.empty(nenergies,float) 

142 v_21 = dagger(v_12) 

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

144 gft1 = self.gft1_emm[e] 

145 if v_11_2!=None: 

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

147 gf1 += gft1 #eq. 16 

148 else: 

149 gf1 = gft1 

150 

151 gft2 = self.gft2_emm[e] 

152 if v_22_1!=None: 

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

154 gf2 += gft2 #eq. 17 

155 else: 

156 gf2 = gft2 

157 

158 a1 = (gf1 - dagger(gf1)) 

159 a2 = (gf2 - dagger(gf2)) 

160 self.v_12 = v_12 

161 self.a2 = a2 

162 self.v_21 = v_21 

163 self.a1 = a1 

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

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

166 self.v12_a2 = v12_a2 

167 self.v21_a1 = v21_a1 

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

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

170 T_e[e] = T.real 

171 self.T_e = T_e 

172 return T_e 

173 

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)*np.trapz(x=energies[i1:i2:step], y=T_e[i1:i2:step]) 

207 

208 

209 

210 

211 

212