From 07bad339a7953bd42e6eec2607956bdb718574f2 Mon Sep 17 00:00:00 2001 From: Colton Hicks Date: Tue, 1 Oct 2024 13:11:30 -0700 Subject: [PATCH] Added encoders and parsers for CREST to support energy, gradient, hessian, and optimization calculations. --- .vscode/settings.json | 2 + CHANGELOG.md | 4 + qcparse/encoders/crest.py | 60 ++++- qcparse/parsers/crest.py | 153 ++++++++++++- tests/data/crest_output/crest.engrad | 20 ++ tests/data/crest_output/crestopt.log | 65 ++++++ tests/data/crest_output/numhess1 | 20 ++ tests/data/crest_output/optstdout.txt | 304 ++++++++++++++++++++++++++ tests/test_crest.py | 171 +++++++++++++++ tests/test_main.py | 2 +- tests/test_models.py | 1 - 11 files changed, 794 insertions(+), 8 deletions(-) create mode 100644 tests/data/crest_output/crest.engrad create mode 100644 tests/data/crest_output/crestopt.log create mode 100644 tests/data/crest_output/numhess1 create mode 100644 tests/data/crest_output/optstdout.txt diff --git a/.vscode/settings.json b/.vscode/settings.json index 0772d54..10d735c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -20,6 +20,8 @@ "qcparse", "rotamer", "rotamers", + "runtypes", + "singlepoint", "spinmult", "tcin", "tcout", diff --git a/CHANGELOG.md b/CHANGELOG.md index f9c3b13..2053acc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ## [unreleased] +### Added + +- Encoders and parsers for CREST to support `energy`, `gradient`, `hessian`, and `optimization` calculations. + ## [0.6.3] - 2024-09-12 ### Added diff --git a/qcparse/encoders/crest.py b/qcparse/encoders/crest.py index c329c60..9ddc21b 100644 --- a/qcparse/encoders/crest.py +++ b/qcparse/encoders/crest.py @@ -8,7 +8,13 @@ from qcparse.exceptions import EncoderError from qcparse.models import NativeInput -SUPPORTED_CALCTYPES = {CalcType.conformer_search} +SUPPORTED_CALCTYPES = { + CalcType.conformer_search, + CalcType.optimization, + CalcType.energy, + CalcType.gradient, + CalcType.hessian, +} def encode(inp_obj: ProgramInput) -> NativeInput: @@ -42,13 +48,47 @@ def validate_input(inp_obj: ProgramInput): """ # These values come from other parts of the ProgramInput and should not be set # in the keywords. - non_allowed_keywords = ["charge", "uhf", "runtype"] + non_allowed_keywords = ["charge", "uhf"] for keyword in non_allowed_keywords: if keyword in inp_obj.keywords: raise EncoderError( f"{keyword} should not be set in keywords for CREST. It is already set " "on the Structure or ProgramInput elsewhere.", ) + if "runtype" in inp_obj.keywords: + _validate_runtype_calctype(inp_obj.keywords["runtype"], inp_obj.calctype) + + +def _validate_runtype_calctype(runtype: str, calctype: CalcType): + """Validate that the runtype is supported for the calctype.""" + invalid_runtype = False + valid_runtypes = set() + + if calctype == CalcType.conformer_search: + valid_runtypes = {"imtd-gc", "imtd-smtd", "entropy", "nci", "nci-mtd"} + if runtype not in valid_runtypes: + invalid_runtype = True + + elif calctype == CalcType.optimization: + valid_runtypes = {"optimize", "ancopt"} + if runtype not in valid_runtypes: + invalid_runtype = True + + elif calctype in {CalcType.energy, CalcType.gradient}: + valid_runtypes = {"singlepoint"} + if runtype not in valid_runtypes: + invalid_runtype = True + + elif calctype == CalcType.hessian: + valid_runtypes = {"numhess"} + if runtype not in valid_runtypes: + invalid_runtype = True + + if invalid_runtype: + raise EncoderError( + f"Unsupported runtype {runtype} for calctype {calctype}. Valid runtypes " + f"are: {valid_runtypes}.", + ) def _to_toml_dict(inp_obj: ProgramInput, struct_filename: str) -> Dict[str, Any]: @@ -64,8 +104,20 @@ def _to_toml_dict(inp_obj: ProgramInput, struct_filename: str) -> Dict[str, Any] toml_dict.setdefault("threads", os.cpu_count()) toml_dict["input"] = struct_filename - # TODO: May need to deal with non-covalent mode at some point - toml_dict["runtype"] = "imtd-gc" + # Set default runtype if not already set + if "runtype" not in inp_obj.keywords: + if inp_obj.calctype == CalcType.conformer_search: + toml_dict["runtype"] = "imtd-gc" + elif inp_obj.calctype == CalcType.optimization: + toml_dict["runtype"] = "optimize" + elif inp_obj.calctype in {CalcType.energy, CalcType.gradient}: + toml_dict["runtype"] = "singlepoint" + elif inp_obj.calctype == CalcType.hessian: + toml_dict["runtype"] = "numhess" + else: + raise EncoderError( + f"Unsupported calctype {inp_obj.calctype} for CREST encoder.", + ) # Calculation level keywords calculation = toml_dict.pop("calculation", {}) diff --git a/qcparse/parsers/crest.py b/qcparse/parsers/crest.py index 6f0fa06..630bf67 100644 --- a/qcparse/parsers/crest.py +++ b/qcparse/parsers/crest.py @@ -1,8 +1,18 @@ +import re from pathlib import Path -from typing import List, Optional, Union +from typing import Any, Dict, List, Optional, Union import numpy as np -from qcio import ConformerSearchResults, Structure +from qcio import ( + CalcType, + ConformerSearchResults, + OptimizationResults, + ProgramInput, + ProgramOutput, + Provenance, + SinglePointResults, + Structure, +) from .utils import regex_search @@ -87,3 +97,142 @@ def parse_conformer_search_dir( rotamers=rotamers, rotamer_energies=np.array(rotamer_energies), ) + + +def parse_energy_grad(text: str) -> SinglePointResults: + """Parse the output of a CREST energy and gradient calculation. + + Args: + text: The text of the output file. + + Returns: + The parsed energy and gradient as a SinglePointResults object. + """ + # Parse the energy + energy_regex = r"# Energy \( Eh \)\n#*\n\s*([-\d.]+)" + gradient_regex = r"# Gradient \( Eh/a0 \)\n#\s*\n((?:\s*[-\d.]+\n)+)" + + energy = float(regex_search(energy_regex, text).group(1)) + gradient = np.array( + [float(x) for x in regex_search(gradient_regex, text).group(1).split()] + ) + return SinglePointResults( + energy=energy, + gradient=gradient, + ) + + +def parse_singlepoint_dir( + directory: Union[Path, str], filename: str = "crest.engrad" +) -> SinglePointResults: + """Parse the output directory of a CREST single point calculation. + + Args: + directory: Path to the directory containing the CREST output files. + filename: The name of the file containing the single point results. + Default is 'crest.engrad'. + + Returns: + The parsed single point results as a SinglePointResults object. + """ + directory = Path(directory) + text = (directory / filename).read_text() + + return parse_energy_grad(text) + + +def parse_numhess_dir( + directory: Union[Path, str], + filename: str = "numhess1", + stdout: Optional[str] = None, +) -> SinglePointResults: + """Parse the output directory of a CREST numerical Hessian calculation. + + Args: + directory: Path to the directory containing the CREST output files. + filename: The name of the file containing the numerical Hessian results. + Default is 'numhess1'. + + Returns: + The parsed numerical Hessian results as a SinglePointResults object. + """ + data = (Path(directory) / filename).read_text() + float_regex = r"[-+]?\d*\.\d+|\d+" + numbers = re.findall(float_regex, data) + array = np.array(numbers, dtype=float) + spr_dict: Dict[str, Any] = {"hessian": array} + if stdout: + energy_regex = r"Energy\s*=\s*([-+]?\d+\.\d+)\s*Eh" + energy = float(regex_search(energy_regex, stdout).group(1)) + spr_dict["energy"] = energy + return SinglePointResults(**spr_dict) + + +def parse_optimization_dir( + directory: Union[Path, str], + *, + inp_obj: ProgramInput, + stdout: str, +) -> OptimizationResults: + """Parse the output directory of a CREST optimization calculation. + + Args: + directory: Path to the directory containing the CREST output files. + inp_obj: The qcio ProgramInput object for the optimization. + stdout: The stdout from CREST. + + Returns: + The parsed optimization results as a OptimizationResults object. + """ + # Read in the xyz file containing the trajectory + directory = Path(directory) + xyz_text = (directory / "crestopt.log").read_text() + + # Parse structures and energies from the xyz file + structures = Structure.from_xyz_multi( + xyz_text, + charge=inp_obj.structure.charge, + multiplicity=inp_obj.structure.multiplicity, + ) + energies = [ + float(struct.extras[Structure._xyz_comment_key][1]) for struct in structures + ] + + # Fake gradient for each step because CREST does not output it + fake_gradient = np.zeros(len(inp_obj.structure.symbols) * 3) + + # Parse program version + program_version = parse_version_string(stdout) + + # Collect final gradient if calculation succeeded + try: + final_spr = parse_singlepoint_dir(directory) + except FileNotFoundError: + # Calculation failed, so we don't have the final energy or gradient + final_spr = SinglePointResults(gradient=fake_gradient) + + # Create the optimization trajectory + trajectory: List[ProgramOutput] = [ + ProgramOutput( + input_data=ProgramInput( + calctype=CalcType.gradient, + structure=struct, + model=inp_obj.model, + ), + success=True, + results=SinglePointResults(energy=energy, gradient=fake_gradient), + provenance=Provenance( + program="crest", + program_version=program_version, + ), + ) + for struct, energy in zip(structures, energies) + ] + + # Fill in final gradient + # https://github.com/crest-lab/crest/issues/354 + trajectory[-1].results.gradient[:] = final_spr.gradient + + return OptimizationResults( + trajectory=trajectory, + ) diff --git a/tests/data/crest_output/crest.engrad b/tests/data/crest_output/crest.engrad new file mode 100644 index 0000000..89dd97f --- /dev/null +++ b/tests/data/crest_output/crest.engrad @@ -0,0 +1,20 @@ +# +# Atoms +# + 3 +# +# Energy ( Eh ) +# + -0.335557824179335 +# +# Gradient ( Eh/a0 ) +# + -0.005962071557911 + -0.004419818102026 + 0.003139227894649 + 0.003048425211480 + 0.001982394235964 + -0.001779667371498 + 0.002913646346432 + 0.002437423866062 + -0.001359560523152 diff --git a/tests/data/crest_output/crestopt.log b/tests/data/crest_output/crestopt.log new file mode 100644 index 0000000..4d9345e --- /dev/null +++ b/tests/data/crest_output/crestopt.log @@ -0,0 +1,65 @@ + 3 + Etot= -4.7918798035 + O -0.0934852751 -0.0692099762 0.0488995019 + H 0.4307247549 1.6181264838 0.5296458319 + H 1.0532005349 -0.5195317162 -1.3058451581 + 3 + Etot= -4.9229264187 + O 0.1136777730 0.0841383057 -0.0594737049 + H 0.3707615594 1.3916895508 0.4552315157 + H 0.9060006825 -0.4464430651 -1.1230576351 + 3 + Etot= -5.0483521241 + O 0.2573182669 0.1904912574 -0.1346012710 + H 0.3599743190 1.1289027403 0.3128185924 + H 0.7731474289 -0.2900092063 -0.9055171457 + 3 + Etot= -5.0170597610 + O 0.1917723991 0.1420049385 -0.1002932645 + H 0.4678576011 0.8951550442 0.0741564841 + H 0.7308100147 -0.0077751914 -0.7011630439 + 3 + Etot= -5.0491088307 + O 0.0623819936 0.0460852580 -0.0326872105 + H 0.5219250724 0.9797836792 0.0717912538 + H 0.8061329488 0.0035158541 -0.7664038677 + 3 + Etot= -5.0660326493 + O 0.1987775715 0.1472783377 -0.1039067694 + H 0.4093039819 1.0814532336 0.2382238721 + H 0.7823584615 -0.1993467800 -0.8616169270 + 3 + Etot= -5.0730217268 + O 0.1770402418 0.1302188641 -0.0930982848 + H 0.4376238607 1.0313507655 0.1821152037 + H 0.7757759123 -0.1321848383 -0.8163167432 + 3 + Etot= -5.0723078880 + O 0.1736292273 0.1445536976 -0.0815172819 + H 0.4380419413 1.0059759479 0.1669722835 + H 0.7787688463 -0.1211448541 -0.8127548260 + 3 + Etot= -5.0683269168 + O 0.1851362467 0.1179180760 -0.1079633483 + H 0.4523839201 0.9974657046 0.1483557889 + H 0.7529198480 -0.0859989893 -0.7676922650 + 3 + Etot= -5.0732163707 + O 0.1760746513 0.1281388107 -0.0933864784 + H 0.4401853431 1.0270846083 0.1771945599 + H 0.7741800205 -0.1258386276 -0.8111079058 + 3 + Etot= -5.0733841586 + O 0.1769285244 0.1328773064 -0.0914470659 + H 0.4405581451 1.0164576782 0.1706642632 + H 0.7729533454 -0.1199501932 -0.8065170216 + 3 + Etot= -5.0734021646 + O 0.1782533830 0.1321019587 -0.0931605147 + H 0.4400623060 1.0186302961 0.1723993514 + H 0.7721243259 -0.1213474635 -0.8065386611 + 3 + Etot= -5.0734025156 + O 0.1788766949 0.1323436930 -0.0936142242 + H 0.4397631667 1.0187613429 0.1727606527 + H 0.7718001532 -0.1217202445 -0.8064462529 diff --git a/tests/data/crest_output/numhess1 b/tests/data/crest_output/numhess1 new file mode 100644 index 0000000..dcb39a6 --- /dev/null +++ b/tests/data/crest_output/numhess1 @@ -0,0 +1,20 @@ + $hessian + 0.02040569 -0.00018059 0.02080099 -0.02081319 0.01511689 + 0.00867078 0.00037976 -0.01495837 -0.02946283 + -0.00018059 -0.01341723 -0.03209513 0.01368595 0.03374600 + 0.01874084 -0.01351862 -0.02035995 0.01336374 + 0.02080099 -0.03209513 0.00327178 0.00784908 0.01737681 + -0.01812512 -0.02863169 0.01472059 0.01485103 + -0.02081319 0.01368595 0.00784908 0.01933555 -0.01625843 + -0.00694960 0.00149263 0.00258608 -0.00090575 + 0.01511689 0.03374600 0.01737681 -0.01625843 -0.03409225 + -0.01710500 0.00114214 0.00035657 -0.00027546 + 0.00867078 0.01874084 -0.01812512 -0.00694960 -0.01710500 + 0.01843539 -0.00173455 -0.00164242 -0.00030677 + 0.00037976 -0.01351862 -0.02863169 0.00149263 0.00114214 + -0.00173455 -0.00185964 0.01238496 0.03036359 + -0.01495837 -0.02035995 0.01472059 0.00258608 0.00035657 + -0.00164242 0.01238496 0.02002423 -0.01308397 + -0.02946283 0.01336374 0.01485103 -0.00090575 -0.00027546 + -0.00030677 0.03036359 -0.01308397 -0.01454546 + $end diff --git a/tests/data/crest_output/optstdout.txt b/tests/data/crest_output/optstdout.txt new file mode 100644 index 0000000..53de2e2 --- /dev/null +++ b/tests/data/crest_output/optstdout.txt @@ -0,0 +1,304 @@ + + ╔════════════════════════════════════════════╗ + ║ ___ ___ ___ ___ _____ ║ + ║ / __| _ \ __/ __|_ _| ║ + ║ | (__| / _|\__ \ | | ║ + ║ \___|_|_\___|___/ |_| ║ + ║ ║ + ║ Conformer-Rotamer Ensemble Sampling Tool ║ + ║ based on the xTB methods ║ + ║ ║ + ╚════════════════════════════════════════════╝ + Version 3.0.2, Sun, 25 August 20:02:44, 08/25/2024 + commit (af7eb99) compiled by 'usr@fv-az732-492' + + Cite work conducted with this code as + + • P.Pracht, F.Bohle, S.Grimme, PCCP, 2020, 22, 7169-7192. + • S.Grimme, JCTC, 2019, 15, 2847-2862. + • P.Pracht, S.Grimme, C.Bannwarth, F.Bohle, S.Ehlert, + G.Feldmann, J.Gorges, M.Müller, T.Neudecker, C.Plett, + S.Spicher, P.Steinbach, P.Wesołowski, F.Zeller, + J. Chem. Phys., 2024, 160, 114110. + + for works involving QCG cite + + • S.Spicher, C.Plett, P.Pracht, A.Hansen, S.Grimme, + JCTC, 2022, 18 (5), 3174-3189. + • C.Plett, S. Grimme, + Angew. Chem. Int. Ed. 2023, 62, e202214477. + + for works involving MECP screening cite + + • P.Pracht, C.Bannwarth, JCTC, 2022, 18 (10), 6370-6385. + + Original code + P.Pracht, S.Grimme, Universität Bonn, MCTC + with help from (alphabetical order): + C.Bannwarth, F.Bohle, S.Ehlert, G.Feldmann, J.Gorges, + S.Grimme, C.Plett, P.Pracht, S.Spicher, P.Steinbach, + P.Wesolowski, F.Zeller + + Online documentation is available at + https://crest-lab.github.io/crest-docs/ + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License (LGPL) for more details. + + Command line input: + $ crest input.toml + + reading input.toml +******************************************************************************** +* INPUT FILE input.toml content (without comments): +******************************************************************************** +* threads = 16 +* input = "structure.xyz" +* runtype = "optimize" +* [calculation] +* [[calculation.level]] +* alpb = "ch2cl2" +* method = "gfn2" +* charge = 0 +* uhf = 0 +******************************************************************************** + ============================= + # threads = 16 + ============================= + + ---------------- + Input structure: + ---------------- + 3 + + O -0.0934852751 -0.0692099762 0.0488995019 + H 0.4307247549 1.6181264838 0.5296458319 + H 1.0532005349 -0.5195317162 -1.3058451581 + + ---------------- + Calculation info + ---------------- +> User-defined calculation level: + : xTB calculation via tblite lib + : GFN2-xTB level + : Molecular charge : 0 + : Solvation model : alpb + : Solvent : ch2cl2 + : Fermi temperature : 300.00000 + : Accuracy : 1.00000 + : max SCC cycles : 500 + +-------------------------------------------------------------------------------- + + ┍━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┑ + │ GEOMETRY OPTIMIZATION SETUP │ + ┝━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┥ + │ algorithm ANCOPT │ + │ optimization level normal (0) │ + │ max. optcycles 200 │ + │ ANC micro-cycles 20 │ + │ degrees of freedom 3 │ + ├─────────────────────────────────────────────────┤ + │ RF solver davidson │ + │ Hessian update bfgs │ + │ write crestopt.log True │ + │ linear? False │ + │ energy convergence 0.5000000E-05 Eh │ + │ grad. convergence 0.1000000E-02 Eh/α │ + │ maximium RF displ. 1.0000000 │ + │ Hlow (freq-cutoff) 0.1000000E-01 │ + │ Hmax (freq-cutoff) 5.0000000 │ + │ S6 in model hess. 20.0000000 │ + └─────────────────────────────────────────────────┘ + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 0 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -4.7918798 Eh change ΔE 0.0000000E+00 Eh + gradient norm : 0.1394460 Eh/a0 predicted 0.0000000E+00 ( -0.00%) + +generating ANC from model Hessian ... +Using Lindh-Hessian (1995) + Shifting diagonal of input Hessian by 9.4088233643030428E-003 + Lowest eigenvalues of input Hessian + 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 + 0.010000 0.048000 0.058377 + Highest eigenvalues + 0.000000 0.000000 0.000000 0.010000 0.048000 0.058377 + + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 1 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -4.7918798 Eh change ΔE -0.3714957E-09 Eh + gradient norm : 0.1394461 Eh/α predicted 0.0000000E+00 (-100.00%) + displ. norm : 0.8430058 α lambda -0.1175421E+00 + maximum displ.: 0.8414127 α in ANC's #2, #1, #3, ... + converged δE/grad : True / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 2 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -4.9229264 Eh change ΔE -0.1310466E+00 Eh + gradient norm : 0.1639971 Eh/α predicted -0.1005373E+00 ( -23.28%) + displ. norm : 0.8792564 α lambda -0.1438084E+00 + maximum displ.: 0.8094865 α in ANC's #2, #1, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 3 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0483521 Eh change ΔE -0.1254257E+00 Eh + gradient norm : 0.0961023 Eh/α predicted -0.1274927E+00 ( 1.65%) + displ. norm : 0.9528974 α lambda -0.9010465E-01 + maximum displ.: 0.8709945 α in ANC's #1, #2, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 4 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0170598 Eh change ΔE 0.3129236E-01 Eh + gradient norm : 0.4380398 Eh/α predicted -0.8596044E-01 (-374.70%) + displ. norm : 0.4253409 α lambda -0.1732278E+00 + maximum displ.: 0.4039886 α in ANC's #2, #1, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 5 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0491088 Eh change ΔE -0.3204907E-01 Eh + gradient norm : 0.1338961 Eh/α predicted -0.1022837E+00 ( 219.15%) + displ. norm : 0.6954090 α lambda -0.5818773E-01 + maximum displ.: 0.6951811 α in ANC's #1, #2, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 6 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0660326 Eh change ΔE -0.1692382E-01 Eh + gradient norm : 0.0688184 Eh/α predicted -0.4316352E-01 ( 155.05%) + displ. norm : 0.2231716 α lambda -0.1146472E-01 + maximum displ.: 0.2117156 α in ANC's #1, #2, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 7 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0730217 Eh change ΔE -0.6989077E-02 Eh + gradient norm : 0.0221586 Eh/α predicted -0.6017861E-02 ( -13.90%) + displ. norm : 0.0698977 α lambda -0.9518713E-03 + maximum displ.: 0.0475382 α in ANC's #3, #1, #2, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 8 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0723079 Eh change ΔE 0.7138388E-03 Eh + gradient norm : 0.0486061 Eh/α predicted -0.4782619E-03 (-167.00%) + displ. norm : 0.1475922 α lambda -0.4204657E-02 + maximum displ.: 0.0991230 α in ANC's #3, #1, #2, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 9 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0683269 Eh change ΔE 0.3980971E-02 Eh + gradient norm : 0.1092617 Eh/α predicted -0.2148123E-02 (-153.96%) + displ. norm : 0.1485682 α lambda -0.1221386E-01 + maximum displ.: 0.1033676 α in ANC's #1, #2, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 10 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0732164 Eh change ΔE -0.4889454E-02 Eh + gradient norm : 0.0167703 Eh/α predicted -0.6241725E-02 ( 27.66%) + displ. norm : 0.0292702 α lambda -0.4059214E-03 + maximum displ.: 0.0201853 α in ANC's #2, #1, #3, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 11 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0733842 Eh change ΔE -0.1677879E-03 Eh + gradient norm : 0.0056902 Eh/α predicted -0.2031358E-03 ( 21.07%) + displ. norm : 0.0075374 α lambda -0.3330981E-04 + maximum displ.: 0.0056458 α in ANC's #1, #3, #2, ... + converged δE/grad : False / False + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 12 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0734022 Eh change ΔE -0.1800607E-04 Eh + gradient norm : 0.0007549 Eh/α predicted -0.1665551E-04 ( -7.50%) + displ. norm : 0.0020200 α lambda -0.7920827E-06 + maximum displ.: 0.0016363 α in ANC's #1, #2, #3, ... + converged δE/grad : False / True + +┌────────────────────────────────────────────────────────────────────────────┐ +│ CYCLE 13 │ +└────────────────────────────────────────────────────────────────────────────┘ + * total energy : -5.0734025 Eh change ΔE -0.3509186E-06 Eh + gradient norm : 0.0003506 Eh/α predicted -0.3957971E-06 ( 12.79%) + displ. norm : 0.0006684 α lambda -0.1013195E-06 + maximum displ.: 0.0004324 α in ANC's #1, #2, #3, ... + converged δE/grad : True / True + + *** GEOMETRY OPTIMIZATION CONVERGED AFTER 13 ITERATIONS *** + +------------------------------------------------------------------------ + total energy gain : -0.2815227 Eh -176.6582 kcal/mol + total RMSD : 1.1486983 a0 0.6079 Å +------------------------------------------------------------------------ + SUCCESS! + geometry successfully optimized! + +-------------------------------------------------------------------------------- + + FINAL CALCULATION SUMMARY + ========================= + +> Final molecular gradient ( Eh/a0 ): + ∂E/∂x ∂E/∂y ∂E/∂z + 0.00017271 -0.00003056 -0.00019257 + -0.00001068 0.00000438 0.00002687 + -0.00016204 0.00002618 0.00016570 +> Gradient norm: 0.00035086 Eh/a0 + +> Writing crest.engrad ... done. + +-------------------------------------------------------------------------------- + + ----------------- + Output structure: + ----------------- + 3 + + O 0.17866812 0.13231044 -0.09343472 + H 0.43983852 1.01871471 0.17266172 + H 0.77193337 -0.12164035 -0.80652682 + +> Optimized geometry written to crestopt.xyz + +========================================== + TOTAL ENERGY -5.0734025156 Eh + GRADIENT NORM 0.0003508607 Eh/a0 +========================================== + + ----------------- + Wall Time Summary + ----------------- + CREST runtime (total) 0 d, 0 h, 0 min, 0.145 sec + ------------------------------------------------------------------ + Geometry optimization ... 0 min, 0.069 sec ( 47.829%) + I/O and setup ... 0 min, 0.075 sec ( 52.171%) + ------------------------------------------------------------------ + * wall-time: 0 d, 0 h, 0 min, 0.145 sec + * cpu-time: 0 d, 0 h, 0 min, 2.048 sec + * ratio c/w: 14.155 speedup + ------------------------------------------------------------------ + * Total number of energy+grad calls: 14 + + CREST terminated normally. diff --git a/tests/test_crest.py b/tests/test_crest.py index 7178adf..dc9ef71 100644 --- a/tests/test_crest.py +++ b/tests/test_crest.py @@ -1,3 +1,4 @@ +import numpy as np import pytest from qcio import ProgramInput, Structure from qcio.utils import water @@ -6,6 +7,10 @@ from qcparse.exceptions import EncoderError from qcparse.parsers.crest import ( parse_conformer_search_dir, + parse_energy_grad, + parse_numhess_dir, + parse_optimization_dir, + parse_singlepoint_dir, parse_structures, parse_version_string, ) @@ -110,3 +115,169 @@ def test_parse_conformer_search_charge_multiplicity(test_data_dir): for struct in getattr(csr, struct_type): assert struct.charge == -2 assert struct.multiplicity == 3 + + +def test_parse_energy_grad(test_data_dir): + text = (test_data_dir / "crest_output" / "crest.engrad").read_text() + spr = parse_energy_grad(text) + assert spr.energy == -0.335557824179335 + np.testing.assert_array_equal( + spr.gradient, + [ + [-0.005962071557911, -0.004419818102026, 0.003139227894649], + [0.003048425211480, 0.001982394235964, -0.001779667371498], + [0.002913646346432, 0.002437423866062, -0.001359560523152], + ], + ) + + +def test_parse_singlepoint_dir(test_data_dir): + spr = parse_singlepoint_dir(test_data_dir / "crest_output") + assert spr.energy == -0.335557824179335 + np.testing.assert_array_equal( + spr.gradient, + [ + [-0.005962071557911, -0.004419818102026, 0.003139227894649], + [0.003048425211480, 0.001982394235964, -0.001779667371498], + [0.002913646346432, 0.002437423866062, -0.001359560523152], + ], + ) + + +def test_parse_optimization_dir(test_data_dir, prog_inp): + prog_input = prog_inp("optimization") + stdout = (test_data_dir / "crest_output" / "optstdout.txt").read_text() + opt_res = parse_optimization_dir( + test_data_dir / "crest_output", inp_obj=prog_input, stdout=stdout + ) + assert len(opt_res.trajectory) == 13 + + # Fills in the final gradient correctly with values from crest.engrad files + np.testing.assert_array_equal( + opt_res.trajectory[-1].results.gradient, + [ + [-0.005962071557911, -0.004419818102026, 0.003139227894649], + [0.003048425211480, 0.001982394235964, -0.001779667371498], + [0.002913646346432, 0.002437423866062, -0.001359560523152], + ], + ) + + +def test_parse_numhess_dir(test_data_dir): + stdout = """ + :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + Initial singlpoint calculation ... + ------------------------------------- + + Energy = -110.490788960984773 Eh + ------------------------------------- + :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + Calculating numerical Hessian ... done. + """ + spr = parse_numhess_dir(test_data_dir / "crest_output", stdout=stdout) + assert spr.energy == -110.490788960984773 + np.testing.assert_array_equal( + spr.hessian, + [ + [ + 0.02040569, + -0.00018059, + 0.02080099, + -0.02081319, + 0.01511689, + 0.00867078, + 0.00037976, + -0.01495837, + -0.02946283, + ], + [ + -0.00018059, + -0.01341723, + -0.03209513, + 0.01368595, + 0.03374600, + 0.01874084, + -0.01351862, + -0.02035995, + 0.01336374, + ], + [ + 0.02080099, + -0.03209513, + 0.00327178, + 0.00784908, + 0.01737681, + -0.01812512, + -0.02863169, + 0.01472059, + 0.01485103, + ], + [ + -0.02081319, + 0.01368595, + 0.00784908, + 0.01933555, + -0.01625843, + -0.00694960, + 0.00149263, + 0.00258608, + -0.00090575, + ], + [ + 0.01511689, + 0.03374600, + 0.01737681, + -0.01625843, + -0.03409225, + -0.01710500, + 0.00114214, + 0.00035657, + -0.00027546, + ], + [ + 0.00867078, + 0.01874084, + -0.01812512, + -0.00694960, + -0.01710500, + 0.01843539, + -0.00173455, + -0.00164242, + -0.00030677, + ], + [ + 0.00037976, + -0.01351862, + -0.02863169, + 0.00149263, + 0.00114214, + -0.00173455, + -0.00185964, + 0.01238496, + 0.03036359, + ], + [ + -0.01495837, + -0.02035995, + 0.01472059, + 0.00258608, + 0.00035657, + -0.00164242, + 0.01238496, + 0.02002423, + -0.01308397, + ], + [ + -0.02946283, + 0.01336374, + 0.01485103, + -0.00090575, + -0.00027546, + -0.00030677, + 0.03036359, + -0.01308397, + -0.01454546, + ], + ], + ) diff --git a/tests/test_main.py b/tests/test_main.py index e0179ae..78d9f42 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -11,7 +11,7 @@ def test_main_terachem_energy(terachem_energy_stdout): def test_encode_raises_error_with_invalid_calctype(prog_inp): - prog_inp = prog_inp("optimization") # Not currently supported by terachem encoder + prog_inp = prog_inp("transition_state") # Not currently supported by crest encoder with pytest.raises(EncoderError): encode(prog_inp, "crest") diff --git a/tests/test_models.py b/tests/test_models.py index 4a5add9..b7ec82b 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -25,7 +25,6 @@ def test_pdc_dict_serialization_method(): pdc.nested_list[1].value = 3 pdc_dict = pdc.dict() - # import pdb; pdb.set_trace() assert pdc_dict["value"] == 1 assert pdc_dict["nested"]["energy"] == -74.964814 assert pdc_dict["nested_list"][0]["value"] == 2