almost ready for packaging

pyproject.toml

@@ -6,8 +6,8 @@ build-backend = "hatchling.build"
 name = "grogupy"
 version = "0.2.0"
 authors = [
-  { name="Daniel Pozsar", email="danielpozsar@student.elte.hu" },
   { name="Laszlo Oroszlany", email="laszlo.oroszlany@ttk.elte.hu" },
+  { name="Daniel Pozsar", email="danielpozsar@student.elte.hu" },
 ]
 description = "Relativistic magnetic interactions from non-orthogonal basis sets"
 readme = "README.md"

requirements-dev.txt

@@ -10,11 +10,15 @@ black[jupyter]
 
 # dev-tools
 pre-commit
-pytest
-pytest-randomly
-pytest-cov
 isort
-tqdm
+
+# plotting
 matplotlib
 plotly
 nbformat
+
+# testing
+pytest
+pytest-randomly
+pytest-cov
+hypothesis

src/grogupy/__init__.py

@@ -21,4 +21,4 @@
 from grogupy.core import *
 from grogupy.io import *
 from grogupy.magnetism import *
-from grogupy.utils import *
+from grogupy.utilities import *

src/grogupy/core.py

@@ -22,7 +22,7 @@ import numpy as np
 from numpy.linalg import inv
 
 from grogupy.magnetism import blow_up_orbindx, parse_magnetic_entity
-from grogupy.utils import commutator
+from grogupy.utilities import commutator
 
 
 def parallel_Gk(HK, SK, eran, eset):

tests/test_core.py

@@ -0,0 +1,171 @@
+# Copyright (c) [2024] []
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import numpy as np
+import pytest
+from hypothesis import given
+from hypothesis import strategies as st
+from numpy.testing import assert_array_almost_equal
+
+from grogupy.core import (
+    build_hh_ss,
+    calc_Vu,
+    commutator,
+    onsite_projection,
+    parallel_Gk,
+    sequential_GK,
+)
+
+
+# Helper function to generate complex arrays
+def complex_arrays(shape):
+    return st.complex_numbers(
+        min_magnitude=0.0, max_magnitude=1e6, allow_infinity=False, allow_nan=False
+    ).map(lambda x: np.full(shape, x))
+
+
+# Test commutator function
+def test_commutator_hermitian():
+    """Test that commutator of Hermitian matrices is anti-Hermitian"""
+    a = np.array([[1, 1j], [-1j, 2]], dtype=np.complex128)
+    b = np.array([[3, -2j], [2j, 4]], dtype=np.complex128)
+
+    result = commutator(a, b)
+
+    # For Hermitian matrices, commutator should be anti-Hermitian
+    assert_array_almost_equal(result, -result.conj().T)
+
+
+@given(st.integers(min_value=2, max_value=10))
+def test_commutator_zero_identity(n):
+    """Test that commutator of identity matrix with any matrix is zero"""
+    # Create random complex matrix
+    rng = np.random.default_rng(42)
+    a = rng.random((n, n)) + 1j * rng.random((n, n))
+    identity = np.eye(n)
+
+    result = commutator(a, identity)
+
+    assert_array_almost_equal(result, np.zeros_like(result))
+
+
+# Test parallel_Gk and sequential_GK implementations
+@pytest.mark.parametrize("size", [2, 4, 8])
+def test_green_function_implementations_match(size):
+    """Test that parallel and sequential Green's function calculations match"""
+    rng = np.random.default_rng(42)
+
+    # Generate test inputs
+    HK = rng.random((size, size)) + 1j * rng.random((size, size))
+    HK = (HK + HK.conj().T) / 2  # Make Hermitian
+    SK = np.eye(size) + 0.1 * (rng.random((size, size)) + 1j * rng.random((size, size)))
+    SK = (SK + SK.conj().T) / 2  # Make Hermitian
+    eran = rng.random(5) + 1j * rng.random(5)
+    eset = len(eran)
+
+    # Calculate using both methods
+    G_parallel = parallel_Gk(HK, SK, eran, eset)
+    G_sequential = sequential_GK(HK, SK, eran, eset)
+
+    assert_array_almost_equal(G_parallel, G_sequential)
+
+
+# Test calc_Vu function
+def test_calc_Vu_hermiticity():
+    """Test that Vu1 and Vu2 maintain expected Hermiticity properties"""
+    rng = np.random.default_rng(42)
+    size = 4
+
+    # Create Hermitian Hamiltonian
+    H = rng.random((size, size)) + 1j * rng.random((size, size))
+    H = (H + H.conj().T) / 2
+
+    # Create unitary rotation matrix
+    Tu = rng.random((size, size)) + 1j * rng.random((size, size))
+    Tu = (Tu + Tu.conj().T) / 2
+
+    Vu1, Vu2 = calc_Vu(H, Tu)
+
+    # Vu1 should be anti-Hermitian (from commutator properties)
+    assert_array_almost_equal(Vu1, -Vu1.conj().T)
+
+    # Vu2 should be Hermitian
+    assert_array_almost_equal(Vu2, Vu2.conj().T)
+
+
+# Test onsite_projection function
+def test_onsite_projection():
+    """Test basic properties of onsite projection"""
+    size = 4
+    idx1 = np.array([0, 1])
+    idx2 = np.array([2, 3])
+
+    # Create test matrix
+    matrix = np.arange(size * size).reshape(size, size)
+
+    result = onsite_projection(matrix, idx1, idx2)
+
+    # Check shape
+    assert result.shape == (len(idx1), len(idx2))
+
+    # Check values
+    expected = matrix[np.ix_(idx1, idx2)]
+    assert_array_almost_equal(result, expected)
+
+
+@given(st.integers(min_value=2, max_value=5))
+def test_build_hh_ss_hermiticity(size):
+    """Test that built Hamiltonians maintain Hermiticity"""
+    from unittest.mock import MagicMock
+
+    # Create mock DFT Hamiltonian class
+    class MockDH:
+        def __init__(self, size):
+            self.no = size
+            self.n_s = 1
+            self.M11r = np.eye(size)
+            self.M11i = np.zeros((size, size))
+            self.M22r = np.eye(size)
+            self.M22i = np.zeros((size, size))
+            self.M12r = np.zeros((size, size))
+            self.M12i = np.zeros((size, size))
+            self.M21r = np.zeros((size, size))
+            self.M21i = np.zeros((size, size))
+            self.S_idx = np.eye(size)
+
+            self.lattice = MagicMock()
+            self.lattice.nsc.prod.return_value = 1
+
+        def tocsr(self, matrix):
+            return matrix
+
+    dh = MockDH(size)
+    hh, ss = build_hh_ss(dh)
+
+    # Check Hermiticity of Hamiltonian and overlap matrices
+    for h in hh:
+        assert_array_almost_equal(h, h.conj().T)
+
+    for s in ss:
+        assert_array_almost_equal(s, s.conj().T)
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])

