From 712c4bd99458c41903f52c727d20d8c211240dfa Mon Sep 17 00:00:00 2001
From: Paul Saxe <paul@thesaxes.net>
Date: Mon, 6 Mar 2023 17:35:21 -0800
Subject: [PATCH 1/4] Added correct molecule numbers, and handle ReaxFF and EAM
 correctly.

---
 lammps_step/lammps.py | 29 ++++++++++++++++++++---------
 1 file changed, 20 insertions(+), 9 deletions(-)

diff --git a/lammps_step/lammps.py b/lammps_step/lammps.py
index e69b2722..be1eb124 100644
--- a/lammps_step/lammps.py
+++ b/lammps_step/lammps.py
@@ -833,15 +833,26 @@ def structure_data(self, eex, triclinic=False):
         lines.append("")
 
         if "charges" in eex:
-            for i, xyz_index, q in zip(
-                range(1, eex["n_atoms"] + 1), eex["atoms"], eex["charges"]
-            ):
-                x, y, z, index = xyz_index
-                # The '1' is molecule ID ... should correct at some point!
-                lines.append(
-                    f"{i:6d}      1 {index:6d} {q:6.3f} {x:12.7f} {y:12.7f} "
-                    f"{z:12.7f}"
-                )
+            if "molecule" in eex:
+                for i, mol, xyz_index, q in zip(
+                    range(1, eex["n_atoms"] + 1),
+                    eex["molecule"],
+                    eex["atoms"],
+                    eex["charges"],
+                ):
+                    x, y, z, index = xyz_index
+                    lines.append(
+                        f"{i:6d} {mol + 1:6d} {index:6d} {q:6.3f} {x:12.7f} {y:12.7f} "
+                        f"{z:12.7f}"
+                    )
+            else:
+                for i, xyz_index, q in zip(
+                    range(1, eex["n_atoms"] + 1), eex["atoms"], eex["charges"]
+                ):
+                    x, y, z, index = xyz_index
+                    lines.append(
+                        f"{i:6d} {index:6d} {q:6.3f} {x:12.7f} {y:12.7f} {z:12.7f}"
+                    )
         else:
             for i, xyz_index in enumerate(eex["atoms"]):
                 x, y, z, index = xyz_index

From 81205163838cd6c56d9092a94dc03100ea2c13c3 Mon Sep 17 00:00:00 2001
From: Paul Saxe <paul@thesaxes.net>
Date: Mon, 6 Mar 2023 17:35:43 -0800
Subject: [PATCH 2/4] Handling for ReaxFF

---
 lammps_step/initialization.py | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/lammps_step/initialization.py b/lammps_step/initialization.py
index b9cc0f6d..77e16c07 100644
--- a/lammps_step/initialization.py
+++ b/lammps_step/initialization.py
@@ -494,7 +494,15 @@ def OpenKIM_input(self):
         lines.append("newton              on")
         lines.append("")
         potential = self.get_variable("_OpenKIM_Potential")
-        lines.append(f"kim_init            {potential} metal " "unit_conversion_mode")
+        if "reax" in potential.lower():
+            lines.append("atom_style          charge")
+            lines.append(
+                f"kim                 init {potential} real unit_conversion_mode"
+            )
+        else:
+            lines.append(
+                f"kim                 init {potential} metal unit_conversion_mode"
+            )
         lines.append("")
         periodicity = configuration.periodicity
         if periodicity == 0:
@@ -511,7 +519,7 @@ def OpenKIM_input(self):
         lines.append("")
         lines.append("read_data           structure.dat")
         lines.append("")
-        lines.append(f'kim_interactions    {" ".join(eex["atom types"])}')
+        lines.append(f'kim                 interactions {" ".join(eex["atom types"])}')
 
         # Set up standard variables
         for variable in thermo_variables:
@@ -558,7 +566,15 @@ def OpenKIM_energy_expression(self):
             x, y, z = xyz
             result.append((x, y, z, index))
 
-        eex["n_atoms"] = len(result)
+        eex["n_atoms"] = n_atoms = len(result)
         eex["n_atom_types"] = len(atom_types)
 
+        # For ReaxFF we need charges
+        potential = self.get_variable("_OpenKIM_Potential")
+        if "reax" in potential.lower():
+            if "charge" in atoms:
+                eex["charges"] = atoms.get_column_data["charges"]
+            else:
+                eex["charges"] = [0.0] * n_atoms
+
         return eex

From c6ec949ce3aa81e2afc76ee01acab21a1f96b43f Mon Sep 17 00:00:00 2001
From: Paul Saxe <paul@thesaxes.net>
Date: Mon, 6 Mar 2023 17:38:48 -0800
Subject: [PATCH 3/4] Added empty icon links to avoid bug in pydata theme.

---
 docs/conf.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/conf.py b/docs/conf.py
index cfdb1eb8..205c20df 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -128,6 +128,7 @@
 html_theme_options = {
     "github_url": "https://github.com/molssi-seamm/lammps_step",
     "twitter_url": "https://twitter.com/MolSSI_NSF",
+    "icon_links": [],
     "logo": {
         "image_light": "SEAMM logo.png",
         "image_dark": "SEAMM Inverted 288x181.png",

From d7afcf508ddeb6cf8ea87efef73945fb784109f7 Mon Sep 17 00:00:00 2001
From: Paul Saxe <paul@thesaxes.net>
Date: Mon, 6 Mar 2023 17:38:58 -0800
Subject: [PATCH 4/4] Preparing for release.

---
 HISTORY.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/HISTORY.rst b/HISTORY.rst
index e8638a94..35931646 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -1,6 +1,10 @@
 =======
 History
 =======
+2023.4.6 -- Better forcefield handling.
+   * Added correct molecule numbers for valence forcefields.
+   * Correctly handle ReaxFF from OpenKim
+   * Updated for some minor changes in OpenKim
 
 2023.2.6 -- Added handling of OPLS-AA forcefield
    * Added handling of the OPLS-AA forcefield