From ee4ade310c5c97ac15add8b1e4e14b027748ecce Mon Sep 17 00:00:00 2001
From: giovanni-br <giovannibenedettidarosa@gmail.com>
Date: Sat, 1 Jun 2024 16:13:12 +0200
Subject: [PATCH 1/2] Fixes #117: Cleaning rauk module

---
 moha/api.py              | 19 ++++-------
 moha/hamiltonians.py     |  7 ++--
 moha/rauk/PariserParr.py | 25 ++++++--------
 moha/rauk/rauk.py        | 48 ++++++++++++---------------
 moha/rauk/utils.py       | 72 ++++++++++++++++++++++++++++++++++++++++
 moha/utils.py            | 45 -------------------------
 6 files changed, 112 insertions(+), 104 deletions(-)
 create mode 100644 moha/rauk/utils.py

diff --git a/moha/api.py b/moha/api.py
index f3244a8..5555e90 100644
--- a/moha/api.py
+++ b/moha/api.py
@@ -8,7 +8,10 @@
 
 from scipy.sparse import csr_matrix, lil_matrix, diags
 
-from .utils import convert_indices, get_atom_type
+from .utils import convert_indices
+
+from moha.rauk.utils import get_atom_type, get_atoms_list
+
 
 __all__ = [
     "HamiltonianAPI",
@@ -28,9 +31,6 @@ def generate_connectivity_matrix(self):
         tuple
             (dictionary, np.ndarray)
         """
-        max_site = 0
-        atoms_sites_lst = []
-
         # check if self.connectivity is a matrix
         # if so, put assign it to self.connectivity_matrix
         # and set the atom_types to None
@@ -41,15 +41,8 @@ def generate_connectivity_matrix(self):
 
             return None, self.connectivity_matrix
 
-        for atom1, atom2, bond in self.connectivity:
-            atom1_name, site1 = get_atom_type(atom1)
-            atom2_name, site2 = get_atom_type(atom2)
-            for pair in [(atom1_name, site1), (atom2_name, site2)]:
-                if pair not in atoms_sites_lst:
-                    atoms_sites_lst.append(pair)
-            if max_site < max(site1, site2):  # finding the max index of site
-                max_site = max(site1, site2)
-        self.n_sites = len(atoms_sites_lst)
+        atoms_sites_lst = get_atoms_list(self.connectivity)
+        max_site = max([site for _, site in atoms_sites_lst])
 
         if self.atom_types is None:
             # Initialize atom_types with None, and adjust size for 0-based
diff --git a/moha/hamiltonians.py b/moha/hamiltonians.py
index 6c7cd0c..c4d6a00 100644
--- a/moha/hamiltonians.py
+++ b/moha/hamiltonians.py
@@ -137,8 +137,8 @@ def generate_one_body_integral(self, basis: str, dense: bool):
                 diags([self.alpha for _ in range(self.n_sites)], format="csr")
                 + self.beta * self.connectivity_matrix
             )
-        # check if alpha and beta are different from the default or a carbon
-        # chain
+        # check if alpha and beta are different from the default or
+        # all atom types are the same
         elif (
             self.alpha != -0.414 and self.beta != -0.0533
         ) or len(np.unique(self.atom_types)) == 1:
@@ -150,9 +150,6 @@ def generate_one_body_integral(self, basis: str, dense: bool):
         elif np.all([isinstance(k, int) for _, _, k in self.connectivity]):
             one_body_term = assign_rauk_parameters(
                 self.connectivity,
-                self.atom_types,
-                self.atoms_num,
-                self.n_sites,
                 self.atom_dictionary,
                 self.bond_dictionary
             )
diff --git a/moha/rauk/PariserParr.py b/moha/rauk/PariserParr.py
index 943f3a6..0926389 100644
--- a/moha/rauk/PariserParr.py
+++ b/moha/rauk/PariserParr.py
@@ -17,7 +17,7 @@
 import json
 
 
-def generate_alpha_beta(distance, atom1_name, atom2_name, ionization, ev_H):
+def populate_PP_dct(distance, atom1_name, atom2_name, ionization):
     r"""
     Calculate the beta for a bond based on atomic ionization and distance.
 
@@ -31,20 +31,19 @@ def generate_alpha_beta(distance, atom1_name, atom2_name, ionization, ev_H):
         Name of the second atom.
     ionization : dict
         Dictionary containing ionization energies of atoms.
-    ev_H : float
-        Conversion factor from electron volts to hartrees.
 
     Returns
     -------
     float
         Calculated beta parameter for the bond.
     """
+    ev_H = constants.value('electron volt-hartree relationship')
     alpha_x = float(-ionization[atom1_name]) * ev_H
     alpha_y = float(-ionization[atom2_name]) * ev_H
-    Rxy = distance
-    p = -(((alpha_x) + (alpha_y)) * Rxy) / 2
+    Rxy = float(distance)
+    p = - 0.5 * Rxy * (alpha_x + alpha_y)
     t = abs((alpha_x - alpha_y) / (alpha_x + alpha_y))
-    beta_xy = 1.75 * (Sxy(t, p)) * ((alpha_x + alpha_y) / 2)
+    beta_xy = 1.75 * Sxy(t, p) * (alpha_x + alpha_y) * 0.5
     return beta_xy
 
 
@@ -114,8 +113,7 @@ def Bn(n, t, p):
     """
     if t == 0:
         return 2 / (n + 1)
-    else:
-        return -np.exp(-p * t) * an(n, p * t) - np.exp(p * t) * bn(n, p * t)
+    return -np.exp(-p * t) * an(n, p * t) - np.exp(p * t) * bn(n, p * t)
 
 
 def An(n, p):
@@ -202,22 +200,19 @@ def compute_param_dist_overlap(
             if atom not in atom_dictionary.keys():
                 atom_dictionary[atom] = -ionization[atom] * ev_H
     if bond_dictionary is None:
-        ev_H = constants.value('electron volt-hartree relationship')
         ionization_path = Path(__file__).parent / "ionization.json"
         ionization = json.load(open(ionization_path, "rb"))
         bond_dictionary = {}
-        for trip in atoms_dist:
-            beta_xy = generate_alpha_beta(
-                trip[2], trip[0], trip[1], ionization, ev_H)
-            bond_key_forward = ','.join([trip[0], trip[1]])
-            bond_key_reverse = ','.join([trip[1], trip[0]])
+        for atom1, atom2, dist in atoms_dist:
+            beta_xy = populate_PP_dct(dist, atom1, atom2, ionization)
+            bond_key_forward = ','.join([atom1, atom2])
+            bond_key_reverse = ','.join([atom2, atom1])
             bond_dictionary[bond_key_forward] = beta_xy
             bond_dictionary[bond_key_reverse] = beta_xy
 
     one_body = build_one_body(
         connectivity,
         atom_dictionary,
-        n_sites,
         bond_dictionary)
 
     return one_body
diff --git a/moha/rauk/rauk.py b/moha/rauk/rauk.py
index f7b71f7..a90fc44 100644
--- a/moha/rauk/rauk.py
+++ b/moha/rauk/rauk.py
@@ -1,7 +1,7 @@
 r"""Rauk Module."""
 import numpy as np
 
-from moha.utils import get_atom_type
+from moha.rauk.utils import get_atom_type
 
 from pathlib import Path
 
@@ -13,28 +13,31 @@
 def build_one_body(
         connectivity,
         atom_dictionary,
-        n_sites,
         bond_dictionary):
     r"""
     Construct the one-body matrix for a compound.
 
-    Parameters:
-    - connectivity (list of tuples): List of tuples where each tuple
-    represents a bondbetween two atoms.
-    - atom_dictionary (dict): Dictionary mapping atom names and properties
-    - atoms_num (list): List of tuples (atom_name, quantity).
-    - n_sites (int): Total number of sites (atoms) in the molecule.
-    - bond_dictionary (dict): Dictionary mapping pairs of atom names
-    to properties.
-
-    Returns:
-    - scipy.sparse.csr_matrix: one-body matrix.
+    Parameters
+    ----------
+    connectivity : list of tuples
+        List of tuples where each tuple represents a bond between two atoms.
+    atom_dictionary : dict
+        Dictionary mapping atom names to properties.
+    bond_dictionary : dict
+        Dictionary mapping pairs of atom names to properties.
+
+    Returns
+    -------
+    scipy.sparse.csr_matrix
+        One-body matrix.
     """
     # Populate diagonal and non-diagonal matrix elements
     # Create a sparse diagonal matrix
     diagonal_values = []
     atom_indices = {}
 
+    n_sites = len(atom_dictionary)
+
     # Initialize matrices and helper dictionaries
     diagonal_values = np.zeros(n_sites)
     param_nodiag_mtrx = np.zeros((n_sites, n_sites))
@@ -73,9 +76,6 @@ def build_one_body(
 
 def assign_rauk_parameters(
         connectivity,
-        atom_types,
-        atoms_num,
-        n_sites,
         atom_dictionary,
         bond_dictionary):
     r"""
@@ -88,12 +88,6 @@ def assign_rauk_parameters(
     ----------
     connectivity : list
         Connections between atoms.
-    atom_types : list
-        Atom types in the molecule.
-    atoms_num : list
-        Tuples of number and types of atoms.
-    n_sites : int
-        Total number of atomic sites.
     atom_dictionary : dict, optional
         Atom parameters dictionary. If None, loaded from JSON.
     bond_dictionary : dict, optional
@@ -113,8 +107,11 @@ def assign_rauk_parameters(
         beta_c = -0.0533  # Value for sp2 orbitals of Carbon atom.
         # Create atom dictionary using predefined values without overlap
         # parameters
-        for atom in atom_types:
-            atom_dictionary[atom] = alpha_c + hx_dictionary[atom] * abs(beta_c)
+        for atom1, atom2, _ in connectivity:
+            atom1_name, site1 = get_atom_type(atom1)
+            atom2_name, site2 = get_atom_type(atom2)
+            hx_value = hx_dictionary[atom1_name] * abs(beta_c)
+            atom_dictionary[atom1_name] = alpha_c + hx_value
     if bond_dictionary is None:
         # Paths to the JSON files
         hx_dictionary_path = Path(__file__).parent / "hx_dictionary.json"
@@ -143,13 +140,12 @@ def assign_rauk_parameters(
             bond_key = ','.join([atom1_name, atom2_name])
             bond_value = kxy_matrix[index1, index2] * abs(beta_c)
             bond_dictionary[bond_key] = bond_value
-            bond_dictionary[','.join([atom2_name, atom1_name])] = bond_value
             # Ensure symmetry
+            bond_dictionary[','.join([atom2_name, atom1_name])] = bond_value
 
     one_body = build_one_body(
         connectivity,
         atom_dictionary,
-        n_sites,
         bond_dictionary)
 
     return one_body
diff --git a/moha/rauk/utils.py b/moha/rauk/utils.py
new file mode 100644
index 0000000..eb0c13f
--- /dev/null
+++ b/moha/rauk/utils.py
@@ -0,0 +1,72 @@
+r"""Utils for Rauk Module."""
+
+import re
+
+
+def get_atom_type(atom):
+    r"""
+    Extract the atom type and position index from the atom string.
+
+    Parameters
+    ----------
+    atom : str
+        String representing an atom, following
+        the pattern "typeposition(site)".
+
+    Returns
+    -------
+    tuple
+        A tuple containing the atom type (including site index) as a string and
+        the position index as an integer.
+    """
+    # The pattern matches an initial letter sequence for the atom type,
+    # followed by a number for the position, and numbers in parentheses for
+    # site index to be appended.
+    pattern = r"([A-Za-z]+)(\d+)\((\d+)\)"
+    match = re.match(pattern, atom)
+
+    if match:
+        # If the pattern matches, append the site index to the atom type
+        atom_type = match.group(1) + match.group(3)
+        position_index = int(match.group(2))
+    else:
+        # Fallback for handling cases without parentheses
+        i = 1
+        while i <= len(atom) and atom[-i].isdigit():
+            i += 1
+        i -= 1
+        if i == 0:
+            raise ValueError(f"Invalid atom format: {atom}")
+        atom_type = atom[:-i]
+        position_index = int(atom[-i:])
+
+    return atom_type, position_index
+
+
+def get_atoms_list(connectivity):
+    r"""
+    Process the connectivity information of a molecular compound.
+
+    Parameters
+    ----------
+    connectivity : list of tuples
+        List of tuples representing bonds between atoms (atom1, atom2, order).
+
+    Returns
+    -------
+    int
+        Total number of unique atomic sites in the molecule.
+    """
+    atoms_sites_lst = []
+    max_site = 0
+
+    for atom1, atom2, _ in connectivity:
+        atom1_name, site1 = get_atom_type(atom1)
+        atom2_name, site2 = get_atom_type(atom2)
+        for pair in [(atom1_name, site1), (atom2_name, site2)]:
+            if pair not in atoms_sites_lst:
+                atoms_sites_lst.append(pair)
+        if max_site < max(site1, site2):  # finding the max index of site
+            max_site = max(site1, site2)
+
+    return atoms_sites_lst
diff --git a/moha/utils.py b/moha/utils.py
index 7d0b33b..d46a936 100644
--- a/moha/utils.py
+++ b/moha/utils.py
@@ -13,51 +13,6 @@
 ]
 
 
-def get_atom_type(atom):
-    r"""
-    Construct the one-body matrix for a molecular compound.
-
-    Parameters
-    ----------
-    connectivity : list of tuples
-        List of tuples representing bonds between atoms (atom1, atom2, order).
-    atom_dictionary : dict
-        Dictionary mapping atom types to properties for  matrix elements.
-    n_sites : int
-        Total number of unique atomic sites in the molecule.
-    bond_dictionary : dict
-        Dictionary mapping 'type1,type2' pairs to properties for off-diagonal
-        elements.
-
-    Returns
-    -------
-    scipy.sparse.csr_matrix
-        Compressed sparse
-    """
-    # The pattern matches an initial letter sequence for the atom type,
-    # followed by a number for the position, and numbers in parentheses for
-    # site index to be appended.
-    pattern = r"([A-Za-z]+)(\d+)\((\d+)\)"
-    match = re.match(pattern, atom)
-
-    if match:
-        # If the pattern matches, append the site index to the atom type
-        atom_type = match.group(1) + match.group(3)
-        position_index = int(match.group(2))
-    else:
-        # Fallback for handling cases without parentheses
-        i = 1
-        while i <= len(atom) and atom[-i].isdigit():
-            i += 1
-        i -= 1
-        if i == 0:
-            raise ValueError(f"Invalid atom format: {atom}")
-        atom_type = atom[:-i]
-        position_index = int(atom[-i:])
-
-    return atom_type, position_index
-
-
 def convert_indices(N, *args):
     r"""
     Convert indices from 4d array to 2d numpy array and vice-versa.

From e4a794747ae5938bed3d62a8b8f242fbf279de31 Mon Sep 17 00:00:00 2001
From: giovanni-br <giovannibenedettidarosa@gmail.com>
Date: Tue, 4 Jun 2024 18:18:35 +0200
Subject: [PATCH 2/2] #118 Fixes #117: Cleaning rauk module

---
 examples/rauk.ipynb      |  50 ++++++++++++++--
 moha/api.py              |   4 +-
 moha/hamiltonians.py     |  22 +++----
 moha/rauk/PariserParr.py | 123 +++++++++++++++++++++++++++++++--------
 moha/rauk/rauk.py        |   7 ++-
 moha/rauk/utils.py       |  15 +++--
 6 files changed, 172 insertions(+), 49 deletions(-)

diff --git a/examples/rauk.ipynb b/examples/rauk.ipynb
index 71afbe2..31d1c1d 100644
--- a/examples/rauk.ipynb
+++ b/examples/rauk.ipynb
@@ -99,10 +99,10 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "[[-0.41380839 -0.73648465  0.         -0.56460276]\n",
-      " [-0.73648465 -0.47655051 -0.84086439  0.        ]\n",
-      " [ 0.         -0.84086439 -0.64027609 -0.54414171]\n",
-      " [-0.56460276  0.         -0.54414171 -0.29956945]]\n"
+      "[[-0.41380839 -0.73648465  0.         -0.56689357]\n",
+      " [-0.73648465 -0.47655051 -0.81351185  0.        ]\n",
+      " [ 0.         -0.81351185 -0.64027609 -0.60447297]\n",
+      " [-0.56689357  0.         -0.60447297 -0.29956945]]\n"
      ]
     }
    ],
