feat(qcmanybody): support fragment-scoped external potentials

pyproject.toml

@@ -86,10 +86,10 @@ default-tag = "0.0.99" # useful for CI/shallow clones
 minversion = "7.0"
 addopts = "--import-mode=importlib"
 pythonpath = [
-    "tests",
+    "qcmanybody/tests",
 ]
 testpaths = [
-    "tests",
+    "qcmanybody/tests",
 ]
 markers = [
     "addon: tests require external non-required software",

qcmanybody/core.py

@@ -526,6 +526,8 @@ def _analyze(
     def analyze(
         self,
         component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],
+        *,
+        external_potential: bool = False,
     ):
         """
 
@@ -548,6 +550,10 @@ def analyze(
                '["mp2", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},
                '["mp2", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},
               }
+        external_potential
+            Whether the component calculations use fragment-scoped external potentials. Because the component
+            Hamiltonians differ, interaction energies are unavailable in the printed summary and interaction-energy
+            properties are omitted from the returned results. This classification is distinct from charge embedding.
 
         """
 
@@ -602,6 +608,7 @@ def analyze(
                 is_embedded,
                 self.supersystem_ie_only,
                 self.max_nbody if self.has_supersystem else None,
+                external_potential=external_potential,
             )
 
         for property_label in available_properties:
@@ -615,6 +622,7 @@ def analyze(
                         is_embedded,
                         self.supersystem_ie_only,
                         self.has_supersystem,
+                        external_potential=external_potential,
                     )
                 )
 

qcmanybody/tests/test_core_singlelevel.py

@@ -33,6 +33,26 @@ def test_h2o_trimer_single(levels, component_file, ref_file):
     compare_results(nbody_results, ref_data, levels)
 
 
+def test_external_potential_reporting():
+    levels = {1: "e_scf", 2: "e_scf", 3: "e_scf"}
+    component_results = load_component_data("h2o_trimer_single_energy_1")
+    mc = ManyBodyCore(
+        mol_h2o_3_dict,
+        [BsseEnum.cp],
+        levels,
+        embedding_charges=None,
+        supersystem_ie_only=False,
+        return_total_data=True,
+    )
+
+    nbody_results = mc.analyze(component_results, external_potential=True)
+
+    assert "cp_corrected_2_body_contribution_to_energy" in nbody_results["results"]
+    assert not any("interaction" in name for name in nbody_results["results"])
+    assert "N/A" in nbody_results["stdout"]
+    assert "nan" not in nbody_results["stdout"].lower()
+
+
 def test_core_mol_error():
     # check sensible error from internal Molecule construction
     odd_mol = mol_h2o_3_dict.copy()

qcmanybody/tests/test_utils.py

@@ -136,6 +136,92 @@ def test_resize_gradient(gin, bas, reverse, gans):
     assert compare_values(gans, gout, atol=1e-5, label="resize_gradient")
 
 
