removed degenerate argument, re-inserted - sign for SHC

After more careful checking, it corresponds well
to litterature: 10.1103/PhysRevMaterials.6.095001.

Also removed the degenerate argument in derivative.

This should be the users responsibility.

Signed-off-by: Nick Papior <[email protected]>

CHANGELOG.md

@@ -101,6 +101,7 @@ we hit release version 1.0.0.
   be done for assigning matrix elements (it fills with 0's).
 
 ### Removed
+- `degenerate` argument in `velocity`/`derivative`, they do not belong there
 - `xvSileSiesta.read_geometry(species_as_Z)`, deprecated in favor of `atoms=`
 - `structSileSiesta.read_geometry(species_as_Z)`, deprecated in favor of `atoms=`
 - `Atom.radii` is removed, `Atom.radius` is the correct invocation

src/sisl/physics/_ufuncs_electron.py

@@ -32,10 +32,6 @@ def velocity(state: StateCElectron, *args, **kwargs):
 
     Notes
     -----
-    The states and energies for the states *may* have changed after calling this routine.
-    This is because of the velocity un-folding for degenerate modes. I.e. calling
-    `PDOS` after this method *may* change the result.
-
     The velocities are calculated without the Berry curvature contribution see Eq. (2) in :cite:`Wang2006`.
     It is thus typically denoted as the *effective velocity operater* (see Ref. 21 in :cite:`Wang2006`.
     The missing contribution may be added in later editions, for completeness sake, it is:
@@ -45,6 +41,11 @@ def velocity(state: StateCElectron, *args, **kwargs):
 
     where :math:`\Omega_i` is the Berry curvature for state :math:`i`.
 
+    Parameters
+    ----------
+    *args, **kwargs:
+        arguments passed directly to `derivative`, see that method
+        for argument details.
 
     See Also
     --------
@@ -154,11 +155,7 @@ def berry_curvature(
     else:
         # different operators
         dA = state.derivative(order=1, operator=opA, matrix=True, **derivative_kwargs)
-
-        dB_kwargs = {**derivative_kwargs}
-        # we shouldn't change the states *again*
-        dB_kwargs.pop("degenerate", None)
-        dB = state.derivative(order=1, operator=opB, matrix=True, **dB_kwargs)
+        dB = state.derivative(order=1, operator=opB, matrix=True, **derivative_kwargs)
 
     ieta2 = 1j * eta**2
     energy = state.c

src/sisl/physics/electron.py

@@ -973,9 +973,9 @@ def noop(M, d):
     # simply multiply by: -1/2
     shc_idx = [i for i in map(direction, J_axes) if i in axes]
     if k_average:
-        cond[shc_idx] *= 0.5
+        cond[shc_idx] *= -0.5
     else:
-        cond[:, shc_idx] *= 0.5
+        cond[:, shc_idx] *= -0.5
 
     return cond
 

src/sisl/physics/state.py

@@ -61,6 +61,49 @@ def degenerate_decouple(state, M):
     return state
 
 
+"""
+        degenerate : float or list of array_like, optional
+           If a float is passed it is regarded as the degeneracy tolerance used to calculate the degeneracy
+           levels. Defaults to 1e-5 eV.
+           If a list, it contains the indices of degenerate states. In that case a prior diagonalization
+           is required to decouple them. See `degenerate_dir` for the sum of directions.
+        degenerate_dir : (3,), optional
+           a direction used for degenerate decoupling. The decoupling based on the velocity along this direction
+
+        # Now parse the degeneracy handling
+        if degenerate is not None:
+            if isinstance(degenerate, Real):
+                degenerate = self.degenerate(degenerate)
+
+            # normalize direction
+            degenerate_dir = _a.asarrayd(degenerate_dir)
+            degenerate_dir /= (degenerate_dir @ degenerate_dir) ** 0.5
+
+            # de-coupling is only done for the 1st derivative
+
+            # create the degeneracy decoupling projector
+            deg_dPk = sum(d * dh for d, dh in zip(degenerate_dir, dPk))
+
+            if is_orthogonal:
+                for deg in degenerate:
+                    # Set the average energy
+                    energy[deg] = np.average(energy[deg])
+                    # Now diagonalize to find the contributions from individual states
+                    # then re-construct the seperated degenerate states
+                    # Since we do this for all directions we should decouple them all
+                    state[deg] = degenerate_decouple(state[deg], deg_dPk)
+            else:
+                for deg in degenerate:
+                    e = np.average(energy[deg])
+                    energy[deg] = e
+                    deg_dSk = sum((d * e) * ds for d, ds in zip(degenerate_dir, dSk))
+                    state[deg] = degenerate_decouple(state[deg], deg_dPk - deg_dSk)
+                    del deg_dSk
+
+            del deg_dPk
+"""
+
+
 class _FakeMatrix:
     """Replacement object which superseedes a matrix"""
 
@@ -230,12 +273,12 @@ def shape(self):
         "0.15",
         "0.16",
     )
-    def degenerate(self, atol: float = 1e-8):
+    def degenerate(self, atol: float):
         """Find degenerate coefficients with a specified precision
 
         Parameters
         ----------
-        atol : float, optional
+        atol :
            the precision above which coefficients are not considered degenerate
 
         Returns
@@ -402,7 +445,7 @@ def __getitem__(self, key):
         """
         return self.sub(key)
 
-    def iter(self, asarray=False):
+    def iter(self, asarray: bool = False):
         """An iterator looping over the states in this system
 
         Parameters
@@ -1098,8 +1141,6 @@ def sort(self, ascending: bool = True):
     def derivative(
         self,
         order: Literal[1, 2] = 1,
-        degenerate=None,
-        degenerate_dir=(1, 1, 1),
         matrix: bool = False,
         axes: CartesianAxes = "xyz",
         operator: _dM_Operator = lambda M, d=None: M,
@@ -1141,13 +1182,6 @@ def derivative(
         ----------
         order :
            an integer specifying which order of the derivative is being calculated.
-        degenerate : float or list of array_like, optional
-           If a float is passed it is regarded as the degeneracy tolerance used to calculate the degeneracy
-           levels. Defaults to 1e-5 eV.
-           If a list, it contains the indices of degenerate states. In that case a prior diagonalization
-           is required to decouple them. See `degenerate_dir` for the sum of directions.
-        degenerate_dir : (3,), optional
-           a direction used for degenerate decoupling. The decoupling based on the velocity along this direction
         matrix :
            whether the full matrix or only the diagonal components are returned
         axes:
@@ -1297,38 +1331,6 @@ def add_keys(opt, *keys):
         # in case the state is not an Eigenstate*
         energy = self.c
 
-        # Now parse the degeneracy handling
-        if degenerate is not None:
-            if isinstance(degenerate, Real):
-                degenerate = self.degenerate(degenerate)
-
-            # normalize direction
-            degenerate_dir = _a.asarrayd(degenerate_dir)
-            degenerate_dir /= (degenerate_dir @ degenerate_dir) ** 0.5
-
-            # de-coupling is only done for the 1st derivative
-
-            # create the degeneracy decoupling projector
-            deg_dPk = sum(d * dh for d, dh in zip(degenerate_dir, dPk))
-
-            if is_orthogonal:
-                for deg in degenerate:
-                    # Set the average energy
-                    energy[deg] = np.average(energy[deg])
-                    # Now diagonalize to find the contributions from individual states
-                    # then re-construct the seperated degenerate states
-                    # Since we do this for all directions we should decouple them all
-                    state[deg] = degenerate_decouple(state[deg], deg_dPk)
-            else:
-                for deg in degenerate:
-                    e = np.average(energy[deg])
-                    energy[deg] = e
-                    deg_dSk = sum((d * e) * ds for d, ds in zip(degenerate_dir, dSk))
-                    state[deg] = degenerate_decouple(state[deg], deg_dPk - deg_dSk)
-                    del deg_dSk
-
-            del deg_dPk
-
         # States have been decoupled and we can calculate things now
         # number of states
         nstate, lstate = state.shape

src/sisl/physics/tests/test_hamiltonian.py

@@ -756,29 +756,32 @@ def test_gauge_velocity(self, setup):
         k = [0.1] * 3
         # This is a degenerate eigenstate:
         #  2, 2, 4, 4, 2, 2
-        # and hence a decoupling is necessary
+        # It would be nice to check decoupling.
+        # Currently this is disabled, but should be used.
+
+        # Some comments from the older code which enabled automatic
+        # decoupling:
         # This test is the reason why a default degenerate=1e-5 is used
         # since the gauge='cell' yields the correct *decoupled* states
         # where as gauge='orbital' mixes them in a bad way.
         es1 = H.eigenstate(k, gauge="cell")
         es2 = H.eigenstate(k, gauge="orbital")
-        assert np.allclose(es1.velocity(), es2.velocity())
-        assert np.allclose(es1.velocity(), es2.velocity(degenerate_dir=(1, 1, 0)))
+        assert not np.allclose(es1.velocity(), es2.velocity())
 
         es2.change_gauge("cell")
-        assert np.allclose(es1.velocity(), es2.velocity())
+        assert not np.allclose(es1.velocity(), es2.velocity())
 
         es2.change_gauge("orbital")
         es1.change_gauge("orbital")
         v1 = es1.velocity()
         v2 = es2.velocity()
-        assert np.allclose(v1, v2)
+        assert not np.allclose(v1, v2)
 
         # Projected velocity
         vv1 = es1.velocity(matrix=True)
         vv2 = es2.velocity(matrix=True)
-        assert np.allclose(np.diagonal(vv1, axis1=1, axis2=2), v2)
-        assert np.allclose(np.diagonal(vv2, axis1=1, axis2=2), v1)
+        assert not np.allclose(np.diagonal(vv1, axis1=1, axis2=2), v2)
+        assert not np.allclose(np.diagonal(vv2, axis1=1, axis2=2), v1)
 
     def test_derivative_orthogonal(self, setup):
         R, param = [0.1, 1.5], [1.0, 0.1]
@@ -971,7 +974,7 @@ def test_berry_curvature(self, setup):
         k = [0.1] * 3
         ie1 = H.eigenstate(k, gauge="cell").berry_curvature()
         ie2 = H.eigenstate(k, gauge="orbital").berry_curvature()
-        assert np.allclose(ie1, ie2)
+        assert not np.allclose(ie1, ie2)
 
     def test_spin_berry_curvature(self, setup):
         R, param = [0.1, 1.5], [[1.0, 0.1, 0, 0], [0.4, 0.2, 0.3, 0.2]]