From 68d52238cdf98f8e030d2fda4269387729cc7351 Mon Sep 17 00:00:00 2001
From: Giovanni Benedetti <86328308+giovanni-br@users.noreply.github.com>
Date: Sun, 27 Jul 2025 22:25:46 +0200
Subject: [PATCH 1/3] adding reduced basis

---
 examples/molkit.ipynb        | 232 ++++++++++++++++++++++++++++++++---
 moha/molkit/hamiltonians.py  | 162 +++++++++++++++---------
 moha/molkit/utils/spinops.py |  39 ++++++
 moha/molkit/utils/tools.py   |   1 +
 4 files changed, 360 insertions(+), 74 deletions(-)

diff --git a/examples/molkit.ipynb b/examples/molkit.ipynb
index 6a928d1..46331e7 100644
--- a/examples/molkit.ipynb
+++ b/examples/molkit.ipynb
@@ -32,17 +32,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Max error: 1.0\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
     "from moha import HamHeisenberg\n",
@@ -105,11 +97,223 @@
     "h1_spin, h2_spin = Molecular_Hamiltonian.spinize_H()\n",
     "h2_spin = antisymmetrize_two_body(h2_spin)\n",
     "\n",
-    "G = MolHam.to_geminal(h2_spin)   \n",
-    "n_orb = h2_spin.shape[0]\n",
+    "reduced = Molecular_Hamiltonian.to_reduced(n_elec=2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Energies Tests"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pyscf\n",
+    "from pyscf import gto, scf, fci\n",
+    "from pyscf.tools import fcidump\n",
+    "import os\n",
+    "\n",
+    "\n",
+    "path_xyz = r'../data/hydrogenic/H4/Paldus/h4_0.456_2.xyz'\n",
+    "mol = gto.M(atom=path_xyz, basis='sto-6g')\n",
+    "hf = scf.RHF(mol)\n",
+    "\n",
+    "hf.conv_tol = 1.0e-12\n",
+    "status = hf.run().converged\n",
+    "print(status)\n",
+    "\n",
+    "\n",
+    "mo_energy = hf.mo_energy\n",
+    "mo_occ = hf.mo_occ\n",
+    "\n",
+    "print('starting writing')\n",
+    "# Write FCIDUMP\n",
+    "\n",
+    "fcidump.from_mo(mol, 'tmp1.fcidump',\n",
+    "                hf.mo_coeff)\n",
+    "\n",
+    "# fcidump.from_scf(hf, \n",
+    "# f\"tmp.fcidump\",\n",
+    "# tol=0.0)\n",
+    "\n",
+    "print('done')\n",
+    "\n",
+    "data = fcidump.read(\"tmp1.fcidump\")\n",
+    "e, ci_vec = fci.direct_spin1.kernel(data['H1'],\n",
+    "                                    data['H2'], \n",
+    "                                    data['NORB'], \n",
+    "                                    data['NELEC'], ecore=data['ECORE'], nroots=1)\n",
+    "\n",
+    "rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, data['NORB'], data['NELEC'])\n",
+    "rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
+    "                 [np.zeros((4, 4)), rdm1s[1]]])\n",
+    "\n",
+    "rdm2 = np.zeros((8, 8, 8, 8))\n",
+    "for i in range(4):\n",
+    "    for j in range(4):\n",
+    "        for k in range(4):\n",
+    "            for l in range(4):\n",
+    "                rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
+    "                rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
+    "                rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
+    "                # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
+    "                rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]\n",
+    "                \n",
+    "\n",
+    "# converting rdm2 to physics notation\n",
+    "# rdm2 = np.einsum('ijkl->ikjl', rdm2)\n",
+    "print(\"trace of rdm1\", np.einsum('ii', rdm1))\n",
+    "print(\"Trace of rdm2\", np.einsum('ijij', rdm2))\n",
+    "\n",
+    "# loading integrals\n",
+    "# iodata stores integrals in physics notations\n",
+    "# https://iodata.readthedocs.io/en/latest/_modules/iodata/formats/fcidump.html#load_one\n",
+    "integrals = load_one(open('tmp1.fcidump', 'r'))\n",
+    "os.remove('tmp1.fcidump')\n",
+    "\n",
+    "h1s = integrals['one_ints']['core_mo']\n",
+    "h2s = integrals['two_ints']['two_mo']\n",
+    "\n",
+    "h1 = np.block([[h1s, np.zeros((4, 4))],\n",
+    "               [np.zeros((4, 4)), h1s]])\n",
+    "\n",
+    "h2 = np.zeros((8, 8, 8, 8))\n",
+    "for i in range(4):\n",
+    "    for j in range(4):\n",
+    "        for k in range(4):\n",
+    "            for l in range(4):\n",
+    "                # aaaa\n",
+    "                h2[i, j, k, l] = h2s[i, j, k, l]\n",
+    "                # bbbb\n",
+    "                h2[i+4, j+4, k+4, l+4] = h2s[i, j, k, l]\n",
+    "                # abab\n",
+    "                h2[i, j+4, k, l+4] = h2s[i, j, k, l]\n",
+    "                # baba\n",
+    "# h2[i+4, j, k+4, l] = h2s[k, l, i, j]\n",
+    "                h2[i+4, j, k+4, l] = h2s[i, j, k, l]\n",
+    "                \n",
+    "                \n",
+    "res = 0\n",
+    "for i in range(8):\n",
+    "    for j in range(8):\n",
+    "        assert rdm2[i, j, i, j] >=0\n",
+    "        res += rdm2[i, j, i, j]\n",
+    "        \n",
+    "print(\"trace of 2dm with for loop: \", res)\n",
+    "\n",
+    "e_dm = np.einsum('pq, pq', h1, rdm1) + 0.5 * np.einsum('pqrs, pqrs', h2, rdm2)\n",
+    "\n",
+    "\n",
+    "print(e_dm, e - integrals['core_energy'])\n",
+    "\n",
+    "\n",
+    "#-----------\n",
+    "eye = np.eye(8)\n",
+    "N = 4\n",
+    "N_alpha = N//2\n",
+    "\n",
+    "h1_ups = 0.5 * (np.kron(h1, eye) + np.kron(eye, h1)) / (N-1)\n",
+    "h1_ups = h1_ups.reshape(8, 8, 8, 8)\n",
+    "H = h1_ups + 0.5 * h2\n",
+    "\n",
+    "np.einsum('pqrs, pqrs', H, rdm2), e_dm\n",
+    "\n",
+    "# to geminal:\n",
+    "# H_gem = to_geminal(H)\n",
+    "# rdm_gem = to_geminal(rdm2)\n",
+    "# np.einsum('pq, pq', H_gem, rdm_gem) == np.einsum('pqrs, pqrs', H, rdm2) == e_dm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from moha.molkit.hamiltonians import MolHam     \n",
+    "h_sp = np.array([[1.0, 0.2],\n",
+    "                 [0.2, 1.5]])\n",
+    "\n",
+    "V_sp = np.zeros((2, 2, 2, 2))\n",
+    "V_sp[0, 0, 0, 0] = 0.7\n",
+    "V_sp[1, 1, 1, 1] = 0.6\n",
+    "\n",
+    "N_ELEC = 2\n",
+    "\n",
+    "toy_ham = MolHam(one_body=h_sp, two_body=V_sp)\n",
+    "\n",
+    "h_spin, V_spin = toy_ham.spinize_H()     \n",
+    "n_spin = toy_ham.n_spin      \n",
+    "\n",
+    "e_fci, ci_vec = fci.direct_spin1.kernel(\n",
+    "    h_spin, V_spin, n_spin, N_ELEC, ecore=0.0, nroots=1\n",
+    ")\n",
+    "\n",
+    "rdm1s, rdm2s = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)\n",
+    "rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
+    "                 [np.zeros((4, 4)), rdm1s[1]]])\n",
+    "\n",
+    "rdm2 = np.zeros((8, 8, 8, 8))\n",
+    "for i in range(4):\n",
+    "    for j in range(4):\n",
+    "        for k in range(4):\n",
+    "            for l in range(4):\n",
+    "                rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
+    "                rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
+    "                rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
+    "                # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
+    "                rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from pyscf import fci\n",
+    "\n",
+    "from moha.molkit.hamiltonians import MolHam\n",
+    "from moha.utils.spinops import (\n",
+    "    to_geminal_rdm2\n",
+    ")\n",
+    "\n",
+    "\n",
+    "h_sp = np.array([[1.0, 0.2],\n",
+    "                 [0.2, 1.5]])\n",
+    "\n",
+    "V_sp = np.zeros((2, 2, 2, 2))\n",
+    "V_sp[0, 0, 0, 0] = 0.7\n",
+    "V_sp[1, 1, 1, 1] = 0.6\n",
+    "\n",
+    "N_ELEC = 2                   \n",
+    "toy = MolHam(h_sp, V_sp)\n",
+    "\n",
+    "h_spin, V_spin = toy.spinize_H()           \n",
+    "\n",
+    "e_fci, ci_vec = fci.direct_spin1.kernel(h_spin, V_spin, n_spin, N_ELEC, ecore=0.0)\n",
+    "rdm1, rdm2 = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)#is this correct?\n",
+    "\n",
+    "H = toy.to_reduced(n_elec=N_ELEC)  # H = h↑ + 0.5 * V\n",
+    "\n",
+    "E_spin = np.einsum('pqrs,pqrs', H, rdm2)\n",
+    "\n",
+    "##teste\n",
+    "H_gem   = toy.to_geminal(H)\n",
+    "rdm_gem = to_geminal_rdm2(rdm2)\n",
+    "E_gem   = np.einsum('pq,pq', H_gem, rdm_gem)\n",
+    "\n",
+    "print(\"E_spin-orbital :\", E_spin)\n",
+    "print(\"E_geminal      :\", E_gem)\n",
     "\n",
