Merge pull request #118 from giovanni-br/issue-117-fix

Fixes #117: Cleaning rauk module

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": {},

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
+
 
 from typing import Union
 
@@ -30,9 +33,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
@@ -43,15 +43,9 @@ 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])
+        self.n_sites = max_site
 
         if self.atom_types is None:
             # Initialize atom_types with None, and adjust size for 0-based
@@ -63,7 +57,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))

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,25 +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_types,
-                self.atoms_num,
-                self.n_sites,
                 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")
@@ -283,6 +280,7 @@ def __init__(
             atom_types=None,
             atom_dictionary=None,
             bond_dictionary=None,
+            orbital_overlap=None,
             Bz=None,
             gamma=None,
     ):
@@ -327,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
         )

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,28 +204,22 @@ 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.
 
     Parameters
     ----------
     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,24 +231,61 @@ 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,
         atom_dictionary,
-        n_sites,
         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