Merge pull request #34 from molssi-seamm/dev

Enhanced results and fixed problem with psi4.ini

HISTORY.rst

@@ -1,7 +1,12 @@
 =======
 History
 =======
-2024.5.23.3 -- Added auxialliary codes that help Psi4
+2024.7.30 -- Enhanced results and fixed problem with psi4.ini
+   * Added to the results by using cclib to parse the output.
+   * If ~/SEAMM/psi4.ini did not exist the version that was automatically created was
+     not complete, causing calculations to fail.
+
+2024.5.23.3 -- Added auxiliary codes that help Psi4
    * Psi4 uses codes like DFTD4 and DFTD4 for parts of e.g. dispersion
      calculations. This release adds them to the Conda install SEAMM uses for Psi4.
 

devtools/conda-envs/test_env.yaml

@@ -9,6 +9,7 @@ dependencies:
   - pip
 
   # Of the code
+  - cclib
   - numpy
   - optking
   - scipy

psi4_step/accelerated_optimization.py

@@ -1123,7 +1123,7 @@ def _bfgs(
         jac=None,
         callback=None,
         gtol=1e-5,
-        norm=np.Inf,
+        norm=np.inf,
         eps=_epsilon,
         maxiter=None,
         disp=False,

psi4_step/energy.py

@@ -5,8 +5,11 @@
 import json
 import logging
 from pathlib import Path
+import pprint
 
 from openbabel import openbabel
+import cclib
+import numpy as np
 
 import psi4_step
 import seamm
@@ -359,12 +362,25 @@ def analyze(self, indent="", data={}, out=[]):
         """Parse the output and generating the text output and store the
         data in variables for other stages to access
         """
+        directory = Path(self.directory)
+
+        # Use cclib to get results
+        path = directory / "output.dat"
+        if path.exists():
+            data = vars(cclib.io.ccread(path))
+            pprint.pformat(data)
+            data = self.process_data(data)
+            pprint.pformat(data)
+        else:
+            data = {}
+
         # Read in the results from json
         directory = Path(self.directory)
         json_file = directory / "properties.json"
         if json_file.exists():
             with json_file.open() as fd:
-                data = json.load(fd)
+                tmp = json.load(fd)
+            data.update(**tmp)
 
             # Put any requested results into variables or tables
             self.store_results(
@@ -498,3 +514,102 @@ def plot_input(self):
         )
 
         return "\n".join(lines)
+
+    def process_data(self, data):
+        """Massage the cclib data to a more easily used form."""
+        logger.debug(pprint.pformat(data))
+        # Convert numpy arrays to Python lists
+        new = {}
+        for key, value in data.items():
+            if isinstance(value, np.ndarray):
+                new[key] = value.tolist()
+            elif isinstance(value, list):
+                if len(value) > 0 and isinstance(value[0], np.ndarray):
+                    new[key] = [i.tolist() for i in value]
+                else:
+                    new[key] = value
+            elif isinstance(value, dict):
+                for k, v in value.items():
+                    newkey = f"{key}/{k}"
+                    if isinstance(v, np.ndarray):
+                        new[newkey] = v.tolist()
+                    else:
+                        new[newkey] = v
+            else:
+                new[key] = value
+
+        for key in ("metadata/cpu_time", "metadata/wall_time"):
+            if key in new:
+                time = new[key][0]
+                for tmp in new[key][1:]:
+                    time += tmp
+                new[key] = str(time).lstrip("0:")
+                if "." in new[key]:
+                    new[key] = new[key].rstrip("0")
+
+        # Pull out the HOMO and LUMO energies as scalars
+        if "homos" in new and "moenergies" in new:
+            homos = new["homos"]
+            if len(homos) == 2:
+                for i, letter in enumerate(["α", "β"]):
+                    Es = new["moenergies"][i]
+                    homo = homos[i]
+                    new[f"N({letter}-homo)"] = homo + 1
+                    new[f"E({letter}-homo)"] = Es[homo]
+                    if homo > 0:
+                        new[f"E({letter}-homo-1)"] = Es[homo - 1]
+                    if homo + 1 < len(Es):
+                        new[f"E({letter}-lumo)"] = Es[homo + 1]
+                        new[f"E({letter}-gap)"] = Es[homo + 1] - Es[homo]
+                    if homo + 2 < len(Es):
+                        new[f"E({letter}-lumo+1)"] = Es[homo + 2]
+                    if "mosyms" in new:
+                        syms = new["mosyms"][i]
+                        new[f"Sym({letter}-homo)"] = syms[homo]
+                        if homo > 0:
+                            new[f"Sym({letter}-homo-1)"] = syms[homo - 1]
+                        if homo + 1 < len(syms):
+                            new[f"Sym({letter}-lumo)"] = syms[homo + 1]
+                        if homo + 2 < len(syms):
+                            new[f"Sym({letter}-lumo+1)"] = syms[homo + 2]
+            else:
+                Es = new["moenergies"][0]
+                homo = homos[0]
+                new["N(homo)"] = homo + 1
+                new["E(homo)"] = Es[homo]
+                if homo > 0:
+                    new["E(homo-1)"] = Es[homo - 1]
+                if homo + 1 < len(Es):
+                    new["E(lumo)"] = Es[homo + 1]
+                    new["E(gap)"] = Es[homo + 1] - Es[homo]
+                if homo + 2 < len(Es):
+                    new["E(lumo+1)"] = Es[homo + 2]
+                if "mosyms" in new:
+                    syms = new["mosyms"][0]
+                    new["Sym(homo)"] = syms[homo]
+                    if homo > 0:
+                        new["Sym(homo-1)"] = syms[homo - 1]
+                    if homo + 1 < len(syms):
+                        new["Sym(lumo)"] = syms[homo + 1]
+                    if homo + 2 < len(syms):
+                        new["Sym(lumo+1)"] = syms[homo + 2]
+
+        # moments
+        if "moments" in new:
+            moments = new["moments"]
+            new["multipole_reference"] = moments[0]
+            new["dipole_moment"] = moments[1]
+            new["dipole_moment_magnitude"] = np.linalg.norm(moments[1])
+            if len(moments) > 2:
+                new["quadrupole_moment"] = moments[2]
+            if len(moments) > 3:
+                new["octapole_moment"] = moments[3]
+            if len(moments) > 4:
+                new["hexadecapole_moment"] = moments[4]
+            del new["moments"]
+
+        for key in ("metadata/symmetry_detected", "metadata/symmetry_used"):
+            if key in new:
+                new[key] = new[key].capitalize()
+
+        return new

psi4_step/psi4.py

@@ -7,13 +7,15 @@
 import importlib
 import json
 import logging
+import os
 from pathlib import Path
 import pprint
 import shutil
 
 import psi4_step
 import seamm
 import seamm_exec
+from seamm_util import Configuration
 import seamm_util.printing as printing
 from seamm_util.printing import FormattedText as __
 
@@ -462,14 +464,19 @@ def run(self):
             ini_dir = Path(seamm_options["root"]).expanduser()
             path = ini_dir / "psi4.ini"
 
-            if path.exists():
-                full_config.read(ini_dir / "psi4.ini")
-
-            # If the section we need doesn't exists, get the default
-            if not path.exists() or executor_type not in full_config:
+            # If the config file doesn't exists, get the default
+            if not path.exists():
                 resources = importlib.resources.files("psi4_step") / "data"
                 ini_text = (resources / "psi4.ini").read_text()
-                full_config.read_string(ini_text)
+                txt_config = Configuration(path)
+                txt_config.from_string(ini_text)
+
+                # Work out the conda info needed
+                txt_config.set_value("local", "conda", os.environ["CONDA_EXE"])
+                txt_config.set_value("local", "conda-environment", "seamm-psi4")
+                txt_config.save()
+
+            full_config.read(ini_dir / "psi4.ini")
 
             # Getting desperate! Look for an executable in the path
             if executor_type not in full_config:
@@ -481,10 +488,12 @@ def run(self):
                         "in the path!"
                     )
                 else:
-                    full_config[executor_type] = {
-                        "installation": "local",
-                        "code": f"{path} -n {{NTASKS}}",
-                    }
+                    txt_config = Configuration(path)
+                    txt_config.add_section(executor_type)
+                    txt_config.set_value(executor_type, "installation", "local")
+                    txt_config.set_value(executor_type, "code", f"{path} -n {{NTASKS}}")
+                    txt_config.save()
+                    full_config.read(ini_dir / "psi4.ini")
 
             config = dict(full_config.items(executor_type))
             # Use the matching version of the seamm-psi4 image by default.

psi4_step/psi4_metadata.py

@@ -2228,7 +2228,7 @@
     },
     "SCF ITERATIONS": {
         "calculation": ["energy", "optimization", "thermochemistry", "vibrations"],
-        "description": "number of itereations in the SCF",
+        "description": "number of iterations in the SCF",
         "dimensionality": "scalar",
         "methods": ["ccsd", "ccsd(t)", "dft", "hf", "lccd", "mp2", "mp3", "mp4"],
         "type": "integer",
@@ -2776,4 +2776,118 @@
         "type": "integer",
         "units": "",
     },
+    # From cclib
+    "moenergies": {
+        "calculation": ["energy", "optimization"],
+        "description": "The energy of the molecular orbitals (MOs)",
+        "dimensionality": "[[nmo]*nspin]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "eV",
+        "format": ".3f",
+    },
+    "moments": {
+        "calculation": ["energy", "optimization"],
+        "description": "",
+        "dimensionality": "scalar",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "",
+        "format": ".6f",
+    },
+    "multipole_reference": {
+        "calculation": ["energy", "optimization"],
+        "description": "The reference point for the multipole expansions",
+        "dimensionality": "[3]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "Å",
+        "format": ".6f",
+    },
+    "dipole_moment": {
+        "calculation": ["energy", "optimization"],
+        "description": "The dipole moment",
+        "dimensionality": "[3]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "debye",
+        "format": ".6f",
+    },
+    "dipole_moment_magnitude": {
+        "calculation": ["energy", "optimization"],
+        "description": "The magnitude of the dipole moment",
+        "dimensionality": "scalar",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "debye",
+        "format": ".6f",
+    },
+    "quadrupole_moment": {
+        "calculation": ["energy", "optimization"],
+        "description": "The quadrupole moment",
+        "dimensionality": "[6]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "debye*Å",
+        "format": ".6f",
+    },
+    "octapole_moment": {
+        "calculation": ["energy", "optimization"],
+        "description": "The octapole moment",
+        "dimensionality": "[10]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "debye*Å^2",
+        "format": ".6f",
+    },
+    "hexadecapole_moment": {
+        "calculation": ["energy", "optimization"],
+        "description": "The hexadecapole moment",
+        "dimensionality": "[15]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "debye*Å^3",
+        "format": ".6f",
+    },
+    "mosysm": {
+        "calculation": ["energy", "optimization"],
+        "description": "The symmetry of the MOs",
+        "dimensionality": "[nmo]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "string",
+        "units": "",
+        "format": "",
+    },
+    "nbasis": {
+        "description": "Number of basis functions",
+        "dimensionality": "scalar",
+        "type": "integer",
+        "units": "",
+        "format": "",
+    },
+    "nmo": {
+        "description": "Number of molecular orbitals",
+        "dimensionality": "scalar",
+        "type": "integer",
+        "units": "",
+        "format": "",
+    },
+    "rotconstants": {
+        "calculation": ["energy", "optimization"],
+        "description": "Rotational constants",
+        "dimensionality": "[[3]*nsteps]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "GHz",
+        "format": ".3f",
+    },
+    "scfenergies": {
+        "calculation": ["energy", "optimization"],
+        "description": "The energies per step",
+        "dimensionality": "[nsteps]",
+        "methods": ["DFT", "HF", "MP2", "MP3", "MP4"],
+        "type": "float",
+        "units": "eV",
+        "format": ".3f",
+    },
 }

requirements.txt

@@ -1,3 +1,4 @@
+cclib
 numpy
 # optking
 scipy