diff --git a/pyproject.toml b/pyproject.toml index ccb10fa..2be54db 100644 --- a/pyproject.toml +++ b/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" diff --git a/requirements-dev.txt b/requirements-dev.txt index 3acf141..4f75850 100644 --- a/requirements-dev.txt +++ b/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 diff --git a/src/grogupy/__init__.py b/src/grogupy/__init__.py index 3864f4b..b6939fc 100644 --- a/src/grogupy/__init__.py +++ b/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 * diff --git a/src/grogupy/core.py b/src/grogupy/core.py index 291a659..125bfae 100644 --- a/src/grogupy/core.py +++ b/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): diff --git a/src/grogupy/utils.py b/src/grogupy/utilities.py similarity index 100% rename from src/grogupy/utils.py rename to src/grogupy/utilities.py diff --git a/tests/test_core.py b/tests/test_core.py new file mode 100644 index 0000000..350d949 --- /dev/null +++ b/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__]) diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..5caa549 --- /dev/null +++ b/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__]) diff --git a/tests/test_jij.py b/tests/test_jij.py deleted file mode 100644 index 93d0ca4..0000000 --- a/tests/test_jij.py +++ /dev/null @@ -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. diff --git a/tests/test_magnetism.py b/tests/test_magnetism.py new file mode 100644 index 0000000..b0a0272 --- /dev/null +++ b/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__]) diff --git a/tests/test_useful.py b/tests/test_useful.py deleted file mode 100644 index 90016ab..0000000 --- a/tests/test_useful.py +++ /dev/null @@ -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() diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 0000000..5ddd221 --- /dev/null +++ b/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__])