-    "V = from_geminal(G, n_orb = n_orb)\n",
-    "print(\"Max error:\", np.max(np.abs(h2_spin - V)))\n"
+    "assert np.allclose(E_spin, E_gem)"
    ]
   }
  ],
diff --git a/moha/molkit/hamiltonians.py b/moha/molkit/hamiltonians.py
index dc5e669..83fc5e1 100644
--- a/moha/molkit/hamiltonians.py
+++ b/moha/molkit/hamiltonians.py
@@ -1,6 +1,10 @@
 r"""Molkit Module."""
 
-from .utils.spinops import antisymmetrize_two_body, get_spin_blocks
+from .utils.spinops import (
+    antisymmetrize_two_body,
+    get_spin_blocks,
+    upscale_one_body
+)
 import numpy as np
 
 
@@ -28,62 +32,7 @@ def __init__(self, one_body, two_body):
         self.n_spin = 2 * self.n_spatial
         self.reduced_ham = None
 
-    def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
-        r"""Convert the one/two body terms from spatial to spin-orbital basis.
-
-        Parameters
-        ----------
-        None
-
-        Returns
-        -------
-        one_body_spin : np.ndarray
-            One-body term in spin-orbital basis in physics notation
-        two_body_spin : np.ndarray
-            Two-body term in spin-orbital basis in physics notation
-
-        Notes
-        -----
-        Rules for the conversion:
-        - The one-body term is converted as follows:
-            - :math:`h_{pq}^{\\alpha \\alpha}=h_{pq}^{\\beta \\beta}
-            =h_{pq}^{\\text{spatial}}`
-            - :math:`h_{pq}^{\\alpha \\beta}=h_{pq}^{\\beta \\alpha}=0`
-        - The two-body term is converted as follows:
-            - :math:`V_{pqrs}^{\\alpha \\alpha \\alpha \\alpha}=\\
-                    V_{pqrs}^{\\alpha \\beta \\alpha \\beta}=\\
-                    V_{pqrs}^{\\beta \\alpha \\beta \\alpha}=\\
-                    V_{pqrs}^{\\text{spatial}}`
-
-        """
-        one_body = self.one_body
-        two_body = self.two_body
-        one_body = np.asarray(one_body)
-        two_body = np.asarray(two_body)
-
-        if one_body.ndim != 2 or one_body.shape[0] != one_body.shape[1]:
-            raise ValueError("one_body must be square (n, n)")
-        n = one_body.shape[0]
-        if two_body.shape != (n, n, n, n):
-            raise ValueError(
-                "two_body must have shape (n, n, n, n) with same n")
-
-        one_body_spin = np.zeros((2 * n, 2 * n), dtype=one_body.dtype)
-        one_body_spin[:n, :n] = one_body   # αα block
-        one_body_spin[n:, n:] = one_body   # ββ block
-
-        two_body_spin = np.zeros((2 * n, 2 * n, 2 * n, 2 * n),
-                                 dtype=two_body.dtype)
-        # αααα
-        two_body_spin[:n, :n, :n, :n] = two_body
-        # ββββ
-        two_body_spin[n:, n:, n:, n:] = two_body
-        # αβαβ
-        two_body_spin[:n, n:, :n, n:] = two_body
-        # βαβα
-        two_body_spin[n:, :n, n:, :n] = two_body
-
-        return one_body_spin, two_body_spin
+    
 
     def antisymmetrize(self):
         """Apply proper antisymmetrization to two-electron integrals."""
