Constant energy bug (#107)

Author: MIWdlB

README.md

@@ -51,11 +51,13 @@ For more details, see the [interop documentation](https://ferrmion.readthedocs.i
 
 
 ## Hamiltonians
-Functions to produce Hamiltonian templates and enumerated encodings are available in `ferrmion.hamiltonians`
+Arbitrary fermionic Hamiltonians can be built with the `FermionHamiltonian` class in `ferrmion.hamiltonians`. Additionally there are named functions for:
 
 - Molecular Hamiltonian (Born-Oppenheimer Approximation)
 - Hubbard Hamiltonian (3D Square lattice)
 
+See the [example notebook](https://ferrmion.readthedocs.io/en/latest/notebooks/encoding_hamiltonians.html#) for more information.
+
 ## Development
 For information on development, check out the project [documentation](https://ferrmion.readthedocs.io/en/latest/development.html).
 

python/ferrmion/core.pyi

@@ -78,7 +78,7 @@ def encode_standard(
     n_qubits: int,
     signatures: list[str],
     coeffs: list[np.ndarray],
-    constant_energy:float,
+    constant_energy: float,
 ) -> dict: ...
 def encode(
     ipowers: npt.NDArray[np.uint8],

python/tests/test_hamiltonians.py

@@ -34,6 +34,7 @@ def test_encode_standard_water_eigvals_equal_expected(encoding, water_data):
     e_nuc = water_data["constant_energy"]
 
     qham = encode_standard(encoding, 14,14, ["+-","++--"], [ones, twos], e_nuc)
+    assert np.isclose(qham["I"*14], -46.465600781952176)
 
     ofop = QubitOperator()
     for k, v in qham.items():

python/tests/test_ternary_tree.py

@@ -1,3 +1,4 @@
+from ferrmion import FermionHamiltonian
 import numpy as np
 import pytest
 import scipy as sp
@@ -17,6 +18,7 @@
 from openfermion.transforms import jordan_wigner
 from ferrmion.hamiltonians import molecular_hamiltonian
 from ferrmion.core import standard_symplectic_matrix
+from scipy.sparse.linalg import eigsh
 
 
 @pytest.fixture
@@ -573,3 +575,50 @@ def test_core_standard_encodings(n_modes,encoding,name):
     ci, cs = standard_symplectic_matrix(name,20)
     assert np.all(i==ci)
     assert np.all(s==cs)
+
+@pytest.mark.parametrize("encoding", [JW, BK, ParityEncoding, JKMN])
+def test_encode_standard_water_eigvals_equal_expected(encoding, water_data):
+    ones = water_data["ones"]
+    twos = water_data["twos"]
+    e_nuc = water_data["constant_energy"]
+
+    fham = FermionHamiltonian(terms = {"+-":ones,"++--":twos}, constant_energy=e_nuc)
+    qham = encoding(fham.n_modes).encode(fham)
+    assert np.isclose(qham["I"*14], -46.465600781952176)
+
+    ofop = QubitOperator()
+    for k, v in qham.items():
+        string = " ".join(
+            [
+                f"{char.upper()}{pos}" if char != "I" else ""
+                for pos, char in enumerate(k)
+            ]
+        )
+        ofop+= QubitOperator(term=string, coefficient=v)
+    print(expected:=water_data["eigvals"])
+    diag, _ = eigsh(get_sparse_operator(ofop), k=2, which="SA")
+    print(diag)
+    assert np.allclose(np.sort(diag), np.sort(expected)[:2])
+
+@pytest.mark.parametrize("encoding", [JW, BK, ParityEncoding, JKMN])
+def test_encode_standard_h2_eigvals_equal_expected(encoding, h2_mol_data_sets):
+    ones = h2_mol_data_sets["ones"]
+    twos = h2_mol_data_sets["twos"]
+    e_nuc = h2_mol_data_sets["constant_energy"]
+    n_modes = ones.shape[0]
+    fham = FermionHamiltonian(terms = {"+-":ones,"++--":twos}, constant_energy=e_nuc)
+    qham = encoding(n_modes).encode(fham)
+
+    ofop = QubitOperator()
+    for k, v in qham.items():
+        string = " ".join(
+            [
+                f"{char.upper()}{pos}" if char != "I" else ""
+                for pos, char in enumerate(k)
+            ]
+        )
+        ofop+= QubitOperator(term=string, coefficient=v)
+    print(expected:=h2_mol_data_sets["eigvals"])
+    diag, _ = eigsh(get_sparse_operator(ofop), k=2*n_modes, which="SA")
+    print(diag)
+    assert np.allclose(np.sort(diag), np.sort(h2_mol_data_sets["eigvals"]))

src/encoding.rs

@@ -269,7 +269,7 @@ impl Encode<&MajoranaSparse> for MajoranaEncoding {
             )
             .or_insert(c64(0., 0.)) += input.constant;
         // qham.into_iter().filter(|(_, v)| v.abs() > 1e-12).collect()
-        qham.into_iter().filter(|(k, v)| v.norm() > 1e-16).collect()
+        qham.into_iter().filter(|(_, v)| v.norm() > 1e-16).collect()
     }
 }
 

src/lib.rs

@@ -1,7 +1,8 @@
 use ::core::panic;
 use log::{debug, info};
 use numpy::{
-    Complex64, IntoPyArray, PyArray1, PyArray2, PyArray3, PyReadonlyArray1, PyReadonlyArray2, PyReadonlyArrayDyn,
+    Complex64, IntoPyArray, PyArray1, PyArray2, PyArray3, PyReadonlyArray1, PyReadonlyArray2,
+    PyReadonlyArrayDyn,
 };
 use pyo3::types::{IntoPyDict, PyComplex, PyDict, PyInt, PyString};
 use pyo3::{prelude::*, pymodule, Bound};
@@ -362,7 +363,8 @@ fn core(m: &Bound<'_, PyModule>) -> PyResult<()> {
             "Must have at least as many qubits as modes."
         );
 
-        let hamiltonian: MajoranaSparse = build_majorana_sparse(signatures, coeffs, constant_energy);
+        let hamiltonian: MajoranaSparse =
+            build_majorana_sparse(signatures, coeffs, constant_energy);
 
         debug!("Got MSparse");
         debug!("Got Hamiltonian");
@@ -379,7 +381,7 @@ fn core(m: &Bound<'_, PyModule>) -> PyResult<()> {
         let encoding = tree.build_encoding(n_qubits).unwrap();
         debug!("Got encoding {:?}", encoding);
         debug!("Got encoding {:?}", encoding);
-        
+
         let qham: QubitHamiltonian = encoding.encode(&hamiltonian);
 
         debug!("Got qham");
@@ -453,7 +455,7 @@ fn core(m: &Bound<'_, PyModule>) -> PyResult<()> {
         debug!("Starting TOPPHATT");
         // let flatpack: TTFlatPack = node_map.extract::<TTFlatPack>()?;
         let hamiltonian: MajoranaSparse = build_majorana_sparse(signatures, coeffs, 0.);
-        
+
         debug!("Got MSparse");
         debug!("Got Hamiltonian");
         let mut tree: TernaryTree = match encoding.as_str() {

src/operators.rs

@@ -1,18 +1,16 @@
 use crate::ternarytree::Edge;
 use ndarray::Dimension;
-use tinyvec::{ArrayVec};
+use tinyvec::ArrayVec;
 /*
 Shared Types.
 */
 use itertools::Itertools;
 use log::debug;
 use num_complex::{c64, ComplexFloat};
-use numpy::ndarray::{
-    arr1, arr2, Array1, Array2, ArrayD, ArrayView1, Axis, IntoDimension, Zip,
-};
+use numpy::ndarray::{arr1, arr2, Array1, Array2, ArrayD, ArrayView1, Axis, IntoDimension, Zip};
 use numpy::Complex64;
 use std::collections::BTreeMap;
-use std::iter::{repeat_n};
+use std::iter::repeat_n;
 use std::{result::Result, str::FromStr};
 
 const MAX_MAJORANAS: usize = 4;
@@ -598,23 +596,17 @@ impl MajoranaSparse {
         if coefficients.len() != indices.len() {
             return Err(MajoranaSparseError);
         };
-        let identity_terms_constant: Complex64 = indices
-            .iter()
-            .zip(&coefficients)
-            .filter(|&(ind, _)| ind == &ArrayVec::new())
-            .map(|(_, coeff)| coeff)
-            .sum();
 
         let (i, c) = indices
             .iter()
             .zip(&coefficients)
-            .filter(|&(&inds, &coeff)| (coeff != Complex64::ZERO) && (inds != ArrayVec::new()))
+            .filter(|&(_, &coeff)| (coeff != Complex64::ZERO))
             .unzip();
 
         Ok(Self {
             indices: i,
             coefficients: c,
-            constant: constant + identity_terms_constant.norm(),
+            constant: constant,
         })
     }
 }