Merge branch 'issue_236' into 'master'

allow computations even if a discrete symmetry is broken

Closes #236 and #242

See merge request kwant/kwant!257

kwant/physics/symmetry.py

@@ -141,10 +141,11 @@ class DiscreteSymmetry:
 
         Returns
         -------
-        broken_symmetry : string or ``None``
-            One of "Conservation law", "Time reversal", "Particle-hole",
-            "Chiral": the symmetry broken by the matrix. If the matrix breaks
-            more than one symmetry, returns only the first failed check.
+        broken_symmetries : list
+            List of strings, the names of symmetries broken by the
+            matrix: any combination of "Conservation law", "Time reversal",
+            "Particle-hole", "Chiral". If no symmetries are broken, returns
+            an empty list.
         """
         # Extra transposes are to enforse sparse dot product in case matrix is
         # dense.
@@ -157,19 +158,22 @@ class DiscreteSymmetry:
                 matrix = new_matrix
             else:
                 matrix = hstack([matrix, csr_matrix((n, n-m))])
+        broken_symmetries = []
         if self.projectors is not None:
             for proj in self.projectors:
                 full = proj.dot(proj.T.conj())
                 commutator = full.dot(matrix) - (full.T.dot(matrix.T)).T
                 if np.linalg.norm(commutator.data) > 1e-8:
-                    return 'Conservation law'
+                    broken_symmetries.append('Conservation law')
+                    break
         for symm, conj, sign, name in zip(self[1:], _conj, _signs, _names):
             if symm is None:
                 continue
             commutator = symm.T.conj().dot((symm.T.dot(matrix.T)).T)
             commutator = commutator - sign * cond_conj(matrix, conj)
             if np.linalg.norm(commutator.data) > 1e-8:
-                return name
+                broken_symmetries.append(name)
+        return broken_symmetries
 
     def __getitem__(self, item):
         return (self.projectors, self.time_reversal,

kwant/physics/tests/test_symmetry.py

@@ -151,22 +151,22 @@ def test_validate():
     csr = sparse.csr_matrix
     sym = DiscreteSymmetry(projectors=[csr(np.array([[1], [0]])),
                                        csr(np.array([[0], [1]]))])
-    assert sym.validate(csr(np.array([[0], [1]]))) == 'Conservation law'
-    assert sym.validate(np.array([[1], [0]])) is None
-    assert sym.validate(np.eye(2)) is None
-    assert sym.validate(1 - np.eye(2)) == 'Conservation law'
+    assert sym.validate(csr(np.array([[0], [1]]))) == ['Conservation law']
+    assert sym.validate(np.array([[1], [0]])) == []
+    assert sym.validate(np.eye(2)) == []
+    assert sym.validate(1 - np.eye(2)) == ['Conservation law']
 
     sym = DiscreteSymmetry(particle_hole=sparse.identity(2))
-    assert sym.validate(1j * sparse.identity(2)) is None
-    assert sym.validate(sparse.identity(2)) == 'Particle-hole'
+    assert sym.validate(1j * sparse.identity(2)) == []
+    assert sym.validate(sparse.identity(2)) == ['Particle-hole']
 
     sym = DiscreteSymmetry(time_reversal=sparse.identity(2))
-    assert sym.validate(sparse.identity(2)) is None
-    assert sym.validate(1j * sparse.identity(2)) == 'Time reversal'
+    assert sym.validate(sparse.identity(2)) == []
+    assert sym.validate(1j * sparse.identity(2)) == ['Time reversal']
 
     sym = DiscreteSymmetry(chiral=csr(np.diag((1, -1))))
-    assert sym.validate(np.eye(2)) == 'Chiral'
-    assert sym.validate(1 - np.eye(2)) is None
+    assert sym.validate(np.eye(2)) == ['Chiral']
+    assert sym.validate(1 - np.eye(2)) == []
 
 
 def random_onsite_hop(n, rng=0):
@@ -178,7 +178,13 @@ def random_onsite_hop(n, rng=0):
 
 
 def test_validate_commutator():
-    symm_class = ['AI', 'AII', 'D', 'C', 'AIII']
+    symm_class = ['AI', 'AII', 'D', 'C', 'AIII', 'BDI']
+    sym_dict = {'AI': ['Time reversal'],
+                'AII': ['Time reversal'],
+                'D': ['Particle-hole'],
+                'C': ['Particle-hole'],
+                'AIII': ['Chiral'],
+                'BDI': ['Time reversal', 'Particle-hole', 'Chiral']}
     n = 10
     rng = 10
     for sym in symm_class:
@@ -201,6 +207,7 @@ def test_validate_commutator():
         disc_symm = DiscreteSymmetry(particle_hole=p_mat,
                                      time_reversal=t_mat,
                                      chiral=c_mat)
-        assert disc_symm.validate(h) == None
+        assert disc_symm.validate(h) == []
         a = random_onsite_hop(n, rng=rng)[1]
-        assert disc_symm.validate(a) in ['Time reversal', 'Particle-hole', 'Chiral']
+        for symmetry in disc_symm.validate(a):
+            assert symmetry in sym_dict[sym]

kwant/system.py

@@ -11,6 +11,7 @@
 __all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
 
 import abc
+import warnings
 from copy import copy
 from . import _system
 
@@ -163,6 +164,17 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
         result.leads = new_leads
         return result
 
+    def validate_symmetries(self, args=(), *, params=None):
+        """Check that the Hamiltonian satisfies discrete symmetries.
+
+        Applies `~kwant.physics.DiscreteSymmetry.validate` to the
+        Hamiltonian, see its documentation for details on the return
+        format.
+        """
+        symmetries = self.discrete_symmetry(args=args, params=params)
+        ham = self.hamiltonian_submatrix(args, sparse=True, params=params)
+        return symmetries.validate(ham)
+
 
 class InfiniteSystem(System, metaclass=abc.ABCMeta):
     """Abstract infinite low-level system.
