Improved generation of denominators for examples

examples/numerical/examples_helpers.py

@@ -23,6 +23,29 @@ def generate_equation(mbeq, nocc, nvir):
     return funct
 
 
+def update_cc_amplitudes(T, R, invD, rank: int):
+    """
+    A function to updated the CCSD amplitudes
+
+    Parameters
+    ----------
+    T : dict[np.ndarray]
+        The cluster amplitudes
+    R : dict[np.ndarray]
+        The CC residual
+    invD : dict[np.ndarray]
+        The inverse MP denominators
+    rank : int
+        The rank of the CC equations (e.g., CCSD : rank = 2)
+    """
+    if rank >= 1:
+        T["ov"] += np.einsum("ia,ia->ia", R["ov"], invD["ov"])
+    if rank >= 2:
+        T["oovv"] += np.einsum("ijab,ijab->ijab", R["oovv"], invD["oovv"])
+    if rank >= 3:
+        T["ooovvv"] += np.einsum("ijkabc,ijkabc->ijkabc", R["ooovvv"], invD["ooovvv"])
+
+
 def antisymmetrize_residual_2_2(Roovv, nocc, nvir):
     # antisymmetrize the residual
     Roovv_anti = np.zeros((nocc, nocc, nvir, nvir))
@@ -33,11 +56,6 @@ def antisymmetrize_residual_2_2(Roovv, nocc, nvir):
     return Roovv_anti
 
 
-def update_amplitudes_ccsd(T, R, invD):
-    T["ov"] += np.einsum("ia,ia->ia", R["ov"], invD["ov"])
-    T["oovv"] += np.einsum("ijab,ijab->ijab", R["oovv"], invD["oovv"])
-
-
 def antisymmetrize_residual_3_3(Rooovvv, nocc, nvir):
     # antisymmetrize the residual
     Rooovvv_anti = np.zeros((nocc, nocc, nocc, nvir, nvir, nvir))
@@ -80,68 +98,34 @@ def antisymmetrize_residual_3_3(Rooovvv, nocc, nvir):
     return Rooovvv_anti
 
 
-def update_amplitudes_ccsdt(T, R, invD):
-    T["ov"] += np.einsum("ia,ia->ia", R["ov"], invD["ov"])
-    T["oovv"] += np.einsum("ijab,ijab->ijab", R["oovv"], invD["oovv"])
-    T["ooovvv"] += np.einsum("ijkabc,ijkabc->ijkabc", R["ooovvv"], invD["ooovvv"])
-
-
-def compute_inverse_denominators(H, nocc, nvir, rank):
+def compute_inverse_denominators(H: dict, nocc: list[int], nvir: list[int], rank: int):
+    """
+    A function to compute the inverse of Møller-Plesset denominators
+    """
     fo = np.diag(H["oo"])
     fv = np.diag(H["vv"])
 
     D = {}
 
-    if rank < 1:
-        return D
-
-    d1 = np.zeros((nocc, nvir))
-    for i in range(nocc):
-        for a in range(nvir):
-            si = i % 2
-            sa = a % 2
-            if si == sa:
-                d1[i][a] = 1.0 / (fo[i] - fv[a])
-    D["ov"] = d1
-
-    if rank < 2:
-        return D
-
-    d2 = np.zeros((nocc, nocc, nvir, nvir))
-    for i in range(nocc):
-        for j in range(nocc):
-            for a in range(nvir):
-                for b in range(nvir):
-                    si = i % 2
-                    sj = j % 2
-                    sa = a % 2
-                    sb = b % 2
-                    if si == sj == sa == sb:
-                        d2[i][j][a][b] = 1.0 / (fo[i] + fo[j] - fv[a] - fv[b])
-                    if si == sa and sj == sb and si != sj:
-                        d2[i][j][a][b] = 1.0 / (fo[i] + fo[j] - fv[a] - fv[b])
-                    if si == sb and sj == sa and si != sj:
-                        d2[i][j][a][b] = 1.0 / (fo[i] + fo[j] - fv[a] - fv[b])
-    D["oovv"] = d2
-
-    if rank < 3:
-        return D
-
-    d3 = np.zeros((nocc, nocc, nocc, nvir, nvir, nvir))
-    for i in range(nocc):
-        for j in range(nocc):
-            for k in range(nocc):
-                for a in range(nvir):
-                    for b in range(nvir):
-                        for c in range(nvir):
-                            si = i % 2
-                            sj = j % 2
-                            sk = k % 2
-                            sa = a % 2
-                            sb = b % 2
-                            sc = c % 2
-                            d3[i][j][k][a][b][c] = 1.0 / (
-                                fo[i] + fo[j] + fo[k] - fv[a] - fv[b] - fv[c]
-                            )
-    D["ooovvv"] = d3
+    if rank >= 1:
+        D["ov"] = 1.0 / (fo.reshape(-1, 1) - fv)
+
+    if rank >= 2:
+        D["oovv"] = 1.0 / (
+            fo.reshape(-1, 1, 1, 1) + fo.reshape(-1, 1, 1) - fv.reshape(-1, 1) - fv
+        )
+
+    if rank >= 3:
+        D["ooovvv"] = 1.0 / (
+            fo.reshape(-1, 1, 1, 1, 1, 1)
+            + fo.reshape(-1, 1, 1, 1, 1)
+            + fo.reshape(-1, 1, 1, 1)
+            - fv.reshape(-1, 1, 1)
+            - fv.reshape(-1, 1)
+            - fv
+        )
+    if rank > 3:
+        raise ValueError(
+            f"compute_inverse_denominators() supports rank up to 3, but was called with rank = {rank}"
+        )
     return D

