Reduced Hamiltonian and Energy Testing

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,

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