diff --git a/pyproject.toml b/pyproject.toml index 9c05dec..ea6f23c 100644 --- a/pyproject.toml +++ b/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", diff --git a/qcmanybody/core.py b/qcmanybody/core.py index a89af67..91353b2 100644 --- a/qcmanybody/core.py +++ b/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, ) ) diff --git a/qcmanybody/tests/test_core_singlelevel.py b/qcmanybody/tests/test_core_singlelevel.py index b5a93a2..b6d4299 100644 --- a/qcmanybody/tests/test_core_singlelevel.py +++ b/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() diff --git a/qcmanybody/tests/test_utils.py b/qcmanybody/tests/test_utils.py index 75cdd06..5d1e1b3 100644 --- a/qcmanybody/tests/test_utils.py +++ b/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 diff --git a/qcmanybody/utils.py b/qcmanybody/utils.py index 41677f2..18afd93 100644 --- a/qcmanybody/utils.py +++ b/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