Coverage for kwant/qsymm.py : 92%
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-2018 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.
10__all__ = ['builder_to_model', 'model_to_builder', 'find_builder_symmetries']
12import itertools as it
13from collections import OrderedDict, defaultdict
15import numpy as np
16import tinyarray as ta
17import scipy.linalg as la
19try:
20 import sympy
21 import sympy.matrices.matrices
22 import qsymm
23 from qsymm.model import Model, BlochModel, BlochCoeff
24 from qsymm.groups import PointGroupElement, ContinuousGroupGenerator
25 from qsymm.symmetry_finder import bravais_point_group
26 from qsymm.linalg import allclose
27 from qsymm.hamiltonian_generator import hamiltonian_from_family
29 one = sympy.S.One
30except ImportError as error:
31 msg = ("'kwant.qsymm' is not available because one or more of its "
32 "dependencies is not installed.")
33 raise ImportError(msg) from error
35from kwant import lattice, builder
36from kwant._common import get_parameters
39def builder_to_model(syst, momenta=None, real_space=True,
40 params=None):
41 """Make a qsymm.BlochModel out of a `~kwant.builder.Builder`.
43 Parameters
44 ----------
45 syst : `~kwant.builder.Builder`
46 May have translational symmetries.
47 momenta : list of strings or None
48 Names of momentum variables. If None, 'k_x', 'k_y', ... is used.
49 real_space : bool (default True)
50 If False, use the unit cell convention for Bloch basis, the
51 exponential has the difference in the unit cell coordinates and
52 k is expressed in the reciprocal lattice basis. This is consistent
53 with `kwant.wraparound`.
54 If True, the difference in the real space coordinates is used
55 and k is given in an absolute basis.
56 Only the default choice guarantees that qsymm is able to find
57 nonsymmorphic symmetries.
58 params : dict, optional
59 Dictionary of parameter names and their values; used when
60 evaluating the Hamiltonian matrix elements.
62 Returns
63 -------
64 model : qsymm.BlochModel
65 Model representing the tight-binding Hamiltonian.
67 Notes
68 -----
69 The sites in the the builder are in lexicographical order, i.e. ordered
70 first by their family and then by their tag. This is the same ordering that
71 is used in finalized kwant systems.
72 """
73 def term_to_model(d, par, matrix):
74 if allclose(matrix, 0):
75 result = BlochModel({}, shape=matrix.shape, format=np.ndarray)
76 else:
77 result = BlochModel({BlochCoeff(d, qsymm.sympify(par)): matrix},
78 momenta=momenta)
79 return result
81 def hopping_to_model(hop, value, proj, params):
82 site1, site2 = hop
83 if real_space:
84 d = proj @ np.array(site2.pos - site1.pos)
85 else:
86 # site in the FD
87 d = np.array(syst.symmetry.which(site2))
89 slice1, slice2 = slices[to_fd(site1)], slices[to_fd(site2)]
90 if callable(value):
91 return sum(term_to_model(d, par, set_block(slice1, slice2, val))
92 for par, val in function_to_terms(hop, value, params))
93 else:
94 matrix = set_block(slice1, slice2, value)
95 return term_to_model(d, '1', matrix)
97 def onsite_to_model(site, value, params):
98 d = np.zeros((dim, ))
99 slice1 = slices[to_fd(site)]
100 if callable(value):
101 return sum(term_to_model(d, par, set_block(slice1, slice1, val))
102 for par, val in function_to_terms(site, value, params))
103 else:
104 return term_to_model(d, '1', set_block(slice1, slice1, value))
106 def function_to_terms(site_or_hop, value, fixed_params):
107 assert callable(value)
108 parameters = get_parameters(value)
109 # remove site or site1, site2 parameters
110 if isinstance(site_or_hop, builder.Site):
111 parameters = parameters[1:]
112 site_or_hop = (site_or_hop,)
113 else:
114 parameters = parameters[2:]
115 free_parameters = (par for par in parameters
116 if par not in fixed_params.keys())
117 # first set all free parameters to 0
118 args = ((fixed_params[par] if par in fixed_params.keys() else 0)
119 for par in parameters)
120 h_0 = value(*site_or_hop, *args)
121 # set one of the free parameters to 1 at a time, the rest 0
122 terms = []
123 for p in free_parameters:
124 args = ((fixed_params[par] if par in fixed_params.keys() else
125 (1 if par == p else 0)) for par in parameters)
126 terms.append((p, value(*site_or_hop, *args) - h_0))
127 return terms + [('1', h_0)]
129 def orbital_slices(syst):
130 orbital_slices = {}
131 start_orb = 0
133 for site in sorted(syst.sites()):
134 n = site.family.norbs
135 if n is None: 135 ↛ 136line 135 didn't jump to line 136, because the condition on line 135 was never true
136 raise ValueError('norbs must be provided for every lattice.')
137 orbital_slices[site] = slice(start_orb, start_orb + n)
138 start_orb += n
139 return orbital_slices, start_orb
141 def set_block(slice1, slice2, val):
142 matrix = np.zeros((N, N), dtype=complex)
143 matrix[slice1, slice2] = val
144 return matrix
146 if params is None:
147 params = dict()
149 periods = np.array(syst.symmetry.periods)
150 dim = len(periods)
151 to_fd = syst.symmetry.to_fd
152 if momenta is None:
153 momenta = ['k_x', 'k_y', 'k_z'][:dim]
154 # If the system is higher dimensional than the number of translation
155 # vectors, we need to project onto the subspace spanned by the
156 # translation vectors.
157 if dim == 0:
158 proj = np.empty((0, len(list(syst.sites())[0].pos)))
159 elif dim < len(list(syst.sites())[0].pos):
160 proj, r = la.qr(np.array(periods).T, mode='economic')
161 sign = np.diag(np.diag(np.sign(r)))
162 proj = sign @ proj.T
163 else:
164 proj = np.eye(dim)
166 slices, N = orbital_slices(syst)
168 one_way_hoppings = [hopping_to_model(hop, value, proj, params)
169 for hop, value in syst.hopping_value_pairs()]
170 other_way_hoppings = [term.T().conj() for term in one_way_hoppings]
171 hoppings = one_way_hoppings + other_way_hoppings
173 onsites = [onsite_to_model(site, value, params)
174 for site, value in syst.site_value_pairs()]
176 result = sum(onsites) + sum(hoppings)
178 return result
181def model_to_builder(model, norbs, lat_vecs, atom_coords, *, coeffs=None):
182 """Make a `~kwant.builder.Builder` out of qsymm.Models or qsymm.BlochModels.
184 Parameters
185 ----------
186 model : qsymm.Model, qsymm.BlochModel, or an iterable thereof
187 The Hamiltonian (or terms of the Hamiltonian) to convert to a
188 Builder.
189 norbs : OrderedDict or sequence of pairs
190 Maps sites to the number of orbitals per site in a unit cell.
191 lat_vecs : list of arrays
192 Lattice vectors of the underlying tight binding lattice.
193 atom_coords : list of arrays
194 Positions of the sites (or atoms) within a unit cell.
195 The ordering of the atoms is the same as in norbs.
196 coeffs : list of sympy.Symbol, default None.
197 Constant prefactors for the individual terms in model, if model
198 is a list of multiple objects. If model is a single Model or BlochModel
199 object, this argument is ignored. By default assigns the coefficient
200 c_n to element model[n].
202 Returns
203 -------
204 syst : `~kwant.builder.Builder`
205 The unfinalized Kwant system representing the qsymm Model(s).
207 Notes
208 -----
209 Onsite terms that are not provided in the input model are set
210 to zero by default.
212 The input model(s) representing the tight binding Hamiltonian in
213 Bloch form should follow the convention where the difference in the real
214 space atomic positions appear in the Bloch factors.
215 """
217 def make_int(R):
218 # If close to an integer array convert to integer tinyarray, else
219 # return None
220 R_int = ta.array(np.round(R), int)
221 if qsymm.linalg.allclose(R, R_int): 221 ↛ 224line 221 didn't jump to line 224, because the condition on line 221 was never false
222 return R_int
223 else:
224 return None
226 def term_onsite(onsites_dict, hopping_dict, hop_mat, atoms,
227 sublattices, coords_dict):
228 """Find the Kwant onsites and hoppings in a qsymm.BlochModel term
229 that has no lattice translation in the Bloch factor.
230 """
231 for atom1, atom2 in it.product(atoms, atoms):
232 # Subblock within the same sublattice is onsite
233 hop = hop_mat[ranges[atom1], ranges[atom2]]
234 if sublattices[atom1] == sublattices[atom2]:
235 onsites_dict[atom1] += Model({coeff: hop}, momenta=momenta)
236 # Blocks between sublattices are hoppings between sublattices
237 # at the same position.
238 # Only include nonzero hoppings
239 elif not allclose(hop, 0):
240 if not allclose(np.array(coords_dict[atom1]), 240 ↛ 242line 240 didn't jump to line 242, because the condition on line 240 was never true
241 np.array(coords_dict[atom2])):
242 raise ValueError(
243 "Position of sites not compatible with qsymm model.")
244 lat_basis = np.array(zer)
245 hop = Model({coeff: hop}, momenta=momenta)
246 hop_dir = builder.HoppingKind(-lat_basis, sublattices[atom1],
247 sublattices[atom2])
248 hopping_dict[hop_dir] += hop
249 return onsites_dict, hopping_dict
251 def term_hopping(hopping_dict, hop_mat, atoms,
252 sublattices, coords_dict):
253 """Find Kwant hoppings in a qsymm.BlochModel term that has a lattice
254 translation in the Bloch factor.
255 """
256 # Iterate over combinations of atoms, set hoppings between each
257 for atom1, atom2 in it.product(atoms, atoms):
258 # Take the block from atom1 to atom2
259 hop = hop_mat[ranges[atom1], ranges[atom2]]
260 # Only include nonzero hoppings
261 if allclose(hop, 0):
262 continue
263 # Adjust hopping vector to Bloch form basis
264 r_lattice = (
265 r_vec
266 + np.array(coords_dict[atom1])
267 - np.array(coords_dict[atom2])
268 )
269 # Bring vector to basis of lattice vectors
270 lat_basis = np.linalg.solve(np.vstack(lat_vecs).T, r_lattice)
271 lat_basis = make_int(lat_basis)
272 # Should only have hoppings that are integer multiples of
273 # lattice vectors
274 if lat_basis is not None: 274 ↛ 281line 274 didn't jump to line 281, because the condition on line 274 was never false
275 hop_dir = builder.HoppingKind(-lat_basis,
276 sublattices[atom1],
277 sublattices[atom2])
278 # Set the hopping as the matrix times the hopping amplitude
279 hopping_dict[hop_dir] += Model({coeff: hop}, momenta=momenta)
280 else:
281 raise RuntimeError('A nonzero hopping not matching a '
282 'lattice vector was found.')
283 return hopping_dict
285 # Disambiguate single model instances from iterables thereof. Because
286 # Model is itself iterable (subclasses dict) this is a bit cumbersome.
287 if isinstance(model, Model):
288 # BlochModel can't yet handle getting a Blochmodel as input
289 if not isinstance(model, BlochModel):
290 model = BlochModel(model)
291 else:
292 model = BlochModel(hamiltonian_from_family(
293 model, coeffs=coeffs, nsimplify=False, tosympy=False))
296 # 'momentum' and 'zer' are used in the closures defined above, so don't
297 # move these declarations down.
298 momenta = model.momenta
299 if len(momenta) != len(lat_vecs): 299 ↛ 300line 299 didn't jump to line 300, because the condition on line 299 was never true
300 raise ValueError("Dimension of the lattice and number of "
301 "momenta do not match.")
302 zer = [0] * len(momenta)
305 # Subblocks of the Hamiltonian for different atoms.
306 N = 0
307 if not any([isinstance(norbs, OrderedDict), isinstance(norbs, list), 307 ↛ 309line 307 didn't jump to line 309, because the condition on line 307 was never true
308 isinstance(norbs, tuple)]):
309 raise ValueError('norbs must be OrderedDict, tuple, or list.')
310 else:
311 norbs = OrderedDict(norbs)
312 ranges = dict()
313 for a, n in norbs.items():
314 ranges[a] = slice(N, N + n)
315 N += n
317 # Extract atoms and number of orbitals per atom,
318 # store the position of each atom
319 atoms, orbs = zip(*norbs.items())
320 coords_dict = dict(zip(atoms, atom_coords))
322 # Make the kwant lattice
323 lat = lattice.general(lat_vecs, atom_coords, norbs=orbs)
324 # Store sublattices by name
325 sublattices = dict(zip(atoms, lat.sublattices))
327 # Keep track of the hoppings and onsites by storing those
328 # which have already been set.
329 hopping_dict = defaultdict(lambda: 0)
330 onsites_dict = defaultdict(lambda: 0)
332 # Iterate over all terms in the model.
333 for key, hop_mat in model.items():
334 # Determine whether this term is an onsite or a hopping, extract
335 # overall symbolic coefficient if any, extract the exponential
336 # part describing the hopping if present.
337 r_vec, coeff = key
338 # Onsite term; modifies onsites_dict and hopping_dict in-place
339 if allclose(r_vec, 0):
340 term_onsite(
341 onsites_dict, hopping_dict, hop_mat,
342 atoms, sublattices, coords_dict)
343 # Hopping term; modifies hopping_dict in-place
344 else:
345 term_hopping(hopping_dict, hop_mat, atoms,
346 sublattices, coords_dict)
348 # If some onsite terms are not set, we set them to zero.
349 for atom in atoms:
350 if atom not in onsites_dict:
351 onsites_dict[atom] = Model(
352 {one: np.zeros((norbs[atom], norbs[atom]))},
353 momenta=momenta)
355 # Make the Kwant system, and set all onsites and hoppings.
357 sym = lattice.TranslationalSymmetry(*lat_vecs)
358 syst = builder.Builder(sym)
360 # Iterate over all onsites and set them
361 for atom, onsite in onsites_dict.items():
362 syst[sublattices[atom](*zer)] = onsite.lambdify(onsite=True)
364 # Finally, iterate over all the hoppings and set them
365 for direction, hopping in hopping_dict.items():
366 syst[direction] = hopping.lambdify(hopping=True)
368 return syst
371# This may be useful in the future, so we'll keep it as internal for now,
372# and can make it part of the API in the future if we wish.
373def _get_builder_symmetries(builder):
374 """Extract the declared symmetries of a Kwant builder.
376 Parameters
377 ----------
378 builder : `~kwant.builder.Builder`
380 Returns
381 -------
382 builder_symmetries : dict
383 Dictionary of the discrete symmetries that the builder has.
384 The symmetries can be particle-hole, time-reversal or chiral,
385 which are returned as qsymm.PointGroupElements, or
386 a conservation law, which is returned as a
387 qsymm.ContinuousGroupGenerators.
388 """
390 dim = len(np.array(builder.symmetry.periods))
391 symmetry_names = ['time_reversal', 'particle_hole', 'chiral',
392 'conservation_law']
393 builder_symmetries = {name: getattr(builder, name)
394 for name in symmetry_names
395 if getattr(builder, name) is not None}
396 for name, symmetry in builder_symmetries.items():
397 if name == 'time_reversal': 397 ↛ 398line 397 didn't jump to line 398, because the condition on line 397 was never true
398 builder_symmetries[name] = PointGroupElement(np.eye(dim),
399 True, False, symmetry)
400 elif name == 'particle_hole':
401 builder_symmetries[name] = PointGroupElement(np.eye(dim),
402 True, True, symmetry)
403 elif name == 'chiral': 403 ↛ 404line 403 didn't jump to line 404, because the condition on line 403 was never true
404 builder_symmetries[name] = PointGroupElement(np.eye(dim),
405 False, True, symmetry)
406 elif name == 'conservation_law': 406 ↛ 410line 406 didn't jump to line 410, because the condition on line 406 was never false
407 builder_symmetries[name] = ContinuousGroupGenerator(R=None,
408 U=symmetry)
409 else:
410 raise ValueError("Invalid symmetry name.")
411 return builder_symmetries
414def find_builder_symmetries(builder, momenta=None, params=None,
415 spatial_symmetries=True, prettify=True,
416 sparse=None):
417 """Finds the symmetries of a Kwant system using qsymm.
419 Parameters
420 ----------
421 builder : `~kwant.builder.Builder`
422 momenta : list of strings or None
423 Names of momentum variables, if None 'k_x', 'k_y', ... is used.
424 params : dict, optional
425 Dictionary of parameter names and their values; used when
426 evaluating the Hamiltonian matrix elements.
427 spatial_symmetries : bool (default True)
428 If True, search for all symmetries.
429 If False, only searches for the symmetries that are declarable in
430 `~kwant.builder.Builder` objects, i.e. time-reversal symmetry,
431 particle-hole symmetry, chiral symmetry, or conservation laws.
432 This can save computation time.
433 prettify : bool (default True)
434 Whether to carry out sparsification of the continuous symmetry
435 generators, in general an arbitrary linear combination of the
436 symmetry generators is returned.
437 sparse : bool, or None (default None)
438 Whether to use sparse linear algebra in the calculation.
439 Can give large performance gain in large systems.
440 If None, uses sparse or dense computation depending on
441 the size of the Hamiltonian.
444 Returns
445 -------
446 symmetries : list of qsymm.PointGroupElements and/or qsymm.ContinuousGroupElement
447 The symmetries of the Kwant system.
448 """
450 if params is None: 450 ↛ 453line 450 didn't jump to line 453, because the condition on line 450 was never false
451 params = dict()
453 ham = builder_to_model(builder, momenta=momenta,
454 real_space=spatial_symmetries, params=params)
456 # Use dense or sparse computation depending on Hamiltonian size
457 if sparse is None:
458 sparse = list(ham.values())[0].shape[0] > 20
460 dim = len(np.array(builder.symmetry.periods))
462 if spatial_symmetries:
463 candidates = bravais_point_group(builder.symmetry.periods, tr=True,
464 ph=True, generators=False,
465 verbose=False)
466 else:
467 candidates = [
468 qsymm.PointGroupElement(np.eye(dim), True, False, None), # T
469 qsymm.PointGroupElement(np.eye(dim), True, True, None), # P
470 qsymm.PointGroupElement(np.eye(dim), False, True, None)] # C
471 sg, cg = qsymm.symmetries(ham, candidates, prettify=prettify,
472 continuous_rotations=False,
473 sparse_linalg=sparse)
474 return list(sg) + list(cg)