diff --git a/qcmanybody/manybody.py b/qcmanybody/manybody.py index 1b6b3de..2a512da 100644 --- a/qcmanybody/manybody.py +++ b/qcmanybody/manybody.py @@ -35,10 +35,9 @@ def __init__( levels: Mapping[Union[int, Literal["supersystem"]], str], return_total_data: bool, supersystem_ie_only: bool, + embedding_charges: Mapping[int, Sequence[float]], ): - # TODO - self.embedding_charges = {} - + self.embedding_charges = embedding_charges self.molecule = molecule self.bsse_type = [BsseEnum(x) for x in bsse_type] self.return_total_data = return_total_data @@ -197,13 +196,13 @@ def iterate_molecules(self) -> Iterable[Tuple[str, str, Molecule]]: mol = self.molecule.get_fragment(real_atoms_0, ghost_atoms_0, orient=False, group_fragments=False) mol = mol.copy(update={"fix_com": True, "fix_orientation": True}) - # if self.embedding_charges: - # embedding_frags = list(set(range(1, self.nfragments + 1)) - set(pair[1])) - # charges = [] - # for frag in embedding_frags: - # positions = self.molecule.extract_subsets(frag).geometry().np.tolist() - # charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[frag])]) - # data['keywords']['function_kwargs'].update({'external_potentials': charges}) + if self.embedding_charges: + embedding_frags = list(set(range(1, self.nfragments + 1)) - set(basis_atoms)) + charges = [] + for ifr in embedding_frags: + positions = self.molecule.get_fragment(ifr-1).geometry.tolist() + charges.extend([[chg, i] for i, chg in zip(positions, self.embedding_charges[ifr])]) + mol.extras["embedding_charges"] = charges done_molecules.add(label) yield mc, label, mol diff --git a/qcmanybody/models/manybody_pydv1.py b/qcmanybody/models/manybody_pydv1.py index 5bb7e11..d43f14b 100644 --- a/qcmanybody/models/manybody_pydv1.py +++ b/qcmanybody/models/manybody_pydv1.py @@ -99,6 +99,8 @@ class ManyBodyKeywords(ProtoModel): {}, 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.", + # TODO embedding charges should sum to fragment charge, right? enforce? + # TODO embedding charges irrelevant to CP (basis sets always present)? json_schema_extra={ "shape": ["nfr", ""], }, diff --git a/qcmanybody/models/qcng_computer.py b/qcmanybody/models/qcng_computer.py index 0fe8c0e..4ce6a1c 100644 --- a/qcmanybody/models/qcng_computer.py +++ b/qcmanybody/models/qcng_computer.py @@ -212,10 +212,12 @@ def set_bsse_type(cls, v: Any) -> List[BsseEnum]: @classmethod # v2: def set_embedding_charges(cls, v: Any, info: FieldValidationInfo) -> Dict[int, List[float]]: def set_embedding_charges(cls, v, values): # -> Dict[int, List[float]]: + print(f"hit embedding_charges validator with {v}", end="") # v2: if len(v) != info.data["nfragments"]: - if len(v) != values["nfragments"]: - raise ValueError("embedding_charges dict should have entries for each 1-indexed fragment.") + if len(v) != len(values["molecule"].fragments): + raise ValueError(f"embedding_charges dict should have entries for each 1-indexed fragment ({len(values['molecule'].fragments)}).") + print(f" ... setting embedding_charges={v}") return v # v2: @field_validator("return_total_data") @@ -397,6 +399,7 @@ def from_qcschema_ben(cls, input_model: ManyBodyInput, build_tasks: bool = True) specifications[mtd] = {} specifications[mtd]["program"] = spec.pop("program") specifications[mtd]["specification"] = spec + specifications[mtd]["specification"]["driver"] = computer_model.driver # overrides atomic driver with mb driver specifications[mtd]["specification"].pop("schema_name", None) specifications[mtd]["specification"].pop("protocols", None) @@ -406,6 +409,7 @@ def from_qcschema_ben(cls, input_model: ManyBodyInput, build_tasks: bool = True) comp_levels, computer_model.return_total_data, computer_model.supersystem_ie_only, + computer_model.embedding_charges, ) if not build_tasks: @@ -416,6 +420,15 @@ def from_qcschema_ben(cls, input_model: ManyBodyInput, build_tasks: bool = True) for chem, label, imol in calculator_cls.iterate_molecules(): inp = AtomicInput(molecule=imol, **specifications[chem]["specification"]) + if imol.extras.get("embedding_charges"): # or test on self.embedding_charges ? + if specifications[chem]["program"] == "psi4": + charges = imol.extras["embedding_charges"] + fkw = inp.keywords.get("function_kwargs", {}) + fkw.update({"external_potentials": charges}) + inp.keywords["function_kwargs"] = fkw + else: + raise RuntimeError(f"Don't know how to handle external charges in {specifications[chem]['program']}") + _, real, bas = delabeler(label) result = qcng.compute(inp, specifications[chem]["program"]) diff --git a/qcmanybody/tests/test_mbe_het4_grad.py b/qcmanybody/tests/test_mbe_het4_grad.py new file mode 100644 index 0000000..2a04f0f --- /dev/null +++ b/qcmanybody/tests/test_mbe_het4_grad.py @@ -0,0 +1,311 @@ +import re +import pprint + +import pytest +import numpy as np + +from qcelemental.models import Molecule +from qcelemental.testing import compare_values + +from qcmanybody.models.manybody_pydv1 import ManyBodyKeywords, ManyBodyInput +from qcmanybody.models.qcng_computer import ManyBodyComputerQCNG, qcvars_to_manybodyproperties + +from .addons import uusing + +def skprop(qcvar): + # qcng: return qcng.procedures.manybody.qcvars_to_manybodyproperties[qcvar] + return qcvars_to_manybodyproperties[qcvar] + + +@pytest.fixture(scope="function") +def mbe_data_grad_dtz(): + # note that spherical/cartesian irrelevant for He & 6-31G, and fc/ae irrelevant for He + p4_kwds = {"scf_type": "pk", "mp2_type": "conv"} + + protocols = {"stdout": False} + return { + "specification": { + "specification": { + "p4-hf": { + "model": { + "method": "hf", + "basis": "cc-pvdz", + }, + "driver": "energy", + "program": "psi4", + "keywords": p4_kwds, + "protocols": protocols, + }, + "p4-mp2": { + "model": { + "method": "mp2", + "basis": "cc-pvdz", + }, + "driver": "energy", + "program": "psi4", + "keywords": p4_kwds, + "protocols": protocols, + }, + "p4-ccsd": { + "model": { + "method": "ccsd", + "basis": "cc-pvdz", + }, + "driver": "energy", + "program": "psi4", + "keywords": p4_kwds, + "protocols": protocols, + }, + }, + "keywords": None, + "driver": "energy", + }, + "molecule": None, + } + + +het4_refs_conv_grad_dtz = { + # 4: hf/cc-pvdz ; copied from psi4 - not computed from pieces + "4b_nocp_rtd": { + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -203.05040290887501, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -203.35196718523923, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -203.32671800908514, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -203.32531316984924, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": -0.3015642763642177, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": -0.27631510021012673, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": -0.2749102609742238, + "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": -0.3015642763642177, + "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": 0.025249176154090947, + "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.001404839235902955, + + "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY": np.array([ + [-7.40362091e-02, 2.84130342e-03, 1.00812419e-02], + [ 7.57644607e-02, -1.48383966e-04, -1.73598306e-03], + [ 2.60812740e-02, 4.40153410e-04, 2.02181413e-03], + [-8.52081479e-03, -2.26666899e-02, -1.10830218e-01], + [-6.03761361e-03, -1.54080944e-02, 1.00510650e-01], + [-1.32510973e-02, 3.49417115e-02, -4.75050095e-05]]), + }, + "4b_nocp_rtd_ee": { + # through 3b from psi4. 4b from qcmb! + # TODO IE should be present" + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -203.4340927642202, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -203.33800560627736, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -203.32602370522727, + #"NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.09608715794283285, + #"NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.10806905899292474, + "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.09608715794283285, + "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": 0.011981901050091892, + #"NOCP-CORRECTED TOTAL GRADIENT THROUGH 3-BODY": np.array([ # fine but "wrongly present" + # [-7.36025774e-02, 3.36991350e-03, 1.00686673e-02], + # [ 7.55655033e-02, -5.37981728e-04, -2.19873195e-03], + # [ 2.63105970e-02, 4.12475495e-04, 2.10686911e-03], + # [-8.78129292e-03, -2.22190993e-02, -1.10826802e-01], + # [-6.14121135e-03, -1.53249058e-02, 1.00569730e-01], + # [-1.32579782e-02, 3.44008392e-02, 3.46216252e-05]]), + + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -203.325313169849, + "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000710535378, + "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY": np.array([ + [-7.40362107e-02, 2.84130177e-03, 1.00812417e-02], + [ 7.57644630e-02, -1.48382170e-04, -1.73598308e-03], + [ 2.60812737e-02, 4.40153586e-04, 2.02181449e-03], + [-8.52081509e-03, -2.26666902e-02, -1.10830218e-01], + [-6.03761335e-03, -1.54080947e-02, 1.00510650e-01], + [-1.32510976e-02, 3.49417117e-02, -4.75050664e-05]]), + }, +} + + + +# only here for keys +he4_refs_conv = { + "CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": None, + "CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": None, + "CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": None, + "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": None, + "CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": None, + "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": None, + "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": None, + "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": None, + "CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": None, + "CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": None, + "CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": None, + + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": None, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": None, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": None, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": None, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": None, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": None, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": None, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": None, + "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": None, + "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": None, + "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": None, + + "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY": None, + + "VMFC-CORRECTED TOTAL ENERGY THROUGH 1-BODY": None, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 2-BODY": None, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY": None, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 4-BODY": None, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 1-BODY": None, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": None, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": None, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": None, + "VMFC-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": None, + "VMFC-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": None, + "VMFC-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": None, +} + + +sumstr = { + "4b_nocp_rtd": r""" + ==> N-Body: Non-Counterpoise Corrected \(NoCP\) energies <== + +^\s+n-Body Total Energy Interaction Energy N-body Contribution to Interaction Energy +^\s+\[Eh\] \[Eh\] \[kcal/mol\] \[Eh\] \[kcal/mol\] +^\s+1\s+-203.0504029\d+ 0.0000000\d+ 0.00\d+ 0.0000000\d+ 0.00\d+ +^\s+2\s+-203.3519671\d+ -0.3015642\d+ -189.23\d+ -0.3015642\d+ -189.23\d+ +^\s+3\s+-203.3267180\d+ -0.2763151\d+ -173.39\d+ 0.0252491\d+ 15.84\d+ +^\s+FULL/RTN\s+4 -203.3253131\d+ -0.2749102\d+ -172.50\d+ 0.0014048\d+ 0.88\d+ +""", + "4b_nocp_rtd_ee": r""" +^\s+==> N-Body: Non-Counterpoise Corrected \(NoCP\) energies <== + +^\s+n-Body Total Energy Interaction Energy N-body Contribution to Interaction Energy +^\s+\[Eh\] \[Eh\] \[kcal/mol\] \[Eh\] \[kcal/mol\] +^\s+1\s+-203.4340927\d+\s+N/A\s+N/A\s+0.0000000\d+ 0.00\d+ +^\s+2\s+-203.3380056\d+\s+N/A\s+N/A\s+0.0960871\d+ 60.29\d+ +^\s+3\s+-203.3260237\d+\s+N/A\s+N/A\s+0.0119819\d+ 7.51\d+ +^\s+FULL/RTN\s+4\s+-203.3253131\d+\s+N/A\s+N/A\s+0.0007105\d+\s+0.44\d+ +""", +} + +sumdict_grad = { + "4b_all": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + "CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + "VMFC-CORRECTED TOTAL ENERGY": "VMFC-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "VMFC-CORRECTED INTERACTION ENERGY": "VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + + "NOCP-CORRECTED TOTAL GRADIENT": "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + #"NOCP-CORRECTED INTERACTION GRADIENT": "NOCP-CORRECTED INTERACTION GRADIENT THROUGH 4-BODY", + "CP-CORRECTED TOTAL GRADIENT": "CP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + #"CP-CORRECTED INTERACTION GRADIENT": "CP-CORRECTED INTERACTION GRADIENT THROUGH 4-BODY", + "VMFC-CORRECTED TOTAL GRADIENT": "VMFC-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + #"VMFC-CORRECTED INTERACTION GRADIENT": "VMFC-CORRECTED INTERACTION GRADIENT THROUGH 4-BODY", + }, + "4b_nocp_rtd": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + + "NOCP-CORRECTED TOTAL GRADIENT": "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + #"NOCP-CORRECTED INTERACTION GRADIENT": "NOCP-CORRECTED INTERACTION GRADIENT THROUGH 4-BODY", + }, + "4b_nocp_rtd_ee": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + #"NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + + "NOCP-CORRECTED TOTAL GRADIENT": "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + #"NOCP-CORRECTED INTERACTION GRADIENT": "NOCP-CORRECTED INTERACTION GRADIENT THROUGH 4-BODY", + }, +} + + +@pytest.fixture +def het_tetramer(): + return Molecule( + symbols=["F", "H", "F", "H", "H", "He"], + fragments=[[0], [1, 2], [3, 4], [5]], + fragment_charges=[-1, 0, 0, 0], + geometry=[-2, 0, 0, 0, 0, 0, 4, 0, 0, 0, 3, 0, 0, 3, 2, 0, -3, 0], + ) + + +@uusing("psi4") +@pytest.mark.parametrize("mbe_keywords,anskeyE,anskeyG,bodykeys,outstrs,calcinfo_nmbe", [ + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "levels": {4: "p4-hf"}}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + [k for k in he4_refs_conv if k.startswith("NOCP-")], + ["4b_nocp_rtd"], + 15, + #All,4b QCVariable: 1_((1, 2, 3, 4), (1, 2, 3, 4)) -203.32531316984924 + id="4b_nocp_rtd"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "levels": {4: "p4-hf"}, "embedding_charges": {1: [-1.0], 2: [0.5, -0.5], 3: [-0.5, 0.5], 4: [0]}}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED TOTAL GRADIENT THROUGH 4-BODY", + [k for k in he4_refs_conv if k.startswith("NOCP-")], + ["4b_nocp_rtd_ee"], + 15, + id="4b_nocp_rtd_ee"), +]) +def test_nbody_het4_grad(mbe_keywords, anskeyE, anskeyG, bodykeys, outstrs, calcinfo_nmbe, het_tetramer, request, mbe_data_grad_dtz): + _inner = request.node.name.split("[")[1].split("]")[0] + kwdsln, pattern, progln = _inner, "", "psi4" + print("LANE", kwdsln, pattern, progln) + + mbe_keywords = ManyBodyKeywords(**mbe_keywords) + mbe_data_grad_dtz["molecule"] = het_tetramer + mbe_data_grad_dtz["specification"]["driver"] = "gradient" + mbe_data_grad_dtz["specification"]["keywords"] = mbe_keywords + mbe_model = ManyBodyInput(**mbe_data_grad_dtz) + + # qcng: ret = qcng.compute_procedure(mbe_model, "manybody", raise_error=True) + ret = ManyBodyComputerQCNG.from_qcschema_ben(mbe_model) + print(f"MMMMMMM {request.node.name}") + pprint.pprint(ret.dict(), width=200) + + refs = het4_refs_conv_grad_dtz[kwdsln] + ansE = refs[anskeyE] + ansG = refs[anskeyG] + ref_nmbe = calcinfo_nmbe + ref_bodykeys = bodykeys + ref_sumdict = sumdict_grad[kwdsln] + atol = 2.5e-8 + + for qcv, ref in refs.items(): + skp = skprop(qcv) + if qcv in ref_bodykeys: + assert compare_values(ref, ret.extras["qcvars"]["nbody"][qcv], atol=atol, label=f"[a] qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=atol, label=f"[b] skprop {skp}") + else: + assert qcv not in ret.extras["qcvars"]["nbody"], f"[z] {qcv=} wrongly present" + assert getattr(ret.properties, skp) is None + + for qcv in sumdict_grad["4b_all"]: + skp = skprop(qcv) + if qcv in ref_sumdict: + ref = refs[ref_sumdict[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: + assert qcv not in ret.extras["qcvars"]["nbody"], f"[y] {qcv=} wrongly present" + assert getattr(ret.properties, skp) is None + + for qcv, ref in { + "CURRENT ENERGY": ansE, + }.items(): + skp = skprop(qcv) + assert compare_values(ref, ret.extras["qcvars"][qcv], atol=atol, label=f"[e] qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=atol, label=f"[f] skprop {skp}") + + for qcv, ref in { + "CURRENT GRADIENT": ansG, + }.items(): + skp = skprop(qcv) + assert compare_values(ref, ret.extras["qcvars"][qcv], atol=atol, label=f"[e] G qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=atol, label=f"[f] G skprop {skp}") + assert compare_values(ansG, ret.return_result, atol=atol, label=f"[g] G ret") + + assert ret.properties.calcinfo_nmbe == ref_nmbe, f"{ret.properties.calcinfo_nmbe=} != {ref_nmbe}" + + if outstrs: + for stdoutkey in outstrs: + assert re.search(sumstr[stdoutkey], ret.stdout, re.MULTILINE), f"[j] N-Body pattern not found: {sumstr[stdoutkey]}" diff --git a/qcmanybody/tests/test_mbe_keywords.py b/qcmanybody/tests/test_mbe_keywords.py index df695c9..05e4346 100644 --- a/qcmanybody/tests/test_mbe_keywords.py +++ b/qcmanybody/tests/test_mbe_keywords.py @@ -272,7 +272,7 @@ def test_noncontiguous_fragments_evaded(): neon_in_hcl = Molecule(symbols=["H", "Ne", "Cl"], geometry=[0, 0, 0, 2, 0, 0, 0, 2, 0], fragments=[[0, 2], [1]], validated=True) with pytest.raises(ValueError) as e: - ManyBodyCalculator(neon_in_hcl, ["cp"], {2: "mp2", 1: "mp2"}, False, False) + ManyBodyCalculator(neon_in_hcl, ["cp"], {2: "mp2", 1: "mp2"}, False, False, None) assert "QCManyBody: non-contiguous fragments could be implemented but aren't at present" in str(e.value) diff --git a/qcmanybody/tests/test_multi.py b/qcmanybody/tests/test_multi.py index 67e9c8d..a08ed55 100644 --- a/qcmanybody/tests/test_multi.py +++ b/qcmanybody/tests/test_multi.py @@ -25,6 +25,6 @@ def test_h2o_trimer_multi(levels, component_file, ref_file): component_results = load_component_data(component_file) ref_data = load_ref_data(ref_file) - mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False, None) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels) diff --git a/qcmanybody/tests/test_multi_ss.py b/qcmanybody/tests/test_multi_ss.py index 872f35a..6807b11 100644 --- a/qcmanybody/tests/test_multi_ss.py +++ b/qcmanybody/tests/test_multi_ss.py @@ -28,6 +28,6 @@ def test_h2o_trimer_multi_ss(levels, component_file, ref_file): component_results = load_component_data(component_file) ref_data = load_ref_data(ref_file) - mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False, None) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels) diff --git a/qcmanybody/tests/test_single.py b/qcmanybody/tests/test_single.py index 9e81748..ffeee76 100644 --- a/qcmanybody/tests/test_single.py +++ b/qcmanybody/tests/test_single.py @@ -25,6 +25,6 @@ def test_h2o_trimer_single(levels, component_file, ref_file): component_results = load_component_data(component_file) ref_data = load_ref_data(ref_file) - mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False, None) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels) diff --git a/qcmanybody/tests/utils.py b/qcmanybody/tests/utils.py index b66a2f2..5aabde8 100644 --- a/qcmanybody/tests/utils.py +++ b/qcmanybody/tests/utils.py @@ -1,7 +1,7 @@ import json import math import os -from typing import Mapping, Union, Literal, Any, Iterable +from typing import Mapping, Union, Literal, Any, Iterable, Optional import numpy import qcengine as qcng @@ -168,10 +168,11 @@ def run_qcengine( specifications: Mapping[str, Mapping[str, Any]], bsse_type: Iterable[BsseEnum], return_total_data: bool, - supersystem_ie_only: bool + supersystem_ie_only: bool, + embedding_charges: Optional[Mapping[int, list]], ): - mc = ManyBodyCalculator(molecule, bsse_type, levels, return_total_data, supersystem_ie_only) + mc = ManyBodyCalculator(molecule, bsse_type, levels, return_total_data, supersystem_ie_only, embedding_charges) component_results = {} diff --git a/qcmanybody/utils.py b/qcmanybody/utils.py index 123e133..cb4f862 100644 --- a/qcmanybody/utils.py +++ b/qcmanybody/utils.py @@ -306,7 +306,10 @@ def print_nbody_energy( info += f""" {lbl:>11} {nb:3} {"N/A":14} {int_e:20.12f} {int_e_kcal:20.12f} {"N/A":20} {"N/A":20}\n""" else: if tot_e: - info += f""" {lbl:>11} {nb:3} {energy_body_dict[nb]:20.12f} {int_e:20.12f} {int_e_kcal:20.12f} {delta_e:20.12f} {delta_e_kcal:20.12f}\n""" + if embedding: + info += f""" {lbl:>11} {nb:3} {energy_body_dict[nb]:20.12f} {"N/A":20} {"N/A":14} {delta_e:20.12f} {delta_e_kcal:20.12f}\n""" + else: + info += f""" {lbl:>11} {nb:3} {energy_body_dict[nb]:20.12f} {int_e:20.12f} {int_e_kcal:20.12f} {delta_e:20.12f} {delta_e_kcal:20.12f}\n""" else: info += f""" {lbl:>11} {nb:3} {"N/A":14} {int_e:20.12f} {int_e_kcal:20.12f} {delta_e:20.12f} {delta_e_kcal:20.12f}\n""" previous_e = energy_body_dict[nb] @@ -355,8 +358,6 @@ def collect_vars( nbody_range = list(body_dict) nbody_range.sort() res = {} - if embedding: - return res if tot_e: res[f"{bsse}-CORRECTED TOTAL {prop}"] = body_dict[max_nbody] @@ -393,6 +394,9 @@ def collect_vars( for nb in nbody_range: res[f"{bsse}-CORRECTED TOTAL {prop} THROUGH {nb}-BODY"] = body_dict[nb] + if embedding: + res = {k: v for k, v in res.items() if "INTERACTION" not in k} + return res