@@ -113,8 +113,8 @@
     "from moha import HamHub\n",
     "\n",
     "\n",
-    "system = [('C1', 'Cl2', 1.5), ('Cl2', 'F3', 1.8),\n",
-    "                    ('F3', 'Si4', 1.8), ('Si4', 'C1', 1.7)]\n",
+    "system = [('C1', 'Cl2', 1.5, 'pi'), ('Cl2', 'F3', 1.8, 'sigma'),\n",
+    "                    ('F3', 'Si4', 1.8, 'sigma'), ('Si4', 'C1', 1.7, 'sigma')]\n",
     "hubbard = HamHub(system, u_onsite=np.array([1, 1]))\n",
     "\n",
     "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))"
@@ -204,6 +204,44 @@
     "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[[-0.41380839 -0.77906404  0.         -1.87261684]\n",
+      " [-0.77906404 -0.47655051 -1.95444655  0.        ]\n",
+      " [ 0.         -1.95444655 -0.64027609 -1.64472969]\n",
+      " [-1.87261684  0.         -1.64472969 -0.29956945]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "import sys  \n",
+    "sys.path.insert(0, '../')\n",
+    "import numpy as np\n",
+    "from moha import HamHub\n",
+    "\n",
+    "# First way to define the Hamiltonian\n",
+    "# two site Hubbard model\n",
+    "\n",
+    "# system = [('C1', 'C2', 1)] is a list of tuples, where each tuple represents a bond\n",
+    "# between two atoms and the third element is the type of bond (singe or double).\n",
+    "# For now, we only support single bonds between carbon atoms. \n",
+    "# For this type of bonds the default values of alpha and beta are -0.414 and -0.0533, respectively.\n",
+    "# In the future we are planning to support different types of bonds for different atoms.\n",
+    "system = [('C1', 'Cl2', 1.1), ('Cl2', 'F3', 1.0),\n",
+    "                    ('F3', 'Si4', 1.0), ('Si4', 'C1', 1.0)]\n",
+    "hubbard = HamHub(system, u_onsite=np.array([1, 1]),\n",
+    "                 orbital_overlap={'C,Cl': 1, 'Cl,F': 2, 'F,Si': 2, 'Si,C': 3})\n",
+    "\n",
+    "print(hubbard.generate_one_body_integral(dense=True, basis='spatial basis'))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
diff --git a/moha/api.py b/moha/api.py
index 5555e90..2b0b9ef 100644
--- a/moha/api.py
+++ b/moha/api.py
@@ -43,6 +43,7 @@ def generate_connectivity_matrix(self):
 
         atoms_sites_lst = get_atoms_list(self.connectivity)
         max_site = max([site for _, site in atoms_sites_lst])
