Skip to content

Instantly share code, notes, and snippets.

@loriab

loriab/mbe_patch

Created February 12, 2024 21:53
Show Gist options
  • Save loriab/efce9ed6e298da622132f496e10f98dc to your computer and use it in GitHub Desktop.
Save loriab/efce9ed6e298da622132f496e10f98dc to your computer and use it in GitHub Desktop.
Download ZIP
tweaks to psi4 nbody
Raw
mbe_patch
diff --git a/psi4/driver/driver_nbody.py b/psi4/driver/driver_nbody.py
index b1c1b2eb5..a21767068 100644
--- a/psi4/driver/driver_nbody.py
+++ b/psi4/driver/driver_nbody.py
@@ -399,6 +399,7 @@ def build_nbody_compute_list(
nbodies: List[Union[int, Literal["supersystem"]]],
nfragments: int,
return_total_data: bool,
+ supersystem_ie_only: bool,
verbose: int = 1,
) -> Dict[str, Dict[int, Set[FragBasIndex]]]:
"""Generates lists of N-Body computations needed for requested BSSE treatments.
@@ -417,6 +418,8 @@ def build_nbody_compute_list(
Formerly max_frag
return_total_data
Whether the total data (True; energy/gradient/Hessian) of the molecular system has been requested, as opposed to interaction data (False).
+ supersystem_ie_only
+ ????
verbose
Control volume of printing.
@@ -465,20 +468,41 @@ def build_nbody_compute_list(
# Everything is in full n-mer basis
basis_tuple = tuple(fragment_range)
- for nb in nbodies:
- if nb > 1:
- for sublevel in range(1, nb + 1):
- for x in itertools.combinations(fragment_range, sublevel):
- # below was `nbodies`, which would never hit. present is closest to pre-DDD. purpose unclear to me.
- # if self.max_nbody == 1: break
- cp_compute_list[nb].add((x, basis_tuple))
-
- if 'nocp' in bsse_type or return_total_data:
+ if supersystem_ie_only:
+ nb = nfragments
+ for sublevel in [1, nfragments]:
+ for x in itertools.combinations(fragment_range, sublevel):
+ cp_compute_list[nb].add((x, basis_tuple))
+ print(f"ADDING CP: {nb=} {sublevel=} {x=} {basis_tuple=}")
+ else:
+ for nb in nbodies:
+ if nb > 1:
+ for sublevel in range(1, nb + 1):
+ for x in itertools.combinations(fragment_range, sublevel):
+ # below was `nbodies`, which would never hit. present is closest to pre-DDD. purpose unclear to me.
+ # if self.max_nbody == 1: break
+ cp_compute_list[nb].add((x, basis_tuple))
+ print(f"adding CP: {nb=} {sublevel=} {x=} {basis_tuple=}")
+
+ if return_total_data and 1 in nbodies:
+ # Monomers in monomer basis
+ nb = 1
+ sublevel = 1
+ for x in itertools.combinations(fragment_range, sublevel):
+ nocp_compute_list[nb].add((x, x))
+
+ if 'nocp' in bsse_type:
# Everything in monomer basis
- for nb in nbodies:
- for sublevel in range(1, nb + 1):
+ if supersystem_ie_only:
+ nb = nfragments
+ for sublevel in [1, nb]:
for x in itertools.combinations(fragment_range, sublevel):
nocp_compute_list[nb].add((x, x))
+ else:
+ for nb in nbodies:
+ for sublevel in range(1, nb + 1):
+ for x in itertools.combinations(fragment_range, sublevel):
+ nocp_compute_list[nb].add((x, x))
if 'vmfc' in bsse_type:
# Like a CP for all combinations of pairs or greater
@@ -487,9 +511,8 @@ def build_nbody_compute_list(
basis_tuple = tuple(cp_combos)
for interior_nbody in range(1, nb + 1):
for x in itertools.combinations(cp_combos, interior_nbody):
- combo_tuple = (x, basis_tuple)
- vmfc_compute_list[nb].add(combo_tuple)
- vmfc_level_list[len(basis_tuple)].add(combo_tuple)
+ vmfc_compute_list[nb].add((x, basis_tuple))
+ vmfc_level_list[len(basis_tuple)].add((x, basis_tuple))
# Build a comprehensive compute range
# * do not use list length to count number of {nb}-body computations
@@ -545,6 +568,8 @@ def assemble_nbody_components(
Formerly max_frag
return_total_data : bool
See class field. Whether the total data (e/g/H) of the molecular system has been requested, as opposed to interaction data.
+ supersystem_ie_only: bool
+ ???
max_nbody : int
See class field. Maximum number of bodies to include in the many-body treatment."
embedding_charges : bool
@@ -581,7 +606,7 @@ def assemble_nbody_components(
# regenerate per-bsse required calcs list
compute_dict = build_nbody_compute_list(
- metadata['bsse_type'], nbodies, metadata['nfragments'], metadata["return_total_data"], verbose=0
+ metadata['bsse_type'], nbodies, metadata['nfragments'], metadata["return_total_data"], metadata["supersystem_ie_only"], verbose=0
)
# Build size and slices dictionaries
@@ -589,7 +614,11 @@ def assemble_nbody_components(
fragment_slice_dict = {}
iat = 0
for ifr in range(1, metadata["nfragments"] + 1):
- nat = metadata["molecule"].extract_subsets(ifr).natom()
+ try:
+ nat = metadata["molecule"].extract_subsets(ifr).natom()
+ except AttributeError:
+ # qcelemental.models.Molecule
+ nat = len(metadata["molecule"].fragments[ifr - 1])
fragment_size_dict[ifr] = nat
fragment_slice_dict[ifr] = slice(iat, iat + nat)
iat += nat
@@ -727,11 +756,21 @@ def assemble_nbody_components(
nbody_dict["CP-CORRECTED TOTAL ENERGY"] = cp_body_dict[metadata['max_nbody']]
nbody_dict["CP-CORRECTED INTERACTION ENERGY"] = cp_body_dict[metadata['max_nbody']] - cp_body_dict[1]
- for nb in nbodies[1:]:
- nbody_dict[f"CP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb] - cp_body_dict[1]
- nbody_dict[f"CP-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = cp_body_dict[nb] - cp_body_dict[nb-1]
- for nb in nbodies:
- nbody_dict[f"CP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb]
+ if metadata["supersystem_ie_only"]:
+ nfr = metadata["nfragments"]
+ for nb in [nfr]:
+ nbody_dict[f"CP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb] - cp_body_dict[1]
+ # define for 2mer nbody_dict[f"CP-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = cp_body_dict[nb] - cp_body_dict[nb-1]
+ if monomer_sum != 0.0:
+ for nb in [1, nfr]:
+ nbody_dict[f"CP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb]
+ else:
+ for nb in range(2, max(nbodies) + 1):
+ nbody_dict[f"CP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb] - cp_body_dict[1]
+ nbody_dict[f"CP-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = cp_body_dict[nb] - cp_body_dict[nb-1]
+ if monomer_sum != 0.0:
+ for nb in nbodies:
+ nbody_dict[f"CP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = cp_body_dict[nb]
# Compute nocp
if 'nocp' in metadata['bsse_type']:
@@ -755,11 +794,18 @@ def assemble_nbody_components(
nbody_dict['NOCP-CORRECTED TOTAL ENERGY'] = nocp_body_dict[metadata['max_nbody']]
nbody_dict['NOCP-CORRECTED INTERACTION ENERGY'] = nocp_body_dict[metadata['max_nbody']] - nocp_body_dict[1]
- for nb in nbodies[1:]:
- nbody_dict[f"NOCP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb] - nocp_body_dict[1]
- nbody_dict[f"NOCP-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = nocp_body_dict[nb] - nocp_body_dict[nb-1]
- for nb in nbodies:
- nbody_dict[f"NOCP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb]
+ if metadata["supersystem_ie_only"]:
+ nfr = metadata["nfragments"]
+ for nb in [nfr]:
+ nbody_dict[f"NOCP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb] - nocp_body_dict[1]
+ for nb in [1, nfr]:
+ nbody_dict[f"NOCP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb]
+ else:
+ for nb in nbodies[1:]:
+ nbody_dict[f"NOCP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb] - nocp_body_dict[1]
+ nbody_dict[f"NOCP-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = nocp_body_dict[nb] - nocp_body_dict[nb-1]
+ for nb in nbodies:
+ nbody_dict[f"NOCP-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = nocp_body_dict[nb]
# Compute vmfc
if 'vmfc' in metadata['bsse_type']:
@@ -862,6 +908,9 @@ class ManyBodyComputer(BaseComputer):
embedding_charges: Dict[int, List[float]] = Field({}, description="Atom-centered point charges to be used on molecule fragments whose basis sets are not included in the computation. Keys: 1-based index of fragment. Values: list of atom charges for that fragment.")
return_total_data: Optional[bool] = Field(None, validate_default=True, description="When True, returns the total data (energy/gradient/Hessian) of the system, otherwise returns interaction data. Default is False for energies, True for gradients and Hessians. Note that the calculation of total counterpoise corrected energies implies the calculation of the energies of monomers in the monomer basis, hence specifying ``return_total_data = True`` may carry out more computations than ``return_total_data = False``.")
+
+ supersystem_ie_only: Optional[bool] = Field(False) # TODO check max_nbody=nfr
+
quiet: bool = Field(False, description="Whether to print/log formatted n-body energy analysis. Presently used by multi to suppress output. Candidate for removal from class once in-class/out-of-class functions sorted.")
task_list: Dict[str, MBETaskComputers] = {}
@@ -892,10 +941,13 @@ class ManyBodyComputer(BaseComputer):
@field_validator("max_nbody")
@classmethod
def set_max_nbody(cls, v: int, info: FieldValidationInfo) -> int:
+ nfragments = info.data["nfragments"]
if v == -1:
- return info.data["nfragments"]
+ return nfragments
+ elif v < 0 or v > nfragments:
+ raise ValueError(f"max_nbody={v} should be between 1 and {nfragments}.")
else:
- return min(v, info.data["nfragments"])
+ return min(v, nfragments)
@field_validator("embedding_charges")
@classmethod
@@ -974,10 +1026,12 @@ class ManyBodyComputer(BaseComputer):
count += 1
compute_dict = build_nbody_compute_list(
- ["nocp"], list(range(1, self.max_nbody + 1)), self.nfragments, self.return_total_data
- )
+ ["nocp"], list(range(1, self.max_nbody + 1)),
+ self.nfragments, self.return_total_data, self.supersystem_ie_only)
else:
- compute_dict = build_nbody_compute_list(self.bsse_type, nbodies, self.nfragments, self.return_total_data)
+ compute_dict = build_nbody_compute_list(
+ self.bsse_type, nbodies,
+ self.nfragments, self.return_total_data, self.supersystem_ie_only)
def labeler(item) -> str:
mc_level_lbl = mc_level_idx + 1
@@ -1300,6 +1354,7 @@ class ManyBodyComputer(BaseComputer):
"bsse_type": self.bsse_type,
"nfragments": self.nfragments,
"return_total_data": self.return_total_data,
+ "supersystem_ie_only": self.supersystem_ie_only,
"molecule": self.molecule,
"embedding_charges": bool(self.embedding_charges),
"max_nbody": self.max_nbody,
diff --git a/psi4/driver/driver_nbody_multilevel.py b/psi4/driver/driver_nbody_multilevel.py
index 219541ecd..5bfeb8e44 100644
--- a/psi4/driver/driver_nbody_multilevel.py
+++ b/psi4/driver/driver_nbody_multilevel.py
@@ -43,8 +43,13 @@ def prepare_results(self, client: Optional["qcportal.FractalClient"] = None) ->
from psi4.driver.driver_nbody import _print_nbody_energy
ptype = self.driver.name
- natoms = self.molecule.natom()
+ try:
+ natoms = self.molecule.natom()
+ except AttributeError:
+ # qcelemental.models.Molecule
+ natoms = len(self.molecule.symbols)
supersystem = {k: v for k, v in self.task_list.items() if k.startswith('supersystem')}
+ is_supersystem = bool(supersystem)
# Initialize with zeros
energy_result, gradient_result, hessian_result = 0, None, None
@@ -54,6 +59,7 @@ def prepare_results(self, client: Optional["qcportal.FractalClient"] = None) ->
gradient_result = np.zeros((natoms, 3))
if ptype == 'hessian':
hessian_result = np.zeros((natoms * 3, natoms * 3))
+ hold_max_nbody = self.max_nbody
# Get numerical label (index) for supersystem tasks
sup_level = 0
@@ -123,6 +129,7 @@ def prepare_results(self, client: Optional["qcportal.FractalClient"] = None) ->
gradient_result += np.array(supersystem_result.return_result).reshape((-1, 3)) - components[f'{ptype}_body_dict'][self.max_nbody]
+ nbody_dict = {}
for b in self.bsse_type:
for n in energy_body_contribution[b]:
energy_body_dict[b][n] = sum(
@@ -132,6 +139,8 @@ def prepare_results(self, client: Optional["qcportal.FractalClient"] = None) ->
for b in self.bsse_type:
_print_nbody_energy(energy_body_dict[b], f"{b.upper()}-corrected multilevel many-body expansion",
self.nfragments, is_embedded)
+ if not is_supersystem: # skipped levels?
+ nbody_dict.update(seteqcvars(b.upper(), energy_body_dict[b], hold_max_nbody, is_embedded))
if not self.return_total_data:
# Remove monomer cotribution for interaction data
@@ -147,5 +156,30 @@ def prepare_results(self, client: Optional["qcportal.FractalClient"] = None) ->
"ret_energy": energy_result,
"ret_ptype": locals()[ptype + '_result'],
"energy_body_dict": energy_body_dict,
+ "nbody": nbody_dict, # ptype == "energy" only?
}
return nbody_results
+
+def seteqcvars(bsse, body_dict, max_nbody: int, embedding: bool = False):
+ previous_e = body_dict[1]
+ tot_e = (previous_e != 0.0)
+ nbody_range = list(body_dict)
+ nbody_range.sort()
+ qcvars = {}
+ if embedding:
+ return qcvars
+
+ if tot_e:
+ qcvars[f"{bsse}-CORRECTED TOTAL ENERGY"] = body_dict[max_nbody]
+ qcvars[f"{bsse}-CORRECTED INTERACTION ENERGY"] = body_dict[max_nbody] - body_dict[1]
+
+ for nb in range(2, max(nbody_range) + 1):
+ print(f"{nb=} {body_dict[nb]=} {body_dict=}")
+ qcvars[f"{bsse}-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = body_dict[nb] - body_dict[1]
+ qcvars[f"{bsse}-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = body_dict[nb] - body_dict[nb-1]
+ if tot_e:
+ for nb in nbody_range:
+ qcvars[f"{bsse}-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = body_dict[nb]
+
+ return qcvars
+
diff --git a/psi4/driver/task_planner.py b/psi4/driver/task_planner.py
index e92e98c05..103389c51 100644
--- a/psi4/driver/task_planner.py
+++ b/psi4/driver/task_planner.py
@@ -159,6 +159,7 @@ def task_planner(driver: DriverEnum, method: str, molecule: core.Molecule, **kwa
plan = ManyBodyComputer(**packet, **kwargs)
original_molecule = packet.pop("molecule")
+ user_max_nbody = kwargs.get("max_nbody", None)
# Add tasks for every nbody level requested
if levels is None:
@@ -169,13 +170,6 @@ def task_planner(driver: DriverEnum, method: str, molecule: core.Molecule, **kwa
# We define cp as being a correction to only interaction energies
# If only doing cp, we need to ignore any user-specified 1st (monomer) level
- if 'cp' in kwargs.get("bsse_type", None) and 'nocp' not in kwargs.get("bsse_type", None):
- if 1 in levels.keys():
- removed_level = levels.pop(1)
- logger.info("NOTE: User specified exclusively 'cp' correction, but provided level 1 details")
- logger.info(f"NOTE: Removing level {removed_level}")
- logger.info("NOTE: For total energies, add 'nocp' to bsse_list")
-
# Organize nbody calculations into modelchem levels
# * expand keys of `levels` into full lists of nbodies covered. save to plan, resetting max_nbody accordingly
@@ -194,7 +188,10 @@ def task_planner(driver: DriverEnum, method: str, molecule: core.Molecule, **kwa
nbodies_per_mc_level.append(nbodies)
prev_body += 1
- plan.max_nbody = max(nb for nb in levels if nb != "supersystem")
+ new_max_nbody = max(nb for nb in levels if nb != "supersystem")
+ if user_max_nbody and (new_max_nbody != plan.max_nbody):
+ raise ValueError(f"levels contradicts user max_nbody={user_max_nbody}.")
+ plan.max_nbody = new_max_nbody
plan.nbodies_per_mc_level = nbodies_per_mc_level
for mc_level_idx, mtd in enumerate(levels.values()):
diff --git a/pytest.ini b/pytest.ini
index cdc4f32b4..0e3e8ac60 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -39,7 +39,7 @@ markers =
cubeprop
d2ints: "tests that fail when 2nd derivative ERI missing from Libint2, either because dertype 2 demanded or tested too tightly for 3-pt FD."
dct
- # ddd: "findif and cbs and nbody" (not defined, so use equiv. expansion)
+ # ddd: "-m 'findif or cbs or nbody'" (not defined, so use equiv. expansion)
df: "tests that employ density-fitting"
dfccd
dfccd-grad
diff --git a/tests/nbody-cp-gradient/input.dat b/tests/nbody-cp-gradient/input.dat
index a4ffd62af..beb8d527d 100644
--- a/tests/nbody-cp-gradient/input.dat
+++ b/tests/nbody-cp-gradient/input.dat
@@ -176,7 +176,7 @@ for result in [ #TEST
]: #TEST
compare_values(ref_2cp_grad, result, 8, 'CP-Corrected Gradient') #TEST
-compare(12, wfn.variable("NBODY NUMBER"), "rtd=T nbody number") #TEST
+compare(9, wfn.variable("NBODY NUMBER"), "rtd=T nbody number") #TEST
# compute nbody again with return_total_data=False
diff --git a/tests/nbody-he-4b/input.dat b/tests/nbody-he-4b/input.dat
index 4a0b99d27..6f372a2e3 100644
--- a/tests/nbody-he-4b/input.dat
+++ b/tests/nbody-he-4b/input.dat
@@ -30,7 +30,7 @@ else:
e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type=["nocp", "cp", "vmfc"], return_wfn=True)
for k, v in sorted(psi4.core.variables().items()):
- print("QCVariable:", k, v)
+ print("[1] All,4b QCVariable:", k, v)
he4_refs = {
"CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530751941948,
@@ -82,6 +82,7 @@ for qcv, ref in {
}.items():
compare_values(ref, variable(qcv), 8, "[1] core " + qcv)
compare_values(ref, wfn.variable(qcv), 8, "[1] wfn " + qcv)
+compare(65, variable("NBODY NUMBER"), "[1] nmbe")
clean_variables()
@@ -94,6 +95,9 @@ if distributed:
else:
e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_wfn=True, max_nbody=3)
+for k, v in sorted(psi4.core.variables().items()):
+ print("[2] NOCP,3b QCVariable:", k, v)
+
for qcv, ref in he4_refs.items(): #TEST
if "4-BODY" not in qcv and qcv.startswith("NOCP-"): #TEST
compare_values(ref, variable(qcv), 8, "[2] core " + qcv) #TEST
@@ -106,6 +110,7 @@ for qcv, ref in {
}.items():
compare_values(ref, variable(qcv), 8, "[2] core " + qcv)
compare_values(ref, wfn.variable(qcv), 8, "[2] wfn " + qcv)
+compare(14, variable("NBODY NUMBER"), "[2] nmbe")
clean_variables()
@@ -118,6 +123,9 @@ if distributed:
else:
e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_wfn=True, max_nbody=3, return_total_data=True)
+for k, v in sorted(psi4.core.variables().items()):
+ print("[3] CP,3b,rtd=T QCVariable:", k, v)
+
for qcv, ref in he4_refs.items(): #TEST
if "4-BODY" not in qcv and qcv.startswith("CP-"): #TEST
compare_values(ref, variable(qcv), 8, "[3] core " + qcv) #TEST
@@ -130,6 +138,8 @@ for qcv, ref in {
}.items():
compare_values(ref, variable(qcv), 8, "[3] core " + qcv) #TEST
compare_values(ref, wfn.variable(qcv), 8, "[3] wfn " + qcv) #TEST
+#compare(28, variable("NBODY NUMBER"), "[3] nmbe") # was
+compare(18, variable("NBODY NUMBER"), "[3] nmbe") # bugfix
clean_variables()
@@ -142,6 +152,9 @@ if distributed:
else:
e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_wfn=True, max_nbody=3, return_total_data=False)
+for k, v in sorted(psi4.core.variables().items()):
+ print("[4] CP,3b,rtd=F QCVariable:", k, v)
+
for qcv, ref in he4_refs.items(): #TEST
if "4-BODY" not in qcv and qcv.startswith("CP-") and "TOTAL ENERGY" not in qcv: #TEST
compare_values(ref, variable(qcv), 8, "[4] core " + qcv) #TEST
@@ -153,17 +166,21 @@ for qcv, ref in {
}.items():
compare_values(ref, variable(qcv), 8, "[4] core " + qcv) #TEST
compare_values(ref, wfn.variable(qcv), 8, "[4] wfn " + qcv) #TEST
+compare(14, variable("NBODY NUMBER"), "[4] nmbe")
clean_variables()
###
if distributed:
- plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_plan=True, max_nbody=3)
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_plan=True, max_nbody=3, return_total_data=True)
plan.compute(client)
snowflake.await_results()
e, wfn = plan.get_psi_results(client, return_wfn=True)
else:
- e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_wfn=True, max_nbody=3)
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_wfn=True, max_nbody=3, return_total_data=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[5] VMFC,3b,rtd=T QCVariable:", k, v)
for qcv, ref in he4_refs.items():
if "4-BODY" not in qcv and qcv.startswith("VMFC-"):
@@ -173,10 +190,177 @@ for qcv, ref in he4_refs.items():
for qcv, ref in {
"VMFC-CORRECTED TOTAL ENERGY": he4_refs["VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY"],
"VMFC-CORRECTED INTERACTION ENERGY": he4_refs["VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
- "CURRENT ENERGY": he4_refs["VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+ "CURRENT ENERGY": he4_refs["VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY"],
}.items():
compare_values(ref, variable(qcv), 8, "[5] core " + qcv) #TEST
compare_values(ref, wfn.variable(qcv), 8, "[5] wfn " + qcv) #TEST
+compare(50, variable("NBODY NUMBER"), "[5] nmbe")
if distributed:
snowflake.stop()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_plan=True, max_nbody=3, return_total_data=False)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="vmfc", return_wfn=True, max_nbody=3, return_total_data=False)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[6] VMFC,3b,rtd=F QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items():
+ if "4-BODY" not in qcv and qcv.startswith("VMFC-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[6] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[6] wfn " + qcv) #TEST
+
+for qcv, ref in {
+ "VMFC-CORRECTED TOTAL ENERGY": he4_refs["VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY"],
+ "VMFC-CORRECTED INTERACTION ENERGY": he4_refs["VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+ "CURRENT ENERGY": he4_refs["VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+}.items():
+ compare_values(ref, variable(qcv), 8, "[6] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[6] wfn " + qcv) #TEST
+compare(50, variable("NBODY NUMBER"), "[6] nmbe")
+
+if distributed:
+ snowflake.stop()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_plan=True, max_nbody=3, return_total_data=True)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_wfn=True, max_nbody=3, return_total_data=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[7] NOCP,3b,rtd=T QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if "4-BODY" not in qcv and qcv.startswith("NOCP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[7] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[7] wfn " + qcv) #TEST
+
+for qcv, ref in {
+ "NOCP-CORRECTED TOTAL ENERGY": he4_refs["NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY"],
+ "NOCP-CORRECTED INTERACTION ENERGY": he4_refs["NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+ "CURRENT ENERGY": he4_refs["NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY"],
+}.items():
+ compare_values(ref, variable(qcv), 8, "[7] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[7] wfn " + qcv) #TEST
+compare(14, variable("NBODY NUMBER"), "[7] nmbe")
+
+clean_variables()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_plan=True, max_nbody=3, return_total_data=False)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_wfn=True, max_nbody=3, return_total_data=False)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[8] NOCP,3b,rtd=F QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if "4-BODY" not in qcv and qcv.startswith("NOCP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[8] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[8] wfn " + qcv) #TEST
+
+for qcv, ref in {
+ "NOCP-CORRECTED INTERACTION ENERGY": he4_refs["NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+ "CURRENT ENERGY": he4_refs["NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY"],
+}.items():
+ compare_values(ref, variable(qcv), 8, "[8] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[8] wfn " + qcv) #TEST
+compare(14, variable("NBODY NUMBER"), "[8] nmbe")
+
+clean_variables()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_total_data=False, supersystem_ie_only=True, return_plan=True)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_total_data=False, supersystem_ie_only=True, return_wfn=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[9] NOCP,4b,rtd=F,sio=T QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if "THROUGH 4-BODY" in qcv and qcv.startswith("NOCP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[9] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[9] wfn " + qcv) #TEST
+compare(5, variable("NBODY NUMBER"), "[9] nmbe")
+
+clean_variables()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_total_data=True, supersystem_ie_only=True, return_plan=True)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="nocp", return_total_data=True, supersystem_ie_only=True, return_wfn=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[10] NOCP,4b,rtd=T,sio=T QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if ("THROUGH 4-BODY" in qcv or "THROUGH 1-BODY" in qcv) and qcv.startswith("NOCP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[10] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[10] wfn " + qcv) #TEST
+compare(5, variable("NBODY NUMBER"), "[10] nmbe")
+
+clean_variables()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_total_data=False, supersystem_ie_only=True, return_plan=True)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_total_data=False, supersystem_ie_only=True, return_wfn=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[11] CP,4b,rtd=F,sio=T QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if "INTERACTION ENERGY THROUGH 4-BODY" in qcv and qcv.startswith("CP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[11] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[11] wfn " + qcv) #TEST
+compare(5, variable("NBODY NUMBER"), "[11] nmbe")
+
+clean_variables()
+
+###
+if distributed:
+ plan = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_total_data=True, supersystem_ie_only=True, return_plan=True)
+ plan.compute(client)
+ snowflake.await_results()
+ e, wfn = plan.get_psi_results(client, return_wfn=True)
+else:
+ e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, bsse_type="cp", return_total_data=True, supersystem_ie_only=True, return_wfn=True)
+
+for k, v in sorted(psi4.core.variables().items()):
+ print("[12] CP,4b,rtd=T,sio=T QCVariable:", k, v)
+
+for qcv, ref in he4_refs.items(): #TEST
+ if ("THROUGH 4-BODY" in qcv or "THROUGH 1-BODY" in qcv) and qcv.startswith("CP-"): #TEST
+ compare_values(ref, variable(qcv), 8, "[12] core " + qcv) #TEST
+ compare_values(ref, wfn.variable(qcv), 8, "[12] wfn " + qcv) #TEST
+compare(9, variable("NBODY NUMBER"), "[12] nmbe")
+
+clean_variables()
+
+
diff --git a/tests/nbody-multi-level-2/input.dat b/tests/nbody-multi-level-2/input.dat
index ead3b3722..274c309b0 100644
--- a/tests/nbody-multi-level-2/input.dat
+++ b/tests/nbody-multi-level-2/input.dat
@@ -19,9 +19,6 @@ He 2.0 2.0 2.0
vars_1 = ["NOCP-CORRECTED TOTAL ENERGY", "CP-CORRECTED TOTAL ENERGY","VMFC-CORRECTED TOTAL ENERGY",
"NOCP-CORRECTED INTERACTION ENERGY", "CP-CORRECTED INTERACTION ENERGY","VMFC-CORRECTED INTERACTION ENERGY"]
-vars_2 = ["4NOCP", "4CP", "4VMFC"]
-vars_3 = ["3NOCP", "3CP", "3VMFC"]
-
levels_1 = {4:'hf/6-31g'}
levels_2 = {1:'ccsd/6-31g', 4:'hf/6-31g'}
levels_3 = {1:'ccsd/6-31g',2:'mp2/6-31g', 3:'hf/6-31g'}
@@ -38,6 +35,41 @@ ref2 = {"4NOCP" : -11.47436885006691,
"4CP" : -11.47390053367237,
"4VMFC" : -11.47396874657488}
+ref2_qcv = {
+ "CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555739,
+ "CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.473787566291,
+ "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.006860989448,
+ "CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.006860989448,
+ "CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.473902077326,
+ "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.006746478413,
+ "CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000114511036,
+ "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.473900533672,
+ "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.006748022067,
+ "CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000001543654,
+
+ "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555739,
+ "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.474336562701,
+ "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.006311993038,
+ "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.006311993038,
+ "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.474364757208,
+ "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.006283798531,
+ "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000028194507,
+ "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.474368850067,
+ "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.006279705672,
+ "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": -0.000004092859,
+
+ "VMFC-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555739,
+ "VMFC-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.473856767634,
+ "VMFC-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.006791788105,
+ "VMFC-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.006791788105,
+ "VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.473970290229,
+ "VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.006678265510,
+ "VMFC-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000113522595,
+ "VMFC-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.473968746575,
+ "VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.006679809164,
+ "VMFC-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000001543654,
+}
+
ref3 = {"3NOCP" : -11.47411506642999,
"3CP" : -11.47349141626227,
"3VMFC" : -11.47357018907939}
@@ -53,27 +85,44 @@ ref4 = {"4NOCP" : -11.47411915928889,
e, wfn = psi4.energy("", levels=levels_1, bsse_type=['cp','nocp','vmfc'], return_wfn=True, return_total_data=True)
for var in vars_1:
- energy = wfn.variable(var)
- compare_values(energy, ref1[var], 8, var + " 4-body HF-6-31g") #TEST
+ compare_values(wfn.variable(var), ref1[var], 8, "[1] " + var + " 4-body HF-6-31g") #TEST
+
+clean_variables()
+
+e, wfn = psi4.energy("", levels=levels_2, bsse_type=['nocp'], return_wfn=True, return_total_data=True)
+for var in ["4NOCP"]: #TEST
+ compare_values(wfn.variable(var), ref2[var], 8, "[2a] " + var) #TEST
+for var in ref2_qcv:
+ if var.startswith("NOCP-"):
+ compare_values(wfn.variable(var), ref2_qcv[var], 8, "[2a] " + var) #TEST
+
+clean_variables()
+
+e, wfn = psi4.energy("", levels=levels_2, bsse_type=['cp'], return_wfn=True, return_total_data=True)
+for var in ["4CP"]: #TEST
+ compare_values(wfn.variable(var), ref2[var], 8, "[2b] " + var) #TEST
+for var in ref2_qcv:
+ if var.startswith("CP-"):
+ compare_values(wfn.variable(var), ref2_qcv[var], 8, "[2b] " + var) #TEST
clean_variables()
-e, wfn = psi4.energy("", levels=levels_2, bsse_type=['cp','nocp','vmfc'], return_wfn=True, return_total_data=True)
-for var in vars_2:
- energy = wfn.variable(var)
- compare_values(energy, ref2[var], 8, var) #TEST
+e, wfn = psi4.energy("", levels=levels_2, bsse_type=['vmfc'], return_wfn=True, return_total_data=True)
+for var in ["4VMFC"]: #TEST
+ compare_values(wfn.variable(var), ref2[var], 8, "[2c] " + var) #TEST
+for var in ref2_qcv:
+ if var.startswith("VMFC-"):
+ compare_values(wfn.variable(var), ref2_qcv[var], 8, "[2c] " + var) #TEST
clean_variables()
e, wfn = psi4.energy("", levels=levels_3, bsse_type=['cp','nocp','vmfc'], return_wfn=True, return_total_data=True)
-for var in vars_3:
- energy = wfn.variable(var)
- compare_values(energy, ref3[var], 8, var) #TEST
+for var in ["3NOCP", "3CP", "3VMFC"]: #TEST
+ compare_values(wfn.variable(var), ref3[var], 8, "[3] " + var) #TEST
clean_variables()
e, wfn = psi4.energy("", levels=levels_4, bsse_type=['cp','nocp','vmfc'], return_wfn=True, return_total_data=True)
-for var in vars_2:
- energy = wfn.variable(var)
- compare_values(energy, ref4[var], 8, var) #TEST
+for var in ["4NOCP", "4CP", "4VMFC"]: #TEST
+ compare_values(wfn.variable(var), ref4[var], 8, "[4] " + var) #TEST
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment