Bugfix: Add ECP Contributions

cpp/src/qdk/chemistry/algorithms/microsoft/cholesky_hamiltonian.cpp

@@ -646,9 +646,8 @@ std::shared_ptr<data::Hamiltonian> CholeskyHamiltonianConstructor::_run_impl(
 
   // Create internal Molecule
   auto structure = basis_set->get_structure();
-  auto mol = utils::microsoft::convert_to_molecule(*structure, 0, 1);
 
-  // Create internal BasisSet
+  // Create internal BasisSet (includes ECP-adjusted nuclear charges)
   auto internal_basis_set =
       utils::microsoft::convert_basis_set_from_qdk(*basis_set);
   // Create dummy SCFConfig
@@ -666,7 +665,7 @@ std::shared_ptr<data::Hamiltonian> CholeskyHamiltonianConstructor::_run_impl(
   // Create Integral Instance
   auto eri = qcs::ERIMultiplexer::create(*internal_basis_set, *scf_config, 0.0);
   auto int1e = std::make_unique<qcs::OneBodyIntegral>(
-      internal_basis_set.get(), mol.get(), scf_config->mpi);
+      internal_basis_set.get(), internal_basis_set->mol.get(), scf_config->mpi);
 
   // Compute Core Hamiltonian in AO basis
   Eigen::MatrixXd T_full(num_atomic_orbitals, num_atomic_orbitals),
@@ -675,6 +674,14 @@ std::shared_ptr<data::Hamiltonian> CholeskyHamiltonianConstructor::_run_impl(
   int1e->nuclear_integral(V_full.data());
   Eigen::MatrixXd H_full = T_full + V_full;
 
+  // Add ECP integrals if present
+  if (internal_basis_set->ecp_shells.size() > 0) {
+    Eigen::MatrixXd ECP_full =
+        Eigen::MatrixXd::Zero(num_atomic_orbitals, num_atomic_orbitals);
+    int1e->ecp_integral(ECP_full.data());
+    H_full += ECP_full;
+  }
+
   // Build active coefficient matrices for alpha and beta (can have different
   // sizes)
   Eigen::MatrixXd Ca_active(num_atomic_orbitals, nactive_alpha);

cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp

@@ -138,9 +138,8 @@ std::shared_ptr<data::Hamiltonian> HamiltonianConstructor::_run_impl(
 
   // Create internal Molecule
   auto structure = basis_set->get_structure();
-  auto mol = utils::microsoft::convert_to_molecule(*structure, 0, 1);
 
-  // Create internal BasisSet
+  // Create internal BasisSet (includes ECP-adjusted nuclear charges)
   auto internal_basis_set =
       utils::microsoft::convert_basis_set_from_qdk(*basis_set);
   // Create dummy SCFConfig
@@ -169,7 +168,7 @@ std::shared_ptr<data::Hamiltonian> HamiltonianConstructor::_run_impl(
   // Create Integral Instance
   auto eri = qcs::ERIMultiplexer::create(*internal_basis_set, *scf_config, 0.0);
   auto int1e = std::make_unique<qcs::OneBodyIntegral>(
-      internal_basis_set.get(), mol.get(), scf_config->mpi);
+      internal_basis_set.get(), internal_basis_set->mol.get(), scf_config->mpi);
 
   // Compute Core Hamiltonian in AO basis
   Eigen::MatrixXd T_full(num_atomic_orbitals, num_atomic_orbitals),
@@ -178,6 +177,14 @@ std::shared_ptr<data::Hamiltonian> HamiltonianConstructor::_run_impl(
   int1e->nuclear_integral(V_full.data());
   Eigen::MatrixXd H_full = T_full + V_full;
 
+  // Add ECP integrals if present
+  if (internal_basis_set->ecp_shells.size() > 0) {
+    Eigen::MatrixXd ECP_full =
+        Eigen::MatrixXd::Zero(num_atomic_orbitals, num_atomic_orbitals);
+    int1e->ecp_integral(ECP_full.data());
+    H_full += ECP_full;
+  }
+
   // Build active coefficient matrices for alpha and beta (can have different
   // sizes)
   Eigen::MatrixXd Ca_active(num_atomic_orbitals, nactive_alpha);

python/tests/test_ecp_hamiltonian.py

@@ -0,0 +1,48 @@
+"""Verify that ECP contributions are included in Hamiltonian one-electron integrals.
+
+Bug: HamiltonianConstructor.run() computes h_core = T + V_nuc but omits V_ecp.
+This causes incorrect one-electron integrals for any atom using an ECP basis
+(e.g., Mo with def2-svp).
+"""
+
+# --------------------------------------------------------------------------------------------
+# Copyright (c) Microsoft Corporation. All rights reserved.
+# Licensed under the MIT License. See LICENSE.txt in the project root for license information.
+# --------------------------------------------------------------------------------------------
+
+import numpy as np
+import pytest
+
+from qdk_chemistry.algorithms import create
+from qdk_chemistry.data import Structure
+
+from .reference_tolerances import scf_energy_tolerance
+
+# Mo ROHF/def2-svp reference: trace of h1e (alpha = beta for ROHF).
+# Trace is unitarily invariant, so it does not depend on the MO rotation
+# within degenerate subspaces (d-shell degeneracy causes different rotations
+# across thread counts / platforms).
+# Without ECP (bare Z=42), the trace is ~ -1800.  With ECP (Z_eff=14), ~ -134.
+_MO_ROHF_H1E_TRACE = -133.9900221637
+
+
+@pytest.fixture(scope="module")
+def mo_rohf_orbitals():
+    """ROHF orbitals for Mo atom with def2-svp (shared across both parameterized tests)."""
+    structure = Structure(np.array([[0.0, 0.0, 0.0]]), ["Mo"])
+    scf = create("scf_solver", "qdk")
+    scf.settings()["method"] = "hf"
+    scf.settings()["scf_type"] = "restricted"
+    scf.settings()["enable_gdm"] = False
+    _, wfn = scf.run(structure, 0, 7, "def2-svp")
+    return wfn.get_orbitals()
+
+
+@pytest.mark.parametrize("constructor_name", ["qdk", "qdk_cholesky"])
+def test_ecp_included_in_hamiltonian_h1e(mo_rohf_orbitals, constructor_name):
+    """h1e from HamiltonianConstructor must include ECP terms for Mo/def2-svp."""
+    ham = create("hamiltonian_constructor", constructor_name).run(mo_rohf_orbitals)
+    h1e_a, h1e_b = ham.get_one_body_integrals()
+
+    np.testing.assert_allclose(np.trace(h1e_a), _MO_ROHF_H1E_TRACE, atol=scf_energy_tolerance)
+    np.testing.assert_allclose(np.trace(h1e_b), _MO_ROHF_H1E_TRACE, atol=scf_energy_tolerance)