Coverage for kwant/lattice.py : 84%
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.
9__all__ = ['TranslationalSymmetry', 'general', 'Polyatomic', 'Monatomic',
10 'chain', 'square', 'cubic', 'triangular', 'honeycomb', 'kagome']
12from math import sqrt
13from itertools import product
14import numpy as np
15import tinyarray as ta
16from . import builder, system
17from .linalg import lll
18from ._common import ensure_isinstance
21def general(prim_vecs, basis=None, name='', norbs=None):
22 """
23 Create a Bravais lattice of any dimensionality, with any number of sites.
25 Parameters
26 ----------
27 prim_vecs : 2d array-like of floats
28 The primitive vectors of the Bravais lattice
29 basis : 2d array-like of floats
30 The coordinates of the basis sites inside the unit cell
31 name : string or sequence of strings
32 Name of the lattice, or sequence of names of all of the sublattices.
33 If the name of the lattice is given, the names of sublattices (if any)
34 are obtained by appending their number to the name of the lattice.
35 norbs : int or sequence of ints, optional
36 The number of orbitals per site on the lattice, or a sequence
37 of the number of orbitals of sites on each of the sublattices.
39 Returns
40 -------
41 lattice : either `Monatomic` or `Polyatomic`
42 Resulting lattice.
44 Notes
45 -----
46 This function is largely an alias to the constructors of corresponding
47 lattices.
48 """
49 if basis is None:
50 return Monatomic(prim_vecs, name=name, norbs=norbs)
51 else:
52 return Polyatomic(prim_vecs, basis, name=name, norbs=norbs)
55def _check_prim_vecs(prim_vecs):
56 """Check constraints to ensure that prim_vecs is correct."""
57 if prim_vecs.ndim != 2:
58 raise ValueError('``prim_vecs`` must be a 2d array-like object.')
60 if prim_vecs.shape[0] > prim_vecs.shape[1]:
61 raise ValueError('Number of primitive vectors exceeds '
62 'the space dimensionality.')
64 if np.linalg.matrix_rank(prim_vecs) < len(prim_vecs):
65 raise ValueError('"prim_vecs" must be linearly independent.')
68class Polyatomic:
69 """
70 A Bravais lattice with an arbitrary number of sites in the basis.
72 Contains `Monatomic` sublattices. Note that an instance of ``Polyatomic`` is
73 not itself a `~kwant.system.SiteFamily`, only its sublattices are.
75 Parameters
76 ----------
77 prim_vecs : 2d array-like of floats
78 The primitive vectors of the Bravais lattice
79 basis : 2d array-like of floats
80 The coordinates of the basis sites inside the unit cell.
81 name : string or sequence of strings, optional
82 The name of the lattice, or a sequence of the names of all the
83 sublattices. If the name of the lattice is given, the names of
84 sublattices are obtained by appending their number to the name of the
85 lattice.
86 norbs : int or sequence of ints, optional
87 The number of orbitals per site on the lattice, or a sequence
88 of the number of orbitals of sites on each of the sublattices.
90 Raises
91 ------
92 ValueError
93 If dimensionalities do not match.
94 """
95 def __init__(self, prim_vecs, basis, name='', norbs=None):
96 prim_vecs = ta.array(prim_vecs, float)
97 _check_prim_vecs(prim_vecs)
99 dim = prim_vecs.shape[1]
100 if name is None: 100 ↛ 101line 100 didn't jump to line 101, because the condition on line 100 was never true
101 name = ''
102 if isinstance(name, str): 102 ↛ 105line 102 didn't jump to line 105, because the condition on line 102 was never false
103 name = [name + str(i) for i in range(len(basis))]
105 basis = ta.array(basis, float)
106 if basis.ndim != 2: 106 ↛ 107line 106 didn't jump to line 107, because the condition on line 106 was never true
107 raise ValueError('`basis` must be a 2d array-like object.')
108 if basis.shape[1] != dim: 108 ↛ 109line 108 didn't jump to line 109, because the condition on line 108 was never true
109 raise ValueError('Basis dimensionality does not match '
110 'the space dimensionality.')
112 try:
113 norbs = list(norbs)
114 if len(norbs) != len(basis):
115 raise ValueError('Length of `norbs` is not the same as '
116 'the number of basis vectors')
117 except TypeError:
118 norbs = [norbs] * len(basis)
120 self.sublattices = [Monatomic(prim_vecs, offset, sname, norb)
121 for offset, sname, norb in zip(basis, name, norbs)]
122 # Sequence of primitive vectors of the lattice.
123 self._prim_vecs = prim_vecs
124 # Precalculation of auxiliary arrays for real space calculations.
125 self.reduced_vecs, self.transf = lll.lll(prim_vecs)
126 self.voronoi = ta.dot(lll.voronoi(self.reduced_vecs), self.transf)
128 def __str__(self):
129 sl_names = ', '.join(str(sl.name) for sl in self.sublattices)
130 return '<Polyatomic lattice with sublattices {0}>'.format(sl_names)
132 def shape(self, function, start):
133 """Return a key for all the lattice sites inside a given shape.
135 The object returned by this method is primarily meant to be used as a
136 key for indexing `~kwant.builder.Builder` instances. See example below.
138 Parameters
139 ----------
140 function : callable
141 A function of real space coordinates that returns a truth value:
142 true for coordinates inside the shape, and false otherwise.
143 start : 1d array-like
144 The real-space origin for the flood-fill algorithm.
146 Returns
147 -------
148 shape_sites : function
150 Notes
151 -----
152 When the function returned by this method is called, a flood-fill
153 algorithm finds and yields all the lattice sites inside the specified
154 shape starting from the specified position.
156 A `~kwant.system.Symmetry` or `~kwant.builder.Builder` may be passed as
157 sole argument when calling the function returned by this method. This
158 will restrict the flood-fill to the fundamental domain of the symmetry
159 (or the builder's symmetry). Note that unless the shape function has
160 that symmetry itself, the result may be unexpected.
162 Examples
163 --------
164 >>> def circle(pos):
165 ... x, y = pos
166 ... return x**2 + y**2 < 100
167 ...
168 >>> lat = kwant.lattice.honeycomb()
169 >>> syst = kwant.Builder()
170 >>> syst[lat.shape(circle, (0, 0))] = 0
171 >>> syst[lat.neighbors()] = 1
172 """
173 def shape_sites(symmetry=None):
174 Site = system.Site
176 if symmetry is None:
177 symmetry = system.NoSymmetry()
178 elif not isinstance(symmetry, system.Symmetry):
179 symmetry = symmetry.symmetry
181 def fd_site(lat, tag):
182 return symmetry.to_fd(Site(lat, tag, True))
184 dim = len(start)
185 if dim != self._prim_vecs.shape[1]: 185 ↛ 186line 185 didn't jump to line 186, because the condition on line 185 was never true
186 raise ValueError('Dimensionality of start position does not '
187 'match the space dimensionality.')
188 lats = self.sublattices
189 deltas = list(self.voronoi)
191 #### Flood-fill ####
192 sites = []
193 for tag in set(lat.closest(start) for lat in lats):
194 for lat in lats:
195 site = fd_site(lat, tag)
196 if function(site.pos):
197 sites.append(site)
198 if not sites:
199 msg = 'No sites close to {0} are inside the desired shape.'
200 raise ValueError(msg.format(start))
201 tags = set(s.tag for s in sites)
203 while sites:
204 old_tags = tags
205 tags = set()
206 for site in sites:
207 yield site
208 tags.add(site.tag)
209 old_tags |= tags
211 new_tags = set()
212 for tag in tags:
213 for delta in deltas:
214 new_tag = tag + delta
215 if new_tag not in old_tags:
216 new_tags.add(new_tag)
218 sites = set()
219 for tag in new_tags:
220 for lat in lats:
221 site = fd_site(lat, tag)
222 if site.tag not in old_tags and function(site.pos):
223 sites.add(site)
225 return shape_sites
227 def wire(self, center, radius):
228 """Return a key for all the lattice sites inside an infinite cylinder.
230 This method makes it easy to define cylindrical (2d: rectangular) leads
231 that point in any direction. The object returned by this method is
232 primarily meant to be used as a key for indexing `~kwant.builder.Builder`
233 instances. See example below.
235 Parameters
236 ----------
237 center : 1d array-like of floats
238 A point belonging to the axis of the cylinder.
239 radius : float
240 The radius of the cylinder.
242 Notes
243 -----
244 The function returned by this method is to be called with a
245 `~kwant.builder.TranslationalSymmetry` instance (or a
246 `~kwant.builder.Builder` instance whose symmetry is used then) as sole
247 argument. All the lattice sites (in the fundamental domain of the
248 symmetry) inside the specified infinite cylinder are yielded. The
249 direction of the cylinder is determined by the symmetry.
251 Examples
252 --------
253 >>> lat = kwant.lattice.honeycomb()
254 >>> sym = kwant.TranslationalSymmetry(lat.a.vec((-2, 1)))
255 >>> lead = kwant.Builder(sym)
256 >>> lead[lat.wire((0, -5), 5)] = 0
257 >>> lead[lat.neighbors()] = 1
258 """
259 center = ta.array(center, float)
261 def wire_sites(sym):
262 if not isinstance(sym, system.Symmetry): 262 ↛ 264line 262 didn't jump to line 264, because the condition on line 262 was never false
263 sym = sym.symmetry
264 if not isinstance(sym, TranslationalSymmetry): 264 ↛ 265line 264 didn't jump to line 265, because the condition on line 264 was never true
265 raise ValueError('wire shape only works with '
266 'translational symmetry.')
267 if len(sym.periods) != 1: 267 ↛ 268line 267 didn't jump to line 268, because the condition on line 267 was never true
268 raise ValueError('wire shape only works with one-dimensional '
269 'translational symmetry.')
270 period = np.array(sym.periods)[0]
271 direction = ta.array(period / np.linalg.norm(period))
272 r_squared = radius**2
274 def wire_shape(pos):
275 rel_pos = pos - center
276 projection = rel_pos - direction * ta.dot(direction, rel_pos)
277 return sum(projection * projection) <= r_squared
279 return self.shape(wire_shape, center)(sym)
281 return wire_sites
283 def neighbors(self, n=1, eps=1e-8):
284 """Return n-th nearest neighbor hoppings.
286 Parameters
287 ----------
288 n : integer
289 Order of the hoppings to return. Note that the zeroth neighbor is
290 the site itself or any other sites with the same position.
291 eps : float
292 Tolerance relative to the length of the shortest lattice vector for
293 when to consider lengths to be approximately equal.
295 Returns
296 -------
297 hoppings : list of kwant.builder.HoppingKind objects
298 The n-th nearest neighbor hoppings.
300 Notes
301 -----
302 The hoppings are ordered lexicographically according to sublattice from
303 which they originate, sublattice on which they end, and their lattice
304 coordinates. Out of the two equivalent hoppings (a hopping and its
305 reverse) only the lexicographically larger one is returned.
306 """
307 # This algorithm is not designed to be fast and can be improved,
308 # however there is no real need.
309 Site = system.Site
310 sls = self.sublattices
311 shortest_hopping = sls[0].n_closest(
312 sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
313 rtol = eps
314 eps *= np.linalg.norm(self.vec(shortest_hopping))
315 nvec = len(self._prim_vecs)
316 sublat_pairs = [(i, j) for (i, j) in product(sls, sls)
317 if sls.index(j) >= sls.index(i)]
318 def first_nonnegative(tag):
319 for i in tag:
320 if i < 0:
321 return False
322 elif i > 0:
323 return True
324 else:
325 continue
326 return True
328 # Find the `n` closest neighbors (with multiplicity) for each
329 # pair of lattices, this surely includes the `n` closest neighbors overall.
330 sites = []
331 for i, j in sublat_pairs:
332 origin = Site(j, ta.zeros(nvec)).pos
333 tags = i.n_closest(origin, n=n+1, group_by_length=True, rtol=rtol)
334 ij_dist = [np.linalg.norm(Site(i, tag).pos - origin)
335 for tag in tags]
336 sites.append((tags, (j, i), ij_dist))
337 distances = np.r_[tuple((i[2] for i in sites))]
338 distances = np.sort(distances)
339 group_boundaries = np.where(np.diff(distances) > eps)[0]
340 # Find distance in `n`-th group.
341 if len(group_boundaries) == n:
342 n_dist = distances[-1]
343 else:
344 n_dist = distances[group_boundaries[n]]
346 # We now have all the required sites, select the ones in `n`-th group.
347 result = []
348 for group in sites:
349 tags, distance = group[0], group[2]
350 i, j = group[1]
351 tags = np.array([tag for tag, dist in zip(tags, distance)
352 if abs(dist - n_dist) < eps])
353 if len(tags):
354 # Sort the tags.
355 tags = tags[np.lexsort(tags.T[::-1])][::-1]
356 # Throw away equivalent hoppings if
357 # two sublattices are the same.
358 if i == j and len(tags) > 1:
359 tags = tags[: len(tags) // 2]
360 for tag in tags:
361 result.append(builder.HoppingKind(tag, j, i))
362 return result
364 @property
365 def prim_vecs(self):
366 """(sequence of vectors) Primitive vectors
368 `prim_vecs[i]`` is the `i`-th primitive basis vector of the lattice
369 displacement of the lattice origin from the real space coordinates
370 origin.
371 """
372 return np.array(self._prim_vecs)
374 def vec(self, int_vec):
375 """
376 Return the coordinates of a Bravais lattice vector in real space.
378 Parameters
379 ----------
380 vec : integer vector
382 Returns
383 -------
384 output : real vector
385 """
386 return ta.dot(int_vec, self._prim_vecs)
389def short_array_repr(array):
390 full = ' '.join([i.lstrip() for i in repr(array).split('\n')])
391 return full[6 : -1]
394def short_array_str(array):
395 full = ', '.join([i.lstrip() for i in str(array).split('\n')])
396 return full[1 : -1]
399class Monatomic(system.SiteFamily, Polyatomic):
400 """
401 A Bravais lattice with a single site in the basis.
403 Instances of this class provide the `~kwant.system.SiteFamily` interface.
404 Site tags (see `~kwant.system.SiteFamily`) are sequences of integers and
405 describe the lattice coordinates of a site.
407 ``Monatomic`` instances are used as site families on their own or as
408 sublattices of `Polyatomic` lattices.
410 Parameters
411 ----------
412 prim_vecs : 2d array-like of floats
413 Primitive vectors of the Bravais lattice.
415 offset : vector of floats
416 Displacement of the lattice origin from the real space
417 coordinates origin.
419 Attributes
420 ----------
421 ``offset`` : vector
422 Displacement of the lattice origin from the real space coordinates origin
423 """
425 def __init__(self, prim_vecs, offset=None, name='', norbs=None):
426 prim_vecs = ta.array(prim_vecs, float)
427 _check_prim_vecs(prim_vecs)
429 dim = prim_vecs.shape[1]
430 if name is None: 430 ↛ 431line 430 didn't jump to line 431, because the condition on line 430 was never true
431 name = ''
433 if offset is None:
434 offset = ta.zeros(dim)
435 else:
436 offset = ta.array(offset, float)
437 if offset.shape != (dim,): 437 ↛ 438line 437 didn't jump to line 438, because the condition on line 437 was never true
438 raise ValueError('Dimensionality of offset does not match '
439 'that of the space.')
441 msg = '{0}({1}, {2}, {3}, {4})'
442 cl = self.__module__ + '.' + self.__class__.__name__
443 canonical_repr = msg.format(cl, short_array_repr(prim_vecs),
444 short_array_repr(offset),
445 repr(name), repr(norbs))
446 super().__init__(canonical_repr, name, norbs)
448 self.sublattices = [self]
449 self._prim_vecs = prim_vecs
450 self.inv_pv = ta.array(np.linalg.pinv(prim_vecs))
451 self.offset = offset
453 # Precalculation of auxiliary arrays for real space calculations.
454 self.reduced_vecs, self.transf = lll.lll(prim_vecs)
455 self.voronoi = ta.dot(lll.voronoi(self.reduced_vecs), self.transf)
457 self.dim = dim
458 self.lattice_dim = len(prim_vecs)
460 if name != '':
461 msg = "<Monatomic lattice {0}{1}>"
462 orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
463 self.cached_str = msg.format(name, orbs)
464 else:
465 msg = "<unnamed Monatomic lattice, vectors {0}, origin [{1}]{2}>"
466 orbs = ', with {0} orbitals'.format(norbs) if norbs else ''
467 self.cached_str = msg.format(short_array_str(self._prim_vecs),
468 short_array_str(self.offset), orbs)
470 def __str__(self):
471 return self.cached_str
473 def normalize_tags(self, tags):
474 tags = np.asarray(tags, int)
475 tags.flags.writeable = False
476 if tags.shape[1] != self.lattice_dim: 476 ↛ 477line 476 didn't jump to line 477, because the condition on line 476 was never true
477 raise ValueError("Dimensionality mismatch.")
478 return tags
480 def normalize_tag(self, tag):
481 tag = ta.array(tag, int)
482 if len(tag) != self.lattice_dim:
483 raise ValueError("Dimensionality mismatch.")
484 return tag
486 def n_closest(self, pos, n=1, group_by_length=False, rtol=1e-9):
487 """Find n sites closest to position `pos`.
489 Returns
490 -------
491 sites : numpy array
492 An array with sites coordinates.
493 """
494 # TODO (Anton): transform to tinyarrays, once ta indexing is better.
495 return np.dot(lll.cvp(pos - self.offset, self.reduced_vecs,
496 n=n, group_by_length=group_by_length, rtol=rtol),
497 self.transf.T)
499 def closest(self, pos):
500 """
501 Find the lattice coordinates of the site closest to position ``pos``.
502 """
503 return ta.array(self.n_closest(pos)[0])
505 def positions(self, tags):
506 """Return the real-space positions of the sites with the given tags."""
507 return tags @ self._prim_vecs + self.offset
509 def pos(self, tag):
510 """Return the real-space position of the site with a given tag."""
511 return ta.dot(tag, self._prim_vecs) + self.offset
514# The following class is designed such that it should avoid floating
515# point precision issues.
517class TranslationalSymmetry(system.Symmetry):
518 """A translational symmetry defined in real space.
520 An alias exists for this common name: ``kwant.TranslationalSymmetry``.
522 Group elements of this symmetry are integer tuples of appropriate length.
524 Parameters
525 ----------
526 p0, p1, p2, ... : sequences of real numbers
527 The symmetry periods in real space.
529 Notes
530 -----
531 This symmetry automatically chooses the fundamental domain for each new
532 `SiteFamily` it encounters. If this site family does not correspond to a
533 Bravais lattice, or if it does not have a commensurate period, an error is
534 produced. A certain flexibility in choice of the fundamental domain can be
535 achieved by calling manually the `add_site_family` method and providing it
536 the `other_vectors` parameter.
538 The fundamental domain for hoppings are all hoppings ``(a, b)`` with site
539 `a` in fundamental domain of sites.
540 """
541 def __init__(self, *periods):
542 self.periods = ta.array(periods)
543 if self.periods.ndim != 2: 543 ↛ 544line 543 didn't jump to line 544, because the condition on line 543 was never true
544 msg = ("TranslationalSymmetry takes 1d sequences as parameters.")
545 raise ValueError(msg)
546 if np.linalg.matrix_rank(periods) < len(periods):
547 raise ValueError("Translational symmetry periods must be "
548 "linearly independent")
549 # A dictionary containing cached data required for applying the
550 # symmetry to different site families.
551 self.site_family_data = {}
552 self.is_reversed = False
554 def subgroup(self, *generators):
555 """Return the subgroup generated by a sequence of group elements.
557 Parameters
558 ----------
559 *generators: sequence of int
560 Each generator must have length ``self.num_directions``.
561 """
562 generators = ta.array(generators)
563 if generators.dtype != int:
564 raise ValueError('Generators must be sequences of integers.')
565 return TranslationalSymmetry(*ta.dot(generators, self.periods))
567 def has_subgroup(self, other):
568 if isinstance(other, system.NoSymmetry):
569 return True
570 elif not isinstance(other, TranslationalSymmetry): 570 ↛ 571line 570 didn't jump to line 571, because the condition on line 570 was never true
571 raise ValueError("Unknown symmetry type.")
573 if other.periods.shape[1] != self.periods.shape[1]: 573 ↛ 574line 573 didn't jump to line 574, because the condition on line 573 was never true
574 return False # Mismatch of spatial dimensionalities.
576 inv = np.linalg.pinv(self.periods)
577 factors = np.dot(other.periods, inv)
578 # Absolute tolerance is correct in the following since we want an error
579 # relative to the closest integer.
580 return (np.allclose(factors, np.round(factors), rtol=0, atol=1e-8) and
581 np.allclose(ta.dot(factors, self.periods), other.periods))
583 def add_site_family(self, fam, other_vectors=None):
584 """
585 Select a fundamental domain for site family and cache associated data.
587 Parameters
588 ----------
589 fam : `SiteFamily`
590 the site family which has to be processed. Be sure to delete the
591 previously processed site families from `site_family_data` if you
592 want to modify the cache.
594 other_vectors : 2d array-like of integers
595 Bravais lattice vectors used to complement the periods in forming
596 a basis. The fundamental domain consists of all the lattice sites
597 for which the zero coefficients corresponding to the symmetry
598 periods in the basis formed by the symmetry periods and
599 `other_vectors`. If an insufficient number of `other_vectors` is
600 provided to form a basis, the missing ones are selected
601 automatically.
603 Raises
604 ------
605 KeyError
606 If `fam` is already stored in `site_family_data`.
607 ValueError
608 If lattice `fam` is incompatible with given periods.
609 """
610 ensure_isinstance(fam, Monatomic)
612 dim = self.periods.shape[1]
613 if fam in self.site_family_data: 613 ↛ 614line 613 didn't jump to line 614, because the condition on line 613 was never true
614 raise KeyError('Family already processed, delete it from '
615 'site_family_data first.')
616 inv = np.linalg.pinv(fam.prim_vecs)
617 try:
618 bravais_periods = np.dot(self.periods, inv)
619 except ValueError:
620 fam_space_dim = fam.prim_vecs.shape[1]
621 if dim == fam_space_dim:
622 raise
623 msg = ("{0}-d-embedded lattice is incompatible with "
624 "{1}-d translational symmetry.")
625 raise ValueError(msg.format(fam_space_dim, dim))
626 # Absolute tolerance is correct in the following since we want an error
627 # relative to the closest integer.
628 if (not np.allclose(bravais_periods, np.round(bravais_periods),
629 rtol=0, atol=1e-8) or
630 not np.allclose([fam.vec(i) for i in bravais_periods],
631 self.periods)):
632 msg = ('Site family {0} does not have commensurate periods with '
633 'symmetry {1}.')
634 raise ValueError(msg.format(fam, self))
635 bravais_periods = np.array(np.round(bravais_periods), dtype='int')
636 (num_dir, lat_dim) = bravais_periods.shape
637 if other_vectors is None:
638 other_vectors = np.zeros((0, lat_dim), dtype=int)
639 else:
640 other_vectors = np.array(other_vectors)
641 if other_vectors.ndim != 2: 641 ↛ 642line 641 didn't jump to line 642, because the condition on line 641 was never true
642 raise ValueError(
643 '`other_vectors` must be a 2d array-like object.')
644 if np.any(np.round(other_vectors) - other_vectors): 644 ↛ 645line 644 didn't jump to line 645, because the condition on line 644 was never true
645 raise ValueError('Only integer other_vectors are allowed.')
646 other_vectors = np.array(np.round(other_vectors), dtype=int)
648 m = np.zeros((lat_dim, lat_dim), dtype=int)
650 m.T[:num_dir] = bravais_periods
651 num_vec = num_dir + len(other_vectors)
652 m.T[num_dir:num_vec] = other_vectors
654 if np.linalg.matrix_rank(m) < num_vec: 654 ↛ 655line 654 didn't jump to line 655, because the condition on line 654 was never true
655 raise ValueError('other_vectors and symmetry periods are not '
656 'linearly independent.')
658 # To define the fundamental domain of the new site family we now need to
659 # choose `lat_dim - num_vec` extra lattice vectors that are not
660 # linearly dependent on the vectors we already have. To do so we
661 # continuously add the lattice basis vectors one by one such that they
662 # are not linearly dependent on the existent vectors
663 while num_vec < lat_dim:
664 vh = np.linalg.svd(np.dot(m[:, :num_vec].T, fam.prim_vecs),
665 full_matrices=False)[2]
666 projector = np.identity(dim) - np.dot(vh.T, vh)
668 residuals = np.dot(fam.prim_vecs, projector)
669 residuals = np.apply_along_axis(np.linalg.norm, 1, residuals)
670 m[np.argmax(residuals), num_vec] = 1
671 num_vec += 1
673 det_m = int(round(np.linalg.det(m)))
674 if det_m == 0: 674 ↛ 675line 674 didn't jump to line 675, because the condition on line 674 was never true
675 raise RuntimeError('Adding site family failed.')
677 det_x_inv_m = np.array(np.round(det_m * np.linalg.inv(m)), dtype=int)
678 assert (np.dot(m, det_x_inv_m) // det_m == np.identity(lat_dim)).all()
680 det_x_inv_m_part = det_x_inv_m[:num_dir, :]
681 m_part = m[:, :num_dir]
682 self.site_family_data[fam] = (ta.array(m_part),
683 ta.array(det_x_inv_m_part), det_m)
685 @property
686 def num_directions(self):
687 return len(self.periods)
689 def _get_site_family_data(self, family):
690 try:
691 return self.site_family_data[family]
692 except KeyError:
693 self.add_site_family(family)
694 return self.site_family_data[family]
696 def which(self, site):
697 det_x_inv_m_part, det_m = self._get_site_family_data(site.family)[-2:]
698 if isinstance(site, system.Site): 698 ↛ 700line 698 didn't jump to line 700, because the condition on line 698 was never false
699 result = ta.dot(det_x_inv_m_part, site.tag) // det_m
700 elif isinstance(site, system.SiteArray):
701 result = np.dot(det_x_inv_m_part, site.tags.transpose()) // det_m
702 else:
703 raise TypeError("'site' must be a Site or a SiteArray")
705 return -result if self.is_reversed else result
707 def act(self, element, a, b=None):
708 is_site = isinstance(a, system.Site)
709 # Tinyarray for small arrays (single site) else numpy
710 array_mod = ta if is_site else np
711 element = array_mod.array(element)
712 if not np.issubdtype(element.dtype, np.integer):
713 raise ValueError("group element must be a tuple of integers")
714 if (len(element.shape) == 2 and is_site): 714 ↛ 715line 714 didn't jump to line 715, because the condition on line 714 was never true
715 raise ValueError("must provide a single group element when "
716 "acting on single sites.")
717 if (len(element.shape) == 1 and not is_site):
718 # We can act on a whole SiteArray with a single group element
719 # using numpy broadcasting.
720 element = element.reshape(-1, 1)
721 m_part = self._get_site_family_data(a.family)[0]
722 try:
723 delta = array_mod.dot(m_part, element)
724 except ValueError:
725 msg = 'Expecting a {0}-tuple group element, but got `{1}` instead.'
726 raise ValueError(msg.format(self.num_directions, element))
727 if self.is_reversed:
728 delta = -delta
729 if b is None:
730 if is_site:
731 return system.Site(a.family, a.tag + delta, True)
732 else:
733 return system.SiteArray(a.family, a.tags + delta.transpose())
734 elif b.family == a.family:
735 if is_site: 735 ↛ 739line 735 didn't jump to line 739, because the condition on line 735 was never false
736 return (system.Site(a.family, a.tag + delta, True),
737 system.Site(b.family, b.tag + delta, True))
738 else:
739 return (system.SiteArray(a.family, a.tags + delta.transpose()),
740 system.SiteArray(b.family, b.tags + delta.transpose()))
741 else:
742 m_part = self._get_site_family_data(b.family)[0]
743 try:
744 delta2 = ta.dot(m_part, element)
745 except ValueError:
746 msg = ('Expecting a {0}-tuple group element, '
747 'but got `{1}` instead.')
748 raise ValueError(msg.format(self.num_directions, element))
749 if self.is_reversed: 749 ↛ 750line 749 didn't jump to line 750, because the condition on line 749 was never true
750 delta2 = -delta2
751 if is_site: 751 ↛ 755line 751 didn't jump to line 755, because the condition on line 751 was never false
752 return (system.Site(a.family, a.tag + delta, True),
753 system.Site(b.family, b.tag + delta2, True))
754 else:
755 return (system.SiteArray(a.family, a.tags + delta.transpose()),
756 system.SiteArray(b.family, b.tags + delta2.transpose()))
758 def reversed(self):
759 """Return a reversed copy of the symmetry.
761 The resulting symmetry has all the period vectors opposite to the
762 original and an identical fundamental domain.
763 """
764 result = TranslationalSymmetry(*self.periods)
765 result.site_family_data = self.site_family_data
766 result.is_reversed = not self.is_reversed
767 result.periods = -self.periods
768 return result
772################ Library of lattices
774def chain(a=1, name='', norbs=None):
775 """Make a one-dimensional lattice."""
776 return Monatomic(((a,),), name=name, norbs=norbs)
779def square(a=1, name='', norbs=None):
780 """Make a square lattice."""
781 return Monatomic(((a, 0), (0, a)), name=name, norbs=norbs)
784def cubic(a=1, name='', norbs=None):
785 """Make a cubic lattice."""
786 return Monatomic(((a, 0, 0), (0, a, 0), (0, 0, a)),
787 name=name, norbs=norbs)
790tri = ta.array(((1, 0), (0.5, 0.5 * sqrt(3))))
792def triangular(a=1, name='', norbs=None):
793 """Make a triangular lattice."""
794 return Monatomic(a * tri, name=name, norbs=norbs)
797def honeycomb(a=1, name='', norbs=None):
798 """Make a honeycomb lattice."""
799 lat = Polyatomic(a * tri, ((0, 0), (0, a / sqrt(3))),
800 name=name, norbs=norbs)
801 lat.a, lat.b = lat.sublattices
802 return lat
805def kagome(a=1, name='', norbs=None):
806 """Make a kagome lattice."""
807 lat = Polyatomic(a * tri, ((0, 0),) + tuple(0.5 * a * tri),
808 name=name, norbs=norbs)
809 lat.a, lat.b, lat.c = lat.sublattices
810 return lat
813# TODO (Anton): unhide _prim_vecs, once tinyarray supports indexing of
814# sub-arrays.