fixing errors

examples/rauk.ipynb

@@ -40,7 +40,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -92,7 +92,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -131,7 +131,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
@@ -167,7 +167,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
@@ -206,7 +206,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [
     {
@@ -251,7 +251,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [
     {
@@ -321,7 +321,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [
     {

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_overlap, compute_u, compute_gamma
+from moha.rauk.PariserParr import compute_overlap, compute_gamma, compute_u
 import warnings
 
 warnings.simplefilter('ignore',
@@ -40,7 +40,9 @@ def __init__(
             sym=1,
             atom_dictionary=None,
             bond_dictionary=None,
-            orbital_overlap=None
+            orbital_overlap=None,
+            affinity_dct=None,
+            Rxy_list=None
     ):
         r"""
         Initialize Pariser-Parr-Pople Hamiltonian.
@@ -103,6 +105,8 @@ def __init__(
         self.bond_dictionary = bond_dictionary
         self.atom_dictionary = atom_dictionary
         self.orbital_overlap = orbital_overlap
+        self.affinity_dct = affinity_dct
+        self.Rxy_list = Rxy_list
 
     def generate_zero_body_integral(self):
         r"""Generate zero body integral.
@@ -701,4 +705,5 @@ def __init__(self,
             mu=mu,
             J_eq=J_eq,
             J_ax=J_ax,
-            connectivity=connectivity)
+            connectivity=connectivity
+        )

moha/rauk/PariserParr.py

@@ -274,23 +274,43 @@ def compute_overlap(
     return one_body
 
 
-def compute_gamma(Uxy_bar, Rxy):
-    """
-    Calculate the gamma value based on Uxy and Rxy.
+def compute_gamma(u_onsite, Rxy_list, connectivity):
+    r"""
+    Calculate the gamma values for each pair of sites.
 
     Parameters
     ----------
-    Uxy_bar (float): Represents the potential energy
-    Rxy (float): Represents the distance or a related measure.
+    u_onsite (list of float): List of potential energies for sites.
+    Rxy_list (list of float): List of distances
+    connectivity (list of tuples): Each tuple contains indices of two
+                                   sites (atom1, atom2), indicating that a
+                                   gamma value should be computed for
+                                   this pair.
 
     Returns
     ----------
-    float: Computed gamma value based on the given parameters.
+    np.ndarray: A matrix of computed gamma values for each pair.
+                Non-connected pairs have a gamma value of zero.
     """
-    # 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
+    num_sites = len(u_onsite)
+    gamma_matrix = np.zeros((num_sites, num_sites))
+
+    for tpl in connectivity:
+        atom1, atom2 = tpl[0], tpl[1]
+        atom1_name, site1 = get_atom_type(atom1)
+        atom2_name, site2 = get_atom_type(atom2)
+
+        if site1 < num_sites and site2 < num_sites:
+            Ux = u_onsite[site1]
+            Uy = u_onsite[site2]
+            Rxy = Rxy_list[site1]
+            Uxy_bar = 0.5 * (Ux + Uy)
+            gamma = Uxy_bar / \
+                (Uxy_bar * Rxy + np.exp(-0.5 * Uxy_bar**2 * Rxy**2))
+            gamma_matrix[site1][site2] = gamma
+
+    return gamma_matrix
+
 
 def compute_u(connectivity, atom_dictionary, affinity_dictionary):
     r"""