tests/test_io.py

@@ -0,0 +1,254 @@
+# Copyright (c) [2024] []
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import io
+import os
+import tempfile
+from contextlib import redirect_stdout
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from grogupy.io import (
+    default_args,
+    load_pickle,
+    print_atoms_and_pairs,
+    print_job_description,
+    print_parameters,
+    print_runtime_information,
+    save_pickle,
+)
+
+
+# Fixtures for common test data
+@pytest.fixture
+def simulation_parameters():
+    """Create sample simulation parameters for testing"""
+    params = default_args.copy()
+    params.update(
+        {
+            "infile": "test.fdf",
+            "outfile": "test_output.pickle",
+            "parallel_size": 4,
+            "cell": np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]),
+            "automatic_ebot": True,
+        }
+    )
+    return params
+
+
+@pytest.fixture
+def magnetic_entities():
+    """Create sample magnetic entities for testing"""
+    return [
+        {
+            "tags": ["[0]Fe(d)"],
+            "xyz": [np.array([0.0, 0.0, 0.0])],
+            "K": np.array([1.0, 2.0, 3.0]),
+            "K_consistency": 0.001,
+        },
+        {
+            "tags": ["[1]Fe(d)"],
+            "xyz": [np.array([1.0, 1.0, 1.0])],
+            "K": np.array([2.0, 3.0, 4.0]),
+            "K_consistency": 0.002,
+        },
+    ]
+
+
+@pytest.fixture
+def pairs():
+    """Create sample pairs for testing"""
+    return [
+        {
+            "tags": ["[0]Fe(d)", "[1]Fe(d)"],
+            "Ruc": np.array([0, 0, 0]),
+            "dist": 1.732,
+            "J_iso": 10.0,
+            "D": np.array([0.1, 0.2, 0.3]),
+            "J_S": np.array([0.5, 0.6, 0.7, 0.8, 0.9]),
+            "J": np.array([[1.0, 0.1, 0.2], [0.1, 1.1, 0.3], [0.2, 0.3, 1.2]]),
+        }
+    ]
+
+
+@pytest.fixture
+def runtime_info():
+    """Create sample runtime information for testing"""
+    return {
+        "start_time": 0.0,
+        "setup_time": 1.0,
+        "H_and_XCF_time": 2.0,
+        "site_and_pair_dictionaries_time": 3.0,
+        "k_set_time": 4.0,
+        "reference_rotations_time": 5.0,
+        "green_function_inversion_time": 6.0,
+        "end_time": 7.0,
+    }
+
+
+# Test pickle save/load functionality
+def test_pickle_save_load(simulation_parameters):
+    """Test saving and loading data using pickle"""
+    with tempfile.NamedTemporaryFile(suffix=".pickle", delete=False) as tmp:
+        temp_path = tmp.name
+
+    try:
+        # Test saving
+        save_pickle(temp_path, simulation_parameters)
+        assert os.path.exists(temp_path)
+        assert os.path.getsize(temp_path) > 0
+
+        # Test loading
+        loaded_data = load_pickle(temp_path)
+        assert loaded_data == simulation_parameters
+
+    finally:
+        # Cleanup
+        if os.path.exists(temp_path):
+            os.unlink(temp_path)
+
+
+# Test print functions
+def test_print_parameters(simulation_parameters):
+    """Test parameters printing function"""
+    output = io.StringIO()
+    with redirect_stdout(output):
+        print_parameters(simulation_parameters)
+
+    printed_output = output.getvalue()
+
+    # Check key elements are present in output
+    assert "Input file:" in printed_output
+    assert "Cell [Ang]:" in printed_output
+    assert "DFT axis:" in printed_output
+    assert str(simulation_parameters["kset"]) in printed_output
+
+
+def test_print_atoms_and_pairs(magnetic_entities, pairs):
+    """Test atoms and pairs printing function"""
+    output = io.StringIO()
+    with redirect_stdout(output):
+        print_atoms_and_pairs(magnetic_entities, pairs)
+
+    printed_output = output.getvalue()
+
+    # Check key elements are present in output
+    assert "Atomic information:" in printed_output
+    assert "Anisotropy [meV]" in printed_output
+    assert "Exchange [meV]" in printed_output
+
+    # Check if magnetic entity information is printed
+    for entity in magnetic_entities:
+        assert entity["tags"][0] in printed_output
+
+    # Check if pair information is printed
+    for pair in pairs:
+        assert pair["tags"][0] in printed_output
+        assert pair["tags"][1] in printed_output
+
+
+def test_print_runtime_information(runtime_info):
+    """Test runtime information printing function"""
+    output = io.StringIO()
+    with redirect_stdout(output):
+        print_runtime_information(runtime_info)
+
+    printed_output = output.getvalue()
+
+    # Check key elements are present in output
+    assert "Runtime information:" in printed_output
+    assert "Total runtime:" in printed_output
+    assert "Initial setup:" in printed_output
+
+
+def test_print_job_description(simulation_parameters):
+    """Test job description printing function"""
+    output = io.StringIO()
+    with redirect_stdout(output):
+        print_job_description(simulation_parameters)
+
+    printed_output = output.getvalue()
+
+    # Check key elements are present in output
+    assert "Input file:" in printed_output
+    assert "Number of nodes in the parallel cluster:" in printed_output
+    assert "Cell [Ang]:" in printed_output
+    assert "Parameters for the contour integral:" in printed_output
+
+
+# Test default arguments
+def test_default_args_structure():
+    """Test the structure and values of default arguments"""
+    assert isinstance(default_args, dict)
+    assert "infile" in default_args
+    assert "outfile" in default_args
+    assert "kset" in default_args
+    assert "kdirs" in default_args
+    assert "eset" in default_args
+    assert "esetp" in default_args
+
+    # Test specific default values
+    assert default_args["kset"] == 2
+    assert default_args["kdirs"] == "xyz"
+    assert default_args["eset"] == 42
+    assert default_args["esetp"] == 1000
+    assert default_args["parallel_solver_for_Gk"] is False
+    assert default_args["padawan_mode"] is True
+
+
+def test_simulation_parameters_validation(simulation_parameters):
+    """Test validation of simulation parameters"""
+    # Test required fields
+    required_fields = ["infile", "outfile", "kset", "kdirs", "eset", "esetp"]
+    for field in required_fields:
+        assert field in simulation_parameters
+
+    # Test numeric parameters are positive
+    assert simulation_parameters["kset"] > 0
+    assert simulation_parameters["eset"] > 0
+    assert simulation_parameters["esetp"] > 0
+
+    # Test cell is 3x3 array
+    assert simulation_parameters["cell"].shape == (3, 3)
+
+
+@pytest.mark.parametrize("test_dir", ["x", "y", "z", "xy", "yz", "xz", "xyz"])
+def test_kdirs_validation(simulation_parameters, test_dir):
+    """Test validation of kdirs parameter"""
+    simulation_parameters["kdirs"] = test_dir
+    # Should not raise any exceptions
+    print_parameters(simulation_parameters)
+
+
+def test_ebot_automatic_detection(simulation_parameters):
+    """Test automatic ebot detection flag"""
+    assert "automatic_ebot" in simulation_parameters
+    output = io.StringIO()
+    with redirect_stdout(output):
+        print_job_description(simulation_parameters)
+
+    printed_output = output.getvalue()
+    assert "WARNING: This was automatically determined!" in printed_output
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])

