Skip to content

Commit

Permalink
1-body fixes (#7)
Browse files Browse the repository at this point in the history 
* fix 1-body

* another note
  • Loading branch information
@loriab
loriab authored Apr 1, 2024
1 parent ff6fd99 commit b8f0645
Show file tree
Hide file tree
Showing 5 changed files with 128 additions and 93 deletions.
46 changes: 29 additions & 17 deletions qcmanybody/manybody.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,13 @@ def _assemble_nbody_components(
property_shape = find_shape(component_results[first_key])

# Final dictionaries
# * {bsse_type}_by_level is filled by sum_cluster_data to contain for NOCP
# & CP the summed total energies (or other property) of each nb-body. That is:
# * NOCP: {1: 1b@1b, 2: 2b@2b, ..., max_nbody: max_nbody-b@max_nbody-b} and
# * CP: {1: 1b@nfr-b, 2: 2b@nfr-b, ..., max_nbody: max_nbody-b@nfr-b}.
# VMFC bookkeeping is different. For key 1 it contains the summed 1b total energies.
# But for higher keys, it contains each nb-body (non-additive) contribution to the energy.
# * VMFC: {1: 1b@1b, 2: 2b contrib, ..., max_nbody: max_nbody-b contrib}
cp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
nocp_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
vmfc_by_level = {n: shaped_zero(property_shape) for n in range(1, nbodies[-1] + 1)}
Expand Down Expand Up @@ -312,8 +319,8 @@ def _analyze(
property_shape = find_shape(property_results[first_key])

property_result = shaped_zero(property_shape)
property_body_dict = {b.value: {} for b in self.bsse_type}
property_body_contribution = {b.value: {} for b in self.bsse_type}
property_body_dict = {bt.value: {} for bt in self.bsse_type}
property_body_contribution = {bt.value: {} for bt in self.bsse_type}

# results per model chemistry
mc_results = {}
Expand All @@ -339,22 +346,22 @@ def _analyze(
mc_results[mc_label] = nb_component_results

for n in nbody_list[::-1]:
property_bsse_dict = {b.value: shaped_zero(property_shape) for b in self.bsse_type}
property_bsse_dict = {bt.value: shaped_zero(property_shape) for bt in self.bsse_type}

for m in range(n - 1, n + 1):
if m == 0:
continue

# Subtract the (n-1)-body contribution from the n-body contribution to get the n-body effect
sign = (-1) ** (1 - m // n)
for b in self.bsse_type:
property_bsse_dict[b.value] += (
sign * mc_results[mc_label][f"{b.value}_{property_label}_body_dict"][m]
for bt in self.bsse_type:
property_bsse_dict[bt.value] += (
sign * mc_results[mc_label][f"{bt.value}_{property_label}_body_dict"][m]
)

property_result += property_bsse_dict[self.return_bsse_type]
for b in self.bsse_type:
property_body_contribution[b.value][n] = property_bsse_dict[b.value]
for bt in self.bsse_type:
property_body_contribution[bt.value][n] = property_bsse_dict[bt.value]

if self.has_supersystem:
# Get the MC label for supersystem tasks
Expand All @@ -372,13 +379,13 @@ def _analyze(
supersystem_result = property_results[ss_label]
property_result += supersystem_result - ss_component_results[f"{property_label}_body_dict"][self.max_nbody]

for b in self.bsse_type:
property_body_contribution[b][self.nfragments] = (
for bt in self.bsse_type:
property_body_contribution[bt][self.nfragments] = (
supersystem_result - ss_component_results[f"{property_label}_body_dict"][self.max_nbody]
)

for b in self.bsse_type:
bstr = b.value
for bt in self.bsse_type:
bstr = bt.value
for n in property_body_contribution[bstr]:
property_body_dict[bstr][n] = sum(
[
Expand Down Expand Up @@ -418,6 +425,11 @@ def analyze(

# Remove any missing data
component_results_inv = {k: v for k, v in component_results_inv.items() if v}
if not component_results_inv:
# Note B: Rarely, "no results" is expected, like for CP-only,
# rtd=False, and max_nbody=1. We'll add a dummy entry so
# processing can continue.
component_results_inv["energy"] = {'["dummy", [1000], [1000]]': 0.0}

# Actually analyze
all_results = {}
Expand All @@ -437,23 +449,23 @@ def analyze(

is_embedded = bool(self.embedding_charges)

for b in self.bsse_type:
for bt in self.bsse_type:
print_nbody_energy(
all_results["energy_body_dict"][b],
f"{b.upper()}-corrected multilevel many-body expansion",
all_results["energy_body_dict"][bt],
f"{bt.upper()}-corrected multilevel many-body expansion",
self.nfragments,
is_embedded,
)

if not self.has_supersystem: # skipped levels?
nbody_dict.update(
collect_vars(b.upper(), all_results["energy_body_dict"][b], self.max_nbody, is_embedded, self.supersystem_ie_only)
collect_vars(bt.upper(), all_results["energy_body_dict"][bt], self.max_nbody, is_embedded, self.supersystem_ie_only)
)

all_results["results"] = nbody_dict

# Make dictionary with "1cp", "2cp", etc
ebd = all_results["energy_body_dict"]
all_results["energy_body_dict"] = {str(k) + b: v for b in ebd for k, v in ebd[b].items()}
all_results["energy_body_dict"] = {str(k) + bt: v for bt in ebd for k, v in ebd[bt].items()}

return all_results
18 changes: 9 additions & 9 deletions qcmanybody/models/manybody_v1.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,8 +130,8 @@ class ManyBodyKeywords(ProtoModel):
"{max_nbody}-BODY`` will be available depending on ``return_total_data``; and ``{max_nbody}-BODY "
"CONTRIBUTION TO {driver}`` won't be available (except for dimers). This keyword produces no savings for a "
"two-fragment molecule. But for the interaction energy of a three-fragment molecule, for example, 2-body "
"subsystems can be skipped with ``supersystem_ie_only=True`` Do not use with ``vmfc`` in ``bsse_type``; it "
"is not implemented as ``cp`` is equivalent."
"subsystems can be skipped with ``supersystem_ie_only=True`` Do not use with ``vmfc`` in ``bsse_type``"
"as it cannot produce savings."
)

# v2: @field_validator("bsse_type", mode="before")
Expand Down Expand Up @@ -342,13 +342,13 @@ class Config:
json_schema_extra=jse,
))

# CP-CORRECTED INTERATION ENERGY THROUGH {nb}-BODY
for nb in range(2, MAX_NBODY):
# CP-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY
for nb in range(1, MAX_NBODY):
mbprop[f"cp_corrected_interaction_{egh}_through_{nb}_body"] = (
Optional[typ],
Field(
None,
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} with cp treatment. Available when when cp in bsse_type & max_nbody>={nb}{availability_of_derivative}.",
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} with cp treatment. Available when when cp in bsse_type & max_nbody>={nb}{availability_of_derivative}. The 1-body quantity is zero by definition.",
json_schema_extra=jse,
))

Expand Down Expand Up @@ -393,12 +393,12 @@ class Config:
))

# NOCP-CORRECTED INTERATION ENERGY THROUGH {nb}-BODY
for nb in range(2, MAX_NBODY):
for nb in range(1, MAX_NBODY):
mbprop[f"nocp_corrected_interaction_{egh}_through_{nb}_body"] = (
Optional[typ],
Field(
None,
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} without cp treatment. Available when when nocp in bsse_type & max_nbody>={nb}{availability_of_derivative}.",
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} without cp treatment. Available when when nocp in bsse_type & max_nbody>={nb}{availability_of_derivative}. The 1-body quantity is zero by definition.",
json_schema_extra=jse,
))

Expand Down Expand Up @@ -446,12 +446,12 @@ class Config:
))

# VMFC-CORRECTED INTERATION ENERGY THROUGH {nb}-BODY
for nb in range(2, MAX_NBODY):
for nb in range(1, MAX_NBODY):
mbprop[f"vmfc_corrected_interaction_{egh}_through_{nb}_body"] = (
Optional[typ],
Field(
None,
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} w/ vmfc treatment. Available when when vmfc in bsse_type & max_nbody>={nb}{availability_of_derivative}.",
description=f"{nb}-body total data less 1-body total data for cumulative IE; inputs are total {plural} w/ vmfc treatment. Available when when vmfc in bsse_type & max_nbody>={nb}{availability_of_derivative}. The 1-body quantity is zero by definition.",
json_schema_extra=jse,
))