+def test_print_nbody_energy_external_potential_without_total_energy():
+    stdout = qcmb.utils.print_nbody_energy(
+        {1: 0.0, 2: 0.005150000059},
+        "Counterpoise Corrected (CP)",
+        2,
+        {1: ("model", "§A"), 2: ("model", "§A")},
+        embedding=False,
+        supersystem_ie_only=False,
+        supersystem_beyond=None,
+        external_potential=True,
+    )
+
+    rows = [line for line in stdout.splitlines() if "§A" in line and "N/A" in line]
+    assert rows == [
+        "              §A  1        N/A                              N/A                   N/A"
+        "        0.000000000000        0.000000000000",
+        "     FULL/RTN §A  2        N/A                              N/A                   N/A"
+        "        0.005150000059        3.231673826977",
+    ]
+    assert "nan" not in stdout.lower()
+
+    properties = qcmb.utils.collect_vars(
+        "cp", "energy", {1: 0.0, 2: 0.005150000059}, 2, external_potential=True
+    )
+    assert "cp_corrected_2_body_contribution_to_energy" in properties
+    assert not any("interaction" in name for name in properties)
+
+
+def test_print_nbody_energy_preserves_legacy_spacing():
+    labels = {1: ("model", "§A"), 2: ("model", "§A")}
+
+    interaction_only = qcmb.utils.print_nbody_energy(
+        {1: 0.0, 2: 0.005150000059},
+        "Counterpoise Corrected (CP)",
+        2,
+        labels,
+        embedding=False,
+        supersystem_ie_only=False,
+        supersystem_beyond=None,
+    )
+    interaction_rows = [line for line in interaction_only.splitlines() if "§A" in line and "Legend" not in line]
+    assert interaction_rows == [
+        "              §A  1        N/A                   0.000000000000        0.000000000000"
+        "        0.000000000000        0.000000000000",
+        "     FULL/RTN §A  2        N/A                   0.005150000059        3.231673826977"
+        "        0.005150000059        3.231673826977",
+    ]
+
+    supersystem = qcmb.utils.print_nbody_energy(
+        {1: -1.0, 2: -0.9},
+        "Counterpoise Corrected (CP)",
+        2,
+        labels,
+        embedding=False,
+        supersystem_ie_only=True,
+        supersystem_beyond=None,
+    )
+    supersystem_rows = [line for line in supersystem.splitlines() if "§A" in line and "Legend" not in line]
+    assert supersystem_rows == [
+        "              §A  1       -1.000000000000        0.000000000000        0.000000000000"
+        "        0.000000000000        0.000000000000",
+        "     FULL/RTN §A  2       -0.900000000000        0.100000000000       62.750947377754"
+        "        N/A                   N/A                 ",
+    ]
+
+    external_supersystem = qcmb.utils.print_nbody_energy(
+        {1: -1.0, 2: -0.9},
+        "Counterpoise Corrected (CP)",
+        2,
+        labels,
+        embedding=False,
+        supersystem_ie_only=True,
+        supersystem_beyond=None,
+        external_potential=True,
+    )
+    external_supersystem_rows = [
+        line for line in external_supersystem.splitlines() if "§A" in line and "Legend" not in line
+    ]
+    assert external_supersystem_rows == [
+        "              §A  1       -1.000000000000        N/A                   N/A"
+        "                   0.000000000000        0.000000000000",
+        "     FULL/RTN §A  2       -0.900000000000                   N/A                   N/A"
+        "        N/A                   N/A                 ",
+    ]
+
+
 @pytest.mark.parametrize("hin,bas,reverse,hans", [
     pytest.param(f3hesses["full"], [1, 2, 3], False, f3hesses["full"]),  # idempotent
     pytest.param(f3hesses["full"], [1, 2, 3], True, f3hesses["full"]),  # idempotent

qcmanybody/utils.py

@@ -322,6 +322,7 @@ def print_nbody_energy(
     embedding: bool,
     supersystem_ie_only: bool,
     supersystem_beyond: Optional[int],
+    external_potential: bool = False,
 ) -> str:
     """Format summary string for energies of a single bsse_type. Logs and returns output.
 
@@ -338,11 +339,13 @@ def print_nbody_energy(
         full model chemistry key and second element a short label. A suitable
         dictionary is `modelchem_labels(manybodycore_instance.nbodies_per_mc_level)`.
     embedding
-        Whether charge embedding present suppress printing, usually False
+        Whether charge embedding is present and suppresses interaction energies, usually False.
     supersystem_ie_only
         Whether only 1-body and nfragments-body levels are available, usually False.
     supersystem_beyond
         If not None, the number of nbody-levels computed by MBE explicitly. Beyond this gets supersystem SS label.
+    external_potential
+        Whether fragment-scoped external potentials are present and suppress interaction energies, usually False.
 
     Returns
     -------
@@ -366,6 +369,10 @@ def print_nbody_energy(
     info += f"""                         [Eh]                    [Eh]                  [kcal/mol]            [Eh]                  [kcal/mol]\n"""
     previous_e = energy_body_dict[1]
     tot_e = previous_e != 0.0
+    interaction_data_unavailable = embedding or external_potential
+    # Preserve the historical numeric-field offsets while replacing only the
+    # literal unavailable ``nan`` token with ``N/A``.
+    unavailable_interaction = f"{np.nan:20.12f}".replace("nan", "N/A")
     nbody_range = list(energy_body_dict)
     if supersystem_ie_only:
         nbody_range = [1, nfragments]
@@ -385,25 +392,31 @@ def print_nbody_energy(
         if nb in nbody_range:
             delta_e = energy_body_dict[nb] - previous_e
             delta_e_kcal = delta_e * constants.hartree2kcalmol
-            if embedding:
-                int_e = np.nan
-                int_e_kcal = np.nan
-            else:
+            if not interaction_data_unavailable:
                 int_e = energy_body_dict[nb] - energy_body_dict[1]
                 int_e_kcal = int_e * constants.hartree2kcalmol
             if supersystem_ie_only and nb == nfragments:
                 if tot_e:
-                    info += f"""  {lbl:>11} {mclbl:2} {nb:2}  {energy_body_dict[nb]:20.12f}  {int_e:20.12f}  {int_e_kcal:20.12f}        {"N/A":20}  {"N/A":20}\n"""
+                    if interaction_data_unavailable:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}  {energy_body_dict[nb]:20.12f}  {unavailable_interaction}  {unavailable_interaction}        {"N/A":20}  {"N/A":20}\n"""
+                    else:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}  {energy_body_dict[nb]:20.12f}  {int_e:20.12f}  {int_e_kcal:20.12f}        {"N/A":20}  {"N/A":20}\n"""
                 else:
-                    info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {int_e:20.12f}  {int_e_kcal:20.12f}        {"N/A":20}  {"N/A":20}\n"""
+                    if interaction_data_unavailable:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {unavailable_interaction}  {unavailable_interaction}        {"N/A":20}  {"N/A":20}\n"""
+                    else:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {int_e:20.12f}  {int_e_kcal:20.12f}        {"N/A":20}  {"N/A":20}\n"""
             else:
                 if tot_e:
-                    if embedding:
+                    if interaction_data_unavailable:
                         info += f"""  {lbl:>11} {mclbl:2} {nb:2}  {energy_body_dict[nb]:20.12f}        {"N/A":20}  {"N/A":14}  {delta_e:20.12f}  {delta_e_kcal:20.12f}\n"""
                     else:
                         info += f"""  {lbl:>11} {mclbl:2} {nb:2}  {energy_body_dict[nb]:20.12f}  {int_e:20.12f}  {int_e_kcal:20.12f}  {delta_e:20.12f}  {delta_e_kcal:20.12f}\n"""
                 else:
-                    info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {int_e:20.12f}  {int_e_kcal:20.12f}  {delta_e:20.12f}  {delta_e_kcal:20.12f}\n"""
+                    if interaction_data_unavailable:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {unavailable_interaction}  {unavailable_interaction}  {delta_e:20.12f}  {delta_e_kcal:20.12f}\n"""
+                    else:
+                        info += f"""  {lbl:>11} {mclbl:2} {nb:2}        {"N/A":14}  {int_e:20.12f}  {int_e_kcal:20.12f}  {delta_e:20.12f}  {delta_e_kcal:20.12f}\n"""
             previous_e = energy_body_dict[nb]
         else:
             info += f"""  {lbl:>11} {"":2} {nb:2}        {"N/A":20}  {"N/A":20}  {"N/A":20}  {"N/A":20}  {"N/A":20}\n"""
@@ -422,6 +435,7 @@ def collect_vars(
     embedding: bool = False,
     supersystem_ie_only: bool = False,
     has_supersystem: bool = False,
+    external_potential: bool = False,
 ) -> Dict:
     """From *body_dict*, construct data for ManyBodyResultProperties/ManyBodyProperties.
 
@@ -445,11 +459,15 @@ def collect_vars(
         By default False: data is available for consecutive levels, up to max_nbody-body.
     has_supersystem
         Whether contributions higher than max_nbody are a summary correction.
+    external_potential
+        Whether fragment-scoped external potentials are present. Interaction properties are unavailable because the
+        component Hamiltonians differ.
 
     Returns
     -------
     dict
-        _description_. Empty return if *embedding* enabled.
+        Many-body result properties. Interaction properties are omitted when *embedding* or *external_potential* is
+        enabled; valid total and n-body contribution properties remain.
     """
     bsse = bsse.lower()
     prop = prop.lower()
@@ -495,7 +513,7 @@ def collect_vars(
             for nb in nbody_range:
                 res[f"{bsse}_corrected_total_{prop}_through_{nb}_body"] = body_dict[nb]
 
-    if embedding:
+    if embedding or external_potential:
         res = {k: v for k, v in res.items() if "interaction" not in k}
 
     return res