diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst index 86d48913b..a480b42fe 100644 --- a/docs/source/changelog.rst +++ b/docs/source/changelog.rst @@ -28,6 +28,38 @@ Changelog .. - UNSOLVED (:issue:`397`) extras failed +v0.31.0 / 2024-MM-DD (Unreleased) +-------------------- + +Breaking Changes +++++++++++++++++ + +New Features +++++++++++++ +- UNMERGED (:pr:`448`) QCManyBody - new procedure for computing interaction energies or truncation or full + many-body expansions. @loriab +- UNMERGED (:pr:`448`) GeomeTRIC, OptKing - new derived procedures running `GeneralizedOptimizationInput`/ + `Results` that allow using a `ManyBodySpecification` instead of `QCSpecification` to drive counterpoise- + corrected or partial-body optimizations. May require explicit use of `schema_name` field. Call with + `qcengine.compute_procedure(qcsk, "gengeometric")` or `genoptking`. This is temporary until QCSchema v2. @loriab + +Enhancements +++++++++++++ + +Bug Fixes ++++++++++ +- UNMERGED (:pr:`448`) CFour - fix error collecting molecule when it's a single atom with two-letter symbol. @loriab + +Misc. ++++++ + +MUST (Unmerged) ++++++++++++++++ + +WIP (Unmerged) +++++++++++++++ + + v0.30.0 / 2024-06-25 -------------------- diff --git a/qcengine/procedures/base.py b/qcengine/procedures/base.py index 2e9979d4a..0c2c23e0c 100644 --- a/qcengine/procedures/base.py +++ b/qcengine/procedures/base.py @@ -6,9 +6,10 @@ from ..exceptions import InputError, ResourceError from .berny import BernyProcedure -from .geometric import GeometricProcedure +from .geometric import GeometricProcedure, GenGeometricProcedure +from .qcmanybody import QCManyBodyProcedure from .nwchem_opt import NWChemDriverProcedure -from .optking import OptKingProcedure +from .optking import OptKingProcedure, GenOptKingProcedure from .torsiondrive import TorsionDriveProcedure from .model import ProcedureHarness @@ -67,7 +68,10 @@ def list_available_procedures() -> Set[str]: register_procedure(GeometricProcedure()) +register_procedure(GenGeometricProcedure()) register_procedure(OptKingProcedure()) +register_procedure(GenOptKingProcedure()) register_procedure(BernyProcedure()) +register_procedure(QCManyBodyProcedure()) register_procedure(NWChemDriverProcedure()) register_procedure(TorsionDriveProcedure()) diff --git a/qcengine/procedures/geometric.py b/qcengine/procedures/geometric.py index 59bc26c86..69a3e6bf2 100644 --- a/qcengine/procedures/geometric.py +++ b/qcengine/procedures/geometric.py @@ -68,3 +68,68 @@ def compute(self, input_model: "OptimizationInput", config: "TaskConfig") -> "Op output_data = OptimizationResult(**output_data) return output_data + + +class GenGeometricProcedure(GeometricProcedure): + + # note that "procedure" value below not used + _defaults = {"name": "GenGeometric", "procedure": "genoptimization"} + + version_cache: Dict[str, str] = {} + + def found(self, raise_error: bool = False) -> bool: + qc = which_import( + "geometric", + return_bool=True, + raise_error=raise_error, + raise_msg="Please install via `conda install geometric -c conda-forge`.", + ) + dep = which_import( + "qcmanybody", + return_bool=True, + raise_error=raise_error, + raise_msg="For GenGeometric harness, please install via `conda install qcmanybody -c conda-forge`.", + ) + + return qc and dep + + def build_input_model( + self, data: Union[Dict[str, Any], "GeneralizedOptimizationInput"] + ) -> "GeneralizedOptimizationInput": + from qcmanybody.models.generalized_optimization import GeneralizedOptimizationInput + + return self._build_model(data, GeneralizedOptimizationInput) + + def compute( + self, input_model: "GeneralizedOptimizationInput", config: "TaskConfig" + ) -> "GeneralizedOptimizationResult": + self.found(raise_error=True) + + import geometric + from qcmanybody.models.generalized_optimization import GeneralizedOptimizationResult + + input_data = input_model.dict() + + # Temporary patch for geomeTRIC + input_data["initial_molecule"]["symbols"] = list(input_data["initial_molecule"]["symbols"]) + + # Set retries to two if zero while respecting local_config + local_config = config.dict() + local_config["retries"] = local_config.get("retries", 2) or 2 + input_data["input_specification"]["extras"]["_qcengine_local_config"] = local_config + + # Run the program + output_data = geometric.run_json.geometric_run_json(input_data) + + output_data["provenance"] = { + "creator": "geomeTRIC", + "routine": "geometric.run_json.geometric_run_json", + "version": geometric.__version__, + } + + output_data["schema_name"] = "qcschema_generalizedoptimizationresult" + output_data["input_specification"]["extras"].pop("_qcengine_local_config", None) + if output_data["success"]: + output_data = GeneralizedOptimizationResult(**output_data) + + return output_data diff --git a/qcengine/procedures/optking.py b/qcengine/procedures/optking.py index 34139a3f3..09ef41e36 100644 --- a/qcengine/procedures/optking.py +++ b/qcengine/procedures/optking.py @@ -1,10 +1,14 @@ -from typing import Any, Dict, Union +from typing import TYPE_CHECKING, Any, Dict, Union from qcelemental.models import OptimizationInput, OptimizationResult from qcelemental.util import safe_version, which_import from .model import ProcedureHarness +if TYPE_CHECKING: + from qcengine.config import TaskConfig + from qcmanybody.models.generalized_optimization import GeneralizedOptimizationInput, GeneralizedOptimizationResult + class OptKingProcedure(ProcedureHarness): @@ -37,7 +41,7 @@ def get_version(self) -> str: return self.version_cache[which_prog] - def compute(self, input_model: "OptimizationInput", config: "TaskConfig") -> "Optimization": + def compute(self, input_model: "OptimizationInput", config: "TaskConfig") -> "OptimizationResult": if self.found(raise_error=True): import optking @@ -57,3 +61,59 @@ def compute(self, input_model: "OptimizationInput", config: "TaskConfig") -> "Op output_data = OptimizationResult(**output_data) return output_data + + +class GenOptKingProcedure(OptKingProcedure): + + # note that "procedure" value below not used + _defaults = {"name": "GenOptKing", "procedure": "genoptimization"} + + version_cache: Dict[str, str] = {} + + def found(self, raise_error: bool = False) -> bool: + qc = which_import( + "optking", + return_bool=True, + raise_error=raise_error, + raise_msg="Please install via `conda install optking -c conda-forge`.", + ) + dep = which_import( + "qcmanybody", + return_bool=True, + raise_error=raise_error, + raise_msg="For GenOptKing harness, please install via `conda install qcmanybody -c conda-forge`.", + ) + + return qc and dep + + def build_input_model( + self, data: Union[Dict[str, Any], "GeneralizedOptimizationInput"] + ) -> "GeneralizedOptimizationInput": + from qcmanybody.models.generalized_optimization import GeneralizedOptimizationInput + + return self._build_model(data, GeneralizedOptimizationInput) + + def compute( + self, input_model: "GeneralizedOptimizationInput", config: "TaskConfig" + ) -> "GeneralizedOptimizationResult": + self.found(raise_error=True) + + import optking + from qcmanybody.models.generalized_optimization import GeneralizedOptimizationResult + + input_data = input_model.dict() + + # Set retries to two if zero while respecting local_config + local_config = config.dict() + local_config["retries"] = local_config.get("retries", 2) or 2 + input_data["input_specification"]["extras"]["_qcengine_local_config"] = local_config + + # Run the program + output_data = optking.optimize_qcengine(input_data) + + output_data["schema_name"] = "qcschema_generalizedoptimizationresult" + output_data["input_specification"]["extras"].pop("_qcengine_local_config", None) + if output_data["success"]: + output_data = GeneralizedOptimizationResult(**output_data) + + return output_data diff --git a/qcengine/procedures/qcmanybody.py b/qcengine/procedures/qcmanybody.py new file mode 100644 index 000000000..b663bb242 --- /dev/null +++ b/qcengine/procedures/qcmanybody.py @@ -0,0 +1,51 @@ +from typing import TYPE_CHECKING, Any, Dict, Union + +from qcelemental.util import safe_version, which_import + +from .model import ProcedureHarness + +if TYPE_CHECKING: + from ..config import TaskConfig + from qcmanybody.models import ManyBodyInput, ManyBodyResult + + +class QCManyBodyProcedure(ProcedureHarness): + + # v2: ClassVar[Dict[str, Any]] + _defaults: Dict[str, Any] = {"name": "QCManyBody", "procedure": "manybody"} + + version_cache: Dict[str, str] = {} + + class Config(ProcedureHarness.Config): + pass + + def found(self, raise_error: bool = False) -> bool: + return which_import( + "qcmanybody", + return_bool=True, + raise_error=raise_error, + raise_msg="Please install via `conda install qcmanybody -c conda-forge`.", + ) + + def build_input_model(self, data: Union[Dict[str, Any], "ManyBodyInput"]) -> "ManyBodyInput": + from qcmanybody.models import ManyBodyInput + + return self._build_model(data, ManyBodyInput) + + def get_version(self) -> str: + self.found(raise_error=True) + + which_prog = which_import("qcmanybody") + if which_prog not in self.version_cache: + import qcmanybody + + self.version_cache[which_prog] = safe_version(qcmanybody.__version__) + + return self.version_cache[which_prog] + + def compute(self, input_model: "ManyBodyInput", config: "TaskConfig") -> "ManyBodyResult": + from qcmanybody import ManyBodyComputer + + output_model = ManyBodyComputer.from_manybodyinput(input_model) + + return output_model diff --git a/qcengine/programs/adcc.py b/qcengine/programs/adcc.py index ff641c5ef..c3fefa334 100644 --- a/qcengine/programs/adcc.py +++ b/qcengine/programs/adcc.py @@ -56,7 +56,7 @@ def found(raise_error: bool = False) -> bool: "psi4", return_bool=True, raise_error=raise_error, - raise_msg="Please install psi4 for adcc harness via `conda install psi4 -c conda-forge/label/libint_dev -c conda-forge`.", + raise_msg="Please install psi4 for adcc harness via `conda install psi4 -c conda-forge`.", ) return found_adcc and found_psi4 diff --git a/qcengine/programs/cfour/harvester.py b/qcengine/programs/cfour/harvester.py index 773d6d3d3..8e436e386 100644 --- a/qcengine/programs/cfour/harvester.py +++ b/qcengine/programs/cfour/harvester.py @@ -1004,7 +1004,9 @@ def harvest_outfile_pass(outtext): # Process atom geometry mobj = re.search(r"^\s+" + r"@GETXYZ-I, 1 atoms read from ZMAT." + r"\s*$", outtext, re.MULTILINE) mobj2 = re.search( - r"^([A-Z]+)#1" + r"\s+" + NUMBER + r"\s+" + NUMBER + r"\s+" + NUMBER + r"\s*$", outtext, re.MULTILINE + r"^\s*([A-Z]+)" + r"\s*" + r"#1" + r"\s+" + NUMBER + r"\s+" + NUMBER + r"\s+" + NUMBER + r"\s*$", + outtext, + re.MULTILINE, ) if mobj and mobj2: logger.debug("matched atom2") # unsavory for when atom never printed except for basis file diff --git a/qcengine/programs/psi4.py b/qcengine/programs/psi4.py index 1ba3d7fac..97aa88056 100644 --- a/qcengine/programs/psi4.py +++ b/qcengine/programs/psi4.py @@ -92,7 +92,7 @@ def found(raise_error: bool = False) -> bool: "psi4", return_bool=True, raise_error=raise_error, - raise_msg="Please install via `conda install psi4 -c conda-forge/label/libint_dev -c conda-forge`. Check it's in your PATH with `which psi4`." + raise_msg="Please install via `conda install psi4 -c conda-forge`. Check it's in your PATH with `which psi4`." + error_msg, ) diff --git a/qcengine/testing.py b/qcengine/testing.py index c7a3ecf11..e6dafb423 100644 --- a/qcengine/testing.py +++ b/qcengine/testing.py @@ -57,9 +57,14 @@ def is_program_new_enough(program, version_feature_introduced): form is equal to or later than `version_feature_introduced`. """ - if program not in qcng.list_available_programs(): - return False - candidate_version = qcng.get_program(program).get_version() + if program in qcng.list_all_procedures(): + if program not in qcng.list_available_procedures(): + return False + candidate_version = qcng.get_procedure(program).get_version() + else: + if program not in qcng.list_available_programs(): + return False + candidate_version = qcng.get_program(program).get_version() return parse_version(candidate_version) >= parse_version(version_feature_introduced) @@ -164,6 +169,7 @@ def get_job(self): "mctc-gcp": is_program_new_enough("mctc-gcp", "2.3.0"), "gcp": which("gcp", return_bool=True), "geometric": which_import("geometric", return_bool=True), + "geometric_genopt": is_program_new_enough("gengeometric", "1.1.0"), "berny": which_import("berny", return_bool=True), "mdi": is_mdi_new_enough("1.2"), "molpro": is_program_new_enough("molpro", "2018.1"), @@ -171,12 +177,14 @@ def get_job(self): "mp2d": which("mp2d", return_bool=True), "nwchem": which("nwchem", return_bool=True), "optking": which_import("optking", return_bool=True), + "optking_genopt": is_program_new_enough("genoptking", "0.3.0"), "psi4": is_program_new_enough("psi4", "1.2"), "psi4_runqcsk": is_program_new_enough("psi4", "1.4a2.dev160"), "psi4_mp2qcsk": is_program_new_enough("psi4", "1.4a2.dev580"), "psi4_derqcsk": is_program_new_enough("psi4", "1.5a1.dev117"), "qcdb": which_import("qcdb", return_bool=True), "qchem": is_program_new_enough("qchem", "5.1"), + "qcmanybody": which_import("qcmanybody", return_bool=True), "rdkit": which_import("rdkit", return_bool=True), "terachem": which("terachem", return_bool=True), "terachem_pbs": is_program_new_enough("terachem_pbs", "0.7.2"), diff --git a/qcengine/tests/test_qcmanybody.py b/qcengine/tests/test_qcmanybody.py new file mode 100644 index 000000000..ead2354ab --- /dev/null +++ b/qcengine/tests/test_qcmanybody.py @@ -0,0 +1,448 @@ +import pprint + +import pytest + +from qcelemental import constants +from qcelemental.models import Molecule +from qcelemental.testing import compare, compare_values + +import qcengine as qcng +from qcengine.testing import using + + +@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], + ) + + +@using("qcmanybody") +@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.0e-10, "d_convergence": 1.0e-10, "scf_type": "pk", "mp2_type": "conv"}, + id="psi4_conv", + marks=using("psi4"), + ), + ], +) +@pytest.mark.parametrize( + "mbe_keywords,anskey,calcinfo_nmbe", + [ + pytest.param( + {"bsse_type": "nocp", "return_total_data": False, "max_nbody": 3}, + "nocp_corrected_interaction_energy_through_3_body", + 14, + id="3b_nocp", + ), + ], +) +def test_nbody_he4_single(program, basis, keywords, mbe_keywords, anskey, calcinfo_nmbe, he_tetramer, request): + from qcmanybody.models import AtomicSpecification, ManyBodyInput + + atomic_spec = AtomicSpecification( + model={"method": "mp2", "basis": basis}, + program=program, + driver="energy", + keywords=keywords, + protocols={"stdout": False}, + ) + mbe_model = ManyBodyInput( + specification={ + "specification": atomic_spec, + "keywords": mbe_keywords, + "driver": "energy", + "protocols": {"component_results": "all"}, + }, + molecule=he_tetramer, + ) + + ret = qcng.compute_procedure(mbe_model, "qcmanybody", raise_error=True) + print(f"SSSSSSS {request.node.name}") + pprint.pprint(ret.dict(), width=200) + + assert ret.extras == {}, f"[w] extras wrongly present: {ret.extras.keys()}" + + refs = { + "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_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, + "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, + } + + ans = refs[anskey] + ref_nmbe = calcinfo_nmbe + atol = 1.0e-8 + + for skp, ref in refs.items(): + if not "4_body" in skp: + assert compare_values(ref, getattr(ret.properties, skp), atol=atol, label=f"[b] skprop {skp}") + else: + assert getattr(ret.properties, skp) is None + + assert compare_values( + refs["nocp_corrected_total_energy_through_3_body"], + ret.properties.nocp_corrected_total_energy, + atol=atol, + label=f"[d] skprop tot", + ) + assert compare_values( + refs["nocp_corrected_interaction_energy_through_3_body"], + ret.properties.nocp_corrected_interaction_energy, + atol=atol, + label=f"[d] skprop IE", + ) + + assert compare_values(ans, ret.properties.return_energy, atol=atol, label=f"[f] skprop {skp}") + assert compare_values(ans, ret.return_result, atol=atol, label=f"[g] ret") + + assert ret.properties.calcinfo_nmbe == ref_nmbe, f"[i] {ret.properties.calcinfo_nmbe} != {ref_nmbe}" + assert ( + len(ret.component_results) == ref_nmbe + ), f"[k] {len(ret.component_results)} != {ref_nmbe}; mbe protocol did not take" + if ref_nmbe > 0: + an_atres = next(iter(ret.component_results.values())) + assert an_atres.stdout is None, f"[l] atomic protocol did not take" + + +@using("qcmanybody") +@pytest.mark.parametrize( + "qcprog", + [ + pytest.param("cfour", marks=using("cfour")), + pytest.param("gamess", marks=using("gamess")), + pytest.param("nwchem", marks=using("nwchem")), + pytest.param("psi4", marks=using("psi4")), + ], +) +def test_bsse_ene_tu6_cp_ne2(qcprog): + """ + from https://github.com/psi4/psi4/blob/master/tests/tu6-cp-ne2/input.dat + Example potential energy surface scan and CP-correction for Ne2 + """ + from qcmanybody.models import ManyBodyInput + + tu6_ie_scan = {2.5: 0.757717, 3.0: 0.015685, 4.0: -0.016266} # Ang: kcal/mol IE + + keywords = { + "cfour": {"frozen_core": True}, + "gamess": {"contrl__ispher": 1}, + "nwchem": {"ccsd__freeze__atomic": True, "basis__spherical": True}, + "psi4": {"freeze_core": True}, + } + basis = { + "cfour": "aug-pvdz", + "gamess": "accd", + "nwchem": "aug-cc-pvdz", + "psi4": "aug-cc-pvdz", + } + + mbe_data = { + "specification": { + "specification": { + "model": { + "method": "ccsd(t)", + "basis": basis[qcprog], + }, + "driver": "energy", + "program": qcprog, + "keywords": keywords[qcprog], + }, + "keywords": { + "bsse_type": "cp", + }, + "driver": "energy", + }, + "molecule": None, + } + + for R in tu6_ie_scan: + # TODO fix_symmetry='c1' propagate + nene = Molecule( + symbols=["Ne", "Ne"], fragments=[[0], [1]], geometry=[0, 0, 0, 0, 0, R / constants.bohr2angstroms] + ) + mbe_data["molecule"] = nene + + mbe_model = ManyBodyInput(**mbe_data) + print("IIIIIII") + pprint.pprint(mbe_model.dict(), width=200) + if qcprog == "gamess": + with pytest.raises(RuntimeError) as exe: + qcng.compute_procedure(mbe_model, "qcmanybody", raise_error=True) + assert "GAMESS+QCEngine can't handle ghost atoms yet" in str(exe.value) + pytest.xfail("GAMESS can't do ghosts") + + ret = qcng.compute_procedure(mbe_model, "qcmanybody", raise_error=True) + # print("SSSSSSS") + # pprint.pprint(ret.dict(), width=200) + + assert compare_values( + tu6_ie_scan[R], ret.return_result * constants.hartree2kcalmol, atol=1.0e-4, label=f"CP-CCSD(T) [{R:3.1f}]" + ) + assert compare(3, ret.properties.calcinfo_nmbe, label="nmbe") + assert compare(1, ret.properties.calcinfo_nmc, label="nmc") + assert compare(2, ret.properties.calcinfo_nfr, label="nfr") + assert compare(2, ret.properties.calcinfo_natom, label="nat") + + +@using("qcmanybody") +def test_mbe_error(): + from qcmanybody.models import ManyBodyInput + + mbe_data = { + "specification": { + "specification": { + "model": { + "method": "nonsense", + "basis": "nonsense", + }, + "driver": "energy", + "program": "cms", + }, + "keywords": {}, + "driver": "energy", + }, + "molecule": None, + } + + nene = Molecule( + symbols=["Ne", "Ne"], fragments=[[0], [1]], geometry=[0, 0, 0, 0, 0, 3.0 / constants.bohr2angstroms] + ) + mbe_data["molecule"] = nene + + mbe_model = ManyBodyInput(**mbe_data) + + # test 1 + with pytest.raises(RuntimeError) as exc: + qcng.compute_procedure(mbe_model, "qcmanybody", raise_error=True) + + assert "Program cms is not registered to QCEngine" in str(exc.value) + + # test 2 + ret = qcng.compute_procedure(mbe_model, "qcmanybody") + assert ret.success is False + assert "Program cms is not registered to QCEngine" in ret.error.error_message + + +@using("psi4") +@using("qcmanybody") +@pytest.mark.parametrize( + "optimizer,bsse_type,sio", + [ + pytest.param("optking", "none", None, marks=using("optking")), + pytest.param("genoptking", "none", None, marks=using("optking_genopt")), + pytest.param("genoptking", "nocp", True, marks=using("optking_genopt")), + pytest.param("genoptking", "cp", False, marks=using("optking_genopt")), + pytest.param("geometric", "none", None, marks=using("geometric")), + pytest.param("gengeometric", "none", None, marks=using("geometric_genopt")), + pytest.param("gengeometric", "nocp", False, marks=using("geometric_genopt")), + pytest.param("gengeometric", "cp", True, marks=using("geometric_genopt")), + ], +) +def test_bsse_opt_hf_trimer(optimizer, bsse_type, sio): + + initial_molecule = Molecule.from_data( + """ +F -0.04288 2.78905 0.00000 +H 0.59079 2.03435 0.00000 +-- +F -1.94320 -0.70822 0.00000 +H -1.60642 0.21789 -0.00000 +-- +F 2.03569 -0.60531 -0.00000 +H 1.06527 -0.77673 0.00000 +units ang +""" + ) + + at_spec = { + # schema_name needed for differentiation in genopt + "schema_name": "qcschema_input", + "model": { + "method": "hf", + "basis": "6-31g", + }, + "keywords": { + "scf_type": "df", + }, + } + + mbe_spec = { + # schema_name needed for differentiation in genopt + "schema_name": "qcschema_manybodyspecification", + "specification": { + "model": { + "method": "hf", + "basis": "6-31g", + }, + "driver": "energy", + "program": "psi4", + "keywords": {}, + "protocols": { + "stdout": False, + }, + "extras": { + "psiapi": True, + }, + }, + "keywords": { + "bsse_type": bsse_type, + "supersystem_ie_only": sio, + }, + "driver": "energy", + "protocols": { + "component_results": "all", + }, + } + + opt_data = { + "initial_molecule": initial_molecule, + "input_specification": at_spec if (bsse_type == "none") else mbe_spec, + "keywords": { + "program": "psi4", + "g_convergence": "nwchem_loose", + }, + "protocols": { + "trajectory": "initial_and_final", + }, + } + # from qcmanybody.models.generalized_optimization import GeneralizedOptimizationInput + # opt_data = GeneralizedOptimizationInput(**opt_data) + + ret = qcng.compute_procedure(opt_data, optimizer, raise_error=True) + + print("FFFFFFFFFF") + pprint.pprint(ret.dict(), width=200) + + r_fh_hb_xptd = { + "none": 2.18 / constants.bohr2angstroms, + "nocp": 2.18 / constants.bohr2angstroms, + "cp": 2.27 / constants.bohr2angstroms, + }[bsse_type] + r_fh_computed = ret.final_molecule.measure([1, 3]) + assert ( + pytest.approx(r_fh_computed, 1.0e-2) == r_fh_hb_xptd + ), f"hydrogen bond length computed ({r_fh_computed}) != expected ({r_fh_hb_xptd})" + assert ( + len(ret.trajectory) == 2 + ), f"trajectory protocol did not take. len(ret.trajectory)={len(ret.trajectory)} != 2 (initial_and_final)" + if bsse_type != "none": + xptd_nmbe = { + ("nocp", False): 7, + ("nocp", True): 4, + ("cp", False): 10, + ("cp", True): 7, + }[(bsse_type, sio)] + assert xptd_nmbe == len(ret.trajectory[-1].component_results), f"mbe protocol did not take" + assert ( + ret.trajectory[-1].component_results['["(auto)", [1, 2, 3], [1, 2, 3]]'].stdout is None + ), f"atomic protocol did not take" + + +@using("qcmanybody") +@pytest.mark.parametrize("bsse_type", ["mbe", "ssfc"]) # aka nocp, cp +@pytest.mark.parametrize( + "qcprog,qc_keywords", + [ + pytest.param("psi4", {}, marks=using("psi4")), + pytest.param("cfour", {}, marks=using("cfour")), + pytest.param("nwchem", {}, marks=using("nwchem")), + ], +) +@pytest.mark.parametrize( + "optimizer,opt_keywords", + [ + pytest.param("optking", {}, marks=using("optking_genopt")), + pytest.param( + "geometric", {"convergence_set": "interfrag_tight", "maxiter": 15}, marks=using("geometric_genopt") + ), + ], +) +def test_bsse_opt_lif_dimer(optimizer, opt_keywords, bsse_type, qcprog, qc_keywords): + # in geomeTRIC: test_lif_bsse + + lif = { + "symbols": ["Li", "F"], + "geometry": [0, 0, 0, 0, 0, 3], + "fragments": [[0], [1]], + "fragment_charges": [+1, -1], + } + + mbe_spec = { + "schema_name": "qcschema_manybodyspecification", + "specification": { + "model": { + "method": "hf", + "basis": "6-31G", + }, + "driver": "energy", + "program": qcprog, + "keywords": qc_keywords, + "protocols": { + "stdout": False, + }, + }, + "keywords": { + "bsse_type": bsse_type, + "supersystem_ie_only": True, + }, + "protocols": { + "component_results": "all", + }, + "driver": "energy", + } + + opt_data = { + "initial_molecule": lif, + "input_specification": mbe_spec, + "keywords": { + "program": "nonsense", + **opt_keywords, + }, + "protocols": { + "trajectory": "final", + }, + } + + ret = qcng.compute_procedure(opt_data, "gen" + optimizer, raise_error=True) + + # printing will show up if job fails + pprint.pprint(ret.dict(), width=200) + + assert ret.success + + assert len(ret.trajectory[0].component_properties) == (5 if bsse_type == "ssfc" else 3) + + assert ret.provenance.creator.lower() == optimizer + assert ret.trajectory[0].provenance.creator == "QCManyBody" + atres = list(ret.trajectory[0].component_results.values())[0] + assert atres.provenance.creator.lower() == qcprog + + Rlif = ret.final_molecule.measure([0, 1]) + Rref = 3.016 if bsse_type == "ssfc" else 2.969 + assert compare_values(Rlif, Rref, "bond length", atol=1.0e-3)