Merge pull request #3 from molssi-seamm/dev

Updated to be compatible with general trajectory analysis

HISTORY.rst

@@ -2,6 +2,8 @@
 History
 =======
 
+2023.5.29 -- Converged with general approach for trajectory analysis
+
 2023.5.6 -- Bugfix
     * Fixed an error handling Nan's and Inf's that caused a crash
     * Added the predictions from the derivatives of the Helfand moments to the output to

thermal_conductivity_step/thermal_conductivity.py

@@ -3,6 +3,7 @@
 """Non-graphical part of the Thermal Conductivity step in a SEAMM flowchart
 """
 
+import json
 import logging
 from math import log10, ceil
 from pathlib import Path
@@ -11,6 +12,7 @@
 import traceback
 
 import numpy as np
+import pandas
 from tabulate import tabulate
 
 from .analysis import (
@@ -29,7 +31,7 @@
 import molsystem
 import seamm
 import seamm_util
-from seamm_util import ureg, Q_  # noqa: F401
+from seamm_util import Q_
 import seamm_util.printing as printing
 from seamm_util.printing import FormattedText as __
 
@@ -792,29 +794,75 @@ def process_run(self, run, run_dir):
         run_dir : pathlib.Path
             The toplevel directory of the run.
         """
-        paths = sorted(run_dir.glob("**/HeatFlux.npz"))
+        paths = sorted(run_dir.glob("**/heat_flux.trj"))
 
         if len(paths) == 0:
-            raise RuntimeError(f"There is no heat flux data for run {run}.")
+            raise RuntimeError(
+                f"There is no heat flux data for run {run} in {run_dir}."
+            )
         elif len(paths) > 1:
             raise NotImplementedError(
-                f"Cannot handle multiple HeatFlux files from run {run}."
+                f"Cannot handle multiple HeatFlux files from run {run} in {run_dir}."
+            )
+
+        # Process the trajectory data
+        control_properties = lambda x: x not in ["tstep"]  # noqa: E731
+        with paths[0].open() as fd:
+            # Get the initial header line and check
+            line = fd.readline()
+            tmp = line.split()
+            if len(tmp) < 3:
+                raise RuntimeError(f"Bad header for {paths[0]}: {line}")
+            if tmp[0] != "!MolSSI":
+                raise RuntimeError(f"Not a MolSSI file? {paths[0]}: {line}")
+            if tmp[1] != "trajectory":
+                raise RuntimeError(f"Not a trajectory file? {paths[0]}: {line}")
+            if tmp[2][0] != "2":
+                raise RuntimeError(
+                    f"Can only handle version 2 trajectory files. {paths[0]}: {line}"
+                )
+            metadata = json.loads(" ".join(tmp[3:]))
+
+            data = pandas.read_csv(
+                fd,
+                sep=" ",
+                header=0,
+                comment="!",
+                usecols=control_properties,
+                index_col=None,
             )
-        with np.load(paths[0]) as data:
-            V = float(data["V"])
-            T = float(data["T"])
-            timestep = float(data["timestep"])
 
-            J = data["J"]
+        dt = Q_(metadata["dt"], metadata["tunits"])
+        data = data.reset_index(drop=True)
+        data.index *= dt.m_as("fs")
 
-        self.V.append(V)
-        self.T.append(T)
-        self.timestep.append(timestep)
+        Jx = data["Jx"].to_numpy()
+        Jy = data["Jy"].to_numpy()
+        Jz = data["Jz"].to_numpy()
+
+        J = np.stack((Jx, Jy, Jz))
+
+        # Get the state function V & T
+        path = paths[0].parent / "state.json"
+        if path.exists():
+            with path.open() as fd:
+                state = json.load(fd)
+            T = state["T"]["mean"]
+            if "V" in state:
+                V = state["V"]["mean"]
+            else:
+                _, configuration = self.get_system_configuration()
+                V = configuration.volume
+            self.V.append(V)
+            self.T.append(T)
+        else:
+            self.V.append(None)
+            self.T.append(None)
+        self.timestep.append(dt.magnitude)
 
         # We need properties like the temperature and volume.
         T = Q_(T, "K")
         V = Q_(V, "Å^3")
-        dt = Q_(timestep, "fs")
         k_B = Q_("k_B")
         Jsq = Q_("W^2/m^4")
 

thermal_conductivity_step/tk_thermal_conductivity.py

@@ -255,6 +255,6 @@ def right_click(self, event):
         """
 
         super().right_click(event)
-        self.popup_menu.add_command(label="Edit..", command=self.edit)
+        self.popup_menu.add_command(label="Edit...", command=self.edit)
 
         self.popup_menu.tk_popup(event.x_root, event.y_root, 0)