Skip to content

Commit 43fc1b9

Browse files
zerothizerothi
committed
fixed precision checks for dtype
1 parent 6534168 commit 43fc1b9

File tree

2 files changed

+45
-36
lines changed

2 files changed

+45
-36
lines changed

‎src/sisl/physics/hamiltonian.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -367,7 +367,7 @@ def eigenvalue(self, k: KPoint = (0, 0, 0), gauge: GaugeType = "lattice", **kwar
367367
if kwargs.pop("sparse", False):
368368
e = self.eigsh(k, gauge=gauge, eigvals_only=True, **kwargs)
369369
else:
370-
e = self.eigh(k, gauge, eigvals_only=True, **kwargs)
370+
e = self.eigh(k, gauge=gauge, eigvals_only=True, **kwargs)
371371
info = {"k": k, "gauge": gauge}
372372
for name in ("spin",):
373373
if name in kwargs:
@@ -408,7 +408,7 @@ def eigenstate(self, k: KPoint = (0, 0, 0), gauge: GaugeType = "lattice", **kwar
408408
if kwargs.pop("sparse", False):
409409
e, v = self.eigsh(k, gauge=gauge, eigvals_only=False, **kwargs)
410410
else:
411-
e, v = self.eigh(k, gauge, eigvals_only=False, **kwargs)
411+
e, v = self.eigh(k, gauge=gauge, eigvals_only=False, **kwargs)
412412
info = {"k": k, "gauge": gauge}
413413
for name in ("spin",):
414414
if name in kwargs:

‎src/sisl/physics/tests/test_hamiltonian.py‎

Lines changed: 43 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313
else:
1414
from numpy import ComplexWarning
1515
import pytest
16-
from scipy.linalg import block_diag
16+
from scipy.linalg import block_diag, eigh
1717
from scipy.sparse import SparseEfficiencyWarning, issparse
1818

1919
from sisl import (
@@ -1584,6 +1584,7 @@ def test_spin_contamination(self, setup, k):
15841584
def test_non_colinear_orthogonal(self, setup, sisl_tolerance):
15851585
atol, rtol = sisl_tolerance[np.complex64]
15861586
allclose = partial(np.allclose, atol=atol, rtol=rtol)
1587+
eigh_kwargs = dict(driver="evd")
15871588

15881589
g = Geometry(
15891590
[[i, 0, 0] for i in range(10)],
@@ -1604,9 +1605,9 @@ def test_non_colinear_orthogonal(self, setup, sisl_tolerance):
16041605
H[i, i + 1, 0] = 1.0
16051606
H[i, i + 1, 1] = 1.0
16061607

1607-
eig1 = H.eigh(dtype=np.complex64)
1608-
assert allclose(H.eigh(dtype=np.complex128), eig1)
1609-
assert allclose(H.eigh(gauge="atom", dtype=np.complex128), eig1)
1608+
eig1 = H.eigh(dtype=np.complex64, **eigh_kwargs)
1609+
assert allclose(H.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
1610+
assert allclose(H.eigh(gauge="atom", dtype=np.complex128, **eigh_kwargs), eig1)
16101611
assert len(eig1) == len(H)
16111612

16121613
H1 = Hamiltonian(g, dtype=np.float64, spin=Spin("non-collinear"))
@@ -1624,38 +1625,38 @@ def test_non_colinear_orthogonal(self, setup, sisl_tolerance):
16241625
H1[i, i + 1, 1] = 1.0
16251626
assert H1.spsame(H)
16261627

1627-
eig1 = H1.eigh(dtype=np.complex64)
1628-
assert allclose(H1.eigh(dtype=np.complex128), eig1)
1629-
assert np.allclose(H.eigh(), H1.eigh())
1628+
eig1 = H1.eigh(dtype=np.complex64, **eigh_kwargs)
1629+
assert allclose(H1.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
1630+
assert allclose(H.eigh(**eigh_kwargs), H1.eigh(**eigh_kwargs))
16301631

16311632
# Create the block matrix for expectation
16321633
SZ = block_diag(*([H1.spin.Z] * H1.no))
16331634

16341635
for dtype in (np.complex64, np.complex128):
1635-
es = H1.eigenstate(dtype=dtype)
1636+
es = H1.eigenstate(dtype=dtype, **eigh_kwargs)
16361637
assert allclose(es.eig, eig1)
1637-
assert np.allclose(es.inner(), 1)
1638+
assert allclose(es.inner(), 1)
16381639

16391640
# Perform spin-moment calculation
16401641
sm = es.spin_moment()
16411642
sm2 = es.inner(matrix=SZ).real
16421643
sm3 = np.diag(np.dot(np.conj(es.state), SZ).dot(es.state.T)).real
1643-
assert np.allclose(sm[2], sm2)
1644-
assert np.allclose(sm[2], sm3)
1644+
assert allclose(sm[2], sm2)
1645+
assert allclose(sm[2], sm3)
16451646

16461647
om = es.spin_moment(projection=True)
1647-
assert np.allclose(sm, om.sum(-1))
1648+
assert allclose(sm, om.sum(-1))
16481649

16491650
PDOS = es.PDOS(np.linspace(-1, 1, 21))
16501651
DOS = es.DOS(np.linspace(-1, 1, 21))
1651-
assert np.allclose(PDOS.sum(1)[0, :], DOS)
1652+
assert allclose(PDOS.sum(1)[0, :], DOS, atol=1e-4)
16521653
es.velocity(matrix=True)
16531654

16541655
# Check the velocities
16551656
# But only compare for np.float64, we need the precision
16561657
v = es.velocity()
16571658
vv = es.velocity(matrix=True)
1658-
assert np.allclose(np.diagonal(vv).T, v)
1659+
assert allclose(np.diagonal(vv).T, v)
16591660

16601661
# Ensure we can change gauge for NC stuff
16611662
es.change_gauge("cell")
@@ -1664,6 +1665,7 @@ def test_non_colinear_orthogonal(self, setup, sisl_tolerance):
16641665
def test_non_colinear_non_orthogonal(self, sisl_tolerance):
16651666
atol, rtol = sisl_tolerance[np.complex64]
16661667
allclose = partial(np.allclose, atol=atol * 10, rtol=rtol * 100)
1668+
eigh_kwargs = dict(driver="gvd")
16671669

16681670
g = Geometry(
16691671
[[i, 0, 0] for i in range(10)],
@@ -1685,8 +1687,8 @@ def test_non_colinear_non_orthogonal(self, sisl_tolerance):
16851687
H[i, i + 1, 1] = 1.0
16861688
H.S[i, i] = 1.0
16871689

1688-
eig1 = H.eigh(dtype=np.complex64)
1689-
assert allclose(H.eigh(dtype=np.complex128), eig1)
1690+
eig1 = H.eigh(dtype=np.complex64, **eigh_kwargs)
1691+
assert allclose(H.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
16901692
assert len(eig1) == len(H)
16911693

16921694
H1 = Hamiltonian(
@@ -1707,29 +1709,32 @@ def test_non_colinear_non_orthogonal(self, sisl_tolerance):
17071709
H1.S[i, i] = 1.0
17081710
assert H1.spsame(H)
17091711

1710-
eig1 = H1.eigh(dtype=np.complex64)
1711-
assert allclose(H1.eigh(dtype=np.complex128), eig1)
1712-
assert allclose(H.eigh(dtype=np.complex64), H1.eigh(dtype=np.complex128))
1712+
eig1 = H1.eigh(dtype=np.complex64, **eigh_kwargs)
1713+
assert allclose(H1.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
1714+
assert allclose(
1715+
H.eigh(dtype=np.complex64, **eigh_kwargs),
1716+
H1.eigh(dtype=np.complex128, **eigh_kwargs),
1717+
)
17131718

17141719
for dtype in (np.complex64, np.complex128):
1715-
es = H1.eigenstate(dtype=dtype)
1720+
es = H1.eigenstate(dtype=dtype, **eigh_kwargs)
17161721
assert allclose(es.eig, eig1)
17171722

17181723
sm = es.spin_moment()
17191724

17201725
om = es.spin_moment(projection=True)
1721-
assert np.allclose(sm, om.sum(-1))
1726+
assert allclose(sm, om.sum(-1))
17221727

17231728
PDOS = es.PDOS(np.linspace(-1, 1, 21))
17241729
DOS = es.DOS(np.linspace(-1, 1, 21))
1725-
assert np.allclose(PDOS.sum(1)[0, :], DOS)
1730+
assert allclose(PDOS.sum(1)[0, :], DOS)
17261731
es.velocity(matrix=True)
17271732

17281733
# Check the velocities
17291734
# But only compare for np.float64, we need the precision
17301735
v = es.velocity()
17311736
vv = es.velocity(matrix=True)
1732-
assert np.allclose(np.diagonal(vv).T, v)
1737+
assert allclose(np.diagonal(vv).T, v)
17331738

17341739
# Ensure we can change gauge for NC stuff
17351740
es.change_gauge("cell")
@@ -1738,6 +1743,7 @@ def test_non_colinear_non_orthogonal(self, sisl_tolerance):
17381743
def test_spin_orbit_orthogonal(self, sisl_tolerance):
17391744
atol, rtol = sisl_tolerance[np.complex64]
17401745
allclose = partial(np.allclose, atol=atol * 10, rtol=rtol * 100)
1746+
eigh_kwargs = dict(driver="evd")
17411747

17421748
g = Geometry(
17431749
[[i, 0, 0] for i in range(10)],
@@ -1762,8 +1768,8 @@ def test_spin_orbit_orthogonal(self, sisl_tolerance):
17621768
H[i, i + 1, 0] = 1.0
17631769
H[i, i + 1, 1] = 1.0
17641770

1765-
eig1 = H.eigh(dtype=np.complex64)
1766-
assert allclose(H.eigh(dtype=np.complex128), eig1)
1771+
eig1 = H.eigh(dtype=np.complex64, **eigh_kwargs)
1772+
assert allclose(H.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
17671773
assert len(H.eigh()) == len(H)
17681774

17691775
H1 = Hamiltonian(g, dtype=np.float64, spin=Spin("spin-orbit"))
@@ -1785,36 +1791,39 @@ def test_spin_orbit_orthogonal(self, sisl_tolerance):
17851791
H1[i, i + 1, 1] = 1.0
17861792
assert H1.spsame(H)
17871793

1788-
eig1 = H1.eigh(dtype=np.complex64)
1789-
assert allclose(H1.eigh(dtype=np.complex128), eig1)
1790-
assert allclose(H.eigh(dtype=np.complex64), H1.eigh(dtype=np.complex128))
1794+
eig1 = H1.eigh(dtype=np.complex64, **eigh_kwargs)
1795+
assert allclose(H1.eigh(dtype=np.complex128, **eigh_kwargs), eig1)
1796+
assert allclose(
1797+
H.eigh(dtype=np.complex64, **eigh_kwargs),
1798+
H1.eigh(dtype=np.complex128, **eigh_kwargs),
1799+
)
17911800

17921801
# Create the block matrix for expectation
17931802
SZ = block_diag(*([H1.spin.Z] * H1.no))
17941803

17951804
for dtype in (np.complex64, np.complex128):
1796-
es = H.eigenstate(dtype=dtype)
1805+
es = H.eigenstate(dtype=dtype, **eigh_kwargs)
17971806
assert allclose(es.eig, eig1)
17981807

17991808
sm = es.spin_moment()
18001809
sm2 = es.inner(matrix=SZ).real
18011810
sm3 = np.diag(np.dot(np.conj(es.state), SZ).dot(es.state.T)).real
1802-
assert np.allclose(sm[2], sm2)
1803-
assert np.allclose(sm[2], sm3)
1811+
assert allclose(sm[2], sm2)
1812+
assert allclose(sm[2], sm3)
18041813

18051814
om = es.spin_moment(projection=True)
1806-
assert np.allclose(sm, om.sum(-1))
1815+
assert allclose(sm, om.sum(-1))
18071816

18081817
PDOS = es.PDOS(np.linspace(-1, 1, 21))
18091818
DOS = es.DOS(np.linspace(-1, 1, 21))
1810-
assert np.allclose(PDOS.sum(1)[0, :], DOS)
1819+
assert allclose(PDOS.sum(1)[0, :], DOS)
18111820
es.velocity(matrix=True)
18121821

18131822
# Check the velocities
18141823
# But only compare for np.float64, we need the precision
18151824
v = es.velocity()
18161825
vv = es.velocity(matrix=True)
1817-
assert np.allclose(np.diagonal(vv).T, v)
1826+
assert allclose(np.diagonal(vv).T, v)
18181827

18191828
# Ensure we can change gauge for SO stuff
18201829
es.change_gauge("cell")

0 commit comments

Comments
 (0)