Skip to content

Commit

Permalink
added more tests for dtype conversion
Browse files Browse the repository at this point in the history
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 <[email protected]>
  • Loading branch information
zerothi committed Nov 26, 2024
1 parent f6df3fb commit cc750c6
Show file tree
Hide file tree
Showing 9 changed files with 141 additions and 165 deletions.
4 changes: 2 additions & 2 deletions src/sisl/_core/sparse.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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.
Expand Down
8 changes: 5 additions & 3 deletions src/sisl/_core/sparse_geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 6 additions & 4 deletions src/sisl/io/siesta/_help.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
6 changes: 3 additions & 3 deletions src/sisl/io/siesta/binaries.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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)
Expand Down
4 changes: 2 additions & 2 deletions src/sisl/io/siesta/siesta_nc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Expand Down
110 changes: 110 additions & 0 deletions src/sisl/io/siesta/tests/test_matrices.py
Original file line number Diff line number Diff line change
@@ -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)
121 changes: 0 additions & 121 deletions src/sisl/io/siesta/tests/test_siesta.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand Down Expand Up @@ -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"))
Expand Down
28 changes: 1 addition & 27 deletions src/sisl/io/siesta/tests/test_tshs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand Down
Loading

0 comments on commit cc750c6

Please sign in to comment.