+        self.n_sites = max_site
 
         if self.atom_types is None:
             # Initialize atom_types with None, and adjust size for 0-based
@@ -54,7 +55,8 @@ def generate_connectivity_matrix(self):
             self.atom_types = atom_types
         connectivity_mtrx = np.zeros((max_site, max_site))
         atoms_dist = []
-        for atom1, atom2, bond in self.connectivity:
+        for tpl in self.connectivity:
+            atom1, atom2, bond = tpl[0], tpl[1], tpl[2]
             atom1_name, site1 = get_atom_type(atom1)
             atom2_name, site2 = get_atom_type(atom2)
             atoms_dist.append((atom1_name, atom2_name, bond))
diff --git a/moha/hamiltonians.py b/moha/hamiltonians.py
index c4d6a00..21c1227 100644
--- a/moha/hamiltonians.py
+++ b/moha/hamiltonians.py
@@ -12,7 +12,7 @@
 from typing import Union
 
 from moha.rauk.rauk import assign_rauk_parameters
-from moha.rauk.PariserParr import compute_param_dist_overlap
+from moha.rauk.PariserParr import compute_overlap
 import warnings
 
 warnings.simplefilter('ignore',
@@ -39,7 +39,8 @@ def __init__(
             charges=None,
             sym=1,
             atom_dictionary=None,
-            bond_dictionary=None
+            bond_dictionary=None,
+            orbital_overlap=None
     ):
         r"""
         Initialize Pariser-Parr-Pople Hamiltonian.
@@ -101,6 +102,7 @@ def __init__(
         self.two_body = None
         self.bond_dictionary = bond_dictionary
         self.atom_dictionary = atom_dictionary
+        self.orbital_overlap = orbital_overlap
 
     def generate_zero_body_integral(self):
         r"""Generate zero body integral.
@@ -147,22 +149,20 @@ def generate_one_body_integral(self, basis: str, dense: bool):
                 + self.beta * self.connectivity_matrix
             )
         # check if elements in connectivity matrix are integer