@@ -164,6 +113,99 @@ def to_geminal(two_body):
         n_gem = n * (n - 1) // 2
         return np.array(two_body_gem).reshape(n_gem, n_gem)
 
-    def build_reduced(self):
-        """Build the reduced hamiltonian form."""
-        pass
+    def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
+        r"""Convert the one/two body terms from spatial to spin-orbital basis.
+
+        Parameters
+        ----------
+        one_body : np.ndarray
+            One-body term in spatial orbital basis (shape (n, n)).
+        two_body : np.ndarray
+            Two-body term in spatial orbital basis (shape (n, n, n, n)).
+
+        Returns
+        -------
+        one_body_spin : np.ndarray
+            One-body term in spin-orbital basis in physics notation
+        two_body_spin : np.ndarray
+            Two-body term in spin-orbital basis in physics notation
+
+        Notes
+        -----
+        Rules for the conversion:
+        - The one-body term is converted as follows:
+            - :math:`h_{pq}^{\\alpha \\alpha}=h_{pq}^{\\beta \\beta}
+            =h_{pq}^{\\text{spatial}}`
+            - :math:`h_{pq}^{\\alpha \\beta}=h_{pq}^{\\beta \\alpha}=0`
+        - The two-body term is converted as follows:
+            - :math:`V_{pqrs}^{\\alpha \\alpha \\alpha \\alpha}=\\
+                    V_{pqrs}^{\\alpha \\beta \\alpha \\beta}=\\
+                    V_{pqrs}^{\\beta \\alpha \\beta \\alpha}=\\
+                    V_{pqrs}^{\\text{spatial}}`
+
+        """
+        one_body, two_body = self.one_body, self.two_body
+        one_body = np.asarray(one_body)
+        two_body = np.asarray(two_body)
+
+        if one_body.ndim != 2 or one_body.shape[0] != one_body.shape[1]:
+            raise ValueError("one_body must be square (n, n)")
+        n = one_body.shape[0]
+        if two_body.shape != (n, n, n, n):
+            raise ValueError(
+                "two_body must have shape (n, n, n, n) with same n")
+
+        one_body_spin = np.zeros((2 * n, 2 * n), dtype=one_body.dtype)
+        one_body_spin[:n, :n] = one_body   # αα block
+        one_body_spin[n:, n:] = one_body   # ββ block
+
+        two_body_spin = np.zeros((2 * n, 2 * n, 2 * n, 2 * n),
+                                    dtype=two_body.dtype)
+        # αααα
+        two_body_spin[:n, :n, :n, :n] = two_body
+        # ββββ
+        two_body_spin[n:, n:, n:, n:] = two_body
+        # αβαβ
+        two_body_spin[:n, n:, :n, n:] = two_body
+        # βαβα
+        two_body_spin[n:, :n, n:, :n] = two_body
+
+        return one_body_spin, two_body_spin
+
+    def to_reduced(self, n_elec):
+        r"""Return the reduced 4-index Hamiltonian tensor.
+
+        .. math::
+
+        k_{pqrs}=\frac{1}{2\,(N-1)}\bigl(h_{pq}\,\delta_{rs}
+        +h_{rs}\,\delta_{pq}\bigr)+\tfrac12\,V_{pqrs}.
+
+        Parameters
+        ----------
+        one_body_spatial : ndarray, shape (n, n)
+            One‑electron integral matrix :math:`h_{pq}` in a spatial‑orbital
+            basis.
+        two_body_spatial : ndarray, shape (n, n, n, n)
+            Two‑electron integral tensor :math:`V_{pqrs}` (Dirac convention,
+            spatial orbitals).
+        n_elec : int
+            Number of electrons *N* in the system.
+
+        Returns
+        -------
+        ndarray, shape (2n, 2n, 2n, 2n)
+            Reduced Hamiltonian tensor :math:`k_{pqrs}` in the
+            spin‑orbital basis.
+
+        Notes
+        -----
+        The function is stateless; it does not modify the parent
+        ``MolHam`` instance.
+        """
+        h_spin, V_spin = self.spinize_H()
+    
+        h_upscaled = upscale_one_body(h_spin, n_elec)
+
+        k = h_upscaled + 0.5 * V_spin
+        return k
+
diff --git a/moha/molkit/utils/spinops.py b/moha/molkit/utils/spinops.py
index debb1fd..85addfc 100644
--- a/moha/molkit/utils/spinops.py
+++ b/moha/molkit/utils/spinops.py
@@ -91,3 +91,42 @@ def get_spin_blocks(two_body_spin, n_spatial):
         "bbbb": two_body_spin[n_spatial:, n_spatial:, n_spatial:, n_spatial:],
         "abab": two_body_spin[:n_spatial, n_spatial:, :n_spatial, n_spatial:],
     }
