adapt s-dftd3 harness for psi4 (#385)

* adapt s-dftd3 harness for psi4

* retrofit s-dftd3 adaptations to old psi

* tidying

* fix pairwise in classic dftd3

* fix sapt pairwise and format

* work on tests

docs/source/changelog.rst

@@ -17,20 +17,36 @@ Changelog
 .. +++++++++
 
 
-v0.26.0 / 2022-MM-DD (Unreleased)
+v0.26.0 / 2022-11-30
 --------------------
 
-.. Breaking Changes
-.. ++++++++++++++++
+Breaking Changes
+++++++++++++++++
 
-.. New Features
-.. ++++++++++++
+- (:pr:`385`) Dispersion - the dispersion parameters resources file has been altered so that for D3 variants there's a
+  2b set (e.g., d3bj2b) that is pure 2-body and doesn't accept s9 (effectively fixed at 0.0) and a atm set (e.g.,
+  d3zeroatm) that does accept s9 (by default 1.0 but user-variable). Previous D3 levels are aliased to 2b. Only
+  downstreams that call the dispersion resources directly should be affected, and retrofits are in place for the known
+  victim/instigator (Psi4). @loriab
 
-.. Enhancements
-.. ++++++++++++
+New Features
+++++++++++++
 
-.. Bug Fixes
-.. +++++++++
+Enhancements
+++++++++++++
+- (:pr:`380`) MRChem - added gradient and thus geometry optimizations support. @robertodr
+- (:pr:`385`) dftd3 - the classic interface now accepts e.g., ``d3mbj2b`` as a level hint. @loriab
+- (:pr:`385`) s-dftd3 - added keyword ``apply_qcengine_aliases`` that when True and ``level_hint`` present allows the
+  levels and aliases in the dispersion resources (e.g., ``d3``, ``d3atm``, ``d32b``) to be given as ``level_hint``. The
+  resource parameters are passed to s-dftd3 as param_tweaks. @loriab
+
+Bug Fixes
++++++++++
+- (:pr:`383`) yaml - uses safe loading. @mbanck, @loriab
+- (:pr:`385`) dftd3 - the pairwise analysis requested through ``AtomicInput.keywords["pair_resolved"] = True`` and
+  returned in ``AtomicResult.extras["qcvars"]["2-BODY PAIRWISE DISPERSION CORRECTION ANALYSIS"]`` was elementwise too
+  large by a factor of 2. It now matches the ``s-dftd3`` harness and fulfills that the sum of the array equals the
+  2-body dispersion energy. @loriab
 
 
 v0.25.0 / 2022-11-11
@@ -39,20 +55,20 @@ v0.25.0 / 2022-11-11
 Breaking Changes
 ++++++++++++++++
 - (:pr:`376`) GAMESS - slight breaking changes of (1) ROHF MP2 ZAPT quantities now stored in "ZAPT" variables, not "MP2"
-  variables; and (2) "HF TOTAL ENERGY" no longer stores DFT energy in DFT computation.
+  variables; and (2) "HF TOTAL ENERGY" no longer stores DFT energy in DFT computation. @loriab
 - (:pr:`376`) testing - reference quantities now indexed by "standard" or "semicanonical" orbitals since program defaults
   differ (mostly in CCSD ROHF FC). Downstream projects using the stdsuite interface will need to add an extra argument to query
-  reference data.
+  reference data. @loriab
 
 New Features
 ++++++++++++
 
 Enhancements
 ++++++++++++
-- (:pr:`376`) Cfour - added parsing for BCCD and BCCD(T) methods.
-- (:pr:`376`) NWChem - B2PLYP double-hybrid can now be run and parsed. Added CC2 parsing.
-- (:pr:`376`) testing - added parsing contracts for ZAPT2, CEPA(1), CEPA(3), ACPF, AQCC, BCCD, BCCD(T), CC2, CC3, and DH-DFT. Added conventional references for most.
-- (:pr:`378`) OpenFF - Support OpenFF Toolkit v0.11+
+- (:pr:`376`) Cfour - added parsing for BCCD and BCCD(T) methods. @loriab
+- (:pr:`376`) NWChem - B2PLYP double-hybrid can now be run and parsed. Added CC2 parsing. @loriab
+- (:pr:`376`) testing - added parsing contracts for ZAPT2, CEPA(1), CEPA(3), ACPF, AQCC, BCCD, BCCD(T), CC2, CC3, and DH-DFT. Added conventional references for most. @loriab
+- (:pr:`378`) OpenFF - Support OpenFF Toolkit v0.11+. @Yoshanuikabundi
 
 Bug Fixes
 +++++++++
@@ -63,7 +79,7 @@ v0.24.1 / 2022-08-16
 
 Enhancements
 ++++++++++++
-- (:pr:`375`) testing - in standard suite, add reference values for occ, a-ccsd(t), olccd grad, remp2, omp2, omp2.5, omp3, oremp2, density fitted ccsd, ccsd(t), a-ccsd(t).
+- (:pr:`375`) testing - in standard suite, add reference values for occ, a-ccsd(t), olccd grad, remp2, omp2, omp2.5, omp3, oremp2, density fitted ccsd, ccsd(t), a-ccsd(t). @loriab
 
 
 v0.24.0 / 2022-07-08

