Coverage for /builds/kinetik161/ase/ase/dft/stm.py: 82.72%
162 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import numpy as np
3from ase.io.jsonio import read_json, write_json
6class STM:
7 def __init__(self, atoms, symmetries=None, use_density=False):
8 """Scanning tunneling microscope.
10 atoms: Atoms object or filename
11 Atoms to scan or name of file to read LDOS from.
12 symmetries: list of int
13 List of integers 0, 1, and/or 2 indicating which surface
14 symmetries have been used to reduce the number of k-points
15 for the DFT calculation. The three integers correspond to
16 the following three symmetry operations::
18 [-1 0] [ 1 0] [ 0 1]
19 [ 0 1] [ 0 -1] [ 1 0]
21 use_density: bool
22 Use the electron density instead of the LDOS.
23 """
25 self.use_density = use_density
27 if isinstance(atoms, str):
28 with open(atoms) as fd:
29 self.ldos, self.bias, self.cell = read_json(fd,
30 always_array=False)
31 self.atoms = None
32 else:
33 self.atoms = atoms
34 self.cell = atoms.cell
35 self.bias = None
36 self.ldos = None
37 assert not self.cell[2, :2].any() and not self.cell[:2, 2].any()
39 self.symmetries = symmetries or []
41 def calculate_ldos(self, bias):
42 """Calculate local density of states for given bias."""
43 if self.ldos is not None and bias == self.bias:
44 return
46 self.bias = bias
48 calc = self.atoms.calc
50 if self.use_density:
51 self.ldos = calc.get_pseudo_density()
52 return
54 if bias < 0:
55 emin = bias
56 emax = 0.0
57 else:
58 emin = 0
59 emax = bias
61 nbands = calc.get_number_of_bands()
62 weights = calc.get_k_point_weights()
63 nkpts = len(weights)
64 nspins = calc.get_number_of_spins()
65 eigs = np.array([[calc.get_eigenvalues(k, s)
66 for k in range(nkpts)]
67 for s in range(nspins)])
68 eigs -= calc.get_fermi_level()
69 ldos = np.zeros(calc.get_pseudo_wave_function(0, 0, 0).shape)
71 for s in range(nspins):
72 for k in range(nkpts):
73 for n in range(nbands):
74 e = eigs[s, k, n]
75 if emin < e < emax:
76 psi = calc.get_pseudo_wave_function(n, k, s)
77 ldos += weights[k] * (psi * np.conj(psi)).real
79 if 0 in self.symmetries:
80 # (x,y) -> (-x,y)
81 ldos[1:] += ldos[:0:-1].copy()
82 ldos[1:] *= 0.5
84 if 1 in self.symmetries:
85 # (x,y) -> (x,-y)
86 ldos[:, 1:] += ldos[:, :0:-1].copy()
87 ldos[:, 1:] *= 0.5
89 if 2 in self.symmetries:
90 # (x,y) -> (y,x)
91 ldos += ldos.transpose((1, 0, 2)).copy()
92 ldos *= 0.5
94 self.ldos = ldos
96 def write(self, filename):
97 """Write local density of states to JSON file."""
98 write_json(filename, (self.ldos, self.bias, self.cell))
100 def get_averaged_current(self, bias, z):
101 """Calculate avarage current at height z (in Angstrom).
103 Use this to get an idea of what current to use when scanning."""
105 self.calculate_ldos(bias)
106 nz = self.ldos.shape[2]
108 # Find grid point:
109 n = z / self.cell[2, 2] * nz
110 dn = n - np.floor(n)
111 n = int(n) % nz
113 # Average and do linear interpolation:
114 return ((1 - dn) * self.ldos[:, :, n].mean() +
115 dn * self.ldos[:, :, (n + 1) % nz].mean())
117 def scan(self, bias, current, z0=None, repeat=(1, 1)):
118 """Constant current 2-d scan.
120 Returns three 2-d arrays (x, y, z) containing x-coordinates,
121 y-coordinates and heights. These three arrays can be passed to
122 matplotlibs contourf() function like this:
124 >>> import matplotlib.pyplot as plt
125 >>> plt.contourf(x, y, z)
126 >>> plt.show()
128 """
130 self.calculate_ldos(bias)
132 L = self.cell[2, 2]
133 nz = self.ldos.shape[2]
134 h = L / nz
136 ldos = self.ldos.reshape((-1, nz))
138 heights = np.empty(ldos.shape[0])
139 for i, a in enumerate(ldos):
140 heights[i] = find_height(a, current, h, z0)
142 s0 = heights.shape = self.ldos.shape[:2]
143 heights = np.tile(heights, repeat)
144 s = heights.shape
146 ij = np.indices(s, dtype=float).reshape((2, -1)).T
147 x, y = np.dot(ij / s0, self.cell[:2, :2]).T.reshape((2,) + s)
149 return x, y, heights
151 def scan2(self, bias, z, repeat=(1, 1)):
152 """Constant height 2-d scan.
154 Returns three 2-d arrays (x, y, I) containing x-coordinates,
155 y-coordinates and currents. These three arrays can be passed to
156 matplotlibs contourf() function like this:
158 >>> import matplotlib.pyplot as plt
159 >>> plt.contourf(x, y, I)
160 >>> plt.show()
162 """
164 self.calculate_ldos(bias)
166 nz = self.ldos.shape[2]
167 ldos = self.ldos.reshape((-1, nz))
169 current = np.empty(ldos.shape[0])
171 zp = z / self.cell[2, 2] * nz
172 zp = int(zp) % nz
174 for i, a in enumerate(ldos):
175 current[i] = self.find_current(a, zp)
177 s0 = current.shape = self.ldos.shape[:2]
178 current = np.tile(current, repeat)
179 s = current.shape
181 ij = np.indices(s, dtype=float).reshape((2, -1)).T
182 x, y = np.dot(ij / s0, self.cell[:2, :2]).T.reshape((2,) + s)
184 # Returing scan with axes in Angstrom.
185 return x, y, current
187 def linescan(self, bias, current, p1, p2, npoints=50, z0=None):
188 """Constant current line scan.
190 Example::
192 stm = STM(...)
193 z = ... # tip position
194 c = stm.get_averaged_current(-1.0, z)
195 stm.linescan(-1.0, c, (1.2, 0.0), (1.2, 3.0))
196 """
198 heights = self.scan(bias, current, z0)[2]
200 p1 = np.asarray(p1, float)
201 p2 = np.asarray(p2, float)
202 d = p2 - p1
203 s = np.dot(d, d)**0.5
205 cell = self.cell[:2, :2]
206 shape = np.array(heights.shape, float)
207 M = np.linalg.inv(cell)
208 line = np.empty(npoints)
209 for i in range(npoints):
210 p = p1 + i * d / (npoints - 1)
211 q = np.dot(p, M) * shape
212 line[i] = interpolate(q, heights)
213 return np.linspace(0, s, npoints), line
215 def pointcurrent(self, bias, x, y, z):
216 """Current for a single x, y, z position for a given bias."""
218 self.calculate_ldos(bias)
220 nx = self.ldos.shape[0]
221 ny = self.ldos.shape[1]
222 nz = self.ldos.shape[2]
224 # Find grid point:
225 xp = x / np.linalg.norm(self.cell[0]) * nx
226 dx = xp - np.floor(xp)
227 xp = int(xp) % nx
229 yp = y / np.linalg.norm(self.cell[1]) * ny
230 dy = yp - np.floor(yp)
231 yp = int(yp) % ny
233 zp = z / np.linalg.norm(self.cell[2]) * nz
234 dz = zp - np.floor(zp)
235 zp = int(zp) % nz
237 # 3D interpolation of the LDOS at point (x,y,z) at given bias.
238 xyzldos = (((1 - dx) + (1 - dy) + (1 - dz)) * self.ldos[xp, yp, zp] +
239 dx * self.ldos[(xp + 1) % nx, yp, zp] +
240 dy * self.ldos[xp, (yp + 1) % ny, zp] +
241 dz * self.ldos[xp, yp, (zp + 1) % nz])
243 return dos2current(bias, xyzldos)
245 def sts(self, x, y, z, bias0, bias1, biasstep):
246 """Returns the dI/dV curve for position x, y at height z (in Angstrom),
247 for bias from bias0 to bias1 with step biasstep."""
249 biases = np.arange(bias0, bias1 + biasstep, biasstep)
250 current = np.zeros(biases.shape)
252 for b in np.arange(len(biases)):
253 print(b, biases[b])
254 current[b] = self.pointcurrent(biases[b], x, y, z)
256 dIdV = np.gradient(current, biasstep)
258 return biases, current, dIdV
260 def find_current(self, ldos, z):
261 """ Finds current for given LDOS at height z."""
262 nz = self.ldos.shape[2]
264 zp = z / self.cell[2, 2] * nz
265 dz = zp - np.floor(zp)
266 zp = int(zp) % nz
268 ldosz = (1 - dz) * ldos[zp] + dz * ldos[(zp + 1) % nz]
270 return dos2current(self.bias, ldosz)
273def dos2current(bias, dos):
274 # Borrowed from gpaw/analyse/simple_stm.py:
275 # The connection between density n and current I
276 # n [e/Angstrom^3] = 0.0002 sqrt(I [nA])
277 # as given in Hofer et al., RevModPhys 75 (2003) 1287
278 return 5000. * dos**2 * (1 if bias > 0 else -1)
281def interpolate(q, heights):
282 qi = q.astype(int)
283 f = q - qi
284 g = 1 - f
285 qi %= heights.shape
286 n0, m0 = qi
287 n1, m1 = (qi + 1) % heights.shape
288 z = (g[0] * g[1] * heights[n0, m0] +
289 f[0] * g[1] * heights[n1, m0] +
290 g[0] * f[1] * heights[n0, m1] +
291 f[0] * f[1] * heights[n1, m1])
292 return z
295def find_height(ldos, current, h, z0=None):
296 if z0 is None:
297 n = len(ldos) - 2
298 else:
299 n = int(z0 / h)
300 while n >= 0:
301 if ldos[n] > current:
302 break
303 n -= 1
304 else:
305 return 0.0
307 c2, c1 = ldos[n:n + 2]
308 return (n + 1 - (current - c1) / (c2 - c1)) * h
311def delta(biases, bias, width):
312 """Return a delta-function centered at 'bias'"""
313 x = -((biases - bias) / width)**2
314 return np.exp(x) / (np.sqrt(np.pi) * width)