tests/test_jij.py

@@ -1,19 +0,0 @@
-# Copyright (c) [2024] [Daniel Pozsar]
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.

tests/test_magnetism.py

@@ -0,0 +1,240 @@
+# Copyright (c) [2024] []
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from unittest.mock import Mock
+
+import numpy as np
+import pytest
+from hypothesis import given
+from hypothesis import strategies as st
+from numpy.testing import assert_allclose, assert_array_almost_equal
+
+from grogupy.magnetism import (
+    blow_up_orbindx,
+    calculate_anisotropy_tensor,
+    calculate_exchange_tensor,
+    parse_magnetic_entity,
+    spin_tracer,
+)
+
+
+# Utility to create mock sisl atoms with orbitals
+def create_mock_atom(orbitals):
+    atom = Mock()
+    atom.orbitals = []
+    for l in orbitals:
+        orbital = Mock()
+        orbital.l = l
+        atom.orbitals.append(orbital)
+    return atom
+
+
+# Test blow_up_orbindx function
+def test_blow_up_orbindx_basic():
+    """Test basic orbital index expansion"""
+    orb_indices = np.array([0, 1, 2])
+    result = blow_up_orbindx(orb_indices)
+    expected = np.array([0, 1, 2, 3, 4, 5])  # Each index expands to two indices
+    assert_array_almost_equal(result, expected)
+
+
+@given(st.lists(st.integers(min_value=0, max_value=100), min_size=1, max_size=10))
+def test_blow_up_orbindx_properties(indices):
+    """Property-based test for orbital index expansion"""
+    orb_indices = np.array(indices)
+    result = blow_up_orbindx(orb_indices)
+
+    # Result should be twice as long as input
+    assert len(result) == 2 * len(orb_indices)
+
+    # Each pair should be consecutive numbers
+    for i in range(0, len(result), 2):
+        assert result[i + 1] == result[i] + 1
+
+
+# Test spin_tracer function
+def test_spin_tracer_basic():
+    """Test basic spin tracing functionality"""
+    # Create a test matrix in SPIN-BOX representation
+    size = 2
+    M = np.zeros((2 * size, 2 * size), dtype=complex)
+    M[0::2, 0::2] = 1.0  # M11
+    M[1::2, 1::2] = -1.0  # M22
+    M[0::2, 1::2] = 0.5j  # M12
+    M[1::2, 0::2] = -0.5j  # M21
+
+    result = spin_tracer(M)
+
+    # Check components
+    assert_array_almost_equal(result["z"], np.eye(size) * 2.0)  # M11 - M22
+    assert_array_almost_equal(result["x"], np.eye(size))  # M12 + M21
+    assert_array_almost_equal(result["y"], np.eye(size) * 1.0j)  # i(M12 - M21)
+    assert_array_almost_equal(result["c"], np.zeros((size, size)))  # M11 + M22
+
+
+def test_spin_tracer_hermiticity():
+    """Test that spin tracer preserves Hermiticity"""
+    size = 2
+    M = np.random.random((2 * size, 2 * size)) + 1j * np.random.random(
+        (2 * size, 2 * size)
+    )
+    M = M + M.conj().T  # Make Hermitian
+
+    result = spin_tracer(M)
+
+    # Check Hermiticity of components
+    for component in ["x", "y", "z", "c"]:
+        assert_array_almost_equal(result[component], result[component].conj().T)
+
+
+# Test parse_magnetic_entity function
+def test_parse_magnetic_entity_single_atom():
+    """Test parsing of single atom magnetic entity"""
+    # Create mock geometry
+    geometry = Mock()
+    atom = create_mock_atom([0, 2, 2, 2])  # s and d orbitals
+    geometry.atoms = [atom]
+    geometry.a2o.return_value = np.array([0, 1, 2, 3])
+
+    dh = Mock()
+    dh.geometry = geometry
+
+    # Test with d orbitals only
+    entity = {"atom": 0, "l": 2}
+    result = parse_magnetic_entity(dh, **entity)
+    assert len(result) == 3  # Should only include d orbitals
+
+    # Test with all orbitals
+    entity = {"atom": 0, "l": None}
+    result = parse_magnetic_entity(dh, **entity)
+    assert len(result) == 4  # Should include all orbitals
+
+
+def test_parse_magnetic_entity_multiple_atoms():
+    """Test parsing of multiple atom magnetic entity"""
+    # Create mock geometry with two atoms
+    geometry = Mock()
+    atom1 = create_mock_atom([0, 2])
+    atom2 = create_mock_atom([0, 2])
+    geometry.atoms = [atom1, atom2]
+    geometry.a2o.side_effect = [np.array([0, 1]), np.array([2, 3])]
+
+    dh = Mock()
+    dh.geometry = geometry
+
+    entity = {"atom": [0, 1], "l": 2}
+    result = parse_magnetic_entity(dh, **entity)
+    assert len(result) == 2  # Should include d orbitals from both atoms
+
+
+# Test calculate_anisotropy_tensor function
+def test_calculate_anisotropy_tensor():
+    """Test calculation of anisotropy tensor"""
+    # Create mock magnetic entity with energies
+    mag_ent = {
+        "energies": np.array(
+            [
+                [0.1, 0.2],  # First rotation axis
+                [0.3, 0.4],  # Second rotation axis
+                [0.5, 0.6],  # Third rotation axis
+            ]
+        )
+    }
+
+    Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent)
+
+    # Check basic properties
+    assert isinstance(Kxx, float)
+    assert isinstance(Kyy, float)
+    assert isinstance(Kzz, float)
+    assert isinstance(consistency, float)
+
+    # Verify calculations
+    assert_allclose(Kxx, mag_ent["energies"][1, 1] - mag_ent["energies"][1, 0])
+    assert_allclose(Kyy, mag_ent["energies"][0, 1] - mag_ent["energies"][0, 0])
+    assert_allclose(Kzz, 0)
+
+
+# Test calculate_exchange_tensor function
+def test_calculate_exchange_tensor():
+    """Test calculation of exchange tensor"""
+    # Create mock pair with energies
+    pair = {
+        "energies": np.zeros((3, 4))  # 3 quantization axes, 4 rotation combinations
+    }
+    # Set some test values
+    pair["energies"][0] = [0.1, 0.2, 0.3, 0.4]  # First quantization axis
+    pair["energies"][1] = [0.2, 0.3, 0.4, 0.5]  # Second quantization axis
+    pair["energies"][2] = [0.3, 0.4, 0.5, 0.6]  # Third quantization axis
+
+    J_iso, J_S, D, J = calculate_exchange_tensor(pair)
+
+    # Check dimensions
+    assert isinstance(J_iso, float)
+    assert J_S.shape == (5,)  # Symmetric exchange has 5 independent components
+    assert D.shape == (3,)  # DM vector has 3 components
+    assert J.shape == (3, 3)  # Full exchange tensor is 3x3
+
+    # Check properties
+    assert_allclose(np.trace(J) / 3, J_iso)  # Isotropic exchange
+    assert_array_almost_equal(J, J.T)  # Symmetric part
+
+    # Verify DM vector components
+    assert_allclose(D[2], (pair["energies"][2, 1] - pair["energies"][2, 2]) / 2)
+
+
+def test_exchange_tensor_symmetries():
+    """Test symmetry properties of exchange tensor calculations"""
+    pair = {"energies": np.zeros((3, 4))}
+    # Set symmetric energies
+    pair["energies"][0] = [0.1, 0.1, 0.1, 0.1]
+    pair["energies"][1] = [0.2, 0.2, 0.2, 0.2]
+    pair["energies"][2] = [0.3, 0.3, 0.3, 0.3]
+
+    J_iso, J_S, D, J = calculate_exchange_tensor(pair)
+
+    # For symmetric energies, DM vector should be zero
+    assert_array_almost_equal(D, np.zeros(3))
+
+    # Exchange tensor should be symmetric
+    assert_array_almost_equal(J, J.T)
+
+
+@pytest.mark.parametrize(
+    "test_energies,expected_D",
+    [
+        (
+            np.array(
+                [[0.1, 0.2, -0.2, 0.1], [0.1, 0.2, -0.2, 0.1], [0.1, 0.2, -0.2, 0.1]]
+            ),
+            np.array([0.2, 0.2, 0.2]),
+        ),
+        (np.zeros((3, 4)), np.zeros(3)),
+    ],
+)
+def test_exchange_tensor_dm_vector(test_energies, expected_D):
+    """Test DM vector calculation with different energy configurations"""
+    pair = {"energies": test_energies}
+    _, _, D, _ = calculate_exchange_tensor(pair)
+    assert_array_almost_equal(D, expected_D)
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])

