Implement coefficient thresholding in Hamiltonian generation

.github/workflows/python-app.yml

@@ -66,6 +66,12 @@ jobs:
             script: hamgen.py
             config_file: ''
 
+          - name: hamgen-unit-tests
+            directory: hamiltonian_generator
+            script: pytest
+            config_file: ''
+            config_setup: ''
+
           - name: analysis-trotter
             directory: analysis
             script: driver.py

hamiltonian_generator/README.md

@@ -74,6 +74,25 @@ Currently allowed values are
 - for Bravyi-Kitaev: "Bravyi-Kitaev", "Bravyi Kitaev", "BK", "bravyi-kitaev", "bravyi\_kitaev",
   "bravyi kitaev", or "bk"
 
+**`hamiltonian.coefficient_threshold`** controls the absolute cutoff, in Hartrees, used when QHAT
+expands active-space spatial integrals into spin-orbital tensors.  QHAT sets tensor coefficients to
+zero when `abs(value) < coefficient_threshold`, so coefficients exactly equal to the threshold are
+retained.  Its default is `1e-8`; set it to `0.0` to disable filtering at this stage:
+
+```python
+hamiltonian.coefficient_threshold = 0.0
+```
+
+The value must be finite and non-negative.  This setting affects only tensor construction; QHAT
+does not apply an additional coefficient cutoff after OpenFermion's standard JW/BK mapping.  The
+tensor threshold is recorded in the generated `.dat` file metadata.
+
+Active-space pickles record the threshold used to construct their tensors.  If the configured value
+does not match `[filestub]_[astag].pickle`, QHAT stops rather than silently reusing incompatible
+tensors.  Remove that active-space pickle and rerun; QHAT will rebuild it and overwrite
+`[filestub]_[astag].tensors.npz`.  Pickles created before this setting existed are interpreted as
+using the historical `1e-8` cutoff.
+
 ## Generated Files
 
 This script will generate various files for intermediate and final stages of the process.  These

hamiltonian_generator/config.py

@@ -24,6 +24,8 @@
 
 hamiltonian.basis = params.get("basis", "sto-3g") # Select the atomic basis functions
 
+# Absolute cutoff applied while building spin-orbital tensors; the default is 1e-8.
+hamiltonian.coefficient_threshold = 1e-8
+
 hamiltonian.num_active_occupied = 4               # Specify the active space
 hamiltonian.num_active_vacant = 6
-

hamiltonian_generator/hamgen.py

@@ -33,6 +33,10 @@
     HamiltonianConfiguration,
     State,
 )
+from qhat.hamiltonian_generator.thresholding import (
+    DEFAULT_COEFFICIENT_THRESHOLD,
+    spinorb_from_spatial,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -287,9 +291,21 @@ def apply_active_space(state, ham1_HartreeFock):
 
     # Build an InteractionOperator (from openfermion) with one- and two-body integrals computed
     state.log_verbose("Build instance of InteractionOperator.")
-    ham2_ActiveSpace = ham1_HartreeFock.get_molecular_hamiltonian(
-        occupied_indices=range(idx_act_occ),
-        active_indices=range(idx_act_occ, idx_frz_vac),
+    core_adjustment, one_body_integrals, two_body_integrals = (
+        ham1_HartreeFock.get_active_space_integrals(
+            occupied_indices=range(idx_act_occ),
+            active_indices=range(idx_act_occ, idx_frz_vac),
+        )
+    )
+    one_body_coefficients, two_body_coefficients = spinorb_from_spatial(
+        one_body_integrals,
+        two_body_integrals,
+        state.config_hamiltonian.coefficient_threshold,
+    )
+    ham2_ActiveSpace = InteractionOperator(
+        ham1_HartreeFock.nuclear_repulsion + core_adjustment,
+        one_body_coefficients,
+        0.5 * two_body_coefficients,
     )
     tstop = time.time()
     # Pass along metadata that needs to be preserved
@@ -299,6 +315,7 @@ def apply_active_space(state, ham1_HartreeFock):
     ham2_ActiveSpace.separation = ham1_HartreeFock.separation
     ham2_ActiveSpace.hf_time = ham1_HartreeFock.hf_time
     ham2_ActiveSpace.as_time = tstop - tstart
+    ham2_ActiveSpace.coefficient_threshold = state.config_hamiltonian.coefficient_threshold
     return ham2_ActiveSpace
 
 # -------------------------------------------------------------------------------------------------
@@ -315,6 +332,7 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace):
     else:
         mapping = state.config_hamiltonian.f2q_mapping()
         raise NotImplementedError(f"invalid fermion-to-qubit mapping \"{mapping}\"")
+
     tstop = time.time()
     # Propagate previously-computed metadata that we need to preserve across all run modes
     ham3_Fermion2Qubit.hf_energy = ham2_ActiveSpace.hf_energy
@@ -325,6 +343,7 @@ def map_fermions_to_qubits(state, ham2_ActiveSpace):
     ham3_Fermion2Qubit.hf_time = ham2_ActiveSpace.hf_time
     ham3_Fermion2Qubit.as_time = ham2_ActiveSpace.as_time
     ham3_Fermion2Qubit.f2q_time = tstop - tstart
+
     # Return result
     return ham3_Fermion2Qubit
 
@@ -437,6 +456,18 @@ def get_ham2(state):
                             two_body=ham2_ActiveSpace.two_body_tensor)
         return ham2_ActiveSpace
     else:
