diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 71f3c62c..920a5ce3 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -163,12 +163,12 @@ jobs:
         # note: any activate/deactivate use the conda cmd. other cmds use mamba cmd.
 
     - name: Special Config - NWChem
-      if: "(matrix.cfg.label == 'NWChem70')"
+      if: matrix.cfg.label == 'NWChem70'
       run: |
         sudo apt-get -y install nwchem
 
     - name: Special Config - QCore
-      if: "(matrix.cfg.label == 'QCore')"
+      if: matrix.cfg.label == 'QCore'
       run: |
         qcore --accept-license
 
@@ -180,18 +180,18 @@ jobs:
       #if: false
       run: |
         conda remove qcelemental --force
-        python -m pip install 'git+https://github.com/loriab/QCElemental.git@csse_layout_536d' --no-deps
+        python -m pip install 'git+https://github.com/loriab/QCElemental.git@csse_layout_536f_rb2025_take2' --no-deps
 
       # note: conda remove --force, not mamba remove --force b/c https://github.com/mamba-org/mamba/issues/412
       #       alt. is micromamba but not yet ready for setup-miniconda https://github.com/conda-incubator/setup-miniconda/issues/75
     - name: Special Config - QCEngine Dep
-      if: "(startsWith(matrix.cfg.label, 'Psi4')) || (matrix.cfg.label == 'ADCC') || (matrix.cfg.label == 'optimization-dispersion')"
+      if: (startsWith(matrix.cfg.label, 'Psi4')) || (matrix.cfg.label == 'ADCC') || (matrix.cfg.label == 'optimization-dispersion')
       run: |
         conda remove qcengine --force
 
       # QCEngine CI and Psi4 are circularly dependent, so a hack is in order
     - name: Special Config - Faux Pydantic Upgrade
