NEB update

.github/workflows/core_tests.yml

@@ -47,7 +47,7 @@ jobs:
         shell: bash -l {0}
         run: |
           pip install -e ./qcportal -e ./qcfractalcompute -e ./qcfractal[services,geoip,snowflake] -e ./qcarchivetesting
-          pip install scipy "geometric @ git+https://github.com/hjnpark/geomeTRIC"
+          pip install scipy "geometric @ git+https://github.com/leeping/geomeTRIC"
 
       - name: Run tests
         shell: bash -l {0}

qcarchivetesting/conda-envs/fulltest_server.yaml

@@ -39,5 +39,5 @@ dependencies:
   - torsiondrive
 
   - pip:
-      - "geometric @ git+https://github.com/hjnpark/geomeTRIC"
+      - "geometric @ git+https://github.com/leeping/geomeTRIC"
       - scipy # for geometric

qcarchivetesting/conda-envs/fulltest_snowflake.yaml

@@ -51,5 +51,5 @@ dependencies:
   - pytest
 
   - pip:
-      - "geometric @ git+https://github.com/hjnpark/geomeTRIC"
+      - "geometric @ git+https://github.com/leeping/geomeTRIC"
       - scipy # for geometric

qcarchivetesting/conda-envs/fulltest_worker.yaml

@@ -36,4 +36,4 @@ dependencies:
   # Geometric service
   - pip:
     - scipy
-    - "geometric @ git+https://github.com/hjnpark/geomeTRIC"
+    - "geometric @ git+https://github.com/leeping/geomeTRIC"

qcarchivetesting/qcarchivetesting/test_full_neb.py

@@ -14,24 +14,107 @@
 
 
 def test_neb_full_1(fulltest_client: PortalClient):
+    """
+    Full NEB test with aligning images.
+    """
     chain = [load_molecule_data("neb/neb_HCN_%i" % i) for i in range(11)]
+
+    # "maximum_force" and "average_force" are set to be high in orger to converge quickly.
+    neb_keywords = NEBKeywords(
+        images=7,
+        spring_constant=1,
+        spring_type=0,
+        maximum_force=3,
+        average_force=2,
+        maximum_cycle=100,
+        optimize_ts=False,
+        optimize_endpoints=False,
+        align=True,
+        epsilon=1e-6,
+        hessian_reset=True,
+    )
+
+    sp_spec = QCSpecification(
+        program="psi4",
+        driver="gradient",
+        method="hf",
+        basis="3-21g",
+        keywords={},
+    )
+
+    meta, ids = fulltest_client.add_nebs(
+        initial_chains=[chain],
+        program="geometric",
+        singlepoint_specification=sp_spec,
+        optimization_specification=None,
+        keywords=neb_keywords,
+    )
+
+    for i in range(600):
+        time.sleep(15)
+        rec = fulltest_client.get_nebs(ids[0])
+        if rec.status not in [RecordStatusEnum.running, RecordStatusEnum.waiting]:
+            break
+    else:
+        raise RuntimeError("Did not finish calculation in time")
+
+    ts_guess = rec.result
+    initial_chain = rec.initial_chain  # List[Molecule]
+    final_chain = rec.final_chain  # List[Singlepoints]
+    optimizations = rec.optimizations
+    singlepoints = rec.singlepoints
+    ts_optimization = rec.ts_optimization
+    hessian = rec.ts_hessian
+
+    # Finding the highest energy image from the last iteration SinglepointRecords.
+    last_sps = singlepoints[max(singlepoints.keys())]
+    sp_energies = [(sp.properties["current energy"], sp.molecule_id) for sp in last_sps]
+    energy, neb_result_id = max(sp_energies)
+
+    # Completed?
+    assert rec.status == RecordStatusEnum.complete
+    # Initial chain length and number of singlepoint records from the last iteration should be 11.
+    assert len(initial_chain) == 7 and len(initial_chain) == len(final_chain)
+    # SinglepointRecords ids of final chain should be the same as the last iteration SinglepointRecords from rec.singlepoints.
+    assert sum([1 if i.id == j.id else 0 for i, j in zip(final_chain, singlepoints[max(singlepoints.keys())])]) == 7
+    # There shouldn't be any OptimizationRecords.
+    assert len(optimizations) == 0
+    # There shouldn't be ts_optimization record.
+    assert ts_optimization is None
+    # When optimize_ts is False, there should not ba a record for the Hessian.
+    assert hessian is None
+    # SinglepointRecords should not have 'return_hessian' either.
+    assert sum(1 if singlepoints[0][i].properties["return_hessian"] is None else 0 for i in range(7)) == 7
+    # Result of the neb and the highest energy image of the last iteration should have the same molecule id.
+    assert ts_guess.id == neb_result_id
+
+
+def test_neb_full_2(fulltest_client: PortalClient):
+    """
+    Full NEB test without aligning images, with endpoint optimizations and TS optimization.
+    This test takes longer than the previous one.
+    """
+    chain = [load_molecule_data("neb/neb_HCN_%i" % i) for i in range(11)]
+
     neb_keywords = NEBKeywords(
         images=11,
         spring_constant=1,
-        optimize_endpoints=True,
-        maximum_force=0.02,
-        average_force=0.02,
+        spring_type=0,
+        maximum_force=0.05,
+        average_force=0.025,
+        maximum_cycle=100,
         optimize_ts=True,
+        optimize_endpoints=True,
+        align=False,
         epsilon=1e-6,
         hessian_reset=True,
-        spring_type=0,
     )
 
     sp_spec = QCSpecification(
         program="psi4",
         driver="gradient",
         method="hf",
-        basis="6-31g",
+        basis="3-21g",
         keywords={},
     )
 
