Merge pull request #139 from molssi-seamm/dev

Add the force constants as a property.

HISTORY.rst

@@ -1,6 +1,10 @@
 =======
 History
 =======
+2024.12.9 -- Add the force constants as a property.
+   * Add the force constants as a property of the configuration when running
+     thermodynamics, IR, or force constants calculations.
+
 2024.10.15 -- Bugfix: error if used in a loop and previous directories deleted.
    * The code crashed if called with a loop in the flowchart, and the last directory of
      a previous loop iteration was deleted before running the next iteration.

mopac_step/data/properties.csv

@@ -10,7 +10,8 @@ electronic energy#MOPAC#{model},float,eV,"The {model} electronic energy.",
 enthalpy of formation#MOPAC#{model},float,kJ/mol,"The {model} enthalpy of formation, DHf.",https://en.wikipedia.org/wiki/Standard_enthalpy_of_formation
 entropy#MOPAC#{model},float,J/mol/K,"The {model} standard entropy.",https://en.wikipedia.org/wiki/Entropy
 ionization energy#MOPAC#{model},float,eV,"The {model} ionization energy.",https://en.wikipedia.org/wiki/Ionization_energy
-surface area#MOPAC,float,Å^2,"The molecular surface area.",
+surface area#MOPAC#{model},float,Å^2,"The molecular surface area.",
 total energy#MOPAC#{model},float,eV,"The {model} total energy (electronic + nuclear).",
-volume#MOPAC,float,Å^3,"The molecular volume.",
+volume#MOPAC#{model},float,Å^3,"The molecular volume.",
 zero point energy#MOPAC#{model},float,kJ/mol,"The {model} zero-point vibrational energy.",
+force constants#MOPAC#{model},json,kcal/mol/Å^2,The {model} force constants,

mopac_step/energy.py

@@ -5,6 +5,7 @@
 import copy
 import csv
 import logging
+from math import sqrt
 from pathlib import Path
 import pprint  # noqa: F401
 import textwrap
@@ -1153,6 +1154,25 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
                         tmp -= delta
                         data["gradients"] = tmp.tolist()
 
+        # Handle the force constant matrix (Hessian) if it exists
+        if "HESSIAN_MATRIX" in data:
+            # It is mass weighted so we need to remove the weighting
+            if "ISOTOPIC_MASSES" not in data:
+                raise RuntimeError("Found no atomic masses")
+            # Expand the mass array for x, y, z
+            mass = [v for v in data["ISOTOPIC_MASSES"] for j in range(3)]
+
+            # Get the atom part of the force constant matrix.
+            hessian = data["HESSIAN_MATRIX"]
+            ij = 0
+            factor = Q_(1.0, "mdyne/Å").m_as("kcal/mol/Å^2")
+            tmp = []
+            for i, mass_i in enumerate(mass):
+                for mass_j in mass[0 : i + 1]:
+                    tmp.append(factor * hessian[ij] * sqrt(mass_i * mass_j))
+                    ij += 1
+            data["force constants"] = tmp
+
         self.store_results(
             configuration=configuration,
             data=data,

mopac_step/forceconstants.py

@@ -26,7 +26,7 @@ def __init__(self, flowchart=None, title="Force Constants", extension=None):
 
         super().__init__(flowchart=flowchart, title=title, extension=extension)
 
-        self._calculation = "vibrations"
+        self._calculation = "force constants"
         self._model = None
         self._metadata = mopac_step.metadata
         self.parameters = mopac_step.ForceconstantsParameters()
@@ -221,7 +221,7 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
             The forceconstant calculation
 
         Since the forceconstant matrix is symmetric (we hope!) we will work with the
-        lower triangle as a lineary array.
+        lower triangle as a linear array.
         """
         P = self.parameters.current_values_to_dict(
             context=seamm.flowchart_variables._data
@@ -466,6 +466,11 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
                     ij += 1
                 fd.write("\n")
 
+        # Save the force constant matrix (Hessian) in the data sections
+        data = data_sections[0]
+        fac = Q_(1.0, P["atom_units"]).m_as("kJ/mol/Å^2")
+        data["force constants"] = [v * fac for v in result]
+
         # Let the energy module do its thing
         super().analyze(
             indent=indent,