diff --git a/.github/workflows/core_tests.yml b/.github/workflows/core_tests.yml index 1b0a30d07..ceda1fe55 100644 --- a/.github/workflows/core_tests.yml +++ b/.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} diff --git a/qcarchivetesting/conda-envs/fulltest_server.yaml b/qcarchivetesting/conda-envs/fulltest_server.yaml index 0353d600a..3d7526751 100644 --- a/qcarchivetesting/conda-envs/fulltest_server.yaml +++ b/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 diff --git a/qcarchivetesting/conda-envs/fulltest_snowflake.yaml b/qcarchivetesting/conda-envs/fulltest_snowflake.yaml index 0304c32ac..8458a7716 100644 --- a/qcarchivetesting/conda-envs/fulltest_snowflake.yaml +++ b/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 diff --git a/qcarchivetesting/conda-envs/fulltest_worker.yaml b/qcarchivetesting/conda-envs/fulltest_worker.yaml index 609fdfb9f..d1690defe 100644 --- a/qcarchivetesting/conda-envs/fulltest_worker.yaml +++ b/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" diff --git a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_b3lyp_opt3.json.xz b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_b3lyp_opt3.json.xz index 5f40d7fdf..4969453bb 100644 Binary files a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_b3lyp_opt3.json.xz and b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_b3lyp_opt3.json.xz differ diff --git a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe.json.xz b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe.json.xz index 0948e32c1..2d515cf47 100644 Binary files a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe.json.xz and b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe.json.xz differ diff --git a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe0_opt1.json.xz b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe0_opt1.json.xz index 3b2334689..b9b9d4f61 100644 Binary files a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe0_opt1.json.xz and b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe0_opt1.json.xz differ diff --git a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt2.json.xz b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt2.json.xz index 7d679b991..1c016a3d5 100644 Binary files a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt2.json.xz and b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt2.json.xz differ diff --git a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt_diff.json.xz b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt_diff.json.xz index a0100821d..00fa07f62 100644 Binary files a/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt_diff.json.xz and b/qcarchivetesting/qcarchivetesting/procedure_data/neb_HCN_psi4_pbe_opt_diff.json.xz differ diff --git a/qcarchivetesting/qcarchivetesting/test_full_neb.py b/qcarchivetesting/qcarchivetesting/test_full_neb.py index e71a98a00..f06897f34 100644 --- a/qcarchivetesting/qcarchivetesting/test_full_neb.py +++ b/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 diff --git a/qcfractal/qcfractal/components/neb/record_socket.py b/qcfractal/qcfractal/components/neb/record_socket.py index cd0e25266..aa162e618 100644 --- a/qcfractal/qcfractal/components/neb/record_socket.py +++ b/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() diff --git a/qcfractal/qcfractal/components/neb/test_record_socket.py b/qcfractal/qcfractal/components/neb/test_record_socket.py index cef335abe..f0a07b6ae 100644 --- a/qcfractal/qcfractal/components/neb/test_record_socket.py +++ b/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() diff --git a/qcportal/qcportal/neb/record_models.py b/qcportal/qcportal/neb/record_models.py index 193d6a491..7a752605b 100644 --- a/qcportal/qcportal/neb/record_models.py +++ b/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,6 +266,10 @@ 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: @@ -264,7 +277,7 @@ def singlepoints(self) -> Dict[int, List[SinglepointRecord]]: 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_ diff --git a/qcportal/qcportal/neb/test_record_models.py b/qcportal/qcportal/neb/test_record_models.py index e24f78998..cd561ce1c 100644 --- a/qcportal/qcportal/neb/test_record_models.py +++ b/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