diff --git a/api/index.html b/api/index.html index cdcabe5..d77545c 100644 --- a/api/index.html +++ b/api/index.html @@ -1543,7 +1543,7 @@
????
+Target the supersystem total/interaction energy (IE) data over the many-body expansion (MBE) " +analysis, thereby omitting intermediate-body calculations.
supersystem_max_nbody
Optional[int]
+ Maximum n-body to use for a supersystem calculation. Must be specified if "supersystem" is in nbodies
None
+ def build_nbody_compute_list(
+175
+176
+177
+178
def build_nbody_compute_list(
bsse_type: Iterable[BsseEnum],
nfragments: int,
nbodies: Iterable[Union[int, Literal["supersystem"]]],
@@ -2355,7 +2373,10 @@ Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested,
as opposed to interaction data (False).
supersystem_ie_only
- ????
+ Target the supersystem total/interaction energy (IE) data over the many-body expansion (MBE) "
+ analysis, thereby omitting intermediate-body calculations.
+ supersystem_max_nbody
+ Maximum n-body to use for a supersystem calculation. Must be specified if "supersystem" is in `nbodies`
Returns
-------
@@ -2388,8 +2409,8 @@ raise ValueError("supersystem_max_nbody must be provided if 'supersystem' contains nbodies")
include_supersystem = True
- nbodies = list(nbodies)
- nbodies.remove("supersystem")
+
+ nbodies: List[int] = [x for x in nbodies if x != "supersystem"]
# What levels do we need?
fragment_range = range(1, nfragments + 1)
diff --git a/search/search_index.json b/search/search_index.json
index 0a6657a..7e1a769 100644
--- a/search/search_index.json
+++ b/search/search_index.json
@@ -1 +1 @@
-{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Welcome to MkDocs","text":"
For full documentation visit mkdocs.org.
"},{"location":"#commands","title":"Commands","text":" mkdocs new [dir-name]
- Create a new project. mkdocs serve
- Start the live-reloading docs server. mkdocs build
- Build the documentation site. mkdocs -h
- Print help message and exit.
"},{"location":"#table-of-contents","title":"Table Of Contents","text":" - Tutorials
- How-To Guides
- Reference
- Explanation
"},{"location":"api/","title":"API Documentation","text":"Source code in qcmanybody/manybody.py
class ManyBodyCalculator:\n def __init__(\n self,\n molecule: Molecule,\n bsse_type: Sequence[BsseEnum],\n levels: Mapping[Union[int, Literal[\"supersystem\"]], str],\n return_total_data: bool,\n supersystem_ie_only: bool,\n ):\n # TODO\n self.embedding_charges = {}\n\n self.molecule = molecule\n self.bsse_type = [BsseEnum(x) for x in bsse_type]\n self.return_total_data = return_total_data\n self.supersystem_ie_only = supersystem_ie_only\n self.nfragments = len(molecule.fragments)\n\n self.levels = levels\n\n # Levels without supersystem\n self.levels_no_ss = {int(k): v for k, v in levels.items() if k != \"supersystem\"}\n\n # Just a set of all the modelchems\n self.mc_levels = set(self.levels.values())\n\n self.max_nbody = max(self.levels_no_ss.keys())\n\n if len(self.bsse_type) == 0:\n raise ValueError(\"No BSSE correction specified\")\n\n if BsseEnum.vmfc in self.bsse_type and len(set(self.levels.values())) == 1:\n # For single-modelchem VMFC, NOCP & sometimes CP are produced for free\n if BsseEnum.nocp not in self.bsse_type:\n self.bsse_type.append(BsseEnum.nocp)\n if BsseEnum.cp not in self.bsse_type and self.max_nbody == self.nfragments:\n self.bsse_type.append(BsseEnum.cp)\n\n self.return_bsse_type = self.bsse_type[0]\n\n ###############################\n # Build nbodies_per_mc_level\n # TODO - use Lori's code\n # TODO - dict to list of lists to handle non-contiguous levels\n # TODO multilevel and supersystem_ie_only=T not allowed together\n # TODO supersystem in levels is not to be trusted -- nfrag only and skips levels\n max_level = max(self.levels_no_ss.keys())\n\n if set(range(1, max_level + 1)) != set(self.levels_no_ss.keys()):\n raise ValueError(f\"Levels must be contiguous from 1 to {max_level}\")\n\n self.nbodies_per_mc_level: Dict[str, list] = {mc_level: [] for mc_level in self.mc_levels}\n for k, v in self.levels_no_ss.items():\n self.nbodies_per_mc_level[v].append(k)\n\n self.nbodies_per_mc_level = {k: sorted(v) for k, v in self.nbodies_per_mc_level.items()}\n\n # Supersystem is always at the end\n if \"supersystem\" in levels:\n ss_mc = levels[\"supersystem\"]\n self.nbodies_per_mc_level[ss_mc].append(\"supersystem\")\n\n # To be built on the fly\n self.mc_compute_dict = None\n\n if not np.array_equal(np.concatenate(self.molecule.fragments), np.arange(len(self.molecule.symbols))):\n raise ValueError(\"\"\"QCManyBody: non-contiguous fragments could be implemented but aren't at present\"\"\")\n\n # Build size and slices dictionaries. Assumes fragments are contiguous\n self.fragment_size_dict = {}\n self.fragment_slice_dict = {}\n iat = 0\n for ifr in range(1, self.nfragments + 1):\n nat = len(self.molecule.fragments[ifr - 1])\n self.fragment_size_dict[ifr] = nat\n self.fragment_slice_dict[ifr] = slice(iat, iat + nat)\n iat += nat\n\n @property\n def has_supersystem(self) -> bool:\n return \"supersystem\" in self.levels\n\n @property\n def compute_map(self) -> Dict[str, Dict[str, Dict[int, Set[FragBasIndex]]]]:\n if self.mc_compute_dict is not None:\n return self.mc_compute_dict\n\n # Build the compute lists\n self.mc_compute_dict = {}\n\n for mc in self.mc_levels:\n nbodies = self.nbodies_per_mc_level[mc]\n self.mc_compute_dict[mc] = build_nbody_compute_list(\n self.bsse_type,\n self.nfragments,\n nbodies,\n self.return_total_data,\n self.supersystem_ie_only,\n self.max_nbody,\n )\n\n return self.mc_compute_dict\n\n def resize_gradient(self, grad: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_gradient(grad, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def resize_hessian(self, hess: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_hessian(hess, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n\n def _assemble_nbody_components(\n self,\n property_label: str,\n component_results: Dict[str, Union[float, np.ndarray]],\n ) -> Dict[str, Any]:\n \"\"\"Assembles N-body components for a single derivative level and a single model chemistry level\n into interaction quantities according to requested BSSE treatment(s).\n \"\"\"\n\n # which level are we assembling?\n delabeled = [delabeler(k) for k in component_results.keys()]\n mc_level_labels = {x[0] for x in delabeled}\n\n if len(mc_level_labels) != 1:\n raise RuntimeError(f\"Multiple model chemistries passed into _assemble_nbody_components: {mc_level_labels}\")\n\n mc_level = mc_level_labels.pop()\n if mc_level not in self.mc_levels:\n raise RuntimeError(f\"Model chemistry {mc_level} not found in {self.mc_levels}\")\n\n # get the range of nbodies and the required calculations for this level\n bsse_type = self.bsse_type\n return_bsse_type = self.return_bsse_type\n nbodies = self.nbodies_per_mc_level[mc_level]\n if \"supersystem\" in nbodies:\n nbodies = list(range(1, self.max_nbody + 1))\n bsse_type = [BsseEnum.nocp]\n return_bsse_type = BsseEnum.nocp\n\n max_nbody = max(nbodies)\n compute_dict = self.compute_map[mc_level]\n\n if not all_same_shape(component_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(component_results.keys()))\n property_shape = find_shape(component_results[first_key])\n\n # Accumulation dictionaries\n # * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP\n # & CP the summed total energies (or other property) of each nb-body. That is:\n # * NOCP: {1: 1b@1b, 2: 2b@2b, ..., max_nbody: max_nbody-b@max_nbody-b} and\n # * CP: {1: 1b@nfr-b, 2: 2b@nfr-b, ..., max_nbody: max_nbody-b@nfr-b}.\n # VMFC bookkeeping is different. For key 1 it contains the summed 1b total energies.\n # But for higher keys, it contains each nb-body (non-additive) contribution to the energy.\n # * VMFC: {1: 1b@1b, 2: 2b contrib, ..., max_nbody: max_nbody-b contrib}\n cp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # * {bsse_type}_body_dict is usually filled with total energies (or other property).\n # Multiple model chemistry levels may be involved.\n # Generally, all consecutive keys between 1 and max_nbody will be present in the body_dict,\n # but if supersystem_ie_only=T, only 1b and nfr-b are present, or if \"supersystem\" in levels, ???\n # * TOT: {1: 1b@1b, 2: 2b tot prop with bsse_type treatment, ..., max_nbody: max_nbody-b tot prop with bsse_type treatment}\n # If 1b@1b (monomers in monomer basis) aren't available, which can happen when return_total_data=F\n # and 1b@1b aren't otherwise needed, body_dict contains interaction energies (or other property).\n # * IE: {1: shaped_zero, 2: 2b interaction prop using bsse_type, ..., max_nbody: max_nbody-b interaction prop using bsse_type}\n # For both TOT and IE cases, body_dict values are cummulative, not additive. For TOT, total,\n # interaction, and contribution data in ManyBodyResultProperties can be computed in\n # collect_vars. For IE, interaction and contribution data can be computed.\n cp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # Sum up all of the levels\n # * compute_dict[bt][nb] holds all the computations needed to compute nb\n # *not* all the nb-level computations, so build the latter\n cp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n nocp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n\n for nb in nbodies:\n for v in compute_dict[\"cp\"][nb]:\n if len(v[1]) != 1:\n cp_compute_list[len(v[0])].add(v)\n for w in compute_dict[\"nocp\"][nb]:\n nocp_compute_list[len(w[0])].add(w)\n\n for nb in range(1, nbodies[-1] + 1):\n cp_by_level[nb] = sum_cluster_data(component_results, cp_compute_list[nb], mc_level)\n nocp_by_level[nb] = sum_cluster_data(component_results, nocp_compute_list[nb], mc_level)\n if nb in compute_dict[\"vmfc_levels\"]:\n vmfc_by_level[nb] = sum_cluster_data(\n component_results, compute_dict[\"vmfc_levels\"][nb], mc_level, vmfc=True, nb=nb\n )\n\n # Extract data for monomers in monomer basis for CP total data\n if 1 in nbodies:\n monomers_in_monomer_basis = [v for v in compute_dict[\"nocp\"][1] if len(v[1]) == 1]\n monomer_sum = sum_cluster_data(component_results, set(monomers_in_monomer_basis), mc_level)\n else:\n monomer_sum = shaped_zero(property_shape)\n\n # Compute cp\n if BsseEnum.cp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n cp_body_dict[nb] = cp_by_level[nb] - bsse\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n cp_body_dict[nb] += take_nk * sign * cp_by_level[k]\n\n if nb == 1:\n bsse = cp_body_dict[nb] - monomer_sum\n cp_body_dict[nb] = copy_value(monomer_sum)\n else:\n cp_body_dict[nb] -= bsse\n\n # Compute nocp\n if BsseEnum.nocp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n nocp_body_dict[nb] = nocp_by_level[nb]\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n nocp_body_dict[nb] += take_nk * sign * nocp_by_level[k]\n\n # Compute vmfc\n if BsseEnum.vmfc in bsse_type:\n for nb in nbodies:\n for k in range(1, nb + 1):\n vmfc_body_dict[nb] += vmfc_by_level[k]\n\n # Collect specific and generalized returns\n results = {\n f\"cp_{property_label}_body_dict\": cp_body_dict,\n f\"nocp_{property_label}_body_dict\": nocp_body_dict,\n f\"vmfc_{property_label}_body_dict\": vmfc_body_dict,\n }\n\n # Overall return body dict & value for this property\n results[f\"{property_label}_body_dict\"] = results[f\"{return_bsse_type.value}_{property_label}_body_dict\"]\n results[f\"ret_{property_label}\"] = copy_value(results[f\"{property_label}_body_dict\"][max_nbody])\n\n if not self.return_total_data:\n results[f\"ret_{property_label}\"] -= results[f\"{property_label}_body_dict\"][1]\n\n return results\n\n def _analyze(\n self,\n property_label: str,\n property_results: Dict[str, Union[float, np.ndarray]], # Label to results\n ):\n # Initialize with zeros\n if not all_same_shape(property_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(property_results.keys()))\n property_shape = find_shape(property_results[first_key])\n\n property_result = shaped_zero(property_shape)\n property_body_dict = {bt.value: {} for bt in self.bsse_type}\n property_body_contribution = {bt.value: {} for bt in self.bsse_type}\n\n # results per model chemistry\n mc_results = {}\n species_results = {}\n\n # sort by nbody level, ignore supersystem\n sorted_nbodies = [(k, v) for k, v in self.nbodies_per_mc_level.items() if v != [\"supersystem\"]]\n sorted_nbodies = sorted(sorted_nbodies, reverse=True, key=lambda x: x[1])\n for mc_label, nbody_list in sorted_nbodies:\n # filter to only one model chemistry\n filtered_results = {k: v for k, v in property_results.items() if delabeler(k)[0] == mc_label}\n\n if not filtered_results:\n if nbody_list == [1]:\n # Note A.2: Note A.1 holds, but for the special case of CP-only\n # and rtd=False and multilevel with a separate level for\n # 1-body, the skipped tasks run afoul of sanity checks, so\n # we'll add a dummy result.\n filtered_results = {labeler(mc_label, [1000], [1000]): shaped_zero(property_shape)}\n else:\n raise RuntimeError(f\"No data found for model chemistry {mc_label}\")\n\n nb_component_results = self._assemble_nbody_components(property_label, filtered_results)\n mc_results[mc_label] = nb_component_results\n\n for n in nbody_list[::-1]:\n property_bsse_dict = {bt.value: shaped_zero(property_shape) for bt in self.bsse_type}\n\n for m in range(n - 1, n + 1):\n if m == 0:\n continue\n\n # Subtract the (n-1)-body contribution from the n-body contribution to get the n-body effect\n sign = (-1) ** (1 - m // n)\n for bt in self.bsse_type:\n property_bsse_dict[bt.value] += (\n sign * mc_results[mc_label][f\"{bt.value}_{property_label}_body_dict\"][m]\n )\n\n property_result += property_bsse_dict[self.return_bsse_type]\n for bt in self.bsse_type:\n property_body_contribution[bt.value][n] = property_bsse_dict[bt.value]\n\n if self.has_supersystem:\n # Get the MC label for supersystem tasks\n supersystem_mc_level = self.levels[\"supersystem\"]\n\n # Super system recovers higher order effects at a lower level\n frag_range = tuple(range(1, self.nfragments + 1))\n\n ss_cresults = {k: v for k, v in property_results.items() if delabeler(k)[0] == supersystem_mc_level}\n ss_component_results = self._assemble_nbody_components(property_label, ss_cresults)\n mc_results[supersystem_mc_level] = ss_component_results\n\n # Compute components at supersystem level of theory\n ss_label = labeler(supersystem_mc_level, frag_range, frag_range)\n supersystem_result = property_results[ss_label]\n property_result += supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n\n for bt in self.bsse_type:\n property_body_contribution[bt][self.nfragments] = (\n supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n )\n\n for bt in self.bsse_type:\n bstr = bt.value\n for n in property_body_contribution[bstr]:\n property_body_dict[bstr][n] = sum(\n [\n property_body_contribution[bstr][i]\n for i in range(1, n + 1)\n if i in property_body_contribution[bstr]\n ]\n )\n\n if not self.return_total_data:\n # Remove monomer contribution for interaction data\n property_result -= property_body_dict[self.return_bsse_type][1]\n\n nbody_results = {\n f\"ret_{property_label}\": property_result,\n f\"{property_label}_body_dict\": property_body_dict,\n \"mc_results\": mc_results,\n }\n return nbody_results\n\n def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties = set([\"energy\"])\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"api/#qcmanybody.ManyBodyCalculator.analyze","title":"analyze(component_results)
","text":"Parameters:
Name Type Description Default component_results
Dict[str, Dict[str, Union[float, ndarray]]]
Nested dictionary with results from all individual molecular system calculations, including all subsystem/basis combinations, all model chemistries, and all properties (e.g., e/g/h).
For example, the below is the format for a nocp gradient run on a helium dimer with 1-body at CCSD and 2-body at MP2. The outer string key can be generated with the qcmanybody.utils.labeler
function. The inner string key is any property; QCManyBody presently knows how to process energy/gradient/Hessian.
{'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])}, '[\"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]])}, }
required Return Source code in qcmanybody/manybody.py
def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties = set([\"energy\"])\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"api/#qcmanybody.ManyBodyCalculator.iterate_molecules","title":"iterate_molecules()
","text":"Iterate over all the molecules needed for the computation.
Yields model chemistry, label, and molecule.
Source code in qcmanybody/manybody.py
def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n
"},{"location":"explanation/","title":"Explanation","text":"This part of the project documentation focuses on an understanding-oriented approach. You'll get a chance to read about the background of the project, as well as reasoning about how it was implemented.
Note: Expand this section by considering the following points:
- Give context and background on your library
- Explain why you created it
- Provide multiple examples and approaches of how to work with it
- Help the reader make connections
- Avoid writing instructions or technical descriptions here
"},{"location":"reference/","title":"Reference","text":"This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the QCManyBody
project code.
Source code in qcmanybody/manybody.py
class ManyBodyCalculator:\n def __init__(\n self,\n molecule: Molecule,\n bsse_type: Sequence[BsseEnum],\n levels: Mapping[Union[int, Literal[\"supersystem\"]], str],\n return_total_data: bool,\n supersystem_ie_only: bool,\n ):\n # TODO\n self.embedding_charges = {}\n\n self.molecule = molecule\n self.bsse_type = [BsseEnum(x) for x in bsse_type]\n self.return_total_data = return_total_data\n self.supersystem_ie_only = supersystem_ie_only\n self.nfragments = len(molecule.fragments)\n\n self.levels = levels\n\n # Levels without supersystem\n self.levels_no_ss = {int(k): v for k, v in levels.items() if k != \"supersystem\"}\n\n # Just a set of all the modelchems\n self.mc_levels = set(self.levels.values())\n\n self.max_nbody = max(self.levels_no_ss.keys())\n\n if len(self.bsse_type) == 0:\n raise ValueError(\"No BSSE correction specified\")\n\n if BsseEnum.vmfc in self.bsse_type and len(set(self.levels.values())) == 1:\n # For single-modelchem VMFC, NOCP & sometimes CP are produced for free\n if BsseEnum.nocp not in self.bsse_type:\n self.bsse_type.append(BsseEnum.nocp)\n if BsseEnum.cp not in self.bsse_type and self.max_nbody == self.nfragments:\n self.bsse_type.append(BsseEnum.cp)\n\n self.return_bsse_type = self.bsse_type[0]\n\n ###############################\n # Build nbodies_per_mc_level\n # TODO - use Lori's code\n # TODO - dict to list of lists to handle non-contiguous levels\n # TODO multilevel and supersystem_ie_only=T not allowed together\n # TODO supersystem in levels is not to be trusted -- nfrag only and skips levels\n max_level = max(self.levels_no_ss.keys())\n\n if set(range(1, max_level + 1)) != set(self.levels_no_ss.keys()):\n raise ValueError(f\"Levels must be contiguous from 1 to {max_level}\")\n\n self.nbodies_per_mc_level: Dict[str, list] = {mc_level: [] for mc_level in self.mc_levels}\n for k, v in self.levels_no_ss.items():\n self.nbodies_per_mc_level[v].append(k)\n\n self.nbodies_per_mc_level = {k: sorted(v) for k, v in self.nbodies_per_mc_level.items()}\n\n # Supersystem is always at the end\n if \"supersystem\" in levels:\n ss_mc = levels[\"supersystem\"]\n self.nbodies_per_mc_level[ss_mc].append(\"supersystem\")\n\n # To be built on the fly\n self.mc_compute_dict = None\n\n if not np.array_equal(np.concatenate(self.molecule.fragments), np.arange(len(self.molecule.symbols))):\n raise ValueError(\"\"\"QCManyBody: non-contiguous fragments could be implemented but aren't at present\"\"\")\n\n # Build size and slices dictionaries. Assumes fragments are contiguous\n self.fragment_size_dict = {}\n self.fragment_slice_dict = {}\n iat = 0\n for ifr in range(1, self.nfragments + 1):\n nat = len(self.molecule.fragments[ifr - 1])\n self.fragment_size_dict[ifr] = nat\n self.fragment_slice_dict[ifr] = slice(iat, iat + nat)\n iat += nat\n\n @property\n def has_supersystem(self) -> bool:\n return \"supersystem\" in self.levels\n\n @property\n def compute_map(self) -> Dict[str, Dict[str, Dict[int, Set[FragBasIndex]]]]:\n if self.mc_compute_dict is not None:\n return self.mc_compute_dict\n\n # Build the compute lists\n self.mc_compute_dict = {}\n\n for mc in self.mc_levels:\n nbodies = self.nbodies_per_mc_level[mc]\n self.mc_compute_dict[mc] = build_nbody_compute_list(\n self.bsse_type,\n self.nfragments,\n nbodies,\n self.return_total_data,\n self.supersystem_ie_only,\n self.max_nbody,\n )\n\n return self.mc_compute_dict\n\n def resize_gradient(self, grad: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_gradient(grad, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def resize_hessian(self, hess: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_hessian(hess, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n\n def _assemble_nbody_components(\n self,\n property_label: str,\n component_results: Dict[str, Union[float, np.ndarray]],\n ) -> Dict[str, Any]:\n \"\"\"Assembles N-body components for a single derivative level and a single model chemistry level\n into interaction quantities according to requested BSSE treatment(s).\n \"\"\"\n\n # which level are we assembling?\n delabeled = [delabeler(k) for k in component_results.keys()]\n mc_level_labels = {x[0] for x in delabeled}\n\n if len(mc_level_labels) != 1:\n raise RuntimeError(f\"Multiple model chemistries passed into _assemble_nbody_components: {mc_level_labels}\")\n\n mc_level = mc_level_labels.pop()\n if mc_level not in self.mc_levels:\n raise RuntimeError(f\"Model chemistry {mc_level} not found in {self.mc_levels}\")\n\n # get the range of nbodies and the required calculations for this level\n bsse_type = self.bsse_type\n return_bsse_type = self.return_bsse_type\n nbodies = self.nbodies_per_mc_level[mc_level]\n if \"supersystem\" in nbodies:\n nbodies = list(range(1, self.max_nbody + 1))\n bsse_type = [BsseEnum.nocp]\n return_bsse_type = BsseEnum.nocp\n\n max_nbody = max(nbodies)\n compute_dict = self.compute_map[mc_level]\n\n if not all_same_shape(component_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(component_results.keys()))\n property_shape = find_shape(component_results[first_key])\n\n # Accumulation dictionaries\n # * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP\n # & CP the summed total energies (or other property) of each nb-body. That is:\n # * NOCP: {1: 1b@1b, 2: 2b@2b, ..., max_nbody: max_nbody-b@max_nbody-b} and\n # * CP: {1: 1b@nfr-b, 2: 2b@nfr-b, ..., max_nbody: max_nbody-b@nfr-b}.\n # VMFC bookkeeping is different. For key 1 it contains the summed 1b total energies.\n # But for higher keys, it contains each nb-body (non-additive) contribution to the energy.\n # * VMFC: {1: 1b@1b, 2: 2b contrib, ..., max_nbody: max_nbody-b contrib}\n cp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # * {bsse_type}_body_dict is usually filled with total energies (or other property).\n # Multiple model chemistry levels may be involved.\n # Generally, all consecutive keys between 1 and max_nbody will be present in the body_dict,\n # but if supersystem_ie_only=T, only 1b and nfr-b are present, or if \"supersystem\" in levels, ???\n # * TOT: {1: 1b@1b, 2: 2b tot prop with bsse_type treatment, ..., max_nbody: max_nbody-b tot prop with bsse_type treatment}\n # If 1b@1b (monomers in monomer basis) aren't available, which can happen when return_total_data=F\n # and 1b@1b aren't otherwise needed, body_dict contains interaction energies (or other property).\n # * IE: {1: shaped_zero, 2: 2b interaction prop using bsse_type, ..., max_nbody: max_nbody-b interaction prop using bsse_type}\n # For both TOT and IE cases, body_dict values are cummulative, not additive. For TOT, total,\n # interaction, and contribution data in ManyBodyResultProperties can be computed in\n # collect_vars. For IE, interaction and contribution data can be computed.\n cp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # Sum up all of the levels\n # * compute_dict[bt][nb] holds all the computations needed to compute nb\n # *not* all the nb-level computations, so build the latter\n cp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n nocp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n\n for nb in nbodies:\n for v in compute_dict[\"cp\"][nb]:\n if len(v[1]) != 1:\n cp_compute_list[len(v[0])].add(v)\n for w in compute_dict[\"nocp\"][nb]:\n nocp_compute_list[len(w[0])].add(w)\n\n for nb in range(1, nbodies[-1] + 1):\n cp_by_level[nb] = sum_cluster_data(component_results, cp_compute_list[nb], mc_level)\n nocp_by_level[nb] = sum_cluster_data(component_results, nocp_compute_list[nb], mc_level)\n if nb in compute_dict[\"vmfc_levels\"]:\n vmfc_by_level[nb] = sum_cluster_data(\n component_results, compute_dict[\"vmfc_levels\"][nb], mc_level, vmfc=True, nb=nb\n )\n\n # Extract data for monomers in monomer basis for CP total data\n if 1 in nbodies:\n monomers_in_monomer_basis = [v for v in compute_dict[\"nocp\"][1] if len(v[1]) == 1]\n monomer_sum = sum_cluster_data(component_results, set(monomers_in_monomer_basis), mc_level)\n else:\n monomer_sum = shaped_zero(property_shape)\n\n # Compute cp\n if BsseEnum.cp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n cp_body_dict[nb] = cp_by_level[nb] - bsse\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n cp_body_dict[nb] += take_nk * sign * cp_by_level[k]\n\n if nb == 1:\n bsse = cp_body_dict[nb] - monomer_sum\n cp_body_dict[nb] = copy_value(monomer_sum)\n else:\n cp_body_dict[nb] -= bsse\n\n # Compute nocp\n if BsseEnum.nocp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n nocp_body_dict[nb] = nocp_by_level[nb]\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n nocp_body_dict[nb] += take_nk * sign * nocp_by_level[k]\n\n # Compute vmfc\n if BsseEnum.vmfc in bsse_type:\n for nb in nbodies:\n for k in range(1, nb + 1):\n vmfc_body_dict[nb] += vmfc_by_level[k]\n\n # Collect specific and generalized returns\n results = {\n f\"cp_{property_label}_body_dict\": cp_body_dict,\n f\"nocp_{property_label}_body_dict\": nocp_body_dict,\n f\"vmfc_{property_label}_body_dict\": vmfc_body_dict,\n }\n\n # Overall return body dict & value for this property\n results[f\"{property_label}_body_dict\"] = results[f\"{return_bsse_type.value}_{property_label}_body_dict\"]\n results[f\"ret_{property_label}\"] = copy_value(results[f\"{property_label}_body_dict\"][max_nbody])\n\n if not self.return_total_data:\n results[f\"ret_{property_label}\"] -= results[f\"{property_label}_body_dict\"][1]\n\n return results\n\n def _analyze(\n self,\n property_label: str,\n property_results: Dict[str, Union[float, np.ndarray]], # Label to results\n ):\n # Initialize with zeros\n if not all_same_shape(property_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(property_results.keys()))\n property_shape = find_shape(property_results[first_key])\n\n property_result = shaped_zero(property_shape)\n property_body_dict = {bt.value: {} for bt in self.bsse_type}\n property_body_contribution = {bt.value: {} for bt in self.bsse_type}\n\n # results per model chemistry\n mc_results = {}\n species_results = {}\n\n # sort by nbody level, ignore supersystem\n sorted_nbodies = [(k, v) for k, v in self.nbodies_per_mc_level.items() if v != [\"supersystem\"]]\n sorted_nbodies = sorted(sorted_nbodies, reverse=True, key=lambda x: x[1])\n for mc_label, nbody_list in sorted_nbodies:\n # filter to only one model chemistry\n filtered_results = {k: v for k, v in property_results.items() if delabeler(k)[0] == mc_label}\n\n if not filtered_results:\n if nbody_list == [1]:\n # Note A.2: Note A.1 holds, but for the special case of CP-only\n # and rtd=False and multilevel with a separate level for\n # 1-body, the skipped tasks run afoul of sanity checks, so\n # we'll add a dummy result.\n filtered_results = {labeler(mc_label, [1000], [1000]): shaped_zero(property_shape)}\n else:\n raise RuntimeError(f\"No data found for model chemistry {mc_label}\")\n\n nb_component_results = self._assemble_nbody_components(property_label, filtered_results)\n mc_results[mc_label] = nb_component_results\n\n for n in nbody_list[::-1]:\n property_bsse_dict = {bt.value: shaped_zero(property_shape) for bt in self.bsse_type}\n\n for m in range(n - 1, n + 1):\n if m == 0:\n continue\n\n # Subtract the (n-1)-body contribution from the n-body contribution to get the n-body effect\n sign = (-1) ** (1 - m // n)\n for bt in self.bsse_type:\n property_bsse_dict[bt.value] += (\n sign * mc_results[mc_label][f\"{bt.value}_{property_label}_body_dict\"][m]\n )\n\n property_result += property_bsse_dict[self.return_bsse_type]\n for bt in self.bsse_type:\n property_body_contribution[bt.value][n] = property_bsse_dict[bt.value]\n\n if self.has_supersystem:\n # Get the MC label for supersystem tasks\n supersystem_mc_level = self.levels[\"supersystem\"]\n\n # Super system recovers higher order effects at a lower level\n frag_range = tuple(range(1, self.nfragments + 1))\n\n ss_cresults = {k: v for k, v in property_results.items() if delabeler(k)[0] == supersystem_mc_level}\n ss_component_results = self._assemble_nbody_components(property_label, ss_cresults)\n mc_results[supersystem_mc_level] = ss_component_results\n\n # Compute components at supersystem level of theory\n ss_label = labeler(supersystem_mc_level, frag_range, frag_range)\n supersystem_result = property_results[ss_label]\n property_result += supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n\n for bt in self.bsse_type:\n property_body_contribution[bt][self.nfragments] = (\n supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n )\n\n for bt in self.bsse_type:\n bstr = bt.value\n for n in property_body_contribution[bstr]:\n property_body_dict[bstr][n] = sum(\n [\n property_body_contribution[bstr][i]\n for i in range(1, n + 1)\n if i in property_body_contribution[bstr]\n ]\n )\n\n if not self.return_total_data:\n # Remove monomer contribution for interaction data\n property_result -= property_body_dict[self.return_bsse_type][1]\n\n nbody_results = {\n f\"ret_{property_label}\": property_result,\n f\"{property_label}_body_dict\": property_body_dict,\n \"mc_results\": mc_results,\n }\n return nbody_results\n\n def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties = set([\"energy\"])\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
Generates lists of N-Body computations needed for requested BSSE treatments.
Parameters:
Name Type Description Default bsse_type
Iterable[BsseEnum]
Requested BSSE treatments.
required nfragments
int
Number of distinct fragments comprising the full molecular supersystem.
required nbodies
Iterable[Union[int, Literal['supersystem']]]
List of n-body levels (e.g., [2]
or [1, 2]
or [\"supersystem\"]
) for which to generate tasks. Note the natural 1-indexing, so [1]
covers one-body contributions.
required return_total_data
bool
Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested, as opposed to interaction data (False).
required supersystem_ie_only
bool
????
required Returns:
Type Description compute_dict
Dictionary containing subdicts enumerating compute lists for each possible BSSE treatment. Subdict keys are n-body levels and values are sets of all the mc_(frag, bas)
indices needed to compute that n-body level. A given index can appear multiple times within a subdict and among subdicts.
compute_dict[\"cp\"] = {\n 1: set(),\n 2: {((1,), (1, 2)),\n ((2,), (1, 2)),\n ((1, 2), (1, 2))}\n}\n
Subdicts below are always returned. Any may be empty if not requested through bsse_type.
'all'
|w---w| full list of computations required 'cp'
|w---w| list of computations required for CP procedure 'nocp'
|w---w| list of computations required for non-CP procedure 'vmfc_compute'
|w---w| list of computations required for VMFC procedure 'vmfc_levels'
|w---w| list of levels required for VMFC procedure
Source code in qcmanybody/builder.py
def build_nbody_compute_list(\n bsse_type: Iterable[BsseEnum],\n nfragments: int,\n nbodies: Iterable[Union[int, Literal[\"supersystem\"]]],\n return_total_data: bool,\n supersystem_ie_only: bool,\n supersystem_max_nbody: Optional[int] = None,\n) -> Dict[str, Dict[int, Set[FragBasIndex]]]:\n \"\"\"Generates lists of N-Body computations needed for requested BSSE treatments.\n\n Parameters\n ----------\n bsse_type\n Requested BSSE treatments.\n nfragments\n Number of distinct fragments comprising the full molecular supersystem.\n nbodies\n List of n-body levels (e.g., `[2]` or `[1, 2]` or `[\"supersystem\"]`) for which to generate tasks.\n Note the natural 1-indexing, so `[1]` covers one-body contributions.\n return_total_data\n Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested,\n as opposed to interaction data (False).\n supersystem_ie_only\n ????\n\n Returns\n -------\n compute_dict\n Dictionary containing subdicts enumerating compute lists for each possible BSSE treatment.\n Subdict keys are n-body levels and values are sets of all the `mc_(frag, bas)` indices\n needed to compute that n-body level. A given index can appear multiple times within a\n subdict and among subdicts.\n\n compute_dict[\"cp\"] = {\n 1: set(),\n 2: {((1,), (1, 2)),\n ((2,), (1, 2)),\n ((1, 2), (1, 2))}\n }\n\n Subdicts below are always returned. Any may be empty if not requested through *bsse_type*.\n\n * ``'all'`` |w---w| full list of computations required\n * ``'cp'`` |w---w| list of computations required for CP procedure\n * ``'nocp'`` |w---w| list of computations required for non-CP procedure\n * ``'vmfc_compute'`` |w---w| list of computations required for VMFC procedure\n * ``'vmfc_levels'`` |w---w| list of levels required for VMFC procedure\n\n \"\"\"\n\n include_supersystem = False\n if \"supersystem\" in nbodies:\n if supersystem_max_nbody is None:\n raise ValueError(\"supersystem_max_nbody must be provided if 'supersystem' contains nbodies\")\n\n include_supersystem = True\n nbodies = list(nbodies)\n nbodies.remove(\"supersystem\")\n\n # What levels do we need?\n fragment_range = range(1, nfragments + 1)\n\n # Need nbodies and all lower-body in full basis\n cp_compute_list = {x: set() for x in nbodies}\n nocp_compute_list = {x: set() for x in nbodies}\n vmfc_compute_list = {x: set() for x in nbodies}\n vmfc_level_list = {x: set() for x in nbodies} # Need to sum something slightly different\n\n # Verify proper passing of bsse_type. already validated in Computer\n bsse_type_remainder = set(bsse_type) - {e.value for e in BsseEnum}\n if bsse_type_remainder:\n raise RuntimeError(f\"Unrecognized BSSE type(s): {bsse_type_remainder}\")\n\n # Build up compute sets\n if \"cp\" in bsse_type:\n # Everything is in counterpoise/nfr-mer basis\n basis_tuple = tuple(fragment_range)\n\n if supersystem_ie_only:\n for sublevel in [1, nfragments]:\n for x in itertools.combinations(fragment_range, sublevel):\n cp_compute_list[nfragments].add((x, basis_tuple))\n else:\n for nb in nbodies:\n # Note A.1: nb=1 is skipped because the nfr-mer-basis monomer\n # contributions cancel at 1-body. These skipped tasks will be\n # ordered anyways if higher bodies are requested. Monomers for\n # the purpose of total energies use monomer basis, not these\n # skipped tasks. See coordinating Note A.2 .\n if nb > 1:\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n cp_compute_list[nb].add((x, basis_tuple))\n\n if \"nocp\" in bsse_type:\n # Everything in natural/n-mer basis\n if supersystem_ie_only:\n for sublevel in [1, nfragments]:\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nfragments].add((x, x))\n else:\n for nb in nbodies:\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nb].add((x, x))\n\n if \"vmfc\" in bsse_type:\n # Like a CP for all combinations of pairs or greater\n for nb in nbodies:\n for cp_combos in itertools.combinations(fragment_range, nb):\n basis_tuple = tuple(cp_combos)\n for interior_nbody in range(1, nb + 1):\n for x in itertools.combinations(cp_combos, interior_nbody):\n combo_tuple = (x, basis_tuple)\n vmfc_compute_list[nb].add(combo_tuple)\n vmfc_level_list[len(basis_tuple)].add(combo_tuple)\n\n if return_total_data and 1 in nbodies:\n # Monomers in monomer basis\n nocp_compute_list.setdefault(1, set())\n for ifr in fragment_range:\n nocp_compute_list[1].add(((ifr,), (ifr,)))\n\n if include_supersystem:\n # Add supersystem info to the compute list (nocp only)\n for nb in range(1, supersystem_max_nbody + 1):\n cp_compute_list.setdefault(nb, set())\n nocp_compute_list.setdefault(nb, set())\n vmfc_compute_list.setdefault(nb, set())\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nb].add((x, x))\n\n # Add the total supersystem (nfragments@nfragments)\n nocp_compute_list.setdefault(nfragments, set())\n nocp_compute_list[nfragments].add((tuple(fragment_range), tuple(fragment_range)))\n\n # Build a comprehensive compute range\n # * do not use list length to count number of {nb}-body computations\n compute_list = {x: set() for x in nbodies}\n for nb in nbodies:\n compute_list[nb] |= cp_compute_list[nb]\n compute_list[nb] |= nocp_compute_list[nb]\n compute_list[nb] |= vmfc_compute_list[nb]\n\n if include_supersystem:\n for nb, lst in nocp_compute_list.items():\n compute_list.setdefault(nb, set())\n compute_list[nb] |= lst\n\n # Rearrange compute_list from key nb having values to compute all of that nb\n # to key nb including values of that nb. Use for counting.\n compute_list_count = {x: set() for x in nbodies}\n for nb in nbodies:\n for nbset in compute_list.values():\n for item in nbset:\n if len(item[0]) == nb:\n compute_list_count[nb].add(item)\n\n compute_dict = {\n \"all\": compute_list,\n \"cp\": cp_compute_list,\n \"nocp\": nocp_compute_list,\n \"vmfc_compute\": vmfc_compute_list,\n \"vmfc_levels\": vmfc_level_list,\n }\n return compute_dict\n
"},{"location":"reference/#qcmanybody.ManyBodyCalculator.analyze","title":"analyze(component_results)
","text":"Parameters:
Name Type Description Default component_results
Dict[str, Dict[str, Union[float, ndarray]]]
Nested dictionary with results from all individual molecular system calculations, including all subsystem/basis combinations, all model chemistries, and all properties (e.g., e/g/h).
For example, the below is the format for a nocp gradient run on a helium dimer with 1-body at CCSD and 2-body at MP2. The outer string key can be generated with the qcmanybody.utils.labeler
function. The inner string key is any property; QCManyBody presently knows how to process energy/gradient/Hessian.
{'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])}, '[\"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]])}, }
required Return Source code in qcmanybody/manybody.py
def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties = set([\"energy\"])\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"reference/#qcmanybody.ManyBodyCalculator.iterate_molecules","title":"iterate_molecules()
","text":"Iterate over all the molecules needed for the computation.
Yields model chemistry, label, and molecule.
Source code in qcmanybody/manybody.py
def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n
"},{"location":"tutorials/","title":"Tutorials","text":"This part of the project documentation focuses on a learning-oriented approach. You'll learn how to get started with the code in this project.
Note: Expand this section by considering the following points:
- Help newcomers with getting started
- Teach readers about your library by making them write code
- Inspire confidence through examples that work for everyone, repeatably
- Give readers an immediate sense of achievement
- Show concrete examples, no abstractions
- Provide the minimum necessary explanation
- Avoid any distractions
"}]}
\ No newline at end of file
+{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Welcome to MkDocs","text":"For full documentation visit mkdocs.org.
"},{"location":"#commands","title":"Commands","text":" mkdocs new [dir-name]
- Create a new project. mkdocs serve
- Start the live-reloading docs server. mkdocs build
- Build the documentation site. mkdocs -h
- Print help message and exit.
"},{"location":"#table-of-contents","title":"Table Of Contents","text":" - Tutorials
- How-To Guides
- Reference
- Explanation
"},{"location":"api/","title":"API Documentation","text":"Source code in qcmanybody/manybody.py
class ManyBodyCalculator:\n def __init__(\n self,\n molecule: Molecule,\n bsse_type: Sequence[BsseEnum],\n levels: Mapping[Union[int, Literal[\"supersystem\"]], str],\n return_total_data: bool,\n supersystem_ie_only: bool,\n ):\n # TODO\n self.embedding_charges = {}\n\n self.molecule = molecule\n self.bsse_type = [BsseEnum(x) for x in bsse_type]\n self.return_total_data = return_total_data\n self.supersystem_ie_only = supersystem_ie_only\n self.nfragments = len(molecule.fragments)\n\n self.levels = levels\n\n # Levels without supersystem\n self.levels_no_ss = {int(k): v for k, v in levels.items() if k != \"supersystem\"}\n\n # Just a set of all the modelchems\n self.mc_levels = set(self.levels.values())\n\n self.max_nbody = max(self.levels_no_ss.keys())\n\n if len(self.bsse_type) == 0:\n raise ValueError(\"No BSSE correction specified\")\n\n if BsseEnum.vmfc in self.bsse_type and len(set(self.levels.values())) == 1:\n # For single-modelchem VMFC, NOCP & sometimes CP are produced for free\n if BsseEnum.nocp not in self.bsse_type:\n self.bsse_type.append(BsseEnum.nocp)\n if BsseEnum.cp not in self.bsse_type and self.max_nbody == self.nfragments:\n self.bsse_type.append(BsseEnum.cp)\n\n self.return_bsse_type = self.bsse_type[0]\n\n ###############################\n # Build nbodies_per_mc_level\n # TODO - use Lori's code\n # TODO - dict to list of lists to handle non-contiguous levels\n # TODO multilevel and supersystem_ie_only=T not allowed together\n # TODO supersystem in levels is not to be trusted -- nfrag only and skips levels\n max_level = max(self.levels_no_ss.keys())\n\n if set(range(1, max_level + 1)) != set(self.levels_no_ss.keys()):\n raise ValueError(f\"Levels must be contiguous from 1 to {max_level}\")\n\n self.nbodies_per_mc_level: Dict[str, list] = {mc_level: [] for mc_level in self.mc_levels}\n for k, v in self.levels_no_ss.items():\n self.nbodies_per_mc_level[v].append(k)\n\n self.nbodies_per_mc_level = {k: sorted(v) for k, v in self.nbodies_per_mc_level.items()}\n\n # Supersystem is always at the end\n if \"supersystem\" in levels:\n ss_mc = levels[\"supersystem\"]\n self.nbodies_per_mc_level[ss_mc].append(\"supersystem\")\n\n # To be built on the fly\n self.mc_compute_dict = None\n\n if not np.array_equal(np.concatenate(self.molecule.fragments), np.arange(len(self.molecule.symbols))):\n raise ValueError(\"\"\"QCManyBody: non-contiguous fragments could be implemented but aren't at present\"\"\")\n\n # Build size and slices dictionaries. Assumes fragments are contiguous\n self.fragment_size_dict = {}\n self.fragment_slice_dict = {}\n iat = 0\n for ifr in range(1, self.nfragments + 1):\n nat = len(self.molecule.fragments[ifr - 1])\n self.fragment_size_dict[ifr] = nat\n self.fragment_slice_dict[ifr] = slice(iat, iat + nat)\n iat += nat\n\n @property\n def has_supersystem(self) -> bool:\n return \"supersystem\" in self.levels\n\n @property\n def compute_map(self) -> Dict[str, Dict[str, Dict[int, Set[FragBasIndex]]]]:\n if self.mc_compute_dict is not None:\n return self.mc_compute_dict\n\n # Build the compute lists\n self.mc_compute_dict = {}\n\n for mc in self.mc_levels:\n nbodies = self.nbodies_per_mc_level[mc]\n self.mc_compute_dict[mc] = build_nbody_compute_list(\n self.bsse_type,\n self.nfragments,\n nbodies,\n self.return_total_data,\n self.supersystem_ie_only,\n self.max_nbody,\n )\n\n return self.mc_compute_dict\n\n def resize_gradient(self, grad: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_gradient(grad, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def resize_hessian(self, hess: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_hessian(hess, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n\n def _assemble_nbody_components(\n self,\n property_label: str,\n component_results: Dict[str, Union[float, np.ndarray]],\n ) -> Dict[str, Any]:\n \"\"\"Assembles N-body components for a single derivative level and a single model chemistry level\n into interaction quantities according to requested BSSE treatment(s).\n \"\"\"\n\n # which level are we assembling?\n delabeled = [delabeler(k) for k in component_results.keys()]\n mc_level_labels = {x[0] for x in delabeled}\n\n if len(mc_level_labels) != 1:\n raise RuntimeError(f\"Multiple model chemistries passed into _assemble_nbody_components: {mc_level_labels}\")\n\n mc_level = mc_level_labels.pop()\n if mc_level not in self.mc_levels:\n raise RuntimeError(f\"Model chemistry {mc_level} not found in {self.mc_levels}\")\n\n # get the range of nbodies and the required calculations for this level\n bsse_type = self.bsse_type\n return_bsse_type = self.return_bsse_type\n nbodies = self.nbodies_per_mc_level[mc_level]\n if \"supersystem\" in nbodies:\n nbodies = list(range(1, self.max_nbody + 1))\n bsse_type = [BsseEnum.nocp]\n return_bsse_type = BsseEnum.nocp\n\n max_nbody = max(nbodies)\n compute_dict = self.compute_map[mc_level]\n\n if not all_same_shape(component_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(component_results.keys()))\n property_shape = find_shape(component_results[first_key])\n\n # Accumulation dictionaries\n # * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP\n # & CP the summed total energies (or other property) of each nb-body. That is:\n # * NOCP: {1: 1b@1b, 2: 2b@2b, ..., max_nbody: max_nbody-b@max_nbody-b} and\n # * CP: {1: 1b@nfr-b, 2: 2b@nfr-b, ..., max_nbody: max_nbody-b@nfr-b}.\n # VMFC bookkeeping is different. For key 1 it contains the summed 1b total energies.\n # But for higher keys, it contains each nb-body (non-additive) contribution to the energy.\n # * VMFC: {1: 1b@1b, 2: 2b contrib, ..., max_nbody: max_nbody-b contrib}\n cp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # * {bsse_type}_body_dict is usually filled with total energies (or other property).\n # Multiple model chemistry levels may be involved.\n # Generally, all consecutive keys between 1 and max_nbody will be present in the body_dict,\n # but if supersystem_ie_only=T, only 1b and nfr-b are present, or if \"supersystem\" in levels, ???\n # * TOT: {1: 1b@1b, 2: 2b tot prop with bsse_type treatment, ..., max_nbody: max_nbody-b tot prop with bsse_type treatment}\n # If 1b@1b (monomers in monomer basis) aren't available, which can happen when return_total_data=F\n # and 1b@1b aren't otherwise needed, body_dict contains interaction energies (or other property).\n # * IE: {1: shaped_zero, 2: 2b interaction prop using bsse_type, ..., max_nbody: max_nbody-b interaction prop using bsse_type}\n # For both TOT and IE cases, body_dict values are cummulative, not additive. For TOT, total,\n # interaction, and contribution data in ManyBodyResultProperties can be computed in\n # collect_vars. For IE, interaction and contribution data can be computed.\n cp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # Sum up all of the levels\n # * compute_dict[bt][nb] holds all the computations needed to compute nb\n # *not* all the nb-level computations, so build the latter\n cp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n nocp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n\n for nb in nbodies:\n for v in compute_dict[\"cp\"][nb]:\n if len(v[1]) != 1:\n cp_compute_list[len(v[0])].add(v)\n for w in compute_dict[\"nocp\"][nb]:\n nocp_compute_list[len(w[0])].add(w)\n\n for nb in range(1, nbodies[-1] + 1):\n cp_by_level[nb] = sum_cluster_data(component_results, cp_compute_list[nb], mc_level)\n nocp_by_level[nb] = sum_cluster_data(component_results, nocp_compute_list[nb], mc_level)\n if nb in compute_dict[\"vmfc_levels\"]:\n vmfc_by_level[nb] = sum_cluster_data(\n component_results, compute_dict[\"vmfc_levels\"][nb], mc_level, vmfc=True, nb=nb\n )\n\n # Extract data for monomers in monomer basis for CP total data\n if 1 in nbodies:\n monomers_in_monomer_basis = [v for v in compute_dict[\"nocp\"][1] if len(v[1]) == 1]\n monomer_sum = sum_cluster_data(component_results, set(monomers_in_monomer_basis), mc_level)\n else:\n monomer_sum = shaped_zero(property_shape)\n\n # Compute cp\n if BsseEnum.cp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n cp_body_dict[nb] = cp_by_level[nb] - bsse\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n cp_body_dict[nb] += take_nk * sign * cp_by_level[k]\n\n if nb == 1:\n bsse = cp_body_dict[nb] - monomer_sum\n cp_body_dict[nb] = copy_value(monomer_sum)\n else:\n cp_body_dict[nb] -= bsse\n\n # Compute nocp\n if BsseEnum.nocp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n nocp_body_dict[nb] = nocp_by_level[nb]\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n nocp_body_dict[nb] += take_nk * sign * nocp_by_level[k]\n\n # Compute vmfc\n if BsseEnum.vmfc in bsse_type:\n for nb in nbodies:\n for k in range(1, nb + 1):\n vmfc_body_dict[nb] += vmfc_by_level[k]\n\n # Collect specific and generalized returns\n results = {\n f\"cp_{property_label}_body_dict\": cp_body_dict,\n f\"nocp_{property_label}_body_dict\": nocp_body_dict,\n f\"vmfc_{property_label}_body_dict\": vmfc_body_dict,\n }\n\n # Overall return body dict & value for this property\n results[f\"{property_label}_body_dict\"] = results[f\"{return_bsse_type.value}_{property_label}_body_dict\"]\n results[f\"ret_{property_label}\"] = copy_value(results[f\"{property_label}_body_dict\"][max_nbody])\n\n if not self.return_total_data:\n results[f\"ret_{property_label}\"] -= results[f\"{property_label}_body_dict\"][1]\n\n return results\n\n def _analyze(\n self,\n property_label: str,\n property_results: Dict[str, Union[float, np.ndarray]], # Label to results\n ):\n # Initialize with zeros\n if not all_same_shape(property_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(property_results.keys()))\n property_shape = find_shape(property_results[first_key])\n\n property_result = shaped_zero(property_shape)\n property_body_dict = {bt.value: {} for bt in self.bsse_type}\n property_body_contribution = {bt.value: {} for bt in self.bsse_type}\n\n # results per model chemistry\n mc_results = {}\n species_results = {}\n\n # sort by nbody level, ignore supersystem\n sorted_nbodies = [(k, v) for k, v in self.nbodies_per_mc_level.items() if v != [\"supersystem\"]]\n sorted_nbodies = sorted(sorted_nbodies, reverse=True, key=lambda x: x[1])\n for mc_label, nbody_list in sorted_nbodies:\n # filter to only one model chemistry\n filtered_results = {k: v for k, v in property_results.items() if delabeler(k)[0] == mc_label}\n\n if not filtered_results:\n if nbody_list == [1]:\n # Note A.2: Note A.1 holds, but for the special case of CP-only\n # and rtd=False and multilevel with a separate level for\n # 1-body, the skipped tasks run afoul of sanity checks, so\n # we'll add a dummy result.\n filtered_results = {labeler(mc_label, [1000], [1000]): shaped_zero(property_shape)}\n else:\n raise RuntimeError(f\"No data found for model chemistry {mc_label}\")\n\n nb_component_results = self._assemble_nbody_components(property_label, filtered_results)\n mc_results[mc_label] = nb_component_results\n\n for n in nbody_list[::-1]:\n property_bsse_dict = {bt.value: shaped_zero(property_shape) for bt in self.bsse_type}\n\n for m in range(n - 1, n + 1):\n if m == 0:\n continue\n\n # Subtract the (n-1)-body contribution from the n-body contribution to get the n-body effect\n sign = (-1) ** (1 - m // n)\n for bt in self.bsse_type:\n property_bsse_dict[bt.value] += (\n sign * mc_results[mc_label][f\"{bt.value}_{property_label}_body_dict\"][m]\n )\n\n property_result += property_bsse_dict[self.return_bsse_type]\n for bt in self.bsse_type:\n property_body_contribution[bt.value][n] = property_bsse_dict[bt.value]\n\n if self.has_supersystem:\n # Get the MC label for supersystem tasks\n supersystem_mc_level = self.levels[\"supersystem\"]\n\n # Super system recovers higher order effects at a lower level\n frag_range = tuple(range(1, self.nfragments + 1))\n\n ss_cresults = {k: v for k, v in property_results.items() if delabeler(k)[0] == supersystem_mc_level}\n ss_component_results = self._assemble_nbody_components(property_label, ss_cresults)\n mc_results[supersystem_mc_level] = ss_component_results\n\n # Compute components at supersystem level of theory\n ss_label = labeler(supersystem_mc_level, frag_range, frag_range)\n supersystem_result = property_results[ss_label]\n property_result += supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n\n for bt in self.bsse_type:\n property_body_contribution[bt][self.nfragments] = (\n supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n )\n\n for bt in self.bsse_type:\n bstr = bt.value\n for n in property_body_contribution[bstr]:\n property_body_dict[bstr][n] = sum(\n [\n property_body_contribution[bstr][i]\n for i in range(1, n + 1)\n if i in property_body_contribution[bstr]\n ]\n )\n\n if not self.return_total_data:\n # Remove monomer contribution for interaction data\n property_result -= property_body_dict[self.return_bsse_type][1]\n\n nbody_results = {\n f\"ret_{property_label}\": property_result,\n f\"{property_label}_body_dict\": property_body_dict,\n \"mc_results\": mc_results,\n }\n return nbody_results\n\n def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties: Set[str] = {\"energy\"}\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"api/#qcmanybody.ManyBodyCalculator.analyze","title":"analyze(component_results)
","text":"Parameters:
Name Type Description Default component_results
Dict[str, Dict[str, Union[float, ndarray]]]
Nested dictionary with results from all individual molecular system calculations, including all subsystem/basis combinations, all model chemistries, and all properties (e.g., e/g/h).
For example, the below is the format for a nocp gradient run on a helium dimer with 1-body at CCSD and 2-body at MP2. The outer string key can be generated with the qcmanybody.utils.labeler
function. The inner string key is any property; QCManyBody presently knows how to process energy/gradient/Hessian.
{'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])}, '[\"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]])}, }
required Return Source code in qcmanybody/manybody.py
def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties: Set[str] = {\"energy\"}\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"api/#qcmanybody.ManyBodyCalculator.iterate_molecules","title":"iterate_molecules()
","text":"Iterate over all the molecules needed for the computation.
Yields model chemistry, label, and molecule.
Source code in qcmanybody/manybody.py
def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n
"},{"location":"explanation/","title":"Explanation","text":"This part of the project documentation focuses on an understanding-oriented approach. You'll get a chance to read about the background of the project, as well as reasoning about how it was implemented.
Note: Expand this section by considering the following points:
- Give context and background on your library
- Explain why you created it
- Provide multiple examples and approaches of how to work with it
- Help the reader make connections
- Avoid writing instructions or technical descriptions here
"},{"location":"reference/","title":"Reference","text":"This part of the project documentation focuses on an information-oriented approach. Use it as a reference for the technical implementation of the QCManyBody
project code.
Source code in qcmanybody/manybody.py
class ManyBodyCalculator:\n def __init__(\n self,\n molecule: Molecule,\n bsse_type: Sequence[BsseEnum],\n levels: Mapping[Union[int, Literal[\"supersystem\"]], str],\n return_total_data: bool,\n supersystem_ie_only: bool,\n ):\n # TODO\n self.embedding_charges = {}\n\n self.molecule = molecule\n self.bsse_type = [BsseEnum(x) for x in bsse_type]\n self.return_total_data = return_total_data\n self.supersystem_ie_only = supersystem_ie_only\n self.nfragments = len(molecule.fragments)\n\n self.levels = levels\n\n # Levels without supersystem\n self.levels_no_ss = {int(k): v for k, v in levels.items() if k != \"supersystem\"}\n\n # Just a set of all the modelchems\n self.mc_levels = set(self.levels.values())\n\n self.max_nbody = max(self.levels_no_ss.keys())\n\n if len(self.bsse_type) == 0:\n raise ValueError(\"No BSSE correction specified\")\n\n if BsseEnum.vmfc in self.bsse_type and len(set(self.levels.values())) == 1:\n # For single-modelchem VMFC, NOCP & sometimes CP are produced for free\n if BsseEnum.nocp not in self.bsse_type:\n self.bsse_type.append(BsseEnum.nocp)\n if BsseEnum.cp not in self.bsse_type and self.max_nbody == self.nfragments:\n self.bsse_type.append(BsseEnum.cp)\n\n self.return_bsse_type = self.bsse_type[0]\n\n ###############################\n # Build nbodies_per_mc_level\n # TODO - use Lori's code\n # TODO - dict to list of lists to handle non-contiguous levels\n # TODO multilevel and supersystem_ie_only=T not allowed together\n # TODO supersystem in levels is not to be trusted -- nfrag only and skips levels\n max_level = max(self.levels_no_ss.keys())\n\n if set(range(1, max_level + 1)) != set(self.levels_no_ss.keys()):\n raise ValueError(f\"Levels must be contiguous from 1 to {max_level}\")\n\n self.nbodies_per_mc_level: Dict[str, list] = {mc_level: [] for mc_level in self.mc_levels}\n for k, v in self.levels_no_ss.items():\n self.nbodies_per_mc_level[v].append(k)\n\n self.nbodies_per_mc_level = {k: sorted(v) for k, v in self.nbodies_per_mc_level.items()}\n\n # Supersystem is always at the end\n if \"supersystem\" in levels:\n ss_mc = levels[\"supersystem\"]\n self.nbodies_per_mc_level[ss_mc].append(\"supersystem\")\n\n # To be built on the fly\n self.mc_compute_dict = None\n\n if not np.array_equal(np.concatenate(self.molecule.fragments), np.arange(len(self.molecule.symbols))):\n raise ValueError(\"\"\"QCManyBody: non-contiguous fragments could be implemented but aren't at present\"\"\")\n\n # Build size and slices dictionaries. Assumes fragments are contiguous\n self.fragment_size_dict = {}\n self.fragment_slice_dict = {}\n iat = 0\n for ifr in range(1, self.nfragments + 1):\n nat = len(self.molecule.fragments[ifr - 1])\n self.fragment_size_dict[ifr] = nat\n self.fragment_slice_dict[ifr] = slice(iat, iat + nat)\n iat += nat\n\n @property\n def has_supersystem(self) -> bool:\n return \"supersystem\" in self.levels\n\n @property\n def compute_map(self) -> Dict[str, Dict[str, Dict[int, Set[FragBasIndex]]]]:\n if self.mc_compute_dict is not None:\n return self.mc_compute_dict\n\n # Build the compute lists\n self.mc_compute_dict = {}\n\n for mc in self.mc_levels:\n nbodies = self.nbodies_per_mc_level[mc]\n self.mc_compute_dict[mc] = build_nbody_compute_list(\n self.bsse_type,\n self.nfragments,\n nbodies,\n self.return_total_data,\n self.supersystem_ie_only,\n self.max_nbody,\n )\n\n return self.mc_compute_dict\n\n def resize_gradient(self, grad: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_gradient(grad, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def resize_hessian(self, hess: np.ndarray, bas: Tuple[int, ...], *, reverse: bool = False) -> np.ndarray:\n return resize_hessian(hess, bas, self.fragment_size_dict, self.fragment_slice_dict, reverse=reverse)\n\n def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n\n def _assemble_nbody_components(\n self,\n property_label: str,\n component_results: Dict[str, Union[float, np.ndarray]],\n ) -> Dict[str, Any]:\n \"\"\"Assembles N-body components for a single derivative level and a single model chemistry level\n into interaction quantities according to requested BSSE treatment(s).\n \"\"\"\n\n # which level are we assembling?\n delabeled = [delabeler(k) for k in component_results.keys()]\n mc_level_labels = {x[0] for x in delabeled}\n\n if len(mc_level_labels) != 1:\n raise RuntimeError(f\"Multiple model chemistries passed into _assemble_nbody_components: {mc_level_labels}\")\n\n mc_level = mc_level_labels.pop()\n if mc_level not in self.mc_levels:\n raise RuntimeError(f\"Model chemistry {mc_level} not found in {self.mc_levels}\")\n\n # get the range of nbodies and the required calculations for this level\n bsse_type = self.bsse_type\n return_bsse_type = self.return_bsse_type\n nbodies = self.nbodies_per_mc_level[mc_level]\n if \"supersystem\" in nbodies:\n nbodies = list(range(1, self.max_nbody + 1))\n bsse_type = [BsseEnum.nocp]\n return_bsse_type = BsseEnum.nocp\n\n max_nbody = max(nbodies)\n compute_dict = self.compute_map[mc_level]\n\n if not all_same_shape(component_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(component_results.keys()))\n property_shape = find_shape(component_results[first_key])\n\n # Accumulation dictionaries\n # * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP\n # & CP the summed total energies (or other property) of each nb-body. That is:\n # * NOCP: {1: 1b@1b, 2: 2b@2b, ..., max_nbody: max_nbody-b@max_nbody-b} and\n # * CP: {1: 1b@nfr-b, 2: 2b@nfr-b, ..., max_nbody: max_nbody-b@nfr-b}.\n # VMFC bookkeeping is different. For key 1 it contains the summed 1b total energies.\n # But for higher keys, it contains each nb-body (non-additive) contribution to the energy.\n # * VMFC: {1: 1b@1b, 2: 2b contrib, ..., max_nbody: max_nbody-b contrib}\n cp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # * {bsse_type}_body_dict is usually filled with total energies (or other property).\n # Multiple model chemistry levels may be involved.\n # Generally, all consecutive keys between 1 and max_nbody will be present in the body_dict,\n # but if supersystem_ie_only=T, only 1b and nfr-b are present, or if \"supersystem\" in levels, ???\n # * TOT: {1: 1b@1b, 2: 2b tot prop with bsse_type treatment, ..., max_nbody: max_nbody-b tot prop with bsse_type treatment}\n # If 1b@1b (monomers in monomer basis) aren't available, which can happen when return_total_data=F\n # and 1b@1b aren't otherwise needed, body_dict contains interaction energies (or other property).\n # * IE: {1: shaped_zero, 2: 2b interaction prop using bsse_type, ..., max_nbody: max_nbody-b interaction prop using bsse_type}\n # For both TOT and IE cases, body_dict values are cummulative, not additive. For TOT, total,\n # interaction, and contribution data in ManyBodyResultProperties can be computed in\n # collect_vars. For IE, interaction and contribution data can be computed.\n cp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n nocp_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n vmfc_body_dict = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}\n\n # Sum up all of the levels\n # * compute_dict[bt][nb] holds all the computations needed to compute nb\n # *not* all the nb-level computations, so build the latter\n cp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n nocp_compute_list = {nb: set() for nb in range(1, nbodies[-1] + 1)}\n\n for nb in nbodies:\n for v in compute_dict[\"cp\"][nb]:\n if len(v[1]) != 1:\n cp_compute_list[len(v[0])].add(v)\n for w in compute_dict[\"nocp\"][nb]:\n nocp_compute_list[len(w[0])].add(w)\n\n for nb in range(1, nbodies[-1] + 1):\n cp_by_level[nb] = sum_cluster_data(component_results, cp_compute_list[nb], mc_level)\n nocp_by_level[nb] = sum_cluster_data(component_results, nocp_compute_list[nb], mc_level)\n if nb in compute_dict[\"vmfc_levels\"]:\n vmfc_by_level[nb] = sum_cluster_data(\n component_results, compute_dict[\"vmfc_levels\"][nb], mc_level, vmfc=True, nb=nb\n )\n\n # Extract data for monomers in monomer basis for CP total data\n if 1 in nbodies:\n monomers_in_monomer_basis = [v for v in compute_dict[\"nocp\"][1] if len(v[1]) == 1]\n monomer_sum = sum_cluster_data(component_results, set(monomers_in_monomer_basis), mc_level)\n else:\n monomer_sum = shaped_zero(property_shape)\n\n # Compute cp\n if BsseEnum.cp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n cp_body_dict[nb] = cp_by_level[nb] - bsse\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n cp_body_dict[nb] += take_nk * sign * cp_by_level[k]\n\n if nb == 1:\n bsse = cp_body_dict[nb] - monomer_sum\n cp_body_dict[nb] = copy_value(monomer_sum)\n else:\n cp_body_dict[nb] -= bsse\n\n # Compute nocp\n if BsseEnum.nocp in bsse_type:\n for nb in range(1, nbodies[-1] + 1):\n if nb == self.nfragments:\n nocp_body_dict[nb] = nocp_by_level[nb]\n continue\n\n for k in range(1, nb + 1):\n take_nk = math.comb(self.nfragments - k - 1, nb - k)\n sign = (-1) ** (nb - k)\n nocp_body_dict[nb] += take_nk * sign * nocp_by_level[k]\n\n # Compute vmfc\n if BsseEnum.vmfc in bsse_type:\n for nb in nbodies:\n for k in range(1, nb + 1):\n vmfc_body_dict[nb] += vmfc_by_level[k]\n\n # Collect specific and generalized returns\n results = {\n f\"cp_{property_label}_body_dict\": cp_body_dict,\n f\"nocp_{property_label}_body_dict\": nocp_body_dict,\n f\"vmfc_{property_label}_body_dict\": vmfc_body_dict,\n }\n\n # Overall return body dict & value for this property\n results[f\"{property_label}_body_dict\"] = results[f\"{return_bsse_type.value}_{property_label}_body_dict\"]\n results[f\"ret_{property_label}\"] = copy_value(results[f\"{property_label}_body_dict\"][max_nbody])\n\n if not self.return_total_data:\n results[f\"ret_{property_label}\"] -= results[f\"{property_label}_body_dict\"][1]\n\n return results\n\n def _analyze(\n self,\n property_label: str,\n property_results: Dict[str, Union[float, np.ndarray]], # Label to results\n ):\n # Initialize with zeros\n if not all_same_shape(property_results.values()):\n raise ValueError(\"All values in data dictionary must have the same shape.\")\n\n # Use first data value to determine shape\n first_key = next(iter(property_results.keys()))\n property_shape = find_shape(property_results[first_key])\n\n property_result = shaped_zero(property_shape)\n property_body_dict = {bt.value: {} for bt in self.bsse_type}\n property_body_contribution = {bt.value: {} for bt in self.bsse_type}\n\n # results per model chemistry\n mc_results = {}\n species_results = {}\n\n # sort by nbody level, ignore supersystem\n sorted_nbodies = [(k, v) for k, v in self.nbodies_per_mc_level.items() if v != [\"supersystem\"]]\n sorted_nbodies = sorted(sorted_nbodies, reverse=True, key=lambda x: x[1])\n for mc_label, nbody_list in sorted_nbodies:\n # filter to only one model chemistry\n filtered_results = {k: v for k, v in property_results.items() if delabeler(k)[0] == mc_label}\n\n if not filtered_results:\n if nbody_list == [1]:\n # Note A.2: Note A.1 holds, but for the special case of CP-only\n # and rtd=False and multilevel with a separate level for\n # 1-body, the skipped tasks run afoul of sanity checks, so\n # we'll add a dummy result.\n filtered_results = {labeler(mc_label, [1000], [1000]): shaped_zero(property_shape)}\n else:\n raise RuntimeError(f\"No data found for model chemistry {mc_label}\")\n\n nb_component_results = self._assemble_nbody_components(property_label, filtered_results)\n mc_results[mc_label] = nb_component_results\n\n for n in nbody_list[::-1]:\n property_bsse_dict = {bt.value: shaped_zero(property_shape) for bt in self.bsse_type}\n\n for m in range(n - 1, n + 1):\n if m == 0:\n continue\n\n # Subtract the (n-1)-body contribution from the n-body contribution to get the n-body effect\n sign = (-1) ** (1 - m // n)\n for bt in self.bsse_type:\n property_bsse_dict[bt.value] += (\n sign * mc_results[mc_label][f\"{bt.value}_{property_label}_body_dict\"][m]\n )\n\n property_result += property_bsse_dict[self.return_bsse_type]\n for bt in self.bsse_type:\n property_body_contribution[bt.value][n] = property_bsse_dict[bt.value]\n\n if self.has_supersystem:\n # Get the MC label for supersystem tasks\n supersystem_mc_level = self.levels[\"supersystem\"]\n\n # Super system recovers higher order effects at a lower level\n frag_range = tuple(range(1, self.nfragments + 1))\n\n ss_cresults = {k: v for k, v in property_results.items() if delabeler(k)[0] == supersystem_mc_level}\n ss_component_results = self._assemble_nbody_components(property_label, ss_cresults)\n mc_results[supersystem_mc_level] = ss_component_results\n\n # Compute components at supersystem level of theory\n ss_label = labeler(supersystem_mc_level, frag_range, frag_range)\n supersystem_result = property_results[ss_label]\n property_result += supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n\n for bt in self.bsse_type:\n property_body_contribution[bt][self.nfragments] = (\n supersystem_result - ss_component_results[f\"{property_label}_body_dict\"][self.max_nbody]\n )\n\n for bt in self.bsse_type:\n bstr = bt.value\n for n in property_body_contribution[bstr]:\n property_body_dict[bstr][n] = sum(\n [\n property_body_contribution[bstr][i]\n for i in range(1, n + 1)\n if i in property_body_contribution[bstr]\n ]\n )\n\n if not self.return_total_data:\n # Remove monomer contribution for interaction data\n property_result -= property_body_dict[self.return_bsse_type][1]\n\n nbody_results = {\n f\"ret_{property_label}\": property_result,\n f\"{property_label}_body_dict\": property_body_dict,\n \"mc_results\": mc_results,\n }\n return nbody_results\n\n def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties: Set[str] = {\"energy\"}\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
Generates lists of N-Body computations needed for requested BSSE treatments.
Parameters:
Name Type Description Default bsse_type
Iterable[BsseEnum]
Requested BSSE treatments.
required nfragments
int
Number of distinct fragments comprising the full molecular supersystem.
required nbodies
Iterable[Union[int, Literal['supersystem']]]
List of n-body levels (e.g., [2]
or [1, 2]
or [\"supersystem\"]
) for which to generate tasks. Note the natural 1-indexing, so [1]
covers one-body contributions.
required return_total_data
bool
Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested, as opposed to interaction data (False).
required supersystem_ie_only
bool
Target the supersystem total/interaction energy (IE) data over the many-body expansion (MBE) \" analysis, thereby omitting intermediate-body calculations.
required supersystem_max_nbody
Optional[int]
Maximum n-body to use for a supersystem calculation. Must be specified if \"supersystem\" is in nbodies
None
Returns:
Type Description compute_dict
Dictionary containing subdicts enumerating compute lists for each possible BSSE treatment. Subdict keys are n-body levels and values are sets of all the mc_(frag, bas)
indices needed to compute that n-body level. A given index can appear multiple times within a subdict and among subdicts.
compute_dict[\"cp\"] = {\n 1: set(),\n 2: {((1,), (1, 2)),\n ((2,), (1, 2)),\n ((1, 2), (1, 2))}\n}\n
Subdicts below are always returned. Any may be empty if not requested through bsse_type.
'all'
|w---w| full list of computations required 'cp'
|w---w| list of computations required for CP procedure 'nocp'
|w---w| list of computations required for non-CP procedure 'vmfc_compute'
|w---w| list of computations required for VMFC procedure 'vmfc_levels'
|w---w| list of levels required for VMFC procedure
Source code in qcmanybody/builder.py
def build_nbody_compute_list(\n bsse_type: Iterable[BsseEnum],\n nfragments: int,\n nbodies: Iterable[Union[int, Literal[\"supersystem\"]]],\n return_total_data: bool,\n supersystem_ie_only: bool,\n supersystem_max_nbody: Optional[int] = None,\n) -> Dict[str, Dict[int, Set[FragBasIndex]]]:\n \"\"\"Generates lists of N-Body computations needed for requested BSSE treatments.\n\n Parameters\n ----------\n bsse_type\n Requested BSSE treatments.\n nfragments\n Number of distinct fragments comprising the full molecular supersystem.\n nbodies\n List of n-body levels (e.g., `[2]` or `[1, 2]` or `[\"supersystem\"]`) for which to generate tasks.\n Note the natural 1-indexing, so `[1]` covers one-body contributions.\n return_total_data\n Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested,\n as opposed to interaction data (False).\n supersystem_ie_only\n Target the supersystem total/interaction energy (IE) data over the many-body expansion (MBE) \"\n analysis, thereby omitting intermediate-body calculations.\n supersystem_max_nbody\n Maximum n-body to use for a supersystem calculation. Must be specified if \"supersystem\" is in `nbodies`\n\n Returns\n -------\n compute_dict\n Dictionary containing subdicts enumerating compute lists for each possible BSSE treatment.\n Subdict keys are n-body levels and values are sets of all the `mc_(frag, bas)` indices\n needed to compute that n-body level. A given index can appear multiple times within a\n subdict and among subdicts.\n\n compute_dict[\"cp\"] = {\n 1: set(),\n 2: {((1,), (1, 2)),\n ((2,), (1, 2)),\n ((1, 2), (1, 2))}\n }\n\n Subdicts below are always returned. Any may be empty if not requested through *bsse_type*.\n\n * ``'all'`` |w---w| full list of computations required\n * ``'cp'`` |w---w| list of computations required for CP procedure\n * ``'nocp'`` |w---w| list of computations required for non-CP procedure\n * ``'vmfc_compute'`` |w---w| list of computations required for VMFC procedure\n * ``'vmfc_levels'`` |w---w| list of levels required for VMFC procedure\n\n \"\"\"\n\n include_supersystem = False\n if \"supersystem\" in nbodies:\n if supersystem_max_nbody is None:\n raise ValueError(\"supersystem_max_nbody must be provided if 'supersystem' contains nbodies\")\n\n include_supersystem = True\n\n nbodies: List[int] = [x for x in nbodies if x != \"supersystem\"]\n\n # What levels do we need?\n fragment_range = range(1, nfragments + 1)\n\n # Need nbodies and all lower-body in full basis\n cp_compute_list = {x: set() for x in nbodies}\n nocp_compute_list = {x: set() for x in nbodies}\n vmfc_compute_list = {x: set() for x in nbodies}\n vmfc_level_list = {x: set() for x in nbodies} # Need to sum something slightly different\n\n # Verify proper passing of bsse_type. already validated in Computer\n bsse_type_remainder = set(bsse_type) - {e.value for e in BsseEnum}\n if bsse_type_remainder:\n raise RuntimeError(f\"Unrecognized BSSE type(s): {bsse_type_remainder}\")\n\n # Build up compute sets\n if \"cp\" in bsse_type:\n # Everything is in counterpoise/nfr-mer basis\n basis_tuple = tuple(fragment_range)\n\n if supersystem_ie_only:\n for sublevel in [1, nfragments]:\n for x in itertools.combinations(fragment_range, sublevel):\n cp_compute_list[nfragments].add((x, basis_tuple))\n else:\n for nb in nbodies:\n # Note A.1: nb=1 is skipped because the nfr-mer-basis monomer\n # contributions cancel at 1-body. These skipped tasks will be\n # ordered anyways if higher bodies are requested. Monomers for\n # the purpose of total energies use monomer basis, not these\n # skipped tasks. See coordinating Note A.2 .\n if nb > 1:\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n cp_compute_list[nb].add((x, basis_tuple))\n\n if \"nocp\" in bsse_type:\n # Everything in natural/n-mer basis\n if supersystem_ie_only:\n for sublevel in [1, nfragments]:\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nfragments].add((x, x))\n else:\n for nb in nbodies:\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nb].add((x, x))\n\n if \"vmfc\" in bsse_type:\n # Like a CP for all combinations of pairs or greater\n for nb in nbodies:\n for cp_combos in itertools.combinations(fragment_range, nb):\n basis_tuple = tuple(cp_combos)\n for interior_nbody in range(1, nb + 1):\n for x in itertools.combinations(cp_combos, interior_nbody):\n combo_tuple = (x, basis_tuple)\n vmfc_compute_list[nb].add(combo_tuple)\n vmfc_level_list[len(basis_tuple)].add(combo_tuple)\n\n if return_total_data and 1 in nbodies:\n # Monomers in monomer basis\n nocp_compute_list.setdefault(1, set())\n for ifr in fragment_range:\n nocp_compute_list[1].add(((ifr,), (ifr,)))\n\n if include_supersystem:\n # Add supersystem info to the compute list (nocp only)\n for nb in range(1, supersystem_max_nbody + 1):\n cp_compute_list.setdefault(nb, set())\n nocp_compute_list.setdefault(nb, set())\n vmfc_compute_list.setdefault(nb, set())\n for sublevel in range(1, nb + 1):\n for x in itertools.combinations(fragment_range, sublevel):\n nocp_compute_list[nb].add((x, x))\n\n # Add the total supersystem (nfragments@nfragments)\n nocp_compute_list.setdefault(nfragments, set())\n nocp_compute_list[nfragments].add((tuple(fragment_range), tuple(fragment_range)))\n\n # Build a comprehensive compute range\n # * do not use list length to count number of {nb}-body computations\n compute_list = {x: set() for x in nbodies}\n for nb in nbodies:\n compute_list[nb] |= cp_compute_list[nb]\n compute_list[nb] |= nocp_compute_list[nb]\n compute_list[nb] |= vmfc_compute_list[nb]\n\n if include_supersystem:\n for nb, lst in nocp_compute_list.items():\n compute_list.setdefault(nb, set())\n compute_list[nb] |= lst\n\n # Rearrange compute_list from key nb having values to compute all of that nb\n # to key nb including values of that nb. Use for counting.\n compute_list_count = {x: set() for x in nbodies}\n for nb in nbodies:\n for nbset in compute_list.values():\n for item in nbset:\n if len(item[0]) == nb:\n compute_list_count[nb].add(item)\n\n compute_dict = {\n \"all\": compute_list,\n \"cp\": cp_compute_list,\n \"nocp\": nocp_compute_list,\n \"vmfc_compute\": vmfc_compute_list,\n \"vmfc_levels\": vmfc_level_list,\n }\n return compute_dict\n
"},{"location":"reference/#qcmanybody.ManyBodyCalculator.analyze","title":"analyze(component_results)
","text":"Parameters:
Name Type Description Default component_results
Dict[str, Dict[str, Union[float, ndarray]]]
Nested dictionary with results from all individual molecular system calculations, including all subsystem/basis combinations, all model chemistries, and all properties (e.g., e/g/h).
For example, the below is the format for a nocp gradient run on a helium dimer with 1-body at CCSD and 2-body at MP2. The outer string key can be generated with the qcmanybody.utils.labeler
function. The inner string key is any property; QCManyBody presently knows how to process energy/gradient/Hessian.
{'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])}, '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])}, '[\"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]])}, }
required Return Source code in qcmanybody/manybody.py
def analyze(\n self,\n component_results: Dict[str, Dict[str, Union[float, np.ndarray]]],\n ):\n \"\"\"\n\n Parameters\n ----------\n component_results\n Nested dictionary with results from all individual molecular system\n calculations, including all subsystem/basis combinations, all model\n chemistries, and all properties (e.g., e/g/h).\n\n For example, the below is the format for a nocp gradient run on a\n helium dimer with 1-body at CCSD and 2-body at MP2. The outer string\n key can be generated with the ``qcmanybody.utils.labeler`` function.\n The inner string key is any property; QCManyBody presently knows how\n to process energy/gradient/Hessian.\n\n {'[\"ccsd\", [1], [1]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"ccsd\", [2], [2]]': {'energy': -2.87, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1], [1]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [2], [2]]': {'energy': -2.86, 'gradient': array([[0., 0., 0.]])},\n '[\"mp2\", [1, 2], [1, 2]]': {'energy': -5.73, 'gradient': array([[ 0., 0., 0.0053], [ 0., 0., -0.0053]])},\n }\n\n Return\n ------\n\n \"\"\"\n\n # All properties that were passed to us\n # * seed with \"energy\" so free/no-op jobs can process\n available_properties: Set[str] = {\"energy\"}\n for property_data in component_results.values():\n available_properties.update(property_data.keys())\n\n # reorganize to component_results_inv[property][label] = 1.23\n component_results_inv = {k: {} for k in available_properties}\n\n for cluster_label, property_data in component_results.items():\n for property_label, property_value in property_data.items():\n component_results_inv[property_label][cluster_label] = property_value\n\n # Remove any missing data\n component_results_inv = {k: v for k, v in component_results_inv.items() if v}\n if not component_results_inv:\n # Note B: Rarely, \"no results\" is expected, like for CP-only,\n # rtd=False, and max_nbody=1. We'll add a dummy entry so\n # processing can continue.\n component_results_inv[\"energy\"] = {'[\"dummy\", [1000], [1000]]': 0.0}\n\n # Actually analyze\n is_embedded = bool(self.embedding_charges)\n component_properties = defaultdict(dict)\n all_results = {}\n nbody_dict = {}\n# all_results[\"energy_body_dict\"] = {\"cp\": {1: 0.0}}\n\n for property_label, property_results in component_results_inv.items():\n # Expand gradient and hessian\n if property_label == \"gradient\":\n property_results = {k: self.resize_gradient(v, delabeler(k)[2]) for k, v in property_results.items()}\n if property_label == \"hessian\":\n property_results = {k: self.resize_hessian(v, delabeler(k)[2]) for k, v in property_results.items()}\n\n r = self._analyze(property_label, property_results)\n for k, v in property_results.items():\n component_properties[k][\"calcinfo_natom\"] = len(self.molecule.symbols)\n component_properties[k][f\"return_{property_label}\"] = v\n all_results.update(r)\n\n for bt in self.bsse_type:\n print_nbody_energy(\n all_results[\"energy_body_dict\"][bt],\n f\"{bt.upper()}-corrected multilevel many-body expansion\",\n self.nfragments,\n is_embedded,\n )\n\n for property_label in available_properties:\n for bt in self.bsse_type:\n if not self.has_supersystem: # skipped levels?\n nbody_dict.update(\n collect_vars(\n bt.upper(),\n property_label.upper(),\n all_results[f\"{property_label}_body_dict\"][bt],\n self.max_nbody,\n is_embedded,\n self.supersystem_ie_only,\n )\n )\n\n all_results[\"results\"] = nbody_dict\n all_results[\"component_properties\"] = component_properties\n\n # Make dictionary with \"1cp\", \"2cp\", etc\n ebd = all_results[\"energy_body_dict\"]\n all_results[\"energy_body_dict\"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}\n\n return all_results\n
"},{"location":"reference/#qcmanybody.ManyBodyCalculator.iterate_molecules","title":"iterate_molecules()
","text":"Iterate over all the molecules needed for the computation.
Yields model chemistry, label, and molecule.
Source code in qcmanybody/manybody.py
def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]:\n \"\"\"Iterate over all the molecules needed for the computation.\n\n Yields model chemistry, label, and molecule.\n \"\"\"\n\n done_molecules = set()\n\n for mc, compute_dict in self.compute_map.items():\n # TODO - this is a bit of a hack. Lots of duplication when reaching higher nbody\n for compute_list in compute_dict[\"all\"].values():\n for real_atoms, basis_atoms in compute_list:\n label = labeler(mc, real_atoms, basis_atoms)\n if label in done_molecules:\n continue\n\n ghost_atoms = list(set(basis_atoms) - set(real_atoms))\n\n # Shift to zero-indexing\n real_atoms_0 = [x - 1 for x in real_atoms]\n ghost_atoms_0 = [x - 1 for x in ghost_atoms]\n mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False)\n mol = mol.copy(update={\"fix_com\": True, \"fix_orientation\": True})\n\n # if self.embedding_charges:\n # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1]))\n # charges = []\n # for frag in embedding_frags:\n # positions = self.molecule.extract_subsets(frag).geometry().np.tolist()\n # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])])\n # data['keywords']['function_kwargs'].update({'external_potentials': charges})\n\n done_molecules.add(label)\n yield mc, label, mol\n
"},{"location":"tutorials/","title":"Tutorials","text":"This part of the project documentation focuses on a learning-oriented approach. You'll learn how to get started with the code in this project.
Note: Expand this section by considering the following points:
- Help newcomers with getting started
- Teach readers about your library by making them write code
- Inspire confidence through examples that work for everyone, repeatably
- Give readers an immediate sense of achievement
- Show concrete examples, no abstractions
- Provide the minimum necessary explanation
- Avoid any distractions
"}]}
\ No newline at end of file
diff --git a/sitemap.xml b/sitemap.xml
index 5f26d53..224d229 100644
--- a/sitemap.xml
+++ b/sitemap.xml
@@ -2,32 +2,32 @@
https://molssi.github.io/QCManyBody/
- 2024-04-19
+ 2024-04-21
daily
https://molssi.github.io/QCManyBody/api/
- 2024-04-19
+ 2024-04-21
daily
https://molssi.github.io/QCManyBody/explanation/
- 2024-04-19
+ 2024-04-21
daily
https://molssi.github.io/QCManyBody/how-to-guides/
- 2024-04-19
+ 2024-04-21
daily
https://molssi.github.io/QCManyBody/reference/
- 2024-04-19
+ 2024-04-21
daily
https://molssi.github.io/QCManyBody/tutorials/
- 2024-04-19
+ 2024-04-21
daily
\ No newline at end of file
diff --git a/sitemap.xml.gz b/sitemap.xml.gz
index 1454439..5a7846c 100644
Binary files a/sitemap.xml.gz and b/sitemap.xml.gz differ