diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e6b6d037..a67247017 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/sisl/physics/_ufuncs_electron.py b/src/sisl/physics/_ufuncs_electron.py index 936bc98e0..5984a95dd 100644 --- a/src/sisl/physics/_ufuncs_electron.py +++ b/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 diff --git a/src/sisl/physics/electron.py b/src/sisl/physics/electron.py index d1ebb018e..caa07ae34 100644 --- a/src/sisl/physics/electron.py +++ b/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 diff --git a/src/sisl/physics/state.py b/src/sisl/physics/state.py index 85111c9e6..6a08e1878 100644 --- a/src/sisl/physics/state.py +++ b/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 diff --git a/src/sisl/physics/tests/test_hamiltonian.py b/src/sisl/physics/tests/test_hamiltonian.py index b22b0dbd2..a77cc67b1 100644 --- a/src/sisl/physics/tests/test_hamiltonian.py +++ b/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]]