Clarified mapping functions, added the unittest for the function to compute overlap of Slater determinants. Since now, the convention is to use the doubled overlap matrix - alpha and beta blocks (user_notation = 1) in sd2indx and to do the sorting, so the overlap doesn't depend on the particular ordering of spin-orbitals entered by user, but only on their identities

Author: alexvakimov

src/libra_py/workflows/nbra/mapping.py

@@ -112,7 +112,7 @@ def sd2indx(inp, nbasis=0, do_sort=False, user_notation=0):
 
 
 
-def ovlp_arb(_SD1, _SD2, S, active_space=None, do_sort=False, user_notation=0, verbose=False):
+def ovlp_arb(_SD1, _SD2, S, active_space=None, verbose=False):
     """Compute the overlap of two generic SDs: <SD1|SD2>
     This functions translates the explicit use of two Python `for` loops
     into a more efficient way using numpy
@@ -129,15 +129,6 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, do_sort=False, user_notation=0, v
         active_space ( list of ints ): indices of the orbitals (starting from 1) to
             include into consideration. If None - all the orbitals will be used [ default: None ]
 
-        do_sort ( Boolean ): the flag to tell whether the indices should be sorted in
-            a particular order (according to canonic ordering defined in the `citools.slatdet`):
-                - True - for new scheme (needed for the SAC) [default]
-                - use with False for Pyxaid mapping!
-
-        user_notation (int):
-            - 0 : short-hand - in this case, the mapping goes into [0, N) range [default - because it has been used for a while]
-            - 1 : extended - in this case, the mapping goes into [0, 2*N) range
-
         verbose ( Bool ): whether to print some extra info [ default: False] 
 
     Returns:
@@ -169,33 +160,27 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, do_sort=False, user_notation=0, v
     # Apply the determinant formula
     det_size = len(SD1)
     s = np.zeros( (det_size, det_size), dtype=np.float64); 
-    """ Commented on 7/22/2025
-    # ============================== The numpy version of the above double-for-loops
-    # Find the tensor product of the SD1 and SD2
-    SD_tensor_product = np.tensordot(SD1, SD2, axes=0)
-    # Next, find the sign of these elementwise multiplications
-    SD_tensor_product_sign = np.sign(SD_tensor_product)
-    # Now, find where we have alpha-beta indices so that the 
-    negative_sign_indices = np.where(SD_tensor_product_sign < 0)
-    s[negative_sign_indices] = 0
-    """
+    
     # Let's build the matrix related to sd1 and sd2 from the KS orbitals
-    # For this, we reuire to turn each element into matrix indices
-    sd1,_,_ = sd2indx(SD1, nbasis, do_sort, user_notation)
-    sd2,_,_ = sd2indx(SD2, nbasis, do_sort, user_notation)
+    # For this, we require to turn each element into matrix indices
+    sd1, p1 = sd2indx(SD1, nbasis, do_sort=True, user_notation=1)
+    sd2, p2 = sd2indx(SD2, nbasis, do_sort=True, user_notation=1)
 
+    if verbose==True:
+        print(sd1, sd2)
+        print(p1, p2)
 
     """ 
-    ALEXEY: Instead of this trick, use user_notation = 1
+    #ALEXEY on 12/08/2025: Instead of this trick, we use user_notation = 1 above
 
     # What about beta indices?! We should add `nbasis/2` to them. This is added on 7/22/2025
     # With this simple approach, there is no need for tensor product or the `if else` clause brought
     # previously since the elements corresponding to the alpha-beta orbitals are already
     # zero in the two-spinor format of the KS matrices :)
+
     beta_indices = np.where(np.array(SD1) < 0)
     sd1 = np.array(sd1)
     sd1[beta_indices] += int(nbasis/2)
-
     beta_indices = np.where(np.array(SD2) < 0)
     sd2 = np.array(sd2)
     sd2[beta_indices] += int(nbasis/2)
