From 357e109878ac599eb1d26e21f57dce36851775a6 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:06:53 -0800
Subject: [PATCH 01/46] Init lammps feature

---
 src/mpmorph/flows/md_flow.py |  5 +--
 src/mpmorph/jobs/lammps.py   | 59 ++++++++++++++++++++++++++++++++++++
 2 files changed, 62 insertions(+), 2 deletions(-)
 create mode 100644 src/mpmorph/jobs/lammps.py

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index e7bbc016..c30f7ba3 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -12,10 +12,11 @@
 M3GNET_MD_FLOW = "M3GNET_MD_FLOW"
 M3GNET_MD_CONVERGED_VOL_FLOW = "M3GNET_MD_CONVERGED_VOL_FLOW"
 
-def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vol_scale = 1):
+def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
     inputs = M3GNetMDInputs(
         temperature=temp,
-        steps=steps
+        steps=steps,
+        **input_kwargs
     )
 
     m3gnet_maker = M3GNetMDMaker(parameters = inputs)
diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
new file mode 100644
index 00000000..59a5fb99
--- /dev/null
+++ b/src/mpmorph/jobs/lammps.py
@@ -0,0 +1,59 @@
+from jobflow import Maker, job
+from pymatgen.io.lammps.utils import LammpsRunner
+from pymatgen.io.lammps.inputs import La
+import logging
+import os
+
+from pymatgen.io.lammps.inputs import TemplateInputGen, LammpsTemplateGen
+from pymatgen.io.lammps.outputs import LammpsDump, parse_lammps_log
+from pymatgen.io.lammps.data import LammpsData
+from pymatgen.core.structure import Structure
+
+class RunLammpsMaker(Maker):
+    """
+    Run LAMMPS directly (no custodian).
+    Required params:
+        lammsps_cmd (str): lammps command to run sans the input file name.
+            e.g. 'mpirun -n 4 lmp_mpi'
+    """
+
+    @job
+    def make(self, lammps_cmd: str,
+                   script_template_path: str,
+                   script_options: dict,
+                   structure: Structure = None,
+                   data_filename: str = None,
+                   log_filename: str = "log.lammps"):
+
+        if data_filename is None:
+            data = LammpsData.from_structure(structure)
+        else:
+            data = LammpsData.from_file(data_filename)
+        # Write the input files
+        linp = LammpsTemplateGen().get_input_set(script_template=script_template_path,
+                                                 settings=script_options, 
+                                                 data=data,
+                                                 data_filename='data.dump')
+
+        input_name = f'lammps.in'
+        linp.write_input(input_name)
+
+        # Run LAMMPS
+        lmps_runner = LammpsRunner(input_name, lammps_cmd)
+        stdout, stderr = lmps_runner.run()
+        logging.info(f"LAMMPS finished running: {stdout} \n {stderr}")
+
+        dump_files = dump_files or []
+        dump_files = [dump_files] if isinstance(dump_files, str) else dump_files
+
+        # Construct various dumps objects
+        dumps = []
+        if dump_files:
+            for df in dump_files:
+                dumps.append((df, LammpsDump.from_file(df)))
+
+        # Construct log object
+        log = parse_lammps_log(log_filename)
+
+        logging.info(f"Getting task doc for base dir :{path}")
+        d = self.generate_doc(path, lmps_input, log, dumps)        
\ No newline at end of file

From 121fb65069d0a9ca732aa5e78405eb71be7f05ba Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:17:29 -0800
Subject: [PATCH 02/46] add output

---
 src/mpmorph/jobs/lammps.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 59a5fb99..8dade216 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -1,10 +1,8 @@
 from jobflow import Maker, job
 from pymatgen.io.lammps.utils import LammpsRunner
-from pymatgen.io.lammps.inputs import La
 import logging
-import os
 
-from pymatgen.io.lammps.inputs import TemplateInputGen, LammpsTemplateGen
+from pymatgen.io.lammps.inputs import LammpsTemplateGen
 from pymatgen.io.lammps.outputs import LammpsDump, parse_lammps_log
 from pymatgen.io.lammps.data import LammpsData
 from pymatgen.core.structure import Structure
@@ -55,5 +53,7 @@ def make(self, lammps_cmd: str,
         # Construct log object
         log = parse_lammps_log(log_filename)
 
-        logging.info(f"Getting task doc for base dir :{path}")
-        d = self.generate_doc(path, lmps_input, log, dumps)        
\ No newline at end of file
+        return {
+            "dumps": dumps,
+            "log": log
+        }
\ No newline at end of file

From b325272be88efe04767760ab6eeefbbda03fd4a2 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:22:39 -0800
Subject: [PATCH 03/46] Add name to lammps maker

---
 src/mpmorph/jobs/lammps.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 8dade216..d78e545a 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -15,6 +15,8 @@ class RunLammpsMaker(Maker):
             e.g. 'mpirun -n 4 lmp_mpi'
     """
 
+    name = "RUN_LAMMPS"
+
     @job
     def make(self, lammps_cmd: str,
                    script_template_path: str,

From 72bd5ae306ff4964deed4b207f72d8f6097fde01 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:24:32 -0800
Subject: [PATCH 04/46] Add dump files parameter

---
 src/mpmorph/jobs/lammps.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index d78e545a..ff212352 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -2,6 +2,8 @@
 from pymatgen.io.lammps.utils import LammpsRunner
 import logging
 
+from typing import List
+
 from pymatgen.io.lammps.inputs import LammpsTemplateGen
 from pymatgen.io.lammps.outputs import LammpsDump, parse_lammps_log
 from pymatgen.io.lammps.data import LammpsData
@@ -23,7 +25,8 @@ def make(self, lammps_cmd: str,
                    script_options: dict,
                    structure: Structure = None,
                    data_filename: str = None,
-                   log_filename: str = "log.lammps"):
+                   log_filename: str = "log.lammps",
+                   dump_files: List[str] = None):
 
         if data_filename is None:
             data = LammpsData.from_structure(structure)

From bd34ffa1e804fbf3943f6ce0333a7ba6b2ffc37a Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:33:49 -0800
Subject: [PATCH 05/46] fix data filename

---
 src/mpmorph/jobs/lammps.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index ff212352..1b4a7900 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -24,19 +24,16 @@ def make(self, lammps_cmd: str,
                    script_template_path: str,
                    script_options: dict,
                    structure: Structure = None,
-                   data_filename: str = None,
                    log_filename: str = "log.lammps",
+                   data_filename: str = "data.lammps",
                    dump_files: List[str] = None):
 
-        if data_filename is None:
-            data = LammpsData.from_structure(structure)
-        else:
-            data = LammpsData.from_file(data_filename)
+        data = LammpsData.from_structure(structure)
         # Write the input files
         linp = LammpsTemplateGen().get_input_set(script_template=script_template_path,
                                                  settings=script_options, 
                                                  data=data,
-                                                 data_filename='data.dump')
+                                                 data_filename=data_filename)
 
         input_name = f'lammps.in'
         linp.write_input(input_name)

From 3ee4b4b99474b97d0e53e3293877d1f1ac8adecb Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 14:48:03 -0800
Subject: [PATCH 06/46] try print

---
 src/mpmorph/jobs/lammps.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 1b4a7900..95dd9d4e 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -41,7 +41,7 @@ def make(self, lammps_cmd: str,
         # Run LAMMPS
         lmps_runner = LammpsRunner(input_name, lammps_cmd)
         stdout, stderr = lmps_runner.run()
-        logging.info(f"LAMMPS finished running: {stdout} \n {stderr}")
+        print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
         dump_files = dump_files or []
         dump_files = [dump_files] if isinstance(dump_files, str) else dump_files

From 2f6032dbd44b9c011f2e1ba9beec83dff62e961f Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Thu, 19 Jan 2023 15:04:38 -0800
Subject: [PATCH 07/46] Add atom style to lammps data

---
 src/mpmorph/jobs/lammps.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 95dd9d4e..438a3874 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -28,7 +28,7 @@ def make(self, lammps_cmd: str,
                    data_filename: str = "data.lammps",
                    dump_files: List[str] = None):
 
-        data = LammpsData.from_structure(structure)
+        data = LammpsData.from_structure(structure, atom_style='atomic')
         # Write the input files
         linp = LammpsTemplateGen().get_input_set(script_template=script_template_path,
                                                  settings=script_options, 

From 7338ed478fc8d05c9445379d42289ded4c43eda5 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 07:32:10 -0800
Subject: [PATCH 08/46] Improve lammps runner inputs

---
 src/mpmorph/jobs/lammps.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 438a3874..29eeb5e7 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -1,3 +1,4 @@
+from subprocess import PIPE, Popen
 from jobflow import Maker, job
 from pymatgen.io.lammps.utils import LammpsRunner
 import logging
@@ -35,12 +36,17 @@ def make(self, lammps_cmd: str,
                                                  data=data,
                                                  data_filename=data_filename)
 
-        input_name = f'lammps.in'
-        linp.write_input(input_name)
-
+        linp.write_input(directory=".")
+        input_name = "in.lammps"
         # Run LAMMPS
         lmps_runner = LammpsRunner(input_name, lammps_cmd)
         stdout, stderr = lmps_runner.run()
+
+        lammps_cmd = self.lammps_bin + ["-in", input_name]
+        print(f"Running: {' '.join(lammps_cmd)}")
+        with Popen(lammps_cmd, stdout=PIPE, stderr=PIPE) as p:
+            (stdout, stderr) = p.communicate()
+
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
         dump_files = dump_files or []

From e9ccfeb006d097ab4446e7e819bc57e340d01e8a Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 07:36:10 -0800
Subject: [PATCH 09/46] fix lammps invocation

---
 src/mpmorph/jobs/lammps.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index 29eeb5e7..ebce5e5a 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -21,7 +21,7 @@ class RunLammpsMaker(Maker):
     name = "RUN_LAMMPS"
 
     @job
-    def make(self, lammps_cmd: str,
+    def make(self, lammps_bin: str,
                    script_template_path: str,
                    script_options: dict,
                    structure: Structure = None,
@@ -39,10 +39,8 @@ def make(self, lammps_cmd: str,
         linp.write_input(directory=".")
         input_name = "in.lammps"
         # Run LAMMPS
-        lmps_runner = LammpsRunner(input_name, lammps_cmd)
-        stdout, stderr = lmps_runner.run()
 
-        lammps_cmd = self.lammps_bin + ["-in", input_name]
+        lammps_cmd = lammps_bin + ["-in", input_name]
         print(f"Running: {' '.join(lammps_cmd)}")
         with Popen(lammps_cmd, stdout=PIPE, stderr=PIPE) as p:
             (stdout, stderr) = p.communicate()

From b1128867cf4e2a6484e3ff7c722ea2fcff5bd5b6 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 07:38:27 -0800
Subject: [PATCH 10/46] fix input format

---
 src/mpmorph/jobs/lammps.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index ebce5e5a..ab4caba6 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -40,7 +40,7 @@ def make(self, lammps_bin: str,
         input_name = "in.lammps"
         # Run LAMMPS
 
-        lammps_cmd = lammps_bin + ["-in", input_name]
+        lammps_cmd = [lammps_bin, "-in", input_name]
         print(f"Running: {' '.join(lammps_cmd)}")
         with Popen(lammps_cmd, stdout=PIPE, stderr=PIPE) as p:
             (stdout, stderr) = p.communicate()

From 161998d85853601bd73819762b0a80fcb02f65d7 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 12:30:19 -0800
Subject: [PATCH 11/46] Add template to package

---
 src/mpmorph/flows/md_flow.py                  | 12 ++++++
 .../jobs/lammps-templates/template.lammps     | 40 +++++++++++++++++++
 src/mpmorph/jobs/lammps.py                    | 38 +++++++-----------
 src/mpmorph/jobs/pv_from_calc.py              | 16 ++++++++
 src/mpmorph/schemas/pv_data_doc.py            |  2 +-
 5 files changed, 83 insertions(+), 25 deletions(-)
 create mode 100644 src/mpmorph/jobs/lammps-templates/template.lammps

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index c30f7ba3..c417b9d1 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -29,6 +29,18 @@ def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vo
         initial_vol_scale=initial_vol_scale
     )
 
+def get_equil_vol_flow_lammps(structure, temp, steps):
+    inputs = M3GNetMDInputs(
+        temperature=temp,
+        steps=steps
+    )
+
+    pv_md_maker = PVFromM3GNet(parameters=inputs)
+    eq_vol_maker = EquilibriumVolumeSearchMaker(pv_md_maker=pv_md_maker)
+    equil_vol_job = eq_vol_maker.make(structure)
+    flow = Flow([equil_vol_job], output=equil_vol_job.output, name=EQUILIBRATE_VOLUME_FLOW)
+    return flow
+
 def get_equil_vol_flow(structure, temp, steps):
     inputs = M3GNetMDInputs(
         temperature=temp,
diff --git a/src/mpmorph/jobs/lammps-templates/template.lammps b/src/mpmorph/jobs/lammps-templates/template.lammps
new file mode 100644
index 00000000..5bdc3723
--- /dev/null
+++ b/src/mpmorph/jobs/lammps-templates/template.lammps
@@ -0,0 +1,40 @@
+#
+# This is an example of the driver of `M3GNet' <https://github.com/materialsvirtuallab/m3gnet>,
+# which contains a state-of-the-art Graph Neural Network Potential trained with data of Materials Projects.
+# This driver is developed by AdvanceSoft Corp <https://www.advancesoft.jp>.
+# Before you use this driver, you have to install python3 and m3gnet (pip install m3gnet).
+#
+# NOTE:
+#   1) the units must be metal
+#   2) the 3D periodic boundary condition must be used
+#   3) cannot use MPI parallelization, only OpenMP or GPU
+#
+
+units         metal
+boundary      p p p
+atom_style    atomic
+
+read_data       data.lammps 
+
+pair_style    m3gnet $m3gnet_path
+pair_coeff    * *  MP-2021.2.8-EFS  $species
+
+
+# Don't need actual positions
+# dump          myDump all custom 10 xyz.lammpstrj id element x y z
+# dump_modify   myDump sort id element $species
+
+thermo_style  custom step time cpu pe ke etotal temp press vol density
+thermo        10
+
+variable p1 equal "step"
+variable p2 equal "temp"
+variable p3 equal "vol"
+variable p4 equal "density"
+fix def1 all print $print_every_n_step "${p1} ${p2} ${p3} ${p4}" file step_temp_vol_density.txt 
+
+
+velocity      all create $temperature 12345
+fix           myEnse all npt temp $temperature $temperature 0.1 aniso 1.0 1.0 1.0
+timestep      0.002
+run           $total_steps 
diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps.py
index ab4caba6..10009a72 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps.py
@@ -1,16 +1,17 @@
 from subprocess import PIPE, Popen
 from jobflow import Maker, job
-from pymatgen.io.lammps.utils import LammpsRunner
-import logging
 
-from typing import List
+import pandas as pd
 
 from pymatgen.io.lammps.inputs import LammpsTemplateGen
-from pymatgen.io.lammps.outputs import LammpsDump, parse_lammps_log
 from pymatgen.io.lammps.data import LammpsData
 from pymatgen.core.structure import Structure
 
-class RunLammpsMaker(Maker):
+from mpmorph.schemas.pv_data_doc import MDPVDataDoc
+
+from pkg_resources import resource_filename
+
+class LammpsVolMaker(Maker):
     """
     Run LAMMPS directly (no custodian).
     Required params:
@@ -22,16 +23,16 @@ class RunLammpsMaker(Maker):
 
     @job
     def make(self, lammps_bin: str,
-                   script_template_path: str,
                    script_options: dict,
                    structure: Structure = None,
-                   log_filename: str = "log.lammps",
-                   data_filename: str = "data.lammps",
-                   dump_files: List[str] = None):
+                   data_filename: str = "data.lammps"):
+        
+        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
+
 
         data = LammpsData.from_structure(structure, atom_style='atomic')
         # Write the input files
