From 32a380f794c997c4b4a238e6fb58dca4761f0c96 Mon Sep 17 00:00:00 2001
From: Christoph Groth <christoph.groth@cea.fr>
Date: Tue, 25 Jul 2017 15:24:47 +0200
Subject: [PATCH] builder: move finalization into initializers of ll-system
classes
---
kwant/builder.py | 445 ++++++++++++++++++------------------
kwant/tests/test_builder.py | 6 +-
2 files changed, 224 insertions(+), 227 deletions(-)
diff --git a/kwant/builder.py b/kwant/builder.py
index f5c9f023..1e86c67d 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -576,7 +576,7 @@ class BuilderLead(Lead):
The order of interface sites is kept during finalization.
"""
- syst = self.builder._finalized_infinite(self.interface)
+ syst = InfiniteSystem(self.builder, self.interface)
_transfer_symmetry(syst, self.builder)
return syst
@@ -1576,9 +1576,9 @@ class Builder:
`Symmetry` can be finalized.
"""
if self.symmetry.num_directions == 0:
- syst = self._finalized_finite()
+ syst = FiniteSystem(self)
elif self.symmetry.num_directions == 1:
- syst = self._finalized_infinite()
+ syst = InfiniteSystem(self)
else:
raise ValueError('Currently, only builders without or with a 1D '
'translational symmetry can be finalized.')
@@ -1587,11 +1587,113 @@ class Builder:
return syst
- def _finalized_finite(self):
- assert self.symmetry.num_directions == 0
+
+
+################ Finalized systems
+
+def _raise_user_error(exc, func):
+ msg = ('Error occurred in user-supplied value function "{0}".\n'
+ 'See the upper part of the above backtrace for more information.')
+ raise UserCodeError(msg.format(func.__name__)) from exc
+
+
+def discrete_symmetry(self, args=(), *, params=None):
+ if self._cons_law is not None:
+ eigvals, eigvecs = self._cons_law
+ eigvals = eigvals.tocoo(args, params=params)
+ if not np.allclose(eigvals.data, np.round(eigvals.data)):
+ raise ValueError("Conservation law must have integer eigenvalues.")
+ eigvals = np.round(eigvals).tocsr()
+ # Avoid appearance of zero eigenvalues
+ eigvals = eigvals + 0.5 * sparse.identity(eigvals.shape[0])
+ eigvals.eliminate_zeros()
+ eigvecs = eigvecs.tocoo(args, params=params)
+ projectors = [eigvecs.dot(eigvals == val)
+ for val in sorted(np.unique(eigvals.data))]
+ else:
+ projectors = None
+
+ def evaluate(op):
+ return None if op is None else op.tocoo(args, params=params)
+
+ return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
+ self._symmetries))
+
+
+def _transfer_symmetry(syst, builder):
+ """Take a symmetry from builder and transfer it to finalized system."""
+ def operator(op):
+ if op is None:
+ return
+ return Density(syst, op, check_hermiticity=False)
+
+ # Conservation law requires preprocessing to split it into eigenvectors
+ # and eigenvalues.
+ cons_law = builder.conservation_law
+
+ if cons_law is None:
+ syst._cons_law = None
+
+ elif callable(cons_law):
+ @wraps(cons_law)
+ def vals(site, *args, **kwargs):
+ if site.family.norbs == 1:
+ return cons_law(site, *args, **kwargs)
+ return np.diag(np.linalg.eigvalsh(cons_law(site, *args, **kwargs)))
+
+ @wraps(cons_law)
+ def vecs(site, *args, **kwargs):
+ if site.family.norbs == 1:
+ return 1
+ return np.linalg.eigh(cons_law(site, *args, **kwargs))[1]
+
+ syst._cons_law = operator(vals), operator(vecs)
+
+ elif isinstance(cons_law, collections.Mapping):
+ vals = {family: (value if family.norbs == 1 else
+ ta.array(np.diag(np.linalg.eigvalsh(value))))
+ for family, value in cons_law.items()}
+ vecs = {family: (1 if family.norbs == 1 else
+ ta.array(np.linalg.eigh(value)[1]))
+ for family, value in cons_law.items()}
+ syst._cons_law = operator(vals), operator(vecs)
+
+ else:
+ try:
+ vals, vecs = np.linalg.eigh(cons_law)
+ vals = np.diag(vals)
+ except np.linalg.LinAlgError as e:
+ if '0-dimensional' not in e.args[0]:
+ raise e # skip coverage
+ vals, vecs = cons_law, 1
+
+ syst._cons_law = operator(vals), operator(vecs)
+
+ syst._symmetries = [operator(symm) for symm in (builder.time_reversal,
+ builder.particle_hole,
+ builder.chiral)]
+
+
+class FiniteSystem(system.FiniteSystem):
+ """Finalized `Builder` with leads.
+
+ Usable as input for the solvers in `kwant.solvers`.
+
+ Attributes
+ ----------
+ sites : sequence
+ ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+ to the integer-labeled site ``i`` of the low-level system. The sites
+ are ordered first by their family and then by their tag.
+ id_by_site : dict
+ The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
+ """
+
+ def __init__(self, builder):
+ assert builder.symmetry.num_directions == 0
#### Make translation tables.
- sites = tuple(sorted(self.H))
+ sites = tuple(sorted(builder.H))
id_by_site = {}
for site_id, site in enumerate(sites):
id_by_site[site] = site_id
@@ -1599,7 +1701,7 @@ class Builder:
#### Make graph.
g = graph.Graph()
g.num_nodes = len(sites) # Some sites could not appear in any edge.
- for tail, hvhv in self.H.items():
+ for tail, hvhv in builder.H.items():
for head in islice(hvhv, 2, None, 2):
if tail == head:
continue
@@ -1609,7 +1711,7 @@ class Builder:
#### Connect leads.
finalized_leads = []
lead_interfaces = []
- for lead_nr, lead in enumerate(self.leads):
+ for lead_nr, lead in enumerate(builder.leads):
try:
with warnings.catch_warnings(record=True) as ws:
warnings.simplefilter("always")
@@ -1638,9 +1740,9 @@ class Builder:
lead_interfaces.append(np.array(interface))
#### Find parameters taken by all value functions
- hoppings = [self._get_edge(sites[tail], sites[head])
+ hoppings = [builder._get_edge(sites[tail], sites[head])
for tail, head in g]
- onsite_hamiltonians = [self.H[site][1] for site in sites]
+ onsite_hamiltonians = [builder.H[site][1] for site in sites]
_ham_param_map = {}
for hams, skip in [(onsite_hamiltonians, 1), (hoppings, 2)]:
@@ -1653,30 +1755,109 @@ class Builder:
params = params[skip:] # remove site argument(s)
_ham_param_map[ham] = (params, defaults, takes_kwargs)
- #### Assemble and return result.
- result = FiniteSystem()
- result.graph = g
- result.sites = sites
- result.site_ranges = _site_ranges(sites)
- result.id_by_site = id_by_site
- result.leads = finalized_leads
- result.hoppings = hoppings
- result.onsite_hamiltonians = onsite_hamiltonians
- result._ham_param_map = _ham_param_map
- result.lead_interfaces = lead_interfaces
- result.symmetry = self.symmetry
- return result
+ self.graph = g
+ self.sites = sites
+ self.site_ranges = _site_ranges(sites)
+ self.id_by_site = id_by_site
+ self.leads = finalized_leads
+ self.hoppings = hoppings
+ self.onsite_hamiltonians = onsite_hamiltonians
+ self._ham_param_map = _ham_param_map
+ self.lead_interfaces = lead_interfaces
+ self.symmetry = builder.symmetry
+
+
+ def hamiltonian(self, i, j, *args, params=None):
+ if args and params:
+ raise TypeError("'args' and 'params' are mutually exclusive.")
+ if i == j:
+ value = self.onsite_hamiltonians[i]
+ if callable(value):
+ if params:
+ required, defaults, takes_kw = self._ham_param_map[value]
+ invalid_params = set(params).intersection(set(defaults))
+ if invalid_params:
+ raise ValueError("Parameters {} have default values "
+ "and may not be set with 'params'"
+ .format(', '.join(invalid_params)))
+ if not takes_kw:
+ params = {pn: params[pn] for pn in required}
+ try:
+ value = value(self.sites[i], **params)
+ except Exception as exc:
+ _raise_user_error(exc, value)
+ else:
+ try:
+ value = value(self.sites[i], *args)
+ except Exception as exc:
+ _raise_user_error(exc, value)
+ else:
+ edge_id = self.graph.first_edge_id(i, j)
+ value = self.hoppings[edge_id]
+ conj = value is Other
+ if conj:
+ i, j = j, i
+ edge_id = self.graph.first_edge_id(i, j)
+ value = self.hoppings[edge_id]
+ if callable(value):
+ sites = self.sites
+ if params:
+ required, defaults, takes_kw = self._ham_param_map[value]
+ invalid_params = set(params).intersection(set(defaults))
+ if invalid_params:
+ raise ValueError("Parameters {} have default values "
+ "and may not be set with 'params'"
+ .format(', '.join(invalid_params)))
+ if not takes_kw:
+ params = {pn: params[pn] for pn in required}
+ try:
+ value = value(sites[i], sites[j], **params)
+ except Exception as exc:
+ _raise_user_error(exc, value)
+ else:
+ try:
+ value = value(sites[i], sites[j], *args)
+ except Exception as exc:
+ _raise_user_error(exc, value)
+ if conj:
+ value = herm_conj(value)
+ return value
+
+ def pos(self, i):
+ return self.sites[i].pos
+
+ discrete_symmetry = discrete_symmetry
+
+
+class InfiniteSystem(system.InfiniteSystem):
+ """Finalized infinite system, extracted from a `Builder`.
- def _finalized_infinite(self, interface_order=None):
+ Attributes
+ ----------
+ sites : sequence
+ ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+ to the integer-labeled site ``i`` of the low-level system.
+ id_by_site : dict
+ The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
+
+ Notes
+ -----
+ In infinite systems ``sites`` consists of 3 parts: sites in the fundamental
+ domain (FD) with hoppings to neighboring cells, sites in the FD with no
+ hoppings to neighboring cells, and sites in FD+1 attached to the FD by
+ hoppings. Each of these three subsequences is individually sorted.
+ """
+
+ def __init__(self, builder, interface_order=None):
"""
- Finalize this builder instance which has to have exactly a single
+ Finalize a builder instance which has to have exactly a single
symmetry direction.
If interface_order is not set, the order of the interface sites in the
finalized system will be arbitrary. If interface_order is set to a
sequence of interface sites, this order will be kept.
"""
- sym = self.symmetry
+ sym = builder.symmetry
assert sym.num_directions == 1
@@ -1684,8 +1865,8 @@ class Builder:
#### neighbors in the previous domain or not.
lsites_with = [] # Fund. domain sites with neighbors in prev. dom
lsites_without = [] # Remaining sites of the fundamental domain
- for tail in self.H: # Loop over all sites of the fund. domain.
- for head in self._out_neighbors(tail):
+ for tail in builder.H: # Loop over all sites of the fund. domain.
+ for head in builder._out_neighbors(tail):
fd = sym.which(head)[0]
if fd == 1:
# Tail belongs to fund. domain, head to the next domain.
@@ -1755,8 +1936,8 @@ class Builder:
g.num_nodes = len(sites) # Some sites could not appear in any edge.
onsite_hamiltonians = []
for tail_id, tail in enumerate(sites[:cell_size]):
- onsite_hamiltonians.append(self.H[tail][1])
- for head in self._out_neighbors(tail):
+ onsite_hamiltonians.append(builder.H[tail][1])
+ for head in builder._out_neighbors(tail):
head_id = id_by_site.get(head)
if head_id is None:
# Head belongs neither to the fundamental domain nor to the
@@ -1785,7 +1966,7 @@ class Builder:
# The tail belongs to the previous domain. Find the
# corresponding hopping with the tail in the fund. domain.
tail, head = sym.to_fd(tail, head)
- hoppings.append(self._get_edge(tail, head))
+ hoppings.append(builder._get_edge(tail, head))
#### Find parameters taken by all value functions
_ham_param_map = {}
@@ -1799,200 +1980,16 @@ class Builder:
params = params[skip:] # remove site argument(s)
_ham_param_map[ham] = (params, defaults, takes_kwargs)
- #### Assemble and return result.
- result = InfiniteSystem()
- result.cell_size = cell_size
- result.sites = sites
- result.id_by_site = id_by_site
- result.site_ranges = _site_ranges(sites)
- result.graph = g
- result.hoppings = hoppings
- result.onsite_hamiltonians = onsite_hamiltonians
- result._ham_param_map = _ham_param_map
- result.symmetry = self.symmetry
- return result
-
-
-################ Finalized systems
-
-def _raise_user_error(exc, func):
- msg = ('Error occurred in user-supplied value function "{0}".\n'
- 'See the upper part of the above backtrace for more information.')
- raise UserCodeError(msg.format(func.__name__)) from exc
-
-
-def discrete_symmetry(self, args=(), *, params=None):
- if self._cons_law is not None:
- eigvals, eigvecs = self._cons_law
- eigvals = eigvals.tocoo(args, params=params)
- if not np.allclose(eigvals.data, np.round(eigvals.data)):
- raise ValueError("Conservation law must have integer eigenvalues.")
- eigvals = np.round(eigvals).tocsr()
- # Avoid appearance of zero eigenvalues
- eigvals = eigvals + 0.5 * sparse.identity(eigvals.shape[0])
- eigvals.eliminate_zeros()
- eigvecs = eigvecs.tocoo(args, params=params)
- projectors = [eigvecs.dot(eigvals == val)
- for val in sorted(np.unique(eigvals.data))]
- else:
- projectors = None
-
- def evaluate(op):
- return None if op is None else op.tocoo(args, params=params)
-
- return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
- self._symmetries))
-
-
-def _transfer_symmetry(syst, builder):
- """Take a symmetry from builder and transfer it to finalized system."""
- def operator(op):
- if op is None:
- return
- return Density(syst, op, check_hermiticity=False)
-
- # Conservation law requires preprocessing to split it into eigenvectors
- # and eigenvalues.
- cons_law = builder.conservation_law
-
- if cons_law is None:
- syst._cons_law = None
-
- elif callable(cons_law):
- @wraps(cons_law)
- def vals(site, *args, **kwargs):
- if site.family.norbs == 1:
- return cons_law(site, *args, **kwargs)
- return np.diag(np.linalg.eigvalsh(cons_law(site, *args, **kwargs)))
-
- @wraps(cons_law)
- def vecs(site, *args, **kwargs):
- if site.family.norbs == 1:
- return 1
- return np.linalg.eigh(cons_law(site, *args, **kwargs))[1]
-
- syst._cons_law = operator(vals), operator(vecs)
-
- elif isinstance(cons_law, collections.Mapping):
- vals = {family: (value if family.norbs == 1 else
- ta.array(np.diag(np.linalg.eigvalsh(value))))
- for family, value in cons_law.items()}
- vecs = {family: (1 if family.norbs == 1 else
- ta.array(np.linalg.eigh(value)[1]))
- for family, value in cons_law.items()}
- syst._cons_law = operator(vals), operator(vecs)
-
- else:
- try:
- vals, vecs = np.linalg.eigh(cons_law)
- vals = np.diag(vals)
- except np.linalg.LinAlgError as e:
- if '0-dimensional' not in e.args[0]:
- raise e # skip coverage
- vals, vecs = cons_law, 1
-
- syst._cons_law = operator(vals), operator(vecs)
-
- syst._symmetries = [operator(symm) for symm in (builder.time_reversal,
- builder.particle_hole,
- builder.chiral)]
-
-
-class FiniteSystem(system.FiniteSystem):
- """Finalized `Builder` with leads.
-
- Usable as input for the solvers in `kwant.solvers`.
-
- Attributes
- ----------
- sites : sequence
- ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
- to the integer-labeled site ``i`` of the low-level system. The sites
- are ordered first by their family and then by their tag.
- id_by_site : dict
- The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
- """
-
- def hamiltonian(self, i, j, *args, params=None):
- if args and params:
- raise TypeError("'args' and 'params' are mutually exclusive.")
- if i == j:
- value = self.onsite_hamiltonians[i]
- if callable(value):
- if params:
- required, defaults, takes_kw = self._ham_param_map[value]
- invalid_params = set(params).intersection(set(defaults))
- if invalid_params:
- raise ValueError("Parameters {} have default values "
- "and may not be set with 'params'"
- .format(', '.join(invalid_params)))
- if not takes_kw:
- params = {pn: params[pn] for pn in required}
- try:
- value = value(self.sites[i], **params)
- except Exception as exc:
- _raise_user_error(exc, value)
- else:
- try:
- value = value(self.sites[i], *args)
- except Exception as exc:
- _raise_user_error(exc, value)
- else:
- edge_id = self.graph.first_edge_id(i, j)
- value = self.hoppings[edge_id]
- conj = value is Other
- if conj:
- i, j = j, i
- edge_id = self.graph.first_edge_id(i, j)
- value = self.hoppings[edge_id]
- if callable(value):
- sites = self.sites
- if params:
- required, defaults, takes_kw = self._ham_param_map[value]
- invalid_params = set(params).intersection(set(defaults))
- if invalid_params:
- raise ValueError("Parameters {} have default values "
- "and may not be set with 'params'"
- .format(', '.join(invalid_params)))
- if not takes_kw:
- params = {pn: params[pn] for pn in required}
- try:
- value = value(sites[i], sites[j], **params)
- except Exception as exc:
- _raise_user_error(exc, value)
- else:
- try:
- value = value(sites[i], sites[j], *args)
- except Exception as exc:
- _raise_user_error(exc, value)
- if conj:
- value = herm_conj(value)
- return value
-
- def pos(self, i):
- return self.sites[i].pos
-
- discrete_symmetry = discrete_symmetry
-
-
-class InfiniteSystem(system.InfiniteSystem):
- """Finalized infinite system, extracted from a `Builder`.
+ self.cell_size = cell_size
+ self.sites = sites
+ self.id_by_site = id_by_site
+ self.site_ranges = _site_ranges(sites)
+ self.graph = g
+ self.hoppings = hoppings
+ self.onsite_hamiltonians = onsite_hamiltonians
+ self._ham_param_map = _ham_param_map
+ self.symmetry = builder.symmetry
- Attributes
- ----------
- sites : sequence
- ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
- to the integer-labeled site ``i`` of the low-level system.
- id_by_site : dict
- The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
-
- Notes
- -----
- In infinite systems ``sites`` consists of 3 parts: sites in the fundamental
- domain (FD) with hoppings to neighboring cells, sites in the FD with no
- hoppings to neighboring cells, and sites in FD+1 attached to the FD by
- hoppings. Each of these three subsequences is individually sorted.
- """
def hamiltonian(self, i, j, *args, params=None):
if args and params:
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 7c290fb5..dee7ced6 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -432,11 +432,11 @@ def test_finalization():
neighbors)
# test that we cannot finalize a system with a badly sorted interface order
- raises(ValueError, lead._finalized_infinite,
+ raises(ValueError, builder.InfiniteSystem, lead,
[builder.Site(fam, n) for n in reversed(neighbors)])
# site ordering independent of whether interface was specified
- flead_order = lead._finalized_infinite([builder.Site(fam, n)
- for n in neighbors])
+ flead_order = builder.InfiniteSystem(lead, [builder.Site(fam, n)
+ for n in neighbors])
assert flead.sites == flead_order.sites
--
GitLab