From 3fbfbfa2f9682a72cd5deb4073b4e178d297c88b Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Mon, 9 Feb 2026 13:29:32 +0000 Subject: [PATCH 1/6] Add ability to create and parse vibrational frequency jobs with ChemShell --- src/aiida_chemshell/calculations/base.py | 9 +++++++- src/aiida_chemshell/parsers/base.py | 28 +++++++++++++++++++++--- tests/test_calculations.py | 24 ++++++++++++++++++++ 3 files changed, 57 insertions(+), 4 deletions(-) diff --git a/src/aiida_chemshell/calculations/base.py b/src/aiida_chemshell/calculations/base.py index eb0cee7..ddfc427 100644 --- a/src/aiida_chemshell/calculations/base.py +++ b/src/aiida_chemshell/calculations/base.py @@ -3,7 +3,7 @@ from aiida.common import CalcInfo, CodeInfo from aiida.common.folders import Folder from aiida.engine import CalcJob, CalcJobProcessSpec -from aiida.orm import ArrayData, Dict, Float, SinglefileData, StructureData +from aiida.orm import ArrayData, Dict, Float, SinglefileData, Str, StructureData from aiida_chemshell.utils import ChemShellMMTheory, ChemShellQMTheory @@ -131,6 +131,12 @@ def define(cls, spec: CalcJobProcessSpec) -> None: "structure is contained within a ChemShell '.pun' file." ), ) + spec.output( + "vibrational_analysis", + valid_type=Str, + required=False, + help="The vibrational analysis from ChemShell/DL_FIND", + ) ## Metadata spec.inputs["metadata"]["options"]["resources"].default = { @@ -282,6 +288,7 @@ def get_valid_optimisation_parameter_keys(cls) -> tuple[str]: "dimer", "delta", "tsrelative", + "thermal", ) @classmethod diff --git a/src/aiida_chemshell/parsers/base.py b/src/aiida_chemshell/parsers/base.py index 5139214..9f687c9 100644 --- a/src/aiida_chemshell/parsers/base.py +++ b/src/aiida_chemshell/parsers/base.py @@ -5,7 +5,7 @@ import numpy from aiida.common import ModificationNotAllowed from aiida.engine import ExitCode -from aiida.orm import ArrayData, Float, SinglefileData +from aiida.orm import ArrayData, Float, SinglefileData, Str from aiida.parsers.parser import Parser from aiida_chemshell.calculations.base import ChemShellCalculation @@ -23,7 +23,7 @@ def parse(self, **kwargs): # Read the 'json' formatted results file results = json.loads( - self.retrieved.get_object_content(ChemShellCalculation.FILE_RESULTS) + self.retrieved.get_object_content(ChemShellCalculation.FILE_RESULTS, "rb") ) # Extract the final energy @@ -67,7 +67,13 @@ def parse(self, **kwargs): # If the calculation was a geometry optimisation, store the optimised structure if "optimisation_parameters" in self.node.inputs: - if ChemShellCalculation.FILE_DLFIND in self.retrieved.list_object_names(): + if self.node.inputs.optimisation_parameters.get("thermal", False): + self.parse_vibrational_analysis( + self.retrieved.get_object_content( + ChemShellCalculation.FILE_STDOUT, "r" + ) + ) + elif ChemShellCalculation.FILE_DLFIND in self.retrieved.list_object_names(): descrip = "Optimised structure from a ChemShell optimisation" input_pk = self.node.inputs.structure.pk descrip += f" of node {input_pk}" @@ -89,3 +95,19 @@ def parse(self, **kwargs): return self.exit_codes.ERROR_MISSING_OPTIMISED_STRUCTURE_FILE return ExitCode(0) + + def parse_vibrational_analysis(self, stdout: str) -> None: + """Extract the vibrational analysis from ChemShell output log.""" + read = False + thermo_analysis = [] + for line in stdout.split("\n"): + if "Thermochemical analysis" in line: + read = True + elif "total S vib" in line: + thermo_analysis.append(line) + read = False + + if read: + thermo_analysis.append(line) + self.out("vibrational_analysis", Str("\n".join(thermo_analysis))) + return diff --git a/tests/test_calculations.py b/tests/test_calculations.py index d95a06e..89359eb 100644 --- a/tests/test_calculations.py +++ b/tests/test_calculations.py @@ -216,6 +216,30 @@ def test_opt_calculation_dlpoly(chemsh_code, get_test_data_file): ) +def test_vibrational_calculation(chemsh_code, get_test_data_file): + """MM based geometry optimisation test.""" + code = chemsh_code() + builder = code.get_builder() + builder.structure = get_test_data_file() + builder.qm_parameters = Dict({"theory": "PySCF", "method": "hf", "basis": "3-21G"}) + builder.optimisation_parameters = Dict({"thermal": True}) + + results, node = run.get_node(builder) + + assert node.is_finished_ok, "CalcJob failed for `test_OptCalculation_dlpoly`" + + thermo_analysis = results.get("vibrational_analysis").value + + assert "Thermochemical analysis" in thermo_analysis + assert ( + "0.1252886592 1819.534 2617.900 0.0041452033 0.0000013454 -0.0000014996" + in thermo_analysis + ) + assert "total ZPE 58728.93143 J/mol" in thermo_analysis + assert "total E vib 3.53277 J/mol" in thermo_analysis + assert "total S vib 0.01313 J/mol/K" in thermo_analysis + + # def test_opt_calculation_qmmm(chemsh_code, get_test_data_file): # """QM/MM geometry optimisation test.""" # code = chemsh_code() From 967cdff6dffc7d5cf7480a87f1a40d6b3acab233 Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Mon, 9 Feb 2026 16:20:42 +0000 Subject: [PATCH 2/6] Create an overall input validator for checking if necessary corresponding inputs have been given --- src/aiida_chemshell/calculations/base.py | 34 ++++++++++- tests/conftest.py | 2 +- tests/test_error_codes.py | 73 ++++++++++++++++++++++++ 3 files changed, 107 insertions(+), 2 deletions(-) diff --git a/src/aiida_chemshell/calculations/base.py b/src/aiida_chemshell/calculations/base.py index ddfc427..cd642a8 100644 --- a/src/aiida_chemshell/calculations/base.py +++ b/src/aiida_chemshell/calculations/base.py @@ -2,7 +2,7 @@ from aiida.common import CalcInfo, CodeInfo from aiida.common.folders import Folder -from aiida.engine import CalcJob, CalcJobProcessSpec +from aiida.engine import CalcJob, CalcJobProcessSpec, PortNamespace from aiida.orm import ArrayData, Dict, Float, SinglefileData, Str, StructureData from aiida_chemshell.utils import ChemShellMMTheory, ChemShellQMTheory @@ -138,6 +138,22 @@ def define(cls, spec: CalcJobProcessSpec) -> None: help="The vibrational analysis from ChemShell/DL_FIND", ) + # Validate inputs namespace + existing_validator = spec.inputs.validator + + def inputs_validator_wrapper(inputs, namespace): + """Wrap the existing and custom input namespace validators.""" + if existing_validator: + error = existing_validator(inputs, namespace) + if error: + return error + error = cls.validate_inputs_namespace(inputs, namespace) + if error: + return error + return None + + spec.inputs.validator = inputs_validator_wrapper + ## Metadata spec.inputs["metadata"]["options"]["resources"].default = { "num_machines": 1, @@ -184,6 +200,22 @@ def define(cls, spec: CalcJobProcessSpec) -> None: return + @classmethod + def validate_inputs_namespace( + cls, value: Dict, namespace: PortNamespace + ) -> str | None: + """Perform additional validation checks on the total inputs namespace.""" + if "mm_parameters" in value and "force_field_file" not in value: + return "A force field must be specified to use molecular mechanics." + if "force_field_file" in value and "mm_parameters" not in value: + return "A MM theory code must be specified to use molecular mechanics." + if "qmmm_parameters" in value: + if "qm_parameters" not in value: + return "Missing QM parameters for QM/MM calculation" + if "mm_parameters" not in value: + return "Missing MM parameters for QM/MM calculation" + return None + @classmethod def validate_structure_file( cls, value: SinglefileData | StructureData | None, _ diff --git a/tests/conftest.py b/tests/conftest.py index 2c30111..10accce 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -69,7 +69,7 @@ def factory( structure_fname: str | StructureData = "water.cjson", ff_fname: str | None = None, opt: dict | None = None, - ): + ) -> dict: if isinstance(structure_fname, str): structure = get_test_data_file(structure_fname) else: diff --git a/tests/test_error_codes.py b/tests/test_error_codes.py index b4d6656..7a228ba 100644 --- a/tests/test_error_codes.py +++ b/tests/test_error_codes.py @@ -37,6 +37,43 @@ def test_mm_theory_validation(generate_calcjob, generate_inputs): raise AssertionError("No error caught when providing invalid MM theory key.") +def test_missing_force_field_file(generate_calcjob, generate_inputs): + """Test the co-validation of mm_parameters and force_field_file (force field).""" + inputs = generate_inputs(mm={"theory": "DL_POLY"}, structure_fname="butanol.cjson") + try: + generate_calcjob(ChemShellCalculation, inputs) + except ValueError as e: + message = str(e) + assert ( + "A force field must be specified to use molecular mechanics." in message + ), "Wrong error message triggered for missing force field validation." + except Exception as e: + raise AssertionError( + f"Wrong error code generated from MM theory input validation: {str(e)}" + ) from e + else: + raise AssertionError("No error caught when providing invalid MM theory key.") + + +def test_missing_mm_theory_with_force_field_file(generate_calcjob, generate_inputs): + """Test the co-validation of mm_parameters and force_field_file (mm_parameters).""" + inputs = generate_inputs(structure_fname="butanol.cjson", ff_fname="butanol.ff") + del inputs["mm_parameters"] + try: + generate_calcjob(ChemShellCalculation, inputs) + except ValueError as e: + message = str(e) + assert ( + "A MM theory code must be specified to use molecular mechanics." in message + ), "Wrong error message triggered for missing mm_parameters with force field." + except Exception as e: + raise AssertionError( + f"Wrong error code generated from MM theory input validation: {str(e)}" + ) from e + else: + raise AssertionError("No error caught when providing invalid MM theory key.") + + def test_structure_validation(generate_calcjob, get_test_data_file): """Test error catching for when no input structure is provided.""" # Test if no structure is given @@ -219,3 +256,39 @@ def test_qm_input_validation(generate_calcjob, generate_inputs): ) from e else: raise AssertionError("No error caught during QM parameter validation.") + + +def test_missing_qm_for_qmmm_validation(generate_calcjob, generate_inputs): + """Test check for no QM if QM/MM specified.""" + inputs = generate_inputs(mm={"theory": "DL_POLY"}, ff_fname="butanol.ff") + inputs["qmmm_parameters"] = {"qm_region": [0, 1, 2]} + try: + generate_calcjob(ChemShellCalculation, inputs) + except ValueError as e: + assert "Missing QM parameters for QM/MM calculation" in str(e) + except Exception as e: + raise AssertionError( + f"Wrong error caught during optimisation parameter validation: {str(e)}" + ) from e + else: + raise AssertionError( + "No error caught during optimisation parameter validation." + ) + + +def test_missing_mm_for_qmmm_validation(generate_calcjob, generate_inputs): + """Test check for no QM if QM/MM specified.""" + inputs = generate_inputs() + inputs["qmmm_parameters"] = {"qm_region": [0, 1, 2]} + try: + generate_calcjob(ChemShellCalculation, inputs) + except ValueError as e: + assert "Missing MM parameters for QM/MM calculation" in str(e) + except Exception as e: + raise AssertionError( + f"Wrong error caught during optimisation parameter validation: {str(e)}" + ) from e + else: + raise AssertionError( + "No error caught during optimisation parameter validation." + ) From 89b301dbe12f411c63584e8b3a1e27cf74291549 Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Wed, 4 Mar 2026 08:34:57 +0000 Subject: [PATCH 3/6] Create a geometry optimisation workflow which includes options for vibrational analysis --- src/aiida_chemshell/workflows/optimisation.py | 163 ++++++++++++++++++ tests/test_workflows.py | 35 ++++ 2 files changed, 198 insertions(+) create mode 100644 src/aiida_chemshell/workflows/optimisation.py create mode 100644 tests/test_workflows.py diff --git a/src/aiida_chemshell/workflows/optimisation.py b/src/aiida_chemshell/workflows/optimisation.py new file mode 100644 index 0000000..b8fad99 --- /dev/null +++ b/src/aiida_chemshell/workflows/optimisation.py @@ -0,0 +1,163 @@ +"""Workflows for geometry optimisation based taks.""" + +from aiida.engine import ToContext, WorkChain +from aiida.orm import ( + Bool, + Dict, + Float, + SinglefileData, + Str, +) + +from aiida_chemshell.calculations.base import ChemShellCalculation + + +class GeometryOptimisationWorkChain(WorkChain): + """Geometry optimisation calculation with extended optional calculation options.""" + + @classmethod + def define(cls, spec) -> None: + """Define the AiiDA process specification for the WorkChain.""" + super().define(spec) + + ## Inputs ## + spec.expose_inputs(ChemShellCalculation, exclude=("metadata")) + + spec.input( + "vibrational_analysis", + valid_type=Bool, + required=False, + help="Calculate vibrational modes of resulting structure. (default=True)", + ) + spec.input( + "basis_quality", + valid_type=Str, + validator=cls.validate_basis_quality_input, + required=False, + help="Set basis set quality for QM calculation based on defined options.", + ) + + ## Outputs ## + spec.output( + "final_energy", + valid_type=Float, + required=True, + help="The final energy for the optimised structure.", + ) + spec.output( + "optimised_structure", + valid_type=SinglefileData, + required=True, + help="The final optimised geometry of the given structure.", + ) + spec.output( + "vibrational_analysis", + valid_type=Str, + required=False, + help="The vibrational analysis for the optimised structure.", + ) + + ## Workflow ## + spec.outline(cls.optimise, cls.energy, cls.result) + + return + + def optimise(self): + """Perform the geometry optimisation.""" + inputs = self.exposed_inputs(ChemShellCalculation) + if "qm_parameters" not in self.inputs: + inputs["qm_parameters"] = Dict( + { + "theory": "NWChem", + "method": "dft", + "functional": "B3LYP", + "d3": True, + } + ) + if "basis_quality" in self.inputs: + inputs["qm_parameters"]["basis"] = self.get_basis_set_label( + self.inputs.basis_quality.value + ) + elif "basis_quality" in self.inputs: + inputs["qm_parameters"] = inputs.get("qm_parameters").get_dict() + inputs["qm_parameters"]["basis"] = self.get_basis_set_label( + self.inputs.basis_quality.value + ) + if "force_field_file" in inputs: + if "mm_parameters" not in inputs: + inputs["mm_parameters"] = Dict({"theory": "DL_POLY"}) + if "qmmm_parameters" not in inputs: + inputs["qmmm_parameters"] = Dict({"qm_region": []}) + elif "mm_parameters" in inputs: + return None + if "optimisation_parameters" not in inputs: + inputs["optimisation_parameters"] = Dict({}) + + future = self.submit(ChemShellCalculation, **inputs) + return ToContext(optimise=future) + + def energy(self): + """Perform a single point energy calculation on the optimised structure.""" + if self.inputs.get("vibrational_analysis", False): + inputs = { + "code": self.exposed_inputs(ChemShellCalculation)["code"], + "structure": self.ctx.optimise.outputs.optimised_structure, + "qm_parameters": self.ctx.optimise.inputs.qm_parameters, + } + if "force_field_file" in self.ctx.optimise.inputs: + inputs["force_field_file"] = self.ctx.optimise.inputs.force_field_file + inputs["mm_parameters"] = self.ctx.optimise.inputs.mm_parameters + inputs["qmmm_parameters"] = self.ctx.optimise.inputs.qmmm_parameters + inputs["optimisation_parameters"] = ( + self.ctx.optimise.inputs.optimisation_parameters.get_dict() + ) + inputs["optimisation_parameters"]["thermal"] = True + future = self.submit(ChemShellCalculation, **inputs) + return ToContext(energy=future) + return None + + def result(self): + """Extract the final workflow results.""" + self.out("optimised_structure", self.ctx.optimise.outputs.optimised_structure) + if self.inputs.get("vibrational_analysis", False): + self.out("final_energy", self.ctx.energy.outputs.energy) + self.out( + "vibrational_analysis", self.ctx.energy.outputs.vibrational_analysis + ) + else: + self.out("final_energy", self.ctx.optimise.outputs.energy) + return + + @classmethod + def validate_basis_quality_input(cls, value: Str, _) -> str | None: + """ + Validate the basis set quality key input. + + Parameters + ---------- + value : Str | None + The input key for the basis set quality mapping. + + Returns + ------- + str | None + Returns `None` if input is valid otherwise returns an error message. + """ + if value.value.lower() not in ["fast", "balanced", "quality"]: + return ( + "Invalid basis set quality key, valid keys are: 'fast', 'balanced' or " + "'quality'." + ) + return None + + @classmethod + def get_basis_set_label(cls, key: str) -> str: + """Define a custom dictionary of basis set labels for user convenience.""" + match key.lower(): + case "fast": + return "3-21G" + case "balanced": + return "cc-pvdz" + case "quality": + return "aug-cc-pvtz" + return "" diff --git a/tests/test_workflows.py b/tests/test_workflows.py new file mode 100644 index 0000000..8c88b8d --- /dev/null +++ b/tests/test_workflows.py @@ -0,0 +1,35 @@ +"""Tests for carrying out pre-defined workflows with aiida-chemshell.""" + +from aiida.engine import run_get_node + +from aiida_chemshell.workflows.optimisation import GeometryOptimisationWorkChain + + +def test_geometry_optimisation_workflow(chemsh_code, get_test_data_file): + """Test a geometry optimisation workflow with vibrational analysis.""" + inputs = { + "code": chemsh_code(), + "structure": get_test_data_file(), + "qm_parameters": {"theory": "PySCF", "method": "HF"}, + "basis_quality": "fast", + "vibrational_analysis": True, + } + results, node = run_get_node(GeometryOptimisationWorkChain, **inputs) + + assert node.is_finished_ok, f"WorkChain failed with exit status {node.exit_status}" + + assert len(node.called) > 0, "WorkChain did not launch any subprocesses" + + assert node.called[0].inputs.qm_parameters.get( + "basis", "" + ) == GeometryOptimisationWorkChain.get_basis_set_label("fast") + + assert abs(results.get("final_energy") - -75.585959742867) < 1e-10, ( + "Incorrect final energy for geometry optimisation workflow." + ) + + assert "Temperature: 300.00 Kelvin" in str(results.get("vibrational_analysis")) + assert " Mode Eigenvalue Frequency" in str(results.get("vibrational_analysis")) + assert "total vibrational energy correction" in str( + results.get("vibrational_analysis") + ) From 1c06fbd4602d4a1ee4cb5cd6ba27874f7bcbfd5e Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Wed, 4 Mar 2026 09:37:04 +0000 Subject: [PATCH 4/6] Get vibrational properties as a dictionary output --- src/aiida_chemshell/calculations/base.py | 12 +++++++++ src/aiida_chemshell/parsers/base.py | 25 ++++++++++++++++++- src/aiida_chemshell/workflows/optimisation.py | 9 +++++++ tests/test_workflows.py | 8 ++++++ 4 files changed, 53 insertions(+), 1 deletion(-) diff --git a/src/aiida_chemshell/calculations/base.py b/src/aiida_chemshell/calculations/base.py index cd642a8..616c168 100644 --- a/src/aiida_chemshell/calculations/base.py +++ b/src/aiida_chemshell/calculations/base.py @@ -137,6 +137,18 @@ def define(cls, spec: CalcJobProcessSpec) -> None: required=False, help="The vibrational analysis from ChemShell/DL_FIND", ) + spec.output( + "vibrational_energies", + valid_type=Dict, + required=False, + help="The calculated thermochemical properties of the system.", + ) + spec.output( + "vibrational_modes", + valid_type=ArrayData, + required=False, + help="The calculated vibrational modes of the system.", + ) # Validate inputs namespace existing_validator = spec.inputs.validator diff --git a/src/aiida_chemshell/parsers/base.py b/src/aiida_chemshell/parsers/base.py index 9f687c9..66b78d4 100644 --- a/src/aiida_chemshell/parsers/base.py +++ b/src/aiida_chemshell/parsers/base.py @@ -5,7 +5,7 @@ import numpy from aiida.common import ModificationNotAllowed from aiida.engine import ExitCode -from aiida.orm import ArrayData, Float, SinglefileData, Str +from aiida.orm import ArrayData, Dict, Float, SinglefileData, Str from aiida.parsers.parser import Parser from aiida_chemshell.calculations.base import ChemShellCalculation @@ -110,4 +110,27 @@ def parse_vibrational_analysis(self, stdout: str) -> None: if read: thermo_analysis.append(line) self.out("vibrational_analysis", Str("\n".join(thermo_analysis))) + + energies = {} + for line in stdout.split("\n"): + if read: + line_vals = line.split() + if "Temperature:" in line: + energies[f"Temperature / {line_vals[2]}"] = float(line_vals[1]) + elif "E_electronic" in line: + energies[f"E_electronic correction / {line_vals[7]}"] = float( + line_vals[6] + ) + elif "total ZPE" in line: + energies[f"ZPE / {line_vals[3]}"] = float(line_vals[2]) + elif "total E vib" in line: + energies[f"Enthalpy / {line_vals[4]}"] = float(line_vals[3]) + elif "total S vib" in line: + energies[f"Entropy / {line_vals[4]}"] = float(line_vals[3]) + + if "Thermochemical analysis" in line: + read = True + elif "total S vib" in line: + read = False + self.out("vibrational_energies", Dict(energies)) return diff --git a/src/aiida_chemshell/workflows/optimisation.py b/src/aiida_chemshell/workflows/optimisation.py index b8fad99..b3a936a 100644 --- a/src/aiida_chemshell/workflows/optimisation.py +++ b/src/aiida_chemshell/workflows/optimisation.py @@ -56,6 +56,12 @@ def define(cls, spec) -> None: required=False, help="The vibrational analysis for the optimised structure.", ) + spec.output( + "vibrational_energies", + valid_type=Dict, + required=False, + help="The calculated thermochemical properties of the optimised structure", + ) ## Workflow ## spec.outline(cls.optimise, cls.energy, cls.result) @@ -124,6 +130,9 @@ def result(self): self.out( "vibrational_analysis", self.ctx.energy.outputs.vibrational_analysis ) + self.out( + "vibrational_energies", self.ctx.energy.outputs.vibrational_energies + ) else: self.out("final_energy", self.ctx.optimise.outputs.energy) return diff --git a/tests/test_workflows.py b/tests/test_workflows.py index 8c88b8d..dadd7b1 100644 --- a/tests/test_workflows.py +++ b/tests/test_workflows.py @@ -33,3 +33,11 @@ def test_geometry_optimisation_workflow(chemsh_code, get_test_data_file): assert "total vibrational energy correction" in str( results.get("vibrational_analysis") ) + + assert results.get("vibrational_energies").get("Temperature / Kelvin") == 300.0 + assert results.get("vibrational_energies").get("ZPE / J/mol") == 57173.49993 + assert results.get("vibrational_energies").get("Enthalpy / J/mol") == 3.84524 + assert results.get("vibrational_energies").get("Entropy / J/mol/K") == 0.01430 + + # print(results.get("vibrational_analysis")) + # assert False From 511f18721ac0ef0ab10157f8aa876aef8e8e9b74 Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Wed, 4 Mar 2026 09:57:55 +0000 Subject: [PATCH 5/6] Get vibrational modes as an ArrayData object --- src/aiida_chemshell/calculations/base.py | 8 +----- src/aiida_chemshell/parsers/base.py | 27 ++++++++++--------- src/aiida_chemshell/workflows/optimisation.py | 24 ++++++----------- tests/test_workflows.py | 13 ++++----- 4 files changed, 28 insertions(+), 44 deletions(-) diff --git a/src/aiida_chemshell/calculations/base.py b/src/aiida_chemshell/calculations/base.py index 616c168..461e2d9 100644 --- a/src/aiida_chemshell/calculations/base.py +++ b/src/aiida_chemshell/calculations/base.py @@ -3,7 +3,7 @@ from aiida.common import CalcInfo, CodeInfo from aiida.common.folders import Folder from aiida.engine import CalcJob, CalcJobProcessSpec, PortNamespace -from aiida.orm import ArrayData, Dict, Float, SinglefileData, Str, StructureData +from aiida.orm import ArrayData, Dict, Float, SinglefileData, StructureData from aiida_chemshell.utils import ChemShellMMTheory, ChemShellQMTheory @@ -131,12 +131,6 @@ def define(cls, spec: CalcJobProcessSpec) -> None: "structure is contained within a ChemShell '.pun' file." ), ) - spec.output( - "vibrational_analysis", - valid_type=Str, - required=False, - help="The vibrational analysis from ChemShell/DL_FIND", - ) spec.output( "vibrational_energies", valid_type=Dict, diff --git a/src/aiida_chemshell/parsers/base.py b/src/aiida_chemshell/parsers/base.py index 66b78d4..ed41848 100644 --- a/src/aiida_chemshell/parsers/base.py +++ b/src/aiida_chemshell/parsers/base.py @@ -5,7 +5,7 @@ import numpy from aiida.common import ModificationNotAllowed from aiida.engine import ExitCode -from aiida.orm import ArrayData, Dict, Float, SinglefileData, Str +from aiida.orm import ArrayData, Dict, Float, SinglefileData from aiida.parsers.parser import Parser from aiida_chemshell.calculations.base import ChemShellCalculation @@ -99,19 +99,8 @@ def parse(self, **kwargs): def parse_vibrational_analysis(self, stdout: str) -> None: """Extract the vibrational analysis from ChemShell output log.""" read = False - thermo_analysis = [] - for line in stdout.split("\n"): - if "Thermochemical analysis" in line: - read = True - elif "total S vib" in line: - thermo_analysis.append(line) - read = False - - if read: - thermo_analysis.append(line) - self.out("vibrational_analysis", Str("\n".join(thermo_analysis))) - energies = {} + modes = [] for line in stdout.split("\n"): if read: line_vals = line.split() @@ -127,10 +116,22 @@ def parse_vibrational_analysis(self, stdout: str) -> None: energies[f"Enthalpy / {line_vals[4]}"] = float(line_vals[3]) elif "total S vib" in line: energies[f"Entropy / {line_vals[4]}"] = float(line_vals[3]) + elif "Mode" in line or "total" in line: + pass + else: + # All remaining lines should be part of the modes table + modes.append(numpy.array([float(x) for x in line_vals[2:]])) if "Thermochemical analysis" in line: read = True elif "total S vib" in line: read = False self.out("vibrational_energies", Dict(energies)) + modes = numpy.asarray(modes) + modes_data_node = ArrayData( + label="Vibrational Modes", + description="Calculated vibrational modes for the system.", + ) + modes_data_node.set_array("Modes", modes) + self.out("vibrational_modes", modes_data_node) return diff --git a/src/aiida_chemshell/workflows/optimisation.py b/src/aiida_chemshell/workflows/optimisation.py index b3a936a..2a4e340 100644 --- a/src/aiida_chemshell/workflows/optimisation.py +++ b/src/aiida_chemshell/workflows/optimisation.py @@ -1,13 +1,7 @@ """Workflows for geometry optimisation based taks.""" from aiida.engine import ToContext, WorkChain -from aiida.orm import ( - Bool, - Dict, - Float, - SinglefileData, - Str, -) +from aiida.orm import ArrayData, Bool, Dict, Float, SinglefileData, Str from aiida_chemshell.calculations.base import ChemShellCalculation @@ -50,18 +44,18 @@ def define(cls, spec) -> None: required=True, help="The final optimised geometry of the given structure.", ) - spec.output( - "vibrational_analysis", - valid_type=Str, - required=False, - help="The vibrational analysis for the optimised structure.", - ) spec.output( "vibrational_energies", valid_type=Dict, required=False, help="The calculated thermochemical properties of the optimised structure", ) + spec.output( + "vibrational_modes", + valid_type=ArrayData, + required=False, + help="The calculated vibrational modes for the optimised structure.", + ) ## Workflow ## spec.outline(cls.optimise, cls.energy, cls.result) @@ -127,12 +121,10 @@ def result(self): self.out("optimised_structure", self.ctx.optimise.outputs.optimised_structure) if self.inputs.get("vibrational_analysis", False): self.out("final_energy", self.ctx.energy.outputs.energy) - self.out( - "vibrational_analysis", self.ctx.energy.outputs.vibrational_analysis - ) self.out( "vibrational_energies", self.ctx.energy.outputs.vibrational_energies ) + self.out("vibrational_modes", self.ctx.energy.outputs.vibrational_modes) else: self.out("final_energy", self.ctx.optimise.outputs.energy) return diff --git a/tests/test_workflows.py b/tests/test_workflows.py index dadd7b1..e60299a 100644 --- a/tests/test_workflows.py +++ b/tests/test_workflows.py @@ -28,16 +28,13 @@ def test_geometry_optimisation_workflow(chemsh_code, get_test_data_file): "Incorrect final energy for geometry optimisation workflow." ) - assert "Temperature: 300.00 Kelvin" in str(results.get("vibrational_analysis")) - assert " Mode Eigenvalue Frequency" in str(results.get("vibrational_analysis")) - assert "total vibrational energy correction" in str( - results.get("vibrational_analysis") - ) - assert results.get("vibrational_energies").get("Temperature / Kelvin") == 300.0 assert results.get("vibrational_energies").get("ZPE / J/mol") == 57173.49993 assert results.get("vibrational_energies").get("Enthalpy / J/mol") == 3.84524 assert results.get("vibrational_energies").get("Entropy / J/mol/K") == 0.01430 - # print(results.get("vibrational_analysis")) - # assert False + assert results.get("vibrational_modes").get_shape("Modes") == (3, 5) + + modes = results.get("vibrational_modes").get_array("Modes") + assert (modes[0][0] - 1799.584) < 1e-10, "Incorrect frequency reported for mode 1" + assert (modes[2][2] - 0.0089900284) < 1e-10, "Incorrect ZPE reported for mode 3" From fc34e4e2a576c5be01106b7db94bf9a394003395 Mon Sep 17 00:00:00 2001 From: "Benjamin T. Speake" Date: Wed, 4 Mar 2026 10:09:01 +0000 Subject: [PATCH 6/6] Fix vibrational frequency calculation test --- tests/test_calculations.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/tests/test_calculations.py b/tests/test_calculations.py index 89359eb..57517f6 100644 --- a/tests/test_calculations.py +++ b/tests/test_calculations.py @@ -228,16 +228,12 @@ def test_vibrational_calculation(chemsh_code, get_test_data_file): assert node.is_finished_ok, "CalcJob failed for `test_OptCalculation_dlpoly`" - thermo_analysis = results.get("vibrational_analysis").value + assert results.get("vibrational_modes").get_shape("Modes") == (3, 5) - assert "Thermochemical analysis" in thermo_analysis - assert ( - "0.1252886592 1819.534 2617.900 0.0041452033 0.0000013454 -0.0000014996" - in thermo_analysis - ) - assert "total ZPE 58728.93143 J/mol" in thermo_analysis - assert "total E vib 3.53277 J/mol" in thermo_analysis - assert "total S vib 0.01313 J/mol/K" in thermo_analysis + assert results.get("vibrational_energies").get("Temperature / Kelvin") == 300.0 + assert results.get("vibrational_energies").get("ZPE / J/mol") == 58728.93143 + assert results.get("vibrational_energies").get("Enthalpy / J/mol") == 3.53277 + assert results.get("vibrational_energies").get("Entropy / J/mol/K") == 0.01313 # def test_opt_calculation_qmmm(chemsh_code, get_test_data_file):