@@ -56,4 +139,34 @@ def test_neb_full_1(fulltest_client: PortalClient):
     else:
         raise RuntimeError("Did not finish calculation in time")
 
+    ts_guess = rec.result
+    initial_chain = rec.initial_chain  # List[Molecule]
+    final_chain = rec.final_chain  # List[Singlepoints]
+    optimizations = rec.optimizations
+    singlepoints = rec.singlepoints
+    ts_optimization = rec.ts_optimization
+    hessian = rec.ts_hessian
+
+    # Finding the highest energy image from the last iteration SinglepointRecords.
+    last_sps = singlepoints[max(singlepoints.keys())]
+    sp_energies = [(sp.properties["current energy"], sp.molecule_id) for sp in last_sps]
+    energy, neb_result_id = max(sp_energies)
+
+    # Completed?
     assert rec.status == RecordStatusEnum.complete
+    # Initial chain length and number of singlepoint records from the last iteration should be 11.
+    assert len(initial_chain) == 11 and len(initial_chain) == len(final_chain)
+    # SinglepointRecords ids of final chain should be the same as the last iteration SinglepointRecords from rec.singlepoints.
+    assert sum([1 if i.id == j.id else 0 for i, j in zip(final_chain, singlepoints[max(singlepoints.keys())])]) == 11
+    # Total 3 OptimizationRecords: two for the endpoints, one for the TS optimization.
+    assert len(optimizations) == 3
+    # rec.tsoptimization should have the same id as the transition record in rec.optimizations.
+    assert optimizations.get("transition").id == ts_optimization.id
+    # When optimize_ts is True, SinglepointRecord containing the Hessian should exist.
+    assert hessian
+    # The rec.hessian should have the Hessian used for the TS optimization.
+    assert hessian.properties["return_hessian"] is not None
+    # Other SinglepointRecords should not have 'return_hessian'
+    assert sum(1 if singlepoints[0][i].properties["return_hessian"] is None else 0 for i in range(11)) == 11
+    # Result of the neb and the highest energy image of the last iteration should have the same molecule id.
+    assert ts_guess.id == neb_result_id

qcfractal/qcfractal/components/neb/record_socket.py

