Repository: quantumlib/OpenFermion
Pull request #681 (2020-10-27 11:09:45)

Anti-symmetrized integrals, coulomb and exchange integrals
Added functionality to extract anti-symmetrized integrals, coulomb and exchange integrals.

google-cla[bot]: 
Thanks for your pull request. It looks like this may be your first contribution to a Google open source project (if not, look below for help). Before we can look at your pull request, you'll need to sign a Contributor License Agreement (CLA).

:memo: **Please visit <https://cla.developers.google.com/> to sign.**

Once you've signed (or fixed any issues), please reply here with `@googlebot I signed it!` and we'll verify it.

----

#### What to do if you already signed the CLA

##### Individual signers

*   It's possible we don't have your GitHub username or you're using a different email address on your commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).

##### Corporate signers

*   Your company has a Point of Contact who decides which employees are authorized to participate. Ask your POC to be added to the group of authorized contributors. If you don't know who your Point of Contact is, direct the Google project maintainer to [go/cla#troubleshoot](http://go/cla#troubleshoot) ([Public version](https://opensource.google/docs/cla/#troubleshoot)).
*   The email used to register you as an authorized contributor must be the email used for the Git commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).
*   The email used to register you as an authorized contributor must also be [attached to your GitHub account](https://github.com/settings/emails).
		

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F681) for more info**.

<!-- need_sender_cla -->
Emieeel: @googlebot I signed it!
Emieeel: Ah, so, maybe there is some confusion here.
Chemists often mean with the coulomb/exchange matrix: 
$$
J_{\mu,\nu} = \sum_{\lambda, \sigma} (\mu \nu|\lambda \sigma) D_{\lambda,\sigma}\\
K_{\mu,\nu} =  \sum_{\lambda, \sigma} (\mu \lambda|\nu \sigma) D_{\lambda,\sigma}
$$
in the AO basis, where D_{\lambda,\sigma} is the 1-RDM.
I implemented instead
$$
J_{p,q} = (pp|qq)
K_{p,q} = (pq|qp)
$$
in the MO basis.

Do you want to return the first or the second option? In standard quantum chemistry codes as psi4 and pyscf, usually the first case is returned when calling a function as get_jk().
ncrubin: I think as a first pass we'll want spin-orbital antisymmeterized 2-electron integrals.  1) we can formulate the Hamiltonian this way (with the 1/4 coeff out front) 2) we can use it to calculate overlaps of the 2-body operator and computational basis states that differ by 2 or fewer occupations.   this also let's us calculate the fock operator (when it is not supplied by a computation) as `einsum('ijjk', antisymm_spin_orb_tei[:, :occ, :occ, :])`.  So I think having spin-orbital antisymm 2 electron integrals will be nice to have.
Emieeel: The newest commit contains the following changes:
- The anti-symmetrized integrals in the spin-orbital basis (and so documented)
- A separate function to convert spatial orbitals to spin-orbitals, that is also employed in get_molecular_hamiltonian
- Tests to calculate the HF ground state energy with these anti-symmetrized integrals in spin-orbital basis and with coulomb/exchange integrals in the spatial orbital basis
ncrubin: @Emieeel Thanks for adding the tests.  It looks like formatting and linting are failing.  You can use the `./check/pylint` script and `./check/format-incremental` scripts to help pass this check.  `./check/format-incremental --apply` will reformat your code according to our style guide.
google-cla[bot]: All (the pull request submitter and all commit authors) CLAs are signed, **but** one or more commits were authored or co-authored by someone other than the pull request submitter.

We need to confirm that all authors are ok with their commits being contributed to this project.  Please have them confirm that by leaving a comment that contains only `@googlebot I consent.` in this pull request.

*Note to project maintainer:* There may be cases where the author cannot leave a comment, or the comment is not properly detected as consent.  In those cases, you can manually confirm consent of the commit author(s), and set the `cla` label to `yes` (if enabled on your project).

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F681) for more info**.

<!-- need_author_consent -->
googlebot: All (the pull request submitter and all commit authors) CLAs are signed, **but** one or more commits were authored or co-authored by someone other than the pull request submitter.

We need to confirm that all authors are ok with their commits being contributed to this project.  Please have them confirm that by leaving a comment that contains only `@googlebot I consent.` in this pull request.

*Note to project maintainer:* There may be cases where the author cannot leave a comment, or the comment is not properly detected as consent.  In those cases, you can manually confirm consent of the commit author(s), and set the `cla` label to `yes` (if enabled on your project).

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F681) for more info**.

