Fix an edge-case bug in gradient of transformed Hamiltonians

evangelistalab · Nov 19, 2024 · e9992cb · e9992cb
1 parent b77f054
commit e9992cb
Show file tree

Hide file tree

Showing 4 changed files with 16 additions and 12 deletions.
diff --git a/forte/apps/hf.py b/forte/apps/hf.py
@@ -12,7 +12,7 @@ def run_hf(geom, basis, state):
         [
             OptionsFactory(),
             MoleculeFactory(geom),
-            StateFactory(charge=0, multiplicity=1, sym="ag"),
+            StateFactory(charge=charge, multiplicity=multiplicity, sym=sym),
             HF(basis=basis),
         ],
         name="HF Workflow",

diff --git a/forte/modules/mcscf.py b/forte/modules/mcscf.py
@@ -33,12 +33,7 @@ def _run(self, data: ForteData) -> ForteData:
             self.solver_type, state_map, data.scf_info, data.mo_space_info, data.options
         )
         mcscf = make_mcscf_two_step(
-            data.active_space_solver,
-            data.state_weights_map,
-            data.scf_info,
-            data.options,
-            data.mo_space_info,
-            data.ints
+            data.active_space_solver, data.state_weights_map, data.scf_info, data.options, data.mo_space_info, data.ints
         )
         energy = mcscf.compute_energy()
         data.results.add("mcscf energy", energy, "MCSCF energy", "hartree")

diff --git a/forte/sparse_ci/sparse_operator_sim_trans.cc b/forte/sparse_ci/sparse_operator_sim_trans.cc
@@ -91,15 +91,14 @@ void fact_unitary_trans_antiherm_grad(SparseOperator& O, const SparseOperatorLis
             throw std::runtime_error("sim_trans_antiherm_impl: the angle theta must be real.");
         }
 
-        // skip transformation if theta is zero
-        if (std::abs(theta) < screen_threshold) {
+        // skip transformation if theta is zero, except if we are taking the gradient
+        if ((std::abs(theta) < screen_threshold) and index != grad_index) {
             return;
         }
 
         // if T = T^dagger, then the transformation is trivial, so we can skip it
         if (sqop.is_self_adjoint()) {
-            // if the operator is the one for which the gradient is being computed, then set it to
-            // zero
+            // if we are taking the gradient of this operator then the result is zero
             if (index == grad_index) {
                 O = SparseOperator();
             }

diff --git a/tests/pytest/sparse_ci/test_sparse_operator_complex.py b/tests/pytest/sparse_ci/test_sparse_operator_complex.py
@@ -210,14 +210,22 @@ def test_sparse_operator_complex_transform5():
     run_sparse_operator_test(O, theta, op)
 
 
-def test_sparse_operator_complex_transform5():
+def test_sparse_operator_complex_transform6():
     """test exponentiation of a sparse operator numerically"""
     O = make_O()
     op = "[0a+ 1a+ 2a+ 2a- 1a- 0a-]"
     theta = 0.35711
     run_sparse_operator_test(O, theta, op)
 
 
+def test_sparse_operator_complex_transform7():
+    """Test exponentiation of a sparse operator numerically with theta = 0 (edge case)"""
+    O = make_O()
+    op = "[1a+ 0a-]"
+    theta = 0.0
+    run_sparse_operator_test(O, theta, op)
+
+
 def run_sparse_operator_test(O, theta, op):
     """test exponentiation of a single sparse operator numerically"""
     sqop, _ = forte.sqop(op)
@@ -275,3 +283,5 @@ def run_sparse_operator_test(O, theta, op):
     test_sparse_operator_complex_transform3()
     test_sparse_operator_complex_transform4()
     test_sparse_operator_complex_transform5()
+    test_sparse_operator_complex_transform6()
+    test_sparse_operator_complex_transform7()