Coverage for kwant/physics/gauge.py : 97%
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
1# Copyright 2011-2019 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
8"""Functions for fixing the magnetic gauge automatically in a Kwant system.
11The "gauge" module has been included in Kwant on a provisional basis.
12Backwards incompatible changes (up to and including removal of the package)
13may occur if deemed necessary by the core developers.
14"""
16import bisect
17from functools import partial
18from itertools import permutations
20import numpy as np
21import scipy
22from scipy.integrate import dblquad
23from scipy.sparse import csgraph
25from .. import system, builder
27from ..graph.dijkstra import dijkstra_directed
29__all__ = ['magnetic_gauge']
32### Integation
34# Integrate vector field over triangle, for internal use by '_surface_integral'
35# Triangle is (origin, origin + v1, origin + v2), 'n' is np.cross(v1, v2)
37def _quad_triangle(f, origin, v1, v2, n, tol):
38 if np.dot(n, n) < tol**2: # does triangle have significant area? 38 ↛ 39line 38 didn't jump to line 39, because the condition on line 38 was never true
39 return 0
41 def g(x, y):
42 return np.dot(n, f(origin + x * v1 + y * v2))
44 result, *_ = dblquad(g, 0, 1, lambda x: 0, lambda x: 1-x)
45 return result.real
48def _const_triangle(f, origin, v1, v2, n, tol):
49 return np.dot(f, n) / 2
52def _average_triangle(f, origin, v1, v2, n, tol):
53 return np.dot(n, f(origin + 1/3 * (v1 + v2))) / 2
56def _surface_integral(f, loop, tol=1e-8, average=False):
57 """Calculate the surface integral of 'f' over a surface enclosed by 'loop'.
59 This function only works for *divergence free* vector fields, where the
60 surface integral depends only on the boundary.
62 Parameters
63 ----------
64 f : callable or real n-vector
65 The vector field for which to calculate the surface integral.
66 If callable, takes a real n-vector as argument and returns a
67 real n-vector.
68 loop : sequence of vectors
69 Ordered sequence of real n-vectors (positions) that define the
70 vertices of the polygon that encloses the surface to integrate over.
71 tol : float, default: 1e-8
72 Error tolerance on the result.
73 average : bool, default: False
74 If True, approximate the integral over each triangle using a
75 single function evaluation at the centre of the triangle.
76 """
77 if callable(f):
78 integrator = _average_triangle if average else _quad_triangle
79 else:
80 integrator = _const_triangle
82 origin, *points = loop
83 integral = 0
84 # Split loop into triangles with 1 vertex on 'origin', evaluate
85 # the integral over each triangle and sum the result
86 for p1, p2 in zip(points, points[1:]):
87 v1 = p1 - origin
88 v2 = p2 - origin
89 n = np.cross(v1, v2)
90 integral += integrator(f, origin, v1, v2, n, tol)
91 return integral
94### Loop finding graph algorithm
96def _find_loops(graph, subgraph):
97 """
98 Parameters
99 ----------
100 graph : COO matrix
101 The complete undirected graph, where the values of the matrix are
102 the weights of the corresponding graph links.
103 subgraph : COO matrix
104 An subgraph of 'graph', with missing edges denoted by infinities.
105 Must have the same sparsity structure as 'graph'.
107 Returns
108 -------
109 A sequence of paths which are partially contained in the subgraph.
110 The loop is formed by adding a link between the first and last node.
112 The ordering is such that the paths are made of links that belong to
113 the subgraph or to the previous closed loops.
114 """
115 # For each link we do 1 update of 'subgraph' and a call to
116 # 'csgraph.shortest_path'. It is cheaper to update the CSR
117 # matrix rather than convert to LIL and back every iteration.
118 subgraph = subgraph.tocsr()
119 graph = graph.tocsr()
120 assert _same_sparsity_structure(subgraph, graph)
122 # Links in graph, but not in subgraph.
123 links_to_find = scipy.sparse.triu(graph - subgraph).tocoo()
124 links_to_find = np.vstack((links_to_find.row, links_to_find.col)).transpose()
126 links_to_find, min_length = _order_links(subgraph, links_to_find)
128 # Find shortest path between each link in turn, updating the subgraph with
129 # the links as we go.
130 loops = []
131 while len(links_to_find) > 0:
132 frm, to = links_to_find[0]
133 (path,), (path_length,) = dijkstra_directed(subgraph,
134 sources=np.array([frm]),
135 targets=np.array([to]))
137 # Reorder links that are still to find based on the loop length in
138 # the updated graph. We only reorder when the path length for *this*
139 # link is a "little bit" longer that the perviously determined minimum.
140 # The "little bit" is needed so we don't needlessly re-order the links
141 # on amorphous lattices.
142 if path_length > min_length * 1.1:
143 links_to_find, min_length = _order_links(subgraph, links_to_find)
144 else:
145 # Assumes that 'graph' and 'subgraph' have the same sparsity structure.
146 _assign_csr(subgraph, graph, (frm, to))
147 _assign_csr(subgraph, graph, (to, frm))
148 loops.append(path)
149 links_to_find = links_to_find[1:]
151 return loops
154def _order_links(subgraph, links_to_find):
155 if len(links_to_find) == 0:
156 return [], None
157 # Order 'links_to_find' by length of shortest path between the nodes of the link
158 path_lengths = dijkstra_directed(subgraph,
159 sources=links_to_find[:, 0],
160 targets=links_to_find[:, 1],
161 return_paths=False)
162 idxs = np.argsort(path_lengths)
163 return links_to_find[idxs], path_lengths[idxs[0]]
166### Generic sparse matrix utilities
168def _assign_csr(a, b, element):
169 """Assign a single element from a CSR matrix to another.
171 Parameters
172 ----------
173 a : CSR matrix
174 b : CSR matrix or scalar
175 If a CSR matrix, must have the same sparsity structure
176 as 'a'. If a scalar, must be the same dtype as 'a'.
177 element: (int, int)
178 Row and column indices of the element to set.
179 """
180 assert isinstance(a, scipy.sparse.csr_matrix)
181 row, col = element
182 for j in range(a.indptr[row], a.indptr[row + 1]): 182 ↛ 186line 182 didn't jump to line 186, because the loop on line 182 didn't complete
183 if a.indices[j] == col:
184 break
185 else:
186 raise ValueError('{} not in sparse matrix'.format(element))
187 if isinstance(b, scipy.sparse.csr_matrix):
188 a.data[j] = b.data[j]
189 else:
190 a.data[j] = b
193def _same_sparsity_structure(a, b):
194 a = a.tocsr().sorted_indices()
195 b = b.tocsr().sorted_indices()
196 return (np.array_equal(a.indices, b.indices)
197 and np.array_equal(a.indptr, b.indptr))
200def _add_coo_matrices(*mats, shape):
201 """Add a sequence of COO matrices by appending their constituent arrays."""
202 values = np.hstack([mat.data for mat in mats])
203 rows = np.hstack([mat.row for mat in mats])
204 cols = np.hstack([mat.col for mat in mats])
205 return scipy.sparse.coo_matrix((values, (rows, cols)), shape=shape)
208def _shift_diagonally(mat, shift, shape):
209 """Shift the row/column indices of a COO matrix."""
210 return scipy.sparse.coo_matrix(
211 (mat.data, (mat.row + shift, mat.col + shift)),
212 shape=shape)
215def _distance_matrix(links, pos, shape):
216 """Return the distances between the provided links as a COO matrix.
218 Parameters
219 ----------
220 links : sequence of pairs of int
221 The links for which to find the lengths.
222 pos : callable: int -> vector
223 Map from link ends (integers) to realspace position.
224 shape : tuple
225 """
226 if len(links) == 0: # numpy does not like 'if array' 226 ↛ 227line 226 didn't jump to line 227, because the condition on line 226 was never true
227 return scipy.sparse.coo_matrix(shape)
228 links = np.array(links)
229 distances = np.array([pos(i) - pos(j) for i, j in links])
230 distances = np.linalg.norm(distances, axis=1)
231 return scipy.sparse.coo_matrix((distances, links.T), shape=shape)
234### Loop finding
235#
236# All of these functions take a finalized Kwant system and return
237# a sequence of loops. Each loop is a sequence of sites (integers)
238# that one visits when traversing the loop. The first and final sites
239# are assumed to be linked, which closes the loop. The links that one
240# traverses when going round a loop is thus:
241#
242# list(zip(loop, loop[1:])) + [(loop[-1], loop[0])]
243#
244# These loops are later used to fix the magnetic gauge in the system.
245# All of the links except the final one are assumed to have their gauge
246# fixed (i.e. the phase across them is known), and gauge of the final
247# link is the one to be determined.
250def _loops_in_finite(syst):
251 """Find the loops in a finite system with no leads.
253 The site indices in the returned loops are those of the system,
254 so they may be used as indices to 'syst.sites', or with 'syst.pos'.
255 """
256 assert isinstance(syst, system.FiniteSystem) and syst.leads == []
257 nsites = len(syst.sites)
259 # Fix the gauge across the minimum spanning tree of the system graph.
260 graph = _distance_matrix(list(syst.graph),
261 pos=syst.pos, shape=(nsites, nsites))
262 spanning_tree = _shortest_distance_forest(graph)
263 return _find_loops(graph, spanning_tree)
266def _shortest_distance_forest(graph):
267 # Grow a forest of minimum distance trees for all connected components of the graph
268 graph = graph.tocsr()
269 tree = graph.copy()
270 # set every entry in tree to infinity
271 tree.data[:] = np.inf
272 unvisited = set(range(graph.shape[0]))
273 # set the target node to be greater than any node in the graph.
274 # This way we explore the whole graph.
275 end = np.array([graph.shape[0] + 1], dtype=np.int32)
277 while unvisited:
278 # Choose an arbitrary element as the root
279 root = unvisited.pop()
280 root = np.array([root], dtype=np.int32)
281 _, pred = dijkstra_directed(graph, sources=root, targets=end,
282 return_predecessors=True, return_paths=False)
283 for i, p in enumerate(pred):
284 # -1 if node 'i' has no predecessor. Either it is the root node,
285 # or it was not reached.
286 if p != -1:
287 unvisited.remove(i)
288 _assign_csr(tree, graph, (i, p))
289 _assign_csr(tree, graph, (p, i))
290 return tree
293def _loops_in_infinite(syst):
294 """Find the loops in an infinite system.
296 Returns
297 -------
298 loops : sequence of sequences of integers
299 The sites in the returned loops belong to two adjacent unit
300 cells. The first 'syst.cell_size' sites are in the first
301 unit cell, and the next 'sys.cell_size' are in the next
302 (in the direction of the translational symmetry).
303 extended_sites : callable : int -> Site
304 Given a site index in the extended system consisting of
305 two unit cells, returns the associated high-level
306 `kwant.system.Site`.
307 """
308 assert isinstance(syst, system.InfiniteSystem)
309 _check_infinite_syst(syst)
311 cell_size = syst.cell_size
313 unit_cell_links = [(i, j) for i, j in syst.graph
314 if i < cell_size and j < cell_size]
315 unit_cell_graph = _distance_matrix(unit_cell_links,
316 pos=syst.pos,
317 shape=(cell_size, cell_size))
319 # Loops in the interior of the unit cell
320 spanning_tree = _shortest_distance_forest(unit_cell_graph)
321 loops = _find_loops(unit_cell_graph, spanning_tree)
323 # Construct an extended graph consisting of 2 unit cells connected
324 # by the inter-cell links.
325 extended_shape = (2 * cell_size, 2 * cell_size)
326 uc1 = _shift_diagonally(unit_cell_graph, 0, shape=extended_shape)
327 uc2 = _shift_diagonally(unit_cell_graph, cell_size, shape=extended_shape)
328 hop_links = [(i, j) for i, j in syst.graph if j >= cell_size]
329 hop = _distance_matrix(hop_links,
330 pos=syst.pos,
331 shape=extended_shape)
332 graph = _add_coo_matrices(uc1, uc2, hop, hop.T,
333 shape=extended_shape)
335 # Construct a subgraph where only the shortest link between the
336 # 2 unit cells is added. The other links are added with infinite
337 # values, so that the subgraph has the same sparsity structure
338 # as 'graph'.
339 idx = np.argmin(hop.data)
340 data = np.full_like(hop.data, np.inf)
341 data[idx] = hop.data[idx]
342 smallest_edge = scipy.sparse.coo_matrix(
343 (data, (hop.row, hop.col)),
344 shape=extended_shape)
345 subgraph = _add_coo_matrices(uc1, uc2, smallest_edge, smallest_edge.T,
346 shape=extended_shape)
348 # Use these two graphs to find the loops between unit cells.
349 loops.extend(_find_loops(graph, subgraph))
351 def extended_sites(i):
352 unit_cell = np.array([i // cell_size])
353 site = syst.sites[i % cell_size]
354 return syst.symmetry.act(-unit_cell, site)
356 return loops, extended_sites
359def _loops_in_composite(syst):
360 """Find the loops in finite system with leads.
362 Parameters
363 ----------
364 syst : kwant.builder.FiniteSystem
366 Returns
367 -------
368 loops : sequence of sequences of integers
369 The sites in each loop belong to the extended scattering region
370 (see notes). The first and last site in each loop are guaranteed
371 to be in 'syst'.
372 which_patch : callable : int -> int
373 Given a site index in the extended scattering region (see notes),
374 returns the lead patch (see notes) to which the site belongs. Returns
375 -1 if the site is part of the reduced scattering region (see notes).
376 extended_sites : callable : int -> Site
377 Given a site index in the extended scattering region (see notes),
378 returns the associated high-level `kwant.system.Site`.
380 Notes
381 -----
382 extended scattering region
383 The scattering region with a single lead unit cell attached at
384 each interface. This unit cell is added so that we can "see" any
385 loops formed with sites in the lead (see 'check_infinite_syst'
386 for details). The sites for each lead are added in the same
387 order as the leads, and within a given added unit cell the sites
388 are ordered in the same way as the associated lead.
389 lead patch
390 Sites in the extended scattering region that belong to the added
391 unit cell for a given lead, or the lead padding for a given lead
392 are said to be in the "lead patch" for that lead.
393 reduced scattering region
394 The sites of the extended scattering region that are not
395 in a lead patch.
396 """
397 # Check that we can consistently fix the gauge in the scattering region,
398 # given that we have independently fixed gauges in the leads.
399 _check_composite_system(syst)
401 # Get distance matrix for the extended scattering region,
402 # a function that maps sites to their lead patches (-1 for sites
403 # in the reduced scattering region), and a function that maps sites
404 # to high-level 'kwant.system.Site' objects.
405 distance_matrix, which_patch, extended_sites =\
406 _extended_scattering_region(syst)
408 spanning_tree = _spanning_tree_composite(distance_matrix, which_patch).tocsr()
410 # Fill in all links with at least 1 site in a lead patch;
411 # their gauge is fixed by the lead gauge.
412 for i, j, v in zip(distance_matrix.row, distance_matrix.col,
413 distance_matrix.data):
414 if which_patch(i) > -1 or which_patch(j) > -1:
415 _assign_csr(spanning_tree, v, (i, j))
416 _assign_csr(spanning_tree, v, (j, i))
418 loops = _find_loops(distance_matrix, spanning_tree)
420 return loops, which_patch, extended_sites
423def _extended_scattering_region(syst):
424 """Return the distance matrix of a finite system with 1 unit cell
425 added to each lead interface.
427 Parameters
428 ----------
429 syst : kwant.builder.FiniteSystem
431 Returns
432 -------
433 extended_scattering_region: COO matrix
434 Distance matrix between connected sites in the extended
435 scattering region.
436 which_patch : callable : int -> int
437 Given a site index in the extended scattering region, returns
438 the lead patch to which the site belongs. Returns
439 -1 if the site is part of the reduced scattering region.
440 extended_sites : callable : int -> Site
441 Given a site index in the extended scattering region, returns
442 the associated high-level `kwant.system.Site`.
444 Notes
445 -----
446 Definitions of the terms 'extended scatteringr region',
447 'lead patch' and 'reduced scattering region' are given
448 in the notes for `kwant.physics.gauge._loops_in_composite`.
449 """
450 extended_size = (syst.graph.num_nodes
451 + sum(l.cell_size for l in syst.leads))
452 extended_shape = (extended_size, extended_size)
454 added_unit_cells = []
455 first_lead_site = syst.graph.num_nodes
456 for lead, interface in zip(syst.leads, syst.lead_interfaces):
457 # Here we assume that the distance between sites in the added
458 # unit cell and sites in the interface is the same as between sites
459 # in neighboring unit cells.
460 uc = _distance_matrix(list(lead.graph),
461 pos=lead.pos, shape=extended_shape)
462 # Map unit cell lead sites to their indices in the extended scattering,
463 # region and sites in next unit cell to their interface sites.
464 hop_from_syst = uc.row >= lead.cell_size
465 uc.row[~hop_from_syst] = uc.row[~hop_from_syst] + first_lead_site
466 uc.row[hop_from_syst] = interface[uc.row[hop_from_syst] - lead.cell_size]
467 # Same for columns
468 hop_to_syst = uc.col >= lead.cell_size
469 uc.col[~hop_to_syst] = uc.col[~hop_to_syst] + first_lead_site
470 uc.col[hop_to_syst] = interface[uc.col[hop_to_syst] - lead.cell_size]
472 added_unit_cells.append(uc)
473 first_lead_site += lead.cell_size
475 scattering_region = _distance_matrix(list(syst.graph),
476 pos=syst.pos, shape=extended_shape)
478 extended_scattering_region = _add_coo_matrices(scattering_region,
479 *added_unit_cells,
480 shape=extended_shape)
482 lead_starts = np.cumsum([syst.graph.num_nodes,
483 *[lead.cell_size for lead in syst.leads]])
484 # Frozenset to quickly check 'is this site in the lead padding?'
485 extra_sites = [frozenset(sites) for sites in syst.lead_paddings]
488 def which_patch(i):
489 if i < len(syst.sites):
490 # In scattering region
491 for patch_num, sites in enumerate(extra_sites):
492 if i in sites:
493 return patch_num
494 # If not in 'extra_sites' it is in the reduced scattering region.
495 return -1
496 else:
497 # Otherwise it's in an attached lead cell
498 which_lead = bisect.bisect(lead_starts, i) - 1
499 assert which_lead > -1
500 return which_lead
503 # Here we use the fact that all the sites in a lead interface belong
504 # to the same symmetry domain.
505 interface_domains = [lead.symmetry.which(syst.sites[interface[0]])
506 for lead, interface in
507 zip(syst.leads, syst.lead_interfaces)]
509 def extended_sites(i):
510 if i < len(syst.sites):
511 # In scattering region
512 return syst.sites[i]
513 else:
514 # Otherwise it's in an attached lead cell
515 which_lead = bisect.bisect(lead_starts, i) - 1
516 assert which_lead > -1
517 lead = syst.leads[which_lead]
518 domain = interface_domains[which_lead] + 1
519 # Map extended scattering region site index to site index in lead.
520 i = i - lead_starts[which_lead]
521 return lead.symmetry.act(domain, lead.sites[i])
523 return extended_scattering_region, which_patch, extended_sites
526def _interior_links(distance_matrix, which_patch):
527 """Return the indices of the links in 'distance_matrix' that
528 connect interface sites of the scattering region to other
529 sites (interface and non-interface) in the scattering region.
530 """
532 def _is_in_lead(i):
533 return which_patch(i) > -1
535 # Sites that connect to/from sites in a lead patch
536 interface_sites = {
537 (i if not _is_in_lead(i) else j)
538 for i, j in zip(distance_matrix.row, distance_matrix.col)
539 if _is_in_lead(i) ^ _is_in_lead(j)
540 }
542 def _we_want(i, j):
543 return i in interface_sites and not _is_in_lead(j)
545 # Links that connect interface sites to the rest of the scattering region.
546 return np.array([
547 k
548 for k, (i, j) in enumerate(zip(distance_matrix.row, distance_matrix.col))
549 if _we_want(i, j) or _we_want(j, i)
550 ])
553def _make_metatree(graph, links_to_delete):
554 """Make a tree of the components of 'graph' that are
555 disconnected by deleting 'links'. The values of
556 the returned tree are indices of edges in 'graph'
557 that connect components.
558 """
559 # Partition the graph into disconnected components
560 dl = partial(np.delete, obj=links_to_delete)
561 partitioned_graph = scipy.sparse.coo_matrix(
562 (dl(graph.data), (dl(graph.row), dl(graph.col)))
563 )
564 # Construct the "metagraph", where each component is reduced to
565 # a single node, and a representative (smallest) edge is chosen
566 # among the edges that connected the components in the original graph.
567 ncc, labels = csgraph.connected_components(partitioned_graph)
568 metagraph = scipy.sparse.dok_matrix((ncc, ncc), int)
569 for k in links_to_delete:
570 i, j = labels[graph.row[k]], labels[graph.col[k]]
571 if i == j:
572 continue # Discard loop edges
573 # Add a representative (smallest) edge from each graph component.
574 if graph.data[k] < metagraph.get((i, j), np.inf):
575 metagraph[i, j] = k
576 metagraph[j, i] = k
578 return csgraph.minimum_spanning_tree(metagraph).astype(int)
581def _spanning_tree_composite(distance_matrix, which_patch):
582 """Find a spanning tree for a composite system.
584 We cannot use a simple minimum-distance spanning tree because
585 we have the additional constraint that all links with at least
586 one end in a lead patch have their gauge fixed. See the notes
587 for details.
589 Parameters
590 ----------
591 distance_matrix : COO matrix
592 Distance matrix between connected sites in the extended
593 scattering region.
594 which_patch : callable : int -> int
595 Given a site index in the extended scattering region (see notes),
596 returns the lead patch (see notes) to which the site belongs. Returns
597 -1 if the site is part of the reduced scattering region (see notes).
598 Returns
599 -------
600 spanning_tree : CSR matrix
601 A spanning tree with the same sparsity structure as 'distance_matrix',
602 where missing links are denoted with infinite weights.
604 Notes
605 -----
606 Definitions of the terms 'extended scattering region', 'lead patch'
607 and 'reduced scattering region' are given in the notes for
608 `kwant.physics.gauge._loops_in_composite`.
610 We cannot use a simple minimum-distance spanning tree because
611 we have the additional constraint that all links with at least
612 one end in a lead patch have their gauge fixed.
613 Consider the following case using a minimum-distance tree
614 where 'x' are sites in the lead patch::
616 o-o-x o-o-x
617 | | | --> | |
618 o-o-x o-o x
620 The removed link on the lower right comes from the lead, and hence
621 is gauge-fixed, however the vertical link in the center is not in
622 the lead, but *is* in the tree, which means that we will fix its
623 gauge to 0. The loop on the right would thus not have the correct
624 gauge on all links.
626 Instead we first cut all links between *interface* sites and
627 sites in the scattering region (including other interface sites).
628 We then construct a minimum distance forest for these disconnected
629 graphs. Finally we add back links from the ones that were cut,
630 ensuring that we do not form any loops; we do this by contructing
631 a tree of representative links from the "metagraph" of components
632 that were disconnected by the link cutting.
633 """
634 # Links that connect interface sites to other sites in the
635 # scattering region (including other interface sites)
636 links_to_delete = _interior_links(distance_matrix, which_patch)
637 # Make a shortest distance tree for each of the components
638 # obtained by cutting the links.
639 cut_syst = distance_matrix.copy()
640 cut_syst.data[links_to_delete] = np.inf
641 forest = _shortest_distance_forest(cut_syst)
642 # Connect the forest back up with representative links until
643 # we have a single tree (if the original system was not connected,
644 # we get a forest).
645 metatree = _make_metatree(distance_matrix, links_to_delete)
646 for k in np.unique(metatree.data):
647 value = distance_matrix.data[k]
648 i, j = distance_matrix.row[k], distance_matrix.col[k]
649 _assign_csr(forest, value, (i, j))
650 _assign_csr(forest, value, (j, i))
652 return forest
655def _check_infinite_syst(syst):
656 r"""Check that the unit cell is a connected graph.
658 If the unit cell is not connected then we cannot be sure whether
659 there are loops or not just by inspecting the unit cell graph
660 (this may be a solved problem, but we could not find an algorithm
661 to do this).
663 To illustrate this, consider the following unit cell consisting
664 of 3 sites and 4 hoppings::
666 o-
667 \
668 o
669 \
670 o-
672 None of the sites are connected within the unit cell, however if we repeat
673 a few unit cells::
675 o-o-o-o
676 \ \ \
677 o o o o
678 \ \ \
679 o-o-o-o
681 we see that there is a loop crossing 4 unit cells. A connected unit cell
682 is a sufficient condition that all the loops can be found by inspecting
683 the graph consisting of two unit cells glued together.
684 """
685 assert isinstance(syst, system.InfiniteSystem)
686 n = syst.cell_size
687 rows, cols = np.array([(i, j) for i, j in syst.graph
688 if i < n and j < n]).transpose()
689 data = np.ones(len(rows))
690 graph = scipy.sparse.coo_matrix((data, (rows, cols)), shape=(n, n))
691 if csgraph.connected_components(graph, return_labels=False) > 1: 691 ↛ 692line 691 didn't jump to line 692, because the condition on line 691 was never true
692 raise ValueError(
693 'Infinite system unit cell is not connected: we cannot determine '
694 'if there are any loops in the system\n\n'
695 'If there are, then you must define your unit cell so that it is '
696 'connected. If there are not, then you must add zero-magnitude '
697 'hoppings to your system.'
698 )
701def _check_composite_system(syst):
702 """Check that we can consistently fix the gauge in a system with leads.
704 If not, raise an exception with an informative error message.
705 """
706 assert isinstance(syst, system.FiniteSystem) and syst.leads
707 # Frozenset to quickly check 'is this site in the lead padding?'
708 extras = [frozenset(sites) for sites in syst.lead_paddings]
709 interfaces = [set(iface) for iface in syst.lead_interfaces]
710 # Make interfaces between lead patches and the reduced scattering region.
711 for interface, extra in zip(interfaces, extras):
712 extra_interface = set()
713 if extra:
714 extra_interface = set()
715 for i, j in syst.graph:
716 if i in extra and j not in extra:
717 extra_interface.add(j)
718 interface -= extra
719 interface |= extra_interface
720 assert not extra.intersection(interface)
722 pre_msg = (
723 'Attaching leads results in gauge-fixed loops in the extended '
724 'scattering region (scattering region plus one lead unit cell '
725 'from every lead). This does not allow consistent gauge-fixing.\n\n'
726 )
727 solution_msg = (
728 'To avoid this error, attach leads further away from each other.\n\n'
729 'Note: calling `attach_lead()` with `add_cells > 0` will not fix '
730 'this problem, as the added sites inherit the gauge from the lead. '
731 'To extend the scattering region, you must manually add sites '
732 'making sure that they use the scattering region gauge.'
733 )
735 # Check that there is at most one overlapping site between
736 # reduced interface of one lead and extra sites of another
737 num_leads = len(syst.leads)
738 metagraph = scipy.sparse.lil_matrix((num_leads, num_leads))
739 for i, j in permutations(range(num_leads), 2):
740 intersection = len(interfaces[i] & (interfaces[j] | extras[j]))
741 if intersection > 1:
742 raise ValueError(
743 pre_msg
744 + ('There is at least one gauge-fixed loop in the overlap '
745 'of leads {} and {}.\n\n'.format(i, j))
746 + solution_msg
747 )
748 elif intersection == 1:
749 metagraph[i, j] = 1
750 # Check that there is no loop formed by gauge-fixed bonds of multiple leads.
751 num_components = scipy.sparse.csgraph.connected_components(metagraph, return_labels=False)
752 if metagraph.nnz // 2 + num_components != num_leads:
753 raise ValueError(
754 pre_msg
755 + ('There is at least one gauge-fixed loop formed by more than 2 leads. '
756 ' The connectivity matrix of the leads is:\n\n'
757 '{}\n\n'.format(metagraph.A))
758 + solution_msg
759 )
761### Phase calculation
763def _calculate_phases(loops, pos, previous_phase, flux):
764 """Calculate the phase across the terminal links of a set of loops
766 Parameters
767 ----------
768 loops : sequence of sequences of int
769 The loops over which to calculate the flux. We wish to find the phase
770 over the link connecting the first and last sites in each loop.
771 The phase across all other links in a given loop is assumed known.
772 pos : callable : int -> ndarray
773 A map from site (integer) to realspace position.
774 previous_phase : callable
775 Takes a dict that maps from links to phases, and a loop and returns
776 the product of the phases across each link in the loop,
777 *except* the link between the first and last site in the loop.
778 flux : callable
779 Takes a sequence of positions and returns the magnetic flux through the
780 surface defined by the provided loop.
782 Returns
783 -------
784 phases : dict : (int, int) -> float
785 A map from links to the phase across those links.
786 """
787 phases = dict()
788 for loop in loops:
789 tail, head = loop[-1], loop[0]
790 integral = flux([pos(p) for p in loop])
791 phase = np.exp(1j * np.pi * integral)
792 phases[tail, head] = phase / previous_phase(phases, loop)
793 return phases
796# These functions are to be used with '_calculate_phases'.
797# 'phases' always stores *either* the phase across (i, j) *or*
798# (j, i), and never both. If a phase is not present it is assumed to
799# be zero.
802def _previous_phase_finite(phases, path):
803 previous_phase = 1
804 for i, j in zip(path, path[1:]):
805 previous_phase *= phases.get((i, j), 1)
806 previous_phase /= phases.get((j, i), 1)
807 return previous_phase
810def _previous_phase_infinite(cell_size, phases, path):
811 previous_phase = 1
812 for i, j in zip(path, path[1:]):
813 # i and j are only in the fundamental unit cell (0 <= i < cell_size)
814 # or the next one (cell_size <= i < 2 * cell_size).
815 if i >= cell_size and j >= cell_size:
816 assert i // cell_size == j // cell_size
817 i = i % cell_size
818 j = j % cell_size
819 previous_phase *= phases.get((i, j), 1)
820 previous_phase /= phases.get((j, i), 1)
821 return previous_phase
824def _previous_phase_composite(which_patch, extended_sites, lead_phases,
825 phases, path):
826 previous_phase = 1
827 for i, j in zip(path, path[1:]):
828 patch_i = which_patch(i)
829 patch_j = which_patch(j)
830 if patch_i == -1 and patch_j == -1:
831 # Both sites in reduced scattering region.
832 previous_phase *= phases.get((i, j), 1)
833 previous_phase /= phases.get((j, i), 1)
834 else:
835 # At least one site in a lead patch; use the phase from the
836 # associated lead. Check that if both are in a patch, they
837 # are in the same patch.
838 assert patch_i * patch_j <= 0 or patch_i == patch_j
839 patch = max(patch_i, patch_j)
840 a, b = extended_sites(i), extended_sites(j)
841 previous_phase *= lead_phases[patch](a, b)
843 return previous_phase
846### High-level interface
847#
848# These functions glue all the above functionality together.
850# Wrapper for phase dict that takes high-level sites
851def _finite_wrapper(syst, phases, a, b):
852 i = syst.id_by_site[a]
853 j = syst.id_by_site[b]
854 # We only store *either* (i, j) *or* (j, i). If not present
855 # then the phase is unity by definition.
856 if (i, j) in phases:
857 return phases[i, j]
858 elif (j, i) in phases:
859 return phases[j, i].conjugate()
860 else:
861 return 1
864def _infinite_wrapper(syst, phases, a, b):
865 sym = syst.symmetry
866 # Bring link to fundamental domain consistently with how
867 # we store the phases.
868 t = max(sym.which(a), sym.which(b))
869 a, b = sym.act(-t, a, b)
870 i = syst.id_by_site[a]
871 j = syst.id_by_site[b]
872 # We only store *either* (i, j) *or* (j, i). If not present
873 # then the phase is unity by definition.
874 if (i, j) in phases:
875 return phases[i, j]
876 elif (j, i) in phases:
877 return phases[j, i].conjugate()
878 else:
879 return 1
882def _peierls_finite(syst, loops, syst_field, tol, average):
883 integrate = partial(_surface_integral, syst_field,
884 tol=tol, average=average)
885 phases = _calculate_phases(
886 loops,
887 syst.pos,
888 _previous_phase_finite,
889 integrate,
890 )
891 return partial(_finite_wrapper, syst, phases)
894def _peierls_infinite(syst, loops, extended_sites, syst_field, tol, average):
895 integrate = partial(_surface_integral, syst_field,
896 tol=tol, average=average)
897 phases = _calculate_phases(
898 loops,
899 lambda i: extended_sites(i).pos,
900 partial(_previous_phase_infinite, syst.cell_size),
901 integrate,
902 )
903 return partial(_infinite_wrapper, syst, phases)
906def _peierls_composite(syst, loops, which_patch, extended_sites, lead_gauges,
907 syst_field, *lead_fields, tol, average):
908 if len(lead_fields) != len(syst.leads): 908 ↛ 909line 908 didn't jump to line 909, because the condition on line 908 was never true
909 raise ValueError('Magnetic fields must be provided for all leads.')
911 lead_phases = [gauge(B, tol=tol, average=average)
912 for gauge, B in zip(lead_gauges, lead_fields)]
914 flux = partial(_surface_integral, syst_field, tol=tol, average=average)
916 # NOTE: uses the scattering region magnetic field to set the phase
917 # of the inteface hoppings this choice is somewhat arbitrary,
918 # but it is consistent with the position defined in the scattering
919 # region coordinate system. the integrate functions for the leads
920 # may be defined far from the interface.
921 phases = _calculate_phases(
922 loops,
923 lambda i: extended_sites(i).pos,
924 partial(_previous_phase_composite,
925 which_patch, extended_sites, lead_phases),
926 flux,
927 )
929 return (partial(_finite_wrapper, syst, phases), *lead_phases)
932# This class is essentially a closure, but documenting closures is a pain.
933# To emphasise the lack of manipulatable or inspectable state, we name the
934# class as we would for a function.
936class magnetic_gauge:
937 """Fix the magnetic gauge for a finalized system.
939 This can be used to later calculate the Peierls phases that
940 should be applied to each hopping, given a magnetic field.
942 This API is currently provisional. Refer to the documentation
943 for details.
945 Parameters
946 ----------
947 syst : `kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`
949 Examples
950 --------
951 The following illustrates basic usage for a scattering region with
952 a single lead attached:
954 >>> import numpy as np
955 >>> import kwant
956 >>>
957 >>> def hopping(a, b, t, peierls):
958 >>> return -t * peierls(a, b)
959 >>>
960 >>> syst = make_system(hopping)
961 >>> lead = make_lead(hopping)
962 >>> lead = lead.substituted(peierls='peierls_lead')
963 >>> syst.attach_lead(lead)
964 >>> syst = syst.finalized()
965 >>>
966 >>> gauge = kwant.physics.magnetic_gauge(syst)
967 >>>
968 >>> def B_syst(pos):
969 >>> return np.exp(-np.sum(pos * pos))
970 >>>
971 >>> peierls_syst, peierls_lead = gauge(B_syst, 0)
972 >>>
973 >>> params = dict(t=1, peierls=peierls_syst, peierls_lead=peierls_lead)
974 >>> kwant.hamiltonian_submatrix(syst, params=params)
975 """
977 def __init__(self, syst):
978 if isinstance(syst, builder.FiniteSystem):
979 if syst.leads:
980 loops, which_patch, extended_sites = _loops_in_composite(syst)
981 lead_gauges = [magnetic_gauge(lead) for lead in syst.leads]
982 self._peierls = partial(_peierls_composite, syst,
983 loops, which_patch,
984 extended_sites, lead_gauges)
985 else:
986 loops = _loops_in_finite(syst)
987 self._peierls = partial(_peierls_finite, syst, loops)
988 elif isinstance(syst, builder.InfiniteSystem): 988 ↛ 993line 988 didn't jump to line 993, because the condition on line 988 was never false
989 loops, extended_sites = _loops_in_infinite(syst)
990 self._peierls = partial(_peierls_infinite, syst,
991 loops, extended_sites)
992 else:
993 raise TypeError('Expected a finalized Builder')
995 def __call__(self, syst_field, *lead_fields, tol=1E-8, average=False):
996 """Return the Peierls phase for a magnetic field configuration.
998 Parameters
999 ----------
1000 syst_field : scalar, vector or callable
1001 The magnetic field to apply to the scattering region.
1002 If callable, takes a position and returns the
1003 magnetic field at that position. Can be a scalar if
1004 the system is 1D or 2D, otherwise must be a vector.
1005 Magnetic field is expressed in units :math:`φ₀ / l²`,
1006 where :math:`φ₀` is the magnetic flux quantum and
1007 :math:`l` is the unit of length.
1008 *lead_fields : scalar, vector or callable
1009 The magnetic fields to apply to each of the leads, in
1010 the same format as 'syst_field'. In addition, if a callable
1011 is provided, then the magnetic field must have the symmetry
1012 of the associated lead.
1013 tol : float, default: 1E-8
1014 The tolerance to which to calculate the flux through each
1015 hopping loop in the system.
1016 average : bool, default: False
1017 If True, estimate the magnetic flux through each hopping loop
1018 in the system by evaluating the magnetic field at a single
1019 position inside the loop and multiplying it by the area of the
1020 loop. If False, then ``scipy.integrate.quad`` is used to integrate
1021 the magnetic field. This parameter is only used when 'syst_field'
1022 or 'lead_fields' are callable.
1024 Returns
1025 -------
1026 phases : callable, or sequence of callables
1027 The first callable computes the Peierls phase in the scattering
1028 region and the remaining callables compute the Peierls phases
1029 in each of the leads. Each callable takes a pair of
1030 `~kwant.system.Site` (a hopping) and returns a unit complex
1031 number (Peierls phase) that multiplies that hopping.
1032 """
1033 return self._peierls(syst_field, *lead_fields, tol=tol, average=False)