<!-- need_author_consent -->
ncrubin: @googlebot I consent.
google-cla[bot]: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F681) for more info**.

<!-- ok -->
googlebot: CLAs look good, thanks!

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F681) for more info**.

<!-- ok -->
ncrubin: Hi @Emieeel you have to run the ` ./check/format-incremental --apply` for the script to automatically reformat the code.  As for the coverage I checked the lines not covered and it seems that one was a residual if/else in the test and the others were code paths (MissingCalculations) that were not touched.  I added coverage for those lines and resubmitted to your branch.  It looks like everything is passing now.  

Thank you for the commits.  Feel free to open a PR editing the readme with your name and affiliations.  You are now officially an OpenFermion developer!  Welcome!

Patch:
diff --git a/src/openfermion/chem/__init__.py b/src/openfermion/chem/__init__.py
index 2ab8089..367b2b8 100644
--- a/src/openfermion/chem/__init__.py
+++ b/src/openfermion/chem/__init__.py
@@ -15,17 +15,11 @@ from .chemical_series import (
     make_atomic_ring,
 )
 
-from .molecular_data import (
-    angstroms_to_bohr,
-    bohr_to_angstroms,
-    MolecularData,
-    name_molecule,
-    geometry_from_file,
-    load_molecular_hamiltonian,
-    periodic_table,
-    periodic_hash_table,
-    periodic_polarization,
-)
+from .molecular_data import (angstroms_to_bohr, bohr_to_angstroms,
+                             MolecularData, name_molecule, geometry_from_file,
+                             load_molecular_hamiltonian, periodic_table,
+                             periodic_hash_table, periodic_polarization,
+                             antisymtei, j_mat, k_mat)
 
 from .pubchem import geometry_from_pubchem
 
diff --git a/src/openfermion/chem/molecular_data.py b/src/openfermion/chem/molecular_data.py
index d8ef916..e51321e 100644
--- a/src/openfermion/chem/molecular_data.py
+++ b/src/openfermion/chem/molecular_data.py
@@ -174,6 +174,92 @@ def geometry_from_file(file_name):
     return geometry
 
 
