Repository: quantumlib/OpenFermion
Pull request #68 (2017-10-19 22:54:04)

Added hubbard model with particle-hole symmetry
I added an alternative version of the routine for creating the Fermi Hubbard hamiltonian. In this version the repulsion term U n_{i,\uparrow} n_{i,\downarrow} is replaced by  U (n_{i,\uparrow}-1/2) (n_{i,\downarrow}-1/2), which is unchanged by a particle-hole transformation. This version has the advantage that at chemical potential and magnetic field equal zero, the ground state is always the one at half-filling and the energy is symmetric with respect to chemical potential = 0. This is more convenient for benchmarks than the previous version. Look at the energies for a 2x2 lattice:
 
![image](https://user-images.githubusercontent.com/10538281/31797891-ab032710-b4fe-11e7-90a4-c5b9f49db7ce.png)

For now it is a separate routine. 

googlebot: Thanks for your pull request.  It looks like this may be your first contribution to a Google open source project. 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, please reply here (e.g. `I signed it!`) and we'll verify.  Thanks.

---

- If you've already signed a CLA, it's possible we don't have your GitHub username or you're using a different email address.  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/).
- If your company signed a CLA, they designated a Point of Contact who decides which employees are authorized to participate. You may need to contact the Point of Contact for your company and ask to be added to the group of authorized contributors. If you don't know who your Point of Contact is, direct the project maintainer to go/cla#troubleshoot.
- In order to pass this check, please resolve this problem and have the pull request author add another comment and the bot will run again.

<!-- need_sender_cla -->
babbush: This looks like a valuable contribution. However, I think we should aim to add this functionality to the existing Hubbard model code. For instance, have an optional input parameter (False by default) which if True, creates the Hubbard model you are interested in here. The goal would be to duplicate as little code as possible and only add the minimum number of extra lines required to create the functionality you want when this flag is true.

