Merge pull request #37 from molssi-seamm/dev

Improvements for forcefield handling

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

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",

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

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