@@ -235,12 +247,19 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
         ham = self.cell_hamiltonian(args, params=params)
         hop = self.inter_cell_hopping(args, params=params)
         symmetries = self.discrete_symmetry(args, params=params)
-        broken = symmetries.validate(ham)
-        if broken is not None:
-            raise ValueError("Cell Hamiltonian breaks " + broken.lower())
-        broken = symmetries.validate(hop)
-        if broken is not None:
-            raise ValueError("Inter-cell hopping breaks " + broken.lower())
+        # Check whether each symmetry is broken.
+        # If a symmetry is broken, it is ignored in the computation.
+        broken = set(symmetries.validate(ham) + symmetries.validate(hop))
+        attribute_names = {'Conservation law': 'projectors',
+                          'Time reversal': 'time_reversal',
+                          'Particle-hole': 'particle-hole',
+                          'Chiral': 'chiral'}
+        for name in broken:
+            warnings.warn('Hamiltonian breaks ' + name +
+                          ', ignoring the symmetry in the computation.')
+            assert name in attribute_names, 'Inconsistent naming of symmetries'
+            setattr(symmetries, attribute_names[name], None)
+
         shape = ham.shape
         assert len(shape) == 2
         assert shape[0] == shape[1]
@@ -269,6 +288,19 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
         return physics.selfenergy(ham,
                                   self.inter_cell_hopping(args, params=params))
 
+    def validate_symmetries(self, args=(), *, params=None):
+        """Check that the Hamiltonian satisfies discrete symmetries.
+
+        Returns `~kwant.physics.DiscreteSymmetry.validate` applied
+        to the onsite matrix and the hopping. See its documentation for
+        details on the return format.
+        """
+        symmetries = self.discrete_symmetry(args=args, params=params)
+        ham = self.cell_hamiltonian(args=args, sparse=True, params=params)
+        hop = self.inter_cell_hopping(args=args, sparse=True, params=params)
+        broken = set(symmetries.validate(ham) + symmetries.validate(hop))
+        return list(broken)
+
 
 class PrecalculatedLead:
     def __init__(self, modes=None, selfenergy=None):