Add Table potential in LAMMPS

pyiron_contrib/atomistics/lammps/potentials.py

@@ -19,6 +19,7 @@
 import warnings
 import itertools
 from functools import wraps
+from pathlib import Path
 
 
 general_doc = """
@@ -249,7 +250,7 @@ def counter(self):
         def pair_style(self):
             """pair_style to be output only in hybrid"""
             if self.is_hybrid:
-                return self._pair_style
+                return [p.split()[0] for p in self._pair_style]
             else:
                 return len(self._pair_style) * [""]
 
@@ -590,7 +591,68 @@ def __init__(self, pair_style, *chemical_elements, cutoff, **kwargs):
             interacting_species=[self._harmonize_species(chemical_elements)],
             pair_coeff=[" ".join([str(cc) for cc in kwargs.values()]) + f" {cutoff}"],
             cutoff=cutoff,
+            **kwargs
         )
 
 
 CustomPotential.__doc__ += general_doc
+
+
+class Table(LammpsPotentials):
+    def __init__(
+        self, *chemical_elements, x, force=None, energy=None, cutoff=None, keyword="", pair_style="table spline"
+    ):
+        if force is None and energy is None:
+            raise ValueError("Either force or energy must be defined")
+        force, energy = self._harmonize_input(x, force, energy)
+        self._check_consistency(x, force, energy)
+        if cutoff is None:
+            cutoff = x.max()
+        self.table_df = pd.DataFrame({"x": x, "E": energy, "f": force})
+        if keyword == "":
+            keyword = "-".join(self._harmonize_species(chemical_elements))
+        self.keyword = keyword
+        file_name = keyword + ".table"
+        self._file_name = str(Path().cwd() / file_name)
+        self._initialize_df(
+            pair_style=[f"{pair_style} {len(energy)}"],
+            interacting_species=[self._harmonize_species(chemical_elements)],
+            pair_coeff=[" ".join([str(cc) for cc in [file_name, keyword, cutoff]])],
+            filename=[self.file_name],
+            cutoff=None,
+        )
+
+    @property
+    def file_name(self):
+        return self._file_name
+
+    @property
+    def file_content(self):
+        txt = "# potential generated by pyiron\n\n"
+        txt += f"{self.keyword}\n"
+        txt += f"N {len(self.table_df)}\n\n"
+        df = table.table_df.copy()
+        df.index += 1
+        txt += df.to_string(header=False)
+        return txt
+
+    @staticmethod
+    def _check_consistency(x, force, energy):
+        if not (len(x) == len(force) == len(energy)):
+            raise ValueError("Lengths of x, force and energy are inconsistent")
+        if np.abs(np.sum(np.gradient(energy, x) * force) / np.sum(force**2) + 1) > 0.01:
+            warnings.warn("force and energy do not seem to be consistent")
+
+    @staticmethod
+    def _harmonize_input(x, force, energy):
+        if force is None:
+            force = -np.gradient(energy, x)
+        elif energy is None:
+            energy = np.append(0, -np.cumsum(force[1:] + force[:-1]) * np.diff(x) / 2)
+        return force, energy
+
+
+    def get_df(self):
+        with open(self.file_name, "w") as f:
+            f.write(self.file_content)
+        return super().get_df()