-      if: "((matrix.cfg.label == 'Psi4-1.6') || (matrix.cfg.label == 'optimization-dispersion')) && (runner.os != 'Windows')"
+      if: ((matrix.cfg.label == 'Psi4-1.6') || (matrix.cfg.label == 'optimization-dispersion')) && (runner.os != 'Windows')
       run: |
         sed -i s/from\ pydantic\ /from\ pydantic.v1\ /g ${CONDA_PREFIX}/lib/python${{ matrix.cfg.python-version }}/site-packages/psi4/driver/*py
 
@@ -203,7 +203,7 @@ jobs:
         python -m pip install 'git+https://github.com/MolSSI/QCElemental.git@next2025' --no-deps
 
     - name: Special Config - Forced Interface Upgrade
-      if: "(matrix.cfg.label == 'Psi4-1.6')"
+      if: matrix.cfg.label == 'Psi4-1.6'
       run: |
         grep -r "local_options" ${CONDA_PREFIX}/lib/python${{ matrix.cfg.python-version }}/site-packages/psi4/driver/
         sed -i "s/local_options/task_config/g" ${CONDA_PREFIX}/lib/python${{ matrix.cfg.python-version }}/site-packages/psi4/driver/procrouting/*py
@@ -221,7 +221,7 @@ jobs:
         git describe
 
     - name: QCEngineRecords
-      if: "(matrix.cfg.label != 'Psi4-1.6')"
+      if: matrix.cfg.label != 'Psi4-1.6'
       run: |
         qcengine info
         export QCER_VER=`python -c "import qcengine.testing; print(qcengine.testing.QCENGINE_RECORDS_COMMIT)"`
@@ -269,7 +269,7 @@ jobs:
       #if: false
       run: |
         conda remove qcelemental --force
-        python -m pip install 'git+https://github.com/loriab/QCElemental.git@csse_layout_536d' --no-deps
+        python -m pip install 'git+https://github.com/loriab/QCElemental.git@csse_layout_536f_rb2025_take2' --no-deps
 
     - name: Environment Information
       run: |
diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index 5d6b3113..1302a63e 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -70,10 +70,16 @@ Misc.
 
 MUST (Unmerged)
 +++++++++++++++
-torsiondrive rewritten in v2.
-berny harness rewritten in v2. optking and geometric natively speak v1, so adapted as well as can be.
-allow nwchemdriver w/o driver=energy. provenance now nwchemdriver not nwchemrelax
- Optking now fills in ``v2.OptimizationResult.stdout``. Through v2, once can alter gradient protocols in an optimization.
+- (:pr:`462`) 
+- (:pr:`462`) adapt harnesses for TD.initial_molecules -> TD.initial_molecule and TD.optimization_history -> TD.scan_results
+- (:pr:`462`) rdkit, store ``AtomicResult.properties.return_gradient`` and ``calcinfo_natom``. mrchem, store ``AtomicResult.properties.return_gradient``.
+- (:pr:`462`) lightly Adapt harnesses for Mol v3
+- (:pr:`462`) Use packaging instead of setuptools to provide version parsing
+- (:pr:`462`) torsiondrive now accepts protocols. use ``protocols={"scan_results": "all"}`` if going to be converted to v1.
+- (:pr:`461`) torsiondrive rewritten in v2.
+- (:pr:`461`) berny harness rewritten in v2. optking and geometric natively speak v1, so adapted as well as can be.
+- (:pr:`461`) allow nwchemdriver w/o driver=energy. provenance now nwchemdriver not nwchemrelax
+- (:pr:`461`)  Optking now fills in ``v2.OptimizationResult.stdout``. Through v2, once can alter gradient protocols in an optimization.
 - (:pr:`460`) integrate ``AtomicInput.specification`` into harnesses and show what new inputs look like in tests
 - (:pr:`459`) gcp, mp2d several got properties.return_energy, retunr_gradient
 - (:pr:`460`) If you're missing something from AtomicResult.extras, check AtomicResult.input_data.extras in case it was passed in on input
diff --git a/qcengine/procedures/berny.py b/qcengine/procedures/berny.py
index eb1ca4a9..5d032676 100644
--- a/qcengine/procedures/berny.py
+++ b/qcengine/procedures/berny.py
@@ -5,7 +5,7 @@
 from typing import Any, ClassVar, Dict, Union
 
 import numpy as np
-from qcelemental.models.v2 import FailedOperation, OptimizationInput, OptimizationResult
+from qcelemental.models.v2 import FailedOperation, Molecule, OptimizationInput, OptimizationResult
 from qcelemental.util import which_import
 
 import qcengine
@@ -85,11 +85,14 @@ def compute(
         except Exception:
             error = {"error_type": "unknown", "error_message": f"Berny error:\n{traceback.format_exc()}"}
         else:
+            final_molecule = trajectory[-1]["molecule"]
             output = {
                 "input_data": input_model,
-                "final_molecule": trajectory[-1]["molecule"],
+                "final_molecule": final_molecule,
                 "properties": {
+                    "nuclear_repulsion_energy": Molecule(**final_molecule).nuclear_repulsion_energy(),
                     "return_energy": trajectory[-1]["properties"]["return_energy"],
+                    "return_gradient": trajectory[-1]["properties"]["return_gradient"],
                     "optimization_iterations": len(trajectory),
                 },
                 "trajectory_results": trajectory,
diff --git a/qcengine/procedures/optking.py b/qcengine/procedures/optking.py
index f312307f..a2a28a9f 100644
--- a/qcengine/procedures/optking.py
+++ b/qcengine/procedures/optking.py
@@ -1,7 +1,7 @@
 import logging
 import sys
 from io import StringIO
-from typing import Any, ClassVar, Dict, Union
+from typing import TYPE_CHECKING, Any, ClassVar, Dict, Union
 
 from qcelemental.models.v1 import OptimizationResult
 from qcelemental.models.v2 import OptimizationInput
@@ -9,6 +9,9 @@
 
 from .model import ProcedureHarness
 
+if TYPE_CHECKING:
+    from ..config import TaskConfig
+
 
 class OptKingProcedure(ProcedureHarness):
 
diff --git a/qcengine/procedures/torsiondrive.py b/qcengine/procedures/torsiondrive.py
index 33439cde..bb50b225 100644
--- a/qcengine/procedures/torsiondrive.py
+++ b/qcengine/procedures/torsiondrive.py
@@ -54,8 +54,8 @@ def _compute(self, input_model: "TorsionDriveInput", config: "TaskConfig"):
         state = torsiondrive.td_api.create_initial_state(
             dihedrals=dihedrals,
             grid_spacing=grid_spacing,
-            elements=input_model.initial_molecules[0].symbols,
-            init_coords=[molecule.geometry.flatten().tolist() for molecule in input_model.initial_molecules],
+            elements=input_model.initial_molecule[0].symbols,
+            init_coords=[molecule.geometry.flatten().tolist() for molecule in input_model.initial_molecule],
             dihedral_ranges=dihedral_ranges,
             energy_upper_limit=energy_upper_limit,
             energy_decrease_thresh=energy_decrease_thresh,
@@ -128,7 +128,7 @@ def _compute(self, input_model: "TorsionDriveInput", config: "TaskConfig"):
             output_data["final_energies"][grid_point] = final_energy
             output_data["final_molecules"][grid_point] = final_molecule
 
-        output_data["optimization_history"] = optimization_results
+        output_data["scan_results"] = optimization_results
 
         if error is not None:
             output_data["error"] = error
@@ -185,7 +185,7 @@ def _spawn_optimization(
 
         from qcengine import compute
 
-        input_molecule = input_model.initial_molecules[0].copy(deep=True).dict()
+        input_molecule = input_model.initial_molecule[0].copy(deep=True).dict()
         input_molecule["geometry"] = np.array(job).reshape(len(input_molecule["symbols"]), 3)
         input_molecule = Molecule.from_data(input_molecule)
 
diff --git a/qcengine/programs/cfour/harvester.py b/qcengine/programs/cfour/harvester.py
index cf6e9987..45a9d7d1 100644
--- a/qcengine/programs/cfour/harvester.py
+++ b/qcengine/programs/cfour/harvester.py
@@ -997,7 +997,7 @@ def harvest_outfile_pass(outtext):
         psivar_coord = Molecule(
             validate=False,
             **qcel.molparse.to_schema(
-                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
+                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=3
             ),
         )
 
@@ -1015,7 +1015,7 @@ def harvest_outfile_pass(outtext):
         psivar_coord = Molecule(
             validate=False,
             **qcel.molparse.to_schema(
-                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
+                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=3
             ),
         )
 
@@ -1034,7 +1034,7 @@ def harvest_outfile_pass(outtext):
         psivar_coord = Molecule(
             validate=False,
             **qcel.molparse.to_schema(
-                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
+                qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=3
             ),
         )
 
@@ -1374,7 +1374,7 @@ def harvest_GRD(grd):
     mol = Molecule(
         validate=False,
         **qcel.molparse.to_schema(
-            qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
+            qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=3
         ),
     )
 
diff --git a/qcengine/programs/gamess/harvester.py b/qcengine/programs/gamess/harvester.py
index cc7bbb7f..043eab77 100644
--- a/qcengine/programs/gamess/harvester.py
+++ b/qcengine/programs/gamess/harvester.py
@@ -544,7 +544,7 @@ def harvest_outfile_pass(outtext):
             qcvar_coord = Molecule(
                 validate=False,
                 **qcel.molparse.to_schema(
-                    qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=2
+                    qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"], dtype=3
                 ),
             )
 
diff --git a/qcengine/programs/mace.py b/qcengine/programs/mace.py
index 5454e05d..f4fba03c 100644
--- a/qcengine/programs/mace.py
+++ b/qcengine/programs/mace.py
@@ -136,7 +136,7 @@ def compute(self, input_data: "AtomicInput", config: "TaskConfig") -> Union["Ato
         ret_data["input_data"] = input_data
         ret_data["molecule"] = input_data.molecule
         ret_data["provenance"] = Provenance(creator="mace", version=mace.__version__, routine="mace")
-        ret_data["schema_name"] = "qcschema_atomic_output"
+        ret_data["schema_name"] = "qcschema_atomic_result"
         ret_data["success"] = True
 
         # Form up a dict first, then sent to BaseModel to avoid repeat kwargs which don't override each other
diff --git a/qcengine/programs/molpro.py b/qcengine/programs/molpro.py
index f0d9b9e2..12d1e9fb 100644
--- a/qcengine/programs/molpro.py
+++ b/qcengine/programs/molpro.py
@@ -443,7 +443,7 @@ def parse_output(self, outfiles: Dict[str, str], input_model: "AtomicInput") ->
         # Final output_data assignments needed for the AtomicResult object
         output_data["properties"] = properties
         output_data["extras"] = extras
-        output_data["schema_name"] = "qcschema_atomic_output"
+        output_data["schema_name"] = "qcschema_atomic_result"
         output_data["stdout"] = outfiles["dispatch.out"]
         output_data["success"] = True
         output_data["provenance"] = input_model.provenance  # TODO better stamp?
diff --git a/qcengine/programs/mrchem.py b/qcengine/programs/mrchem.py
index 8ee66524..4be0bf9a 100644
--- a/qcengine/programs/mrchem.py
+++ b/qcengine/programs/mrchem.py
@@ -101,7 +101,7 @@ def compute(self, input_model: "AtomicInput", config: "TaskConfig") -> "AtomicRe
         input_data = copy.deepcopy(job_input["mrchem_json"])
 
         output_data = {
-            "schema_name": "qcschema_atomic_output",
+            "schema_name": "qcschema_atomic_result",
             "schema_version": 2,
             "input_data": input_model,
             "molecule": input_model.molecule,
@@ -170,9 +170,9 @@ def compute(self, input_model: "AtomicInput", config: "TaskConfig") -> "AtomicRe
                 if input_model.specification.driver == "energy":
                     output_data["return_result"] = mrchem_output["properties"]["scf_energy"]["E_tot"]
                 elif input_model.specification.driver == "gradient":
-                    output_data["return_result"] = mrchem_output["properties"]["geometric_derivative"]["geom-1"][
-                        "total"
-                    ]
+                    grad = mrchem_output["properties"]["geometric_derivative"]["geom-1"]["total"]
+                    output_data["return_result"] = grad
+                    output_data["properties"]["return_gradient"] = grad
                 elif input_model.specification.driver == "properties":
                     output_data["return_result"] = {
                         f"{ks[1]}": {f"{ks[2]}": _nested_get(mrchem_output, ks)} for ks in computed_rsp_props
diff --git a/qcengine/programs/nwchem/harvester.py b/qcengine/programs/nwchem/harvester.py
index b83d81e3..9500be38 100644
--- a/qcengine/programs/nwchem/harvester.py
+++ b/qcengine/programs/nwchem/harvester.py
@@ -2,7 +2,7 @@
 import logging
 import re
 from decimal import Decimal
-from typing import Optional, Tuple
+from typing import TYPE_CHECKING, Optional, Tuple
 
 import numpy as np
 import qcelemental as qcel
@@ -10,6 +10,9 @@
 
 from ..util import PreservingDict
 
+if TYPE_CHECKING:
+    from qcelemental.models.v2 import Molecule
+
 logger = logging.getLogger(__name__)
 
 
@@ -845,7 +848,7 @@ def harvest_outfile_pass(outtext):
                     validate=False,
                     **qcel.molparse.to_schema(
                         qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"],
-                        dtype=2,
+                        dtype=3,
                     ),
                 )
 
@@ -867,7 +870,7 @@ def harvest_outfile_pass(outtext):
                     validate=False,
                     **qcel.molparse.to_schema(
                         qcel.molparse.from_string(molxyz, dtype="xyz+", fix_com=True, fix_orientation=True)["qm"],
-                        dtype=2,
+                        dtype=3,
                     ),
                 )
 
@@ -921,6 +924,7 @@ def harvest_outfile_pass(outtext):
             logger.debug("matched total dipole")
 
             # UNIT = DEBYE(S)
+            d2au = Decimal(qcel.constants.conversion_factor("debye", "e * bohr"))
             psivar[f"CURRENT DIPOLE"] = d2au * np.array([mobj.group(7), mobj.group(8), mobj.group(9)])
             # total?
 
diff --git a/qcengine/programs/psi4.py b/qcengine/programs/psi4.py
index 441c1ea7..cec84e49 100644
--- a/qcengine/programs/psi4.py
+++ b/qcengine/programs/psi4.py
@@ -208,7 +208,7 @@ def compute(self, input_model: "AtomicInput", config: "TaskConfig") -> "AtomicRe
                     error_type = "execution_error"
 
                 # Reset the schema if required
-                output_data["schema_name"] = "qcschema_atomic_output"
+                output_data["schema_name"] = "qcschema_atomic_result"
                 output_data.pop("memory", None)
                 output_data.pop("nthreads", None)
                 output_data["stdout"] = output_data.pop("raw_output", None)
diff --git a/qcengine/programs/qcore.py b/qcengine/programs/qcore.py
index 4303c248..5629a07b 100644
--- a/qcengine/programs/qcore.py
+++ b/qcengine/programs/qcore.py
@@ -242,7 +242,7 @@ def parse_output(self, output: Dict[str, Any], input_model: "AtomicInput") -> "A
 
         output_data["properties"] = properties
 
-        output_data["schema_name"] = "qcschema_atomic_output"
+        output_data["schema_name"] = "qcschema_atomic_result"
         output_data["success"] = True
 
         return AtomicResult(**output_data)
diff --git a/qcengine/programs/qcvar_identities_resources.py b/qcengine/programs/qcvar_identities_resources.py
index ceb3111b..b568b2d2 100644
--- a/qcengine/programs/qcvar_identities_resources.py
+++ b/qcengine/programs/qcvar_identities_resources.py
@@ -3,7 +3,7 @@
 from typing import Any, Dict, List
 
 import numpy as np
-from qcelemental.models.v2 import AtomicResultProperties
+from qcelemental.models.v2 import AtomicProperties
 
 from .util import PreservingDict
 
@@ -427,7 +427,7 @@ def build_out(rawvars: Dict[str, Any], verbose: int = 1) -> None:
 }
 
 
-def build_atomicproperties(qcvars: "PreservingDict") -> AtomicResultProperties:
+def build_atomicproperties(qcvars: "PreservingDict") -> AtomicProperties:
     """For results extracted from QC output in QCDB terminology, translate to QCSchema terminology.
 
     Parameters
@@ -437,8 +437,8 @@ def build_atomicproperties(qcvars: "PreservingDict") -> AtomicResultProperties:
 
     Returns
     -------
-    atprop : AtomicResultProperties
-        Object of calculation information in QCSchema AtomicResultProperties terminology.
+    atprop : AtomicProperties
+        Object of calculation information in QCSchema AtomicProperties terminology.
 
     """
     atprop = {}
@@ -446,4 +446,4 @@ def build_atomicproperties(qcvars: "PreservingDict") -> AtomicResultProperties:
         if pv in qcvars_to_atomicproperties:
             atprop[qcvars_to_atomicproperties[pv]] = dpv
 
-    return AtomicResultProperties(**atprop)
+    return AtomicProperties(**atprop)
diff --git a/qcengine/programs/rdkit.py b/qcengine/programs/rdkit.py
index 423d5a19..3c453948 100644
--- a/qcengine/programs/rdkit.py
+++ b/qcengine/programs/rdkit.py
@@ -126,13 +126,18 @@ def compute(self, input_data: "AtomicInput", config: "TaskConfig") -> "AtomicRes
 
         ff.Initialize()
 
-        ret_data["properties"] = {"return_energy": ff.CalcEnergy() * ureg.conversion_factor("kJ / mol", "hartree")}
+        ret_data["properties"] = {
+            "return_energy": ff.CalcEnergy() * ureg.conversion_factor("kJ / mol", "hartree"),
+            "calcinfo_natom": len(jmol.symbols),
+        }
+        if input_data.specification.driver == "gradient":
+            coef = ureg.conversion_factor("kJ / mol", "hartree") * ureg.conversion_factor("angstrom", "bohr")
+            ret_data["properties"]["return_gradient"] = [x * coef for x in ff.CalcGrad()]
 
         if input_data.specification.driver == "energy":
             ret_data["return_result"] = ret_data["properties"]["return_energy"]
         elif input_data.specification.driver == "gradient":
-            coef = ureg.conversion_factor("kJ / mol", "hartree") * ureg.conversion_factor("angstrom", "bohr")
-            ret_data["return_result"] = [x * coef for x in ff.CalcGrad()]
+            ret_data["return_result"] = ret_data["properties"]["return_gradient"]
         else:
             raise InputError(f"Driver {input_data.specification.driver} not implemented for RDKit.")
 
diff --git a/qcengine/programs/terachem.py b/qcengine/programs/terachem.py
index 3ca0d12e..974233f8 100644
--- a/qcengine/programs/terachem.py
+++ b/qcengine/programs/terachem.py
@@ -182,7 +182,7 @@ def parse_output(self, outfiles: Dict[str, str], input_model: "AtomicInput") ->
 
         output_data["properties"] = properties
 
-        output_data["schema_name"] = "qcschema_atomic_output"
+        output_data["schema_name"] = "qcschema_atomic_result"
         output_data["stdout"] = outfiles["tc.out"]
         # TODO Should only return True if TeraChem calculation terminated properly
         output_data["success"] = True
diff --git a/qcengine/programs/tests/test_dftd3_mp2d.py b/qcengine/programs/tests/test_dftd3_mp2d.py
index 9a73f6da..c2b8a6cc 100644
--- a/qcengine/programs/tests/test_dftd3_mp2d.py
+++ b/qcengine/programs/tests/test_dftd3_mp2d.py
@@ -1394,7 +1394,7 @@ def eneyne_ne_qcdbmols():
         },
         "ne": {"atom": ne},
     }
-    return mols
+    return {2: mols, 3: mols}
 
 
 def eneyne_ne_psi4mols():
@@ -1414,27 +1414,29 @@ def eneyne_ne_psi4mols():
         },
         "ne": {"atom": ne},
     }
-    return mols
+    return {2: mols, 3: mols}
 
 
 def eneyne_ne_qcschemamols():
 
-    eneyne = qcel.molparse.to_schema(qcel.molparse.from_string(seneyne)["qm"], dtype=2)
-    mA = qcel.molparse.to_schema(qcel.molparse.from_string("\n".join(seneyne.splitlines()[:7]))["qm"], dtype=2)
-    mB = qcel.molparse.to_schema(qcel.molparse.from_string("\n".join(seneyne.splitlines()[-4:]))["qm"], dtype=2)
-    ne = qcel.molparse.to_schema(qcel.molparse.from_string(sne)["qm"], dtype=2)
+    mols = {}
+    for dtype in [2, 3]:
+        eneyne = qcel.molparse.to_schema(qcel.molparse.from_string(seneyne)["qm"], dtype=dtype)
+        mA = qcel.molparse.to_schema(qcel.molparse.from_string("\n".join(seneyne.splitlines()[:7]))["qm"], dtype=dtype)
+        mB = qcel.molparse.to_schema(qcel.molparse.from_string("\n".join(seneyne.splitlines()[-4:]))["qm"], dtype=dtype)
+        ne = qcel.molparse.to_schema(qcel.molparse.from_string(sne)["qm"], dtype=dtype)
 
-    mAgB = qcel.molparse.from_string(seneyne)["qm"]
-    mAgB["real"] = [
-        (iat < mAgB["fragment_separators"][0]) for iat in range(len(mAgB["elem"]))
-    ]  # works b/c chgmult doesn't need refiguring
-    mAgB = qcel.molparse.to_schema(mAgB, dtype=2)
+        mAgB = qcel.molparse.from_string(seneyne)["qm"]
+        mAgB["real"] = [
+            (iat < mAgB["fragment_separators"][0]) for iat in range(len(mAgB["elem"]))
+        ]  # works b/c chgmult doesn't need refiguring
+        mAgB = qcel.molparse.to_schema(mAgB, dtype=dtype)
 
-    gAmB = qcel.molparse.from_string(seneyne)["qm"]
-    gAmB["real"] = [(iat >= gAmB["fragment_separators"][0]) for iat in range(len(gAmB["elem"]))]
-    gAmB = qcel.molparse.to_schema(gAmB, dtype=2)
+        gAmB = qcel.molparse.from_string(seneyne)["qm"]
+        gAmB["real"] = [(iat >= gAmB["fragment_separators"][0]) for iat in range(len(gAmB["elem"]))]
+        gAmB = qcel.molparse.to_schema(gAmB, dtype=dtype)
 
-    mols = {"eneyne": {"dimer": eneyne, "mA": mA, "mB": mB, "mAgB": mAgB, "gAmB": gAmB}, "ne": {"atom": ne}}
+        mols[dtype] = {"eneyne": {"dimer": eneyne, "mA": mA, "mB": mB, "mAgB": mAgB, "gAmB": gAmB}, "ne": {"atom": ne}}
     return mols
 
 
@@ -1564,7 +1566,7 @@ def test_3(schema_versions, request):
         resinp = {
             "schema_name": "qcschema_atomic_input",
             "schema_version": 2,
-            "molecule": qcel.molparse.to_schema(sys, dtype=2),
+            "molecule": qcel.molparse.to_schema(sys, dtype=3),
             "specification": {
                 "driver": "energy",
                 "model": {"method": "b3lyp"},
@@ -1599,7 +1601,7 @@ def test_3(schema_versions, request):
             eneyne_ne_qcdbmols, marks=using("psi4")
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
     ],
-    ids=["qmol", "pmol"],
+    ids=["pmol", "qmol"],
 )
 @pytest.mark.parametrize(
     "inp",
@@ -1625,8 +1627,9 @@ def test_3(schema_versions, request):
         ),
     ],
 )
-def test_molecule__run_dftd3__23body(inp, subjects):
-    subject = subjects()[inp["parent"]][inp["subject"]]
+def test_molecule__run_dftd3__23body(inp, subjects, request):
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
@@ -1679,7 +1682,7 @@ def test_qcdb__energy_d3():
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "inp",
@@ -1695,14 +1698,15 @@ def test_qcdb__energy_d3():
 def test_mp2d__run_mp2d__2body(inp, subjects, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     if from_v2(request.node.name):
         resinp = {
@@ -1754,7 +1758,7 @@ def test_mp2d__run_mp2d__2body(inp, subjects, schema_versions, request):
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "program, inp",
@@ -1799,14 +1803,15 @@ def test_mp2d__run_mp2d__2body(inp, subjects, schema_versions, request):
 def test_dftd3__run_dftd3__2body(inp, program, subjects, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     if from_v2(request.node.name):
         atin = models.AtomicInput(
@@ -1856,7 +1861,7 @@ def test_dftd3__run_dftd3__2body(inp, program, subjects, schema_versions, reques
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "inp",
@@ -1874,14 +1879,15 @@ def test_dftd3__run_dftd3__2body(inp, program, subjects, schema_versions, reques
 def test_dftd3__run_dftd3__2body_error(inp, subjects, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     program = "dftd4" if ("D4(BJ" in inp["lbl"]) else "dftd3"
 
@@ -1915,7 +1921,7 @@ def test_dftd3__run_dftd3__2body_error(inp, subjects, schema_versions, request):
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "inp",
@@ -1931,14 +1937,15 @@ def test_dftd3__run_dftd3__2body_error(inp, subjects, schema_versions, request):
 def test_dftd3__run_dftd3__3body(inp, subjects, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     if from_v2(request.node.name):
         resinp = {
@@ -1987,7 +1994,7 @@ def test_dftd3__run_dftd3__3body(inp, subjects, schema_versions, request):
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "inp, extrakw, program",
@@ -2007,14 +2014,15 @@ def test_dftd3__run_dftd3__3body(inp, subjects, schema_versions, request):
 def test_sapt_pairwise(inp, program, extrakw, subjects, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     expected_pairwise = pref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     if from_v2(request.node.name):
         atin = models.AtomicInput(
@@ -2065,7 +2073,7 @@ def test_sapt_pairwise(inp, program, extrakw, subjects, schema_versions, request
         ),  # needs qcdb.Molecule, presently more common in psi4 than in qcdb
         pytest.param(eneyne_ne_qcschemamols),
     ],
-    ids=["qmol", "pmol", "qcmol"],
+    ids=["pmol", "qmol", "qcmol"],
 )
 @pytest.mark.parametrize(
     "inp",
@@ -2081,14 +2089,15 @@ def test_sapt_pairwise(inp, program, extrakw, subjects, schema_versions, request
 def test_gcp(inp, subjects, program, schema_versions, request):
     models, retver, _ = schema_versions
 
-    subject = subjects()[inp["parent"]][inp["subject"]]
+    vmol = 3 if from_v2(request.node.name) else 2
+    subject = subjects()[vmol][inp["parent"]][inp["subject"]]
     expected = ref[inp["parent"]][inp["lbl"]][inp["subject"]]
     gexpected = gref[inp["parent"]][inp["lbl"]][inp["subject"]]
 
     if "qcmol" in request.node.name:
         mol = subject
     else:
-        mol = subject.to_schema(dtype=2)
+        mol = subject.to_schema(dtype=vmol)
 
     if from_v2(request.node.name):
         resinp = {
diff --git a/qcengine/programs/tests/test_ghost.py b/qcengine/programs/tests/test_ghost.py
index 668945d0..888f6c84 100644
--- a/qcengine/programs/tests/test_ghost.py
+++ b/qcengine/programs/tests/test_ghost.py
@@ -182,7 +182,8 @@ def test_simple_ghost(driver, program, basis, keywords, hene_data, schema_versio
 def test_tricky_ghost(driver, qcprog, subject, basis, keywords, schema_versions, request):
     models, retver, _ = schema_versions
 
-    dmol = eneyne_ne_qcschemamols()["eneyne"][subject]
+    vmol = 3 if from_v2(request.node.name) else 2
+    dmol = eneyne_ne_qcschemamols()[vmol]["eneyne"][subject]
     # Freeze the input orientation so that output arrays are aligned to input
     #   and all programs match gradient.
     dmol["fix_com"] = True
diff --git a/qcengine/programs/tests/test_programs.py b/qcengine/programs/tests/test_programs.py
index 7e7e1a5c..59add297 100644
--- a/qcengine/programs/tests/test_programs.py
+++ b/qcengine/programs/tests/test_programs.py
@@ -817,6 +817,7 @@ def test_psi4_properties_driver(schema_versions, request):
     }
     if from_v2(request.node.name):
         json_data["schema_name"] = "qcschema_atomic_input"
+        json_data["schema_version"] = 2
         json_data["specification"] = {
             "driver": json_data.pop("driver"),
             "model": json_data.pop("model"),
@@ -877,7 +878,7 @@ def test_psi4_properties_driver(schema_versions, request):
     json_ret = checkver_and_convert(json_ret, request.node.name, "post")
 
     assert json_ret.success
-    xptd_schema_name = "qcschema_atomic_output" if "v2" in request.node.name else "qcschema_output"
+    xptd_schema_name = "qcschema_atomic_result" if "v2" in request.node.name else "qcschema_output"
     assert json_ret.schema_name == xptd_schema_name
     for k in expected_return_result.keys():
         assert compare_values(expected_return_result[k], json_ret.return_result[k], atol=1.0e-5)
diff --git a/qcengine/programs/torchani.py b/qcengine/programs/torchani.py
index 1e32f731..89820d11 100644
--- a/qcengine/programs/torchani.py
+++ b/qcengine/programs/torchani.py
@@ -181,7 +181,7 @@ def compute(self, input_data: "AtomicInput", config: "TaskConfig") -> "AtomicRes
             creator="torchani", version="unknown", routine="torchani.builtin.aev_computer"
         )
 
-        ret_data["schema_name"] = "qcschema_atomic_output"
+        ret_data["schema_name"] = "qcschema_atomic_result"
         ret_data["success"] = True
 
         # Form up a dict first, then sent to BaseModel to avoid repeat kwargs which don't override each other
diff --git a/qcengine/testing.py b/qcengine/testing.py
index 53d758fe..fafa26cb 100644
--- a/qcengine/testing.py
+++ b/qcengine/testing.py
@@ -7,7 +7,7 @@
 import numpy as np
 import pytest
 import qcelemental as qcel
-from pkg_resources import parse_version
+from packaging.version import parse as parse_version
 from pydantic import ConfigDict
 from qcelemental.util import which, which_import
 
diff --git a/qcengine/tests/test_mdi.py b/qcengine/tests/test_mdi.py
index 285ecaff..a32ae409 100644
--- a/qcengine/tests/test_mdi.py
+++ b/qcengine/tests/test_mdi.py
@@ -3,6 +3,7 @@
 """
 
 
+import qcelemental as qcel
 from qcelemental.testing import compare_values
 
 import qcengine as qcng
@@ -22,7 +23,7 @@ def test_mdi_water():
     engine = qcng.MDIServer(
         "-role DRIVER -name QCEngine -method TEST",
         "psi4",
-        qcng.get_molecule("water"),
+        qcel.models.v2.Molecule(**qcng.get_molecule("water", return_dict=True)),
         {"method": "SCF", "basis": "sto-3g"},
         {"scf_type": "df"},
     )
diff --git a/qcengine/tests/test_procedures.py b/qcengine/tests/test_procedures.py
index 1f092524..a7714b12 100644
--- a/qcengine/tests/test_procedures.py
+++ b/qcengine/tests/test_procedures.py
@@ -4,6 +4,7 @@
 
 import warnings
 
+import numpy as np
 import pytest
 import qcelemental as qcel
 
@@ -57,6 +58,7 @@ def test_geometric_psi4(input_data, optimizer, ncores, schema_versions, request)
         input_data["specification"]["specification"]["program"] = grad_program
         input_data["specification"]["specification"]["extras"] = {"myqctag": "hello psi4"}
         input_data["specification"]["extras"] = {"myopttag": "hello qcengine"}
+        input_data["specification"]["protocols"] = {"trajectory_results": "all"}
 
     else:
         input_data["input_specification"]["model"] = {"method": "HF", "basis": "sto-3g"}
@@ -142,6 +144,7 @@ def test_geometric_local_options(input_data, schema_versions, request, optimizer
     if from_v2(request.node.name):
         input_data["specification"]["specification"]["model"] = {"method": "HF", "basis": "sto-3g"}
         input_data["specification"]["specification"]["program"] = "psi4"
+        input_data["specification"]["protocols"] = {"trajectory_results": "final"}
     else:
         input_data["input_specification"]["model"] = {"method": "HF", "basis": "sto-3g"}
         input_data["keywords"]["program"] = "psi4"
@@ -194,6 +197,7 @@ def test_optimizer_stdout(optimizer, gradprog, gradmodel, converged, input_data,
         input_data["specification"]["specification"]["model"] = gradmodel
         input_data["specification"]["specification"]["program"] = gradprog
         input_data["specification"]["specification"]["protocols"] = {"stdout": False}
+        input_data["specification"]["protocols"] = {"trajectory_results": "all"}
     else:
         input_data["input_specification"]["model"] = gradmodel
         input_data["keywords"]["program"] = gradprog
@@ -217,6 +221,79 @@ def test_optimizer_stdout(optimizer, gradprog, gradmodel, converged, input_data,
     else:
         assert trajs_tgt[0].stdout is not None
 
+    if "v2" in request.node.name:
+        assert ret.properties.optimization_iterations in [3, 4, 5]
+        assert len(ret.trajectory_properties) == ret.properties.optimization_iterations
+        assert np.allclose(ret.properties.return_gradient, ret.trajectory_results[-1].return_result)
+        assert ret.trajectory_results[-1].properties.return_energy == pytest.approx(
+            ret.properties.return_energy, 1.0e-4
+        )
+        if gradprog == "psi4":
+            assert ret.properties.nuclear_repulsion_energy == pytest.approx(8.906, 1.0e-3)
+            assert ret.properties.return_energy == pytest.approx(-74.96599, 1.0e-4)
+        elif gradprog == "rdkit":
+            assert ret.properties.nuclear_repulsion_energy == pytest.approx(8.888, 1.0e-3)
+            assert ret.properties.return_energy == pytest.approx(0, 1.0e-2)  # 0 seems wrong
+
+
+@pytest.mark.parametrize(
+    "optimizer",
+    [
+        pytest.param("geometric", marks=using("geometric")),
+        pytest.param("optking", marks=using("optking")),
+        pytest.param("berny", marks=using("berny")),
+    ],
+)
+@pytest.mark.parametrize(
+    "gradprog,gradmodel",
+    [
+        pytest.param("rdkit", {"method": "UFF", "basis": ""}, marks=using("rdkit")),
+        # pytest.param("psi4", {"method": "HF", "basis": "sto-3g"}, marks=using("psi4")),
+    ],
+)
+@pytest.mark.parametrize("traj_ptcl", ["none", "final", "all", "default"])
+def test_optimizer_protocols(optimizer, gradprog, gradmodel, input_data, schema_versions, request, traj_ptcl):
+    models, retver, _ = schema_versions
+
+    input_data["initial_molecule"] = models.Molecule(**qcng.get_molecule("water", return_dict=True))
+
+    if from_v2(request.node.name):
+        input_data["specification"]["specification"]["model"] = gradmodel
+        input_data["specification"]["specification"]["program"] = gradprog
+        if traj_ptcl != "default":
+            input_data["specification"]["protocols"] = {"trajectory_results": traj_ptcl}
+    else:
+        input_data["input_specification"]["model"] = gradmodel
+        input_data["keywords"]["program"] = gradprog
+        if traj_ptcl != "default":
+            input_data["protocols"] = {"trajectory": traj_ptcl}
+
+    input_data = models.OptimizationInput(**input_data)
+
+    input_data = checkver_and_convert(input_data, request.node.name, "pre")
+    ret = qcng.compute(input_data, optimizer, raise_error=True, return_version=retver)
+    ret = checkver_and_convert(ret, request.node.name, "post")
+
+    assert ret.success is True
+
+    trajs_tgt = ret.trajectory_results if "v2" in request.node.name else ret.trajectory
+    assert ret.provenance.creator.lower() == optimizer
+
+    if "v2" in request.node.name:
+        assert len(ret.trajectory_properties) > 2
+
+    if traj_ptcl == "none":
+        assert len(trajs_tgt) == 0
+    elif traj_ptcl == "final":
+        assert len(trajs_tgt) == 1
+    elif traj_ptcl == "all":
+        assert len(trajs_tgt) > 2
+    elif traj_ptcl == "default":
+        if from_v2(request.node.name):
+            assert len(trajs_tgt) == 0
+        else:
+            assert len(trajs_tgt) > 2
+
 
 @using("psi4")
 @using("berny")
@@ -293,7 +370,7 @@ def test_optimization_protocols(optimizer, input_data, schema_versions, request)
     if from_v2(request.node.name):
         input_data["specification"]["specification"]["model"] = grad_model
         input_data["specification"]["specification"]["program"] = grad_program
-        input_data["specification"]["protocols"] = {"trajectory": "initial_and_final"}
+        input_data["specification"]["protocols"] = {"trajectory_results": "initial_and_final"}
     else:
         input_data["input_specification"]["model"] = grad_model
         input_data["keywords"]["program"] = grad_program
@@ -344,6 +421,7 @@ def test_geometric_retries(failure_engine, input_data, schema_versions, request)
         input_data["specification"]["keywords"][
             "coordsys"
         ] = "cart"  # needed by geometric v1.0 to play nicely with failure_engine
+        input_data["specification"]["protocols"] = {"trajectory_results": "all"}
     else:
         input_data["input_specification"]["model"] = {"method": "something"}
         input_data["keywords"]["program"] = failure_engine.name
@@ -442,6 +520,7 @@ def test_geometric_generic(input_data, program, model, bench, schema_versions, r
         input_data["specification"]["specification"]["extras"] = {
             "_secret_tags": {"mysecret_tag": "data1"}  # pragma: allowlist secret
         }
+        input_data["specification"]["protocols"] = {"trajectory_results": "final"}
     else:
         input_data["input_specification"]["model"] = model
         input_data["keywords"]["program"] = program
@@ -484,6 +563,7 @@ def test_nwchem_relax(linopt, schema_versions, request):
                     "keywords": {"set__driver:linopt": linopt},
                     "driver": "gradient",
                 },
+                "protocols": {"trajectory_results": "all"},
             },
             "initial_molecule": models.Molecule(**qcng.get_molecule("hydrogen", return_dict=True)),
         }
@@ -553,15 +633,18 @@ def test_nwchem_restart(tmpdir, schema_versions, request):
 
 @using("rdkit")
 @using("torsiondrive")
-def test_torsiondrive_generic(schema_versions, request):
+@pytest.mark.parametrize("scan_ptcl", ["none", "all", "lowest"])
+def test_torsiondrive_generic(schema_versions, request, scan_ptcl):
     models, retver, _ = schema_versions
 
     if from_v2(request.node.name):
         input_data = models.TorsionDriveInput(
-            initial_molecules=[models.Molecule(**qcng.get_molecule("ethane", return_dict=True))] * 2,
+            initial_molecule=[models.Molecule(**qcng.get_molecule("ethane", return_dict=True))] * 2,
             specification=models.TorsionDriveSpecification(
-                keywords=models.TDKeywords(dihedrals=[(2, 0, 1, 5)], grid_spacing=[180]),
+                keywords=models.TorsionDriveKeywords(dihedrals=[(2, 0, 1, 5)], grid_spacing=[180]),
+                protocols=models.TorsionDriveProtocols(scan_results=scan_ptcl),
                 specification=models.OptimizationSpecification(
+                    protocols=models.OptimizationProtocols(trajectory_results="all"),
                     program="geomeTRIC",
                     keywords={
                         "coordsys": "hdlc",
@@ -599,8 +682,10 @@ def test_torsiondrive_generic(schema_versions, request):
     assert ret.success
 
     expected_grid_ids = {"180", "0"}
+    opthist_tgt = ret.scan_results if "v2" in request.node.name else ret.optimization_history
 
-    assert {*ret.optimization_history} == expected_grid_ids
+    if not (from_v2(request.node.name) and scan_ptcl == "none"):
+        assert {*opthist_tgt} == expected_grid_ids
 
     assert {*ret.final_energies} == expected_grid_ids
     assert {*ret.final_molecules} == expected_grid_ids
@@ -612,13 +697,26 @@ def test_torsiondrive_generic(schema_versions, request):
     assert pytest.approx(ret.final_molecules["0"].measure([2, 0, 1, 5]), abs=1.0e-2) == 0.0
 
     assert ret.provenance.creator.lower() == "torsiondrive"
-    assert ret.optimization_history["180"][0].provenance.creator.lower() == "geometric"
+
+    if from_v2(request.node.name):
+        # properly "to_v1" should be always `== 4` but data lost in v2 by default
+        if scan_ptcl == "none":
+            assert not opthist_tgt
+        else:
+            assert len(opthist_tgt["180"]) == {"lowest": 1, "all": 4}[scan_ptcl]
+    else:
+        assert len(opthist_tgt["180"]) == 4
+
+    if not (from_v2(request.node.name) and scan_ptcl == "none"):
+        assert opthist_tgt["180"][0].provenance.creator.lower() == "geometric"
     if "v2" in request.node.name:
-        assert ret.optimization_history["180"][0].trajectory_results[0].provenance.creator.lower() == "rdkit"
-        assert ret.optimization_history["180"][0].trajectory_results[0].schema_version == 2
+        if scan_ptcl != "none":
+            assert opthist_tgt["180"][0].trajectory_results[0].provenance.creator.lower() == "rdkit"
+            assert opthist_tgt["180"][0].trajectory_results[0].schema_version == 2
     else:
-        assert ret.optimization_history["180"][0].trajectory[0].provenance.creator.lower() == "rdkit"
-        assert ret.optimization_history["180"][0].trajectory[0].schema_version == 1
+        if not ("to_v1" in request.node.name and scan_ptcl == "none"):
+            assert opthist_tgt["180"][0].trajectory[0].provenance.creator.lower() == "rdkit"
+            assert opthist_tgt["180"][0].trajectory[0].schema_version == 1
 
     assert ret.stdout == "All optimizations converged at lowest energy. Job Finished!\n"
 
@@ -725,6 +823,7 @@ def test_optimization_mrchem(input_data, optimizer, schema_versions, request):
         input_data["specification"]["specification"]["model"] = {"method": "HF"}
         input_data["specification"]["specification"]["keywords"] = {"world_prec": 1.0e-4}
         input_data["specification"]["specification"]["program"] = "mrchem"
+        input_data["specification"]["protocols"] = {"trajectory_results": "final"}  # to test provenance
     else:
         input_data["input_specification"]["model"] = {"method": "HF"}
         input_data["input_specification"]["keywords"] = {"world_prec": 1.0e-4}
@@ -736,9 +835,13 @@ def test_optimization_mrchem(input_data, optimizer, schema_versions, request):
     ret = qcng.compute(input_data, optimizer, raise_error=True, return_version=retver)
     ret = checkver_and_convert(ret, request.node.name, "post")
 
-    trajs_tgt = ret.trajectory_results if "v2" in request.node.name else ret.trajectory
+    trajs_tgt = ret.trajectory_results if "v2" in request.node.name else ret.trajectory  # for looking at grad jobs
+    props_tgt = ret.trajectory_properties if "v2" in request.node.name else ret.energies  # for counting opt iterations
+
+    assert 10 > len(props_tgt) > 1
+    if "v2" in request.node.name:
+        assert 10 > ret.properties.optimization_iterations > 1
 
-    assert 10 > len(trajs_tgt) > 1
     assert pytest.approx(ret.final_molecule.measure([0, 1]), 1.0e-3) == 1.3860734486984705
     assert ret.provenance.creator.lower() == optimizer
     assert trajs_tgt[0].provenance.creator.lower() == "mrchem"
diff --git a/setup.py b/setup.py
index 381f2ec0..911eb23d 100644
--- a/setup.py
+++ b/setup.py
@@ -47,7 +47,7 @@
                 "numpydoc",
                 "autodoc-pydantic>=2.0",
             ],
-            "tests": ["pytest", "pytest-cov"],
+            "tests": ["pytest", "pytest-cov", "packaging"],
             "lint": ["black>=22.1.0,<23.0a0"],
         },
         tests_require=["pytest", "pytest-cov"],