diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 82cb539..97e1ef0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -27,8 +27,7 @@ jobs: - label: Py-max python-version: "3.12" runs-on: ubuntu-latest - pytest: "" - #pytest: "-k 'not (hes2 or qchem)'" + pytest: "-k 'not (he4 and all)'" name: "🐍 ${{ matrix.cfg.python-version }} • ${{ matrix.cfg.label }} • ${{ matrix.cfg.runs-on }}" runs-on: ${{ matrix.cfg.runs-on }} @@ -58,11 +57,16 @@ jobs: - zstandard - pytest-cov - codecov + # Testing CMS + - nwchem + - networkx EOF if [[ "${{ runner.os }}" == "Linux" ]]; then : if [[ "${{ matrix.cfg.label }}" == "Py-min" ]]; then sed -i "s;pydantic;pydantic=1;g" export.yaml + sed -i "s;- nwchem;#- nwchem;g" export.yaml + sed -i "s;- networkx;#- networkx;g" export.yaml fi fi # model sed for L/W @@ -105,7 +109,7 @@ jobs: - name: PyTest run: | - pytest -rws -v ${{ matrix.cfg.pytest }} --cov=qcmanybody --color=yes --cov-report=xml qcmanybody/ tests/ + pytest -rws -v ${{ matrix.cfg.pytest }} --cov=qcmanybody --color=yes --cov-report=xml --durations 50 --durations-min 5 qcmanybody/ tests/ - name: CodeCov uses: codecov/codecov-action@v3 diff --git a/qcmanybody/builder.py b/qcmanybody/builder.py index 0a5093e..5516db4 100644 --- a/qcmanybody/builder.py +++ b/qcmanybody/builder.py @@ -11,6 +11,7 @@ def build_nbody_compute_list( nfragments: int, nbodies: Iterable[Union[int, Literal["supersystem"]]], return_total_data: bool, + supersystem_ie_only: bool, supersystem_max_nbody: Optional[int] = None, ) -> Dict[str, Dict[int, Set[FragBasIndex]]]: """Generates lists of N-Body computations needed for requested BSSE treatments. @@ -27,6 +28,8 @@ def build_nbody_compute_list( 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 + ???? Returns ------- @@ -81,18 +84,28 @@ 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): - cp_compute_list[nb].add((x, basis_tuple)) + if supersystem_ie_only: + for sublevel in [1, nfragments]: + for x in itertools.combinations(fragment_range, sublevel): + cp_compute_list[nfragments].add((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): + cp_compute_list[nb].add((x, basis_tuple)) if "nocp" in bsse_type: # Everything in monomer basis - for nb in nbodies: - for sublevel in range(1, nb + 1): + if supersystem_ie_only: + for sublevel in [1, nfragments]: for x in itertools.combinations(fragment_range, sublevel): - nocp_compute_list[nb].add((x, x)) + nocp_compute_list[nfragments].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 diff --git a/qcmanybody/manybody.py b/qcmanybody/manybody.py index 8370126..79ee0b9 100644 --- a/qcmanybody/manybody.py +++ b/qcmanybody/manybody.py @@ -33,6 +33,7 @@ def __init__( bsse_type: Sequence[BsseEnum], levels: Mapping[Union[int, Literal["supersystem"]], str], return_total_data: bool, + supersystem_ie_only: bool, ): # TODO self.embedding_charges = {} @@ -40,6 +41,7 @@ def __init__( self.molecule = molecule self.bsse_type = [BsseEnum(x) for x in bsse_type] self.return_total_data = return_total_data + self.supersystem_ie_only = supersystem_ie_only self.nfragments = len(molecule.fragments) self.levels = levels @@ -61,6 +63,8 @@ def __init__( # Build nbodies_per_mc_level # TODO - use Lori's code # TODO - dict to list of lists to handle non-contiguous levels + # TODO multilevel and supersystem_ie_only=T not allowed together + # TODO supersystem in levels is not to be trusted -- nfrag only and skips levels max_level = max(self.levels_no_ss.keys()) if set(range(1, max_level + 1)) != set(self.levels_no_ss.keys()): @@ -109,6 +113,7 @@ def compute_map(self) -> Dict[str, Dict[str, Dict[int, Set[FragBasIndex]]]]: self.nfragments, nbodies, self.return_total_data, + self.supersystem_ie_only, self.max_nbody, ) @@ -435,7 +440,7 @@ def analyze( 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) + collect_vars(b.upper(), all_results["energy_body_dict"][b], self.max_nbody, is_embedded, self.supersystem_ie_only) ) all_results["results"] = nbody_dict diff --git a/qcmanybody/models/manybody_v1.py b/qcmanybody/models/manybody_v1.py index 4215e87..458080a 100644 --- a/qcmanybody/models/manybody_v1.py +++ b/qcmanybody/models/manybody_v1.py @@ -130,7 +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``." + "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." ) # v2: @field_validator("bsse_type", mode="before") diff --git a/qcmanybody/models/qcng_computer.py b/qcmanybody/models/qcng_computer.py index d26df4e..1d6829c 100644 --- a/qcmanybody/models/qcng_computer.py +++ b/qcmanybody/models/qcng_computer.py @@ -1,18 +1,29 @@ import abc import copy + +# printing and logging formatting niceties import pprint +from functools import partial +import numpy as np +pp = pprint.PrettyPrinter(width=120, compact=True, indent=1) +nppp = partial(np.array_str, max_line_width=120, precision=8, suppress_small=True) +nppp10 = partial(np.array_str, max_line_width=120, precision=10, suppress_small=True) + from ast import literal_eval # v2: from typing import TYPE_CHECKING, Any, ClassVar, Dict, List, Tuple, Union, Literal, Optional -from typing import Any, Dict, List, Tuple, Union, Literal, Optional +from typing import Any, Dict, List, Mapping, Tuple, Union, Literal, Optional # v2: from pydantic import ConfigDict, field_validator, FieldValidationInfo, computed_field, BaseModel, Field try: - from pydantic.v1 import create_model, Field, validator + from pydantic.v1 import create_model, Field, validator, BaseModel except ImportError: - from pydantic import create_model, Field, validator + from pydantic import create_model, Field, validator, BaseModel from qcelemental.models import FailedOperation, Molecule, DriverEnum, ProtoModel, AtomicResult, AtomicInput +import qcengine as qcng from .manybody_v1 import BsseEnum, ManyBodyKeywords, ManyBodyInput, ManyBodyResult, ManyBodyResultProperties +from qcmanybody import ManyBodyCalculator +from qcmanybody.utils import delabeler, provenance_stamp class BaseComputerQCNG(ProtoModel): @@ -88,7 +99,7 @@ def compute(self, client: Optional["qcportal.FractalClient"] = None) -> None: NOTE: client logic removed (compared to psi4.driver.AtomicComputer) """ - from ..compute import compute as qcng_compute + from qcengine import compute as qcng_compute if self.computed: return @@ -350,9 +361,74 @@ def set_supersystem_ie_only(cls, v: Optional[bool], values) -> bool: if (sio is True) and (_max_nbody != _nfr): raise ValueError(f"Cannot skip intermediate n-body jobs when max_nbody={_max_nbody} != nfragments={_nfr}.") + if (sio is True) and ("vmfc" in values["bsse_type"]): + raise ValueError(f"Cannot skip intermediate n-body jobs when VMFC in bsse_type={values['bsse_type']}. Use CP instead.") + print(f" ... setting {sio=}") return sio + @classmethod + def from_qcschema_ben(cls, input_model: ManyBodyInput): + + computer_model = cls( + molecule=input_model.molecule, + driver=input_model.specification.driver, + # v2: **input_model.specification.keywords.model_dump(exclude_unset=True), + **input_model.specification.keywords.dict(exclude_unset=True), + input_data=input_model, # storage, to reconstitute ManyBodyResult + ) + nb_per_mc = computer_model.nbodies_per_mc_level + + comp_levels = {} + for mc_level_idx, mtd in enumerate(computer_model.levels.values()): + for lvl1 in nb_per_mc[mc_level_idx]: + comp_levels[int(lvl1)] = mtd + + specifications = {} + for mtd, spec in computer_model.input_data.specification.specification.items(): + spec = spec.dict() + specifications[mtd] = {} + specifications[mtd]["program"] = spec.pop("program") + specifications[mtd]["specification"] = spec + specifications[mtd]["specification"].pop("schema_name", None) + specifications[mtd]["specification"].pop("protocols", None) + + calculator_cls = ManyBodyCalculator( + computer_model.molecule, + computer_model.bsse_type, + comp_levels, + computer_model.return_total_data, + computer_model.supersystem_ie_only, + ) + + component_results = {} + + for chem, label, imol in calculator_cls.iterate_molecules(): + inp = AtomicInput(molecule=imol, **specifications[chem]["specification"]) + + _, real, bas = delabeler(label) + result = qcng.compute(inp, specifications[chem]["program"]) + + if not result.success: + print(result.error.error_message) + raise RuntimeError("Calculation did not succeed! Error:\n" + result.error.error_message) + + # pull out stuff + props = {"energy", "gradient", "hessian"} + + component_results[label] = {} + + for p in props: + if hasattr(result.properties, f"return_{p}"): + v = getattr(result.properties, f"return_{p}") + if v is not None: + component_results[label][p] = v + + analyze_back = calculator_cls.analyze(component_results) + analyze_back["nbody_number"] = len(component_results) + + return computer_model.get_results(external_results=analyze_back) + @classmethod def from_qcschema(cls, input_model: ManyBodyInput, build_tasks: bool = False): @@ -421,7 +497,7 @@ def build_tasks( Formerly, didn't include supersystem in count. """ - from psi4.driver.driver_nbody import build_nbody_compute_list +# qcng: from psi4.driver.driver_nbody import build_nbody_compute_list # TODO method not coming from levels right @@ -456,7 +532,7 @@ def build_tasks( self.bsse_type, nbodies, self.nfragments, self.return_total_data, self.supersystem_ie_only) - def labeler(item) -> str: + def lab_labeler(item) -> str: # note 0-index to 1-index shift for label return f"{mc_level_idx + 1}_{item}" @@ -466,7 +542,7 @@ def labeler(item) -> str: # and second is basis set fragment indices, all 1-indexed for nb in compute_dict["all"]: for pair in compute_dict["all"][nb]: - lbl = labeler(pair) + lbl = lab_labeler(pair) if lbl in self.task_list: continue @@ -498,9 +574,9 @@ def compute(self, client: Optional["qcportal.FractalClient"] = None) -> None: NOTE: client logic removed (compared to psi4.driver.ManyBodyComputer) """ - from psi4.driver.p4util import banner +# qcng: from psi4.driver.p4util import banner - info = "\n" + banner(f" ManyBody Computations ", strNotOutfile=True) + "\n" +# qcng: info = "\n" + banner(f" ManyBody Computations ", strNotOutfile=True) + "\n" #logger.info(info) for t in self.task_list.values(): @@ -539,8 +615,8 @@ def prepare_results( # rtd := return_total_data # """ - from psi4.driver.driver_nbody import assemble_nbody_components - from psi4.driver.driver_nbody_multilevel import prepare_results as multilevel_prepare_results +# qcng: from psi4.driver.driver_nbody import assemble_nbody_components +# qcng: from psi4.driver.driver_nbody_multilevel import prepare_results as multilevel_prepare_results if results is None: print(f"RESULTS setting empty") @@ -600,7 +676,7 @@ def prepare_results( nbody_results["intermediates"] = {} for idx, task in results_list.items(): - mc, frag, bas = delabeler(idx) + mc, frag, bas = lab_delabeler(idx) nbody_results["intermediates"][f"N-BODY ({frag})@({bas}) TOTAL ENERGY"] = task.properties.return_energy nbody_results["intermediates_energy"] = trove["energy"] @@ -629,38 +705,51 @@ def prepare_results( return nbody_results - def get_results(self, client: Optional["qcportal.FractalClient"] = None) -> ManyBodyResult: + def get_results(self, client: Optional["qcportal.FractalClient"] = None, external_results: Optional[Dict] = None) -> ManyBodyResult: """Return results as ManyBody-flavored QCSchema. NOTE: client removed (compared to psi4.driver.ManyBodyComputer) """ - from psi4.driver.p4util import banner +# qcng: from psi4.driver.p4util import banner - info = "\n" + banner(f" ManyBody Results ", strNotOutfile=True) + "\n" +# qcng: info = "\n" + banner(f" ManyBody Results ", strNotOutfile=True) + "\n" #logger.info(info) - results = self.prepare_results(client=client) - ret_energy = results.pop("ret_energy") - ret_ptype = results.pop("ret_ptype") - ret_gradient = results.pop("ret_gradient", None) + if external_results is None: + results = self.prepare_results(client=client) + ret_energy = results.pop("ret_energy") + ret_ptype = results.pop("ret_ptype") + ret_gradient = results.pop("ret_gradient", None) + nbody_number = len(self.task_list) + else: + ret_energy = external_results.pop("ret_energy") + ret_ptype = ret_energy if self.driver == "energy" else external_results.pop(f"ret_{self.driver.name}") + ret_gradient = external_results.pop("ret_gradient", None) + nbody_number = external_results.pop("nbody_number") # load QCVariables qcvars = { 'NUCLEAR REPULSION ENERGY': self.molecule.nuclear_repulsion_energy(), - 'NBODY NUMBER': len(self.task_list), + 'NBODY NUMBER': nbody_number, } properties = { "calcinfo_nmc": len(self.nbodies_per_mc_level), "calcinfo_nfr": self.nfragments, # or len(self.molecule.fragments) "calcinfo_natom": len(self.molecule.symbols), - "calcinfo_nmbe": len(self.task_list), + "calcinfo_nmbe": nbody_number, "nuclear_repulsion_energy": self.molecule.nuclear_repulsion_energy(), "return_energy": ret_energy, } - for k, val in results.items(): - qcvars[k] = val + if external_results is None: + for k, val in results.items(): + qcvars[k] = val + else: + for k, val in external_results.items(): + if k == "results": + k = "nbody" + qcvars[k] = val qcvars['CURRENT ENERGY'] = ret_energy if self.driver == 'gradient': @@ -758,7 +847,7 @@ def get_results(self, client: Optional["qcportal.FractalClient"] = None) -> Many -def delabeler(item: str, return_obj: bool = False) -> Union[Tuple[str, str, str], Tuple[int, Tuple[int], Tuple[int]]]: +def lab_delabeler(item: str, return_obj: bool = False) -> Union[Tuple[str, str, str], Tuple[int, Tuple[int], Tuple[int]]]: """Transform labels like string "1_((2,), (1, 2))" into string tuple ("1", "2", "1, 2") or object tuple (1, (2,), (1, 2)). @@ -782,7 +871,7 @@ def delabeler(item: str, return_obj: bool = False) -> Union[Tuple[str, str, str] qcvars_to_manybodyproperties["CURRENT HESSIAN"] = "return_hessian" -def build_manybodyproperties(qcvars: Dict) -> ManyBodyResultProperties: +def build_manybodyproperties(qcvars: Mapping) -> ManyBodyResultProperties: """For results extracted from QC output in QCDB terminology, translate to QCSchema terminology. Parameters diff --git a/qcmanybody/models/test_mbe_he4_singlelevel.py b/qcmanybody/models/test_mbe_he4_singlelevel.py new file mode 100644 index 0000000..6cc402e --- /dev/null +++ b/qcmanybody/models/test_mbe_he4_singlelevel.py @@ -0,0 +1,375 @@ +import pprint + +import pytest + +from qcelemental import constants +from qcelemental.models import Molecule +# v2: from qcelemental.models.procedures_manybody import AtomicSpecification, ManyBodyKeywords, ManyBodyInput +from qcelemental.testing import compare_values + +from qcmanybody.models.manybody_v1 import AtomicSpecification, ManyBodyKeywords, ManyBodyInput +from qcmanybody.models.qcng_computer import ManyBodyComputerQCNG, qcvars_to_manybodyproperties + +import qcengine as qcng +from qcengine.testing import using + +def skprop(qcvar): + # qcng: return qcng.procedures.manybody.qcvars_to_manybodyproperties[qcvar] + return qcvars_to_manybodyproperties[qcvar] + + +he4_refs_conv = { + "CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530668717083888, + "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 2-BODY": 0.008200959993875045, + "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.007965853003739198, + "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008028846432448944, + "CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.008200959993875045, + "CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.00023510699013584713, + "CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 6.299342870974556e-05, + + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530668717083888, + "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 2-BODY": 0.007817510905059777, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.0075734474125397355, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.007630623419519367, + "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.007817510905059777, + "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.00024406349252004134, + "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 5.717600697963121e-05, + + "VMFC-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530668717083888, + "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 2-BODY": 0.00821979538669737, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.007984194798996924, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.00804718822770667, + "VMFC-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.00821979538669737, + "VMFC-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.00023560058770044634, + "VMFC-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 6.299342870974556e-05, +} + +he4_refs_df = { + "CP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530751941948, + "CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.522403579651, + "CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.522640167467, + "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.522576639404, + "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008348362297, + "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.008111774481, + "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008175302544, + "CP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.008348362297, + "CP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000236587816, + "CP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000063528063, + + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530751941948, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.522760073327, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.523005411447, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.522948420000, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.007991868621, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.007746530501, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.007803521948, + "NOCP-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.007991868621, + "NOCP-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000245338120, + "NOCP-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000056991448, + + "VMFC-CORRECTED TOTAL ENERGY THROUGH 1-BODY": -11.530751941948, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 2-BODY": -11.522390319401, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 3-BODY": -11.522627256726, + "VMFC-CORRECTED TOTAL ENERGY THROUGH 4-BODY": -11.522563728663, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 2-BODY": 0.008361622547, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 3-BODY": 0.008124685222, + "VMFC-CORRECTED INTERACTION ENERGY THROUGH 4-BODY": 0.008188213285, + "VMFC-CORRECTED 2-BODY CONTRIBUTION TO ENERGY": 0.008361622547, + "VMFC-CORRECTED 3-BODY CONTRIBUTION TO ENERGY": -0.000236937325, + "VMFC-CORRECTED 4-BODY CONTRIBUTION TO ENERGY": 0.000063528063, + } + +sumdict = { + "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", + }, + "4b-nocp-rtd-sio": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "4b-nocp-sio": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "4b-cp-rtd-sio": { + "CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "4b-cp-sio": { + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY 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", + }, + "4b-nocp": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "4b-cp-rtd": { + "CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "4b-cp": { + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + }, + "3b-nocp-rtd": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + }, + "3b-nocp": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + }, + "3b-cp-rtd": { + "CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY", + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + }, + "3b-cp": { + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + }, + "2b-nocp-rtd": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + }, + "2b-nocp": { + "NOCP-CORRECTED TOTAL ENERGY": "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY", + "NOCP-CORRECTED INTERACTION ENERGY": "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + }, + "2b-cp-rtd": { + "CP-CORRECTED TOTAL ENERGY": "CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY", + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + }, + "2b-cp": { + "CP-CORRECTED INTERACTION ENERGY": "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + }, + "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", + }, + "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", + }, + "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", + }, + "1b-cp": { + "CP-CORRECTED INTERACTION ENERGY": "zero", #"CP-CORRECTED INTERACTION ENERGY THROUGH 1-BODY", + }, +# TODO table defines the general qcvar as 0 even if 1-body qcvar not available. continue? +} + + +@pytest.fixture +def he_tetramer(): + a2 = 2 / constants.bohr2angstroms + return Molecule(symbols=["He", "He", "He", "He"], fragments=[[0], [1], [2], [3]], geometry=[0, 0, 0, 0, 0, a2, 0, a2, 0, 0, a2, a2]) + + +@pytest.mark.parametrize("program,basis,keywords", [ + pytest.param("cfour", "aug-pvdz", {"frozen_core": False}, id="cfour-conv", marks=using("cfour")), + pytest.param("gamess", "accd", {"contrl__ispher": 1, "mp2__nacore": 0}, id="gamess-conv", marks=using("gamess")), + pytest.param("nwchem", "aug-cc-pvdz", {"basis__spherical": True, "scf__thresh": 1.0e-8, "mp2__freeze": False}, id="nwchem-conv", marks=using("nwchem")), + pytest.param("psi4", "aug-cc-pvdz", {"e_convergence": 1.e-10, "d_convergence": 1.e-10, "scf_type": "pk", "mp2_type": "conv"}, id="psi4-conv", marks=using("psi4")), + pytest.param("psi4", "aug-cc-pvdz", {"e_convergence": 1.e-10, "d_convergence": 1.e-10}, id="psi4-df", marks=using("psi4")), +]) +@pytest.mark.parametrize("mbe_keywords,anskey,bodykeys,calcinfo_nmbe", [ + pytest.param( + {"bsse_type": ["nocp", "cp", "vmfc"]}, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv], + 65, + id="4b-all"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "supersystem_ie_only": True}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("THROUGH 4-BODY" in k or "THROUGH 1-BODY" in k))], + 5, + id="4b-nocp-rtd-sio"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": False, "supersystem_ie_only": True}, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("THROUGH 4-BODY" in k or "THROUGH 1-BODY" in k))], + 5, + id="4b-nocp-sio"), + pytest.param( + {"bsse_type": "cp", "return_total_data": True, "supersystem_ie_only": True}, + "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and ("THROUGH 4-BODY" in k or "THROUGH 1-BODY" in k))], + 9, + id="4b-cp-rtd-sio"), + pytest.param( + {"bsse_type": "cp", "return_total_data": False, "supersystem_ie_only": True}, + "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and ("THROUGH 4-BODY" in k or "THROUGH 1-BODY" in k) and "TOTAL ENERGY" not in k)], + 5, + id="4b-cp-sio"), +## TODO add vmfc. 3b nmbe=50 + pytest.param( + {"bsse_type": "nocp", "return_total_data": True}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-"))], + 15, + id="4b-nocp-rtd"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": False}, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-"))], + 15, + id="4b-nocp"), + pytest.param( + {"bsse_type": "cp", "return_total_data": True}, + "CP-CORRECTED TOTAL ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-"))], + 19, + id="4b-cp-rtd"), + pytest.param( + {"bsse_type": "cp", "return_total_data": False}, + "CP-CORRECTED INTERACTION ENERGY THROUGH 4-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and "TOTAL ENERGY" not in k)], + 15, + id="4b-cp"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "max_nbody": 3}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 3-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("4-BODY" not in k))], + 14, + id="3b-nocp-rtd"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": False, "max_nbody": 3}, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and "4-BODY" not in k)], + 14, + id="3b-nocp"), + pytest.param( + {"bsse_type": "cp", "return_total_data": True, "max_nbody": 3}, + "CP-CORRECTED TOTAL ENERGY THROUGH 3-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and "4-BODY" not in k)], + 18, # bugfix: was 28 + id="3b-cp-rtd"), + pytest.param( + {"bsse_type": "cp", "return_total_data": False, "max_nbody": 3}, + "CP-CORRECTED INTERACTION ENERGY THROUGH 3-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and "4-BODY" not in k and "TOTAL ENERGY" not in k)], + 14, + id="3b-cp"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "max_nbody": 2}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 2-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("4-BODY" not in k) and ("3-BODY" not in k))], + 10, + id="2b-nocp-rtd"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": False, "max_nbody": 2}, + "NOCP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("4-BODY" not in k) and ("3-BODY" not in k))], + 10, + id="2b-nocp"), + pytest.param( + {"bsse_type": "cp", "return_total_data": True, "max_nbody": 2}, + "CP-CORRECTED TOTAL ENERGY THROUGH 2-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and ("4-BODY" not in k) and ("3-BODY" not in k))], + 14, + id="2b-cp-rtd"), + pytest.param( + {"bsse_type": "cp", "return_total_data": False, "max_nbody": 2}, + "CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY", + [k for k in he4_refs_conv if (k.startswith("CP-") and ("4-BODY" not in k) and ("3-BODY" not in k) and "TOTAL ENERGY" not in k)], + 10, + id="2b-cp"), + pytest.param( + {"bsse_type": "nocp", "return_total_data": True, "max_nbody": 1}, + "NOCP-CORRECTED TOTAL ENERGY THROUGH 1-BODY", + [k for k in he4_refs_conv if (k.startswith("NOCP-") and ("1-BODY" in k))], + 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))], +# 10, +# 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))], + 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)], +# 4, +# id="1b-cp"), +]) +def test_nbody_he4_single(program, basis, keywords, mbe_keywords, anskey, bodykeys, calcinfo_nmbe, he_tetramer, request): + #! MP2/aug-cc-pvDZ many body energies of an arbitrary Helium complex, + # addressing 4-body formulas + # e, wfn = energy('MP2/aug-cc-pVDZ', molecule=he_tetramer, ...) + + atomic_spec = AtomicSpecification(model={"method": "mp2", "basis": basis}, program=program, driver="energy", keywords=keywords) + mbe_model = ManyBodyInput(specification={"specification": atomic_spec, "keywords": mbe_keywords, "driver": "energy"}, molecule=he_tetramer) + + if program == "gamess": + with pytest.raises(ValueError) as exe: + # qcng: qcng.compute_procedure(mbe_model, "manybody", raise_error=True) + qcng_compute_procedure(mbe_model) #, "manybody", raise_error=True) + assert "GAMESS+QCEngine can't handle ghost atoms yet" in str(exe.value) + pytest.xfail("GAMESS can't do ghosts") + + # qcng: ret = qcng.compute_procedure(mbe_model, "manybody", raise_error=True) + ret = ManyBodyComputerQCNG.from_qcschema_ben(mbe_model) + print(f"SSSSSSS {request.node.name}") + # v2: pprint.pprint(ret.model_dump(), width=200) + pprint.pprint(ret.dict(), width=200) + + _inner = request.node.name.split("[")[1].split("]")[0] + kwdsln, progln = "-".join(_inner.split("-")[:-2]), "-".join(_inner.split("-")[-2:]) + refs = he4_refs_df if progln == "psi4-df" else he4_refs_conv + ans = refs[anskey] + + for qcv, ref in refs.items(): + skp = skprop(qcv) + if qcv in bodykeys: + assert compare_values(ref, ret.extras["qcvars"]["nbody"][qcv], atol=1.0e-8, label=f"[a] qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=1.0e-8, 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["4b-all"]: + skp = skprop(qcv) + if qcv in sumdict[kwdsln]: + refkey = sumdict[kwdsln][qcv] + ref = 0.0 if refkey == "zero" else refs[refkey] + assert compare_values(ref, ret.extras["qcvars"]["nbody"][qcv], atol=1.0e-8, label=f"[c] qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=1.0e-8, 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": ans, + }.items(): + skp = skprop(qcv) + assert compare_values(ref, ret.extras["qcvars"][qcv], atol=1.0e-8, label=f"[e] qcvars {qcv}") + assert compare_values(ref, getattr(ret.properties, skp), atol=1.0e-8, label=f"[f] skprop {skp}") + assert compare_values(ans, ret.return_result, atol=1.0e-8, label=f"[g] ret") + assert ret.properties.calcinfo_nmbe == calcinfo_nmbe, f"{ret.properties.calcinfo_nmbe=} != {calcinfo_nmbe}" diff --git a/qcmanybody/models/test_mbe_keywords.py b/qcmanybody/models/test_mbe_keywords.py index 696314b..71ed8d1 100644 --- a/qcmanybody/models/test_mbe_keywords.py +++ b/qcmanybody/models/test_mbe_keywords.py @@ -237,6 +237,10 @@ def test_mbe_bsse_type(mbe_data, kws, ans): pytest.param({"supersystem_ie_only": False}, False), pytest.param({"max_nbody": 1, "supersystem_ie_only": True}, "Cannot skip intermediate n-body jobs when max_nbody=1 != nfragments"), + pytest.param({"bsse_type": "vmFC", "supersystem_ie_only": True}, + "Cannot skip intermediate n-body jobs when VMFC in bsse_type"), + pytest.param({"bsse_type": ["cp", "Vmfc"], "supersystem_ie_only": True}, + "Cannot skip intermediate n-body jobs when VMFC in bsse_type"), ]) def test_mbe_sie(mbe_data, kws, ans): mbe_data["specification"]["keywords"] = kws diff --git a/qcmanybody/utils.py b/qcmanybody/utils.py index 3293b89..be63afe 100644 --- a/qcmanybody/utils.py +++ b/qcmanybody/utils.py @@ -146,7 +146,7 @@ def print_nbody_energy( print(info) -def collect_vars(bsse, body_dict, max_nbody: int, embedding: bool = False): +def collect_vars(bsse, body_dict, max_nbody: int, embedding: bool = False, supersystem_ie_only: bool = False): previous_e = body_dict[1] tot_e = previous_e != 0.0 nbody_range = list(body_dict) @@ -159,11 +159,32 @@ def collect_vars(bsse, body_dict, max_nbody: int, embedding: bool = False): res[f"{bsse}-CORRECTED TOTAL ENERGY"] = body_dict[max_nbody] res[f"{bsse}-CORRECTED INTERACTION ENERGY"] = body_dict[max_nbody] - body_dict[1] - for nb in range(2, max(nbody_range) + 1): - res[f"{bsse}-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = body_dict[nb] - body_dict[1] - res[f"{bsse}-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = body_dict[nb] - body_dict[nb - 1] - if tot_e: - for nb in nbody_range: - res[f"{bsse}-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = body_dict[nb] + if supersystem_ie_only: + nfr = nbody_range[-1] + for nb in [nfr]: + res[f"{bsse}-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = body_dict[nb] - body_dict[1] + if nb == 2: + res[f"{bsse}-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = body_dict[nb] - body_dict[nb - 1] + if tot_e: + for nb in [1, nfr]: + res[f"{bsse}-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = body_dict[nb] + else: + for nb in range(2, max(nbody_range) + 1): + res[f"{bsse}-CORRECTED INTERACTION ENERGY THROUGH {nb}-BODY"] = body_dict[nb] - body_dict[1] + res[f"{bsse}-CORRECTED {nb}-BODY CONTRIBUTION TO ENERGY"] = body_dict[nb] - body_dict[nb - 1] + if tot_e: + for nb in nbody_range: + res[f"{bsse}-CORRECTED TOTAL ENERGY THROUGH {nb}-BODY"] = body_dict[nb] return res + + +def provenance_stamp(routine: str) -> Dict[str, str]: + """Return dictionary satisfying QCSchema, + https://github.com/MolSSI/QCSchema/blob/master/qcschema/dev/definitions.py#L23-L41 + with QCManyBody's credentials for creator and version. The + generating routine's name is passed in through `routine`. + + """ + import qcmanybody + return {"creator": "QCManyBody", "version": qcmanybody.__version__, "routine": routine} diff --git a/tests/test_multi.py b/tests/test_multi.py index 06099bf..67e9c8d 100644 --- a/tests/test_multi.py +++ b/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) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels) diff --git a/tests/test_multi_ss.py b/tests/test_multi_ss.py index 78f19da..872f35a 100644 --- a/tests/test_multi_ss.py +++ b/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) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels) diff --git a/tests/test_single.py b/tests/test_single.py index 97a714d..9e81748 100644 --- a/tests/test_single.py +++ b/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) + mc = ManyBodyCalculator(mol_h2o_3, [BsseEnum.cp, BsseEnum.nocp, BsseEnum.vmfc], levels, True, False) nbody_results = mc.analyze(component_results) compare_results(nbody_results, ref_data, levels)