-        linp = LammpsTemplateGen().get_input_set(script_template=script_template_path,
+        linp = LammpsTemplateGen().get_input_set(script_template=template_path,
                                                  settings=script_options, 
                                                  data=data,
                                                  data_filename=data_filename)
@@ -47,19 +48,8 @@ def make(self, lammps_bin: str,
 
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
-        dump_files = dump_files or []
-        dump_files = [dump_files] if isinstance(dump_files, str) else dump_files
-
-        # Construct various dumps objects
-        dumps = []
-        if dump_files:
-            for df in dump_files:
-                dumps.append((df, LammpsDump.from_file(df)))
 
-        # Construct log object
-        log = parse_lammps_log(log_filename)
+        filecontents = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", skiprows=1, index_col="step", names=["step", "temp", "vol", "density"])
+        eq_vol = filecontents[["vol"]].iloc[-1].values[0]
 
-        return {
-            "dumps": dumps,
-            "log": log
-        }
\ No newline at end of file
+        return eq_vol
\ No newline at end of file
diff --git a/src/mpmorph/jobs/pv_from_calc.py b/src/mpmorph/jobs/pv_from_calc.py
index 518cced7..57f7bb92 100644
--- a/src/mpmorph/jobs/pv_from_calc.py
+++ b/src/mpmorph/jobs/pv_from_calc.py
@@ -44,6 +44,22 @@ def build_doc(self, m3gnet_calc: M3GNetMDCalculation):
         p_data = m3gnet_calc_to_pressure(m3gnet_calc)
         return MDPVDataDoc(volume=v_data, pressure=p_data)
 
+@dataclass
+class PVFromM3GNetLammps(PVFromCalc):
+    """Generates a MDPVDataDoc using Lammps run with M3gnet and a npt ensemble.
+    """
+
+    name: str = "PV_FROM_M3GNET_LAMMPS"
+    parameters: M3GNetMDInputs = None
+
+    def run_md(self, structure: Structure, **kwargs):
+        calc_doc = run_m3gnet(structure, self.parameters, self.name, **kwargs)
+
+        return calc_doc
+
+    def build_doc(self, pvdoc: MDPVDataDoc):
+        return pvdoc
+
 
 def m3gnet_calc_to_vol(m3gnet_calc: M3GNetMDCalculation):
     volume = m3gnet_calc.trajectory[-1].lattice.volume
diff --git a/src/mpmorph/schemas/pv_data_doc.py b/src/mpmorph/schemas/pv_data_doc.py
index 7219d8ea..bb40b89b 100644
--- a/src/mpmorph/schemas/pv_data_doc.py
+++ b/src/mpmorph/schemas/pv_data_doc.py
@@ -5,4 +5,4 @@ class MDPVDataDoc(BaseModel):
 
     task_label: str = Field(None, description="The name of the task.")
     volume: float = Field(None, description="The volume data from the MD run")
-    pressure: float = Field(None, description="The volume data from the MD run")
+    pressure: float = Field(None, description="The pressure of the MD run")

From c3744dd3403343d5050cbfe2b171a76d97e6370c Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 12:51:29 -0800
Subject: [PATCH 12/46] Try vol workflow w lammps

---
 src/mpmorph/flows/md_flow.py                  | 24 +++++++++-----
 src/mpmorph/flows/vt_flow.py                  | 33 ++++++++++++++++++-
 .../jobs/lammps-templates/template.lammps     |  9 ++---
 .../jobs/{lammps.py => lammps_volume.py}      | 22 +++++++++----
 4 files changed, 65 insertions(+), 23 deletions(-)
 rename src/mpmorph/jobs/{lammps.py => lammps_volume.py} (74%)

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index c417b9d1..3a0cff21 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -3,6 +3,7 @@
 from mpmorph.jobs.core import M3GNetMDMaker
 
 from mpmorph.jobs.equilibrate_volume import EquilibriumVolumeSearchMaker
+from mpmorph.jobs.lammps_volume import LammpsVolMaker
 from pymatgen.core.structure import Structure
 
 from mpmorph.jobs.pv_from_calc import PVFromCalc, PVFromM3GNet, PVFromVasp
@@ -11,6 +12,7 @@
 EQUILIBRATE_VOLUME_FLOW = "EQUILIBRATE_VOLUME_FLOW"
 M3GNET_MD_FLOW = "M3GNET_MD_FLOW"
 M3GNET_MD_CONVERGED_VOL_FLOW = "M3GNET_MD_CONVERGED_VOL_FLOW"
+LAMMPS_VOL_FLOW = "LAMMPS_VOL_FLOW"
 
 def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
     inputs = M3GNetMDInputs(
@@ -29,16 +31,20 @@ def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vo
         initial_vol_scale=initial_vol_scale
     )
 
-def get_equil_vol_flow_lammps(structure, temp, steps):
-    inputs = M3GNetMDInputs(
-        temperature=temp,
-        steps=steps
+def get_equil_vol_flow_lammps(structure,
+                              temp,
+                              steps,
+                              lammps_bin_path,
+                              m3gnet_path):
+    vol_maker = LammpsVolMaker()
+    vol_job = vol_maker.make(
+        lammps_bin_path,
+        temp,
+        m3gnet_path,
+        steps,
+        structure
     )
-
-    pv_md_maker = PVFromM3GNet(parameters=inputs)
-    eq_vol_maker = EquilibriumVolumeSearchMaker(pv_md_maker=pv_md_maker)
-    equil_vol_job = eq_vol_maker.make(structure)
-    flow = Flow([equil_vol_job], output=equil_vol_job.output, name=EQUILIBRATE_VOLUME_FLOW)
+    flow = Flow([vol_job], output=vol_job, name=LAMMPS_VOL_FLOW)
     return flow
 
 def get_equil_vol_flow(structure, temp, steps):
diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index 4b789363..abd54dee 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -1,7 +1,7 @@
 from jobflow import Flow, job
 import json
 
-from .md_flow import get_equil_vol_flow
+from .md_flow import get_equil_vol_flow, get_equil_vol_flow_lammps
 from ..jobs.pv_from_calc import m3gnet_calc_to_vol
 
 VOLUME_TEMPERATURE_SWEEP = "VOLUME_TEMPERATURE_SWEEP"
@@ -33,6 +33,37 @@ def get_vt_sweep_flow(
     new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
     return new_flow
 
+def get_vt_sweep_flow_lammps(
+    structure,
+    lammps_bin_path,
+    m3gnet_path,
+    lower_bound=100,
+    upper_bound=1100,
+    temp_step=100,
+    output_name="vt.out",
+    steps=2000,
+):
+
+    vs = []
+    volume_jobs = []
+    temps = list(range(lower_bound, upper_bound, temp_step))
+
+    for temp in temps:
+        job = get_equil_vol_flow_lammps(
+            structure=structure,
+            temp=temp,
+            steps=steps,
+            lammps_bin_path=lammps_bin_path,
+            m3gnet_path=m3gnet_path
+        )
+        volume_jobs.append(job)
+        vs.append(job.output)
+
+    collect_job = _collect_vt_results(vs, temps, structure, output_name)
+
+    new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
+    return new_flow
+
 
 @job
 def _collect_vt_results(vs, ts, structure, output_fn):
diff --git a/src/mpmorph/jobs/lammps-templates/template.lammps b/src/mpmorph/jobs/lammps-templates/template.lammps
index 5bdc3723..045569be 100644
--- a/src/mpmorph/jobs/lammps-templates/template.lammps
+++ b/src/mpmorph/jobs/lammps-templates/template.lammps
@@ -14,18 +14,13 @@ units         metal
 boundary      p p p
 atom_style    atomic
 
-read_data       data.lammps 
+read_data     data.lammps 
 
 pair_style    m3gnet $m3gnet_path
 pair_coeff    * *  MP-2021.2.8-EFS  $species
 
-
-# Don't need actual positions
-# dump          myDump all custom 10 xyz.lammpstrj id element x y z
-# dump_modify   myDump sort id element $species
-
 thermo_style  custom step time cpu pe ke etotal temp press vol density
-thermo        10
+thermo        $print_every_n_step
 
 variable p1 equal "step"
 variable p2 equal "temp"
diff --git a/src/mpmorph/jobs/lammps.py b/src/mpmorph/jobs/lammps_volume.py
similarity index 74%
rename from src/mpmorph/jobs/lammps.py
rename to src/mpmorph/jobs/lammps_volume.py
index 10009a72..0e8aec72 100644
--- a/src/mpmorph/jobs/lammps.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -13,23 +13,33 @@
 
 class LammpsVolMaker(Maker):
     """
-    Run LAMMPS directly (no custodian).
+    Run LAMMPS directly using m3gnet (no custodian).
     Required params:
         lammsps_cmd (str): lammps command to run sans the input file name.
             e.g. 'mpirun -n 4 lmp_mpi'
     """
 
-    name = "RUN_LAMMPS"
+    name = "LAMMPS_TO_VOLUME"
 
     @job
     def make(self, lammps_bin: str,
-                   script_options: dict,
-                   structure: Structure = None,
-                   data_filename: str = "data.lammps"):
+                   temperature: int,
+                   m3gnet_path: str,
+                   total_steps: int,
+                   structure: Structure = None):
+
+        script_options = {
+            "temperature": temperature,
+            "m3gnet_path": m3gnet_path,
+            "species": structure.composition.chemical_system.replace("-", " "),
+            "total_steps": total_steps,
+            "print_every_n_step": 10
+        }                   
         
-        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
 
+        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
 
+        data_filename: str = "data.lammps"
         data = LammpsData.from_structure(structure, atom_style='atomic')
         # Write the input files
         linp = LammpsTemplateGen().get_input_set(script_template=template_path,

From 28e601298ca91aa1ae1f578449a510684059df0e Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 13:03:39 -0800
Subject: [PATCH 13/46] fix output ref

---
 src/mpmorph/flows/vt_flow.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index abd54dee..d018915f 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -57,7 +57,7 @@ def get_vt_sweep_flow_lammps(
             m3gnet_path=m3gnet_path
         )
         volume_jobs.append(job)
-        vs.append(job.output)
+        vs.append(job.output.output)
 
     collect_job = _collect_vt_results(vs, temps, structure, output_name)
 

From 538320efbbfe2c4152c07afe105bf81e730ed6e3 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 13:13:23 -0800
Subject: [PATCH 14/46] fix chem sys string

---
 src/mpmorph/jobs/lammps_volume.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 0e8aec72..f91b8c85 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -28,10 +28,11 @@ def make(self, lammps_bin: str,
                    total_steps: int,
                    structure: Structure = None):
 
+        chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
         script_options = {
             "temperature": temperature,
             "m3gnet_path": m3gnet_path,
-            "species": structure.composition.chemical_system.replace("-", " "),
+            "species": chem_sys_str,
             "total_steps": total_steps,
             "print_every_n_step": 10
         }                   

From 7406ddeaf177ee5081814a41ccfa44822835442c Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 14:41:25 -0800
Subject: [PATCH 15/46] get bins from path

---
 src/mpmorph/jobs/lammps_volume.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index f91b8c85..3fec543b 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -1,4 +1,6 @@
 from subprocess import PIPE, Popen
+
+import os
 from jobflow import Maker, job
 
 import pandas as pd
@@ -22,12 +24,13 @@ class LammpsVolMaker(Maker):
     name = "LAMMPS_TO_VOLUME"
 
     @job
-    def make(self, lammps_bin: str,
-                   temperature: int,
-                   m3gnet_path: str,
+    def make(self, temperature: int,
                    total_steps: int,
                    structure: Structure = None):
 
+        lammps_bin = os.environ.get("LAMMPS_CMD")
+        m3gnet_path = os.environ.get("M3GNET_PATH")
+
         chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
         script_options = {
             "temperature": temperature,

From e215b089ca1ed50e05fed29c1f259a694889f6a0 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 14:44:22 -0800
Subject: [PATCH 16/46] remove extraneous params

---
 src/mpmorph/flows/md_flow.py | 6 +-----
 src/mpmorph/flows/vt_flow.py | 4 ----
 2 files changed, 1 insertion(+), 9 deletions(-)

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index 3a0cff21..b25a4394 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -33,14 +33,10 @@ def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vo
 
 def get_equil_vol_flow_lammps(structure,
                               temp,
-                              steps,
-                              lammps_bin_path,
-                              m3gnet_path):
+                              steps):
     vol_maker = LammpsVolMaker()
     vol_job = vol_maker.make(
-        lammps_bin_path,
         temp,
-        m3gnet_path,
         steps,
         structure
     )
diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index d018915f..3b50c68a 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -35,8 +35,6 @@ def get_vt_sweep_flow(
 
 def get_vt_sweep_flow_lammps(
     structure,
-    lammps_bin_path,
-    m3gnet_path,
     lower_bound=100,
     upper_bound=1100,
     temp_step=100,
@@ -53,8 +51,6 @@ def get_vt_sweep_flow_lammps(
             structure=structure,
             temp=temp,
             steps=steps,
-            lammps_bin_path=lammps_bin_path,
-            m3gnet_path=m3gnet_path
         )
         volume_jobs.append(job)
         vs.append(job.output.output)

From b7d0da0854f16cfe32f3622a75bf8b4ed0f60f42 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 16:08:18 -0800
Subject: [PATCH 17/46] add mpid and formula to output

---
 src/mpmorph/flows/vt_flow.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index 3b50c68a..7e05ff82 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -40,6 +40,7 @@ def get_vt_sweep_flow_lammps(
     temp_step=100,
     output_name="vt.out",
     steps=2000,
+    mp_id=None
 ):
 
     vs = []
@@ -55,15 +56,21 @@ def get_vt_sweep_flow_lammps(
         volume_jobs.append(job)
         vs.append(job.output.output)
 
-    collect_job = _collect_vt_results(vs, temps, structure, output_name)
+    collect_job = _collect_vt_results(vs, temps, structure, output_name, mp_id)
 
     new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
     return new_flow
 
 
 @job
-def _collect_vt_results(vs, ts, structure, output_fn):
-    result = {"structure": structure.as_dict(), "volumes": vs, "temps": ts}
+def _collect_vt_results(vs, ts, structure, output_fn, mp_id):
+    result = {
+        "structure": structure.as_dict(),
+        "volumes": vs,
+        "temps": ts,
+        "mp_id": mp_id,
+        "formula": structure.composition.reduced_formula
+    }
 
     with open(output_fn, "+w") as f:
         f.write(json.dumps(result))

From 5a91e2265392f653e4f100995794072d162f0f1e Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 16:28:48 -0800
Subject: [PATCH 18/46] use average of last 10% of data points

---
 src/mpmorph/jobs/lammps_volume.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 3fec543b..04d73da1 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -63,7 +63,8 @@ def make(self, temperature: int,
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
 
-        filecontents = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", skiprows=1, index_col="step", names=["step", "temp", "vol", "density"])
-        eq_vol = filecontents[["vol"]].iloc[-1].values[0]
+        avging_window = int(min(total_steps / 100, 50))
+        df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", skiprows=1, names=["step", "temp", "vol", "density"])
+        eq_vol = df.iloc[-avging_window::]['vol'].values.mean()
 
         return eq_vol
\ No newline at end of file

From ba41df6f48ce3e23195209eccb3a384c862d4404 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 18:05:19 -0800
Subject: [PATCH 19/46] take 30%

---
 src/mpmorph/jobs/lammps_volume.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 04d73da1..94b584c3 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -63,7 +63,7 @@ def make(self, temperature: int,
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
 
-        avging_window = int(min(total_steps / 100, 50))
+        avging_window = int(total_steps / 300)
         df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", skiprows=1, names=["step", "temp", "vol", "density"])
         eq_vol = df.iloc[-avging_window::]['vol'].values.mean()
 

From 7248d802073b5666e40e2a0c37e7ab866595ee04 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 18:35:50 -0800
Subject: [PATCH 20/46] Save more data than we have been

---
 src/mpmorph/flows/vt_flow.py      | 25 +++++++++++++++++++------
 src/mpmorph/jobs/lammps_volume.py |  7 +++----
 2 files changed, 22 insertions(+), 10 deletions(-)

diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index 7e05ff82..0ae8cb0a 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -1,8 +1,10 @@
 from jobflow import Flow, job
 import json
+import uuid
 
 from .md_flow import get_equil_vol_flow, get_equil_vol_flow_lammps
 from ..jobs.pv_from_calc import m3gnet_calc_to_vol
+import pandas as pd
 
 VOLUME_TEMPERATURE_SWEEP = "VOLUME_TEMPERATURE_SWEEP"
 
@@ -43,7 +45,7 @@ def get_vt_sweep_flow_lammps(
     mp_id=None
 ):
 
-    vs = []
+    v_outputs = []
     volume_jobs = []
     temps = list(range(lower_bound, upper_bound, temp_step))
 
@@ -54,24 +56,35 @@ def get_vt_sweep_flow_lammps(
             steps=steps,
         )
         volume_jobs.append(job)
-        vs.append(job.output.output)
+        v_outputs.append(job.output.output)
 
-    collect_job = _collect_vt_results(vs, temps, structure, output_name, mp_id)
+    collect_job = _collect_vt_results(v_outputs, temps, structure, output_name, mp_id)
 
     new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
     return new_flow
 
 
 @job
-def _collect_vt_results(vs, ts, structure, output_fn, mp_id):
+def _collect_vt_results(v_outputs, ts, structure, output_fn, mp_id):
     result = {
         "structure": structure.as_dict(),
-        "volumes": vs,
+        "volumes": [get_converged_vol(v) for v in v_outputs],
         "temps": ts,
         "mp_id": mp_id,
-        "formula": structure.composition.reduced_formula
+        "reduced_formula": structure.composition.reduced_formula,
+        "formula": structure.composition.formula,
+        "uuid": str(uuid.uuid4())
     }
 
     with open(output_fn, "+w") as f:
         f.write(json.dumps(result))
     return result
+
+def get_converged_vol(v_output):   
+    df = pd.DataFrame.from_dict(v_output)
+    total_steps = (len(df) - 1) * 10
+    avging_window = int(total_steps / 30)
+    vols = df.iloc[-avging_window::]['vol']
+    eq_vol = vols.values.mean()
+    return float(eq_vol)
+
diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 94b584c3..8b1f2c0d 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -63,8 +63,7 @@ def make(self, temperature: int,
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
 
-        avging_window = int(total_steps / 300)
-        df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", skiprows=1, names=["step", "temp", "vol", "density"])
-        eq_vol = df.iloc[-avging_window::]['vol'].values.mean()
+        
+        df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", index_col="step", skiprows=1, names=["step", "temp", "vol", "density"])
 
-        return eq_vol
\ No newline at end of file
+        return df.to_dict()
\ No newline at end of file

From 9cf017c7080c15ec49f72b277ae7058d4f944070 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 20 Jan 2023 20:06:58 -0800
Subject: [PATCH 21/46] fix flow name

---
 src/mpmorph/flows/vt_flow.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index 0ae8cb0a..2771f747 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -60,7 +60,9 @@ def get_vt_sweep_flow_lammps(
 
     collect_job = _collect_vt_results(v_outputs, temps, structure, output_name, mp_id)
 
-    new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
+
+    flow_name = f'{structure.composition.reduced_formula}-Melting Point'
+    new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=flow_name)
     return new_flow
 
 

From d16f39865b41fd4868463f5c2b6486029b537393 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 27 Jan 2023 15:07:07 -0800
Subject: [PATCH 22/46] Update LAMMPS output

---
 src/mpmorph/analysis/melting_points.py        | 99 +++++++++++++++----
 .../jobs/lammps-templates/template.lammps     |  4 +
 src/mpmorph/jobs/lammps_volume.py             | 39 ++++++--
 src/mpmorph/jobs/volume_temperature_sweep.py  | 25 -----
 src/mpmorph/schemas/lammps_calc.py            | 34 +++++++
 5 files changed, 150 insertions(+), 51 deletions(-)
 create mode 100644 src/mpmorph/schemas/lammps_calc.py

diff --git a/src/mpmorph/analysis/melting_points.py b/src/mpmorph/analysis/melting_points.py
index 2b42e97f..17bd8da5 100644
--- a/src/mpmorph/analysis/melting_points.py
+++ b/src/mpmorph/analysis/melting_points.py
@@ -1,28 +1,51 @@
 import numpy as np
 from scipy.stats import linregress
 import matplotlib.pyplot as plt
-
+from sklearn.cluster import AgglomerativeClustering
+from sklearn.metrics import mean_squared_error
+import numpy as np
+import matplotlib.pyplot as plt
 import math
 
-class MeltingPointAnalyzer():
+class MeltingPointClusterAnalyzer():
 
-    def split_dset(self, pts, split_idx):
-        return pts[0:split_idx], pts[split_idx:]
+    def _get_clusters(self, points):
+        clustering = AgglomerativeClustering(n_clusters=2).fit(points)
+        cluster1 = points[np.argwhere(clustering.labels_ == 1).squeeze()].T
+        cluster2 = points[np.argwhere(clustering.labels_ == 0).squeeze()].T
+        return cluster1, cluster2
 
-    def get_split_fit(self, xs, ys, split_idx):
-        leftx, rightx = self.split_dset(xs, split_idx)
-        lefty, righty = self.split_dset(ys, split_idx)
-        
-        leftfit = linregress(leftx, lefty)
-        lefterr = leftfit.stderr
-        
-        rightfit = linregress(rightx, righty)
-        righterr = rightfit.stderr
-        
-        combined_err = math.sqrt(lefterr ** 2 + righterr ** 2)
-        combined_err = lefterr + righterr
-        return leftfit.slope, leftfit.intercept, rightfit.slope, rightfit.intercept, combined_err
+    def plot_vol_vs_temp(self, ts, vs, plot_title = None):
+        points = np.array(list(zip(ts, vs)))
+        cluster1, cluster2 = self._get_clusters(points)
+        plt.scatter(*cluster1)
+        plt.scatter(*cluster2)
+        plt.xlabel("Temperature (K)")
+        plt.ylabel("Volume (A^3)")
+        Tm = self.estimate_melting_temp(ts, vs)
+        plt.plot([Tm, Tm], [min(vs), max(vs)], color='r')
+      
+        if plot_title is None:
+            plt.title("Volume vs Temperature by Clustering")
+        else:
+            plt.title(plot_title)
+    
+    def estimate_melting_temp(self, temps, vols):
+        points = np.array(list(zip(temps, vols)))
+        cluster1, cluster2 = self._get_clusters(points)
+        if min(cluster1[0]) < min(cluster2[0]):
+            solid_range = cluster1[0]
+            liquid_range = cluster2[0]
+        else:
+            solid_range = cluster2[0]
+            liquid_range = cluster1[0]
+
+        return np.mean([max(solid_range), min(liquid_range)])
+
+class MeltingPointSlopeAnalyzer():
 
+    def split_dset(self, pts, split_idx):
+        return pts[0:split_idx], pts[split_idx:]
 
     def assess_splits(self, xs, ys):
         dset_size = len(xs)
@@ -48,7 +71,7 @@ def plot_split(self, xs, ys, split_idx):
 
         plt.scatter(leftxs, leftys)
         plt.plot(leftxs, left_fit_ys)
-        plt.title("Volume vs Temperature (w/ best fits)")
+        plt.title("Volume vs Temperature (w/ best fits by Slope Method")
         plt.xlabel("Temperature (K)")
         plt.ylabel("Equil. Volume (cubic Angstroms)")
 
@@ -67,7 +90,45 @@ def get_best_split(self, xs, ys):
     def plot_vol_vs_temp(self, temps, vols):
         split_idx = self.get_best_split(temps, vols)
         self.plot_split(temps, vols, split_idx)
+        Tm = self.estimate_melting_temp(temps, vols)
+        print(Tm)
+        plt.plot([Tm, Tm], [min(vols), max(vols)], color='r')        
+
 
     def estimate_melting_temp(self, temps, vols):
         best_split_idx = self.get_best_split(temps, vols)
-        return temps[best_split_idx]
\ No newline at end of file
+        return np.mean([temps[best_split_idx], temps[best_split_idx - 1]])
+
+class MeltingPointSlopeRMSEAnalyzer(MeltingPointSlopeAnalyzer):
+
+    def get_split_fit(self, xs, ys, split_idx):
+        leftx, rightx = self.split_dset(xs, split_idx)
+        lefty, righty = self.split_dset(ys, split_idx)
+        
+        lslope, lintercept, r_value, p_value, std_err = linregress(leftx, lefty)
+        left_y_pred = lintercept + lslope * np.array(leftx)
+        lefterr = mean_squared_error(y_true=lefty, y_pred=left_y_pred, squared=False)
+
+        rslope, rintercept, r_value, p_value, std_err = linregress(rightx, righty)
+        right_y_pred = rintercept + rslope * np.array(rightx)
+        righterr = mean_squared_error(y_true=righty, y_pred=right_y_pred, squared=False)
+        
+        combined_err = math.sqrt(lefterr ** 2 + righterr ** 2)
+        combined_err = lefterr + righterr
+        return lslope, lintercept, rslope, rintercept, combined_err
+
+class MeltingPointSlopeStdErrAnalyzer(MeltingPointSlopeAnalyzer):
+
+    def get_split_fit(self, xs, ys, split_idx):
+        leftx, rightx = self.split_dset(xs, split_idx)
+        lefty, righty = self.split_dset(ys, split_idx)
+        
+        leftfit = linregress(leftx, lefty)
+        lefterr = leftfit.stderr
+        
+        rightfit = linregress(rightx, righty)
+        righterr = rightfit.stderr
+        
+        combined_err = math.sqrt(lefterr ** 2 + righterr ** 2)
+        combined_err = lefterr + righterr
+        return leftfit.slope, leftfit.intercept, rightfit.slope, rightfit.intercept, combined_err
diff --git a/src/mpmorph/jobs/lammps-templates/template.lammps b/src/mpmorph/jobs/lammps-templates/template.lammps
index 045569be..ab2a161a 100644
--- a/src/mpmorph/jobs/lammps-templates/template.lammps
+++ b/src/mpmorph/jobs/lammps-templates/template.lammps
@@ -22,6 +22,10 @@ pair_coeff    * *  MP-2021.2.8-EFS  $species
 thermo_style  custom step time cpu pe ke etotal temp press vol density
 thermo        $print_every_n_step
 
+# Record trajectory
+dump          myDump all custom 10 trajectory.lammpstrj id element x y z
+dump_modify   myDump sort id element $species
+
 variable p1 equal "step"
 variable p2 equal "temp"
 variable p3 equal "vol"
diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 8b1f2c0d..ea70d753 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -8,12 +8,15 @@
 from pymatgen.io.lammps.inputs import LammpsTemplateGen
 from pymatgen.io.lammps.data import LammpsData
 from pymatgen.core.structure import Structure
+from ase.io.lammpsrun import read_lammps_dump_text
+from pymatgen.io.ase import AseAtomsAdaptor
+from pymatgen.core.trajectory import Trajectory
 
-from mpmorph.schemas.pv_data_doc import MDPVDataDoc
+from mpmorph.schemas.lammps_calc import LammpsCalc
 
 from pkg_resources import resource_filename
 
-class LammpsVolMaker(Maker):
+class LammpsCalcMaker(Maker):
     """
     Run LAMMPS directly using m3gnet (no custodian).
     Required params:
@@ -21,9 +24,9 @@ class LammpsVolMaker(Maker):
             e.g. 'mpirun -n 4 lmp_mpi'
     """
 
-    name = "LAMMPS_TO_VOLUME"
+    name = "LAMMPS_CALCULATION"
 
-    @job
+    @job(trajectory="trajectory", output_schema=LammpsCalc)
     def make(self, temperature: int,
                    total_steps: int,
                    structure: Structure = None):
@@ -38,7 +41,7 @@ def make(self, temperature: int,
             "species": chem_sys_str,
             "total_steps": total_steps,
             "print_every_n_step": 10
-        }                   
+        }
         
 
         template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
@@ -62,8 +65,30 @@ def make(self, temperature: int,
 
         print(f"LAMMPS finished running: {stdout} \n {stderr}")
 
+        # Build trajectory from LAMMPS output .xyz file
+        with open("trajectory.lammpstrj", "r+") as f:
+            atoms = read_lammps_dump_text(f, index=slice(-1))
+
+        structs = []
+
+        for a in atoms:
+            structs.append(AseAtomsAdaptor().get_structure(a))
+
+        trajectory = Trajectory.from_structures(structs, constant_lattice=False)
 
-        
         df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", index_col="step", skiprows=1, names=["step", "temp", "vol", "density"])
 
-        return df.to_dict()
\ No newline at end of file
+        metadata = {
+            "temperature": temperature,
+            "total_steps": total_steps
+        }
+
+        output = LammpsCalc(
+            dir_name=os.getcwd(),
+            trajectory=trajectory,
+            composition=structure.composition,
+            reduced_formula=structure.composition.reduced_formula,
+            metadata=metadata,
+            dump_data=df.to_dict()
+        )
+        return output
\ No newline at end of file
diff --git a/src/mpmorph/jobs/volume_temperature_sweep.py b/src/mpmorph/jobs/volume_temperature_sweep.py
index c9a9683b..dd8aa53c 100644
--- a/src/mpmorph/jobs/volume_temperature_sweep.py
+++ b/src/mpmorph/jobs/volume_temperature_sweep.py
@@ -1,9 +1,5 @@
 from jobflow import Flow, Maker, Response, job
-import json
-
 from mpmorph.jobs.tasks.m3gnet_input import M3GNetMDInputs
-from ..schemas.vt_sweep_doc import VTSweepDoc
-
 import dataclasses
 
 class VolumeTemperatureSweepMaker(Maker):
@@ -41,24 +37,3 @@ def make(
         return Response(replace=new_flow)
 
 
-@job
-def _collect_vt_results(vs, ts, structure, output_fn = None):
-    filtered_vs = []
-    filtered_ts = []
-    for v, t in zip(vs, ts):
-        if v is not None:
-            filtered_vs.append(v)
-            filtered_ts.append(t)
-
-
-    result = VTSweepDoc(
-        volumes=filtered_vs,
-        temps=filtered_ts,
-        structure=structure
-    )
-
-    if output_fn is not None:
-        with open(output_fn, "+w") as f:
-            f.write(json.dumps(result))
-        
-    return result
diff --git a/src/mpmorph/schemas/lammps_calc.py b/src/mpmorph/schemas/lammps_calc.py
new file mode 100644
index 00000000..918086e3
--- /dev/null
+++ b/src/mpmorph/schemas/lammps_calc.py
@@ -0,0 +1,34 @@
+from pydantic import BaseModel, Field
+from pymatgen.core.composition import Composition
+from pymatgen.core.trajectory import Trajectory as PmgTrajectory
+
+from mpmorph.utils import datetime_str
+
+
+class LammpsCalc(BaseModel):
+    task_label: str = Field(None, description="The name of the task.")
+    dir_name: str = Field(
+        None, description="The directory where the LAMMPS calculation was run"
+    )
+    last_updated: str = Field(
+        default_factory=datetime_str,
+        description="Timestamp of when the document was last updated.",
+    )
+    trajectory: PmgTrajectory = Field(
+        None, description="The pymatgen Trajectory object stored ad dictionary"
+    )
+    composition: Composition = Field(description="The composition of the structure.")
+    reduced_formula: str = Field(
+        description="The reduced formula of the structure's composition."
+    )
+    dump_data: dict = Field(
+        None,
+        description="Any additional data collected via LAMMPS dump files"
+    )
+    metadata: dict = Field(
+        None,
+        description=(
+            "Important info about the calculation, including ensemble type,"
+            " temperature, etc."
+        ),
+    )

From 38843a5ea2e2d53c278742a1ed705a6849b3ea76 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 27 Jan 2023 15:08:00 -0800
Subject: [PATCH 23/46] fix typo"

---
 src/mpmorph/flows/md_flow.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index b25a4394..a193cdde 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -3,7 +3,7 @@
 from mpmorph.jobs.core import M3GNetMDMaker
 
 from mpmorph.jobs.equilibrate_volume import EquilibriumVolumeSearchMaker
-from mpmorph.jobs.lammps_volume import LammpsVolMaker
+from mpmorph.jobs.lammps_volume import LammpsCalcMaker
 from pymatgen.core.structure import Structure
 
 from mpmorph.jobs.pv_from_calc import PVFromCalc, PVFromM3GNet, PVFromVasp
@@ -34,7 +34,7 @@ def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vo
 def get_equil_vol_flow_lammps(structure,
                               temp,
                               steps):
-    vol_maker = LammpsVolMaker()
+    vol_maker = LammpsCalcMaker()
     vol_job = vol_maker.make(
         temp,
         steps,

From 3971037ad5054b8f243b69aaf6c3c78e8025fc7b Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Fri, 27 Jan 2023 15:40:36 -0800
Subject: [PATCH 24/46] include last step of trajectory

---
 src/mpmorph/jobs/lammps_volume.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index ea70d753..2598ef09 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -67,7 +67,7 @@ def make(self, temperature: int,
 
         # Build trajectory from LAMMPS output .xyz file
         with open("trajectory.lammpstrj", "r+") as f:
-            atoms = read_lammps_dump_text(f, index=slice(-1))
+            atoms = read_lammps_dump_text(f, index=slice(0, None))
 
         structs = []
 

From c70678e921c53ceead4c2026ae46fd055dc52397 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 8 Feb 2023 10:09:36 -0800
Subject: [PATCH 25/46] add temperature sweeping lammps input template

---
 .../jobs/lammps-templates/in.melt_template    | 118 ++++++++++++++++++
 1 file changed, 118 insertions(+)
 create mode 100644 src/mpmorph/jobs/lammps-templates/in.melt_template

diff --git a/src/mpmorph/jobs/lammps-templates/in.melt_template b/src/mpmorph/jobs/lammps-templates/in.melt_template
new file mode 100644
index 00000000..0381b7c3
--- /dev/null
+++ b/src/mpmorph/jobs/lammps-templates/in.melt_template
@@ -0,0 +1,118 @@
+# Script originally made by Oscar Guerrero
+# Reference: https://orca.cardiff.ac.uk/id/eprint/101322/1/MRSPaper2.pdf
+units  metal
+atom_style atomic
+boundary p p p
+
+atom_modify map array
+
+#define variables
+
+variable tempstart equal $tempstart
+variable tempstop equal $tempstop 
+variable myseed equal 12345
+variable atomrate equal 1000
+variable time_step equal 0.002
+variable time_eq equal 1000
+#variable tdamp equal 1.
+
+variable tamp equal "v_time_step*1000" # DO NOT CHANGE
+variable pdamp equal "v_time_step*1000" # DO NOT CHANGE
+timestep ${time_step} # DO NOT CHANGE
+
+
+#Create structure
+read_data data.dump
+
+
+#Define Interatomic Potential
+
+pair_style    m3gnet /global/home/users/huizheng/repos/lammps/potentials/M3GNET
+pair_coeff    * *  MP-2021.2.8-EFS   $species
+
+# Equilibration
+reset_timestep 0
+velocity all create ${tempstart} ${myseed} mom yes rot no dist gaussian
+fix equilibration all npt temp ${tempstart} ${tempstart} $(100.0*dt) iso 1 1 ${pdamp} drag 0.2
+
+
+variable eq1 equal "step"
+variable eq2 equal "pxx"
+variable eq3 equal "pyy"
+variable eq4 equal "pzz"
+variable eq5 equal "lx"
+variable eq6 equal "ly"
+variable eq7 equal "lz"
+variable eq8 equal "vol"
+variable eq9 equal "temp"
+variable eq10 equal "etotal"
+
+fix data_equilibration all print 10 "${eq1} ${eq2} ${eq3} ${eq4} ${eq5} ${eq6} ${eq7} ${eq8} ${eq9} ${eq10}" file ${tempstart}K.data
+thermo 1000
+thermo_style custom step pxx pyy pzz lx ly lz temp etotal
+
+# RUN
+run 1000
+
+# store final volume Vo to calculate V/Vo (reduce units)
+variable tmp equal "vol"
+variable Vo equal ${tmp}
+print "Volume initial is , Vo: ${Vo}"
+
+#reset
+unfix equilibration
+unfix data_equilibration
+
+#----------------------------- Increase temperature------------------------------------
+reset_timestep 0
+fix melting all npt temp ${tempstart} ${tempstop} $(100.0*dt) iso 1 1 ${pdamp} drag 0.2
+# fix melting all nvt temp ${tempstart} ${tempstop} ${tdamp} drag 0.2
+
+variable eq1 equal "step"
+variable eq2 equal "pxx"
+variable eq3 equal "pyy"
+variable eq4 equal "pzz"
+variable eq5 equal "lx"
+variable eq6 equal "ly"
+variable eq7 equal "lz"
+variable eq8 equal "temp"
+variable eq9 equal "vol/v_Vo"
+variable eq10 equal "etotal"
+run 0 
+fix data_melting all print $print_every_n_step "${eq8} ${eq9}" file temp_vs_ref_vol.txt screen no
+
+dump 1 all cfg 100 HgF2.step*.cfg mass type xs ys zs id
+dump_modify 1 element  $species
+dump 2 all custom  100  dump.* id type x y z
+# Compute msd command and dump every 10 steps
+compute msd all msd com yes
+fix msd all ave/time 1 1 10 c_msd[4] file msd.txt
+
+
+# use velocity auto-correlation function (VACF) to calculate diffusion coefficient
+compute         2 all vacf
+fix             5 all vector 1 c_2[4]
+variable        diff equal dt*trap(f_5)
+fix    vacf all print 1 "${eq1} ${eq8} ${eq9} ${diff}"
+
+run $total_steps 
+thermo $print_every_n_step
+thermo_style    custom step v_diff
+
+
+# print to screen out
+
+#run 10000
+#reset 
+unfix melting 
+unfix data_melting 
+undump 1
+undump 2
+write_restart restart.equil
+write_data data.*
+# SAVE THE DATA OF THE CALCULATION OR ELSE YOU NEED TO START OVER = ( OUCH !
+# SIMULATION DONE
+clear
+print "You've done great job =)"
+
+

From 9e9af360532c80931b50653917047848c3e9656a Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Fri, 10 Feb 2023 15:37:14 -0800
Subject: [PATCH 26/46] run black

---
 src/mpmorph/analysis/diffusion.py            |  2 +-
 src/mpmorph/analysis/melting_points.py       | 47 +++++++++++---------
 src/mpmorph/analysis/structural_analysis.py  |  3 --
 src/mpmorph/database.py                      |  3 +-
 src/mpmorph/firetasks/dbtasks.py             |  1 -
 src/mpmorph/fireworks/powerups.py            |  6 +--
 src/mpmorph/flows/md_flow.py                 |  6 +++
 src/mpmorph/flows/vt_flow.py                 | 33 +++++++-------
 src/mpmorph/io.py                            |  8 ++--
 src/mpmorph/jobs/equilibrate_volume.py       |  1 -
 src/mpmorph/jobs/lammps_volume.py            | 42 +++++++++--------
 src/mpmorph/jobs/pv_from_calc.py             |  6 +--
 src/mpmorph/jobs/tasks/m3gnet_input.py       |  1 -
 src/mpmorph/jobs/volume_temperature_sweep.py |  4 +-
 src/mpmorph/schemas/lammps_calc.py           |  3 +-
 src/mpmorph/schemas/m3gnet_md_calc.py        |  1 -
 src/mpmorph/schemas/pv_data_doc.py           |  1 -
 src/mpmorph/schemas/vt_sweep_doc.py          |  9 ++--
 src/mpmorph/workflows/quench.py              |  2 +-
 19 files changed, 92 insertions(+), 87 deletions(-)

diff --git a/src/mpmorph/analysis/diffusion.py b/src/mpmorph/analysis/diffusion.py
index 879572ad..bb06bd9a 100644
--- a/src/mpmorph/analysis/diffusion.py
+++ b/src/mpmorph/analysis/diffusion.py
@@ -290,7 +290,7 @@ def plot(self, title=None, annotate=True, el="", **kwargs):
             self.y,
             yerr=self.yerr.T,
             label="Q[{}]: ".format(el) + tx + " K",
-            **kwargs
+            **kwargs,
         )
         plt.ylabel("ln(D cm$^2$/s)", fontsize=15)
         plt.xlabel("1000/T K$^{-1}$", fontsize=15)
diff --git a/src/mpmorph/analysis/melting_points.py b/src/mpmorph/analysis/melting_points.py
index 17bd8da5..a3ff8724 100644
--- a/src/mpmorph/analysis/melting_points.py
+++ b/src/mpmorph/analysis/melting_points.py
@@ -7,15 +7,15 @@
 import matplotlib.pyplot as plt
 import math
 
-class MeltingPointClusterAnalyzer():
 
+class MeltingPointClusterAnalyzer:
     def _get_clusters(self, points):
         clustering = AgglomerativeClustering(n_clusters=2).fit(points)
         cluster1 = points[np.argwhere(clustering.labels_ == 1).squeeze()].T
         cluster2 = points[np.argwhere(clustering.labels_ == 0).squeeze()].T
         return cluster1, cluster2
 
-    def plot_vol_vs_temp(self, ts, vs, plot_title = None):
+    def plot_vol_vs_temp(self, ts, vs, plot_title=None):
         points = np.array(list(zip(ts, vs)))
         cluster1, cluster2 = self._get_clusters(points)
         plt.scatter(*cluster1)
@@ -23,13 +23,13 @@ def plot_vol_vs_temp(self, ts, vs, plot_title = None):
         plt.xlabel("Temperature (K)")
         plt.ylabel("Volume (A^3)")
         Tm = self.estimate_melting_temp(ts, vs)
-        plt.plot([Tm, Tm], [min(vs), max(vs)], color='r')
-      
+        plt.plot([Tm, Tm], [min(vs), max(vs)], color="r")
+
         if plot_title is None:
             plt.title("Volume vs Temperature by Clustering")
         else:
             plt.title(plot_title)
-    
+
     def estimate_melting_temp(self, temps, vols):
         points = np.array(list(zip(temps, vols)))
         cluster1, cluster2 = self._get_clusters(points)
@@ -42,8 +42,8 @@ def estimate_melting_temp(self, temps, vols):
 
         return np.mean([max(solid_range), min(liquid_range)])
 
-class MeltingPointSlopeAnalyzer():
 
+class MeltingPointSlopeAnalyzer:
     def split_dset(self, pts, split_idx):
         return pts[0:split_idx], pts[split_idx:]
 
@@ -56,7 +56,7 @@ def assess_splits(self, xs, ys):
         for idx in pt_idxs:
             _, _, _, _, total_err = self.get_split_fit(xs, ys, idx)
             errs.append(total_err)
-        
+
         return list(zip(pt_idxs, errs))
 
     def get_linear_ys(self, m, b, xs):
@@ -79,32 +79,31 @@ def plot_split(self, xs, ys, split_idx):
 
         plt.scatter(rightxs, rightys)
         plt.plot(rightxs, right_fit_ys)
-        
+
     def get_best_split(self, xs, ys):
         split_errs = self.assess_splits(xs, ys)
         errs = [pt[1] for pt in split_errs]
         idxs = [pt[0] for pt in split_errs]
         best_split_idx = idxs[np.argmin(errs)]
         return best_split_idx
-    
+
     def plot_vol_vs_temp(self, temps, vols):
         split_idx = self.get_best_split(temps, vols)
         self.plot_split(temps, vols, split_idx)
         Tm = self.estimate_melting_temp(temps, vols)
         print(Tm)
-        plt.plot([Tm, Tm], [min(vols), max(vols)], color='r')        
-
+        plt.plot([Tm, Tm], [min(vols), max(vols)], color="r")
 
     def estimate_melting_temp(self, temps, vols):
         best_split_idx = self.get_best_split(temps, vols)
         return np.mean([temps[best_split_idx], temps[best_split_idx - 1]])
 
-class MeltingPointSlopeRMSEAnalyzer(MeltingPointSlopeAnalyzer):
 
+class MeltingPointSlopeRMSEAnalyzer(MeltingPointSlopeAnalyzer):
     def get_split_fit(self, xs, ys, split_idx):
         leftx, rightx = self.split_dset(xs, split_idx)
         lefty, righty = self.split_dset(ys, split_idx)
-        
+
         lslope, lintercept, r_value, p_value, std_err = linregress(leftx, lefty)
         left_y_pred = lintercept + lslope * np.array(leftx)
         lefterr = mean_squared_error(y_true=lefty, y_pred=left_y_pred, squared=False)
@@ -112,23 +111,29 @@ def get_split_fit(self, xs, ys, split_idx):
         rslope, rintercept, r_value, p_value, std_err = linregress(rightx, righty)
         right_y_pred = rintercept + rslope * np.array(rightx)
         righterr = mean_squared_error(y_true=righty, y_pred=right_y_pred, squared=False)
-        
-        combined_err = math.sqrt(lefterr ** 2 + righterr ** 2)
+
+        combined_err = math.sqrt(lefterr**2 + righterr**2)
         combined_err = lefterr + righterr
         return lslope, lintercept, rslope, rintercept, combined_err
 
-class MeltingPointSlopeStdErrAnalyzer(MeltingPointSlopeAnalyzer):
 
+class MeltingPointSlopeStdErrAnalyzer(MeltingPointSlopeAnalyzer):
     def get_split_fit(self, xs, ys, split_idx):
         leftx, rightx = self.split_dset(xs, split_idx)
         lefty, righty = self.split_dset(ys, split_idx)
-        
+
         leftfit = linregress(leftx, lefty)
         lefterr = leftfit.stderr
-        
+
         rightfit = linregress(rightx, righty)
         righterr = rightfit.stderr
-        
-        combined_err = math.sqrt(lefterr ** 2 + righterr ** 2)
+
+        combined_err = math.sqrt(lefterr**2 + righterr**2)
         combined_err = lefterr + righterr
-        return leftfit.slope, leftfit.intercept, rightfit.slope, rightfit.intercept, combined_err
+        return (
+            leftfit.slope,
+            leftfit.intercept,
+            rightfit.slope,
+            rightfit.intercept,
+            combined_err,
+        )
diff --git a/src/mpmorph/analysis/structural_analysis.py b/src/mpmorph/analysis/structural_analysis.py
index 2becc70f..ef2884cf 100644
--- a/src/mpmorph/analysis/structural_analysis.py
+++ b/src/mpmorph/analysis/structural_analysis.py
@@ -47,7 +47,6 @@ def polyhedra_connectivity(structures, pair, cutoff, step_freq=1):
             polyhedra_list.append(set(current_poly))
 
         for polypair in itertools.combinations(polyhedra_list, 2):
-
             polyhedra_pair_type = (len(polypair[0]), len(polypair[1]))
 
             shared_vertices = len(polypair[0].intersection(polypair[1]))
@@ -143,7 +142,6 @@ class BondAngleDistribution(object):
     """
 
     def __init__(self, structures, cutoffs, step_freq=1):
-
         self.bond_angle_distribution = None
         self.structures = structures
         self.step_freq = step_freq
@@ -251,7 +249,6 @@ def get_bond_angle_distribution(self):
 
                 # get all pair combinations of neoghbor sites of i:
                 for p in itertools.combinations(neighbors[i], 2):
-
                     # check if pairs are within the defined cutoffs
                     if self._cutoff_type == "dict":
                         if self._check_skip_triplet(s_index, i, p[0][2], p[1][2]):
diff --git a/src/mpmorph/database.py b/src/mpmorph/database.py
index 20684cba..272fff15 100644
--- a/src/mpmorph/database.py
+++ b/src/mpmorph/database.py
@@ -30,7 +30,7 @@ def __init__(
         collection="tasks",
         user=None,
         password=None,
-        **kwargs
+        **kwargs,
     ):
         super(VaspMDCalcDb, self).__init__(
             host, port, database, collection, user, password, **kwargs
@@ -75,7 +75,6 @@ def insert_task(
 
         # insert structures  at each ionic step into GridFS
         if parse_ionic_steps and "calcs_reversed" in task_doc:
-
             # Convert from ionic steps dictionary to pymatgen.core.trajectory.Trajectory object
             ionic_steps_dict = task_doc["calcs_reversed"][0]["output"]["ionic_steps"]
             time_step = task_doc["input"]["incar"]["POTIM"]
diff --git a/src/mpmorph/firetasks/dbtasks.py b/src/mpmorph/firetasks/dbtasks.py
index 1e07f82d..0f660f99 100644
--- a/src/mpmorph/firetasks/dbtasks.py
+++ b/src/mpmorph/firetasks/dbtasks.py
@@ -222,7 +222,6 @@ def load_trajectories_from_gfs(runs, mmdb, gfs_keys=None):
 
     trajectory = None
     for i, (fs_id, fs) in enumerate(gfs_keys):
-
         if fs == "trajectories_fs" or fs == "rebuild_trajectories_fs":
             # Load stored Trajectory
             print(fs_id, "is stored in trajectories_fs")
diff --git a/src/mpmorph/fireworks/powerups.py b/src/mpmorph/fireworks/powerups.py
index 22f16257..4d244c33 100644
--- a/src/mpmorph/fireworks/powerups.py
+++ b/src/mpmorph/fireworks/powerups.py
@@ -46,7 +46,7 @@ def aggregate_trajectory(fw, **kwargs):
 def add_cont_structure(fw):
     prev_struct_task = PreviousStructureTask()
     insert_i = 2
-    for (i, task) in enumerate(fw.tasks):
+    for i, task in enumerate(fw.tasks):
         if task.fw_name == "{{atomate.vasp.firetasks.run_calc.RunVaspCustodian}}":
             insert_i = i
             break
@@ -68,7 +68,7 @@ def add_pass_pv(fw, **kwargs):
 
 def add_pv_volume_rescale(fw):
     insert_i = 2
-    for (i, task) in enumerate(fw.tasks):
+    for i, task in enumerate(fw.tasks):
         if task.fw_name == "{{atomate.vasp.firetasks.run_calc.RunVaspCustodian}}":
             insert_i = i
             break
@@ -80,7 +80,7 @@ def add_pv_volume_rescale(fw):
 def add_rescale_volume(fw, **kwargs):
     rsv_task = RescaleVolumeTask(**kwargs)
     insert_i = 2
-    for (i, task) in enumerate(fw.tasks):
+    for i, task in enumerate(fw.tasks):
         if task.fw_name == "{{atomate.vasp.firetasks.run_calc.RunVaspCustodian}}":
             insert_i = i
             break
diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index a193cdde..b63aa6c6 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -14,6 +14,10 @@
 M3GNET_MD_CONVERGED_VOL_FLOW = "M3GNET_MD_CONVERGED_VOL_FLOW"
 LAMMPS_VOL_FLOW = "LAMMPS_VOL_FLOW"
 
+
+def get_md_temperature_sweeping(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
+
+
 def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
     inputs = M3GNetMDInputs(
         temperature=temp,
@@ -91,3 +95,5 @@ def _get_converge_flow(structure: Structure, pv_md_maker: PVFromCalc, production
     flow = Flow([equil_vol_job, final_md_job], output=final_md_job.output, name=M3GNET_MD_CONVERGED_VOL_FLOW)
 
     return flow
+
+
diff --git a/src/mpmorph/flows/vt_flow.py b/src/mpmorph/flows/vt_flow.py
index 2771f747..7c1ee123 100644
--- a/src/mpmorph/flows/vt_flow.py
+++ b/src/mpmorph/flows/vt_flow.py
@@ -8,6 +8,7 @@
 
 VOLUME_TEMPERATURE_SWEEP = "VOLUME_TEMPERATURE_SWEEP"
 
+
 def get_vt_sweep_flow(
     structure,
     lower_bound=100,
@@ -16,25 +17,25 @@ def get_vt_sweep_flow(
     output_name="vt.out",
     steps=2000,
 ):
-
     vs = []
     volume_jobs = []
     temps = list(range(lower_bound, upper_bound, temp_step))
 
     for temp in temps:
-        job = get_equil_vol_flow(
-            structure=structure,
-            temp=temp,
-            steps=steps
-        )
+        job = get_equil_vol_flow(structure=structure, temp=temp, steps=steps)
         volume_jobs.append(job)
         vs.append(job.output.volume)
 
     collect_job = _collect_vt_results(vs, temps, structure, output_name)
 
-    new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=VOLUME_TEMPERATURE_SWEEP)
+    new_flow = Flow(
+        [*volume_jobs, collect_job],
+        output=collect_job.output,
+        name=VOLUME_TEMPERATURE_SWEEP,
+    )
     return new_flow
 
+
 def get_vt_sweep_flow_lammps(
     structure,
     lower_bound=100,
@@ -42,9 +43,8 @@ def get_vt_sweep_flow_lammps(
     temp_step=100,
     output_name="vt.out",
     steps=2000,
-    mp_id=None
+    mp_id=None,
 ):
-
     v_outputs = []
     volume_jobs = []
     temps = list(range(lower_bound, upper_bound, temp_step))
@@ -60,9 +60,10 @@ def get_vt_sweep_flow_lammps(
 
     collect_job = _collect_vt_results(v_outputs, temps, structure, output_name, mp_id)
 
-
-    flow_name = f'{structure.composition.reduced_formula}-Melting Point'
-    new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output, name=flow_name)
+    flow_name = f"{structure.composition.reduced_formula}-Melting Point"
+    new_flow = Flow(
+        [*volume_jobs, collect_job], output=collect_job.output, name=flow_name
+    )
     return new_flow
 
 
@@ -75,18 +76,18 @@ def _collect_vt_results(v_outputs, ts, structure, output_fn, mp_id):
         "mp_id": mp_id,
         "reduced_formula": structure.composition.reduced_formula,
         "formula": structure.composition.formula,
-        "uuid": str(uuid.uuid4())
+        "uuid": str(uuid.uuid4()),
     }
 
     with open(output_fn, "+w") as f:
         f.write(json.dumps(result))
     return result
 
-def get_converged_vol(v_output):   
+
+def get_converged_vol(v_output):
     df = pd.DataFrame.from_dict(v_output)
     total_steps = (len(df) - 1) * 10
     avging_window = int(total_steps / 30)
-    vols = df.iloc[-avging_window::]['vol']
+    vols = df.iloc[-avging_window::]["vol"]
     eq_vol = vols.values.mean()
     return float(eq_vol)
-
diff --git a/src/mpmorph/io.py b/src/mpmorph/io.py
index d095ecbe..8a43529c 100644
--- a/src/mpmorph/io.py
+++ b/src/mpmorph/io.py
@@ -17,13 +17,13 @@ def get_string_from_struct(
     ):
         format_str = "{{:.{0}f}}".format(significant_figures)
 
-        for (si, structure) in enumerate(structures):
+        for si, structure in enumerate(structures):
             lines = [system, "1.0", str(structure.lattice)]
             lines.append(" ".join(self.get_site_symbols(structure)))
             lines.append(" ".join([str(x) for x in self.get_natoms(structure)]))
 
             lines.append("Direct configuration=     " + str(si + 1))
-            for (i, site) in enumerate(structure):
+            for i, site in enumerate(structure):
                 coords = site.frac_coords
                 line = " ".join([format_str.format(c) for c in coords])
                 line += " " + site.species_string
@@ -72,9 +72,9 @@ def get_string(self, system="unknown system", significant_figures=6):
         #  positions = np.add(self.trajectory[0].frac_coords, self.trajectory.displacements)
         atoms = [site.specie.symbol for site in self.trajectory[0]]
 
-        for (si, position_array) in enumerate(positions):
+        for si, position_array in enumerate(positions):
             lines.append("Direct configuration=     " + str(si + 1))
-            for (i, coords) in enumerate(position_array):
+            for i, coords in enumerate(position_array):
                 line = " ".join([format_str.format(c) for c in coords])
                 line += " " + atoms[i]
                 lines.append(line)
diff --git a/src/mpmorph/jobs/equilibrate_volume.py b/src/mpmorph/jobs/equilibrate_volume.py
index 646fe719..74533c35 100644
--- a/src/mpmorph/jobs/equilibrate_volume.py
+++ b/src/mpmorph/jobs/equilibrate_volume.py
@@ -29,7 +29,6 @@ class EquilibriumVolumeSearchMaker(Maker):
     def make(
         self, original_structure: Structure, md_pv_data_docs: List[MDPVDataDoc] = None
     ):
-
         if md_pv_data_docs is not None and len(md_pv_data_docs) > MAX_MD_JOBS:
             raise RuntimeError(
                 "Maximum number of jobs for equilibrium volume search exceeded"
diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps_volume.py
index 2598ef09..6efab85c 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps_volume.py
@@ -16,6 +16,7 @@
 
 from pkg_resources import resource_filename
 
+
 class LammpsCalcMaker(Maker):
     """
     Run LAMMPS directly using m3gnet (no custodian).
@@ -27,10 +28,7 @@ class LammpsCalcMaker(Maker):
     name = "LAMMPS_CALCULATION"
 
     @job(trajectory="trajectory", output_schema=LammpsCalc)
-    def make(self, temperature: int,
-                   total_steps: int,
-                   structure: Structure = None):
-
+    def make(self, temperature: int, total_steps: int, structure: Structure = None):
         lammps_bin = os.environ.get("LAMMPS_CMD")
         m3gnet_path = os.environ.get("M3GNET_PATH")
 
@@ -40,19 +38,22 @@ def make(self, temperature: int,
             "m3gnet_path": m3gnet_path,
             "species": chem_sys_str,
             "total_steps": total_steps,
-            "print_every_n_step": 10
+            "print_every_n_step": 10,
         }
-        
 
-        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
+        template_path = resource_filename(
+            "mpmorph", "jobs/lammps-templates/template.lammps"
+        )
 
         data_filename: str = "data.lammps"
-        data = LammpsData.from_structure(structure, atom_style='atomic')
+        data = LammpsData.from_structure(structure, atom_style="atomic")
         # Write the input files
-        linp = LammpsTemplateGen().get_input_set(script_template=template_path,
-                                                 settings=script_options, 
-                                                 data=data,
-                                                 data_filename=data_filename)
+        linp = LammpsTemplateGen().get_input_set(
+            script_template=template_path,
+            settings=script_options,
+            data=data,
+            data_filename=data_filename,
+        )
 
         linp.write_input(directory=".")
         input_name = "in.lammps"
@@ -76,12 +77,15 @@ def make(self, temperature: int,
 
         trajectory = Trajectory.from_structures(structs, constant_lattice=False)
 
-        df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", index_col="step", skiprows=1, names=["step", "temp", "vol", "density"])
+        df = pd.read_csv(
+            "step_temp_vol_density.txt",
+            delimiter=" ",
+            index_col="step",
+            skiprows=1,
+            names=["step", "temp", "vol", "density"],
+        )
 
-        metadata = {
-            "temperature": temperature,
-            "total_steps": total_steps
-        }
+        metadata = {"temperature": temperature, "total_steps": total_steps}
 
         output = LammpsCalc(
             dir_name=os.getcwd(),
@@ -89,6 +93,6 @@ def make(self, temperature: int,
             composition=structure.composition,
             reduced_formula=structure.composition.reduced_formula,
             metadata=metadata,
-            dump_data=df.to_dict()
+            dump_data=df.to_dict(),
         )
-        return output
\ No newline at end of file
+        return output
diff --git a/src/mpmorph/jobs/pv_from_calc.py b/src/mpmorph/jobs/pv_from_calc.py
index 57f7bb92..e82ec028 100644
--- a/src/mpmorph/jobs/pv_from_calc.py
+++ b/src/mpmorph/jobs/pv_from_calc.py
@@ -30,7 +30,6 @@ def make(self, structure, scale_factor=None):
 
 @dataclass
 class PVFromM3GNet(PVFromCalc):
-
     name: str = "PV_FROM_M3GNET"
     parameters: M3GNetMDInputs = None
 
@@ -44,10 +43,10 @@ def build_doc(self, m3gnet_calc: M3GNetMDCalculation):
         p_data = m3gnet_calc_to_pressure(m3gnet_calc)
         return MDPVDataDoc(volume=v_data, pressure=p_data)
 
+
 @dataclass
 class PVFromM3GNetLammps(PVFromCalc):
-    """Generates a MDPVDataDoc using Lammps run with M3gnet and a npt ensemble.
-    """
+    """Generates a MDPVDataDoc using Lammps run with M3gnet and a npt ensemble."""
 
     name: str = "PV_FROM_M3GNET_LAMMPS"
     parameters: M3GNetMDInputs = None
@@ -76,7 +75,6 @@ def m3gnet_calc_to_pressure(m3gnet_calc: M3GNetMDCalculation):
 
 @dataclass
 class PVFromVasp(PVFromCalc):
-
     name: str = "PV_FROM_VASP"
     md_maker: Maker = MDMaker()
 
diff --git a/src/mpmorph/jobs/tasks/m3gnet_input.py b/src/mpmorph/jobs/tasks/m3gnet_input.py
index cf47e6c9..8392eceb 100644
--- a/src/mpmorph/jobs/tasks/m3gnet_input.py
+++ b/src/mpmorph/jobs/tasks/m3gnet_input.py
@@ -11,7 +11,6 @@ def one_atmosphere():
 
 @dataclass
 class M3GNetMDInputs:
-
     ensemble: str = "nvt"
     temperature: float = 2000.0
     pressure: float = 1.01325 * units.bar
diff --git a/src/mpmorph/jobs/volume_temperature_sweep.py b/src/mpmorph/jobs/volume_temperature_sweep.py
index dd8aa53c..83c66693 100644
--- a/src/mpmorph/jobs/volume_temperature_sweep.py
+++ b/src/mpmorph/jobs/volume_temperature_sweep.py
@@ -2,8 +2,8 @@
 from mpmorph.jobs.tasks.m3gnet_input import M3GNetMDInputs
 import dataclasses
 
-class VolumeTemperatureSweepMaker(Maker):
 
+class VolumeTemperatureSweepMaker(Maker):
     name: str = "VOLUME_TEMPERATURE_SWEEP"
     md_parameters: M3GNetMDInputs = None
 
@@ -35,5 +35,3 @@ def make(
 
         new_flow = Flow([*volume_jobs, collect_job], output=collect_job.output)
         return Response(replace=new_flow)
-
-
diff --git a/src/mpmorph/schemas/lammps_calc.py b/src/mpmorph/schemas/lammps_calc.py
index 918086e3..c76c17d1 100644
--- a/src/mpmorph/schemas/lammps_calc.py
+++ b/src/mpmorph/schemas/lammps_calc.py
@@ -22,8 +22,7 @@ class LammpsCalc(BaseModel):
         description="The reduced formula of the structure's composition."
     )
     dump_data: dict = Field(
-        None,
-        description="Any additional data collected via LAMMPS dump files"
+        None, description="Any additional data collected via LAMMPS dump files"
     )
     metadata: dict = Field(
         None,
diff --git a/src/mpmorph/schemas/m3gnet_md_calc.py b/src/mpmorph/schemas/m3gnet_md_calc.py
index 941639a0..586e23c6 100644
--- a/src/mpmorph/schemas/m3gnet_md_calc.py
+++ b/src/mpmorph/schemas/m3gnet_md_calc.py
@@ -47,7 +47,6 @@ def from_directory(
         ),
         **kwargs,
     ):
-
         """
         Create a M3GnetCalculation document from a directory containing output files of
         a M3GNet MD run.
diff --git a/src/mpmorph/schemas/pv_data_doc.py b/src/mpmorph/schemas/pv_data_doc.py
index bb40b89b..b836d8fd 100644
--- a/src/mpmorph/schemas/pv_data_doc.py
+++ b/src/mpmorph/schemas/pv_data_doc.py
@@ -2,7 +2,6 @@
 
 
 class MDPVDataDoc(BaseModel):
-
     task_label: str = Field(None, description="The name of the task.")
     volume: float = Field(None, description="The volume data from the MD run")
     pressure: float = Field(None, description="The pressure of the MD run")
diff --git a/src/mpmorph/schemas/vt_sweep_doc.py b/src/mpmorph/schemas/vt_sweep_doc.py
index 4681ba43..c3b7d61f 100644
--- a/src/mpmorph/schemas/vt_sweep_doc.py
+++ b/src/mpmorph/schemas/vt_sweep_doc.py
@@ -3,8 +3,11 @@
 
 
 class VTSweepDoc(BaseModel):
-
     task_label: str = Field(None, description="The name of the task.")
     volumes: float = Field(description="The volume at each temperature")
-    temps: float = Field(description="The temperatures at which the volume was equilibrated")
-    structure: Structure = Field(description="The original structure for which this sweep was performed")
+    temps: float = Field(
+        description="The temperatures at which the volume was equilibrated"
+    )
+    structure: Structure = Field(
+        description="The original structure for which this sweep was performed"
+    )
diff --git a/src/mpmorph/workflows/quench.py b/src/mpmorph/workflows/quench.py
index 21832e2c..6e8066af 100644
--- a/src/mpmorph/workflows/quench.py
+++ b/src/mpmorph/workflows/quench.py
@@ -33,7 +33,7 @@ def get_quench_wf(
     hold_args = kwargs.get("hold_args", {"md_params": {"nsteps": 500}})
     quench_args = kwargs.get("quench_args", {})
 
-    for (i, structure) in enumerate(structures):
+    for i, structure in enumerate(structures):
         _fw_list = []
         if quench_type == "slow_quench":
             for temp in np.arange(

From f68009e9340dba734cc8e4dcb4d6f98144a761ae Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Mon, 13 Feb 2023 19:44:07 -0800
Subject: [PATCH 27/46] remove the unfinished function

---
 src/mpmorph/flows/md_flow.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/flows/md_flow.py b/src/mpmorph/flows/md_flow.py
index b63aa6c6..376fbf89 100644
--- a/src/mpmorph/flows/md_flow.py
+++ b/src/mpmorph/flows/md_flow.py
@@ -15,7 +15,7 @@
 LAMMPS_VOL_FLOW = "LAMMPS_VOL_FLOW"
 
 
-def get_md_temperature_sweeping(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
+# def get_md_temperature_sweeping(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):
 
 
 def get_md_flow_m3gnet(structure, temp, steps, converge_first = True, initial_vol_scale = 1, **input_kwargs):

From 268601a2ce8e9b3ab2e4128ec5b10e95e2b5b659 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Wed, 22 Feb 2023 10:35:27 -0800
Subject: [PATCH 28/46] Add basic temp sweep code

---
 src/mpmorph/jobs/lammps/helpers.py            | 37 +++++++++++
 .../lammps_basic_const_temp.py}               | 45 ++-----------
 .../jobs/lammps/lammps_basic_temp_sweep.py    | 65 +++++++++++++++++++
 .../templates/basic_constant_temp.lammps}     |  0
 .../templates/basic_temp_sweep.lammps}        | 25 ++++++-
 5 files changed, 131 insertions(+), 41 deletions(-)
 create mode 100644 src/mpmorph/jobs/lammps/helpers.py
 rename src/mpmorph/jobs/{lammps_volume.py => lammps/lammps_basic_const_temp.py} (50%)
 create mode 100644 src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
 rename src/mpmorph/jobs/{lammps-templates/template.lammps => lammps/templates/basic_constant_temp.lammps} (100%)
 rename src/mpmorph/jobs/{lammps-templates/in.melt_template => lammps/templates/basic_temp_sweep.lammps} (77%)

diff --git a/src/mpmorph/jobs/lammps/helpers.py b/src/mpmorph/jobs/lammps/helpers.py
new file mode 100644
index 00000000..b5fa38ab
--- /dev/null
+++ b/src/mpmorph/jobs/lammps/helpers.py
@@ -0,0 +1,37 @@
+from ase.io.lammpsrun import read_lammps_dump_text
+from pymatgen.io.ase import AseAtomsAdaptor
+from pymatgen.core.trajectory import Trajectory
+from pymatgen.io.lammps.inputs import LammpsTemplateGen
+from pymatgen.io.lammps.data import LammpsData
+from subprocess import PIPE, Popen
+
+def trajectory_from_lammps_dump(dump_path):
+    with open(dump_path, "r+") as f:
+        atoms = read_lammps_dump_text(f, index=slice(0, None))
+
+    structs = []
+
+    for a in atoms:
+        structs.append(AseAtomsAdaptor().get_structure(a))
+
+    return Trajectory.from_structures(structs, constant_lattice=False)
+
+def run_lammps(structure, template_path, template_opts, lammps_bin):
+    data_filename: str = "data.lammps"
+    data = LammpsData.from_structure(structure, atom_style='atomic')
+    # Write the input files
+    linp = LammpsTemplateGen().get_input_set(script_template=template_path,
+                                                settings=template_opts, 
+                                                data=data,
+                                                data_filename=data_filename)
+
+    linp.write_input(directory=".")
+    input_name = "in.lammps"
+    # Run LAMMPS
+
+    lammps_cmd = [lammps_bin, "-in", input_name]
+    print(f"Running: {' '.join(lammps_cmd)}")
+    with Popen(lammps_cmd, stdout=PIPE, stderr=PIPE) as p:
+        (stdout, stderr) = p.communicate()
+
+    print(f"LAMMPS finished running: {stdout} \n {stderr}")
\ No newline at end of file
diff --git a/src/mpmorph/jobs/lammps_volume.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
similarity index 50%
rename from src/mpmorph/jobs/lammps_volume.py
rename to src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index 2598ef09..f5bfdecc 100644
--- a/src/mpmorph/jobs/lammps_volume.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -1,24 +1,17 @@
-from subprocess import PIPE, Popen
-
 import os
 from jobflow import Maker, job
 
 import pandas as pd
-
-from pymatgen.io.lammps.inputs import LammpsTemplateGen
-from pymatgen.io.lammps.data import LammpsData
 from pymatgen.core.structure import Structure
-from ase.io.lammpsrun import read_lammps_dump_text
-from pymatgen.io.ase import AseAtomsAdaptor
-from pymatgen.core.trajectory import Trajectory
+from .helpers import run_lammps, trajectory_from_lammps_dump
 
 from mpmorph.schemas.lammps_calc import LammpsCalc
 
 from pkg_resources import resource_filename
 
-class LammpsCalcMaker(Maker):
+class BasicLammpsConstantTempMaker(Maker):
     """
-    Run LAMMPS directly using m3gnet (no custodian).
+    Run LAMMPS directly using m3gnet at a constant temperature.
     Required params:
         lammsps_cmd (str): lammps command to run sans the input file name.
             e.g. 'mpirun -n 4 lmp_mpi'
@@ -44,37 +37,11 @@ def make(self, temperature: int,
         }
         
 
-        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/template.lammps')
-
-        data_filename: str = "data.lammps"
-        data = LammpsData.from_structure(structure, atom_style='atomic')
-        # Write the input files
-        linp = LammpsTemplateGen().get_input_set(script_template=template_path,
-                                                 settings=script_options, 
-                                                 data=data,
-                                                 data_filename=data_filename)
-
-        linp.write_input(directory=".")
-        input_name = "in.lammps"
-        # Run LAMMPS
-
-        lammps_cmd = [lammps_bin, "-in", input_name]
-        print(f"Running: {' '.join(lammps_cmd)}")
-        with Popen(lammps_cmd, stdout=PIPE, stderr=PIPE) as p:
-            (stdout, stderr) = p.communicate()
-
-        print(f"LAMMPS finished running: {stdout} \n {stderr}")
-
-        # Build trajectory from LAMMPS output .xyz file
-        with open("trajectory.lammpstrj", "r+") as f:
-            atoms = read_lammps_dump_text(f, index=slice(0, None))
-
-        structs = []
+        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/basic_constant_temp.lammps')
 
-        for a in atoms:
-            structs.append(AseAtomsAdaptor().get_structure(a))
+        run_lammps(structure, template_path, script_options, lammps_bin)
 
-        trajectory = Trajectory.from_structures(structs, constant_lattice=False)
+        trajectory = trajectory_from_lammps_dump("trajectory.lammpstrj")
 
         df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", index_col="step", skiprows=1, names=["step", "temp", "vol", "density"])
 
diff --git a/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
new file mode 100644
index 00000000..8135d48f
--- /dev/null
+++ b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
@@ -0,0 +1,65 @@
+import os
+from jobflow import Maker, job
+
+import pandas as pd
+from pymatgen.core.structure import Structure
+
+from mpmorph.schemas.lammps_calc import LammpsCalc
+from .helpers import run_lammps, trajectory_from_lammps_dump
+
+from pkg_resources import resource_filename
+
+class BasicLammpsTempSweepMaker(Maker):
+    """
+    Run LAMMPS directly using m3gnet sweeping over a range of temperatures.
+    Required params:
+        lammsps_cmd (str): lammps command to run sans the input file name.
+            e.g. 'mpirun -n 4 lmp_mpi'
+    """
+
+    name = "LAMMPS_CALCULATION"
+
+    @job(trajectory="trajectory", output_schema=LammpsCalc)
+    def make(self, temp_initial: int,
+                   temp_final: int,
+                   total_steps: int,
+                   structure: Structure = None):
+
+        lammps_bin = os.environ.get("LAMMPS_CMD")
+        m3gnet_path = os.environ.get("M3GNET_PATH")
+
+        chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
+
+        script_options = {
+            "tempstart": temp_initial,
+            "tempstop": temp_final,
+            "m3gnet_path": m3gnet_path,
+            "species": chem_sys_str,
+            "total_steps": total_steps,
+            "print_every_n_step": 10
+        }
+        
+
+        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/basic_temp_sweep.lammps')
+
+        run_lammps(structure, template_path, script_options, lammps_bin)
+
+        trajectory = trajectory_from_lammps_dump("trajectory.lammpstrj")
+
+        # df = pd.read_csv("step_temp_vol_density.txt", delimiter=" ", index_col="step", skiprows=1, names=["step", "temp", "vol", "density"])
+
+        metadata = {
+            "temp_initial": temp_initial,
+            "temp_final": temp_final,
+            "total_steps": total_steps
+        }
+
+        output = LammpsCalc(
+            dir_name=os.getcwd(),
+            trajectory=trajectory,
+            composition=structure.composition,
+            reduced_formula=structure.composition.reduced_formula,
+            metadata=metadata,
+            dump_data={}
+        )
+        return output
\ No newline at end of file
diff --git a/src/mpmorph/jobs/lammps-templates/template.lammps b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
similarity index 100%
rename from src/mpmorph/jobs/lammps-templates/template.lammps
rename to src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
diff --git a/src/mpmorph/jobs/lammps-templates/in.melt_template b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
similarity index 77%
rename from src/mpmorph/jobs/lammps-templates/in.melt_template
rename to src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index 0381b7c3..8629693c 100644
--- a/src/mpmorph/jobs/lammps-templates/in.melt_template
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -1,5 +1,23 @@
 # Script originally made by Oscar Guerrero
 # Reference: https://orca.cardiff.ac.uk/id/eprint/101322/1/MRSPaper2.pdf
+
+#####################################################################
+# READ THIS BEFORE MAKING CHANGES
+#
+# NOTICE: ANY VARIABLES BEGINNING WITH _ ARE INTERNAL TO THE LAMMPS SCRIPT
+#         VARIABLES FOR USE REPLACEMENT BY TEMPLATE HAVE NO _ IN NAME
+#
+# Parameters for this template:
+#   tempstart - The starting temp for the simulation
+#   tempstop - The ending temp for the simulation
+#   species - A space-separated list of elements present in the system, e.g. "Y Mn O"
+#   m3gnet_path - The path to the m3gnet potential installation
+#   print_every_n_step - The frequency with which info should be printed to stdout
+#   total_steps - The total number of simulation steps
+##########################################################################
+
+
+
 units  metal
 atom_style atomic
 boundary p p p
@@ -27,7 +45,7 @@ read_data data.dump
 
 #Define Interatomic Potential
 
-pair_style    m3gnet /global/home/users/huizheng/repos/lammps/potentials/M3GNET
+pair_style    m3gnet $m3gnet_path
 pair_coeff    * *  MP-2021.2.8-EFS   $species
 
 # Equilibration
@@ -82,8 +100,11 @@ run 0
 fix data_melting all print $print_every_n_step "${eq8} ${eq9}" file temp_vs_ref_vol.txt screen no
 
 dump 1 all cfg 100 HgF2.step*.cfg mass type xs ys zs id
+# What does this do?
 dump_modify 1 element  $species
-dump 2 all custom  100  dump.* id type x y z
+
+dump 2 all custom 100 trajectory.lammpstrj id element x y z
+
 # Compute msd command and dump every 10 steps
 compute msd all msd com yes
 fix msd all ave/time 1 1 10 c_msd[4] file msd.txt

From deed20e41ac9efc1453df4ce5f05e73e2dfea5b9 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 11:41:25 -0800
Subject: [PATCH 29/46] update the directory name to templates

---
 src/mpmorph/jobs/core.py                           | 1 +
 src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/core.py b/src/mpmorph/jobs/core.py
index 50ef7524..ecb1c828 100644
--- a/src/mpmorph/jobs/core.py
+++ b/src/mpmorph/jobs/core.py
@@ -38,5 +38,6 @@ def make(self, structure: Structure, **kwargs):
         """
 
         calc_doc = run_m3gnet(structure, self.parameters, self.name, **kwargs)
+        
 
         return calc_doc
diff --git a/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
index 8135d48f..37a4e0f5 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
@@ -40,7 +40,7 @@ def make(self, temp_initial: int,
         }
         
 
-        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/basic_temp_sweep.lammps')
+        template_path = resource_filename('mpmorph', 'jobs/templates/basic_temp_sweep.lammps')
 
         run_lammps(structure, template_path, script_options, lammps_bin)
 

From cdc670c6204ec4d00ac5f1ad2f5149e8b5b4c955 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 12:13:22 -0800
Subject: [PATCH 30/46] fix the template path

---
 src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
index 37a4e0f5..943fd6a8 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_temp_sweep.py
@@ -40,7 +40,7 @@ def make(self, temp_initial: int,
         }
         
 
-        template_path = resource_filename('mpmorph', 'jobs/templates/basic_temp_sweep.lammps')
+        template_path = resource_filename('mpmorph', 'jobs/lammps/templates/basic_temp_sweep.lammps')
 
         run_lammps(structure, template_path, script_options, lammps_bin)
 

From a9b25e7bb94362a309e59a0475ca996807678181 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 12:22:14 -0800
Subject: [PATCH 31/46] rename data.lammps to data.dump to match with template

---
 src/mpmorph/jobs/lammps/helpers.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/helpers.py b/src/mpmorph/jobs/lammps/helpers.py
index b5fa38ab..57f91845 100644
--- a/src/mpmorph/jobs/lammps/helpers.py
+++ b/src/mpmorph/jobs/lammps/helpers.py
@@ -17,7 +17,7 @@ def trajectory_from_lammps_dump(dump_path):
     return Trajectory.from_structures(structs, constant_lattice=False)
 
 def run_lammps(structure, template_path, template_opts, lammps_bin):
-    data_filename: str = "data.lammps"
+    data_filename: str = "data.dump"
     data = LammpsData.from_structure(structure, atom_style='atomic')
     # Write the input files
     linp = LammpsTemplateGen().get_input_set(script_template=template_path,

From a91d4db3e71d026eef04f92dc92a053716258610 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 14:45:23 -0800
Subject: [PATCH 32/46] update the lammps template file

---
 .../lammps/templates/basic_temp_sweep.lammps  | 33 +++++++++++++++----
 1 file changed, 27 insertions(+), 6 deletions(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index 8629693c..91915927 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -39,7 +39,7 @@ variable pdamp equal "v_time_step*1000" # DO NOT CHANGE
 timestep ${time_step} # DO NOT CHANGE
 
 
-#Create structure
+# Create structure
 read_data data.dump
 
 
@@ -99,16 +99,37 @@ variable eq10 equal "etotal"
 run 0 
 fix data_melting all print $print_every_n_step "${eq8} ${eq9}" file temp_vs_ref_vol.txt screen no
 
-dump 1 all cfg 100 HgF2.step*.cfg mass type xs ys zs id
-# What does this do?
+dump 1 all cfg 100 step*.cfg mass type xs ys zs id
+# What does this do? so the element str can be written to the cfg file
 dump_modify 1 element  $species
-
-dump 2 all custom 100 trajectory.lammpstrj id element x y z
+dump 2 all custom 100 trajectory.lammpstrj id type x y z
+# dump all the dump files into one file
+# wildcard * is used to dump all the snapshots into invividual files
+dump 2 all custom  100  dump.* id type x y z
 
 # Compute msd command and dump every 10 steps
 compute msd all msd com yes
 fix msd all ave/time 1 1 10 c_msd[4] file msd.txt
 
+group lithium type 1
+# group lanthanum type 2
+# group zirconium type 3
+# group oxygen type 4
+
+compute           mymsd1 lithium msd com yes   
+# compute ID(mymsd) group-ID(lithium) msd keyword(com)values(yes: If the com option is set to yes
+#  then the effect of any drift in the center-of-mass of the group of atoms is subtracted out 
+# before the displacement of each atom is calculated.)
+variable          msdxLi equal "c_mymsd1[1]"
+variable          msdyLi equal "c_mymsd1[2]"
+variable          msdzLi equal "c_mymsd1[3]"
+variable          msdtotLi equal "c_mymsd1[4]"
+fix               msdT1 lithium ave/time 1 1 1000 v_msdxLi v_msdyLi v_msdzLi v_msdtotLi file msd_Li.dat 
+#fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ... 
+# For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on timesteps 90,92,94,96,98,100 will be used to compute the final average on timestep 100. 
+# Similarly for timesteps 190,192,194,196,198,200 on timestep 200, etc.
+# If Nrepeat=1 and Nfreq = 100, then no time averaging is done; values are simply generated on timesteps 100,200,etc. (v_msdx v_msdy v_msdz v_msdtot)
+
 
 # use velocity auto-correlation function (VACF) to calculate diffusion coefficient
 compute         2 all vacf
@@ -134,6 +155,6 @@ write_data data.*
 # SAVE THE DATA OF THE CALCULATION OR ELSE YOU NEED TO START OVER = ( OUCH !
 # SIMULATION DONE
 clear
-print "You've done great job =)"
+print "Simulation done! You have done great job!"
 
 

From 2cc395140140fcb29c331138707da9283c6f6b0e Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 15:09:09 -0800
Subject: [PATCH 33/46] fix template file

---
 src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index 91915927..eedce2fc 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -105,7 +105,7 @@ dump_modify 1 element  $species
 dump 2 all custom 100 trajectory.lammpstrj id type x y z
 # dump all the dump files into one file
 # wildcard * is used to dump all the snapshots into invividual files
-dump 2 all custom  100  dump.* id type x y z
+dump 3 all custom  100  dump.* id type x y z
 
 # Compute msd command and dump every 10 steps
 compute msd all msd com yes
@@ -150,6 +150,7 @@ unfix melting
 unfix data_melting 
 undump 1
 undump 2
+undump 3
 write_restart restart.equil
 write_data data.*
 # SAVE THE DATA OF THE CALCULATION OR ELSE YOU NEED TO START OVER = ( OUCH !

From c7755cd96b5a0790cd27baffd130710ec33c9183 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 15:45:32 -0800
Subject: [PATCH 34/46] fix path of template file

---
 src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
index ab2a161a..dec51563 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
@@ -14,7 +14,7 @@ units         metal
 boundary      p p p
 atom_style    atomic
 
-read_data     data.lammps 
+read_data     data.dump 
 
 pair_style    m3gnet $m3gnet_path
 pair_coeff    * *  MP-2021.2.8-EFS  $species

From 4b4060f78a70cf1a8fae93509722543c79b55c97 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 15:49:04 -0800
Subject: [PATCH 35/46] fix template path

---
 src/mpmorph/jobs/lammps/lammps_basic_const_temp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index 47efd8ea..6ca7d7e4 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -33,7 +33,7 @@ def make(self, temperature: int, total_steps: int, structure: Structure = None):
             "print_every_n_step": 10,
         }
 
-        template_path = resource_filename('mpmorph', 'jobs/lammps-templates/basic_constant_temp.lammps')
+        template_path = resource_filename('mpmorph', 'jobs/lammps/templates/basic_constant_temp.lammps')
 
         run_lammps(structure, template_path, script_options, lammps_bin)
 

From 737b9f0109d4124f492a8133c66356eca42f96d5 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 22 Feb 2023 16:27:58 -0800
Subject: [PATCH 36/46] dump ele str into trajectory.lammpstrj

---
 src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index eedce2fc..b1622557 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -102,8 +102,11 @@ fix data_melting all print $print_every_n_step "${eq8} ${eq9}" file temp_vs_ref_
 dump 1 all cfg 100 step*.cfg mass type xs ys zs id
 # What does this do? so the element str can be written to the cfg file
 dump_modify 1 element  $species
-dump 2 all custom 100 trajectory.lammpstrj id type x y z
-# dump all the dump files into one file
+
+dump 2 all custom 100 trajectory.lammpstrj id element x y z
+dump_modify   2 sort id element $species
+# dump all the dump files into one file trajectory.lammpstrj 
+
 # wildcard * is used to dump all the snapshots into invividual files
 dump 3 all custom  100  dump.* id type x y z
 

From 0ab76f0ec848c09d916221c1bc65fc2e967db03e Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Sat, 25 Feb 2023 11:20:02 -0800
Subject: [PATCH 37/46] add ensemble option for constant temperature run

---
 src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
index dec51563..cb64eee0 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
@@ -34,6 +34,6 @@ fix def1 all print $print_every_n_step "${p1} ${p2} ${p3} ${p4}" file step_temp_
 
 
 velocity      all create $temperature 12345
-fix           myEnse all npt temp $temperature $temperature 0.1 aniso 1.0 1.0 1.0
+fix           myEnse all $ensemble temp $temperature $temperature 0.1 aniso 1.0 1.0 1.0
 timestep      0.002
 run           $total_steps 

From 41619a5c321d70f5d2ba16b94ffa333c6a41df3b Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Sat, 25 Feb 2023 11:25:10 -0800
Subject: [PATCH 38/46] update ensemble option

---
 src/mpmorph/jobs/lammps/lammps_basic_const_temp.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index 6ca7d7e4..94595a38 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -20,12 +20,13 @@ class BasicLammpsConstantTempMaker(Maker):
     name = "LAMMPS_CALCULATION"
 
     @job(trajectory="trajectory", output_schema=LammpsCalc)
-    def make(self, temperature: int, total_steps: int, structure: Structure = None):
+    def make(self, temperature: int, ensemble:int, total_steps: int, structure: Structure = None):
         lammps_bin = os.environ.get("LAMMPS_CMD")
         m3gnet_path = os.environ.get("M3GNET_PATH")
 
         chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
         script_options = {
+            "ensemble": ensemble,
             "temperature": temperature,
             "m3gnet_path": m3gnet_path,
             "species": chem_sys_str,

From b92aac4aa0e4326f491d28207c634efeb72bf72c Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Sat, 25 Feb 2023 11:32:28 -0800
Subject: [PATCH 39/46] minor fix ensemble

---
 src/mpmorph/jobs/lammps/lammps_basic_const_temp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index 94595a38..9e4fafd6 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -20,7 +20,7 @@ class BasicLammpsConstantTempMaker(Maker):
     name = "LAMMPS_CALCULATION"
 
     @job(trajectory="trajectory", output_schema=LammpsCalc)
-    def make(self, temperature: int, ensemble:int, total_steps: int, structure: Structure = None):
+    def make(self, temperature: int, ensemble:str, total_steps: int, structure: Structure = None):
         lammps_bin = os.environ.get("LAMMPS_CMD")
         m3gnet_path = os.environ.get("M3GNET_PATH")
 

From 91ae04491d6b1e59887af0184a8c85e59c24447d Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Tue, 28 Feb 2023 13:07:52 -0800
Subject: [PATCH 40/46] Add ensemble and polyfit estimators

---
 src/mpmorph/analysis/melting_points.py | 197 ++++++++++++++++++-------
 1 file changed, 144 insertions(+), 53 deletions(-)

diff --git a/src/mpmorph/analysis/melting_points.py b/src/mpmorph/analysis/melting_points.py
index a3ff8724..93a8f647 100644
--- a/src/mpmorph/analysis/melting_points.py
+++ b/src/mpmorph/analysis/melting_points.py
@@ -6,31 +6,84 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import math
+from abc import ABC, abstractmethod
+import numpy.polynomial.polynomial as poly
+from numpy.polynomial import Polynomial as P
 
+class AbstractMeltingPointEstimator(ABC):
 
-class MeltingPointClusterAnalyzer:
-    def _get_clusters(self, points):
-        clustering = AgglomerativeClustering(n_clusters=2).fit(points)
-        cluster1 = points[np.argwhere(clustering.labels_ == 1).squeeze()].T
-        cluster2 = points[np.argwhere(clustering.labels_ == 0).squeeze()].T
-        return cluster1, cluster2
+    def plot(self, ts, vs, plot_title=None):
+        fig, axs = self._plot_ts_vs(ts, vs)
+        Tm = self.estimate(ts, vs)
+        axs.plot([Tm, Tm], [min(vs), max(vs)], color="r")
 
-    def plot_vol_vs_temp(self, ts, vs, plot_title=None):
+        if plot_title is None:
+            axs.set_title("Volume vs Temperature by Polynomial Fit")
+        else:
+            axs.set_title(plot_title)
+        
+        return fig, axs
+
+    @abstractmethod
+    def estimate(self, temps, vols):
+        pass
+
+    def _plot_ts_vs(self, ts, vs):
+        fig, axs = plt.subplots()
+        axs.scatter(ts, vs)
+        axs.set_xlabel("Temperature (K)")
+        axs.set_ylabel("Volume (A^3)")
+        return fig, axs        
+
+class MeltingPointEnsembleEstimator(AbstractMeltingPointEstimator):
+
+    def __init__(self, estimators):
+        self._estimators = estimators
+    
+    def estimate(self, temps, vols):
+        tm_estimates = [e.estimate(temps, vols) for e in self._estimators]
+        return np.mean(tm_estimates)
+
+    def plot(self, ts, vs, plot_title = None):
+        _, axs = self._plot_ts_vs(ts, vs)
+        tm = self.estimate(ts, vs)
+        axs.plot([tm, tm], [np.min(vs), np.max(vs)])
+    
+    def plot_all_estimates(self, ts, vs):
+        _, axs = self._plot_ts_vs(ts, vs)
+        for e in self._estimators:
+            tm = e.estimate(ts, vs)
+            axs.plot([tm, tm], [np.min(vs), np.max(vs)], label=e.name)
+
+        avg = self.estimate(ts, vs)
+        axs.plot([avg, avg], [np.min(vs), np.max(vs)], label="Mean Estimate")
+        axs.set_title("Ensemble Tm Estimates")
+        axs.legend()
+
+
+class MeltingPointClusterEstimator(AbstractMeltingPointEstimator):
+
+    name: str = "Clustering"
+
+    def plot_clusters(self, ts, vs, plot_title=None):
         points = np.array(list(zip(ts, vs)))
         cluster1, cluster2 = self._get_clusters(points)
-        plt.scatter(*cluster1)
-        plt.scatter(*cluster2)
-        plt.xlabel("Temperature (K)")
-        plt.ylabel("Volume (A^3)")
-        Tm = self.estimate_melting_temp(ts, vs)
-        plt.plot([Tm, Tm], [min(vs), max(vs)], color="r")
+        fig, axs = plt.subplots()
+        axs.scatter(*cluster1)
+        axs.scatter(*cluster2)
+        axs.set_xlabel("Temperature (K)")
+        axs.set_ylabel("Volume (A^3)")
+        Tm = self.estimate(ts, vs)
+        axs.plot([Tm, Tm], [min(vs), max(vs)], color="r")
 
         if plot_title is None:
-            plt.title("Volume vs Temperature by Clustering")
+            axs.title("Volume vs Temperature by Clustering")
         else:
-            plt.title(plot_title)
+            axs.title(plot_title)
+        
+        return fig, axs
 
-    def estimate_melting_temp(self, temps, vols):
+    def estimate(self, temps, vols):
         points = np.array(list(zip(temps, vols)))
         cluster1, cluster2 = self._get_clusters(points)
         if min(cluster1[0]) < min(cluster2[0]):
@@ -42,67 +95,79 @@ def estimate_melting_temp(self, temps, vols):
 
         return np.mean([max(solid_range), min(liquid_range)])
 
+    def _get_clusters(self, points):
+        clustering = AgglomerativeClustering(n_clusters=2).fit(points)
+        cluster1 = points[np.argwhere(clustering.labels_ == 1).squeeze()].T
+        cluster2 = points[np.argwhere(clustering.labels_ == 0).squeeze()].T
+        return cluster1, cluster2
+
+
+class MeltingPointSlopeEstimator(AbstractMeltingPointEstimator):
 
-class MeltingPointSlopeAnalyzer:
-    def split_dset(self, pts, split_idx):
+    def plot_best_split(self, temps, vols):
+        split_idx = self.get_best_split(temps, vols)
+        fig, axs = self._plot_split(temps, vols, split_idx)
+        Tm = self.estimate(temps, vols)
+        axs.plot([Tm, Tm], [min(vols), max(vols)], color="r")
+        return fig, axs
+
+    def estimate(self, temps, vols):
+        best_split_idx = self.get_best_split(temps, vols)
+        return np.mean([temps[best_split_idx], temps[best_split_idx - 1]])
+
+    def _split_dset(self, pts, split_idx):
         return pts[0:split_idx], pts[split_idx:]
 
-    def assess_splits(self, xs, ys):
+    def _assess_splits(self, xs, ys):
         dset_size = len(xs)
         buffer = max(round(dset_size / 10), 3)
 
         pt_idxs = list(range(buffer + 1, len(xs) - buffer - 1))
         errs = []
         for idx in pt_idxs:
-            _, _, _, _, total_err = self.get_split_fit(xs, ys, idx)
+            _, _, _, _, total_err = self._get_split_fit(xs, ys, idx)
             errs.append(total_err)
 
         return list(zip(pt_idxs, errs))
 
-    def get_linear_ys(self, m, b, xs):
+    def _get_linear_ys(self, m, b, xs):
         return [m * x + b for x in xs]
 
-    def plot_split(self, xs, ys, split_idx):
-        m1, b1, m2, b2, _ = self.get_split_fit(xs, ys, split_idx)
-        leftxs, rightxs = self.split_dset(xs, split_idx)
-        leftys, rightys = self.split_dset(ys, split_idx)
-
-        left_fit_ys = self.get_linear_ys(m1, b1, leftxs)
+    def _plot_split(self, xs, ys, split_idx):
+        m1, b1, m2, b2, _ = self._get_split_fit(xs, ys, split_idx)
+        leftxs, rightxs = self._split_dset(xs, split_idx)
+        leftys, rightys = self._split_dset(ys, split_idx)
 
-        plt.scatter(leftxs, leftys)
-        plt.plot(leftxs, left_fit_ys)
-        plt.title("Volume vs Temperature (w/ best fits by Slope Method")
-        plt.xlabel("Temperature (K)")
-        plt.ylabel("Equil. Volume (cubic Angstroms)")
+        left_fit_ys = self._get_linear_ys(m1, b1, leftxs)
+        fig, axs = plt.subplots()
+        axs.scatter(leftxs, leftys)
+        axs.plot(leftxs, left_fit_ys)
+        axs.title("Volume vs Temperature (w/ best fits by Slope Method")
+        axs.xlabel("Temperature (K)")
+        axs.ylabel("Equil. Volume (cubic Angstroms)")
 
-        right_fit_ys = self.get_linear_ys(m2, b2, rightxs)
+        right_fit_ys = self._get_linear_ys(m2, b2, rightxs)
 
-        plt.scatter(rightxs, rightys)
-        plt.plot(rightxs, right_fit_ys)
+        axs.scatter(rightxs, rightys)
+        axs.plot(rightxs, right_fit_ys)
+        return fig, axs
 
     def get_best_split(self, xs, ys):
-        split_errs = self.assess_splits(xs, ys)
+        split_errs = self._assess_splits(xs, ys)
         errs = [pt[1] for pt in split_errs]
         idxs = [pt[0] for pt in split_errs]
         best_split_idx = idxs[np.argmin(errs)]
         return best_split_idx
 
-    def plot_vol_vs_temp(self, temps, vols):
-        split_idx = self.get_best_split(temps, vols)
-        self.plot_split(temps, vols, split_idx)
-        Tm = self.estimate_melting_temp(temps, vols)
-        print(Tm)
-        plt.plot([Tm, Tm], [min(vols), max(vols)], color="r")
 
-    def estimate_melting_temp(self, temps, vols):
-        best_split_idx = self.get_best_split(temps, vols)
-        return np.mean([temps[best_split_idx], temps[best_split_idx - 1]])
 
+class MeltingPointSlopeRMSEEstimator(MeltingPointSlopeEstimator):
+
+    name: str = "RMSE Bisection"
 
-class MeltingPointSlopeRMSEAnalyzer(MeltingPointSlopeAnalyzer):
-    def get_split_fit(self, xs, ys, split_idx):
-        leftx, rightx = self.split_dset(xs, split_idx)
-        lefty, righty = self.split_dset(ys, split_idx)
+    def _get_split_fit(self, xs, ys, split_idx):
+        leftx, rightx = self._split_dset(xs, split_idx)
+        lefty, righty = self._split_dset(ys, split_idx)
 
         lslope, lintercept, r_value, p_value, std_err = linregress(leftx, lefty)
         left_y_pred = lintercept + lslope * np.array(leftx)
@@ -117,10 +182,13 @@ def get_split_fit(self, xs, ys, split_idx):
         return lslope, lintercept, rslope, rintercept, combined_err
 
 
-class MeltingPointSlopeStdErrAnalyzer(MeltingPointSlopeAnalyzer):
-    def get_split_fit(self, xs, ys, split_idx):
-        leftx, rightx = self.split_dset(xs, split_idx)
-        lefty, righty = self.split_dset(ys, split_idx)
+class MeltingPointSlopeStdErrEstimator(MeltingPointSlopeEstimator):
+
+    name: str = "StdErr Bisection"
+
+    def _get_split_fit(self, xs, ys, split_idx):
+        leftx, rightx = self._split_dset(xs, split_idx)
+        lefty, righty = self._split_dset(ys, split_idx)
 
         leftfit = linregress(leftx, lefty)
         lefterr = leftfit.stderr
@@ -137,3 +205,26 @@ def get_split_fit(self, xs, ys, split_idx):
             rightfit.intercept,
             combined_err,
         )
+
+class MeltingPointPolyFitEstimator(AbstractMeltingPointEstimator):
+
+    name: str = "Polynomial Fit"
+
+    def plot_fit(self, ts, vs, plot_title=None):
+        fig, axs = super().plot(ts, vs, plot_title=plot_title)
+        p = self._polyfit(ts, vs)
+        fit_ts = p(vs)
+        axs.plot(fit_ts, vs, color='r')
+        return fig, axs
+
+    def estimate(self, temps, vols):
+        p = self._polyfit(temps, vols)
+        second = p.deriv(2)
+        melting_v = second.roots()[0]
+        melting_t = p(melting_v)
+        return melting_t
+    
+    def _polyfit(self, temps, vols):
+        coefs = poly.polyfit(vols, temps, 3)
+        p = P(coefs)
+        return p
\ No newline at end of file

From c2bada76e2e1e0330fe4be4225bbf7609b99197b Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Tue, 28 Feb 2023 13:08:00 -0800
Subject: [PATCH 41/46] fix typo in lammps input template

---
 src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index 8629693c..e65d6a41 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -40,7 +40,7 @@ timestep ${time_step} # DO NOT CHANGE
 
 
 #Create structure
-read_data data.dump
+read_data data.lammps
 
 
 #Define Interatomic Potential

From 606f0615ac80102f0bb27ce380cdd49632da24a5 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 1 Mar 2023 12:47:17 -0800
Subject: [PATCH 42/46] update the contant_temp template, enable NVT/NPT

---
 .../jobs/lammps/lammps_basic_const_temp.py    | 21 ++++++++++++++++---
 .../templates/basic_constant_temp.lammps      |  2 +-
 2 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index 9e4fafd6..fdb88a23 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -20,19 +20,34 @@ class BasicLammpsConstantTempMaker(Maker):
     name = "LAMMPS_CALCULATION"
 
     @job(trajectory="trajectory", output_schema=LammpsCalc)
-    def make(self, temperature: int, ensemble:str, total_steps: int, structure: Structure = None):
+    def make(self, temperature: int, ensemble:str, total_steps: int, structure: Structure = None, 
+    barostat:str =None, Pstart:str = None,Pstop:str = None, Pdamp:str = None):
         lammps_bin = os.environ.get("LAMMPS_CMD")
         m3gnet_path = os.environ.get("M3GNET_PATH")
 
         chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
-        script_options = {
+        if ensemble == "npt":
+            script_options = {
             "ensemble": ensemble,
+            "barostat":barostat,
+            "Pstart":Pstart,
+            "Pstop":Pstop,
+            "Pdamp":Pdamp,
             "temperature": temperature,
             "m3gnet_path": m3gnet_path,
             "species": chem_sys_str,
             "total_steps": total_steps,
             "print_every_n_step": 10,
-        }
+        }            
+        elif ensemble =='nvt':
+            script_options = {
+                "ensemble": ensemble,
+                "temperature": temperature,
+                "m3gnet_path": m3gnet_path,
+                "species": chem_sys_str,
+                "total_steps": total_steps,
+                "print_every_n_step": 10,
+            }
 
         template_path = resource_filename('mpmorph', 'jobs/lammps/templates/basic_constant_temp.lammps')
 
diff --git a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
index cb64eee0..e59fc3d5 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_constant_temp.lammps
@@ -34,6 +34,6 @@ fix def1 all print $print_every_n_step "${p1} ${p2} ${p3} ${p4}" file step_temp_
 
 
 velocity      all create $temperature 12345
-fix           myEnse all $ensemble temp $temperature $temperature 0.1 aniso 1.0 1.0 1.0
+fix           myEnse all $ensemble temp $temperature $temperature 0.1 $barostat $Pstart $Pstop $Pdamp
 timestep      0.002
 run           $total_steps 

From d84dbb0a7e39e071bcbc687f39beb8639d5880b2 Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Wed, 1 Mar 2023 13:39:58 -0800
Subject: [PATCH 43/46] reduce the redundancy

---
 .../jobs/lammps/lammps_basic_const_temp.py    | 34 ++++++++-----------
 1 file changed, 15 insertions(+), 19 deletions(-)

diff --git a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
index fdb88a23..566fc189 100644
--- a/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
+++ b/src/mpmorph/jobs/lammps/lammps_basic_const_temp.py
@@ -21,33 +21,29 @@ class BasicLammpsConstantTempMaker(Maker):
 
     @job(trajectory="trajectory", output_schema=LammpsCalc)
     def make(self, temperature: int, ensemble:str, total_steps: int, structure: Structure = None, 
-    barostat:str =None, Pstart:str = None,Pstop:str = None, Pdamp:str = None):
+    barostat:str =None, Pstart:str = None, Pstop:str = None, Pdamp:str = None):
         lammps_bin = os.environ.get("LAMMPS_CMD")
         m3gnet_path = os.environ.get("M3GNET_PATH")
 
-        chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)
-        if ensemble == "npt":
-            script_options = {
+        chem_sys_str = " ".join(el.symbol for el in structure.composition.elements)          
+        script_options ={"temperature": temperature,
+                        "m3gnet_path": m3gnet_path,
+                        "species": chem_sys_str,
+                        "total_steps": total_steps,
+                        "print_every_n_step": 10}
+
+        if ensemble =='nvt':
+            script_options.update({
+                "ensemble": ensemble
+                })
+        elif ensemble =='npt':
+            script_options.update({ 
             "ensemble": ensemble,
             "barostat":barostat,
             "Pstart":Pstart,
             "Pstop":Pstop,
             "Pdamp":Pdamp,
-            "temperature": temperature,
-            "m3gnet_path": m3gnet_path,
-            "species": chem_sys_str,
-            "total_steps": total_steps,
-            "print_every_n_step": 10,
-        }            
-        elif ensemble =='nvt':
-            script_options = {
-                "ensemble": ensemble,
-                "temperature": temperature,
-                "m3gnet_path": m3gnet_path,
-                "species": chem_sys_str,
-                "total_steps": total_steps,
-                "print_every_n_step": 10,
-            }
+        }  )
 
         template_path = resource_filename('mpmorph', 'jobs/lammps/templates/basic_constant_temp.lammps')
 

From 597261222d8e292440c6d3a629cf9317f099ecbd Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Wed, 1 Mar 2023 15:57:38 -0800
Subject: [PATCH 44/46] A few fixes

---
 src/mpmorph/analysis/melting_points.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/mpmorph/analysis/melting_points.py b/src/mpmorph/analysis/melting_points.py
index 93a8f647..250910d8 100644
--- a/src/mpmorph/analysis/melting_points.py
+++ b/src/mpmorph/analysis/melting_points.py
@@ -56,7 +56,7 @@ def plot_all_estimates(self, ts, vs):
             axs.plot([tm, tm], [np.min(vs), np.max(vs)], label=e.name)
 
         avg = self.estimate(ts, vs)
-        axs.plot([avg, avg], [np.min(vs), np.max(vs)], label="Mean Estimate")
+        axs.plot([avg , avg ], [np.min(vs), np.max(vs)], label="Mean Estimate")
         axs.set_title("Ensemble Tm Estimates")
         axs.legend()
 
@@ -138,18 +138,20 @@ def _plot_split(self, xs, ys, split_idx):
         leftxs, rightxs = self._split_dset(xs, split_idx)
         leftys, rightys = self._split_dset(ys, split_idx)
 
-        left_fit_ys = self._get_linear_ys(m1, b1, leftxs)
         fig, axs = plt.subplots()
+
+        left_fit_ys = self._get_linear_ys(m1, b1, leftxs)
         axs.scatter(leftxs, leftys)
         axs.plot(leftxs, left_fit_ys)
-        axs.title("Volume vs Temperature (w/ best fits by Slope Method")
-        axs.xlabel("Temperature (K)")
-        axs.ylabel("Equil. Volume (cubic Angstroms)")
 
-        right_fit_ys = self._get_linear_ys(m2, b2, rightxs)
 
+        right_fit_ys = self._get_linear_ys(m2, b2, rightxs)
         axs.scatter(rightxs, rightys)
         axs.plot(rightxs, right_fit_ys)
+
+        axs.title("Volume vs Temperature (w/ best fits by Slope Method")
+        axs.xlabel("Temperature (K)")
+        axs.ylabel("Equil. Volume (cubic Angstroms)")
         return fig, axs
 
     def get_best_split(self, xs, ys):
@@ -165,6 +167,12 @@ class MeltingPointSlopeRMSEEstimator(MeltingPointSlopeEstimator):
 
     name: str = "RMSE Bisection"
 
+    def _get_fit_error(self, xs, ys):
+        slope, intercept, r_value, p_value, std_err = linregress(xs, ys)
+        y_pred = intercept + slope * np.array(xs)
+        err = mean_squared_error(y_true=ys, y_pred=y_pred, squared=False)        
+        return slope, intercept, err
+
     def _get_split_fit(self, xs, ys, split_idx):
         leftx, rightx = self._split_dset(xs, split_idx)
         lefty, righty = self._split_dset(ys, split_idx)

From 9b8f31185df539ba796d8a64917ce21d989007f0 Mon Sep 17 00:00:00 2001
From: Max Gallant <mcgallant72@gmail.com>
Date: Wed, 1 Mar 2023 16:09:08 -0800
Subject: [PATCH 45/46] Add trisection method

---
 src/mpmorph/analysis/melting_points.py | 104 ++++++++++++++++++++++++-
 1 file changed, 100 insertions(+), 4 deletions(-)

diff --git a/src/mpmorph/analysis/melting_points.py b/src/mpmorph/analysis/melting_points.py
index 250910d8..30cfcbec 100644
--- a/src/mpmorph/analysis/melting_points.py
+++ b/src/mpmorph/analysis/melting_points.py
@@ -61,6 +61,7 @@ def plot_all_estimates(self, ts, vs):
         axs.legend()
 
 
+
 class MeltingPointClusterEstimator(AbstractMeltingPointEstimator):
 
     name: str = "Clustering"
@@ -101,8 +102,104 @@ def _get_clusters(self, points):
         cluster2 = points[np.argwhere(clustering.labels_ == 0).squeeze()].T
         return cluster1, cluster2
 
+class MeltingPointTrisectionEstimator(AbstractMeltingPointEstimator):
+
+    def _get_fit_error_total(self, xs, ys):
+        slope, intercept, r_value, p_value, std_err = linregress(xs, ys)
+        y_pred = intercept + slope * np.array(xs)
+        err = np.sum(np.abs(y_pred - ys))        
+        return slope, intercept, err
+
+    def _unzip_pts(self, pts):
+        ts = [pt[0] for pt in pts]
+        vs = [pt[1] for pt in pts]
+        return ts, vs
+
+    def _plot_pts(self, pts):
+        ts, vs = self._unzip_pts(pts)
+        plt.scatter(ts,vs, color='grey')
+
+
+    def _split_pts(self, pts, x1, x2):
+        set1 = [pt for pt in pts if pt[0] < x1]
+        set2 = [pt for pt in pts if pt[0] > x1 and pt[0] < x2]
+        set3 = [pt for pt in pts if pt[0] > x2]
+        return set1, set2, set3
+
+    def _plot_split(self, pts, x1, x2):
+        set1, set2, set3 = self._plit_pts(pts, x1, x2)
+        self._plot_pts(set1)
+        self._plot_pts(set2)
+        self._plot_pts(set3)
+
+    def _get_linear_ys(self, m, b, xs):
+        return [m * x + b for x in xs]
+
+    def _plot_fit_line(self, m, b, xs):
+        fit_ys = self._get_linear_ys(m, b, xs)
+        plt.plot(xs, fit_ys)
+
+    def _plot_fits(self, pts, x1, x2):
+        set1, set2, set3 = self._split_pts(pts, x1, x2)
+        for dset in [set1, set2, set3]:
+            xs, ys = self._unzip_pts(dset)
+            m, b, err = self._get_fit_error_total(xs, ys)
+            self._plot_pts(dset)
+            self._plot_fit_line(m, b, xs)
 
-class MeltingPointSlopeEstimator(AbstractMeltingPointEstimator):
+
+    def _find_best_trisection(self, points, min_x = None, max_x = None, min_window_size = 100, step_size = 50):
+        lowest_err = math.inf
+        
+        xs, ys = self._unzip_pts(points)
+        if min_x is None:
+            min_x = math.floor(np.min(xs))
+
+        if max_x is None:
+            max_x = math.ceil(np.max(xs))
+        
+        for pt1 in range(min_x + min_window_size, max_x - 2 * min_window_size, step_size):
+            for pt2 in range(pt1 + min_window_size, max_x - min_window_size, step_size):
+                set1, set2, set3 = self._split_pts(points, pt1, pt2)
+                errs_total = 0
+                errs = []
+                try:
+                    for dset in [set1, set2, set3]:
+                        xs, ys = self._unzip_pts(dset)
+                        m, b, err = self._get_fit_error_total(xs, ys)      
+                        errs.append(err)
+                        errs_total += err ** 2
+
+                    errs_total = math.sqrt(errs_total)
+
+                    if errs_total < lowest_err:
+                        lowest_err = errs_total
+                        best_pt1 = pt1
+                        best_pt2 = pt2
+                except:
+                    print("problem encountered")
+
+        return best_pt1, best_pt2
+
+    def _estimate_melting_pt(self, points):
+        pt1, pt2 = self._find_best_trisection(points, step_size=100)
+        # print(f'Coarse guess: {pt1, pt2}')
+        pt1, pt2 = self._find_best_trisection(points, pt1 - 150, pt2 + 150, step_size = 5)
+        # print(f'Fine guess: {pt1, pt2}')
+        return pt1, pt2, (pt2 - pt1) / 2 + pt1
+    
+    def estimate(self, temps, vols):
+        pts = list(zip(temps, vols))
+        pt1, pt2, tm = self._estimate_melting_pt(pts)
+        return tm
+    
+    def plot(self, ts, vs, plot_title=None):
+        pts = list(zip(ts, vs))
+        pt1, pt2, tm = self._estimate_melting_pt(pts)
+        self._plot_fits(pts, pt1, pt2)
+
+
+class MeltingPointBisectionEstimator(AbstractMeltingPointEstimator):
 
     def plot_best_split(self, temps, vols):
         split_idx = self.get_best_split(temps, vols)
@@ -162,8 +259,7 @@ def get_best_split(self, xs, ys):
         return best_split_idx
 
 
-
-class MeltingPointSlopeRMSEEstimator(MeltingPointSlopeEstimator):
+class MeltingPointSlopeRMSEEstimator(MeltingPointBisectionEstimator):
 
     name: str = "RMSE Bisection"
 
@@ -190,7 +286,7 @@ def _get_split_fit(self, xs, ys, split_idx):
         return lslope, lintercept, rslope, rintercept, combined_err
 
 
-class MeltingPointSlopeStdErrEstimator(MeltingPointSlopeEstimator):
+class MeltingPointSlopeStdErrEstimator(MeltingPointBisectionEstimator):
 
     name: str = "StdErr Bisection"
 

From 96ead3de660c8f96d9c74a0e1b7bb2edef8bb7ac Mon Sep 17 00:00:00 2001
From: Hui Zheng <huz071@eng.ucsd.edu>
Date: Thu, 16 Mar 2023 23:10:29 -0700
Subject: [PATCH 46/46] add box tilt large to lammps template

---
 src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
index b15379bb..4262ee5a 100644
--- a/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
+++ b/src/mpmorph/jobs/lammps/templates/basic_temp_sweep.lammps
@@ -21,6 +21,7 @@
 units  metal
 atom_style atomic
 boundary p p p
+box tilt large
 
 atom_modify map array