diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/cholesky_hamiltonian.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/cholesky_hamiltonian.cpp index fba412253..3d47dd6a5 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/cholesky_hamiltonian.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/cholesky_hamiltonian.cpp @@ -646,9 +646,8 @@ std::shared_ptr 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 CholeskyHamiltonianConstructor::_run_impl( // Create Integral Instance auto eri = qcs::ERIMultiplexer::create(*internal_basis_set, *scf_config, 0.0); auto int1e = std::make_unique( - 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 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); diff --git a/cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp b/cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp index 22521f585..a7d3e2018 100644 --- a/cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp +++ b/cpp/src/qdk/chemistry/algorithms/microsoft/hamiltonian.cpp @@ -138,9 +138,8 @@ std::shared_ptr 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 HamiltonianConstructor::_run_impl( // Create Integral Instance auto eri = qcs::ERIMultiplexer::create(*internal_basis_set, *scf_config, 0.0); auto int1e = std::make_unique( - 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 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); diff --git a/python/tests/test_ecp_hamiltonian.py b/python/tests/test_ecp_hamiltonian.py new file mode 100644 index 000000000..e148a0cbd --- /dev/null +++ b/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)