Coverage for kwant/operator.pyx : 94%
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"""Tools for working with operators for acting on wavefunctions."""
10__all__ = ['Density', 'Current', 'Source']
12import cython
13import itertools
14import functools as ft
15import collections
16import warnings
18import numpy as np
19import tinyarray as ta
20from scipy.sparse import coo_matrix
22from libc cimport math
24cdef extern from "complex.h":
25 double cabs(double complex)
27from .graph.core cimport EdgeIterator
28from .graph.core import DisabledFeatureError, NodeDoesNotExistError
29from .graph.defs cimport gint
30from .graph.defs import gint_dtype
31from .system import (
32 is_infinite, is_vectorized, Site, SiteArray, _normalize_matrix_blocks
33)
34from ._common import (
35 UserCodeError, KwantDeprecationWarning, get_parameters, deprecate_args
36)
39################ Generic Utility functions
41@cython.boundscheck(False)
42@cython.wraparound(False)
43cdef gint _bisect(gint[:] a, gint x):
44 "bisect.bisect specialized for searching `site_ranges`"
45 cdef gint mid, lo = 0, hi = a.shape[0]
46 while lo < hi:
47 mid = (lo + hi) // 2
48 if x < a[mid]:
49 hi = mid
50 else:
51 lo = mid + 1
52 return lo
55@cython.boundscheck(False)
56@cython.wraparound(False)
57cdef int _is_hermitian(
58 complex[:, :] a, double atol=1e-300, double rtol=1e-13
59) except -1:
60 "Return True if 'a' is Hermitian"
62 if a.shape[0] != a.shape[1]:
63 return False
65 # compute max(a)
66 cdef double tmp, max_a = 0
67 cdef gint i, j, k
68 for i in range(a.shape[0]):
69 for j in range(a.shape[1]):
70 tmp = cabs(a[i, j])
71 if tmp > max_a:
72 max_a = tmp
73 max_a = math.sqrt(max_a)
75 cdef double tol = rtol * max_a + atol
76 for i in range(a.shape[0]):
77 for j in range(i, a.shape[1]):
78 tmp = cabs(a[i, j] - a[j, i].conjugate())
79 if tmp > tol:
80 return False
81 return True
83@cython.boundscheck(False)
84@cython.wraparound(False)
85cdef int _is_hermitian_3d(
86 complex[:, :, :] a, double atol=1e-300, double rtol=1e-13
87) except -1:
88 "Return True if 'a' is Hermitian"
90 if a.shape[1] != a.shape[2]:
91 return False
93 # compute max(a)
94 cdef double tmp, max_a = 0
95 cdef gint i, j, k
96 for k in range(a.shape[0]):
97 for i in range(a.shape[1]):
98 for j in range(a.shape[2]):
99 tmp = cabs(a[k, i, j])
100 if tmp > max_a:
101 max_a = tmp
102 max_a = math.sqrt(max_a)
104 cdef double tol = rtol * max_a + atol
105 for k in range(a.shape[0]):
106 for i in range(a.shape[1]):
107 for j in range(i, a.shape[2]):
108 tmp = cabs(a[k, i, j] - a[k, j, i].conjugate())
109 if tmp > tol:
110 return False
111 return True
114@cython.boundscheck(False)
115@cython.wraparound(False)
116cdef _select(gint[:, :] arr, gint[:] indexes):
117 ret = np.empty((indexes.shape[0], arr.shape[1]), dtype=gint_dtype)
118 cdef gint[:, :] ret_view = ret
119 cdef gint i, j
121 for i in range(indexes.shape[0]):
122 for j in range(arr.shape[1]):
123 ret_view[i, j] = arr[indexes[i], j]
125 return ret
127################ Helper functions
129_shape_msg = ('{0} matrix dimensions do not match '
130 'the declared number of orbitals')
132_herm_msg = ('{0} matrix is not hermitian, use the option '
133 '`check_hermiticity=True` if this is intentional.')
135cdef int _check_onsite(complex[:, :] M, gint norbs,
136 int check_hermiticity) except -1:
137 "Check onsite matrix for correct shape and hermiticity."
138 if M.shape[0] != M.shape[1]:
139 raise UserCodeError('Onsite matrix is not square')
140 if M.shape[0] != norbs:
141 raise UserCodeError(_shape_msg.format('Onsite'))
142 if check_hermiticity and not _is_hermitian(M):
143 raise ValueError(_herm_msg.format('Onsite'))
144 return 0
147cdef int _check_onsites(complex[:, :, :] M, gint norbs,
148 int check_hermiticity) except -1:
149 "Check onsite matrix for correct shape and hermiticity."
150 if M.shape[1] != M.shape[2]:
151 raise UserCodeError('Onsite matrix is not square')
152 if M.shape[1] != norbs:
153 raise UserCodeError(_shape_msg.format('Onsite'))
154 if check_hermiticity and not _is_hermitian_3d(M):
155 raise ValueError(_herm_msg.format('Onsite'))
156 return 0
159cdef int _check_hams(complex[:, :, :] H, gint to_norbs, gint from_norbs,
160 int check_hermiticity) except -1:
161 if H.shape[1] != to_norbs or H.shape[2] != from_norbs:
162 raise UserCodeError(_shape_msg.format('Hamiltonian'))
163 if check_hermiticity and not _is_hermitian_3d(H):
164 raise ValueError(_herm_msg.format('Hamiltonian'))
165 return 0
168def _check_norbs(syst):
169 # TODO: Update this when it becomes clear how ND systems will be
170 # implemented.
171 if syst.site_ranges is None:
172 raise ValueError('Number of orbitals not defined.\n'
173 'Declare the number of orbitals using the '
174 '`norbs` keyword argument when constructing '
175 'the site families (lattices).')
178@cython.boundscheck(False)
179@cython.wraparound(False)
180cdef void _get_orbs(gint[:, :] site_ranges, gint site,
181 gint *start_orb, gint *norbs):
182 """Return the first orbital of this site and the number of orbitals"""
183 cdef gint run_idx, first_site, norb, orb_offset, orb
184 # Calculate the index of the range that contains the site.
185 run_idx = _bisect(site_ranges[:, 0], site) - 1
186 first_site = site_ranges[run_idx, 0]
187 norb = site_ranges[run_idx, 1]
188 orb_offset = site_ranges[run_idx, 2]
189 # calculate the slice
190 start_orb[0] = orb_offset + (site - first_site) * norb
191 norbs[0] = norb
194@cython.boundscheck(False)
195@cython.wraparound(False)
196def _get_all_orbs(gint[:, :] where, gint[:, :] site_ranges):
197 cdef gint[:, :] offsets = np.empty((where.shape[0], 2), dtype=gint_dtype)
198 cdef gint[:, :] norbs = np.empty((where.shape[0], 2), dtype=gint_dtype)
200 cdef gint w, a, a_offset, a_norbs, b, b_offset, b_norbs
201 for w in range(where.shape[0]):
202 a = where[w, 0]
203 _get_orbs(site_ranges, a, &a_offset, &a_norbs)
204 if where.shape[1] == 1:
205 b, b_offset, b_norbs = a, a_offset, a_norbs
206 else:
207 b = where[w, 1]
208 _get_orbs(site_ranges, b, &b_offset, &b_norbs)
209 offsets[w, 0] = a_offset
210 offsets[w, 1] = b_offset
211 norbs[w, 0] = a_norbs
212 norbs[w, 1] = b_norbs
214 return offsets, norbs
217def _get_tot_norbs(syst):
218 cdef gint _unused, tot_norbs
219 _check_norbs(syst)
220 is_infinite_system = is_infinite(syst)
221 n_sites = syst.cell_size if is_infinite_system else syst.graph.num_nodes
222 _get_orbs(np.asarray(syst.site_ranges, dtype=gint_dtype),
223 n_sites, &tot_norbs, &_unused)
224 return tot_norbs
227def _normalize_site_where(syst, where):
228 """Normalize the format of `where` when `where` contains sites.
230 If `where` is None, then all sites in the system are returned.
231 If it is a general sequence then it is expanded into an array. If `syst`
232 is a finalized Builder then `where` may contain `Site` objects,
233 otherwise it should contain integers.
234 """
235 if where is None:
236 if is_infinite(syst):
237 where = list(range(syst.cell_size))
238 else:
239 where = list(range(syst.graph.num_nodes))
240 elif callable(where):
241 try:
242 where = [syst.id_by_site[s] for s in filter(where, syst.sites)]
243 except AttributeError:
244 if is_infinite(syst):
245 where = [s for s in range(syst.cell_size) if where(s)]
246 else:
247 where = [s for s in range(syst.graph.num_nodes) if where(s)]
248 else:
249 if isinstance(where[0], Site):
250 try:
251 where = [syst.id_by_site[s] for s in where]
252 except AttributeError:
253 raise TypeError("'where' contains Sites, but the system is not "
254 "a finalized Builder.")
256 where = np.asarray(where, dtype=gint_dtype).reshape(-1, 1)
258 if is_infinite(syst) and np.any(where >= syst.cell_size):
259 raise ValueError('Only sites in the fundamental domain may be '
260 'specified using `where`.')
261 if np.any(np.logical_or(where < 0, where >= syst.graph.num_nodes)):
262 raise ValueError('`where` contains sites that are not in the '
263 'system.')
265 return where
268def _normalize_hopping_where(syst, where):
269 """Normalize the format of `where` when `where` contains hoppings.
271 If `where` is None, then all hoppings in the system are returned.
272 If it is a general iterator then it is expanded into an array. If `syst` is
273 a finalized Builder then `where` may contain pairs of `Site` objects,
274 otherwise it should contain pairs of integers.
275 """
276 if where is None:
277 # we cannot extract the hoppings in the same order as they are in the
278 # graph while simultaneously excluding all inter-cell hoppings
279 if is_infinite(syst):
280 raise ValueError('`where` must be provided when calculating '
281 'current in an InfiniteSystem.')
282 where = list(syst.graph)
283 elif callable(where):
284 if hasattr(syst, "sites"):
285 def idxwhere(hop):
286 a, b = hop
287 return where(syst.sites[a], syst.sites[b])
288 where = list(filter(idxwhere, syst.graph))
289 else:
290 where = list(filter(lambda h: where(*h), syst.graph))
291 else:
292 if isinstance(where[0][0], Site):
293 try:
294 where = list((syst.id_by_site[a], syst.id_by_site[b])
295 for a, b in where)
296 except AttributeError:
297 raise TypeError("'where' contains Sites, but the system is not "
298 "a finalized Builder.")
299 # NOTE: if we ever have operators that contain elements that are
300 # not in the system graph, then we should modify this check
301 try:
302 error = ValueError('`where` contains hoppings that are not '
303 'in the system.')
304 if any(not syst.graph.has_edge(*w) for w in where):
305 raise error
306 # If where contains: negative integers, or integers > # of sites
307 except (NodeDoesNotExistError, DisabledFeatureError):
308 raise error
310 where = np.asarray(where, dtype=gint_dtype)
312 if is_infinite(syst) and np.any(where > syst.cell_size):
313 raise ValueError('Only intra-cell hoppings may be specified '
314 'using `where`.')
316 return where
319## These four classes are here to avoid using closures, as these will
320## break pickle support. These are only used inside '_normalize_onsite'.
322def _raise_user_error(exc, func, vectorized):
323 msg = [
324 'Error occurred in user-supplied onsite function "{0}".',
325 'Did you remember to vectorize "{0}"?' if vectorized else '',
326 'See the upper part of the above backtrace for more information.',
327 ]
328 msg = '\n'.join(line for line in msg if line).format(func.__name__)
329 raise UserCodeError(msg) from exc
332class _FunctionalOnsite:
334 def __init__(self, onsite, sites, site_ranges):
335 self.onsite = onsite
336 self.sites = sites
337 self.site_ranges = site_ranges
339 def __call__(self, site_range, site_offsets, *args):
340 sites = self.sites
341 offset, norbs, _ = self.site_ranges[site_range]
342 try:
343 ret = [self.onsite(sites[offset + i], *args) for i in site_offsets]
344 except Exception as exc:
345 _raise_user_error(exc, self.onsite, vectorized=False)
347 expected_shape = (len(site_offsets), norbs, norbs)
348 return _normalize_matrix_blocks(ret, expected_shape,
349 calling_function=self.onsite)
352class _VectorizedFunctionalOnsite:
354 def __init__(self, onsite, site_arrays):
355 self.onsite = onsite
356 self.site_arrays = site_arrays
358 def __call__(self, site_range, site_offsets, *args):
359 site_array = self.site_arrays[site_range]
360 tags = site_array.tags[site_offsets]
361 sites = SiteArray(site_array.family, tags)
362 try:
363 ret = self.onsite(sites, *args)
364 except Exception as exc:
365 _raise_user_error(exc, self.onsite, vectorized=True)
367 expected_shape = (len(sites), sites.family.norbs, sites.family.norbs)
368 return _normalize_matrix_blocks(ret, expected_shape,
369 calling_function=self.onsite)
372class _FunctionalOnsiteNoTransform:
374 def __init__(self, onsite, site_ranges):
375 self.onsite = onsite
376 self.site_ranges = site_ranges
378 def __call__(self, site_range, site_offsets, *args):
379 offset, norbs, _ = self.site_ranges[site_range]
380 site_ids = offset + site_offsets
381 try:
382 ret = [self.onsite(id, *args) for id in site_ids]
383 except Exception as exc:
384 _raise_user_error(exc, self.onsite, vectorized=False)
386 expected_shape = (len(site_offsets), norbs, norbs)
387 return _normalize_matrix_blocks(ret, expected_shape,
388 calling_function=self.onsite)
391class _DictOnsite:
393 def __init__(self, onsite, range_family_map):
394 self.onsite = onsite
395 self.range_family_map = range_family_map
397 def __call__(self, site_range, site_offsets, *args):
398 fam = self.range_family_map[site_range]
399 ret = [self.onsite[fam]] * len(site_offsets)
401 expected_shape = (len(site_offsets), fam.norbs, fam.norbs)
402 return _normalize_matrix_blocks(ret, expected_shape)
405def _normalize_onsite(syst, onsite, check_hermiticity):
406 """Normalize the format of `onsite`.
408 If `onsite` is a function or a mapping (dictionary) then a function
409 is returned.
410 """
411 param_names = ()
413 if callable(onsite):
414 # make 'onsite' compatible with hamiltonian value functions
415 param_names = get_parameters(onsite)[1:]
416 if is_vectorized(syst):
417 _onsite = _VectorizedFunctionalOnsite(onsite, syst.site_arrays)
418 elif hasattr(syst, "sites"): # probably a non-vectorized finalized Builder
419 _onsite = _FunctionalOnsite(onsite, syst.sites, syst.site_ranges)
420 else: # generic old-style system, therefore *not* vectorized.
421 _onsite = _FunctionalOnsiteNoTransform(onsite, syst.site_ranges)
423 elif isinstance(onsite, collections.abc.Mapping):
424 # onsites known; immediately check for correct shape and hermiticity
425 for fam, _onsite in onsite.items():
426 _onsite = ta.matrix(_onsite, complex)
427 _check_onsite(_onsite, fam.norbs, check_hermiticity)
429 if is_vectorized(syst):
430 range_family_map = [arr.family for arr in syst.site_arrays]
431 elif not hasattr(syst, 'sites'):
432 raise TypeError('Provide `onsite` as a value or a function for '
433 'systems that are not finalized Builders.')
434 else:
435 # The last entry in 'site_ranges' is just an end marker, so we remove it
436 range_family_map = [syst.sites[r[0]].family for r in syst.site_ranges[:-1]]
437 _onsite = _DictOnsite(onsite, range_family_map)
439 else:
440 # single onsite; immediately check for correct shape and hermiticity
441 _onsite = ta.matrix(onsite, complex)
442 _check_onsite(_onsite, _onsite.shape[0], check_hermiticity)
443 if _onsite.shape[0] == 1:
444 # NOTE: this is wasteful when many orbitals per site, but it
445 # simplifies the code in `_operate`. If this proves to be a
446 # bottleneck, then we can add a code path for scalar onsites
447 max_norbs = max(norbs for (_, norbs, _) in syst.site_ranges)
448 _onsite = _onsite[0, 0] * ta.identity(max_norbs, complex)
449 elif len(set(norbs for _, norbs, _ in syst.site_ranges[:-1])) == 1:
450 # we have the same number of orbitals everywhere
451 norbs = syst.site_ranges[0][1]
452 if _onsite.shape[0] != norbs:
453 msg = ('Single `onsite` matrix of shape ({0}, {0}) provided '
454 'but there are {1} orbitals per site in the system')
455 raise ValueError(msg.format(_onsite.shape[0], norbs))
456 else:
457 msg = ('Single `onsite` matrix provided, but there are '
458 'different numbers of orbitals on different sites')
459 raise ValueError(msg)
461 return _onsite, param_names
464def _make_onsite_or_hopping_terms(site_ranges, where):
466 terms = dict()
468 cdef gint[:] site_offsets_ = np.asarray(site_ranges, dtype=gint_dtype)[:, 0]
469 cdef gint i
471 if where.shape[1] == 1: # onsite
472 for i in range(where.shape[0]):
473 a = _bisect(site_offsets_, where[i, 0]) - 1
474 terms.setdefault((a, a), []).append(i)
475 else: # hopping
476 for i in range(where.shape[0]):
477 a = _bisect(site_offsets_, where[i, 0]) - 1
478 b = _bisect(site_offsets_, where[i, 1]) - 1
479 terms.setdefault((a, b), []).append(i)
480 return [(a, None, b) for a, b in terms.items()]
483def _vectorized_make_onsite_terms(syst, where):
484 assert is_vectorized(syst)
485 assert where.shape[1] == 1
486 site_offsets = [r[0] for r in syst.site_ranges]
488 terms = {}
489 term_by_id = syst._onsite_term_by_site_id
490 for i in range(where.shape[0]):
491 w = where[i, 0]
492 terms.setdefault(term_by_id[w], []).append(i)
494 ret = []
495 for term_id, which in terms.items():
496 term = syst.terms[term_id]
497 ((term_sa, _), (term_site_offsets, _)) = syst.subgraphs[term.subgraph]
498 term_sites = site_offsets[term_sa] + term_site_offsets
499 which = np.asarray(which, dtype=gint_dtype)
500 sites = _select(where, which).reshape(-1)
501 selector = np.searchsorted(term_sites, sites)
502 ret.append((term_id, selector, which))
504 return ret
507def _vectorized_make_hopping_terms(syst, where):
508 assert is_vectorized(syst)
509 assert where.shape[1] == 2
510 site_offsets = [r[0] for r in syst.site_ranges]
512 terms = {}
513 term_by_id = syst._hopping_term_by_edge_id
514 for i in range(where.shape[0]):
515 a, b = where[i, 0], where[i, 1]
516 edge = syst.graph.first_edge_id(a, b)
517 terms.setdefault(term_by_id[edge], []).append(i)
519 ret = []
520 dtype = np.dtype([('f0', int), ('f1', int)])
521 for term_id, which in terms.items():
522 herm_conj = term_id < 0
523 if herm_conj:
524 real_term_id = -term_id - 1
525 else:
526 real_term_id = term_id
527 which = np.asarray(which, dtype=gint_dtype)
528 # Select out the hoppings and reverse them if we are
529 # dealing with Hermitian conjugate hoppings
530 hops = _select(where, which)
531 if herm_conj:
532 hops = hops[:, ::-1]
533 # Force contiguous array to handle herm conj case also.
534 # Needs to be contiguous to cast to compound dtype
535 hops = np.ascontiguousarray(hops, dtype=int)
536 hops = hops.view(dtype).reshape(-1)
537 term = syst.terms[real_term_id]
538 # Get array of pairs
539 ((to_sa, from_sa), term_hoppings) = syst.subgraphs[term.subgraph]
540 term_hoppings = term_hoppings.transpose() + (site_offsets[to_sa], site_offsets[from_sa])
541 term_hoppings = term_hoppings.view(dtype).reshape(-1)
543 selector = np.recarray.searchsorted(term_hoppings, hops)
545 ret.append((term_id, selector, which))
547 return ret
550def _make_matrix_elements(evaluate_term, terms):
551 matrix_elements = []
552 for (term_id, term_selector, which) in terms:
553 which = np.asarray(which, dtype=gint_dtype)
554 data = evaluate_term(term_id, term_selector, which)
555 matrix_elements.append((which, data))
556 return matrix_elements
559cdef class BlockSparseMatrix:
560 """A sparse matrix stored as dense blocks.
562 Parameters
563 ----------
564 where : gint[:, :]
565 ``Nx2`` matrix or ``Nx1`` matrix: the arguments ``a``
566 and ``b`` to be used when evaluating ``f``. If an
567 ``Nx1`` matrix, then ``b=a``.
568 block_offsets : gint[:, :]
569 The row and column offsets for the start of each block
570 in the sparse matrix: ``(row_offset, col_offset)``.
571 block_shapes : gint[:, :]
572 ``Nx2`` array: the shapes of each block, ``(n_rows, n_cols)``.
573 matrix_elements : sequence of pairs (where_indices, data)
574 'data' is a 3D complex array; a vector of matrices.
575 'where_indices' is a 1D array of indices for 'where';
576 'data[i]' should be placed at 'where[where_indices[i]]'.
578 Attributes
579 ----------
580 block_offsets : gint[:, :]
581 The row and column offsets for the start of each block
582 in the sparse matrix: ``(row_offset, col_offset)``.
583 block_shapes : gint[:, :]
584 The shape of each block: ``(n_rows, n_cols)``
585 data_offsets : gint[:]
586 The offsets of the start of each matrix block in `data`.
587 data : complex[:]
588 The matrix of each block, stored in row-major (C) order.
589 """
591 @cython.embedsignature
592 @cython.boundscheck(False)
593 @cython.wraparound(False)
594 def __init__(self, gint[:, :] where, gint[:, :] block_offsets,
595 gint[:, :] block_shapes, matrix_elements):
596 if (block_offsets.shape[0] != where.shape[0] or
597 block_shapes.shape[0] != where.shape[0]):
598 raise ValueError('Arrays should be the same length along '
599 'the first axis.')
600 self.block_shapes = block_shapes
601 self.block_offsets = block_offsets
602 self.data_offsets = np.empty(where.shape[0], dtype=gint_dtype)
603 ### calculate shapes and data_offsets
604 cdef gint w, data_size = 0
605 for w in range(where.shape[0]):
606 self.data_offsets[w] = data_size
607 data_size += block_shapes[w, 0] * block_shapes[w, 1]
608 ### Populate data array
609 self.data = np.empty((data_size,), dtype=complex)
610 cdef complex[:, :, :] data
611 cdef gint[:] where_indexes
612 cdef gint i, j, k, off, a, b, a_norbs, b_norbs
613 for where_indexes, data in matrix_elements:
614 for i in range(where_indexes.shape[0]):
615 w = where_indexes[i]
616 off = self.data_offsets[w]
617 a_norbs = self.block_shapes[w, 0]
618 b_norbs = self.block_shapes[w, 1]
619 # Copy data
620 for j in range(a_norbs):
621 for k in range(b_norbs):
622 self.data[off + j * b_norbs + k] = data[i, j, k]
624 cdef complex* get(self, gint block_idx):
625 return <complex*> &self.data[0] + self.data_offsets[block_idx]
627 def __getstate__(self):
628 return tuple(map(np.asarray, (
629 self.block_offsets,
630 self.block_shapes,
631 self.data_offsets,
632 self.data
633 )))
635 def __setstate__(self, state):
636 (self.block_offsets,
637 self.block_shapes,
638 self.data_offsets,
639 self.data,
640 ) = state
643################ Local Observables
645# supported operations within the `_operate` method
646ctypedef enum operation:
647 MAT_ELS
648 ACT
651cdef class _LocalOperator:
652 """Base class for operators defined by an on-site matrix and the
653 Hamiltonian.
655 This includes "true" local operators, as well as "currents" and "sources".
657 Attributes
658 ----------
659 syst : `~kwant.system.System`
660 The system for which this operator is defined. Must have the
661 number of orbitals defined for all site families.
662 onsite : complex 2D array, or callable
663 If a complex array, then the same onsite is used everywhere.
664 Otherwise, function that can be called with a single site (integer) and
665 extra arguments, and returns the representation of the operator on
666 that site. This should return either a scalar or a square matrix of the
667 same shape as that returned by the system Hamiltonian evaluated on the
668 same site. The extra arguments must be the same as the extra arguments
669 to ``syst.hamiltonian``.
670 where : 2D array of `int` or `None`
671 where to evaluate the operator. A list of sites for on-site
672 operators (accessed like `where[n, 0]`), otherwise a list of pairs
673 of sites (accessed like `where[n, 0]` and `where[n, 1]`).
674 check_hermiticity : bool
675 If True, checks that ``onsite``, as well as any relevant parts
676 of the Hamiltonian are hermitian.
677 sum : bool, default: False
678 If True, then calling this operator will return a single scalar,
679 otherwise a vector will be returned (see
680 `~kwant.operator._LocalOperator.__call__` for details).
681 """
683 @cython.embedsignature
684 def __init__(self, syst, onsite, where, *,
685 check_hermiticity=True, sum=False):
686 _check_norbs(syst)
687 if is_vectorized(syst) and is_infinite(syst):
688 raise TypeError('Vectorized infinite systems cannot yet be '
689 'used with operators.')
691 self.syst = syst
692 self.onsite, self._onsite_param_names = _normalize_onsite(
693 syst, onsite, check_hermiticity)
694 self.check_hermiticity = check_hermiticity
695 self.sum = sum
696 self._site_ranges = np.asarray(syst.site_ranges, dtype=gint_dtype)
697 self.where = where
698 self._bound_onsite = None
699 self._bound_hamiltonian = None
701 # Here we pre-compute the datastructures that will enable us to evaluate
702 # the Hamiltonian and onsite functions in a vectorized way. Essentially
703 # we group the sites/hoppings in 'where' by what term of 'syst' they are
704 # in (for vectorized systems), or by the site family(s) (for
705 # non-vectorized systems). If the system is vectorized we store a list
706 # of triples:
707 #
708 # (term_id, term_selector, which)
709 #
710 # otherwise
711 #
712 # ((to_site_range, from_site_range), None, which)
713 #
714 # Where:
715 #
716 # 'term_id' → integer: term index, may be negative (indicates herm. conj.)
717 # 'term_selector' → 1D integer array: selects which elements from the
718 # subgraph of term number 'term_id' should be evaluated.
719 # 'which' → 1D integer array: selects which elements of 'where' this
720 # vectorized evaluation corresponds to.
721 # 'to/from_site_range' → integer: the site ranges that the elements of
722 # 'where' indexed by 'which' correspond to.
723 #
724 # Note that all sites/hoppings indicated by 'which' belong to the *same*
725 # pair of site families by construction. This is what allows for
726 # vectorized evaluation.
727 if is_vectorized(syst):
728 if self.where.shape[1] == 1:
729 self._terms = _vectorized_make_onsite_terms(syst, where)
730 else:
731 self._terms = _vectorized_make_hopping_terms(syst, where)
732 else:
733 self._terms = _make_onsite_or_hopping_terms(self._site_ranges, where)
735 @cython.embedsignature
736 def __call__(self, bra, ket=None, args=(), *, params=None):
737 r"""Return the matrix elements of the operator.
739 An operator ``A`` can be called like
741 >>> A(psi)
743 to compute the expectation value :math:`\bra{ψ} A \ket{ψ}`,
744 or like
746 >>> A(phi, psi)
748 to compute the matrix element :math:`\bra{φ} A \ket{ψ}`.
750 If ``sum=True`` was provided when constructing the operator, then
751 a scalar is returned. If ``sum=False``, then a vector is returned.
752 The vector is defined over the sites of the system if the operator
753 is a `~kwant.operator.Density`, or over the hoppings if it is a
754 `~kwant.operator.Current` or `~kwant.operator.Source`. By default,
755 the returned vector is ordered in the same way as the sites
756 (for `~kwant.operator.Density`) or hoppings in the graph of the
757 system (for `~kwant.operator.Current` or `~kwant.operator.Density`).
758 If the keyword parameter ``where`` was provided when constructing
759 the operator, then the returned vector is instead defined only over
760 the sites or hoppings specified, and is ordered in the same way
761 as ``where``.
763 Alternatively stated, for an operator :math:`Q_{iαβ}`, ``bra``
764 :math:`φ_α` and ``ket`` :math:`ψ_β` this computes
765 :math:`q_i = ∑_{αβ} φ^*_α Q_{iαβ} ψ_β` if ``self.sum`` is False,
766 otherwise computes :math:`q = ∑_{iαβ} φ^*_α Q_{iαβ} ψ_β`. where
767 :math:`i` runs over all sites or hoppings, and
768 :math:`α` and :math:`β` run over all the degrees of freedom.
770 Parameters
771 ----------
772 bra, ket : sequence of complex
773 Must have the same length as the number of orbitals
774 in the system. If only one is provided, both ``bra``
775 and ``ket`` are taken as equal.
776 args : tuple, optional
777 The arguments to pass to the system. Used to evaluate
778 the ``onsite`` elements and, possibly, the system Hamiltonian.
779 Deprecated in favor of 'params' (and mutually exclusive with it).
780 params : dict, optional
781 Dictionary of parameter names and their values. Mutually exclusive
782 with 'args'.
784 Returns
785 -------
786 `float` if ``check_hermiticity`` is True, and ``ket`` is ``None``,
787 otherwise `complex`. If this operator was created with ``sum=True``,
788 then a single value is returned, otherwise an array is returned.
789 """
790 if (self._bound_onsite or self._bound_hamiltonian) and (args or params):
791 raise ValueError("Extra arguments are already bound to this "
792 "operator. You should call this operator "
793 "providing neither 'args' nor 'params'.")
794 if args:
795 # deprecate_args does not play nicely with methods of cdef classes,
796 # when used as a decorator, so we manually raise the
797 # deprecation warning here.
798 deprecate_args()
799 if args and params:
800 raise TypeError("'args' and 'params' are mutually exclusive.")
801 if bra is None:
802 raise TypeError('bra must be an array')
803 bra = np.asarray(bra, dtype=complex)
804 ket = bra if ket is None else np.asarray(ket, dtype=complex)
805 tot_norbs = _get_tot_norbs(self.syst)
806 if bra.shape != (tot_norbs,):
807 msg = 'vector is incorrect shape'
808 msg = 'bra ' + msg if ket is None else msg
809 raise ValueError(msg)
810 elif ket.shape != (tot_norbs,):
811 raise ValueError('ket vector is incorrect shape')
813 result = np.zeros((self.where.shape[0],), dtype=complex)
814 self._operate(out_data=result, bra=bra, ket=ket, args=args,
815 params=params, op=MAT_ELS)
816 # if everything is Hermitian then result is real if bra == ket
817 if self.check_hermiticity and bra is ket:
818 result = result.real
819 return np.sum(result) if self.sum else result
821 @cython.embedsignature
822 def act(self, ket, args=(), *, params=None):
823 """Act with the operator on a wavefunction.
825 For an operator :math:`Q_{iαβ}` and ``ket`` :math:`ψ_β`
826 this computes :math:`∑_{iβ} Q_{iαβ} ψ_β`.
828 Parameters
829 ----------
830 ket : sequence of complex
831 Wavefunctions defined over all the orbitals of the system.
832 args : tuple
833 The extra arguments to the Hamiltonian value functions and
834 the operator ``onsite`` function.
835 Deprecated in favor of 'params' (and mutually exclusive with it).
836 params : dict, optional
837 Dictionary of parameter names and their values. Mutually exclusive
838 with 'args'.
840 Returns
841 -------
842 Array of `complex`.
843 """
844 if (self._bound_onsite or self._bound_hamiltonian) and (args or params):
845 raise ValueError("Extra arguments are already bound to this "
846 "operator. You should call this operator "
847 "providing neither 'args' nor 'params'.")
848 if args:
849 # deprecate_args does not play nicely with methods of cdef classes,
850 # when used as a decorator, so we manually raise the
851 # deprecation warning here.
852 deprecate_args()
853 if args and params:
854 raise TypeError("'args' and 'params' are mutually exclusive.")
856 if ket is None:
857 raise TypeError('ket must be an array')
858 ket = np.asarray(ket, dtype=complex)
859 tot_norbs = _get_tot_norbs(self.syst)
860 if ket.shape != (tot_norbs,):
861 raise ValueError('ket vector is incorrect shape')
862 result = np.zeros((tot_norbs,), dtype=complex)
863 self._operate(out_data=result, bra=None, ket=ket, args=args,
864 params=params, op=ACT)
865 return result
867 @cython.embedsignature
868 def bind(self, args=(), *, params=None):
869 """Bind the given arguments to this operator.
871 Returns a copy of this operator that does not need to be passed extra
872 arguments when subsequently called or when using the ``act`` method.
874 Providing positional arguments via 'args' is deprecated,
875 instead provide named parameters as a dictionary via 'params'.
876 """
877 if args:
878 # deprecate_args does not play nicely with methods of cdef classes,
879 # when used as a decorator, so we manually raise the
880 # deprecation warning here.
881 deprecate_args()
882 if args and params:
883 raise TypeError("'args' and 'params' are mutually exclusive.")
884 # generic creation of new instance
885 cls = self.__class__
886 q = cls.__new__(cls)
887 q.syst = self.syst
888 q.onsite = self.onsite
889 q._onsite_param_names = self._onsite_param_names
890 q.where = self.where
891 q.sum = self.sum
892 q._site_ranges = self._site_ranges
893 q.check_hermiticity = self.check_hermiticity
894 if callable(self.onsite):
895 q._bound_onsite = self._eval_onsites(args, params)
896 # NOTE: subclasses should populate `bound_hamiltonian` if needed
897 return q
899 def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
900 args, operation op, *, params=None):
901 """Do an operation with the operator.
903 Parameters
904 ----------
905 out_data : ndarray
906 Output array, zero on entry. On exit should contain the required
907 data. What this means depends on the value of `op`, as does the
908 length of the array.
909 bra, ket : ndarray
910 Wavefunctions defined over all the orbitals of the system.
911 If `op` is `ACT` then `bra` is None.
912 args : tuple
913 The extra arguments to the Hamiltonian value functions and
914 the operator ``onsite`` function.
915 Deprecated in favor of 'params' (and mutually exclusive with it).
916 op : operation
917 The operation to perform.
918 `MAT_ELS`: calculate matrix elements between `bra` and `ket`
919 `ACT`: act on `ket` with the operator
920 params : dict, optional
921 Dictionary of parameter names and their values. Mutually exclusive
922 with 'args'.
923 """
924 raise NotImplementedError()
926 cdef BlockSparseMatrix _eval_onsites(self, args, params):
927 """Evaluate the onsite matrices on all elements of `where`"""
928 assert callable(self.onsite)
929 assert not (args and params)
930 check_hermiticity = self.check_hermiticity
931 syst = self.syst
933 _is_vectorized = is_vectorized(syst)
935 if params:
936 try:
937 args = tuple(params[pn] for pn in self._onsite_param_names)
938 except KeyError:
939 missing = [p for p in self._onsite_param_names
940 if p not in params]
941 msg = ('Operator is missing required arguments: ',
942 ', '.join(map('"{}"'.format, missing)))
943 raise TypeError(''.join(msg))
945 # Evaluate many onsites at once. See _LocalOperator.__init__
946 # for an explanation the parameters.
947 def eval_onsite(term_id, term_selector, which):
948 if _is_vectorized:
949 if term_id >= 0:
950 (sr, _), _ = syst.subgraphs[syst.terms[term_id].subgraph]
951 else:
952 (_, sr), _ = syst.subgraphs[syst.terms[-term_id - 1].subgraph]
953 else:
954 sr, _ = term_id
955 start_site, norbs, _ = self.syst.site_ranges[sr]
956 # All sites selected by 'which' are part of the same site family.
957 site_offsets = _select(self.where, which)[:, 0] - start_site
958 data = self.onsite(sr, site_offsets, *args)
959 _check_onsites(data, norbs, self.check_hermiticity)
960 return data
962 matrix_elements = _make_matrix_elements(eval_onsite, self._terms)
963 offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
964 return BlockSparseMatrix(self.where, offsets, norbs, matrix_elements)
967 cdef BlockSparseMatrix _eval_hamiltonian(self, args, params):
968 """Evaluate the Hamiltonian on all elements of `where`."""
970 where = self.where
971 syst = self.syst
972 is_onsite = self.where.shape[1] == 1
973 check_hermiticity = self.check_hermiticity
975 if is_vectorized(self.syst):
977 # Evaluate many Hamiltonian elements at once.
978 # See _LocalOperator.__init__ for an explanation the parameters.
979 def eval_hamiltonian(term_id, term_selector, which):
980 herm_conj = term_id < 0
981 assert not is_onsite or (is_onsite and not herm_conj) # onsite terms are never hermitian conjugate
982 if herm_conj:
983 term_id = -term_id - 1
984 data = syst.hamiltonian_term(term_id, term_selector,
985 args=args, params=params)
986 if herm_conj:
987 data = data.conjugate().transpose(0, 2, 1)
988 # Checks for data consistency
989 (to_sr, from_sr), _ = syst.subgraphs[syst.terms[term_id].subgraph]
990 to_norbs = syst.site_ranges[to_sr][1]
991 from_norbs = syst.site_ranges[from_sr][1]
992 if herm_conj:
993 to_norbs, from_norbs = from_norbs, to_norbs
994 _check_hams(data, to_norbs, from_norbs, is_onsite and check_hermiticity)
996 return data
998 else:
1000 # Evaluate many Hamiltonian elements at once.
1001 # See _LocalOperator.__init__ for an explanation the parameters.
1002 def eval_hamiltonian(term_id, term_selector, which):
1003 if is_onsite:
1004 data = [
1005 syst.hamiltonian(where[i, 0], where[i, 0], *args, params=params)
1006 for i in which
1007 ]
1008 else:
1009 data = [
1010 syst.hamiltonian(where[i, 0], where[i, 1], *args, params=params)
1011 for i in which
1012 ]
1014 (to_sr, from_sr) = term_id
1015 to_norbs, from_norbs = syst.site_ranges[to_sr][1], syst.site_ranges[from_sr][1]
1016 expected_shape = (len(which), to_norbs, from_norbs)
1017 data = _normalize_matrix_blocks(data, expected_shape)
1018 # Checks for data consistency
1020 _check_hams(data, to_norbs, from_norbs, is_onsite and check_hermiticity)
1022 return data
1024 matrix_elements = _make_matrix_elements(eval_hamiltonian, self._terms)
1025 offsets, norbs = _get_all_orbs(where, self._site_ranges)
1026 return BlockSparseMatrix(where, offsets, norbs, matrix_elements)
1028 def __getstate__(self):
1029 return (
1030 (self.check_hermiticity, self.sum),
1031 (self.syst, self.onsite, self._onsite_param_names),
1032 tuple(map(np.asarray, (self.where, self._site_ranges))),
1033 (self._terms,),
1034 (self._bound_onsite, self._bound_hamiltonian),
1035 )
1037 def __setstate__(self, state):
1038 ((self.check_hermiticity, self.sum),
1039 (self.syst, self.onsite, self._onsite_param_names),
1040 (self.where, self._site_ranges),
1041 (self._terms,),
1042 (self._bound_onsite, self._bound_hamiltonian),
1043 ) = state
1046cdef class Density(_LocalOperator):
1047 """An operator for calculating general densities.
1049 An instance of this class can be called like a function to evaluate the
1050 expectation value with a wavefunction. See
1051 `~kwant.operator.Density.__call__` for details.
1053 Parameters
1054 ----------
1055 syst : `~kwant.system.System`
1056 onsite : scalar or square matrix or dict or callable
1057 The onsite matrix that defines the operator. If a dict is given, it
1058 maps from site families to square matrices. If a function is given it
1059 must take the same arguments as the onsite Hamiltonian functions of the
1060 system.
1061 where : sequence of `int` or `~kwant.system.Site`, or callable, optional
1062 Where to evaluate the operator. If ``syst`` is not a finalized Builder,
1063 then this should be a sequence of integers. If a function is provided,
1064 it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
1065 a finalized builder) and return True or False. If not provided, the
1066 operator will be calculated over all sites in the system.
1067 check_hermiticity: bool
1068 Check whether the provided ``onsite`` is Hermitian. If it is not
1069 Hermitian, then an error will be raised when the operator is
1070 evaluated.
1071 sum : bool, default: False
1072 If True, then calling this operator will return a single scalar,
1073 otherwise a vector will be returned (see
1074 `~kwant.operator.Density.__call__` for details).
1076 Notes
1077 -----
1078 In general, if there is a certain "density" (e.g. charge or spin) that is
1079 represented by a square matrix :math:`M_i` associated with each site
1080 :math:`i` then an instance of this class represents the tensor
1081 :math:`Q_{iαβ}` which is equal to :math:`M_i` when α and β are orbitals on
1082 site :math:`i`, and zero otherwise.
1083 """
1085 @cython.embedsignature
1086 def __init__(self, syst, onsite=1, where=None, *,
1087 check_hermiticity=True, sum=False):
1088 where = _normalize_site_where(syst, where)
1089 super().__init__(syst, onsite, where,
1090 check_hermiticity=check_hermiticity, sum=sum)
1092 @cython.boundscheck(False)
1093 @cython.wraparound(False)
1094 def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
1095 args, operation op, *, params=None):
1096 matrix = ta.matrix
1097 cdef int unique_onsite = not callable(self.onsite)
1098 # prepare onsite matrices
1099 cdef complex[:, :] _tmp_mat
1100 cdef complex *M_a = NULL
1101 cdef BlockSparseMatrix M_a_blocks
1103 if unique_onsite:
1104 _tmp_mat = self.onsite
1105 M_a = <complex*> &_tmp_mat[0, 0]
1106 elif self._bound_onsite:
1107 M_a_blocks = self._bound_onsite
1108 else:
1109 M_a_blocks = self._eval_onsites(args, params)
1111 # loop-local variables
1112 cdef gint a, a_s, a_norbs
1113 cdef gint i, j, w
1114 cdef complex tmp, bra_conj
1115 ### loop over sites
1116 for w in range(self.where.shape[0]):
1117 ### get the next site, start orbital and number of orbitals
1118 a = self.where[w, 0]
1119 _get_orbs(self._site_ranges, a, &a_s, &a_norbs)
1120 ### get the next onsite matrix, if necessary
1121 if not unique_onsite:
1122 M_a = M_a_blocks.get(w)
1123 ### do the actual calculation
1124 if op == MAT_ELS:
1125 tmp = 0
1126 for i in range(a_norbs):
1127 for j in range(a_norbs):
1128 tmp += (bra[a_s + i].conjugate() *
1129 M_a[i * a_norbs + j] * ket[a_s + j])
1130 out_data[w] = tmp
1131 elif op == ACT:
1132 for i in range(a_norbs):
1133 tmp = 0
1134 for j in range(a_norbs):
1135 tmp += M_a[i * a_norbs + j] * ket[a_s + j]
1136 out_data[a_s + i] = out_data[a_s + i] + tmp
1138 @cython.boundscheck(False)
1139 @cython.wraparound(False)
1140 @cython.cdivision(True)
1141 @cython.embedsignature
1142 def tocoo(self, args=(), *, params=None):
1143 """Convert the operator to coordinate format sparse matrix.
1145 Providing positional arguments via 'args' is deprecated,
1146 instead provide named parameters as a dictionary via 'params'.
1147 """
1148 cdef int blk, blk_size, n_blocks, n, k = 0
1149 cdef int [:, :] offsets, shapes
1150 cdef int [:] row, col
1151 if self._bound_onsite and (args or params):
1152 raise ValueError("Extra arguments are already bound to this "
1153 "operator. You should call this operator "
1154 "providing neither 'args' nor 'params'.")
1155 if args:
1156 # deprecate_args does not play nicely with methods of cdef classes,
1157 # when used as a decorator, so we manually raise the
1158 # deprecation warning here.
1159 deprecate_args()
1160 if args and params:
1161 raise TypeError("'args' and 'params' are mutually exclusive.")
1163 if not callable(self.onsite):
1164 offsets = _get_all_orbs(self.where, self._site_ranges)[0]
1165 n_blocks = len(self.where)
1166 shapes = np.asarray(np.resize([self.onsite.shape], (n_blocks, 2)),
1167 gint_dtype)
1168 data = np.asarray(self.onsite).flatten()
1169 data = np.resize(data, [len(data) * n_blocks])
1170 else:
1171 if self._bound_onsite is not None:
1172 onsite_matrix = self._bound_onsite
1173 else:
1174 onsite_matrix = self._eval_onsites(args, params)
1175 data = onsite_matrix.data
1176 offsets = np.asarray(onsite_matrix.block_offsets)
1177 shapes = np.asarray(onsite_matrix.block_shapes)
1179 row = np.empty(len(data), gint_dtype)
1180 col = np.empty(len(data), gint_dtype)
1181 for blk in range(len(offsets)):
1182 blk_size = shapes[blk, 0] * shapes[blk, 1]
1183 for n in range(blk_size):
1184 row[k] = offsets[blk, 0] + n // shapes[blk, 1]
1185 col[k] = offsets[blk, 1] + n % shapes[blk, 1]
1186 k += 1
1188 norbs = _get_tot_norbs(self.syst)
1189 return coo_matrix((np.asarray(data),
1190 (np.asarray(row), np.asarray(col))),
1191 shape=(norbs, norbs))
1194cdef class Current(_LocalOperator):
1195 r"""An operator for calculating general currents.
1197 An instance of this class can be called like a function to evaluate the
1198 expectation value with a wavefunction. See
1199 `~kwant.operator.Current.__call__` for details.
1201 Parameters
1202 ----------
1203 syst : `~kwant.system.System`
1204 onsite : scalar or square matrix or dict or callable
1205 The onsite matrix that defines the density from which this current is
1206 derived. If a dict is given, it maps from site families to square
1207 matrices (scalars are allowed if the site family has 1 orbital per
1208 site). If a function is given it must take the same arguments as the
1209 onsite Hamiltonian functions of the system.
1210 where : sequence of pairs of `int` or `~kwant.system.Site`, or callable, optional
1211 Where to evaluate the operator. If ``syst`` is not a finalized Builder,
1212 then this should be a sequence of pairs of integers. If a function is
1213 provided, it should take a pair of integers or a pair of
1214 `~kwant.system.Site` (if ``syst`` is a finalized builder) and return
1215 True or False. If not provided, the operator will be calculated over
1216 all hoppings in the system.
1217 check_hermiticity : bool
1218 Check whether the provided ``onsite`` is Hermitian. If it
1219 is not Hermitian, then an error will be raised when the
1220 operator is evaluated.
1221 sum : bool, default: False
1222 If True, then calling this operator will return a single scalar,
1223 otherwise a vector will be returned (see
1224 `~kwant.operator.Current.__call__` for details).
1226 Notes
1227 -----
1228 In general, if there is a certain "density" (e.g. charge or spin) that is
1229 represented by a square matrix :math:`M_i` associated with each site
1230 :math:`i` and :math:`H_{ij}` is the hopping Hamiltonian from site :math:`j`
1231 to site `i`, then an instance of this class represents the tensor
1232 :math:`J_{ijαβ}` which is equal to :math:`i\left[(H_{ij})^† M_i - M_i
1233 H_{ij}\right]` when α and β are orbitals on sites :math:`i` and :math:`j`
1234 respectively, and zero otherwise.
1236 The tensor :math:`J_{ijαβ}` will also be referred to as :math:`Q_{nαβ}`,
1237 where :math:`n` is the index of hopping :math:`(i, j)` in ``where``.
1238 """
1240 @cython.embedsignature
1241 def __init__(self, syst, onsite=1, where=None, *,
1242 check_hermiticity=True, sum=False):
1243 where = _normalize_hopping_where(syst, where)
1244 super().__init__(syst, onsite, where,
1245 check_hermiticity=check_hermiticity, sum=sum)
1247 @cython.embedsignature
1248 def bind(self, args=(), *, params=None):
1249 """Bind the given arguments to this operator.
1251 Returns a copy of this operator that does not need to be passed extra
1252 arguments when subsequently called or when using the ``act`` method.
1254 Providing positional arguments via 'args' is deprecated,
1255 instead provide named parameters as a dictionary via 'params'.
1256 """
1257 q = super().bind(args, params=params)
1258 q._bound_hamiltonian = self._eval_hamiltonian(args, params)
1259 return q
1261 @cython.boundscheck(False)
1262 @cython.wraparound(False)
1263 def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
1264 args, operation op, *, params=None):
1265 # prepare onsite matrices and hamiltonians
1266 cdef int unique_onsite = not callable(self.onsite)
1267 cdef complex[:, :] _tmp_mat
1268 cdef complex *M_a = NULL
1269 cdef complex *H_ab = NULL
1270 cdef BlockSparseMatrix M_a_blocks, H_ab_blocks
1272 if unique_onsite:
1273 _tmp_mat = self.onsite
1274 M_a = <complex*> &_tmp_mat[0, 0]
1275 elif self._bound_onsite:
1276 M_a_blocks = self._bound_onsite
1277 else:
1278 M_a_blocks = self._eval_onsites(args, params)
1280 if self._bound_hamiltonian:
1281 H_ab_blocks = self._bound_hamiltonian
1282 else:
1283 H_ab_blocks = self._eval_hamiltonian(args, params)
1285 # main loop
1286 cdef gint a, a_s, a_norbs, b, b_s, b_norbs
1287 cdef gint i, j, k, w
1288 cdef complex tmp
1289 for w in range(self.where.shape[0]):
1290 ### get the next hopping's start orbitals and numbers of orbitals
1291 a_s = H_ab_blocks.block_offsets[w, 0]
1292 b_s = H_ab_blocks.block_offsets[w, 1]
1293 a_norbs = H_ab_blocks.block_shapes[w, 0]
1294 b_norbs = H_ab_blocks.block_shapes[w, 1]
1295 ### get the next onsite and Hamiltonian matrices
1296 H_ab = H_ab_blocks.get(w)
1297 if not unique_onsite:
1298 M_a = M_a_blocks.get(w)
1299 ### do the actual calculation
1300 if op == MAT_ELS:
1301 tmp = 0
1302 for i in range(b_norbs):
1303 for j in range(a_norbs):
1304 for k in range(a_norbs):
1305 tmp += (bra[b_s + i].conjugate() *
1306 H_ab[j * b_norbs + i].conjugate() *
1307 M_a[j * a_norbs + k] * ket[a_s + k]
1308 - bra[a_s + j].conjugate() *
1309 M_a[j * a_norbs + k] *
1310 H_ab[k * b_norbs + i] * ket[b_s + i])
1311 out_data[w] = 1j * tmp
1312 elif op == ACT:
1313 for i in range(b_norbs):
1314 for j in range(a_norbs):
1315 for k in range(a_norbs):
1316 out_data[b_s + i] = (
1317 out_data[b_s + i] +
1318 1j * H_ab[j * b_norbs + i].conjugate() *
1319 M_a[j * a_norbs + k] * ket[a_s + k])
1320 out_data[a_s + j] = (
1321 out_data[a_s + j] -
1322 1j * M_a[j * a_norbs + k] * H_ab[k * b_norbs + i] *
1323 ket[b_s + i])
1326cdef class Source(_LocalOperator):
1327 """An operator for calculating general sources.
1329 An instance of this class can be called like a function to evaluate the
1330 expectation value with a wavefunction. See
1331 `~kwant.operator.Source.__call__` for details.
1333 Parameters
1334 ----------
1335 syst : `~kwant.system.System`
1336 onsite : scalar or square matrix or dict or callable
1337 The onsite matrix that defines the density from which this source is
1338 defined. If a dict is given, it maps from site families to square
1339 matrices (scalars are allowed if the site family has 1 orbital per
1340 site). If a function is given it must take the same arguments as the
1341 onsite Hamiltonian functions of the system.
1342 where : sequence of `int` or `~kwant.system.Site`, or callable, optional
1343 Where to evaluate the operator. If ``syst`` is not a finalized Builder,
1344 then this should be a sequence of integers. If a function is provided,
1345 it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
1346 a finalized builder) and return True or False. If not provided, the
1347 operator will be calculated over all sites in the system.
1348 check_hermiticity : bool
1349 Check whether the provided ``onsite`` is Hermitian. If it is not
1350 Hermitian, then an error will be raised when the operator is
1351 evaluated.
1352 sum : bool, default: False
1353 If True, then calling this operator will return a single scalar,
1354 otherwise a vector will be returned (see
1355 `~kwant.operator.Source.__call__` for details).
1357 Notes
1358 -----
1359 An example of a "source" is a spin torque. In general, if there is a
1360 certain "density" (e.g. charge or spin) that is represented by a square
1361 matrix :math:`M_i` associated with each site :math:`i`, and :math:`H_{i}`
1362 is the onsite Hamiltonian on site site `i`, then an instance of this class
1363 represents the tensor :math:`Q_{iαβ}` which is equal to :math:`(H_{i})^†
1364 M_i - M_i H_{i}` when α and β are orbitals on site :math:`i`, and zero
1365 otherwise.
1366 """
1368 @cython.embedsignature
1369 def __init__(self, syst, onsite=1, where=None, *,
1370 check_hermiticity=True, sum=False):
1371 where = _normalize_site_where(syst, where)
1372 super().__init__(syst, onsite, where,
1373 check_hermiticity=check_hermiticity, sum=sum)
1375 @cython.embedsignature
1376 def bind(self, args=(), *, params=None):
1377 """Bind the given arguments to this operator.
1379 Returns a copy of this operator that does not need to be passed extra
1380 arguments when subsequently called or when using the ``act`` method.
1382 Providing positional arguments via 'args' is deprecated,
1383 instead provide named parameters as a dictionary via 'params'.
1384 """
1385 q = super().bind(args, params=params)
1386 q._bound_hamiltonian = self._eval_hamiltonian(args, params)
1387 return q
1389 @cython.boundscheck(False)
1390 @cython.wraparound(False)
1391 def _operate(self, complex[:] out_data, complex[:] bra, complex[:] ket,
1392 args, operation op, *, params=None):
1393 # prepare onsite matrices and hamiltonians
1394 cdef int unique_onsite = not callable(self.onsite)
1395 cdef complex[:, :] _tmp_mat
1396 cdef complex *M_a = NULL
1397 cdef complex *H_aa = NULL
1398 cdef BlockSparseMatrix M_a_blocks, H_aa_blocks
1400 if unique_onsite:
1401 _tmp_mat = self.onsite
1402 M_a = <complex*> &_tmp_mat[0, 0]
1403 elif self._bound_onsite:
1404 M_a_blocks = self._bound_onsite
1405 else:
1406 M_a_blocks = self._eval_onsites(args, params)
1408 if self._bound_hamiltonian:
1409 H_aa_blocks = self._bound_hamiltonian
1410 else:
1411 H_aa_blocks = self._eval_hamiltonian(args, params)
1413 # main loop
1414 cdef gint a, a_s, a_norbs
1415 cdef gint i, j, k, w
1416 cdef complex tmp, tmp2
1417 for w in range(self.where.shape[0]):
1418 ### get the next site, start orbital and number of orbitals
1419 # row offsets and block size are the same as for columns, as
1420 # we are only dealing with the block-diagonal part of H
1421 a_s = H_aa_blocks.block_offsets[w, 0]
1422 a_norbs = H_aa_blocks.block_shapes[w, 0]
1423 ### get the next onsite and Hamiltonian matrices
1424 H_aa = H_aa_blocks.get(w)
1425 if not unique_onsite:
1426 M_a = M_a_blocks.get(w)
1427 ### do the actual calculation
1428 if op == MAT_ELS:
1429 tmp2 = 0
1430 for i in range(a_norbs):
1431 tmp = 0
1432 for j in range(a_norbs):
1433 for k in range(a_norbs):
1434 tmp += (H_aa[j * a_norbs + i].conjugate() *
1435 M_a[j * a_norbs + k] * ket[a_s + k]
1436 - M_a[i * a_norbs + j] *
1437 H_aa[j * a_norbs + k] * ket[a_s + k])
1438 tmp2 += bra[a_s + i].conjugate() * tmp
1439 out_data[w] = 1j * tmp2
1440 elif op == ACT:
1441 for i in range(a_norbs):
1442 tmp = 0
1443 for j in range(a_norbs):
1444 for k in range(a_norbs):
1445 tmp += (H_aa[j * a_norbs + i].conjugate() *
1446 M_a[j * a_norbs + k] * ket[a_s + k]
1447 - M_a[i * a_norbs + j] *
1448 H_aa[j * a_norbs + k] * ket[a_s + k])
1449 out_data[a_s + i] = out_data[a_s + i] + 1j * tmp