Merge pull request #22 from molssi-seamm/dev

Added thermochemistry and orbital/density plots.
molssi-seamm · Jun 5, 2023 · 80c6e4b · 80c6e4b
2 parents 5d645c1 + ab5634a
commit 80c6e4b
Show file tree

Hide file tree

Showing 18 changed files with 1,589 additions and 181 deletions.
diff --git a/HISTORY.rst b/HISTORY.rst
@@ -2,6 +2,11 @@
 History
 =======
 
+2023.6.4 -- Enhancements
+   * Added thermochemistry substep to compute the vibrational and other corrections for
+     thermochemistry.
+   * Added ability to create .cube files for plotting the density and orbitals.
+
 2023.2.21 -- Added control over SCF convergence
    * Able to set convergence criteria.
    * Options for damping. level shifting and second-order SCF

diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml
@@ -14,6 +14,7 @@ dependencies:
   - psutil
   - scipy
   - seamm
+  - tabulate
 
   # Testing
   - black

diff --git a/psi4_step/__init__.py b/psi4_step/__init__.py
@@ -8,43 +8,48 @@
 # Bring up the classes so that they appear to be directly in
 # the psi4_step package.
 
-from psi4_step.psi4_metadata import methods  # noqa: F401
-from psi4_step.psi4_metadata import dft_functionals  # noqa: F401
-from psi4_step.psi4_metadata import optimization_methods  # noqa: F401
-from psi4_step.psi4_metadata import optimization_convergence  # noqa: F401
-from psi4_step.psi4_metadata import metadata  # noqa: F401
-
-from psi4_step.psi4 import Psi4  # noqa: F401
-from psi4_step.psi4_parameters import Psi4Parameters  # noqa: F401
-from psi4_step.psi4_step import Psi4Step  # noqa: F401
-from psi4_step.tk_psi4 import TkPsi4  # noqa: F401
-
-from psi4_step.initialization import Initialization  # noqa: F401
-from psi4_step.initialization_parameters import InitializationParameters  # noqa: F401
-from psi4_step.initialization_step import InitializationStep  # noqa: F401
-from psi4_step.tk_initialization import TkInitialization  # noqa: F401
-
-from psi4_step.energy import Energy  # noqa: F401
-from psi4_step.energy_parameters import EnergyParameters  # noqa: F401
-from psi4_step.energy_step import EnergyStep  # noqa: F401
-from psi4_step.tk_energy import TkEnergy  # noqa: F401
-
-from psi4_step.optimization import Optimization  # noqa: F401
-from psi4_step.optimization_parameters import OptimizationParameters  # noqa: F401
-from psi4_step.optimization_step import OptimizationStep  # noqa: F401
-from psi4_step.tk_optimization import TkOptimization  # noqa: F401
-
-# from psi4_step.accelerated_optimization import AcceleratedOptimization  # noqa: F401
-# from psi4_step.accelerated_optimization_parameters import (  # noqa: F401
+from .psi4_metadata import methods  # noqa: F401
+from .psi4_metadata import dft_functionals  # noqa: F401
+from .psi4_metadata import optimization_methods  # noqa: F401
+from .psi4_metadata import optimization_convergence  # noqa: F401
+from .psi4_metadata import metadata  # noqa: F401
+
+from .psi4 import Psi4  # noqa: F401
+from .psi4_parameters import Psi4Parameters  # noqa: F401
+from .psi4_step import Psi4Step  # noqa: F401
+from .tk_psi4 import TkPsi4  # noqa: F401
+
+from .initialization import Initialization  # noqa: F401
+from .initialization_parameters import InitializationParameters  # noqa: F401
+from .initialization_step import InitializationStep  # noqa: F401
+from .tk_initialization import TkInitialization  # noqa: F401
+
+from .energy import Energy  # noqa: F401
+from .energy_parameters import EnergyParameters  # noqa: F401
+from .energy_step import EnergyStep  # noqa: F401
+from .tk_energy import TkEnergy  # noqa: F401
+
+from .optimization import Optimization  # noqa: F401
+from .optimization_parameters import OptimizationParameters  # noqa: F401
+from .optimization_step import OptimizationStep  # noqa: F401
+from .tk_optimization import TkOptimization  # noqa: F401
+
+# from .accelerated_optimization import AcceleratedOptimization  # noqa: F401
+# from .accelerated_optimization_parameters import (  # noqa: F401
 #     AcceleratedOptimizationParameters,
 # )
-# from psi4_step.accelerated_optimization_step import (  # noqa: F401
+# from .accelerated_optimization_step import (  # noqa: F401
 #     AcceleratedOptimizationStep,
 # )
