more misc

qcmanybody/manybody.py

@@ -2,6 +2,7 @@
 
 import logging
 import math
+from collections import defaultdict
 from typing import Set, Dict, Tuple, Union, Literal, Mapping, Any, Sequence
 
 import numpy as np
@@ -206,7 +207,7 @@ def _assemble_nbody_components(
         first_key = next(iter(component_results.keys()))
         property_shape = find_shape(component_results[first_key])
 
-        # Final dictionaries
+        # Accumulation dictionaries
         # * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP
         #   & CP the summed total energies (or other property) of each nb-body. That is:
         #   * NOCP: {1: 1b@1b,    2: 2b@2b,      ..., max_nbody: max_nbody-b@max_nbody-b} and
@@ -218,6 +219,17 @@ def _assemble_nbody_components(
         nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
         vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
 
+        # * {bsse_type}_body_dict is usually filled with total energies (or other property).
+        #   Multiple model chemistry levels may be involved.
+        #   Generally, all consecutive keys between 1 and max_nbody will be present in the body_dict,
+        #   but if supersystem_ie_only=T, only 1b and nfr-b are present, or if "supersystem" in levels, ???
+        #   * TOT: {1: 1b@1b, 2: 2b tot prop with bsse_type treatment, ..., max_nbody: max_nbody-b tot prop with bsse_type treatment}
+        #   If 1b@1b (monomers in monomer basis) aren't available, which can happen when return_total_data=F
+        #   and 1b@1b aren't otherwise needed, body_dict contains interaction energies (or other property).
+        #   * IE: {1: shaped_zero, 2: 2b interaction prop using bsse_type, ..., max_nbody: max_nbody-b interaction prop using bsse_type}
+        #   For both TOT and IE cases, body_dict values are cummulative, not additive. For TOT, total,
+        #   interaction, and contribution data in ManyBodyResultProperties can be computed in
+        #   collect_vars. For IE, interaction and contribution data can be computed.
         cp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
         nocp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
         vmfc_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
@@ -321,6 +333,7 @@ def _analyze(
 
         # results per model chemistry
         mc_results = {}
+        species_results = {}
 
         # sort by nbody level, ignore supersystem
         sorted_nbodies = [(k, v) for k, v in self.nbodies_per_mc_level.items() if v != ["supersystem"]]
@@ -455,23 +468,28 @@ def analyze(
             component_results_inv["energy"] = {'["dummy", [1000], [1000]]': 0.0}
 
         # Actually analyze
+        is_embedded = bool(self.embedding_charges)
+        component_results_qq = defaultdict(dict)
         all_results = {}
+        nbody_dict = {}
+#        all_results["energy_body_dict"] = {"cp": {1: 0.0}}
 
         for property_label, property_results in component_results_inv.items():
+            for k, v in property_results.items():
+                # TODO before or after expansion?
+                component_results_qq[k][property_label] = v
             # Expand gradient and hessian
             if property_label == "gradient":
                 property_results = {k: self.expand_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}
             if property_label == "hessian":
                 property_results = {k: self.expand_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}
 
             r = self._analyze(property_label, property_results)
+            # TODO before or after?
+            #for k, v in property_results.items():
+            #    component_results_qq[k][property_label] = v
             all_results.update(r)
 
-        # Analyze the total results
-        nbody_dict = {}
-
-        is_embedded = bool(self.embedding_charges)
-
         for bt in self.bsse_type:
             print_nbody_energy(
                 all_results["energy_body_dict"][bt],
@@ -480,12 +498,22 @@ def analyze(
                 is_embedded,
             )
 
-            if not self.has_supersystem:  # skipped levels?
-                nbody_dict.update(
-                    collect_vars(bt.upper(), all_results["energy_body_dict"][bt], self.max_nbody, is_embedded, self.supersystem_ie_only)
-                )
+        for property_label in available_properties:
+            for bt in self.bsse_type:
+                if not self.has_supersystem:  # skipped levels?
+                    nbody_dict.update(
+                        collect_vars(
+                            bt.upper(),
+                            property_label.upper(),
+                            all_results[f"{property_label}_body_dict"][bt],
+                            self.max_nbody,
+                            is_embedded,
+                            self.supersystem_ie_only,
+                        )
+                    )
 
         all_results["results"] = nbody_dict
+        all_results["component_data"] = component_results_qq
 
         # Make dictionary with "1cp", "2cp", etc
         ebd = all_results["energy_body_dict"]

qcmanybody/utils.py

@@ -99,9 +99,9 @@ def labeler(mc_level_lbl: str, frag: Tuple[int, ...], bas: Tuple[int, ...]) -> s
         Key identifying the model chemistry. May be `"(auto)"`. Often the
         ManyBodyInput.specification.specification keys.
     frag
-        1-indexed list of fragments active in the supersystem.
+        List of 1-indexed fragments active in the supersystem.
     bas
-        1-indexed list of fragments with active basis sets in the supersystem.
+        List of 1-indexed fragments with active basis sets in the supersystem.
         All those in *frag* plus any ghost.
 
     Returns
@@ -164,7 +164,35 @@ def print_nbody_energy(
     print(info)
 
 
-def collect_vars(bsse, body_dict, max_nbody: int, embedding: bool = False, supersystem_ie_only: bool = False):
+def collect_vars(
+        bsse: str,
+        prop: str,
+        body_dict: Mapping[int, Union[float, np.ndarray]],
+        max_nbody: int,
+        embedding: bool = False,
+        supersystem_ie_only: bool = False,
+    ) -> Dict:
+    """From *body_dict*, construct QCVariables.
+
+    Parameters
+    ----------
+    bsse
+        Uppercase label for a single many-body treatment, generally a value of BsseEnum.
+    prop
+        Uppercase label for a single property, generally a value of DriverEnum.
+    body_dict
+        Dictionary of minimal per-body info already specialized for property *prop* and treatment *bsse*. May contain either total data or interaction data, both cummulative not additive, from 1-body to max_nbody-body (see also *supersystem_ie_only*). Interaction data signaled by zero float or array for 1-body. May contain multiple model chemistry levels.
+    max_nbody
+        _description_
+    embedding, optional
+        Is charge embedding enabled, by default False?
+    supersystem_ie_only, optional
+        Is data available in *body_dict* only for 1-body (possibly zero) and nfr-body levels? By default False: data is available for consecutive levels, up to max_nbody-body.
+
+    Returns
+    -------
+        _description_. Empty return if *embedding* enabled.
+    """
     previous_e = body_dict[1]
     property_shape = find_shape(previous_e)
     tot_e = bool(np.count_nonzero(previous_e))