Molham test

docs/molham.rst

@@ -0,0 +1,8 @@
+Molecular Hamiltonian
+=====================
+MolHam class
+-----------------
+... autoclass:: moha.molkit.hamiltonians.MolHam
+    :members:
+
+    .. automethod:: __init__

moha/molkit/hamiltonians.py

@@ -42,7 +42,11 @@ def __init__(self, one_body=None, two_body=None, zero_body=0):
 
     def antisymmetrize(self):
         """Apply proper antisymmetrization to two-electron integrals."""
-        self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
+        if not hasattr(self, "two_body_spin"):
+            raise RuntimeError(
+                "Call .spinize_H() first to compute spin-orbital form.")
+        self.two_body_spin = antisymmetrize_two_body(self.two_body_spin,
+                                                     inplace=True)
 
     def get_spin_blocks(self):
         """Return the main spin blocks of the two-body spin-orbital tensor.

moha/molkit/utils/spinops.py

@@ -20,7 +20,7 @@ def antisymmetrize_two_body(
     Returns
     -------
     np.ndarray
-        Antisymmetrised tensor obeying
+        Antisymmetrised tensor obeys
 
         .. math::
 
@@ -32,18 +32,28 @@ def antisymmetrize_two_body(
 
     Notes
     -----
-    The operation applied is
+    https://vergil.chemistry.gatech.edu/static/content/permsymm.pdf
+
+    The input tensor obbeys the following 8-fold symmetry
 
     .. math::
 
-        V^{\\sigma\\sigma\\sigma\\sigma}_{pqrs}
-            \\;\\;\\leftarrow\\;\\;
-            \\tfrac12\\,\\bigl(V^{\\sigma\\sigma\\sigma\\sigma}_{pqrs}
-                              -V^{\\sigma\\sigma\\sigma\\sigma}_{pqsr}\\bigr),
+        V_{ijkl} = V_{jilk}
+                 = V_{klij}
+                 = V_{lkji}
+                 = V_{kjjl}
+                 = V_{lijk}
+                 = V_{ilkj}
+                 = V_{jklj}
+
+    The operation applied is
 
-    for ``σ = α`` and ``σ = β``.  All other spin sectors are returned
-    untouched.
+    .. math::
+        \langle ij\| kl\rangle=
+        \langle ij\mid kl\rangle-\langle ij\mid lk\rangle
 
+    Keep in mind, that such operation produces terms of the form
+    abba, baab, that are not present in the original tensor.
     """
     if not inplace:
         tensor = tensor.copy()
@@ -54,15 +64,7 @@ def antisymmetrize_two_body(
 
     n = nspin // 2      # number of spatial orbitals
 
-    # αααα block
-    aa = tensor[:n, :n, :n, :n]
-    aa -= np.swapaxes(aa, 2, 3)
-    aa *= 0.5
-
-    # ββββ block
-    bb = tensor[n:, n:, n:, n:]
-    bb -= np.swapaxes(bb, 2, 3)
-    bb *= 0.5
+    tensor = tensor - np.einsum('pqrs->pqsr', tensor)
 
     return tensor
 

moha/molkit/utils/tools.py

@@ -55,16 +55,14 @@ def to_geminal(two_body=None, type='h2'):
                 for s in range(r + 1, n):
                     if type == 'rdm2':
                         two_body_gem.append(
-                            0.5 * 4 * two_body[p, q, r, s]
+                            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

moha/test/test_molham.py

@@ -0,0 +1,70 @@
+"""Testing the MolHam class."""
+
+import numpy as np
+from numpy.testing import assert_allclose, assert_equal
+from moha.molkit import utils
+from moha.molkit.hamiltonians import MolHam
+
+
+def test_antisymmetrize():
+    """Test antisymmetrization of two-electron integrals."""
+    one_body = np.eye(2)
+    two_body = np.zeros((2, 2, 2, 2))
+    two_body[0, 0, 1, 1] = 1.0
+    two_body[1, 1, 0, 0] = 1.0
+
+    mol_ham = MolHam(one_body=one_body, two_body=two_body)
+    mol_ham.spinize_H()
+    mol_ham.antisymmetrize()
+
+    expected_two_body_aa = np.zeros((2, 2, 2, 2))
+    expected_two_body_aa[0, 0, 1, 1] = 1
+    expected_two_body_aa[1, 1, 0, 0] = 1
+    expected_two_body_aa[0, 1, 0, 1] = 0
+    expected_two_body_aa[1, 0, 1, 0] = 0
+
+    assert_allclose(mol_ham.two_body[:2, :2, :2, :2], expected_two_body_aa)
+
+
+def test_to_geminal():
+    """Test conversion to geminal basis."""
+    n_orb = 4
+    two_body = np.zeros((2, 2, 2, 2))
+    two_body[0, 0, 1, 1] = 1.0
+    two_body[1, 1, 0, 0] = 1.0
+    two_body[1, 0, 0, 1] = 1.0
+    two_body[0, 1, 1, 0] = 1.0
+
+    geminal_true = np.array([
+        [-2.0]
+    ])
+    geminal = utils.tools.to_geminal(two_body, type='h2')
+    assert_equal(geminal.shape, (1, 1))
+    assert_allclose(geminal, geminal_true)
+
+
+def test_to_reduced_ham():
+    """Test conversion to reduced Hamiltonian."""
+    n_orb = 2
+    one_body = np.eye(n_orb)
+    two_body = np.zeros((n_orb, n_orb, n_orb, n_orb))
+    two_body[0, 0, 1, 1] = 1.0
+    two_body[1, 1, 0, 0] = 1.0
+
+    mol_ham = MolHam(one_body=one_body, two_body=two_body)
+    reduced_ham = mol_ham.to_reduced(n_elec=2)
+
+    # sum over the spin-orbital indices
+    reduced_ham = reduced_ham[:2, :2, :2, :2] +\
+        reduced_ham[2:, 2:, 2:, 2:] +\
+        reduced_ham[:2, 2:, :2, 2:] +\
+        reduced_ham[2:, :2, 2:, :2]
+    reduced_ham *= 0.25
+
+    reduced_ham_true = 0.5 * two_body
+    reduced_ham_true[0, 0, 0, 0] = 1
+    reduced_ham_true[1, 1, 1, 1] = 1
+    reduced_ham_true[0, 1, 0, 1] = 1
+    reduced_ham_true[1, 0, 1, 0] = 1
+
+    assert_allclose(reduced_ham, reduced_ham_true)

website/_toc.yml

@@ -12,6 +12,7 @@ parts:
       - file: examples/mohagpt.ipynb
       - file: examples/gui.ipynb
       - file: examples/RG.ipynb
+      - file: examples/molkit.ipynb
 
     - caption: API
       chapters: