From cc750c63b9b816130206030fa827cf869795842a Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Tue, 26 Nov 2024 11:19:45 +0100 Subject: [PATCH] added more tests for dtype conversion fixed errors in col creation of dtypes when calling csr.diags() Ensured that sparsegeometry.finalize accepts any arguments the csr.finalize accepts. Removed casting errors in mat2dtype. This may *hide* potential problems when there are non-zero imaginary parts. We should perhaps later revisit the problem considering that TB + Peierls has complex overlap matrices. Fixed issue when a construct was called with Python intrinsic complex values. Signed-off-by: Nick Papior --- src/sisl/_core/sparse.py | 4 +- src/sisl/_core/sparse_geometry.py | 8 +- src/sisl/io/siesta/_help.py | 10 +- src/sisl/io/siesta/binaries.py | 6 +- src/sisl/io/siesta/siesta_nc.py | 4 +- src/sisl/io/siesta/tests/test_matrices.py | 110 ++++++++++++++++++++ src/sisl/io/siesta/tests/test_siesta.py | 121 ---------------------- src/sisl/io/siesta/tests/test_tshs.py | 28 +---- src/sisl/physics/sparse.py | 15 ++- 9 files changed, 141 insertions(+), 165 deletions(-) create mode 100644 src/sisl/io/siesta/tests/test_matrices.py diff --git a/src/sisl/_core/sparse.py b/src/sisl/_core/sparse.py index aa75fdc39b..e8cabbb0a4 100644 --- a/src/sisl/_core/sparse.py +++ b/src/sisl/_core/sparse.py @@ -404,7 +404,7 @@ def diags(self, diagonals, offsets=0, dim=None, dtype=None): shape[2] = dim shape = tuple(shape) - offsets = array_fill_repeat(offsets, shape[0], cls=dtype) + offsets = array_fill_repeat(offsets, shape[0], cls=np.int32) # Create the index-pointer, data and values data = array_fill_repeat(diagonals, shape[0], axis=0, cls=dtype) @@ -488,7 +488,7 @@ def finalized(self): """Whether the contained data is finalized and non-used elements have been removed""" return self._finalized - def finalize(self, sort=True): + def finalize(self, sort: bool = True): """Finalizes the sparse matrix by removing all non-set elements One may still interact with the sparse matrix as one would previously. diff --git a/src/sisl/_core/sparse_geometry.py b/src/sisl/_core/sparse_geometry.py index 0adc7035cb..45c96f081f 100644 --- a/src/sisl/_core/sparse_geometry.py +++ b/src/sisl/_core/sparse_geometry.py @@ -1003,15 +1003,17 @@ def unrepeat( atoms = np.arange(self.geometry.na).reshape(-1, reps).T.ravel() return self.sub(atoms).untile(reps, axis, segment, *args, sym=sym, **kwargs) - def finalize(self): + def finalize(self, *args, **kwargs): """Finalizes the model Finalizes the model so that all non-used elements are removed. I.e. this simply reduces the memory requirement for the sparse matrix. - Note that adding more elements to the sparse matrix is more time-consuming than for a non-finalized sparse matrix due to the + Notes + ----- + Adding more elements to the sparse matrix is more time-consuming than for a non-finalized sparse matrix due to the internal data-representation. """ - self._csr.finalize() + self._csr.finalize(*args, **kwargs) def tocsr(self, dim: int = 0, isc=None, **kwargs): """Return a :class:`~scipy.sparse.csr_matrix` for the specified dimension diff --git a/src/sisl/io/siesta/_help.py b/src/sisl/io/siesta/_help.py index b356f85ec3..d7ac67585d 100644 --- a/src/sisl/io/siesta/_help.py +++ b/src/sisl/io/siesta/_help.py @@ -127,7 +127,8 @@ def toc(D, re, im): csr._D = csr._D.astype(dtype) elif spin.is_noncolinear: D = np.empty(shape[:-1] + (shape[-1] - 1,), dtype=dtype) - D[..., [0, 1]] = csr._D[..., [0, 1]].astype(dtype) + # These should be real only anyways! + D[..., [0, 1]] = csr._D[..., [0, 1]].real.astype(dtype) D[..., 2] = toc(csr._D, 2, 3) if D.shape[-1] > 4: D[..., 3:] = csr._D[..., 4:].astype(dtype) @@ -158,11 +159,12 @@ def toc(D, re, im): csr._D = csr._D.astype(dtype) elif spin.is_noncolinear: D = np.empty(shape[:-1] + (shape[-1] + 1,), dtype=dtype) - D[..., [0, 1]] = csr._D[..., [0, 1]].astype(dtype) + # These should be real only anyways! + D[..., [0, 1]] = csr._D[..., [0, 1]].real.astype(dtype) D[..., 2] = csr._D[..., 2].real.astype(dtype) D[..., 3] = csr._D[..., 2].imag.astype(dtype) if D.shape[-1] > 4: - D[..., 4:] = csr._D[..., 3:].astype(dtype) + D[..., 4:] = csr._D[..., 3:].real.astype(dtype) csr._D = D elif spin.is_spinorbit: D = np.empty(shape[:-1] + (shape[-1] + 4,), dtype=dtype) @@ -175,7 +177,7 @@ def toc(D, re, im): D[..., 6] = csr._D[..., 3].real.astype(dtype) D[..., 7] = csr._D[..., 3].imag.astype(dtype) if D.shape[-1] > 8: - D[..., 8:] = csr._D[..., 4:].astype(dtype) + D[..., 8:] = csr._D[..., 4:].real.astype(dtype) csr._D = D else: raise NotImplementedError diff --git a/src/sisl/io/siesta/binaries.py b/src/sisl/io/siesta/binaries.py index ef8336c0db..fd7224b038 100644 --- a/src/sisl/io/siesta/binaries.py +++ b/src/sisl/io/siesta/binaries.py @@ -721,8 +721,9 @@ def write_density_matrices(self, DM, EDM, Ef: float = 0.0, **kwargs): Ef : fermi-level to be contained """ - DM = DM.transpose(spin=False, sort=False) - EDM = EDM.transpose(spin=False, sort=False) + sort = kwargs.get("sort", True) + DM = DM.transpose(spin=False, sort=sort) + EDM = EDM.transpose(spin=False, sort=sort) DM._csr.align(EDM._csr) EDM._csr.align(DM._csr) @@ -734,7 +735,6 @@ def write_density_matrices(self, DM, EDM, Ef: float = 0.0, **kwargs): _csr_to_siesta(DM.geometry, DM._csr) _csr_to_siesta(DM.geometry, EDM._csr) - sort = kwargs.get("sort", True) DM._csr.finalize(sort=sort) EDM._csr.finalize(sort=sort) _mat_sisl2siesta(DM, dtype=np.float64) diff --git a/src/sisl/io/siesta/siesta_nc.py b/src/sisl/io/siesta/siesta_nc.py index 49a3f4b38b..16d795e53d 100644 --- a/src/sisl/io/siesta/siesta_nc.py +++ b/src/sisl/io/siesta/siesta_nc.py @@ -253,8 +253,8 @@ def read_hamiltonian(self, **kwargs) -> Hamiltonian: _mat_siesta2sisl(H, dtype=kwargs.get("dtype")) # Shift to the Fermi-level - Ef = -self._value("Ef")[:] * Ry2eV - H.shift(Ef) + Ef = self._value("Ef")[:] * Ry2eV + H.shift(-Ef) return H.transpose(spin=False, sort=kwargs.get("sort", True)) diff --git a/src/sisl/io/siesta/tests/test_matrices.py b/src/sisl/io/siesta/tests/test_matrices.py new file mode 100644 index 0000000000..447351055b --- /dev/null +++ b/src/sisl/io/siesta/tests/test_matrices.py @@ -0,0 +1,110 @@ +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at https://mozilla.org/MPL/2.0/. +from __future__ import annotations + +import numpy as np +import pytest + +import sisl +from sisl.io.siesta._help import _mat2dtype + +pytestmark = [pytest.mark.io, pytest.mark.siesta] + +listify = sisl.utils.listify + + +@pytest.mark.parametrize("sort", [True, False]) +@pytest.mark.parametrize( + "matrix,ext", + (map(lambda x: ("Hamiltonian", x), ["nc", "TSHS"]) | listify) + + (map(lambda x: ("DensityMatrix", x), ["nc", "DM"]) | listify) + + (map(lambda x: ("EnergyDensityMatrix", x), ["nc"]) | listify), +) +@pytest.mark.parametrize("read_dtype", [np.float64, np.complex128]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex128]) +def test_non_colinear(sisl_tmp, sort, matrix, ext, dtype, read_dtype): + M = getattr(sisl, matrix)(sisl.geom.graphene(), spin=sisl.Spin("NC"), dtype=dtype) + if np.issubdtype(dtype, np.complexfloating): + onsite = [0.1 + 0j, 0.2 + 0j, 0.3 + 0.4j] + nn = [0.2, 0.3, 0.4 + 0.5j] + else: + onsite = [0.1, 0.2, 0.3, 0.4] + nn = [0.2, 0.3, 0.4, 0.5] + M.construct(([0.1, 1.44], [onsite, nn])) + + f1 = sisl_tmp(f"M1.{ext}") + f2 = sisl_tmp(f"M2.{ext}") + M.write(f1, sort=sort) + M.finalize() + with sisl.get_sile(f1) as sile: + M2 = M.read(sile, dtype=read_dtype) + M2.write(f2, sort=sort) + with sisl.get_sile(f2) as sile: + M3 = M2.read(sile, dtype=read_dtype) + + if sort: + M.finalize(sort=sort) + assert M._csr.spsame(M2._csr) + assert M._csr.spsame(M3._csr) + + from sisl.io.siesta._help import _mat2dtype + + # Convert to the same dtype + _mat2dtype(M2, dtype) + _mat2dtype(M3, dtype) + if M.orthogonal and not M2.orthogonal: + assert np.allclose(M._csr._D, M2._csr._D[..., :-1]) + else: + assert np.allclose(M._csr._D, M2._csr._D) + if M.orthogonal and not M3.orthogonal: + assert np.allclose(M._csr._D, M3._csr._D[..., :-1]) + else: + assert np.allclose(M._csr._D, M3._csr._D) + + +@pytest.mark.parametrize("sort", [True, False]) +@pytest.mark.parametrize( + "matrix,ext", + (map(lambda x: ("Hamiltonian", x), ["nc", "TSHS"]) | listify) + + (map(lambda x: ("DensityMatrix", x), ["nc", "DM"]) | listify) + + (map(lambda x: ("EnergyDensityMatrix", x), ["nc"]) | listify), +) +@pytest.mark.parametrize("read_dtype", [np.float64, np.complex128]) +@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex128]) +def test_spin_orbit(sisl_tmp, sort, matrix, ext, dtype, read_dtype): + M = getattr(sisl, matrix)(sisl.geom.graphene(), spin=sisl.Spin("SO"), dtype=dtype) + if np.issubdtype(dtype, np.complexfloating): + onsite = [0.1 + 0j, 0.2 + 0j, 0.3 + 0.4j, 0.3 - 0.4j] + nn = [0.2 + 0.1j, 0.3 + 0.3j, 0.4 + 0.5j, 0.4 - 0.5j] + else: + onsite = [0.1, 0.2, 0.3, 0.4, 0, 0, 0.3, -0.4] + nn = [0.2, 0.3, 0.4, 0.5, 0.1, 0.3, 0.4, -0.5] + M.construct(([0.1, 1.44], [onsite, nn])) + + f1 = sisl_tmp(f"M1.{ext}") + f2 = sisl_tmp(f"M2.{ext}") + M.write(f1, sort=sort) + M.finalize() + with sisl.get_sile(f1) as sile: + M2 = M.read(sile, dtype=read_dtype) + M2.write(f2, sort=sort) + with sisl.get_sile(f2) as sile: + M3 = M2.read(sile, dtype=read_dtype) + + if sort: + M.finalize(sort=sort) + assert M._csr.spsame(M2._csr) + assert M._csr.spsame(M3._csr) + + # Convert to the same dtype + _mat2dtype(M2, dtype) + _mat2dtype(M3, dtype) + if M.orthogonal and not M2.orthogonal: + assert np.allclose(M._csr._D, M2._csr._D[..., :-1]) + else: + assert np.allclose(M._csr._D, M2._csr._D) + if M.orthogonal and not M3.orthogonal: + assert np.allclose(M._csr._D, M3._csr._D[..., :-1]) + else: + assert np.allclose(M._csr._D, M3._csr._D) diff --git a/src/sisl/io/siesta/tests/test_siesta.py b/src/sisl/io/siesta/tests/test_siesta.py index 6f62ec52ee..e959fe6ed3 100644 --- a/src/sisl/io/siesta/tests/test_siesta.py +++ b/src/sisl/io/siesta/tests/test_siesta.py @@ -161,99 +161,6 @@ def test_nc_density_matrix(sisl_tmp, sisl_system): assert sisl_system.g.atoms.equal(ndm.atoms, R=False) -def test_nc_H_non_colinear(sisl_tmp): - H1 = Hamiltonian(sisl.geom.graphene(), spin=sisl.Spin("NC")) - H1.construct(([0.1, 1.44], [[0.1, 0.2, 0.3, 0.4], [0.2, 0.3, 0.4, 0.5]])) - - f1 = sisl_tmp("H1.nc") - f2 = sisl_tmp("H2.nc") - H1.write(f1) - H1.finalize() - with sisl.get_sile(f1) as sile: - H2 = sile.read_hamiltonian() - H2.write(f2) - with sisl.get_sile(f2) as sile: - H3 = sile.read_hamiltonian() - assert H1._csr.spsame(H2._csr) - assert np.allclose(H1._csr._D, H2._csr._D) - assert H1._csr.spsame(H3._csr) - assert np.allclose(H1._csr._D, H3._csr._D) - - -def test_nc_DM_non_colinear(sisl_tmp): - DM1 = DensityMatrix(sisl.geom.graphene(), spin=sisl.Spin("NC")) - DM1.construct(([0.1, 1.44], [[0.1, 0.2, 0.3, 0.4], [0.2, 0.3, 0.4, 0.5]])) - - f1 = sisl_tmp("DM1.nc") - f2 = sisl_tmp("DM2.nc") - DM1.write(f1) - DM1.finalize() - with sisl.get_sile(f1) as sile: - DM2 = sile.read_density_matrix() - DM2.write(f2) - with sisl.get_sile(f2) as sile: - DM3 = sile.read_density_matrix() - assert DM1._csr.spsame(DM2._csr) - assert DM1._csr.spsame(DM3._csr) - # DM1 is finalized, but DM2 is not finalized - assert np.allclose(DM1._csr._D, DM2._csr._D) - # DM2 and DM3 are the same - assert np.allclose(DM2._csr._D, DM3._csr._D) - DM2.finalize() - assert np.allclose(DM1._csr._D, DM2._csr._D) - - -def test_nc_EDM_non_colinear(sisl_tmp): - EDM1 = EnergyDensityMatrix(sisl.geom.graphene(), spin=sisl.Spin("NC")) - EDM1.construct(([0.1, 1.44], [[0.1, 0.2, 0.3, 0.4], [0.2, 0.3, 0.4, 0.5]])) - - f1 = sisl_tmp("EDM1.nc") - f2 = sisl_tmp("EDM2.nc") - EDM1.write(f1, sort=False) - EDM1.finalize() - with sisl.get_sile(f1) as sile: - EDM2 = sile.read_energy_density_matrix(sort=False) - EDM2.write(f2, sort=False) - with sisl.get_sile(f2) as sile: - EDM3 = sile.read_energy_density_matrix(sort=False) - assert EDM1._csr.spsame(EDM2._csr) - assert EDM1._csr.spsame(EDM3._csr) - # EDM1 is finalized, but EDM2 is not finalized - assert not np.allclose(EDM1._csr._D, EDM2._csr._D) - # EDM2 and EDM3 are the same - assert np.allclose(EDM2._csr._D, EDM3._csr._D) - EDM2.finalize() - assert np.allclose(EDM1._csr._D, EDM2._csr._D) - - -@pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") -def test_nc_H_spin_orbit(sisl_tmp): - H1 = Hamiltonian(sisl.geom.graphene(), spin=sisl.Spin("SO")) - H1.construct( - ( - [0.1, 1.44], - [ - [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8], - [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], - ], - ) - ) - - f1 = sisl_tmp("H1.nc") - f2 = sisl_tmp("H2.nc") - H1.write(f1) - H1.finalize() - with sisl.get_sile(f1) as sile: - H2 = sile.read_hamiltonian() - H2.write(f2) - with sisl.get_sile(f2) as sile: - H3 = sile.read_hamiltonian() - assert H1._csr.spsame(H2._csr) - assert np.allclose(H1._csr._D, H2._csr._D) - assert H1._csr.spsame(H3._csr) - assert np.allclose(H1._csr._D, H3._csr._D) - - @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") def test_nc_H_spin_orbit_nc2tshs2nc(sisl_tmp): H1 = Hamiltonian(sisl.geom.graphene(), spin=sisl.Spin("SO")) @@ -282,34 +189,6 @@ def test_nc_H_spin_orbit_nc2tshs2nc(sisl_tmp): assert np.allclose(H1._csr._D, H3._csr._D) -@pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") -def test_nc_DM_spin_orbit(sisl_tmp): - DM1 = DensityMatrix(sisl.geom.graphene(), spin=sisl.Spin("SO")) - DM1.construct( - ( - [0.1, 1.44], - [ - [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8], - [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], - ], - ) - ) - - f1 = sisl_tmp("DM1.nc") - f2 = sisl_tmp("DM2.nc") - DM1.write(f1) - DM1.finalize() - with sisl.get_sile(f1) as sile: - DM2 = sile.read_density_matrix() - DM2.write(f2) - with sisl.get_sile(f2) as sile: - DM3 = sile.read_density_matrix() - assert DM1._csr.spsame(DM2._csr) - assert np.allclose(DM1._csr._D, DM2._csr._D) - assert DM1._csr.spsame(DM3._csr) - assert np.allclose(DM1._csr._D, DM3._csr._D) - - @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") def test_nc_DM_spin_orbit_nc2dm2nc(sisl_tmp): DM1 = DensityMatrix(sisl.geom.graphene(), orthogonal=False, spin=sisl.Spin("SO")) diff --git a/src/sisl/io/siesta/tests/test_tshs.py b/src/sisl/io/siesta/tests/test_tshs.py index 9d18d50476..20e6fc4759 100644 --- a/src/sisl/io/siesta/tests/test_tshs.py +++ b/src/sisl/io/siesta/tests/test_tshs.py @@ -29,7 +29,7 @@ def test_tshs_si_pdos_kgrid(sisl_files, sisl_tmp): @pytest.mark.filterwarnings("ignore", message="*Casting complex values") -def test_tshs_si_pdos_dtypes(sisl_files, sisl_tmp): +def test_tshs_si_pdos_dtypes_eigs(sisl_files, sisl_tmp): si = sisl.get_sile(sisl_files("siesta", "Si_pdos_k", "Si_pdos.TSHS")) data = [] eigs = None @@ -193,32 +193,6 @@ def test_tshs_si_pdos_kgrid_overlap(sisl_files): assert np.allclose(HS._csr._D[:, HS.S_idx], S._csr._D[:, 0]) -@pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") -def test_tshs_spin_orbit(sisl_tmp): - H1 = sisl.Hamiltonian(sisl.geom.graphene(), spin=sisl.Spin("SO")) - H1.construct( - ( - [0.1, 1.44], - [ - [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8], - [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], - ], - ) - ) - - f1 = sisl_tmp("tmp1.TSHS") - f2 = sisl_tmp("tmp2.TSHS") - H1.write(f1) - H1.finalize() - H2 = sisl.get_sile(f1).read_hamiltonian() - H2.write(f2) - H3 = sisl.get_sile(f2).read_hamiltonian() - assert H1._csr.spsame(H2._csr) - assert np.allclose(H1._csr._D, H2._csr._D) - assert H1._csr.spsame(H3._csr) - assert np.allclose(H1._csr._D, H3._csr._D) - - @pytest.mark.filterwarnings("ignore", message="*is NOT Hermitian for on-site") def test_tshs_spin_orbit_tshs2nc2tshs(sisl_tmp): pytest.importorskip("netCDF4") diff --git a/src/sisl/physics/sparse.py b/src/sisl/physics/sparse.py index c4be7cf32b..c5e8a79fd4 100644 --- a/src/sisl/physics/sparse.py +++ b/src/sisl/physics/sparse.py @@ -979,7 +979,13 @@ def create_construct(self, R, params): # Hermitian parameters # The input order is [uu, dd, ud, du] paramsH = [ - [p[0].conj(), p[1].conj(), p[3].conj(), p[2].conj(), *p[4:]] + [ + p[0].conjugate(), + p[1].conjugate(), + p[3].conjugate(), + p[2].conjugate(), + *p[4:], + ] for p in params ] else: @@ -1009,7 +1015,7 @@ def create_construct(self, R, params): ], dtype_cplx, ) - if not np.allclose(onsite, onsite.T.conj()): + if not np.allclose(onsite, onsite.T.conjugate()): warn( f"{self.__class__.__name__}.create_construct is NOT " "Hermitian for on-site terms. This is your responsibility! " @@ -1020,7 +1026,10 @@ def create_construct(self, R, params): if is_complex: nv = 3 # Hermitian parameters - paramsH = [[p[0].conj(), p[1].conj(), p[2], *p[3:]] for p in params] + paramsH = [ + [p[0].conjugate(), p[1].conjugate(), p[2], *p[3:]] + for p in params + ] else: nv = 4 # Hermitian parameters