@@ -207,13 +192,14 @@ def ovlp_arb(_SD1, _SD2, S, active_space=None, do_sort=False, user_notation=0, v
 
     if verbose==True:
         print(s)
+
     res = np.linalg.det(s)
 
     return res
 
 
 
-def ovlp_mat_arb(SD1, SD2, _S, active_space=None, do_sort=False, user_notation=0, verbose=False):
+def ovlp_mat_arb(SD1, SD2, _S, active_space=None, verbose=False):
     """Compute a matrix of overlaps in the SD basis
 
     Args:
@@ -233,15 +219,6 @@ def ovlp_mat_arb(SD1, SD2, _S, active_space=None, do_sort=False, user_notation=0
         active_space ( list of ints ): indices of the orbitals (starting from 1) to
             include into consideration. If None - all the orbitals will be used [ default: None ]
 
-        do_sort ( Boolean ): the flag to tell whether the indices should be sorted in
-            a particular order (according to canonic ordering defined in the `citools.slatdet`):
-                - True - for new scheme (needed for the SAC) [default]
-                - use with False for Pyxaid mapping!
-
-        user_notation (int):
-            - 0 : short-hand - in this case, the mapping goes into [0, N) range [default - because it has been used for a while]
-            - 1 : extended - in this case, the mapping goes into [0, 2*N) range
-
         verbose ( Bool ): whether to print some extra info [ default: False]
 
     Returns:
@@ -261,7 +238,7 @@ def ovlp_mat_arb(SD1, SD2, _S, active_space=None, do_sort=False, user_notation=0
 
     for n in range(0, N):
         for m in range(0, M):
-            res[n, m] = ovlp_arb(SD1[n], SD2[m], S, active_space, do_sort, user_notation, verbose)
+            res[n, m] = ovlp_arb(SD1[n], SD2[m], S, active_space, verbose)
 
     return res
 

unittests/test_mapping.py

@@ -26,8 +26,10 @@
 elif sys.platform == "linux" or sys.platform == "linux2":
     from liblibra_core import *
 
-from libra_py.workflows.nbra.mapping import sd2indx
+from libra_py.workflows.nbra.mapping import sd2indx, ovlp_arb
 
+
+#=============================== testing sd2indx ============================
 def generate():
     """
     Run this function to see what is the current results 
@@ -43,6 +45,7 @@ def generate():
                 out, p = sd2indx(inp, nbasis=6, do_sort=do_sort, user_notation=user_notation)
                 print(F"{inp} do_sort={do_sort}, user_notation={user_notation} -> {out}, {p}" )
 
+#generate()
 
 # Test cases
 test_data = [
@@ -84,3 +87,80 @@ def test_sd2indx(inp, do_sort, user_notation, expected_out, expected_p):
     assert out == expected_out, f"Output mismatch for inp={inp}, do_sort={do_sort}, user_notation={user_notation}"
     assert p == expected_p, f"Phase mismatch for inp={inp}, do_sort={do_sort}, user_notation={user_notation}"
 
+
+#================================== testing ovlp_arb =================================
+def generate2():
+    S = np.eye(6)
+    S[0,0] = 0.91
+    S[1,1] = -0.98
+    S[3,3] = 0.91
+    S[4,4] = -0.98
+
+    basis = [ [1,-1, 2, -2],
+                 [-1,1, 2, -2],
+                 [1,-1, 2, -3],
+                 [1,-1,-2, 3]
+               ]
+    for inp1 in basis:
+        for inp2 in basis:
+            # We should always sort - to get user-input-independent results
+            # We always use doubled matrix notation - so that the translations work as intended
+            out1, _ = sd2indx(inp1, nbasis=6, do_sort=do_sort, user_notation=user_notation)
+            out2, _ = sd2indx(inp2, nbasis=6, do_sort=do_sort, user_notation=user_notation)
+            x = ovlp_arb(inp1, inp2, S, None, False)
+            print(F"<{inp1}|{inp2}> = <{out1}|{out2}> = {x}")
+
+
+#generate2()
+
+
+basis = [
+    [1, -1, 2, -2],
+    [-1, 1, 2, -2],
+    [1, -1, 2, -3],
+    [1, -1, -2, 3],
+]
+
+
+# Expected results
+expected = {
+    ((1,-1, 2,-2), (1,-1, 2,-2)): 0.79530724,
+    ((1,-1, 2,-2), (-1,1, 2,-2)): 0.79530724,
+    ((1,-1, 2,-2), (1,-1, 2,-3)): 0.0,
+    ((1,-1, 2,-2), (1,-1,-2, 3)): 0.0,
+
+    ((-1,1, 2,-2), (1,-1, 2,-2)): 0.79530724,
+    ((-1,1, 2,-2), (-1,1, 2,-2)): 0.79530724,
+    ((-1,1, 2,-2), (1,-1, 2,-3)): 0.0,
+    ((-1,1, 2,-2), (1,-1,-2, 3)): 0.0,
+
+    ((1,-1, 2,-3), (1,-1, 2,-2)): 0.0,
+    ((1,-1, 2,-3), (-1,1, 2,-2)): 0.0,
+    ((1,-1, 2,-3), (1,-1, 2,-3)): -0.8115380000000001,
+    ((1,-1, 2,-3), (1,-1,-2, 3)): 0.0,
+
+    ((1,-1,-2, 3), (1,-1, 2,-2)): 0.0,
+    ((1,-1,-2, 3), (-1,1, 2,-2)): 0.0,
+    ((1,-1,-2, 3), (1,-1, 2,-3)): 0.0,
+    ((1,-1,-2, 3), (1,-1,-2, 3)): -0.8115380000000001,
+}
+
+@pytest.mark.parametrize("inp1", basis)
+@pytest.mark.parametrize("inp2", basis)
+def test_overlap(inp1, inp2):
+
+    key = (tuple(inp1), tuple(inp2))
+    x_expected = expected[key]
+
+    # Build the S matrix
+    S = np.eye(6)
+    S[0,0] = 0.91
+    S[1,1] = -0.98
+    S[3,3] = 0.91
+    S[4,4] = -0.98
+
+    x = ovlp_arb(inp1, inp2, S, None, False)
+
+    # Numerical tolerance
+    assert np.isclose(x, x_expected, atol=1e-12)
+