Expand Down
100 changes: 52 additions & 48 deletions qcmanybody/models/test_mbe_he4_multilevel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,7 @@ def mbe_data_multilevel_631g():
"NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.472000052247,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.472089645469,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.472068853166,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008648503357,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.008558910134,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008579702437,
Expand All @@ -177,36 +178,39 @@ def mbe_data_multilevel_631g():
"CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.471058574581,
"CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.471324608815,
"CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.471272244751,
"CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.009589981022,
"CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.009323946788,
"CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.009376310852,
"CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.009589981022,
"CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000266034234,
"CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000052364064,
"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.009589981022,
"CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.009323946788,
"CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.009376310852,
"CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.009589981022,
"CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000266034234,
"CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000052364064,
},
# 1,2: ccsd; 3,4: mp2, all 6-31G
"22": {
"NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555603,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.471764016410,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.471853609632,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.471834096023,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008884539193,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.008794945971,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008814459580,
"NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.008884539193,
"NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000089593222,
"NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000019513609,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555603,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.471764016410,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.471853609632,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.471834096023,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008884539193,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.008794945971,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008814459580,
"NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.008884539193,
"NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000089593222,
"NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000019513609,

"CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555603,
"CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.470705938773,
"CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.470971973006,
"CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.470913449084,
"CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.009942616831,
"CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.009676582597,
"CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.009735106519,
"CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.009942616831,
"CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000266034234,
"CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000058523922,
"CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.480648555603,
"CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.470705938773,
"CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.470971973006,
"CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.470913449084,
"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.009942616831,
"CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.009676582597,
"CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.009735106519,
"CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.009942616831,
"CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000266034234,
"CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000058523922,
},
}

