diff --git a/HISTORY.rst b/HISTORY.rst
index 56a410c..53ea868 100644
--- a/HISTORY.rst
+++ b/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.
diff --git a/mopac_step/data/properties.csv b/mopac_step/data/properties.csv
index db6e630..41835ff 100644
--- a/mopac_step/data/properties.csv
+++ b/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,
diff --git a/mopac_step/energy.py b/mopac_step/energy.py
index 1274363..8a77e9c 100644
--- a/mopac_step/energy.py
+++ b/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,
diff --git a/mopac_step/forceconstants.py b/mopac_step/forceconstants.py
index 5768d28..5f23ba9 100644
--- a/mopac_step/forceconstants.py
+++ b/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,
diff --git a/mopac_step/metadata.py b/mopac_step/metadata.py
index ee3fc67..5997877 100644
--- a/mopac_step/metadata.py
+++ b/mopac_step/metadata.py
@@ -1246,102 +1246,194 @@
 """
 metadata["results"] = {
     "energy": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "enthalpy of formation",
         "dimensionality": "scalar",
         "type": "float",
         "units": "kcal/mol",
     },
     "gradients": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "gradients on the atoms",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
         "units": "kcal/mol/Å",
     },
+    "force constants": {
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
+        "description": "the Cartesian force constants",
+        "dimensionality": "[natoms*(natoms+1)]",
+        "property": "force constants#MOPAC#{model}",
+        "type": "float",
+        "units": "kcal/mol/Å^2",
+        "format": ".2f",
+    },
     "ERROR_MESSAGE": {
         "description": "An error message",
         "dimensionality": "scalar",
         "type": "string",
     },
     "AO_ATOMINDEX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "atom for AO",
         "dimensionality": ["n_aos"],
         "type": "integer",
     },
     "AO_CHARGES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "electrons in the AOs",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "AO_SPINS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "spins in the AOs",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "AO_ZETA": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "Slater exponent",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "AREA": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "surface area",
         "dimensionality": "scalar",
-        "property": "surface area#MOPAC",
+        "property": "surface area#MOPAC#{model}",
         "type": "float",
         "units": "Å^2",
     },
     "ATOM_CHARGES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "atom charges",
         "dimensionality": ["n_atoms"],
         "type": "float",
     },
     "ATOM_CORE": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "number of valence electrons",
         "dimensionality": ["n_atoms"],
         "type": "integer",
     },
     "ATOM_EL": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "atomic symbol",
         "dimensionality": ["n_atoms"],
         "type": "string",
     },
     "ATOM_PQN": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "principal quantum number",
         "dimensionality": ["n_aos"],
         "type": "integer",
     },
     "ATOM_SYMTYPE": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "atomic orbital angular shape",
         "dimensionality": ["n_aos"],
         "type": "string",
     },
     "ATOM_X": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "x, y, z coordinates",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
         "units": "Å",
     },
     "ATOM_X_FORCE": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "reoriented coordinates",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
         "units": "Å",
     },
     "ATOM_X_OPT": {
-        "calculation": ["optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "optimized x, y, z coordinates",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
@@ -1355,7 +1447,13 @@
         "units": "Å",
     },
     "BOND_ORDERS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "bond order matrix",
         "dimensionality": ["triangular", "n_atoms", "n_atoms"],
         "type": "float",
@@ -1367,14 +1465,26 @@
         "type": "string",
     },
     "CPU_TIME": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "cpu time for calculation",
         "dimensionality": "scalar",
         "type": "float",
         "units": "s",
     },
     "DATE": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "date and time of calculation",
         "dimensionality": "scalar",
         "type": "date_time",
@@ -1391,7 +1501,13 @@
         "units": "g/mL",
     },
     "DIPOLE": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "dipole moment",
         "dimensionality": "scalar",
         "property": "dipole moment#MOPAC#{model}",
@@ -1399,72 +1515,138 @@
         "units": "debye",
     },
     "DIP_VEC": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "dipole vector",
         "dimensionality": [3],
         "type": "float",
         "units": "debye",
     },
     "EIGENVALUES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital energies",
         "dimensionality": ["n_aos"],
         "type": "float",
         "units": "eV",
     },
     "ALPHA_EIGENVALUES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital energies",
         "dimensionality": ["n_aos"],
         "type": "float",
         "units": "eV",
     },
     "BETA_EIGENVALUES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital energies",
         "dimensionality": ["n_aos"],
         "type": "float",
         "units": "eV",
     },
     "LMO_ENERGY_LEVELS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "Localized MO energies",
         "dimensionality": ["n_aos"],
         "type": "float",
         "units": "eV",
     },
     "EIGENVECTORS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital coefficients",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "ALPHA_EIGENVECTORS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital coefficients",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "BETA_EIGENVECTORS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital coefficients",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "LMO_VECTORS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "localized molecular orbital coefficients",
         "dimensionality": ["n_aos"],
         "type": "float",
     },
     "EMPIRICAL_FORMULA": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "empirical formula",
         "dimensionality": "scalar",
         "type": "string",
     },
     "ENERGY_ELECTRONIC": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "electronic energy",
         "dimensionality": "scalar",
         "property": "electronic energy#MOPAC#{model}",
@@ -1472,7 +1654,13 @@
         "units": "eV",
     },
     "ENERGY_NUCLEAR": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "nuclear repulsion energy",
         "dimensionality": "scalar",
         "property": "nuclear repulsion energy",
@@ -1480,7 +1668,13 @@
         "units": "eV",
     },
     "DIEL_ENER": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "the dielectric energy from COSMO",
         "dimensionality": "scalar",
         "property": "dielectric energy#MOPAC#{model}",
@@ -1488,7 +1682,7 @@
         "units": "eV",
     },
     "ENTHALPY_TOT": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "enthalpy",
         "dimensionality": [
             "n_temps",
@@ -1497,7 +1691,7 @@
         "units": "cal/mol",
     },
     "ENTROPY_TOT": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "entropy",
         "dimensionality": [
             "n_temps",
@@ -1506,14 +1700,25 @@
         "units": "cal/K/mol",
     },
     "GRADIENTS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "gradients on the atoms",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
         "units": "kcal/mol/Å",
     },
     "GRADIENT_NORM": {
-        "calculation": ["optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "final norm of the gradient",
         "dimensionality": "scalar",
         "type": "float",
@@ -1546,7 +1751,7 @@
         "units": "GPa",
     },
     "HEAT_CAPACITY_TOT": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "heat capacity",
         "dimensionality": [
             "n_temps",
@@ -1555,7 +1760,13 @@
         "units": "cal/K/mol",
     },
     "HEAT_OF_FORMATION": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "enthalpy of formation",
         "dimensionality": "scalar",
         "property": "enthalpy of formation#MOPAC#{model}",
@@ -1570,21 +1781,21 @@
         "units": "kcal/mol",
     },
     "HESSIAN_MATRIX": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "Hessian matrix",
         "dimensionality": ["n_dof", "n_dof"],
         "type": "float",
         "units": "mdyne/Å/Da",
     },
     "H_O_F(T)": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "enthalpy of formation vs T",
         "dimensionality": ["n_temps"],
         "type": "float",
         "units": "kcal/mol",
     },
     "INT_FORCE_CONSTS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "force constants for internals",
         "dimensionality": ["n_dofs"],
         "type": "float",
@@ -1601,7 +1812,13 @@
         "units": "Å",
     },
     "IONIZATION_POTENTIAL": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "ionization energy (IE)",
         "dimensionality": "scalar",
         "property": "ionization energy#MOPAC#{model}",
@@ -1609,7 +1826,7 @@
         "units": "eV",
     },
     "ISOTOPIC_MASSES": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "isotopic masses",
         "dimensionality": ["n_atoms"],
         "type": "float",
@@ -1622,123 +1839,219 @@
         "type": "string",
     },
     "M.O.SYMMETRY_LABELS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "molecular orbital symmetries",
         "dimensionality": ["n_mos"],
         "type": "string",
     },
     "ALPHA_M.O.SYMMETRY_LABELS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "molecular orbital symmetries",
         "dimensionality": ["n_mos"],
         "type": "string",
     },
     "BETA_M.O.SYMMETRY_LABELS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "molecular orbital symmetries",
         "dimensionality": ["n_mos"],
         "type": "string",
     },
     "METHOD": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "hamiltonian",
         "dimensionality": "scalar",
         "type": "string",
     },
     "MOLECULAR_ORBITAL_OCCUPANCIES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital occupancies",
         "dimensionality": ["n_mos"],
         "type": "float",
     },
     "ALPHA_MOLECULAR_ORBITAL_OCCUPANCIES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital occupancies",
         "dimensionality": ["n_mos"],
         "type": "float",
     },
     "BETA_MOLECULAR_ORBITAL_OCCUPANCIES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "orbital occupancies",
         "dimensionality": ["n_mos"],
         "type": "float",
     },
     "MOLECULAR_WEIGHT": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "molecular weight",
         "dimensionality": "scalar",
         "type": "float",
         "units": "Da",
     },
     "MOPAC_VERSION": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "MOPAC version",
         "dimensionality": "scalar",
         "type": "string",
     },
     "NORMAL_MODES": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "normal modes",
         "dimensionality": ["n_dof", "n_dof"],
         "type": "float",
     },
     "NORMAL_MODE_SYMMETRY_LABELS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "symmetry labels of normal modes",
         "dimensionality": ["n_dof"],
         "type": "string",
     },
     "NUMBER_SCF_CYCLES": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "number of scf's",
         "dimensionality": "scalar",
         "type": "integer",
     },
     "NUM_ALPHA_ELECTRONS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "number of spin-up electrons",
         "dimensionality": "scalar",
         "type": "integer",
     },
     "NUM_BETA_ELECTRONS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "number of spin-down electrons",
         "dimensionality": "scalar",
         "type": "integer",
     },
     "NUM_ELECTRONS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "number of electrons",
         "dimensionality": "scalar",
         "type": "integer",
     },
     "ORIENTATION_ATOM_X": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "",
         "dimensionality": [3, "n_atoms"],
         "type": "float",
         "units": "Å",
     },
     "OVERLAP_MATRIX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "the AO overlap matrix",
         "dimensionality": ["triangular", "n_aos", "n_aos"],
         "type": "float",
     },
     "DENSITY_MATRIX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "the density matrix",
         "dimensionality": ["triangular", "n_aos", "n_aos"],
         "type": "float",
     },
     "POINT_GROUP": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "the molecular symmetry",
         "dimensionality": "scalar",
         "type": "string",
     },
     "PRI_MOM_OF_I": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "primary moments of inertia",
         "dimensionality": [3],
         "type": "float",
@@ -1755,26 +2068,44 @@
         "units": "GPa",
     },
     "ROTAT_CONSTS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "rotational constants",
         "dimensionality": [3],
         "type": "float",
         "units": "1/cm",
     },
     "SET_OF_MOS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "set of MOs",
         "dimensionality": [2],
         "type": "string",
     },
     "SET_OF_ALPHA_MOS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "set of MOs",
         "dimensionality": [2],
         "type": "string",
     },
     "SET_OF_BETA_MOS": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "set of MOs",
         "dimensionality": [2],
         "type": "string",
@@ -1854,13 +2185,19 @@
         "type": "string",
     },
     "SPIN_COMPONENT": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "spin component",
         "dimensionality": [3],
         "type": "float",
     },
     "THERMODYNAMIC_PROPERTIES_TEMPS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "temperature for thermodynamic properties",
         "dimensionality": ["n_temps"],
         "type": "float",
@@ -1873,28 +2210,52 @@
         "type": "string",
     },
     "TOTAL_DENSITY_MATRIX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "density matrix",
         "dimensionality": ["n_mos", "n_mos"],
         "shape": "triangular",
         "type": "float",
     },
     "ALPHA_DENSITY_MATRIX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "density matrix",
         "dimensionality": ["n_mos", "n_mos"],
         "shape": "triangular",
         "type": "float",
     },
     "BETA_DENSITY_MATRIX": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "density matrix",
         "dimensionality": ["n_mos", "n_mos"],
         "shape": "triangular",
         "type": "float",
     },
     "TOTAL_ENERGY": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "total energy",
         "dimensionality": "scalar",
         "property": "total energy#MOPAC#{model}",
@@ -1902,7 +2263,13 @@
         "units": "eV",
     },
     "HOMO Energy": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "energy of the HOMO",
         "dimensionality": "scalar",
         "property": "HOMO energy#MOPAC#{model}",
@@ -1910,7 +2277,13 @@
         "units": "eV",
     },
     "LUMO Energy": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "energy of the LUMO",
         "dimensionality": "scalar",
         "property": "LUMO energy#MOPAC#{model}",
@@ -1918,7 +2291,13 @@
         "units": "eV",
     },
     "HOMO-LUMO Gap": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "The energy of the HOMO-LUMO gap",
         "dimensionality": "scalar",
         "property": "band gap#MOPAC#{model}",
@@ -1926,7 +2305,13 @@
         "units": "eV",
     },
     "TOTAL_SPIN": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "spin",
         "dimensionality": "scalar",
         "property": "S^2#MOPAC#{model}",
@@ -1953,49 +2338,55 @@
         "units": "Å",
     },
     "VIB._EFF_MASS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "vibrational effective masses",
         "dimensionality": ["n_dof"],
         "type": "float",
         "units": "Da",
     },
     "VIB._FREQ": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "vibrational frequencies",
         "dimensionality": ["n_dof"],
         "type": "float",
         "units": "1/cm",
     },
     "VIB._RED_MASS": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "vibrational reduced masses",
         "dimensionality": ["n_dof"],
         "type": "float",
         "units": "Da",
     },
     "VIB._TRAVEL": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "vibrational lenght of travel",
         "dimensionality": ["n_dof"],
         "type": "float",
         "units": "Å",
     },
     "VIB._T_DIP": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "vibrational transition dipole",
         "dimensionality": ["n_dof"],
         "type": "float",
     },
     "VOLUME": {
-        "calculation": ["energy", "optimization", "thermodynamics", "vibrations"],
+        "calculation": [
+            "energy",
+            "optimization",
+            "thermodynamics",
+            "vibrations",
+            "force constants",
+        ],
         "description": "volume",
         "dimensionality": "scalar",
-        "property": "volume#MOPAC",
+        "property": "volume#MOPAC#{model}",
         "type": "float",
         "units": "Å^3",
     },
     "ZERO_POINT_ENERGY": {
-        "calculation": ["thermodynamics", "vibrations"],
+        "calculation": ["thermodynamics", "vibrations", "force constants"],
         "description": "zero point energy",
         "dimensionality": "scalar",
         "property": "zero point energy#MOPAC#{model}",