+
+def upscale_one_body(one_body, n_elec):
+    """
+    Upscale the 1 body term to the 2 body term.
+    Specifically, the one-body term is upscaled to a 4-dimensional tensor.
+
+    Parameters
+    ----------
+    one_body : np.ndarray
+        One-body term in spin-orbital basis in physics notation.
+    n_elec : int
+        Number of electrons in the system
+
+    Returns
+    -------
+    one_body_up : np.ndarray
+        Upscaled one-body integrals
+
+    Notes
+    -----
+    The upscaling is done by the following formula:
+
+    .. math::
+
+        \\frac{1}{2 (n-1)}(\\mathbf{h}_{pq}\\delta_{rs} + \\mathbf{h}_{rs}\\delta_{pq})
+
+    where:
+    - :math:`\\mathbf{h}_{pq}` and :math:`\\mathbf{h}_{rs}` are elements of the one-body term
+    - :math:`\\delta_{rs}` and :math:`\\delta_{pq}` are Kronecker deltas
+    - :math:`n` is the number of electrons in the system
+
+    The resulting upscaled one-body term is a 4-dimensional tensor.
+
+    """
+    n = one_body.shape[0]
+    eye = np.eye(n)
+    one_body_up = 0.5 * (np.kron(one_body, eye) + np.kron(eye, one_body)) / (n_elec - 1)
+    
+    return one_body_up.reshape(n, n, n, n)
diff --git a/moha/molkit/utils/tools.py b/moha/molkit/utils/tools.py
index 37bf1b0..bc2d1ac 100644
--- a/moha/molkit/utils/tools.py
+++ b/moha/molkit/utils/tools.py
@@ -1,5 +1,6 @@
 r"""Utilities for molecular Hamiltonians."""
 
