Allow constant in jellium; gitignore profiling files (#119)

* Added code and tests for Madelung constant; profiling files to gitignore

* Increase test coverage

Author: idk3

.gitignore

@@ -70,3 +70,5 @@ ENV/
 
 # mkdocs documentation
 /site
+
+*.lprof

src/fermilib/utils/_jellium.py

@@ -249,14 +249,16 @@ def plane_wave_potential(grid, spinless=False):
 
 
 def dual_basis_jellium_model(grid, spinless=False,
-                             kinetic=True, potential=True):
+                             kinetic=True, potential=True,
+                             include_constant=False):
     """Return jellium Hamiltonian in the dual basis of arXiv:1706.00023
 
     Args:
         grid (Grid): The discretization to use.
         spinless (bool): Whether to use the spinless model or not.
         kinetic (bool): Whether to include kinetic terms.
         potential (bool): Whether to include potential terms.
+        include_constant (bool): Whether to include the Madelung constant.
 
     Returns:
         operator (FermionOperator)
@@ -325,6 +327,10 @@ def dual_basis_jellium_model(grid, spinless=False,
                         operator += FermionOperator(operators,
                                                     potential_coefficient)
 
+    # Include the Madelung constant if requested.
+    if include_constant:
+        operator += FermionOperator.identity() * (2.8372 / grid.scale)
+
     # Return.
     return operator
 
@@ -355,14 +361,16 @@ def dual_basis_potential(grid, spinless=False):
     return dual_basis_jellium_model(grid, spinless, False, True)
 
 
-def jellium_model(grid, spinless=False, plane_wave=True):
+def jellium_model(grid, spinless=False, plane_wave=True,
+                  include_constant=False):
     """Return jellium Hamiltonian as FermionOperator class.
 
     Args:
         grid (fermilib.utils.Grid): The discretization to use.
         spinless (bool): Whether to use the spinless model or not.
         plane_wave (bool): Whether to return in momentum space (True)
             or position space (False).
+        include_constant (bool): Whether to include the Madelung constant.
 
     Returns:
         FermionOperator: The Hamiltonian of the model.
@@ -372,15 +380,20 @@ def jellium_model(grid, spinless=False, plane_wave=True):
         hamiltonian += plane_wave_potential(grid, spinless)
     else:
         hamiltonian = dual_basis_jellium_model(grid, spinless)
+    # Include the Madelung constant if requested.
+    if include_constant:
+        hamiltonian += FermionOperator.identity() * (2.8372 / grid.scale)
     return hamiltonian
 
 