tests/test_useful.py

@@ -1,80 +0,0 @@
-# Copyright (c) [2024] [Daniel Pozsar]
-#
-# Permission is hereby granted, free of charge, to any person obtaining a copy
-# of this software and associated documentation files (the "Software"), to deal
-# in the Software without restriction, including without limitation the rights
-# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-# copies of the Software, and to permit persons to whom the Software is
-# furnished to do so, subject to the following conditions:
-#
-# The above copyright notice and this permission notice shall be included in all
-# copies or substantial portions of the Software.
-#
-# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-# SOFTWARE.
-
-import pytest
-from useful import *
-
-
-@pytest.mark.parametrize(
-    "a, b, c",
-    [
-        (tau_x, tau_y, tau_z),
-        (tau_z, tau_x, tau_y),
-        (tau_y, tau_z, tau_x),
-        (tau_z, tau_y, -tau_x),
-    ],
-)
-def test_comm(a, b, c):
-    assert (comm(a, b) == 2j * c).all()
-
-
-@pytest.mark.parametrize(
-    "dirs, NUMK, out",
-    [
-        ("", 10, np.array([[0, 0, 0]])),
-        (
-            "x",
-            10,
-            np.array(
-                [
-                    [0.0, 0.0, 0.0],
-                    [0.1, 0.0, 0.0],
-                    [0.2, 0.0, 0.0],
-                    [0.3, 0.0, 0.0],
-                    [0.4, 0.0, 0.0],
-                    [0.5, 0.0, 0.0],
-                    [0.6, 0.0, 0.0],
-                    [0.7, 0.0, 0.0],
-                    [0.8, 0.0, 0.0],
-                    [0.9, 0.0, 0.0],
-                ]
-            ),
-        ),
-        (
-            "xy",
-            3,
-            np.array(
-                [
-                    [0.0, 0.0, 0.0],
-                    [0.33333333, 0.0, 0.0],
-                    [0.66666667, 0.0, 0.0],
-                    [0.0, 0.33333333, 0.0],
-                    [0.0, 0.66666667, 0.0],
-                    [0.33333333, 0.33333333, 0.0],
-                    [0.33333333, 0.66666667, 0.0],
-                    [0.66666667, 0.33333333, 0.0],
-                    [0.66666667, 0.66666667, 0.0],
-                ]
-            ),
-        ),
-    ],
-)
-def test_make_kset(dirs, NUMK, out):
-    assert (make_kset(dirs, NUMK) == out).all()