@@ -202,7 +202,7 @@ def iterate_service(
                     initial_molecules[-1] = complete_opts[-1].record.final_molecule.to_model(Molecule)
 
                 with capture_all_output("geometric.nifty") as (rdout, _):
-                    respaced_chain = geometric.qcf_neb.arrange(initial_molecules)
+                    respaced_chain = geometric.qcf_neb.arrange(initial_molecules, params.get("align"))
 
                 output += "\n" + rdout.getvalue()
 

qcfractal/qcfractal/components/neb/test_record_socket.py

@@ -220,7 +220,7 @@ def test_neb_socket_run(
 
     time_0 = now_at_utc()
     finished, n_spopt = run_service(
-        storage_socket, activated_manager_name, id_1[0], generate_task_key, result_data_1, 100
+        storage_socket, activated_manager_name, id_1[0], generate_task_key, result_data_1, 150
     )
     time_1 = now_at_utc()
 

qcportal/qcportal/neb/record_models.py

@@ -61,6 +61,8 @@ class Config:
         description="Setting it equal to True will optimize two end points of the initial chain before starting NEB.",
     )
 
+    align: bool = Field(True, description="Align the images before starting the NEB calculation.")
+
     epsilon: float = Field(1e-5, description="Small eigenvalue threshold for resetting Hessian.")
 
     hessian_reset: bool = Field(
@@ -133,6 +135,7 @@ class NEBRecord(BaseRecord):
     initial_chain_molecule_ids_: Optional[List[int]] = None
     singlepoints_: Optional[List[NEBSinglepoint]] = None
     optimizations_: Optional[Dict[str, NEBOptimization]] = None
+    ts_hessian_: Optional[SinglepointRecord] = None
     neb_result_: Optional[Molecule] = None
     initial_chain_: Optional[List[Molecule]] = None
 
@@ -214,6 +217,12 @@ def _fetch_singlepoints(self):
             self._singlepoints_cache.setdefault(sp_info.chain_iteration, list())
             self._singlepoints_cache[sp_info.chain_iteration].append(sp_rec)
 
+        if len(self._singlepoints_cache) > 0:
+            if len(self._singlepoints_cache[max(self._singlepoints_cache)]) == 1:
+                _, temp_list = self._singlepoints_cache.popitem()
+                self.ts_hessian_ = temp_list[0]
+                assert self.ts_hessian_.specification.driver == "hessian"
+
         self.propagate_client(self._client)
 
     def _fetch_initial_chain(self):
@@ -257,14 +266,18 @@ def initial_chain(self) -> List[Molecule]:
             self._fetch_initial_chain()
         return self.initial_chain_
 
+    @property
+    def final_chain(self) -> List[SinglepointRecord]:
+        return self.singlepoints[max(self.singlepoints.keys())]
+
     @property
     def singlepoints(self) -> Dict[int, List[SinglepointRecord]]:
         if self._singlepoints_cache is None:
             self._fetch_singlepoints()
         return self._singlepoints_cache
 
     @property
-    def neb_result(self):
+    def result(self):
         if self.neb_result_ is None and "neb_result_" not in self.__fields_set__:
             self._fetch_neb_result()
         return self.neb_result_
@@ -278,3 +291,9 @@ def optimizations(self) -> Optional[Dict[str, OptimizationRecord]]:
     @property
     def ts_optimization(self) -> Optional[OptimizationRecord]:
         return self.optimizations.get("transition", None)
+
+    @property
+    def ts_hessian(self) -> Optional[SinglepointRecord]:
+        if self._singlepoints_cache is None:
+            self._fetch_singlepoints()
+        return self.ts_hessian_

qcportal/qcportal/neb/test_record_models.py

@@ -40,13 +40,11 @@ def test_neb_record_model(snowflake: QCATestingSnowflake, includes: Optional[Lis
 
     assert len(record.singlepoints) > 0
 
-    # last one is a hessian calculation?
-    all_sps = list(record.singlepoints.values())
-    for sps in all_sps[:-1]:
-        assert len(sps) == len(molecules)
+    assert all(len(sp) == len(molecules) for sp in record.singlepoints.values())
 
-    assert len(all_sps[-1]) == 1
-    assert all_sps[-1][0].specification.driver == "hessian"
+    ts_hessian = record.ts_hessian
+    if ts_hessian is not None:
+        assert ts_hessian.specification.driver == "hessian"
 
     assert "initial" in record.optimizations
     assert "final" in record.optimizations