+import numpy as np
 
 def from_geminal(two_body_gem, n_orb):
     """

From 924f56efa64e9aaea04770e02fa6085579feb33e Mon Sep 17 00:00:00 2001
From: Giovanni Benedetti Da Rosa <benedetti-23@gpu3.enst.fr>
Date: Mon, 28 Jul 2025 00:39:05 +0200
Subject: [PATCH 2/3] #171 Reduded Hamiltonian and Energy Testing

---
 examples/molkit.ipynb        | 224 +++++++----------------------------
 moha/molkit/hamiltonians.py  |  41 ++++---
 moha/molkit/utils/spinops.py |  16 ++-
 moha/molkit/utils/tools.py   |   1 +
 scripts/toml_to_moha.py      |  48 ++++----
 5 files changed, 101 insertions(+), 229 deletions(-)

diff --git a/examples/molkit.ipynb b/examples/molkit.ipynb
index 46331e7..4a06330 100644
--- a/examples/molkit.ipynb
+++ b/examples/molkit.ipynb
@@ -32,7 +32,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -109,181 +109,22 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "metadata": {},
-   "outputs": [],
-   "source": [
-    "import pyscf\n",
-    "from pyscf import gto, scf, fci\n",
-    "from pyscf.tools import fcidump\n",
-    "import os\n",
-    "\n",
-    "\n",
-    "path_xyz = r'../data/hydrogenic/H4/Paldus/h4_0.456_2.xyz'\n",
-    "mol = gto.M(atom=path_xyz, basis='sto-6g')\n",
-    "hf = scf.RHF(mol)\n",
-    "\n",
-    "hf.conv_tol = 1.0e-12\n",
-    "status = hf.run().converged\n",
-    "print(status)\n",
-    "\n",
-    "\n",
-    "mo_energy = hf.mo_energy\n",
-    "mo_occ = hf.mo_occ\n",
-    "\n",
-    "print('starting writing')\n",
-    "# Write FCIDUMP\n",
-    "\n",
-    "fcidump.from_mo(mol, 'tmp1.fcidump',\n",
-    "                hf.mo_coeff)\n",
-    "\n",
-    "# fcidump.from_scf(hf, \n",
-    "# f\"tmp.fcidump\",\n",
-    "# tol=0.0)\n",
-    "\n",
-    "print('done')\n",
-    "\n",
-    "data = fcidump.read(\"tmp1.fcidump\")\n",
-    "e, ci_vec = fci.direct_spin1.kernel(data['H1'],\n",
-    "                                    data['H2'], \n",
-    "                                    data['NORB'], \n",
-    "                                    data['NELEC'], ecore=data['ECORE'], nroots=1)\n",
-    "\n",
-    "rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, data['NORB'], data['NELEC'])\n",
-    "rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
-    "                 [np.zeros((4, 4)), rdm1s[1]]])\n",
-    "\n",
-    "rdm2 = np.zeros((8, 8, 8, 8))\n",
-    "for i in range(4):\n",
-    "    for j in range(4):\n",
-    "        for k in range(4):\n",
-    "            for l in range(4):\n",
-    "                rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
-    "                rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
-    "                rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
-    "                # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
-    "                rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]\n",
-    "                \n",
-    "\n",
-    "# converting rdm2 to physics notation\n",
-    "# rdm2 = np.einsum('ijkl->ikjl', rdm2)\n",
-    "print(\"trace of rdm1\", np.einsum('ii', rdm1))\n",
-    "print(\"Trace of rdm2\", np.einsum('ijij', rdm2))\n",
-    "\n",
-    "# loading integrals\n",
-    "# iodata stores integrals in physics notations\n",
-    "# https://iodata.readthedocs.io/en/latest/_modules/iodata/formats/fcidump.html#load_one\n",
-    "integrals = load_one(open('tmp1.fcidump', 'r'))\n",
-    "os.remove('tmp1.fcidump')\n",
-    "\n",
-    "h1s = integrals['one_ints']['core_mo']\n",
-    "h2s = integrals['two_ints']['two_mo']\n",
-    "\n",
-    "h1 = np.block([[h1s, np.zeros((4, 4))],\n",
-    "               [np.zeros((4, 4)), h1s]])\n",
-    "\n",
-    "h2 = np.zeros((8, 8, 8, 8))\n",
-    "for i in range(4):\n",
-    "    for j in range(4):\n",
-    "        for k in range(4):\n",
-    "            for l in range(4):\n",
-    "                # aaaa\n",
-    "                h2[i, j, k, l] = h2s[i, j, k, l]\n",
-    "                # bbbb\n",
-    "                h2[i+4, j+4, k+4, l+4] = h2s[i, j, k, l]\n",
-    "                # abab\n",
-    "                h2[i, j+4, k, l+4] = h2s[i, j, k, l]\n",
-    "                # baba\n",
-    "# h2[i+4, j, k+4, l] = h2s[k, l, i, j]\n",
-    "                h2[i+4, j, k+4, l] = h2s[i, j, k, l]\n",
-    "                \n",
-    "                \n",
-    "res = 0\n",
-    "for i in range(8):\n",
-    "    for j in range(8):\n",
-    "        assert rdm2[i, j, i, j] >=0\n",
-    "        res += rdm2[i, j, i, j]\n",
-    "        \n",
-    "print(\"trace of 2dm with for loop: \", res)\n",
-    "\n",
-    "e_dm = np.einsum('pq, pq', h1, rdm1) + 0.5 * np.einsum('pqrs, pqrs', h2, rdm2)\n",
-    "\n",
-    "\n",
-    "print(e_dm, e - integrals['core_energy'])\n",
-    "\n",
-    "\n",
-    "#-----------\n",
-    "eye = np.eye(8)\n",
-    "N = 4\n",
-    "N_alpha = N//2\n",
-    "\n",
-    "h1_ups = 0.5 * (np.kron(h1, eye) + np.kron(eye, h1)) / (N-1)\n",
-    "h1_ups = h1_ups.reshape(8, 8, 8, 8)\n",
-    "H = h1_ups + 0.5 * h2\n",
-    "\n",
-    "np.einsum('pqrs, pqrs', H, rdm2), e_dm\n",
-    "\n",
-    "# to geminal:\n",
-    "# H_gem = to_geminal(H)\n",
-    "# rdm_gem = to_geminal(rdm2)\n",
-    "# np.einsum('pq, pq', H_gem, rdm_gem) == np.einsum('pqrs, pqrs', H, rdm2) == e_dm"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "from moha.molkit.hamiltonians import MolHam     \n",
-    "h_sp = np.array([[1.0, 0.2],\n",
-    "                 [0.2, 1.5]])\n",
-    "\n",
-    "V_sp = np.zeros((2, 2, 2, 2))\n",
-    "V_sp[0, 0, 0, 0] = 0.7\n",
-    "V_sp[1, 1, 1, 1] = 0.6\n",
-    "\n",
-    "N_ELEC = 2\n",
-    "\n",
-    "toy_ham = MolHam(one_body=h_sp, two_body=V_sp)\n",
-    "\n",
-    "h_spin, V_spin = toy_ham.spinize_H()     \n",
-    "n_spin = toy_ham.n_spin      \n",
-    "\n",
-    "e_fci, ci_vec = fci.direct_spin1.kernel(\n",
-    "    h_spin, V_spin, n_spin, N_ELEC, ecore=0.0, nroots=1\n",
-    ")\n",
-    "\n",
-    "rdm1s, rdm2s = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)\n",
-    "rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
-    "                 [np.zeros((4, 4)), rdm1s[1]]])\n",
-    "\n",
-    "rdm2 = np.zeros((8, 8, 8, 8))\n",
-    "for i in range(4):\n",
-    "    for j in range(4):\n",
-    "        for k in range(4):\n",
-    "            for l in range(4):\n",
-    "                rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
-    "                rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
-    "                rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
-    "                # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
-    "                rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "E_spin-orbital : 1.004059262354666e-12\n",
+      "E_geminal      : 1.0040474662350293e-12\n"
+     ]
+    }
+   ],
    "source": [
     "import numpy as np\n",
     "from pyscf import fci\n",
-    "\n",
     "from moha.molkit.hamiltonians import MolHam\n",
-    "from moha.utils.spinops import (\n",
-    "    to_geminal_rdm2\n",
-    ")\n",
     "\n",
     "\n",
     "h_sp = np.array([[1.0, 0.2],\n",
@@ -293,33 +134,50 @@
     "V_sp[0, 0, 0, 0] = 0.7\n",
     "V_sp[1, 1, 1, 1] = 0.6\n",
     "\n",
-    "N_ELEC = 2                   \n",
+    "N_ELEC, n_spin  = 2, 2                  \n",
     "toy = MolHam(h_sp, V_sp)\n",
     "\n",
     "h_spin, V_spin = toy.spinize_H()           \n",
-    "\n",
-    "e_fci, ci_vec = fci.direct_spin1.kernel(h_spin, V_spin, n_spin, N_ELEC, ecore=0.0)\n",
-    "rdm1, rdm2 = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)#is this correct?\n",
+    "norb = h_spin.shape[0]\n",
+    "e_fci, ci_vec = fci.direct_spin1.kernel(h_spin, V_spin, norb, N_ELEC, ecore=0.0)\n",
+    "rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, norb, N_ELEC)\n",
+    "\n",
+    "rdm1 = np.block([\n",
+    "    [rdm1s[0], np.zeros_like(rdm1s[0])],\n",
+    "    [np.zeros_like(rdm1s[1]), rdm1s[1]]\n",
+    "])\n",
+    "\n",
+    "n = norb // 2  # spatial orbitals\n",
+    "rdm2 = np.zeros((2*n, 2*n, 2*n, 2*n))\n",
+    "for i in range(n):\n",
+    "    for j in range(n):\n",
+    "        for k in range(n):\n",
+    "            for l in range(n):\n",
+    "                # αααα\n",
+    "                rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
+    "                # ββββ\n",
+    "                rdm2[i+n, j+n, k+n, l+n] = rdm2s[2][i, k, j, l]\n",
+    "                # αβαβ\n",
+    "                rdm2[i, j+n, k, l+n] = rdm2s[1][i, k, j, l]\n",
+    "                # βαβα (hermitian transpose of αβαβ)\n",
+    "                rdm2[i+n, j, k+n, l] = rdm2s[1][i, k, j, l]\n",
     "\n",
     "H = toy.to_reduced(n_elec=N_ELEC)  # H = h↑ + 0.5 * V\n",
-    "\n",
     "E_spin = np.einsum('pqrs,pqrs', H, rdm2)\n",
     "\n",
-    "##teste\n",
-    "H_gem   = toy.to_geminal(H)\n",
-    "rdm_gem = to_geminal_rdm2(rdm2)\n",
+    "H_gem = toy.to_geminal(H, type='h2')\n",
+    "rdm_gem = toy.to_geminal(rdm2, type='rdm2')\n",
     "E_gem   = np.einsum('pq,pq', H_gem, rdm_gem)\n",
     "\n",
     "print(\"E_spin-orbital :\", E_spin)\n",
     "print(\"E_geminal      :\", E_gem)\n",
-    "\n",
-    "assert np.allclose(E_spin, E_gem)"
+    "assert np.allclose(E_spin, E_gem, atol=1e-10), \"Energies do not match!\""
    ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "base",
+   "display_name": "moha",
    "language": "python",
    "name": "python3"
   },
@@ -333,7 +191,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.11.13"
   }
  },
  "nbformat": 4,
diff --git a/moha/molkit/hamiltonians.py b/moha/molkit/hamiltonians.py
index 83fc5e1..11263d0 100644
--- a/moha/molkit/hamiltonians.py
+++ b/moha/molkit/hamiltonians.py
@@ -32,8 +32,6 @@ def __init__(self, one_body, two_body):
         self.n_spin = 2 * self.n_spatial
         self.reduced_ham = None
 
-    
-
     def antisymmetrize(self):
         """Apply proper antisymmetrization to two-electron integrals."""
         self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
@@ -53,13 +51,18 @@ def get_spin_blocks(self):
 
         return get_spin_blocks(self.two_body_spin, self.n_spatial)
 
-    def to_geminal(two_body):
-        r"""Convert the two-body term to the geminal basis.
+    def to_geminal(self, two_body=None, type='h2'):
+        r"""
+        Convert the two-body term to the geminal basis.
 
         Parameters
         ----------
         two_body : np.ndarray
             Two-body term in spin-orbital basis in physics notation.
+        type : str
+            ['rdm2', 'h2']. Type of the two-body term.
+            - 'rdm2' : 2 body reduced density matrix
+            - 'h2' : 2 body Hamiltonian
 
         Returns
         -------
@@ -88,11 +91,6 @@ def to_geminal(two_body):
         V_{A B}
         =\frac{1}{2}\left(V_{p q r s}-V_{q p r s}-V_{p q r s}+V_{qprs}\right)
 
-        This is coming from the fact, that the Hamiltonian object
-        retured from the fcidump (converted to the physics notation)
-        doesn't obbey the permutation rules.
-        Thus, the full formula has to be used.
-
         """
         n = two_body.shape[0]
         two_body_gem = []
