From 4619851f0571bc72c904d8f7de495c011a6b571f Mon Sep 17 00:00:00 2001 From: Heejune Park Date: Fri, 13 Oct 2023 12:08:02 -0700 Subject: [PATCH 1/4] neb.final_chain is added --- qcportal/qcportal/neb/record_models.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/qcportal/qcportal/neb/record_models.py b/qcportal/qcportal/neb/record_models.py index b1827325d..e0f21a84e 100644 --- a/qcportal/qcportal/neb/record_models.py +++ b/qcportal/qcportal/neb/record_models.py @@ -136,6 +136,7 @@ class NEBRecord(BaseRecord): # Caches ######################################## initial_chain_: Optional[List[Molecule]] = None + final_chain_: Optional[List[SinglepointRecord]] = None optimizations_cache_: Optional[Dict[str, OptimizationRecord]] = None singlepoints_cache_: Optional[Dict[int, List[SinglepointRecord]]] = None @@ -213,6 +214,8 @@ def _handle_includes(self, includes: Optional[Iterable[str]]): if "initial_chain" in includes: self._fetch_initial_chain() + if "final_chain" in includes: + self._fetch_singlepoints() if "singlepoints" in includes: self._fetch_singlepoints() if "optimizations" in includes: @@ -224,6 +227,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: From fb75107160145d5a565c73545e02596955c2e985 Mon Sep 17 00:00:00 2001 From: Heejune Park Date: Fri, 13 Oct 2023 17:04:42 -0700 Subject: [PATCH 2/4] neb test is updated --- .../qcarchivetesting/test_full_neb.py | 115 +++++++++++++++++- .../qcfractal/components/neb/record_socket.py | 1 + 2 files changed, 114 insertions(+), 2 deletions(-) diff --git a/qcarchivetesting/qcarchivetesting/test_full_neb.py b/qcarchivetesting/qcarchivetesting/test_full_neb.py index e71a98a00..f0563b1b1 100644 --- a/qcarchivetesting/qcarchivetesting/test_full_neb.py +++ b/qcarchivetesting/qcarchivetesting/test_full_neb.py @@ -14,13 +14,16 @@ def test_neb_full_1(fulltest_client: PortalClient): + """ + Full NEB test with optimizing end points and the guessed TS (result of the NEB method) optimization. + """ 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, + maximum_force=0.05, + average_force=0.025, optimize_ts=True, epsilon=1e-6, hessian_reset=True, @@ -56,4 +59,112 @@ def test_neb_full_1(fulltest_client: PortalClient): else: raise RuntimeError("Did not finish calculation in time") + ts_guess = rec.neb_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 + + # Finding the highest energy image from the last iteration SinglepointRecords. + neb_result_id = 0 + energy = -9999999 + for sp in singlepoints[max(singlepoints.keys())]: + if sp.properties["current energy"] > energy: + energy = sp.properties["current energy"] + neb_result_id = sp.molecule_id + + # 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 OptimizationRecord + 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, rec.singlepoints should have -1 key. + assert -1 in singlepoints + # The singlepoints[-1] should have the Hessian used for the TS optimization. + assert singlepoints[-1][0].properties["return_hessian"] is not None + # And other SP records should not have 'return_hessian' + assert singlepoints[0][0].properties["return_hessian"] is None + # 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): + """ + Identical as the previous test without optimizing endpoints and transition state. + """ + chain = [load_molecule_data("neb/neb_HCN_%i" % i) for i in range(11)] + neb_keywords = NEBKeywords( + images=11, + spring_constant=1, + optimize_endpoints=False, + maximum_force=0.05, + average_force=0.025, + optimize_ts=False, + epsilon=1e-6, + hessian_reset=True, + spring_type=0, + ) + + sp_spec = QCSpecification( + program="psi4", + driver="gradient", + method="hf", + basis="6-31g", + keywords={}, + ) + + opt_spec = OptimizationSpecification( + program="geometric", + qc_specification=sp_spec, + ) + + meta, ids = fulltest_client.add_nebs( + initial_chains=[chain], + program="geometric", + singlepoint_specification=sp_spec, + optimization_specification=opt_spec, + 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.neb_result + initial_chain = rec.initial_chain # List[Molecule] + final_chain = rec.final_chain # List[Singlepoints] + optimizations = rec.optimizations + singlepoints = rec.singlepoints + + # Finding the highest energy image from the last iteration SinglepointRecords. + neb_result_id = 0 + energy = -9999999 + for sp in singlepoints[max(singlepoints.keys())]: + if sp.properties["current energy"] > energy: + energy = sp.properties["current energy"] + neb_result_id = sp.molecule_id + + # 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 + # There shouldn't be any OptimizationRecords. + assert len(optimizations) == 0 + # There should not be ts_optimization record. + assert rec.ts_optimization is None + # When optimize_ts is False, rec.singlepoints should not have -1 key. + assert -1 not in singlepoints + # 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 66444b8af..315ffc76d 100644 --- a/qcfractal/qcfractal/components/neb/record_socket.py +++ b/qcfractal/qcfractal/components/neb/record_socket.py @@ -400,6 +400,7 @@ def submit_singlepoints( opt_spec = neb_orm.specification.optimization_specification.model_dict() qc_spec = opt_spec["qc_specification"] qc_spec["driver"] = "hessian" + service_state.iteration = -1 meta, sp_ids = self.root_socket.records.singlepoint.add( chain, From 803cc51f5e71ef2e1016e0d28bd521aad4268e65 Mon Sep 17 00:00:00 2001 From: Heejune Park Date: Mon, 16 Oct 2023 17:32:33 -0700 Subject: [PATCH 3/4] SinglepointRecord for the TS Hessian can be request by rec.ts_hessian --- .../qcarchivetesting/test_full_neb.py | 38 +++++++++---------- .../qcfractal/components/neb/record_socket.py | 1 - qcportal/qcportal/neb/record_models.py | 18 +++++++-- 3 files changed, 32 insertions(+), 25 deletions(-) diff --git a/qcarchivetesting/qcarchivetesting/test_full_neb.py b/qcarchivetesting/qcarchivetesting/test_full_neb.py index f0563b1b1..27675ad2d 100644 --- a/qcarchivetesting/qcarchivetesting/test_full_neb.py +++ b/qcarchivetesting/qcarchivetesting/test_full_neb.py @@ -65,14 +65,12 @@ def test_neb_full_1(fulltest_client: PortalClient): 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. - neb_result_id = 0 - energy = -9999999 - for sp in singlepoints[max(singlepoints.keys())]: - if sp.properties["current energy"] > energy: - energy = sp.properties["current energy"] - neb_result_id = sp.molecule_id + 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 @@ -80,16 +78,16 @@ def test_neb_full_1(fulltest_client: PortalClient): 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 OptimizationRecord + # Total 3 OptimizationRecords 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, rec.singlepoints should have -1 key. - assert -1 in singlepoints - # The singlepoints[-1] should have the Hessian used for the TS optimization. - assert singlepoints[-1][0].properties["return_hessian"] is not None - # And other SP records should not have 'return_hessian' - assert singlepoints[0][0].properties["return_hessian"] is None + # 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 + # And 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 @@ -140,6 +138,7 @@ def test_neb_full_2(fulltest_client: PortalClient): else: raise RuntimeError("Did not finish calculation in time") + hessian = rec.ts_hessian # Calling the Hessian first ts_guess = rec.neb_result initial_chain = rec.initial_chain # List[Molecule] final_chain = rec.final_chain # List[Singlepoints] @@ -147,12 +146,9 @@ def test_neb_full_2(fulltest_client: PortalClient): singlepoints = rec.singlepoints # Finding the highest energy image from the last iteration SinglepointRecords. - neb_result_id = 0 - energy = -9999999 - for sp in singlepoints[max(singlepoints.keys())]: - if sp.properties["current energy"] > energy: - energy = sp.properties["current energy"] - neb_result_id = sp.molecule_id + 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 @@ -164,7 +160,7 @@ def test_neb_full_2(fulltest_client: PortalClient): assert len(optimizations) == 0 # There should not be ts_optimization record. assert rec.ts_optimization is None - # When optimize_ts is False, rec.singlepoints should not have -1 key. - assert -1 not in singlepoints + # When optimize_ts is False, there should not ba a record for the Hessian. + assert hessian is None # 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 315ffc76d..66444b8af 100644 --- a/qcfractal/qcfractal/components/neb/record_socket.py +++ b/qcfractal/qcfractal/components/neb/record_socket.py @@ -400,7 +400,6 @@ def submit_singlepoints( opt_spec = neb_orm.specification.optimization_specification.model_dict() qc_spec = opt_spec["qc_specification"] qc_spec["driver"] = "hessian" - service_state.iteration = -1 meta, sp_ids = self.root_socket.records.singlepoint.add( chain, diff --git a/qcportal/qcportal/neb/record_models.py b/qcportal/qcportal/neb/record_models.py index e0f21a84e..203963312 100644 --- a/qcportal/qcportal/neb/record_models.py +++ b/qcportal/qcportal/neb/record_models.py @@ -121,7 +121,6 @@ def _convert_basis(cls, v): class NEBRecord(BaseRecord): - record_type: Literal["neb"] = "neb" specification: NEBSpecification @@ -131,6 +130,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 ######################################## # Caches @@ -214,8 +214,6 @@ def _handle_includes(self, includes: Optional[Iterable[str]]): if "initial_chain" in includes: self._fetch_initial_chain() - if "final_chain" in includes: - self._fetch_singlepoints() if "singlepoints" in includes: self._fetch_singlepoints() if "optimizations" in includes: @@ -229,12 +227,26 @@ def initial_chain(self) -> List[Molecule]: @property def final_chain(self) -> List[SinglepointRecord]: + if self.optimizations.get("transition", None) is not None and self.ts_hessian_ is None: + self.singlepoints return self.singlepoints[max(self.singlepoints.keys())] + @property + def ts_hessian(self) -> Optional[SinglepointRecord]: + if self.optimizations.get("transition", None) is not None: + if self.ts_hessian_ is None: + self.singlepoints + return self.ts_hessian_ + else: + return None + @property def singlepoints(self) -> Dict[int, List[SinglepointRecord]]: if self.singlepoints_cache_ is None: self._fetch_singlepoints() + if self.optimizations.get("transition", None) is not None and self.ts_hessian_ is None: + _, temp_list = self.singlepoints_cache_.popitem() + self.ts_hessian_ = temp_list[0] return self.singlepoints_cache_ @property From aecc5b65d5ae529b535c5020ed219a4379e6bed3 Mon Sep 17 00:00:00 2001 From: Heejune Park Date: Fri, 20 Oct 2023 14:37:59 -0700 Subject: [PATCH 4/4] _fetch_singlepoints extracts the hessian SP record when ts optimization is requested --- qcportal/qcportal/neb/record_models.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/qcportal/qcportal/neb/record_models.py b/qcportal/qcportal/neb/record_models.py index 203963312..686d6c87a 100644 --- a/qcportal/qcportal/neb/record_models.py +++ b/qcportal/qcportal/neb/record_models.py @@ -193,6 +193,11 @@ 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_[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): @@ -227,26 +232,18 @@ def initial_chain(self) -> List[Molecule]: @property def final_chain(self) -> List[SinglepointRecord]: - if self.optimizations.get("transition", None) is not None and self.ts_hessian_ is None: - self.singlepoints return self.singlepoints[max(self.singlepoints.keys())] @property def ts_hessian(self) -> Optional[SinglepointRecord]: - if self.optimizations.get("transition", None) is not None: - if self.ts_hessian_ is None: - self.singlepoints - return self.ts_hessian_ - else: - return None + if self.singlepoints_cache_ is None: + self._fetch_singlepoints() + return self.ts_hessian_ @property def singlepoints(self) -> Dict[int, List[SinglepointRecord]]: if self.singlepoints_cache_ is None: self._fetch_singlepoints() - if self.optimizations.get("transition", None) is not None and self.ts_hessian_ is None: - _, temp_list = self.singlepoints_cache_.popitem() - self.ts_hessian_ = temp_list[0] return self.singlepoints_cache_ @property