-# from psi4_step.tk_accelerated_optimization import (  # noqa: F401
+# from .tk_accelerated_optimization import (  # noqa: F401
 #     TkAcceleratedOptimization,
 # )
 
+from .thermochemistry_step import ThermochemistryStep  # noqa: F401
+from .thermochemistry import Thermochemistry  # noqa: F401
+from .thermochemistry_parameters import ThermochemistryParameters  # noqa: F401
+from .tk_thermochemistry import TkThermochemistry  # noqa: F401
+
 # Handle versioneer
 from ._version import get_versions
 

diff --git a/psi4_step/energy.py b/psi4_step/energy.py
@@ -19,7 +19,7 @@
 
 
 class Energy(seamm.Node):
-    def __init__(self, flowchart=None, title="Energy", extension=None):
+    def __init__(self, flowchart=None, title="Energy", extension=None, logger=logger):
         """Initialize the node"""
 
         logger.debug("Creating Energy {}".format(self))
@@ -92,6 +92,14 @@ def description_text(self, P=None, calculation_type="Single-point energy"):
             text += " The spin will not be restricted and may not be a "
             text += "proper eigenstate."
 
+        # Plotting
+        if P["density"]:
+            if P["orbitals"]:
+                text += "\nThe alpha and beta electron, total, and spin densities, "
+                text += f"and orbitals {P['selected orbitals']} will be plotted."
+        elif P["orbitals"]:
+            text += f"\nThe orbitals {P['selected orbitals']} will be plotted."
+
         return self.header + "\n" + __(text, **P, indent=4 * " ").__str__()
 
     def get_input(self, calculation_type="energy", restart=None):
@@ -277,7 +285,9 @@ def get_input(self, calculation_type="energy", restart=None):
             lines.append(
                 "    json.dump(fix_multipoles(variables), fd, sort_keys=True, indent=3)"
             )
-        lines.append("")
+
+        # Orbital plots
+        lines.append(self.plot_input())
 
         return "\n".join(lines)
 
@@ -303,7 +313,7 @@ def analyze(self, indent="", data={}, out=[]):
                 create_tables=self.parameters["create tables"].get(),
             )
 
-            text = "The calculated energy is {Eelec:.6f} Ha."
+            text = "The calculated energy is {Eelec:.6f} E_h."
         else:
             data = {}
             tmp = str(json_file)
@@ -315,3 +325,106 @@ def analyze(self, indent="", data={}, out=[]):
             raise RuntimeError(text)
 
         printer.normal(__(text, **data, indent=self.indent + 4 * " "))
+
+    def plot_input(self):
+        """Generate the input for plotting to cube files."""
+        _, configuration = self.get_system_configuration(None)
+
+        P = self.parameters.current_values_to_dict(
+            context=seamm.flowchart_variables._data
+        )
+
+        tasks = []
+        if P["density"]:
+            tasks.append("density")
+
+        orbitals = []
+        if P["orbitals"]:
+            tasks.append("orbitals")
+            # and work out the orbitals
+            txt = P["selected orbitals"]
+            if txt == "all":
+                pass
+            else:
+                # Which is the HOMO orbital?
+                # This will not work with ECPs.
+                n_electrons = (
+                    sum(configuration.atoms.atomic_numbers) - configuration.charge
+                )
+                multiplicity = configuration.spin_multiplicity
+                homo = (n_electrons - (multiplicity - 1)) // 2 + (multiplicity - 1)
+
+                for chunk in txt.split(","):
+                    chunk = chunk.strip()
+                    if ":" in chunk or ".." in chunk:
+                        if ":" in chunk:
+                            first, last = chunk.split(":")
+                        elif ".." in chunk:
+                            first, last = chunk.split("..")
+                        first = first.strip().upper()
+                        last = last.strip().upper()
+
+                        if first == "HOMO":
+                            first = homo
+                        elif first == "LUMO":
+                            first = homo + 1
+                        else:
+                            first = int(first.removeprefix("HOMO").removeprefix("LUMO"))
+                            if first < 0:
+                                first = homo + first
+                            else:
+                                first = homo + 1 + first
+
+                        if last == "HOMO":
+                            last = homo
+                        elif last == "LUMO":
+                            last = homo + 1
+                        else:
+                            last = int(last.removeprefix("HOMO").removeprefix("LUMO"))
+                            if last < 0:
+                                last = homo + last
+                            else:
+                                last = homo + 1 + last
+
+                        orbitals.extend(range(first, last + 1))
+                    else:
+                        first = chunk.strip().upper()
+
+                        if first == "HOMO":
+                            first = homo
+                        elif first == "LUMO":
+                            first = homo + 1
+                        else:
+                            first = int(first.removeprefix("HOMO").removeprefix("LUMO"))
+                            if first < 0:
+                                first = homo + first
+                            else:
+                                first = homo + 1 + first
+                        orbitals.append(first)
+
+        if len(tasks) == 0:
+            return ""
+
+        lines = []
+        txt = "', '".join(tasks)
+        lines.append("")
+        lines.append("# Cube files for density and orbitals")
+        lines.append(f"set cubeprop_tasks ['{txt}']")
+        if len(orbitals) > 0:
+            txt = ", ".join([f"{i}, {-i}" for i in orbitals])
+            lines.append(f"set cubeprop_orbitals [{txt}]")
+
+        lines.append("")
+        lines.append("cubeprop(wfn)")
+        lines.append(
+            f"""
+# Prefix the files with the substep number
+paths = Path.cwd().glob('*.cube')
+for path in paths:
+    name = path.name
+    newpath = path.with_name('@{self._id[-1]}+' + name)
+    path.rename(newpath)
+"""
+        )
+
+        return "\n".join(lines)
diff --git a/psi4_step/energy_parameters.py b/psi4_step/energy_parameters.py
@@ -272,10 +272,45 @@ class EnergyParameters(seamm.Parameters):
         },
     }
 