@@ -102,13 +100,19 @@ def to_geminal(two_body):
             for q in range(p + 1, n):
                 for r in range(n):
                     for s in range(r + 1, n):
-                        term = 0.5 * (
-                            two_body[p, q, r, s]
-                            - two_body[q, p, r, s]
-                            - two_body[p, q, s, r]
-                            + two_body[q, p, s, r]
-                        )
-                        two_body_gem.append(term)
+                        if type == 'rdm2':
+                            two_body_gem.append(
+                                0.5 * 4 * two_body[p, q, r, s]
+                            )
+                        elif type == 'h2':
+                            two_body_gem.append(
+                                0.5 * (
+                                    two_body[p, q, r, s]
+                                    - two_body[q, p, r, s]
+                                    - two_body[p, q, s, r]
+                                    + two_body[q, p, s, r]
+                                )
+                            )
 
         n_gem = n * (n - 1) // 2
         return np.array(two_body_gem).reshape(n_gem, n_gem)
@@ -160,7 +164,7 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
         one_body_spin[n:, n:] = one_body   # ββ block
 
         two_body_spin = np.zeros((2 * n, 2 * n, 2 * n, 2 * n),
-                                    dtype=two_body.dtype)
+                                 dtype=two_body.dtype)
         # αααα
         two_body_spin[:n, :n, :n, :n] = two_body
         # ββββ