+def antisymtei(two_body_integrals):
+    """
+    Args:
+    two_body_integrals : Numpy array of two-electron integrals with OpenFermion
+        Ordering.
+
+    Returns:
+    antisymints : Numpy array of anti-symmetrized integrals
+        <ij||kl> = <ij|kl> - <ij|lk> (physicist ordering).
+    """
+    symints = numpy.copy(two_body_integrals.transpose(0, 1, 3, 2), order='C')
+    antisymints = symints - two_body_integrals
+    return antisymints
+
+
+def j_mat(two_body_integrals):
+    """
+    Args:
+    two_body_integrals : Numpy array of two-electron integrals with OpenFermion
+        Ordering.
+
+    Returns:
+    j_matr : Numpy array of the coulomb integrals J_{p,q} = (pp|qq)
+        (in chemist notation).
+    """
+    chem_ordering = numpy.copy(two_body_integrals.transpose(0, 3, 1, 2),
+                               order='C')
+    return numpy.einsum('iijj -> ij', chem_ordering)
+
+
+def k_mat(two_body_integrals):
+    """
+    Args:
+    two_body_integrals : Numpy array of two-electron integrals with OpenFermion
+        Ordering.
+
+    Returns:
+    k_matr : Numpy array of the exchange integrals K_{p,q} = (pq|qp)
+        (in chemist notation).
+    """
+    chem_ordering = numpy.copy(two_body_integrals.transpose(0, 3, 1, 2),
+                               order='C')
+    return numpy.einsum('ijji -> ij', chem_ordering)
+
+
+def spinorb_from_spatial(one_body_integrals, two_body_integrals):
+    n_qubits = 2 * one_body_integrals.shape[0]
+
+    # Initialize Hamiltonian coefficients.
+    one_body_coefficients = numpy.zeros((n_qubits, n_qubits))
+    two_body_coefficients = numpy.zeros(
+        (n_qubits, n_qubits, n_qubits, n_qubits))
+    # Loop through integrals.
+    for p in range(n_qubits // 2):
+        for q in range(n_qubits // 2):
+
+            # Populate 1-body coefficients. Require p and q have same spin.
+            one_body_coefficients[2 * p, 2 * q] = one_body_integrals[p, q]
+            one_body_coefficients[2 * p + 1, 2 * q +
+                                  1] = one_body_integrals[p, q]
+            # Continue looping to prepare 2-body coefficients.
+            for r in range(n_qubits // 2):
+                for s in range(n_qubits // 2):
+
+                    # Mixed spin
+                    two_body_coefficients[2 * p, 2 * q + 1, 2 * r + 1, 2 *
+                                          s] = (two_body_integrals[p, q, r, s])
+                    two_body_coefficients[2 * p + 1, 2 * q, 2 * r, 2 * s +
+                                          1] = (two_body_integrals[p, q, r, s])
+
+                    # Same spin
+                    two_body_coefficients[2 * p, 2 * q, 2 * r, 2 *
+                                          s] = (two_body_integrals[p, q, r, s])
+                    two_body_coefficients[2 * p + 1, 2 * q + 1, 2 * r +
+                                          1, 2 * s +
+                                          1] = (two_body_integrals[p, q, r, s])
+
+    # Truncate.
+    one_body_coefficients[
+        numpy.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
+    two_body_coefficients[
+        numpy.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.
+
+    return one_body_coefficients, two_body_coefficients
+
+
 class MolecularData(object):
     """Class for storing molecule data from a fixed basis set at a fixed
     geometry that is obtained from classical electronic structure
@@ -741,7 +827,7 @@ class MolecularData(object):
                 shape of (n_orbitals, n_orbitals, n_orbitals, n_orbitals).
 
         Raises:
-          MisissingCalculationError: If integrals are not calculated.
+          MissingCalculationError: If integrals are not calculated.
         """
         # Make sure integrals have been computed.
         if self.one_body_integrals is None or self.two_body_integrals is None:
@@ -838,51 +924,12 @@ class MolecularData(object):
                 get_active_space_integrals(occupied_indices, active_indices)
             constant = self.nuclear_repulsion + core_adjustment
 
-        n_qubits = 2 * one_body_integrals.shape[0]
-
-        # Initialize Hamiltonian coefficients.
-        one_body_coefficients = numpy.zeros((n_qubits, n_qubits))
-        two_body_coefficients = numpy.zeros(
-            (n_qubits, n_qubits, n_qubits, n_qubits))
-        # Loop through integrals.
-        for p in range(n_qubits // 2):
-            for q in range(n_qubits // 2):
-
-                # Populate 1-body coefficients. Require p and q have same spin.
-                one_body_coefficients[2 * p, 2 * q] = one_body_integrals[p, q]
-                one_body_coefficients[2 * p + 1, 2 * q +
-                                      1] = one_body_integrals[p, q]
-                # Continue looping to prepare 2-body coefficients.
-                for r in range(n_qubits // 2):
-                    for s in range(n_qubits // 2):
-
-                        # Mixed spin
-                        two_body_coefficients[2 * p, 2 * q + 1, 2 * r + 1, 2 *
-                                              s] = (
-                                                  two_body_integrals[p, q, r, s]
-                                                  / 2.)
-                        two_body_coefficients[2 * p + 1, 2 * q, 2 * r, 2 * s +
-                                              1] = (
-                                                  two_body_integrals[p, q, r, s]
-                                                  / 2.)
-
-                        # Same spin
-                        two_body_coefficients[2 * p, 2 * q, 2 * r, 2 * s] = (
-                            two_body_integrals[p, q, r, s] / 2.)
-                        two_body_coefficients[2 * p + 1, 2 * q + 1, 2 * r +
-                                              1, 2 * s + 1] = (
-                                                  two_body_integrals[p, q, r, s]
-                                                  / 2.)
-
-        # Truncate.
-        one_body_coefficients[
-            numpy.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
-        two_body_coefficients[
-            numpy.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.
+        one_body_coefficients, two_body_coefficients = spinorb_from_spatial(
+            one_body_integrals, two_body_integrals)
 
         # Cast to InteractionOperator class and return.
         molecular_hamiltonian = reps.InteractionOperator(
-            constant, one_body_coefficients, two_body_coefficients)
+            constant, one_body_coefficients, 1 / 2 * two_body_coefficients)
 
         return molecular_hamiltonian
 
@@ -926,6 +973,59 @@ class MolecularData(object):
         rdm = reps.InteractionRDM(one_rdm, two_rdm)
         return rdm
 
+    def get_j(self):
+        """Method to return coulomb matrix.
+
+        Returns:
+          k_matr : Numpy array of the coulomb integrals J_{p,q} = (pp|qq)
+            (in chemist notation).
+
+        Raises:
+          MissingCalculationError: If integrals are not calculated.
+        """
+        # Make sure integrals have been computed.
+        if self.two_body_integrals is None:
+            raise MissingCalculationError(
+                'Missing integral calculation in {}, run before loading '
+                'integrals.'.format(self.filename))
+        return j_mat(self.two_body_integrals)
+
+    def get_k(self):
+        """Method to return exchange matrix.
+
+        Returns:
+          k_matr : Numpy array of the coulomb integrals K_{p,q} = (pq|qp)
+            (in chemist notation).
+
+        Raises:
+          MissingCalculationError: If integrals are not calculated.
+        """
+        # Make sure integrals have been computed.
+        if self.two_body_integrals is None:
+            raise MissingCalculationError(
+                'Missing integral calculation in {}, run before loading '
+                'integrals.'.format(self.filename))
+        return k_mat(self.two_body_integrals)
+
+    def get_antisym(self):
+        """Method to return anti-symmetrized integrals in spin-orbital basis.
+
+        Returns:
+          antisymints : Numpy array of anti-symmetrized integrals in spin-
+            orbital basis <ij||kl> = <ij|kl> - <ij|lk> (physicist ordering)
+
+        Raises:
+          MissingCalculationError: If integrals are not calculated.
+        """
+        # Make sure integrals have been computed.
+        if self.two_body_integrals is None:
+            raise MissingCalculationError(
+                'Missing integral calculation in {}, run before loading '
+                'integrals.'.format(self.filename))
+        _, two_body_coefficients = spinorb_from_spatial(self.one_body_integrals,
+                                                        self.two_body_integrals)
+        return antisymtei(two_body_coefficients)
+
 
 def load_molecular_hamiltonian(geometry,
                                basis,


Test patch:
diff --git a/src/openfermion/chem/molecular_data_test.py b/src/openfermion/chem/molecular_data_test.py
index d0afbd3..c6776a7 100644
--- a/src/openfermion/chem/molecular_data_test.py
+++ b/src/openfermion/chem/molecular_data_test.py
@@ -14,6 +14,7 @@ import os
 import unittest
 import numpy.random
 import scipy.linalg
+import numpy as np
 
 from openfermion.config import THIS_DIRECTORY, DATA_DIRECTORY
 from openfermion.hamiltonians import jellium_model
@@ -314,3 +315,42 @@ class MolecularDataTest(unittest.TestCase):
                                                      format(bond_length), None,
                                                      None)
         self.assertEqual(count_qubits(lih_hamiltonian), 12)
+
+    def test_jk_matr(self):
+        h2mol = self.molecule
+        j = h2mol.get_j()
+        k = h2mol.get_k()
+        pyscf_j = [[0.67448877, 0.6634681], [0.6634681, 0.69739377]]
+        pyscf_k = [[0.67448877, 0.18128881], [0.18128881, 0.69739377]]
+        for p in range(j.shape[0]):
+            for q in range(j.shape[1]):
+                self.assertAlmostEqual(j[p][q], pyscf_j[p][q])
+                self.assertAlmostEqual(k[p][q], pyscf_k[p][q])
+        ndocc = h2mol.n_electrons // 2
+        E1bdy = 2 * h2mol.one_body_integrals[:ndocc, :ndocc]
+        E2bdy = (2 * j[:ndocc, :ndocc]) - k[:ndocc, :ndocc]
+        E_test = h2mol.nuclear_repulsion + E1bdy + E2bdy
+        self.assertAlmostEqual(E_test, h2mol.hf_energy)
+
+    def test_antisymint(self):
+        h2mol = self.molecule
+        antisymm_spin_orb_tei = h2mol.get_antisym()
+        nocc = h2mol.n_electrons
+        mol_H = h2mol.get_molecular_hamiltonian()
+        E1bdy = np.sum(np.diag(mol_H.one_body_tensor[:nocc, :nocc]))
+        E2bdy = 1/2.*np.einsum('ijij',\
+                                 antisymm_spin_orb_tei[:nocc,:nocc,:nocc,:nocc])
+        E_test = h2mol.nuclear_repulsion + E1bdy + E2bdy
+        self.assertAlmostEqual(E_test, h2mol.hf_energy)
+
+    def test_missing_calcs_for_integrals(self):
+        geometry = [('H', (0., 0., 0.)), ('H', (0., 0., 0.8764))]
+        basis = 'sto-3g'
+        multiplicity = 1
+        molecule = MolecularData(geometry, basis, multiplicity)
+        with self.assertRaises(MissingCalculationError):
+            molecule.get_j()
+        with self.assertRaises(MissingCalculationError):
+            molecule.get_k()
+        with self.assertRaises(MissingCalculationError):
+            molecule.get_antisym()
\ No newline at end of file