qcengine/programs/dftd3.py

@@ -201,7 +201,7 @@ def parse_output(self, outfiles: Dict[str, str], input_model: "AtomicInput") ->
                     atom1 = int(data[0]) - 1
                     atom2 = int(data[1]) - 1
                     Edisp = Decimal(data[-1])
-                    D3pairs[atom1, atom2] = Edisp / Decimal(Grimme_h2kcal)
+                    D3pairs[atom1, atom2] = Decimal(0.5) * Edisp / Decimal(Grimme_h2kcal)
                     D3pairs[atom2, atom1] = D3pairs[atom1, atom2]
 
             elif re.match(" normal termination of dftd3", ln):
@@ -297,6 +297,9 @@ def parse_output(self, outfiles: Dict[str, str], input_model: "AtomicInput") ->
         output_data["extras"]["local_keywords"] = input_model.extras["info"]
         output_data["extras"]["qcvars"] = calcinfo
         if input_model.keywords.get("pair_resolved", False):
+            assert (
+                abs(D3pairs.sum() - float(retres)) < 1.0e-6
+            ), f"pairwise sum {D3pairs.sum()} != energy {float(retres)}"
             output_data["extras"]["qcvars"]["2-BODY PAIRWISE DISPERSION CORRECTION ANALYSIS"] = D3pairs
         output_data["success"] = True
 
@@ -315,6 +318,8 @@ def dftd3_coeff_formatter(dashlvl: str, dashcoeff: Dict) -> str:
     d3mbj:   s6      a1       s8      a2      alpha6=None version=6
     atmgr:   s6=1.0  sr6=None s8=None a2=None alpha6      version=3 (needs -abc, too)
 
+    2-body variant here. that is, d3zero2b
+
     Parameters
     ----------
     dashlvl : {'d2', 'd3zero', d3bj', 'd3mzero', 'd3mbj', 'atmgr'}
@@ -339,15 +344,15 @@ def dftd3_coeff_formatter(dashlvl: str, dashcoeff: Dict) -> str:
     dashlvl = dashlvl.lower()
     if dashlvl == "d2":
         return dashformatter.format(dashcoeff["s6"], dashcoeff["sr6"], 0.0, 0.0, dashcoeff["alpha6"], 2)
-    elif dashlvl == "d3zero":
+    elif dashlvl == "d3zero2b":
         return dashformatter.format(
             dashcoeff["s6"], dashcoeff["sr6"], dashcoeff["s8"], dashcoeff["sr8"], dashcoeff["alpha6"], 3
         )
-    elif dashlvl == "d3bj":
+    elif dashlvl == "d3bj2b":
         return dashformatter.format(dashcoeff["s6"], dashcoeff["a1"], dashcoeff["s8"], dashcoeff["a2"], 0.0, 4)
-    elif dashlvl == "d3mzero":
+    elif dashlvl == "d3mzero2b":
         return dashformatter.format(dashcoeff["s6"], dashcoeff["sr6"], dashcoeff["s8"], dashcoeff["beta"], 14.0, 5)
-    elif dashlvl == "d3mbj":
+    elif dashlvl == "d3mbj2b":
         return dashformatter.format(dashcoeff["s6"], dashcoeff["a1"], dashcoeff["s8"], dashcoeff["a2"], 0.0, 6)
     elif dashlvl == "atmgr":
         # need to set first four parameters to something other than None, otherwise DFTD3 gets mad or a bit wrong

qcengine/programs/dftd_ng.py

@@ -242,7 +242,7 @@ def compute(self, input_model: AtomicInput, config: TaskConfig) -> AtomicResult:
             input_data["model"]["method"] = method
         qcvkey = method.upper() if method is not None else None
 
-        # send `from_arrays` the dftd3 behavior of functional specification overrides explicit parameters specification
+        # send `from_arrays` the s-dftd3 behavior of functional specification overrides explicit parameters specification
         # * differs from dftd3 harness behavior where parameters extend or override functional
         # * stash the resolved plan in extras or, if errored, leave it for the proper dftd3 api to reject
         param_tweaks = None if method else input_model.keywords.get("params_tweaks", None)
@@ -265,6 +265,18 @@ def compute(self, input_model: AtomicInput, config: TaskConfig) -> AtomicResult:
                 method = method[: -(len(alias) + 1)]
                 input_data["model"]["method"] = method
 
+        # consolidate dispersion level aliases
+        if input_model.keywords.pop("apply_qcengine_aliases", False):
+            level_hint = input_model.keywords.get("level_hint", None)
+            if level_hint:
+                level_hint = get_dispersion_aliases()[level_hint.lower()]
+                if level_hint.endswith("atm"):
+                    level_hint = level_hint[:-3]
+                if level_hint.endswith("2b"):
+                    level_hint = level_hint[:-2]
+                    input_data["keywords"]["params_tweaks"] = {**planinfo["dashparams"], "s9": 0.0}
+                input_data["keywords"]["level_hint"] = level_hint
+
         input_model = AtomicInput(**input_data)
 
         # Run the Harness