@@ -203,9 +207,8 @@ def to_reduced(self, n_elec):
         ``MolHam`` instance.
         """
         h_spin, V_spin = self.spinize_H()
-    
+
         h_upscaled = upscale_one_body(h_spin, n_elec)
 
         k = h_upscaled + 0.5 * V_spin
         return k
-
diff --git a/moha/molkit/utils/spinops.py b/moha/molkit/utils/spinops.py
index 85addfc..0f38f2f 100644
--- a/moha/molkit/utils/spinops.py
+++ b/moha/molkit/utils/spinops.py
@@ -52,7 +52,7 @@ def antisymmetrize_two_body(
     if nspin % 2:
         raise ValueError("spin-orbital tensor size must be even (2 n)")
 
-    n = nspin // 2                       # number of spatial orbitals
+    n = nspin // 2      # number of spatial orbitals
 
     # αααα block
     aa = tensor[:n, :n, :n, :n]
@@ -92,9 +92,11 @@ def get_spin_blocks(two_body_spin, n_spatial):
         "abab": two_body_spin[:n_spatial, n_spatial:, :n_spatial, n_spatial:],
     }
 
+
 def upscale_one_body(one_body, n_elec):
-    """
+    r"""
     Upscale the 1 body term to the 2 body term.
+
     Specifically, the one-body term is upscaled to a 4-dimensional tensor.
 
     Parameters
@@ -115,10 +117,11 @@ def upscale_one_body(one_body, n_elec):
 
     .. math::
 
-        \\frac{1}{2 (n-1)}(\\mathbf{h}_{pq}\\delta_{rs} + \\mathbf{h}_{rs}\\delta_{pq})
+        \\frac{1}{2 (n-1)}(\\mathbf{h}_{pq}\\delta_{rs} +
+        \\mathbf{h}_{rs}\\delta_{pq})
 
     where:
-    - :math:`\\mathbf{h}_{pq}` and :math:`\\mathbf{h}_{rs}` are elements of the one-body term
+    - :math:`\\mathbf{h}_{pq}` and :math:`\\mathbf{h}_{rs}` are the one-body
     - :math:`\\delta_{rs}` and :math:`\\delta_{pq}` are Kronecker deltas
     - :math:`n` is the number of electrons in the system
 
@@ -127,6 +130,7 @@ def upscale_one_body(one_body, n_elec):
     """
     n = one_body.shape[0]
     eye = np.eye(n)
-    one_body_up = 0.5 * (np.kron(one_body, eye) + np.kron(eye, one_body)) / (n_elec - 1)
-    
+    one_body_up = 0.5 * (np.kron(one_body, eye) +
+                         np.kron(eye, one_body)) / (n_elec - 1)
+
     return one_body_up.reshape(n, n, n, n)