examples/numerical/spinorbital-CCSD.ipynb

@@ -254,7 +254,7 @@
       "CCSD correlation energy             -0.107582941496 [Eh]\n",
       "Reference CCSD correlation energy   -0.107582941213 [Eh]\n",
       "Error                               -2.825710915255e-10 [Eh]\n",
-      "Timing                              +1.613265000000e-01 [s]\n"
+      "Timing                              +1.642418330000e-01 [s]\n"
      ]
     }
    ],
@@ -282,7 +282,7 @@
     "    R[\"oovv\"] = antisymmetrize_residual_2_2(Roovv, nocc, nvir)\n",
     "\n",
     "    # 2. amplitude update\n",
-    "    update_amplitudes_ccsd(T, R, invD)\n",
+    "    update_cc_amplitudes(T, R, invD, 2)\n",
     "\n",
     "    # 3. check for convergence\n",
     "    norm_R = np.sqrt(np.linalg.norm(R[\"ov\"]) ** 2 + np.linalg.norm(R[\"oovv\"]) ** 2)\n",

examples/numerical/spinorbital-CCSDT.ipynb

@@ -214,7 +214,7 @@
       "CCSDT correlation energy             -0.108354658566 [Eh]\n",
       "Reference CCSDT correlation energy   -0.108354659115 [Eh]\n",
       "Error                               +5.487632293022e-10 [Eh]\n",
-      "Timing                              +1.901947000000e+00 [s]\n"
+      "Timing                              +1.924095000000e+00 [s]\n"
      ]
     }
    ],
@@ -248,7 +248,7 @@
     "    R[\"ooovvv\"] = antisymmetrize_residual_3_3(Rooovvv, nocc, nvir)\n",
     "\n",
     "    # 2. amplitude update\n",
-    "    update_amplitudes_ccsdt(T, R, invD)\n",
+    "    update_cc_amplitudes(T, R, invD, 3)\n",
     "\n",
     "    # 3. check for convergence\n",
     "    norm_R = np.sqrt(np.linalg.norm(R[\"ov\"]) ** 2 + np.linalg.norm(R[\"oovv\"]) ** 2)\n",

examples/numerical/spinorbital-MP2.ipynb

@@ -125,22 +125,6 @@
     "exec(funct)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "id": "64e8285b",
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "# psi4.set_options({'mp2_type': 'conv','freeze_core': True})\n",
-    "# Emp2_psi4 = psi4.energy('mp2')\n",
-    "# print(f\"RHF energy:       {Escf:.12f} Eh\")\n",
-    "# print(f\"MP2 energy:       {Emp2_psi4:.12f} Eh\")\n",
-    "# print(f\"MP2 corr. energy: {Emp2_psi4 - Escf:.12f} Eh\")"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 6,
@@ -224,7 +208,7 @@
     "print(\"-\" * len(header))\n",
     "\n",
     "print(f\"MP2 corr. energy: {Emp2_wicked:+.12f} Eh\")\n",
-    "print(f\"MP2 corr. energy: {Emp2:+.12f} Eh\")\n",
+    "print(f\"MP2 corr. energy: {Emp2:+.12f} Eh (psi4)\")\n",
     "print(f\"Err corr. energy: {Emp2_wicked - Emp2:+.12f} Eh\")\n",
     "\n",
     "assert np.isclose(Emp2_wicked, Emp2)"