-        elif np.all([isinstance(k, int) for _, _, k in self.connectivity]):
+        elif np.all([isinstance(elem[2], int) for elem in self.connectivity]):
             one_body_term = assign_rauk_parameters(
                 self.connectivity,
                 self.atom_dictionary,
                 self.bond_dictionary
             )
         # check if elements in connectivity matrix are float
-        elif np.all([isinstance(k, float) for _, _, k in self.connectivity]):
-            one_body_term = compute_param_dist_overlap(
+        elif np.all([isinstance(elem[2], float)
+                    for elem in self.connectivity]):
+            one_body_term = compute_overlap(
                 self.connectivity,
-                self.atom_types,
-                self.atoms_num,
-                self.n_sites,
-                self.atoms_dist,
                 self.atom_dictionary,
-                self.bond_dictionary
+                self.bond_dictionary,
+                self.orbital_overlap
             )
         else:
             raise TypeError("Connectivity matrix has wrong type")
@@ -280,6 +280,7 @@ def __init__(
             atom_types=None,
             atom_dictionary=None,
             bond_dictionary=None,
+            orbital_overlap=None,
             Bz=None,
             gamma=None,
     ):
@@ -324,6 +325,7 @@ def __init__(
             gamma=None,
             atom_dictionary=atom_dictionary,
             bond_dictionary=bond_dictionary,
+            orbital_overlap=orbital_overlap,
             charges=np.array(0),
             sym=sym
         )
diff --git a/moha/rauk/PariserParr.py b/moha/rauk/PariserParr.py
index 0926389..8d5d491 100644
--- a/moha/rauk/PariserParr.py
+++ b/moha/rauk/PariserParr.py
@@ -11,13 +11,21 @@
 
 from moha.rauk.rauk import build_one_body
 
+from moha.rauk.utils import get_atom_type
+
 import numpy as np
 from scipy.special import gamma
 from pathlib import Path
 import json
 
 
-def populate_PP_dct(distance, atom1_name, atom2_name, ionization):
+def populate_PP_dct(
+        distance,
+        atom1_name,
+        atom2_name,
+        ionization,
+        Sxy=None,
+        bond_type='sigma'):
     r"""
     Calculate the beta for a bond based on atomic ionization and distance.
 
@@ -40,10 +48,17 @@ def populate_PP_dct(distance, atom1_name, atom2_name, ionization):
     ev_H = constants.value('electron volt-hartree relationship')
     alpha_x = float(-ionization[atom1_name]) * ev_H
     alpha_y = float(-ionization[atom2_name]) * ev_H
-    Rxy = float(distance)
-    p = - 0.5 * Rxy * (alpha_x + alpha_y)
-    t = abs((alpha_x - alpha_y) / (alpha_x + alpha_y))
-    beta_xy = 1.75 * Sxy(t, p) * (alpha_x + alpha_y) * 0.5
+    if Sxy is None:
+        Rxy = float(distance)
+        p = - 0.5 * Rxy * (alpha_x + alpha_y)
+        t = abs((alpha_x - alpha_y) / (alpha_x + alpha_y))
+        if bond_type == 'sigma':
+            Sxy = Sxy_sigma(t, p)
+        elif bond_type == 'pi':
+            Sxy = Sxy_pi(t, p)
+        else:
+            raise ValueError(f"Invalid bond type: {bond_type}")
+    beta_xy = 1.75 * Sxy * (alpha_x + alpha_y) * 0.5
     return beta_xy
 
 
@@ -135,9 +150,35 @@ def An(n, p):
     return np.exp(-p) * an(n, p)
 
 
-def Sxy(t, p):
+def Sxy_sigma(t, p):
     r"""
-    Calculate the overlap integral Sxy for given parameters t and p.
+    Calculate the overlap integral Sxy for a sigma bond.
+
+    Parameters
+    ----------
+    t : float
+        Parameter t in the formula.
+    p : float
+        Parameter p in the formula.
+
+    Returns
+    -------
+    float
+        Calculated Sxy value.
+    """
+    if t == 0:
+        return np.exp(-p) * (1 + p + (1 / 3) * p**2)
+    else:
+        A2 = An(2, p)
+        A0 = An(0, p)
+        B0 = Bn(0, t, p)
+        B2 = Bn(2, t, p)
+        return ((1 - t**2)**(3 / 2) * p**3) * (A2 * B0 - A0 * B2) / 4
+
+
+def Sxy_pi(t, p):
+    r"""
+    Calculate the overlap integral Sxy for a pi bond.
 
     Parameters
     ----------
@@ -163,9 +204,9 @@ def Sxy(t, p):
                ((1 - t**2)**(5 / 2)) * (p**5) / 32
 
 
-def compute_param_dist_overlap(
-        connectivity, atom_types, atoms_num, n_sites, atoms_dist,
-        atom_dictionary, bond_dictionary):
+def compute_overlap(
+        connectivity, atom_dictionary,
+        bond_dictionary, orbital_overlap):
     r"""
     Compute the parameterized distance overlap matrix for a set of atoms.
 
@@ -173,18 +214,12 @@ def compute_param_dist_overlap(
     ----------
     connectivity : list of tuples
         List defining connectivity between atoms (atom1, atom2, order).
-    atom_types : list
-        List of atom types present in the molecule.
-    atoms_num : list
-        List of tuples defining atom types and their quantities.
-    n_sites : int
-        Number of sites (atoms) in the molecule.
-    atoms_dist : list
-        List defining distances between connected atoms.
     atom_dictionary : dict
         Dictionary mapping atom types to properties.
     bond_dictionary : dict
         Dictionary mapping pairs of atom types to bond properties.
+    orbital_overlap : dict
+        Dictionary mapping pairs of atom types to orbital overlap properties.
 
     Returns
     -------
@@ -196,19 +231,38 @@ def compute_param_dist_overlap(
         ionization = json.load(open(ionization_path, "rb"))
         atom_dictionary = {}  # alpha as first ionization potential
         ev_H = constants.value('electron volt-hartree relationship')
-        for atom in atom_types:
+        for tpl in connectivity:
+            atom, _ = get_atom_type(tpl[0])
             if atom not in atom_dictionary.keys():
                 atom_dictionary[atom] = -ionization[atom] * ev_H
     if bond_dictionary is None:
         ionization_path = Path(__file__).parent / "ionization.json"
         ionization = json.load(open(ionization_path, "rb"))
         bond_dictionary = {}
-        for atom1, atom2, dist in atoms_dist:
-            beta_xy = populate_PP_dct(dist, atom1, atom2, ionization)
-            bond_key_forward = ','.join([atom1, atom2])
-            bond_key_reverse = ','.join([atom2, atom1])
-            bond_dictionary[bond_key_forward] = beta_xy
-            bond_dictionary[bond_key_reverse] = beta_xy
+        if orbital_overlap is None:
+            for atom1, atom2, dist, bond_type in connectivity:
+                atom1_name, _ = get_atom_type(atom1)
+                atom2_name, _ = get_atom_type(atom2)
+                bond_key_forward = ','.join([atom1_name, atom2_name])
+                bond_key_reverse = ','.join([atom2_name, atom1_name])
+                beta_xy = populate_PP_dct(
+                    dist, atom1_name, atom2_name, ionization,
+                    bond_type=bond_type)
+                bond_dictionary[bond_key_forward] = beta_xy
+                bond_dictionary[bond_key_reverse] = beta_xy
+        else:
+            for tpl in connectivity:
+                atom1, atom2, dist = tpl[0], tpl[1], tpl[2]
+                atom1_name, _ = get_atom_type(atom1)
+                atom2_name, _ = get_atom_type(atom2)
+                bond_key_forward = ','.join([atom1_name, atom2_name])
+                bond_key_reverse = ','.join([atom2_name, atom1_name])
+                Sxy = orbital_overlap[bond_key_forward]
+
+                beta_xy = populate_PP_dct(
+                    dist, atom1_name, atom2_name, ionization, Sxy)
+                bond_dictionary[bond_key_forward] = beta_xy
+                bond_dictionary[bond_key_reverse] = beta_xy
 
     one_body = build_one_body(
         connectivity,
@@ -216,3 +270,22 @@ def compute_param_dist_overlap(
         bond_dictionary)
 
     return one_body
+
+
+def calculate_gamma(Uxy_bar, Rxy):
+    """
+    Calculate the gamma value based on Uxy and Rxy.
+
+    Parameters
+    ----------
+    Uxy_bar (float): Represents the potential energy
+    Rxy (float): Represents the distance or a related measure.
+
+    Returns
+    ----------
+    float: Computed gamma value based on the given parameters.
+    """
+    # Example formula, needs actual formula to be replaced here
+    # This is just a placeholder formula
+    gamma = Uxy_bar / (Uxy_bar * Rxy + np.exp(-1 / 2 * Uxy_bar**2 * Rxy ** 2))
+    return gamma
diff --git a/moha/rauk/rauk.py b/moha/rauk/rauk.py
index a90fc44..7d36810 100644
--- a/moha/rauk/rauk.py
+++ b/moha/rauk/rauk.py
@@ -1,7 +1,7 @@
 r"""Rauk Module."""
 import numpy as np
 
-from moha.rauk.utils import get_atom_type
+from moha.rauk.utils import get_atom_type, get_atoms_list
 
 from pathlib import Path
 
@@ -36,7 +36,7 @@ def build_one_body(
     diagonal_values = []
     atom_indices = {}
 
-    n_sites = len(atom_dictionary)
+    _, n_sites = get_atoms_list(connectivity, return_nsites=True)
 
     # Initialize matrices and helper dictionaries
     diagonal_values = np.zeros(n_sites)
@@ -44,7 +44,8 @@ def build_one_body(
     atom_indices = {}
 
     # Populate diagonal and non-diagonal matrix elements
-    for atom1, atom2, _ in connectivity:
+    for tpl in connectivity:
+        atom1, atom2 = tpl[0], tpl[1]
         for atom in [atom1, atom2]:
             atom_type, index = get_atom_type(atom)
             if index not in atom_indices:
diff --git a/moha/rauk/utils.py b/moha/rauk/utils.py
index eb0c13f..9708f8d 100644
--- a/moha/rauk/utils.py
+++ b/moha/rauk/utils.py
@@ -43,7 +43,7 @@ def get_atom_type(atom):
     return atom_type, position_index
 
 
-def get_atoms_list(connectivity):
+def get_atoms_list(connectivity, return_nsites=False):
     r"""
     Process the connectivity information of a molecular compound.
 
@@ -51,16 +51,21 @@ def get_atoms_list(connectivity):
     ----------
     connectivity : list of tuples
         List of tuples representing bonds between atoms (atom1, atom2, order).
+    return_nsites : bool, optional
+        Whether to return total number of sites in the system
 
     Returns
     -------
-    int
-        Total number of unique atomic sites in the molecule.
+    atoms_sites_lst: list
+        List with atom types and site indices.
+    max_site : int
+        Maximum site index in the system.
     """
     atoms_sites_lst = []
     max_site = 0
 
-    for atom1, atom2, _ in connectivity:
+    for tpl in connectivity:
+        atom1, atom2 = tpl[0], tpl[1]
         atom1_name, site1 = get_atom_type(atom1)
         atom2_name, site2 = get_atom_type(atom2)
         for pair in [(atom1_name, site1), (atom2_name, site2)]:
@@ -68,5 +73,7 @@ def get_atoms_list(connectivity):
                 atoms_sites_lst.append(pair)
         if max_site < max(site1, site2):  # finding the max index of site
             max_site = max(site1, site2)
+    if return_nsites:
+        return atoms_sites_lst, max_site
 
     return atoms_sites_lst