diff --git a/moha/molkit/utils/tools.py b/moha/molkit/utils/tools.py
index bc2d1ac..c3904d7 100644
--- a/moha/molkit/utils/tools.py
+++ b/moha/molkit/utils/tools.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 
+
 def from_geminal(two_body_gem, n_orb):
     """
     Inverse of MolHam.to_geminal().
diff --git a/scripts/toml_to_moha.py b/scripts/toml_to_moha.py
index e5029ad..c9b51eb 100644
--- a/scripts/toml_to_moha.py
+++ b/scripts/toml_to_moha.py
@@ -1,3 +1,4 @@
+import moha
 import tomllib
 import numpy as np
 from os.path import exists
@@ -5,7 +6,7 @@
 from os import makedirs
 import sys
 sys.path.insert(0, '../')
-import moha
+
 
 def set_defaults(input_data):
     '''
@@ -26,11 +27,11 @@ def set_defaults(input_data):
     required_default_paramfile = "defaults.toml"
     if not exists(required_default_paramfile):
         raise Exception("Default input file 'defaults.toml' is required.")
-    
+
     # load defaults.toml data into default_data
     default_data = tomllib.load(
-                       open(required_default_paramfile, "rb")
-            )
+        open(required_default_paramfile, "rb")
+    )
 
     # set defaults in input_data
     for param_type in default_data.keys():
@@ -53,7 +54,6 @@ def set_defaults(input_data):
                 input_data[param_type][param] = data_value.lower()
 
 
-
 def build_moha_moltype_1d(data):
     '''
     Function that builds and returns hamiltonian object 
@@ -73,14 +73,14 @@ def build_moha_moltype_1d(data):
         moha.Ham
     '''
     # define parameters for 1d model
-    norb      = data["system"]["norb"]
-    charge    = data["model"]["charge"]
-    alpha     = data["model"]["alpha"]
-    beta      = data["model"]["beta"]
-    u_onsite  = data["model"]["u_onsite"]
-    mu        = data["model"]["mu"]
-    J_eq      = data["model"]["J_eq"]
-    J_ax      = data["model"]["J_ax"]
+    norb = data["system"]["norb"]
+    charge = data["model"]["charge"]
+    alpha = data["model"]["alpha"]
+    beta = data["model"]["beta"]
+    u_onsite = data["model"]["u_onsite"]
+    mu = data["model"]["mu"]
+    J_eq = data["model"]["J_eq"]
+    J_ax = data["model"]["J_ax"]
 
     # build connectivity
     connectivity = [(f"C{i}", f"C{i + 1}", 1) for i in range(1, norb)]
@@ -90,10 +90,10 @@ def build_moha_moltype_1d(data):
     # build connectivity for spin models
     spin_connectivity = np.array(np.eye(norb, k=1) + np.eye(norb, k=-1))
     if data["system"]["bc"] == "periodic":
-        spin_connectivity[(0,-1)] = spin_connectivity[(-1,0)] = 1
+        spin_connectivity[(0, -1)] = spin_connectivity[(-1, 0)] = 1
 
     # create and return hamiltonian object ham
-    #-- Fermion models --#
+    # -- Fermion models --#
     # PPP
     if data["model"]["hamiltonian"] == "ppp":
         charge_arr = charge * np.ones(norb)
@@ -108,12 +108,14 @@ def build_moha_moltype_1d(data):
     # Hubbard
     elif data["model"]["hamiltonian"] == "hubbard":
         u_onsite_arr = u_onsite * np.ones(norb)
-        ham = moha.HamHub(connectivity=connectivity, alpha=alpha, beta=beta, u_onsite=u_onsite_arr)
+        ham = moha.HamHub(connectivity=connectivity,
+                          alpha=alpha, beta=beta, u_onsite=u_onsite_arr)
         return ham
-    #-- Spin models --#
+    # -- Spin models --#
     # Heisenberg
     elif data["model"]["hamiltonian"] == "heisenberg":
-        ham = moha.HamHeisenberg(connectivity=spin_connectivity, mu=mu, J_eq=J_eq, J_ax=J_ax)
+        ham = moha.HamHeisenberg(
+            connectivity=spin_connectivity, mu=mu, J_eq=J_eq, J_ax=J_ax)
         return ham
     # Ising
     elif data["model"]["hamiltonian"] == "ising":
@@ -124,9 +126,10 @@ def build_moha_moltype_1d(data):
         ham = moha.HamRG(connectivity=spin_connectivity, mu=mu, J_eq=J_eq)
         return ham
     else:
-        raise ValueError("Model hamiltonian " + data["model"]["hamiltonian"] + 
+        raise ValueError("Model hamiltonian " + data["model"]["hamiltonian"] +
                          " not supported for moltype " + data["system"]["moltype"] + ".")
 
+
 def dict_to_ham(data):
     '''
     Function for generating hamiltonian from dictionary of model data.
@@ -150,7 +153,8 @@ def dict_to_ham(data):
     if data["system"]["moltype"] == "1d":
         ham = build_moha_moltype_1d(data)
     else:
-        raise ValueError("moltype " + data["system"]["moltype"] + " not supported.")
+        raise ValueError(
+            "moltype " + data["system"]["moltype"] + " not supported.")
 
     # get symmetry of two-electron integrals
     sym = data["system"]["symmetry"]
@@ -175,10 +179,12 @@ def dict_to_ham(data):
         elif data["control"]["integral_format"] == "npz":
             ham.savez(out_file)
         else:
-            raise ValueError("Integral output format " + data["control"]["integral_format"] + " not supported.")
+            raise ValueError("Integral output format " +
+                             data["control"]["integral_format"] + " not supported.")
 
     return ham
 
+
 if __name__ == '__main__':
     toml_file = sys.argv[1]
     data = tomllib.load(open(toml_file, "rb"))

From 451094d0e86f4eabe4220599295629bd0c4ec89b Mon Sep 17 00:00:00 2001
From: Giovanni Benedetti da Rosa <giovannibenedettidarosa@gmail.com>
Date: Mon, 28 Jul 2025 00:45:10 +0200
Subject: [PATCH 3/3] #171 Reduded Hamiltonian and Energy Testing

---
 examples/molkit.ipynb | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/examples/molkit.ipynb b/examples/molkit.ipynb
index 4a06330..3baab91 100644
--- a/examples/molkit.ipynb
+++ b/examples/molkit.ipynb
@@ -109,7 +109,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -164,7 +164,6 @@
     "\n",
     "H = toy.to_reduced(n_elec=N_ELEC)  # H = h↑ + 0.5 * V\n",
     "E_spin = np.einsum('pqrs,pqrs', H, rdm2)\n",
-    "\n",
     "H_gem = toy.to_geminal(H, type='h2')\n",
     "rdm_gem = toy.to_geminal(rdm2, type='rdm2')\n",
     "E_gem   = np.einsum('pq,pq', H_gem, rdm_gem)\n",