Otherwise, the contribution looks pretty solid. You might need a few more tests and you do need to sign the CLA if you haven't done so already.
googlebot: We found a Contributor License Agreement for you (the sender of this pull request), but were unable to find agreements for the commit author(s).  If you authored these, maybe you used a different email address in the git commits than was used to sign the CLA ([login here](https://cla.developers.google.com/) to double check)?  If these were authored by someone else, then they will need to sign a CLA as well, and confirm that they're okay with these being contributed to Google.
In order to pass this check, please resolve this problem and have the pull request author add another comment and the bot will run again.

<!-- need_author_cla -->
jromerofontalvo: I updated it to have a single function with a particle_hole_symmetry Boolean option set by default as False.
googlebot: CLAs look good, thanks!

<!-- ok -->
babbush: Thanks for a nice contribution Jhonathan! And welcome to OpenFermion!

Patch:
diff --git a/src/openfermion/hamiltonians/_hubbard.py b/src/openfermion/hamiltonians/_hubbard.py
index 348fbcf..62e4d8f 100644
--- a/src/openfermion/hamiltonians/_hubbard.py
+++ b/src/openfermion/hamiltonians/_hubbard.py
@@ -60,7 +60,8 @@ def down_index(index):
 
 def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
                   chemical_potential=None, magnetic_field=None,
-                  periodic=True, spinless=False):
+                  periodic=True, spinless=False,
+                  particle_hole_symmetry=False):
     """Return symbolic representation of a Fermi-Hubbard Hamiltonian.
 
     Args:
@@ -76,6 +77,13 @@ def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
         spinless: An optional Boolean. If False, each site has spin up
             orbitals and spin down orbitals. If True, return a spinless
             Fermi-Hubbard model.
+        particle_hole_symmetry: an optional Boolean. If False, the repulsion
+            term corresponds to:
+            coulomb sum_{k=1}^{N-1} a_k^dagger a_k a_{k+1}^dagger a_{k+1}
+            If true, the repulsion term is replaced by:
+            coulomb sum_{k=1}^{N-1} (a_k^dagger a_k - 1/2)
+            (a_{k+1}^dagger a_{k+1} - 1/2)
+            which is unchanged under a particle-hole transformation.
         verbose: An optional Boolean. If True, print all second quantized
             terms.
 
@@ -89,6 +97,11 @@ def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
     else:
         n_spin_orbitals = 2 * n_sites
     hubbard_model = FermionOperator((), 0.0)
+    # select particle-hole symmetry
+    if particle_hole_symmetry:
+        coulomb_shift = FermionOperator((), 0.5)
+    else:
+        coulomb_shift = FermionOperator((), 0.0)
 
     # Loop through sites and add terms.
     for site in range(n_sites):
@@ -118,9 +131,11 @@ def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
 
         # Add local pair interaction terms.
         if not spinless:
-            operators = ((up_index(site), 1), (up_index(site), 0),
-                         (down_index(site), 1), (down_index(site), 0))
-            hubbard_model += FermionOperator(operators, coulomb)
+            operator_1 = number_operator(
+                n_spin_orbitals, up_index(site), 1.0) - coulomb_shift
+            operator_2 = number_operator(
+                n_spin_orbitals, down_index(site), 1.0) - coulomb_shift
+            hubbard_model += coulomb * operator_1 * operator_2
 
         # Index coupled orbitals.
         right_neighbor = site + 1
@@ -137,9 +152,11 @@ def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
         if (site + 1) % x_dimension or (periodic and x_dimension > 2):
             if spinless:
                 # Add Coulomb term.
-                operators = ((site, 1), (site, 0),
-                             (right_neighbor, 1), (right_neighbor, 0))
-                hubbard_model += FermionOperator(operators, coulomb)
+                operator_1 = number_operator(
+                    n_spin_orbitals, site, 1.0) - coulomb_shift
+                operator_2 = number_operator(
+                    n_spin_orbitals, right_neighbor, 1.0) - coulomb_shift
+                hubbard_model += coulomb * operator_1 * operator_2
 
                 # Add hopping term.
                 operators = ((site, 1), (right_neighbor, 0))
@@ -163,9 +180,11 @@ def fermi_hubbard(x_dimension, y_dimension, tunneling, coulomb,
         if site + x_dimension + 1 <= n_sites or (periodic and y_dimension > 2):
             if spinless:
                 # Add Coulomb term.
-                operators = ((site, 1), (site, 0),
-                             (bottom_neighbor, 1), (bottom_neighbor, 0))
-                hubbard_model += FermionOperator(operators, coulomb)
+                operator_1 = number_operator(
+                    n_spin_orbitals, site, 1.0) - coulomb_shift
+                operator_2 = number_operator(
+                    n_spin_orbitals, bottom_neighbor, 1.0) - coulomb_shift
+                hubbard_model += coulomb * operator_1 * operator_2
 
                 # Add hopping term.
                 operators = ((site, 1), (bottom_neighbor, 0))


Test patch:
diff --git a/src/openfermion/hamiltonians/_hubbard_test.py b/src/openfermion/hamiltonians/_hubbard_test.py
index a324b91..ad1d2ba 100644
--- a/src/openfermion/hamiltonians/_hubbard_test.py
+++ b/src/openfermion/hamiltonians/_hubbard_test.py
@@ -32,11 +32,11 @@ class FermiHubbardTest(unittest.TestCase):
 
     def test_two_by_two_spinful(self):
 
-        # Initialize the Hamiltonian.
+        # Initialize the Hamiltonians.
         hubbard_model = fermi_hubbard(
             self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            self.periodic, self.spinless)
+            self.periodic, self.spinless, False)
 
         # Check up spin on site terms.
         self.assertAlmostEqual(hubbard_model.terms[((0, 1), (0, 0))], -.75)
@@ -84,29 +84,74 @@ class FermiHubbardTest(unittest.TestCase):
         self.assertAlmostEqual(hubbard_model.terms[((6, 1), (6, 0),
                                                     (7, 1), (7, 0))], 1.)
 
+    def test_two_by_two_spinful_phs(self):
+        hubbard_model = fermi_hubbard(
+            self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
+            self.chemical_potential, self.magnetic_field,
+            self.periodic, self.spinless, False)
+
+        hubbard_model_phs = fermi_hubbard(
+            self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
+            self.chemical_potential, self.magnetic_field,
+            self.periodic, self.spinless, True)
+
+        # Compute difference between the models with and without
+        # spin symmetry.
+        difference = hubbard_model_phs - hubbard_model
+
+        # Check constant term in difference.
+        self.assertAlmostEqual(difference.terms[()], 1.0)
+
+        # Check up spin on site terms in difference.
+        self.assertAlmostEqual(difference.terms[((0, 1), (0, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((2, 1), (2, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((4, 1), (4, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((6, 1), (6, 0))], -.50)
+
+        # Check down spin on site terms in difference.
+        self.assertAlmostEqual(difference.terms[((1, 1), (1, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((3, 1), (3, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((5, 1), (5, 0))], -.50)
+        self.assertAlmostEqual(difference.terms[((7, 1), (7, 0))], -.50)
+
     def test_two_by_two_spinful_periodic_rudimentary(self):
         hubbard_model = fermi_hubbard(
             self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            periodic=True, spinless=False)
+            periodic=True, spinless=False, particle_hole_symmetry=False)
+
+        hubbard_model = fermi_hubbard(
+            self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
+            self.chemical_potential, self.magnetic_field,
+            periodic=True, spinless=False, particle_hole_symmetry=True)
 
     def test_two_by_two_spinful_aperiodic_rudimentary(self):
         hubbard_model = fermi_hubbard(
             self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            periodic=False, spinless=False)
+            periodic=False, spinless=False, particle_hole_symmetry=False)
+
+        hubbard_model = fermi_hubbard(
+            self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
+            self.chemical_potential, self.magnetic_field,
+            periodic=False, spinless=False, particle_hole_symmetry=True)
 
     def test_two_by_two_spinless_periodic_rudimentary(self):
         hubbard_model = fermi_hubbard(
             self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            periodic=True, spinless=True)
+            periodic=True, spinless=True, particle_hole_symmetry=False)
+
+        hubbard_model = fermi_hubbard(
+            self.x_dimension, self.y_dimension, self.tunneling, self.coulomb,
+            self.chemical_potential, self.magnetic_field,
+            periodic=True, spinless=True, particle_hole_symmetry=True)
 
     def test_two_by_three_spinless_periodic_rudimentary(self):
         hubbard_model = fermi_hubbard(
             2, 3, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            periodic=True, spinless=True)
+            periodic=True, spinless=True, particle_hole_symmetry=False)
         # Check up top/bottom hopping terms.
         self.assertAlmostEqual(hubbard_model.terms[((4, 1), (0, 0))],
                                -self.tunneling)
@@ -115,7 +160,7 @@ class FermiHubbardTest(unittest.TestCase):
         hubbard_model = fermi_hubbard(
             3, 2, self.tunneling, self.coulomb,
             self.chemical_potential, self.magnetic_field,
-            periodic=True, spinless=True)
+            periodic=True, spinless=True, particle_hole_symmetry=False)
         # Check up top/bottom hopping terms.
         self.assertAlmostEqual(hubbard_model.terms[((2, 1), (0, 0))],
                                -self.tunneling)