Merge pull request #52 from JohanSchott/sanity_check_all_pickled_hami…

…ltonians

sanity check all pickled non-interacting Hamiltonians in h0 directory

.github/workflows/buildci.yml

@@ -15,7 +15,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: Set up Python 3.10
       uses: actions/setup-python@v2
       with:
@@ -38,5 +38,7 @@ jobs:
       run: |
         # Activate virtual environment, and set PYTHONPATH
         source env.sh
-        pytest
-
+        # test_comparison_with_reference.py is ignored in pytest.ini
+        # Test it separatly since it starts MPI processes in a subprocess.run.
+        pytest 
+        pytest impurityModel/test/test_comparison_with_reference.py

README.md

@@ -121,6 +121,10 @@ Type
 ```bash
 pytest
 ``` 
+and 
+```bash
+pytest impurityModel/test/test_comparison_with_reference.py
+```
 to run all python unit tests in the repository.
 
 ### Documentation

impurityModel/ed/get_spectra.py

@@ -396,7 +396,7 @@ def get_hamiltonian_operator(nBaths, nValBaths, slaterCondon, SOCs, DCinfo, hFie
     hHfieldOperator = finite.gethHfieldop(hx, hy, hz, l=2)
 
     # Read the non-relativistic non-interacting Hamiltonian operator from file.
-    h0_operator = get_h0_operator(h0_filename, nBaths)
+    h0_operator = read_pickled_file(h0_filename)
 
     # Add Hamiltonian terms to one operator.
     hOperator = finite.addOps(
@@ -416,37 +416,10 @@ def get_hamiltonian_operator(nBaths, nValBaths, slaterCondon, SOCs, DCinfo, hFie
     return hOp
 
 
-def get_h0_operator(h0_filename, nBaths):
-    """
-    Return h0 operator.
-
-    Parameters
-    ----------
-    h0_filename : str
-        Filename of non-interacting, non-relativistic operator.
-    nBaths : dict
-        Number of bath states for each angular momentum.
-
-    Returns
-    -------
-    h0_operator : dict
-        The non-relativistic non-interacting Hamiltonian in operator form.
-        Hamiltonian describes 3d orbitals and bath orbitals.
-        tuple : complex,
-        where each tuple describes a process of two steps (annihilation and then creation).
-        Each step is described by a tuple of the form:
-        (spin_orb, 'c') or (spin_orb, 'a'),
-        where spin_orb is a tuple of the form (l, s, m) or (l, b) or ((l_a, l_b), b).
-
-    """
-    with open(h0_filename, "rb") as handle:
-        h0_operator = pickle.loads(handle.read())
-    # Sanity check
-    for process in h0_operator.keys():
-        for event in process:
-            if len(event[0]) == 2:
-                assert nBaths[event[0][0]] > event[0][1]
-    return h0_operator
+def read_pickled_file(filename: str):
+    with open(filename, "rb") as handle:
+        content = pickle.load(handle)
+    return content
 
 
 if __name__ == "__main__":

impurityModel/test/test_comparison_with_reference.py

@@ -1,8 +1,10 @@
 """
 Module with test comparing new simulations with reference data.
+
+test_comparison_with_reference.py is not allowed to import any MPI stuff, except in the subprocess.run.
+Otherwise MPI gets confused, since MPI can't handle that both the parent and the child process use MPI.
 """
 
-import inspect
 import math
 import os
 import subprocess
@@ -11,50 +13,51 @@
 import h5py
 import numpy as np
 
+DIR_PATH = os.path.dirname(os.path.realpath(__file__))
+
+SCRIPT_PATH = os.path.join(DIR_PATH, "../../scripts/run_Ni_NiO_Xbath.sh")
+REFERENCE_SPECTRA_PATH = os.path.join(DIR_PATH, "referenceOutput/Ni_NiO_50bath/spectra.h5")
+
 
 def test_comparison():
     compare_spectra()
 
 
 def compare_spectra(
-    script_file="scripts/run_Ni_NiO_Xbath.sh",
+    script_path=SCRIPT_PATH,
     script_argument=50,
-    reference_file="referenceOutput/Ni_NiO_50bath/spectra.h5",
+    reference_spectra_path=REFERENCE_SPECTRA_PATH,
 ):
-    print("Start comparison of spectra...")
     # Create a temporary directory using the context manager
     with tempfile.TemporaryDirectory() as tmpdirname:
         print("Created temporary directory", tmpdirname)
         os.chdir(tmpdirname)
         print("Current working dir:", os.getcwd())
-        path = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
-        print(f"{path = }")
-        cmd = os.path.join(path[:-19], script_file)
-        print("Run command:", cmd)
-        print("Use command argument:", script_argument)
+        print(f"{script_path = }")
+        print(f"{script_argument = }")
 
-        subprocess.run(args=[cmd, str(script_argument)], stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=True)
+        subprocess.run(args=[script_path, str(script_argument)], check=True)
 
         files_and_dirs = os.listdir()
         print("Files and folders in temporary folder:", files_and_dirs)
+        assert os.path.isfile("spectra.h5")
+        assert os.path.isfile(reference_spectra_path), reference_spectra_path
         # Open spectra file and the reference spectra file
-        file_handle = h5py.File("spectra.h5", "r")
-        ref_file_handle = h5py.File(os.path.join(path, reference_file), "r")
-        # Compare file contents
-        for key in ref_file_handle:
-            print("Compare dataset:", key)
-            x = file_handle[key][()]
-            x_ref = ref_file_handle[key][()]
-            abs_diff = np.abs(x - x_ref)
-            i = np.argmax(abs_diff)
-            print("Max abs diff:", np.ravel(abs_diff)[i])
-            print("Reference value at max diff:", np.ravel(x_ref)[i])
-            np.testing.assert_allclose(x, x_ref, atol=3e-2)
-            np.testing.assert_allclose(x, x_ref, atol=2e-2, rtol=0.1)
-            print("Mean abs diff:", np.mean(abs_diff))
-            assert math.isclose(np.mean(abs_diff), 0, abs_tol=2e-5)
-            print("Median abs diff:", np.median(abs_diff))
-            assert math.isclose(np.median(abs_diff), 0, abs_tol=1e-8)
+        with h5py.File("spectra.h5", "r") as file_handle, h5py.File(reference_spectra_path, "r") as ref_file_handle:
+            # Compare file contents
+            for key in ref_file_handle:
+                print("Compare dataset:", key)
+                x = file_handle[key][()]
+                x_ref = ref_file_handle[key][()]
+                abs_diff = np.abs(x - x_ref)
+                i = np.argmax(abs_diff)
+                print("Reference value at max diff:", np.ravel(x_ref)[i])
+                np.testing.assert_allclose(x, x_ref, atol=3e-2)
+                np.testing.assert_allclose(x, x_ref, atol=2e-2, rtol=0.1)
+                print("Mean abs diff:", np.mean(abs_diff))
+                assert math.isclose(np.mean(abs_diff), 0, abs_tol=2e-5)
+                print("Median abs diff:", np.median(abs_diff))
+                assert math.isclose(np.median(abs_diff), 0, abs_tol=1e-8)
         print("Comparison successful")
 
 

impurityModel/test/test_h0.py

@@ -0,0 +1,58 @@
+import os
+from glob import glob
+
+from impurityModel.ed.get_spectra import read_pickled_file
+
+DIR_PATH = os.path.dirname(os.path.realpath(__file__))
+
+
+def test_read_all_h0_pickle_files():
+    for h0_filename in glob(os.path.join(DIR_PATH, "../../h0/h0*.pickle")):
+        h0 = read_pickled_file(h0_filename)
+
+        string = os.path.basename(h0_filename).split(".")[0].split("_")[-1].split("bath")[0]
+        nBaths_tot = sum([int(nbath) for nbath in string.split("p")])
+        # So far, in all the non-interacting Hamiltonians the angular momentum is equal to two
+        # for the correlated orbitals that are coupled to the bath states.
+        nBaths = {2: nBaths_tot}
+
+        sanity_check_non_interacting_hamiltonian(h0, nBaths)
+
+
+def sanity_check_non_interacting_hamiltonian(h0: dict[tuple, complex], nBaths: dict[int, int]):
+    """
+    Sanity check non-interacting Hamiltonian operator.
+
+    Parameters
+    -------
+    h0:
+        Non-interacting Hamiltonian operator.
+        Describes impurity orbitals and bath orbitals.
+        Each tuple key describes a process of two steps: annihilation followed by creation.
+        Each step is described by a tuple of the form:
+        (spin_orb, 'c') for creation, and
+        (spin_orb, 'a') for annihilation,
+        where spin_orb is a tuple of the form (l, s, m) or (l, b) or ((l_a, l_b), b).
+    nBaths:
+        angular momentum: number of bath states coupled to the correlated orbitals with this angular momentum.
+    """
+    for process, value in h0.items():
+        assert isinstance(value, complex)
+        assert len(process) == 2  # two events in non-interacting Hamiltonian
+        assert process[1][1] == "a"  # First event is annihilation
+        assert process[0][1] == "c"  # Second event is creation
+        for event in process[::-1]:
+            assert len(event) == 2  # spin-orbit and create or remove
+            spin_orbit_info = event[0]
+            if len(spin_orbit_info) == 2:
+                l, bath_index = spin_orbit_info
+                assert 0 <= l <= 3
+                assert 0 <= bath_index < nBaths[l]
+            elif len(spin_orbit_info) == 3:
+                l, s, m = spin_orbit_info
+                assert l in nBaths
+                assert 0 <= l <= 3
+                assert s in (0, 1)
+                assert m in range(-l, l + 1)
+            else:
+                raise ValueError(f"{spin_orbit_info}")

pytest.ini

@@ -1,6 +1,8 @@
 [pytest]
 testpaths = impurityModel
 
+# Test test_comparison_with_reference.py in a separate pytest call since it starts MPI processes in a subprocess.run, and
+# MPI can't handle that both the parent and the child process use MPI.
+addopts = --ignore impurityModel/test/test_comparison_with_reference.py
 
-filterwarnings =
-    error
+filterwarnings = error

requirements-ubuntu.txt

@@ -1,7 +1,8 @@
 build-essential
 gfortran
-openmpi-bin
-libopenmpi-dev
 python3
 python3-dev
 python3-venv
+python3-mpi4py
+openmpi-bin
+libopenmpi-dev

scripts/run_Ni_NiO_Xbath.sh

@@ -30,4 +30,4 @@ echo "Number of bath states: $nBath3d"
 echo "H0 filename: $h0_filename"
 echo "Radial wavefunction filename: $radial_filename"
 
-mpirun -n $ranks python -m impurityModel.ed.get_spectra $h0_filename $radial_filename --nBaths 0 $nBath3d --nValBaths 0 $nBath3d
+mpiexec -n $ranks --verbose python -m impurityModel.ed.get_spectra $h0_filename $radial_filename --nBaths 0 $nBath3d --nValBaths 0 $nBath3d