tests/test_utilities.py

@@ -0,0 +1,264 @@
+# Copyright (c) [2024] []
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import numpy as np
+import pytest
+from hypothesis import given
+from hypothesis import strategies as st
+from numpy.testing import assert_allclose, assert_array_almost_equal
+
+from grogupy.utilities import (
+    RotM,
+    RotMa2b,
+    crossM,
+    hsk,
+    int_de_ke,
+    make_contour,
+    make_kset,
+    read_siesta_emin,
+    tau_0,
+    tau_u,
+    tau_x,
+    tau_y,
+    tau_z,
+)
+
+
+# Test basic Pauli matrix properties
+def test_pauli_matrices_properties():
+    """Test fundamental properties of Pauli matrices"""
+    # Test anticommutation relations
+    assert_array_almost_equal(tau_x @ tau_y + tau_y @ tau_x, np.zeros((2, 2)))
+    assert_array_almost_equal(tau_y @ tau_z + tau_z @ tau_y, np.zeros((2, 2)))
+    assert_array_almost_equal(tau_z @ tau_x + tau_x @ tau_z, np.zeros((2, 2)))
+
+    # Test square of Pauli matrices equals identity
+    assert_array_almost_equal(tau_x @ tau_x, tau_0)
+    assert_array_almost_equal(tau_y @ tau_y, tau_0)
+    assert_array_almost_equal(tau_z @ tau_z, tau_0)
+
+    # Test Hermiticity
+    assert_array_almost_equal(tau_x, tau_x.conj().T)
+    assert_array_almost_equal(tau_y, tau_y.conj().T)
+    assert_array_almost_equal(tau_z, tau_z.conj().T)
+
+
+# Test hsk function
+def test_hsk_basics():
+    """Test basic properties of the hsk function"""
+    size = 4
+    H = np.random.random((3, size, size)) + 1j * np.random.random((3, size, size))
+    ss = np.random.random((3, size, size)) + 1j * np.random.random((3, size, size))
+    sc_off = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]])
+    k = np.array([0.1, 0.2, 0.3])
+
+    HK, SK = hsk(H, ss, sc_off, k)
+
+    # Check shapes
+    assert HK.shape == (size, size)
+    assert SK.shape == (size, size)
+
+    # Test k=0 case gives the central cell contribution
+    HK0, SK0 = hsk(H, ss, sc_off, np.zeros(3))
+    assert_array_almost_equal(HK0, H[0])
+    assert_array_almost_equal(SK0, ss[0])
+
+
+@given(st.integers(min_value=1, max_value=20))
+def test_make_kset_dimensions(num_k):
+    """Test k-point set generation with different dimensions"""
+    # Test 1D
+    kset_x = make_kset("x", num_k)
+    assert kset_x.shape[1] == 3
+    assert len(np.unique(kset_x[:, 1])) == 1  # y coordinates should be same
+    assert len(np.unique(kset_x[:, 2])) == 1  # z coordinates should be same
+
+    # Test 2D
+    kset_xy = make_kset("xy", num_k)
+    assert kset_xy.shape[1] == 3
+    assert len(np.unique(kset_xy[:, 2])) == 1  # z coordinates should be same
+
+    # Test 3D
+    kset_xyz = make_kset("xyz", num_k)
+    assert kset_xyz.shape[1] == 3
+    assert len(kset_xyz) == num_k**3
+
+
+def test_make_contour():
+    """Test contour generation for energy integration"""
+    emin = -20
+    emax = 0.0
+    enum = 42
+    p = 150
+
+    cont = make_contour(emin, emax, enum, p)
+
+    # Check attributes
+    assert hasattr(cont, "R")
+    assert hasattr(cont, "z0")
+    assert hasattr(cont, "ze")
+    assert hasattr(cont, "we")
+    assert hasattr(cont, "enum")
+
+    # Check dimensions
+    assert len(cont.ze) == enum
+    assert len(cont.we) == enum
+
+    # Check basic properties
+    assert cont.z0 == (emax + emin) / 2
+    assert cont.R == (emax - emin) / 2
+    assert cont.enum == enum
+
+
+@given(st.lists(st.floats(min_value=-1, max_value=1), min_size=3, max_size=3))
+def test_tau_u(u):
+    """Test generation of Pauli matrix in arbitrary direction"""
+    u = np.array(u)
+    if np.all(u == 0):
+        u = np.array([1, 0, 0])  # Handle zero vector case
+    u = u / np.linalg.norm(u)
+
+    result = tau_u(u)
+
+    # Check result is 2x2
+    assert result.shape == (2, 2)
+
+    # Check Hermiticity
+    assert_array_almost_equal(result, result.conj().T)
+
+    # Check trace is zero
+    assert abs(np.trace(result)) < 1e-10
+
+    # Check square equals identity
+    assert_array_almost_equal(result @ result, np.eye(2))
+
+
+def test_crossM():
+    """Test cross product matrix generation"""
+    u = np.array([1.0, 2.0, 3.0])
+    v = np.array([4.0, 5.0, 6.0])
+
+    # Test that crossM(u) @ v equals np.cross(u, v)
+    assert_array_almost_equal(crossM(u) @ v, np.cross(u, v))
+
+    # Test antisymmetry
+    assert_array_almost_equal(crossM(u), -crossM(u).T)
+
+
+@pytest.mark.parametrize("theta", [0, np.pi / 4, np.pi / 2, np.pi])
+def test_RotM(theta):
+    """Test rotation matrix generation"""
+    # Test rotation around z-axis
+    u = np.array([0, 0, 1])
+    R = RotM(theta, u)
+
+    # Check orthogonality
+    assert_array_almost_equal(R @ R.T, np.eye(3))
+
+    # Check determinant is 1
+    assert_allclose(np.linalg.det(R), 1.0)
+
+    # Test rotation of vector
+    v = np.array([1, 0, 0])
+    rotated_v = R @ v
+    expected_v = np.array([np.cos(theta), np.sin(theta), 0])
+    assert_array_almost_equal(rotated_v, expected_v)
+
+
+def test_RotMa2b():
+    """Test rotation matrix between two vectors"""
+    # Test simple cases
+    a = np.array([1, 0, 0])
+    b = np.array([0, 1, 0])
+    R = RotMa2b(a, b)
+
+    # Check that R rotates a to b
+    assert_array_almost_equal(R @ a, b)
+
+    # Check orthogonality
+    assert_array_almost_equal(R @ R.T, np.eye(3))
+
+    # Check determinant is 1
+    assert_allclose(np.linalg.det(R), 1.0)
+
+
+def test_read_siesta_emin(tmpdir):
+    """Test reading minimum energy from SIESTA .EIG file"""
+    # Create mock .EIG file
+    eig_file = tmpdir / "test.EIG"
+    with open(eig_file, "w") as f:
+        # Write mock eigenvalues
+        f.write("1.0 2.0 3.0\n-1.0 0.0 1.0\n")
+
+    # Test reading minimum energy
+    emin = read_siesta_emin(str(eig_file))
+    assert emin == -1.0
+
+
+def test_int_de_ke():
+    """Test energy integral calculation"""
+    # Create test data
+    energy = np.linspace(-10, 10, 1000)
+    traced = np.exp(-(energy**2))  # Gaussian function
+    we = np.ones_like(energy)  # Unit weights
+
+    result = int_de_ke(traced, we)
+
+    # Result should be real for real input
+    assert np.imag(result) < 1e-10
+
+    # Test with complex input
+    traced_complex = traced + 1j * traced
+    result_complex = int_de_ke(traced_complex, we)
+    assert isinstance(result_complex, complex)
+
+
+# Property-based tests
+@given(st.integers(min_value=2, max_value=10))
+def test_make_kset_dimensions_property(num_k):
+    """Property-based test for k-point set generation"""
+    kset = make_kset("xyz", num_k)
+
+    # All k-points should be in [0,1)
+    assert np.all(kset >= 0)
+    assert np.all(kset < 1)
+
+    # Shape should be (num_k^3, 3)
+    assert kset.shape == (num_k**3, 3)
+
+
+@given(
+    st.floats(min_value=-20, max_value=-1),
+    st.floats(min_value=0, max_value=1),
+    st.integers(min_value=10, max_value=100),
+)
+def test_make_contour_properties(emin, emax, enum):
+    """Property-based test for contour generation"""
+    cont = make_contour(emin, emax, enum, 150)
+
+    # Check basic properties
+    assert len(cont.ze) == enum
+    assert len(cont.we) == enum
+    assert cont.R > 0
+    assert cont.z0 == (emax + emin) / 2
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])