-def jordan_wigner_dual_basis_jellium(grid, spinless=False):
+def jordan_wigner_dual_basis_jellium(grid, spinless=False,
+                                     include_constant=False):
     """Return the jellium Hamiltonian as QubitOperator in the dual basis.
 
     Args:
         grid (Grid): The discretization to use.
         spinless (bool): Whether to use the spinless model or not.
+        include_constant (bool): Whether to include the Madelung constant.
 
     Returns:
         hamiltonian (QubitOperator)
@@ -471,5 +484,9 @@ def jordan_wigner_dual_basis_jellium(grid, spinless=False):
             hamiltonian += QubitOperator(xzx_operators, term_coefficient)
             hamiltonian += QubitOperator(yzy_operators, term_coefficient)
 
+    # Include the Madelung constant if requested.
+    if include_constant:
+        hamiltonian += QubitOperator((),) * (2.8372 / grid.scale)
+
     # Return Hamiltonian.
     return hamiltonian

src/fermilib/utils/_jellium_test.py

@@ -16,9 +16,11 @@
 
 import numpy
 
+from fermilib.ops import FermionOperator
 from fermilib.transforms import jordan_wigner
 from fermilib.utils import count_qubits, eigenspectrum, Grid
 from fermilib.utils._jellium import (
+    dual_basis_jellium_model,
     dual_basis_kinetic,
     dual_basis_potential,
     jellium_model,
@@ -31,6 +33,8 @@
     position_vector,
 )
 
+from projectq.ops import QubitOperator
+
 
 class JelliumTest(unittest.TestCase):
 
@@ -169,6 +173,31 @@ def test_model_integration(self):
             numpy.absolute(momentum_spectrum - position_spectrum))
         self.assertAlmostEqual(difference, 0.)
 
+    def test_model_integration_with_constant(self):
+        # Compute Hamiltonian in both momentum and position space.
+        length_scale = 0.7
+
+        grid = Grid(dimensions=2, length=3, scale=length_scale)
+        spinless = True
+
+        # Include the Madelung constant in the momentum but not the position
+        # Hamiltonian.
+        momentum_hamiltonian = jellium_model(grid, spinless, True,
+                                             include_constant=True)
+        position_hamiltonian = jellium_model(grid, spinless, False)
+
+        # Diagonalize and confirm the same energy.
+        jw_momentum = jordan_wigner(momentum_hamiltonian)
+        jw_position = jordan_wigner(position_hamiltonian)
+        momentum_spectrum = eigenspectrum(jw_momentum)
+        position_spectrum = eigenspectrum(jw_position)
+
+        # Confirm momentum spectrum is shifted 2.8372 / length_scale higher.
+        max_difference = numpy.amax(momentum_spectrum - position_spectrum)
+        min_difference = numpy.amax(momentum_spectrum - position_spectrum)
+        self.assertAlmostEqual(max_difference, 2.8372 / length_scale)
+        self.assertAlmostEqual(min_difference, 2.8372 / length_scale)
+
     def test_coefficients(self):
 
         # Test that the coefficients post-JW transform are as claimed in paper.
@@ -278,9 +307,11 @@ def test_jordan_wigner_dual_basis_jellium(self):
         grid = Grid(dimensions=2, length=3, scale=1.)
         spinless = True
 
-        # Compute fermionic Hamiltonian.
-        fermion_hamiltonian = jellium_model(grid, spinless, False)
+        # Compute fermionic Hamiltonian. Include then subtract constant.
+        fermion_hamiltonian = dual_basis_jellium_model(
+            grid, spinless, include_constant=True)
         qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)
+        qubit_hamiltonian -= QubitOperator((), 2.8372)
 
         # Compute Jordan-Wigner Hamiltonian.
         test_hamiltonian = jordan_wigner_dual_basis_jellium(grid, spinless)
@@ -296,3 +327,18 @@ def test_jordan_wigner_dual_basis_jellium(self):
         num_nonzeros = sum(1 for coeff in qubit_hamiltonian.terms.values() if
                            coeff != 0.0)
         self.assertTrue(num_nonzeros <= paper_n_terms)
+
+    def test_jordan_wigner_dual_basis_jellium_constant_shift(self):
+        length_scale = 0.6
+        grid = Grid(dimensions=2, length=3, scale=length_scale)
+        spinless = True
+
+        hamiltonian_without_constant = jordan_wigner_dual_basis_jellium(
+            grid, spinless, include_constant=False)
+        hamiltonian_with_constant = jordan_wigner_dual_basis_jellium(
+            grid, spinless, include_constant=True)
+
+        difference = hamiltonian_with_constant - hamiltonian_without_constant
+        expected = FermionOperator.identity() * (2.8372 / length_scale)
+
+        self.assertTrue(expected.isclose(difference))

src/fermilib/utils/_plane_wave_hamiltonian.py

@@ -158,7 +158,8 @@ def plane_wave_external_potential(grid, geometry, spinless):
 
 
 def plane_wave_hamiltonian(grid, geometry=None,
-                           spinless=False, plane_wave=True):
+                           spinless=False, plane_wave=True,
+                           include_constant=False):
     """Returns Hamiltonian as FermionOperator class.
 
     Args:
@@ -169,11 +170,12 @@ def plane_wave_hamiltonian(grid, geometry=None,
         spinless (bool): Whether to use the spinless model or not.
         plane_wave (bool): Whether to return in plane wave basis (True)
             or plane wave dual basis (False).
+        include_constant (bool): Whether to include the Madelung constant.
 
     Returns:
         FermionOperator: The hamiltonian.
     """
-    jellium_op = jellium_model(grid, spinless, plane_wave)
+    jellium_op = jellium_model(grid, spinless, plane_wave, include_constant)
 
     if geometry is None:
         return jellium_op
@@ -283,7 +285,8 @@ def _fourier_transform_helper(hamiltonian,
     return hamiltonian_t
 
 
-def jordan_wigner_dual_basis_hamiltonian(grid, geometry=None, spinless=False):
+def jordan_wigner_dual_basis_hamiltonian(grid, geometry=None, spinless=False,
+                                         include_constant=False):
     """Return the dual basis Hamiltonian as QubitOperator.
 
     Args:
@@ -292,11 +295,13 @@ def jordan_wigner_dual_basis_hamiltonian(grid, geometry=None, spinless=False):
             example is [('H', (0, 0, 0)), ('H', (0, 0, 0.7414))].
             Distances in atomic units. Use atomic symbols to specify atoms.
         spinless (bool): Whether to use the spinless model or not.
+        include_constant (bool): Whether to include the Madelung constant.
 
     Returns:
         hamiltonian (QubitOperator)
     """
-    jellium_op = jordan_wigner_dual_basis_jellium(grid, spinless)
+    jellium_op = jordan_wigner_dual_basis_jellium(
+        grid, spinless, include_constant)
 
     if geometry is None:
         return jellium_op

src/fermilib/utils/_plane_wave_hamiltonian_test.py

@@ -117,21 +117,27 @@ def test_plane_wave_hamiltonian_integration(self):
         length_set = [3, 4]
         spinless_set = [True, False]
         geometry = [('H', (0,)), ('H', (0.8,))]
+        length_scale = 1.1
+
         for l in length_set:
             for spinless in spinless_set:
-                grid = Grid(dimensions=1, scale=1.0, length=l)
-                h_plane_wave = plane_wave_hamiltonian(grid, geometry, spinless,
-                                                      True)
+                grid = Grid(dimensions=1, scale=length_scale, length=l)
+                h_plane_wave = plane_wave_hamiltonian(
+                    grid, geometry, spinless, True, include_constant=True)
                 h_dual_basis = plane_wave_hamiltonian(grid, geometry, spinless,
                                                       False)
                 jw_h_plane_wave = jordan_wigner(h_plane_wave)
                 jw_h_dual_basis = jordan_wigner(h_dual_basis)
                 h_plane_wave_spectrum = eigenspectrum(jw_h_plane_wave)
                 h_dual_basis_spectrum = eigenspectrum(jw_h_dual_basis)
 
-                diff = numpy.amax(numpy.absolute(
-                    h_plane_wave_spectrum - h_dual_basis_spectrum))
-                self.assertAlmostEqual(diff, 0)
+                max_diff = numpy.amax(
+                    h_plane_wave_spectrum - h_dual_basis_spectrum)
+                min_diff = numpy.amin(
+                    h_plane_wave_spectrum - h_dual_basis_spectrum)
+
+                self.assertAlmostEqual(max_diff, 2.8372 / length_scale)
+                self.assertAlmostEqual(min_diff, 2.8372 / length_scale)
 
     def test_plane_wave_hamiltonian_default_to_jellium_with_no_geometry(self):
         grid = Grid(dimensions=1, scale=1.0, length=4)
@@ -154,12 +160,12 @@ def test_jordan_wigner_dual_basis_hamiltonian(self):
         spinless = True
         geometry = [('H', (0, 0)), ('H', (0.5, 0.8))]
 
-        fermion_hamiltonian = plane_wave_hamiltonian(grid, geometry, spinless,
-                                                     False)
+        fermion_hamiltonian = plane_wave_hamiltonian(
+            grid, geometry, spinless, False, include_constant=True)
         qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)
 
-        test_hamiltonian = jordan_wigner_dual_basis_hamiltonian(grid, geometry,
-                                                                spinless)
+        test_hamiltonian = jordan_wigner_dual_basis_hamiltonian(
+            grid, geometry, spinless, include_constant=True)
         self.assertTrue(test_hamiltonian.isclose(qubit_hamiltonian))
 
     def test_jordan_wigner_dual_basis_hamiltonian_default_to_jellium(self):