Expand All @@ -216,6 +220,7 @@ def mbe_data_multilevel_631g():
"CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.522467757090013,
"CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.522702864080149,
"CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.522639870651439,
"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008200959993875045,
"CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.007965853003739198,
"CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008028846432448944,
Expand All @@ -227,6 +232,7 @@ def mbe_data_multilevel_631g():
"NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.522851206178828,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.523095269671348,
"NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.523038093664368,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.007817510905059777,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.0075734474125397355,
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.007630623419519367,
Expand All @@ -238,6 +244,7 @@ def mbe_data_multilevel_631g():
"VMFC-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.52244892169719,
"VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.52268452228489,
"VMFC-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.522621528856181,
"VMFC-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": 0.0,
"VMFC-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.00821979538669737,
"VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.007984194798996924,
"VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.00804718822770667,
Expand Down Expand Up @@ -323,20 +330,19 @@ def mbe_data_multilevel_631g():
},
"1b_nocp_rtd": {
"NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY",
"NOCP-CORRECTED INTERACTION ENERGY": "zero", #"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
"NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
},
"1b_nocp": {
"NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY",
"NOCP-CORRECTED INTERACTION ENERGY": "zero", #"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
"NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
},
"1b_cp_rtd": {
"CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY",
"CP-CORRECTED INTERACTION ENERGY": "zero", #"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
"CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
},
"1b_cp": {
"CP-CORRECTED INTERACTION ENERGY": "zero", #"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
"CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
},
# TODO table defines the general qcvar as 0 even if 1-body qcvar not available. continue?
}


Expand Down Expand Up @@ -498,28 +504,27 @@ def he_tetramer():
{"121": 4, # 4hi
"22": 4},
id="1b_nocp_rtd"),
## TODO fix 1b for rtd=F
## pytest.param(
## {"bsse_type": "nocp", "return_total_data": False, "max_nbody": 1},
## "NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
## [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("1-BODY" in k))],
## {"121": 10,
## "22": 99}, #
## id="1b_nocp"),
pytest.param(
{"bsse_type": "nocp", "return_total_data": False, "max_nbody": 1},
"NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
[k for k in he4_refs_conv if (k.startswith("NOCP-") and ("1-BODY" in k))],
{"121": 4, # 4hi
"22": 4},
id="1b_nocp"),
pytest.param(
{"bsse_type": "cp", "return_total_data": True, "max_nbody": 1},
"CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY",
[k for k in he4_refs_conv if (k.startswith("CP-") and ("1-BODY" in k))],
{"121": 4, # 4hi
"22": 4},
id="1b_cp_rtd"),
## pytest.param(
## {"bsse_type": "cp", "return_total_data": False, "max_nbody": 1},
## "CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
## [k for k in he4_refs_conv if (k.startswith("CP-") and ("1-BODY" in k) and "TOTAL ENERGY" not in k)],
## {"121": 4,
## "22": 99}, #
## id="1b_cp"),
pytest.param(
{"bsse_type": "cp", "return_total_data": False, "max_nbody": 1},
"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY",
[k for k in he4_refs_conv if (k.startswith("CP-") and ("1-BODY" in k) and "TOTAL ENERGY" not in k)],
{"121": 0,
"22": 0},
id="1b_cp"),
])
def test_nbody_he4_multi(levels, mbe_keywords, anskey, bodykeys, calcinfo_nmbe, he_tetramer, request, mbe_data_multilevel_631g):
_inner = request.node.name.split("[")[1].split("]")[0]
Expand Down Expand Up @@ -576,8 +581,7 @@ def test_nbody_he4_multi(levels, mbe_keywords, anskey, bodykeys, calcinfo_nmbe,
for qcv in sumdict["4b_all"]:
skp = skprop(qcv)
if qcv in sumdict[kwdsln]:
refkey = sumdict[kwdsln][qcv]
ref = 0.0 if refkey == "zero" else refs[refkey]
ref = refs[sumdict[kwdsln][qcv]]
assert compare_values(ref, ret.extras["qcvars"]["nbody"][qcv], atol=atol, label=f"[c] qcvars {qcv}")
assert compare_values(ref, getattr(ret.properties, skp), atol=atol, label=f"[d] skprop {skp}")
else:
Expand Down
Loading

0 comments on commit b8f0645

Please sign in to comment.