diff --git a/examples/molkit.ipynb b/examples/molkit.ipynb
index 6a928d1..3baab91 100644
--- a/examples/molkit.ipynb
+++ b/examples/molkit.ipynb
@@ -32,17 +32,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "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,17 +97,86 @@
     "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": 4,
+   "metadata": {},
+   "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",
+    "from moha.molkit.hamiltonians import MolHam\n",
+    "\n",
     "\n",
-    "V = from_geminal(G, n_orb = n_orb)\n",
-    "print(\"Max error:\", np.max(np.abs(h2_spin - V)))\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, n_spin  = 2, 2                  \n",
+    "toy = MolHam(h_sp, V_sp)\n",
+    "\n",
+    "h_spin, V_spin = toy.spinize_H()           \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",
+    "E_spin = np.einsum('pqrs,pqrs', H, 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",
+    "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"
   },
@@ -129,7 +190,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 dc5e669..11263d0 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,12 +32,100 @@ 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)
+
+    def get_spin_blocks(self):
+        """Return the main spin blocks of the two-body spin-orbital tensor.
+
+        Returns
+        -------
+        dict
+            Dictionary with spin blocks: 'aaaa', 'bbbb', 'abab'
+
+        """
+        if not hasattr(self, "two_body_spin"):
+            raise RuntimeError(
+                "Call .spinize() first to compute spin-orbital form.")
+
+        return get_spin_blocks(self.two_body_spin, self.n_spatial)
+
+    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
+        -------
+        two_body_gem : np.ndarray
+            Two-body term in the geminal basis
+
+        Notes
+        -----
+        Assuming that rdm2 obbey the following permutation rules:
+        - :math:`\Gamma_{p q r s}=-\Gamma_{q p r s}=-\Gamma_{p q s r}
+        =\Gamma_{q p s r}`
+        we can convert the two-body term to the geminal basis
+        by the following formula:
+
+        .. math::
+
+            \Gamma_{p q}=0.5 * 4 \Gamma_{p q r s}
+
+        where:
+        - :math:`\Gamma_{p q}` is the two-body term in the geminal basis
+        - :math:`\Gamma_{p q r s}` is the two-body term in the spin-orbital
+        Hamiltonian in the geminal basis is obtained by the following formula:
+
+        .. math::
+
+        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)
+
+        """
+        n = two_body.shape[0]
+        two_body_gem = []
+
+        # i,j,k,l -> pqrs
+        for p in range(n):
+            for q in range(p + 1, n):
+                for r in range(n):
+                    for s in range(r + 1, n):
+                        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)
+
     def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
         r"""Convert the one/two body terms from spatial to spin-orbital basis.
 
         Parameters
         ----------
-        None
+        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
         -------
@@ -56,8 +148,7 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
                     V_{pqrs}^{\\text{spatial}}`
 
         """
-        one_body = self.one_body
-        two_body = self.two_body
+        one_body, two_body = self.one_body, self.two_body
         one_body = np.asarray(one_body)
         two_body = np.asarray(two_body)
 
@@ -85,85 +176,39 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
 
         return one_body_spin, two_body_spin
 
-    def antisymmetrize(self):
-        """Apply proper antisymmetrization to two-electron integrals."""
-        self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
-
-    def get_spin_blocks(self):
-        """Return the main spin blocks of the two-body spin-orbital tensor.
-
-        Returns
-        -------
-        dict
-            Dictionary with spin blocks: 'aaaa', 'bbbb', 'abab'
-
-        """
-        if not hasattr(self, "two_body_spin"):
-            raise RuntimeError(
-                "Call .spinize() first to compute spin-orbital form.")
+    def to_reduced(self, n_elec):
+        r"""Return the reduced 4-index Hamiltonian tensor.
 
-        return get_spin_blocks(self.two_body_spin, self.n_spatial)
+        .. math::
 
-    def to_geminal(two_body):
-        r"""Convert the two-body term to the geminal basis.
+        k_{pqrs}=\frac{1}{2\,(N-1)}\bigl(h_{pq}\,\delta_{rs}
+        +h_{rs}\,\delta_{pq}\bigr)+\tfrac12\,V_{pqrs}.
 
         Parameters
         ----------
-        two_body : np.ndarray
-            Two-body term in spin-orbital basis in physics notation.
+        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
         -------
-        two_body_gem : np.ndarray
-            Two-body term in the geminal basis
+        ndarray, shape (2n, 2n, 2n, 2n)
+            Reduced Hamiltonian tensor :math:`k_{pqrs}` in the
+            spin‑orbital basis.
 
         Notes
         -----
-        Assuming that rdm2 obbey the following permutation rules:
-        - :math:`\Gamma_{p q r s}=-\Gamma_{q p r s}=-\Gamma_{p q s r}
-        =\Gamma_{q p s r}`
-        we can convert the two-body term to the geminal basis
-        by the following formula:
-
-        .. math::
-
-            \Gamma_{p q}=0.5 * 4 \Gamma_{p q r s}
-
-        where:
-        - :math:`\Gamma_{p q}` is the two-body term in the geminal basis
-        - :math:`\Gamma_{p q r s}` is the two-body term in the spin-orbital
-        Hamiltonian in the geminal basis is obtained by the following formula:
-
-        .. math::
-
-        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.
-
+        The function is stateless; it does not modify the parent
+        ``MolHam`` instance.
         """
-        n = two_body.shape[0]
-        two_body_gem = []
+        h_spin, V_spin = self.spinize_H()
 
-        # i,j,k,l -> pqrs
-        for p in range(n):
-            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)
-
-        n_gem = n * (n - 1) // 2
-        return np.array(two_body_gem).reshape(n_gem, n_gem)
+        h_upscaled = upscale_one_body(h_spin, n_elec)
 
-    def build_reduced(self):
-        """Build the reduced hamiltonian form."""
-        pass
+        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..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]
@@ -91,3 +91,46 @@ 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):
+    r"""
+    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 the one-body
+    - :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..c3904d7 100644
--- a/moha/molkit/utils/tools.py
+++ b/moha/molkit/utils/tools.py
@@ -1,5 +1,7 @@
 r"""Utilities for molecular Hamiltonians."""
 
+import numpy as np
+
 
 def from_geminal(two_body_gem, n_orb):
     """
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"))