+    output = {
+        "density": {
+            "default": "no",
+            "kind": "boolean",
+            "default_units": "",
+            "enumeration": ("yes", "no"),
+            "format_string": "",
+            "description": "Plot total density:",
+            "help_text": "Whether to plot the total charge density.",
+        },
+        "orbitals": {
+            "default": "no",
+            "kind": "boolean",
+            "default_units": "",
+            "enumeration": ("yes", "no"),
+            "format_string": "",
+            "description": "Plot orbitals:",
+            "help_text": "Whether to plot orbitals.",
+        },
+        "selected orbitals": {
+            "default": "-1, HOMO, LUMO, +1",
+            "kind": "string",
+            "default_units": "",
+            "enumeration": ("all", "-1, HOMO, LUMO, +1"),
+            "format_string": "",
+            "description": "Selected orbitals:",
+            "help_text": "Which orbitals to plot.",
+        },
+    }
+
     def __init__(self, defaults={}, data=None):
         """Initialize the instance, by default from the default
         parameters given in the class"""
 
         super().__init__(
-            defaults={**EnergyParameters.parameters, **defaults}, data=data
+            defaults={
+                **EnergyParameters.parameters,
+                **EnergyParameters.output,
+                **defaults,
+            },
+            data=data,
         )
diff --git a/psi4_step/initialization.py b/psi4_step/initialization.py
@@ -87,7 +87,7 @@ def get_input(self):
 
         result.append(f"set basis {P['basis']}")
         result.append("")
-        # result.append(f"initial.symmetrize({P['symmetry_tolerance']})")
+        result.append(f"initial.symmetrize({P['symmetry_tolerance']})")
         result.append("point_group = initial.point_group().symbol()")
 
         # Dump the properties to a json file

diff --git a/psi4_step/initialization_parameters.py b/psi4_step/initialization_parameters.py
@@ -18,18 +18,33 @@ class InitializationParameters(seamm.Parameters):
             "default_units": "",
             "enumeration": (
                 "6-31G",
-                "6-31G*",
-                "6-31G**",
+                "6-31G(d)",
+                "6-31G(d,p)",
+                "6-31+G",
+                "6-31+G(d)",
+                "6-31+G(d,p)",
+                "6-311G",
+                "6-311G(d)",
+                "6-311G(d,p)",
+                "6-311+G",
+                "6-311+G(d)",
+                "6-311+G(d,p)",
                 "cc-pVDZ",
                 "cc-pVTZ",
                 "cc-pVQZ",
+                "def2-SV(P)",
+                "def2-SVP",
+                "def2-TZVP",
+                "def2-TZVPP",
+                "def2-QZVP",
+                "def2-QZVPP",
             ),
             "format_string": "s",
             "description": "Basis:",
             "help_text": ("The basis set to use."),
         },
         "symmetry_tolerance": {
-            "default": "0.05",
+            "default": "0.00001",
             "kind": "float",
             "default_units": None,
             "enumeration": tuple(),

diff --git a/psi4_step/optimization.py b/psi4_step/optimization.py
@@ -19,7 +19,13 @@
 
 
 class Optimization(psi4_step.Energy):
-    def __init__(self, flowchart=None, title="Optimization", extension=None):
+    def __init__(
+        self,
+        flowchart=None,
+        title="Optimization",
+        extension=None,
+        logger=logger,
+    ):
         """Initialize the node"""
 
         logger.debug("Creating Optimization {}".format(self))