+        cached_threshold = getattr(
+            ham2_ActiveSpace,
+            "coefficient_threshold",
+            DEFAULT_COEFFICIENT_THRESHOLD,
+        )
+        configured_threshold = state.config_hamiltonian.coefficient_threshold
+        if cached_threshold != configured_threshold:
+            raise ValueError(
+                f"Active-space cache \"{ham2_filename}\" was built with coefficient threshold "
+                f"{cached_threshold}, but the configuration requests {configured_threshold}; "
+                "remove the cache and rerun."
+            )
         state.log(' '.join([
             f"Loaded \"{ham2_cache_path}\".",
             "Continuing from after the active space is applied."]))
@@ -498,6 +529,9 @@ def compute_metadata(state, ham3_Fermion2Qubit):
     state.metadata["interatomic separation (angstroms)"] = ham3_Fermion2Qubit.separation
     # fermion-to-qubit transformation
     state.metadata["fermion-to-qubit operator mapping"] = ham3_Fermion2Qubit.f2q_mapping
+    state.metadata["spin-orbital coefficient threshold (Hartrees)"] = (
+        state.config_hamiltonian.coefficient_threshold
+    )
     # number of terms in sum of Pauli strings
     state.metadata["number of terms in sum of Pauli strings"] = len(ham3_Fermion2Qubit.terms)
     # one-norm of sum of Pauli strings

hamiltonian_generator/hamgen_types.py

@@ -7,6 +7,10 @@
 import sys
 
 from qhat.common.git_utils import get_git_hash
+from qhat.hamiltonian_generator.thresholding import (
+    DEFAULT_COEFFICIENT_THRESHOLD,
+    validate_coefficient_threshold,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -40,9 +44,16 @@ def __init__(self):
         self._geometry = ()
         self.basis = "sto-3g"
         self.f2q_mapping = "Jordan-Wigner"
+        self._coefficient_threshold = DEFAULT_COEFFICIENT_THRESHOLD
         # The number of active occupied and vacant spin (single-occupancy) orbitals in the molecule
         self.num_active_occupied = None
         self.num_active_vacant = None
+    @property
+    def coefficient_threshold(self):
+        return self._coefficient_threshold
+    @coefficient_threshold.setter
+    def coefficient_threshold(self, value):
+        self._coefficient_threshold = validate_coefficient_threshold(value)
     def add_atom(self, element, x, y, z):
         # TODO: Clarify the units involved.  Right now we just pass numbers along, so the user has
         #       to know the right units for the various choices they make.  It would be better to
@@ -208,4 +219,3 @@ def filename_ham3(self):
                 ext=self.config_general.ham3_ext())
 
 # -------------------------------------------------------------------------------------------------
-