Fixes for manually specified potentials

pyiron_lammps/compatibility/file.py

@@ -88,18 +88,24 @@ def lammps_file_interface_function(
         calc_kwargs = {}
 
     os.makedirs(working_directory, exist_ok=True)
-    if isinstance(potential, str):
-        potential_dataframe = get_potential_by_name(
-            potential_name=potential, resource_path=resource_path
-        )
-    elif isinstance(potential, pandas.DataFrame):
-        potential_dataframe = potential.iloc[0]
-    elif isinstance(potential, pandas.Series):
-        potential_dataframe = potential
-    else:
-        raise TypeError()
-    lmp_str_lst = lammps_file_initialization(structure=structure)
-    lmp_str_lst += potential_dataframe["Config"]
+    potential_lst, potential_replace, species = _get_potential(
+        potential=potential, resource_path=resource_path
+    )
+
+    lmp_str_lst = []
+    atom_type = "atomic"
+    for l in lammps_file_initialization(structure=structure):
+        if l.startswith("units") and "units" in potential_replace:
+            lmp_str_lst.append(potential_replace["units"])
+        elif l.startswith("atom_style") and "atom_style" in potential_replace:
+            lmp_str_lst.append(potential_replace["atom_style"])
+            atom_type = potential_replace["atom_style"].split()[-1]
+        elif l.startswith("dimension") and "dimension" in potential_replace:
+            lmp_str_lst.append(potential_replace["dimension"])
+        else:
+            lmp_str_lst.append(l)
+
+    lmp_str_lst += potential_lst
     lmp_str_lst += ["variable dumptime equal {} ".format(calc_kwargs.get("n_print", 1))]
     lmp_str_lst += [
         "dump 1 all custom ${dumptime} dump.out id type xsu ysu zsu fx fy fz vx vy vz",
@@ -161,11 +167,12 @@ def lammps_file_interface_function(
 
     write_lammps_datafile(
         structure=structure,
-        potential_elements=potential_dataframe["Species"],
+        potential_elements=species,
         bond_dict=None,
         units=units,
         file_name="lammps.data",
         working_directory=working_directory,
+        atom_type=atom_type,
     )
 
     shell = subprocess.check_output(
@@ -178,7 +185,7 @@ def lammps_file_interface_function(
     output = parse_lammps_output(
         working_directory=working_directory,
         structure=structure,
-        potential_elements=potential_dataframe["Species"],
+        potential_elements=species,
         units=units,
         prism=None,
         dump_h5_file_name="dump.h5",
@@ -222,3 +229,30 @@ def _modify_input_dict(
         return lmp_tmp_lst
     else:
         return lmp_str_lst
+
+
+def _get_potential(potential, resource_path: Optional[str] = None):
+    if isinstance(potential, str):
+        potential_dataframe = get_potential_by_name(
+            potential_name=potential, resource_path=resource_path
+        )
+    elif isinstance(potential, pandas.DataFrame):
+        potential_dataframe = potential.iloc[0]
+    elif isinstance(potential, pandas.Series):
+        potential_dataframe = potential
+    else:
+        raise TypeError()
+
+    potential_replace = {}
+    potential_lst = []
+    for l in potential_dataframe["Config"]:
+        if l.startswith("units"):
+            potential_replace["units"] = l
+        elif l.startswith("atom_style"):
+            potential_replace["atom_style"] = l
+        elif l.startswith("dimension"):
+            potential_replace["dimension"] = l
+        else:
+            potential_lst.append(l)
+
+    return potential_lst, potential_replace, potential_dataframe["Species"]

pyiron_lammps/potential.py

@@ -1,6 +1,6 @@
 import os
 from pathlib import Path
-from typing import Union
+from typing import Optional, Union
 
 import pandas
 from ase.atoms import Atoms
@@ -346,7 +346,7 @@ def get_potential_dataframe(structure: Atoms, resource_path=None):
     )
 
 
-def get_potential_by_name(potential_name: str, resource_path=None):
+def get_potential_by_name(potential_name: str, resource_path: Optional[str] = None):
     if resource_path is None:
         resource_path = get_resource_path_from_conda()
     df = LammpsPotentialFile(resource_path=resource_path).list()

tests/test_compatibility_file.py

@@ -3,7 +3,10 @@
 import shutil
 from ase.build import bulk
 import pandas
-from pyiron_lammps.compatibility.file import lammps_file_interface_function
+from pyiron_lammps.compatibility.file import (
+    lammps_file_interface_function,
+    _get_potential,
+)
 from pyiron_lammps.potential import get_potential_by_name
 
 
@@ -271,10 +274,43 @@ def test_calc_md_nvt_langevin(self):
 
     def test_calc_md_nve(self):
         calc_kwargs = {"n_print": 100, "langevin": True}
+        potential = [
+            "# Bouhadja et al., J. Chem. Phys. 138, 224510 (2013) \n",
+            "units metal\n",
+            "dimension 3\n",
+            "atom_style charge\n",
+            "\n",
+            "# create groups ###\n",
+            "group Al type 1\n",
+            "group Ca type 2\n",
+            "group O type 3\n",
+            "group Si type 4\n",
+            "\n### set charges ###\n",
+            "set type 1 charge 1.8\n",
+            "set type 2 charge 1.2\n",
+            "set type 3 charge -1.2\n",
+            "set type 4 charge 2.4\n",
+            "\n### Bouhadja Born-Mayer-Huggins + Coulomb Potential Parameters ###\n",
+            "pair_style born/coul/dsf 0.25 8.0\n",
+            "pair_coeff 1 1 0.002900 0.068000 1.570400 14.049800 0.000000\n",
+            "pair_coeff 1 2 0.003200 0.074000 1.957200 17.171000 0.000000\n",
+            "pair_coeff 1 3 0.007500 0.164000 2.606700 34.574700 0.000000\n",
+            "pair_coeff 1 4 0.002500 0.057000 1.505600 18.811600 0.000000\n",
+            "pair_coeff 2 2 0.003500 0.080000 2.344000 20.985600 0.000000\n",
+            "pair_coeff 2 3 0.007700 0.178000 2.993500 42.255600 0.000000\n",
+            "pair_coeff 2 4 0.002700 0.063000 1.892400 22.990700 0.000000\n",
+            "pair_coeff 3 3 0.012000 0.263000 3.643000 85.084000 0.000000\n",
+            "pair_coeff 3 4 0.007000 0.156000 2.541900 46.293000 0.000000\n",
+            "pair_coeff 4 4 0.001200 0.046000 1.440800 25.187300 0.000000\n",
+            "\npair_modify shift yes\n",
+        ]
+        element_lst = ["Al", "Ca", "O", "Si"]
         shell_output, parsed_output, job_crashed = lammps_file_interface_function(
             working_directory=self.working_dir,
             structure=self.structure,
-            potential=self.potential,
+            potential=pandas.DataFrame(
+                {"Config": [potential], "Species": [element_lst]}
+            ),
             calc_mode="md",
             calc_kwargs=calc_kwargs,
             units=self.units,
@@ -292,15 +328,40 @@ def test_calc_md_nve(self):
             "units metal\n",
             "dimension 3\n",
             "boundary p p p\n",
-            "atom_style atomic\n",
+            "atom_style charge\n",
+            "\n",
             "read_data lammps.data\n",
-            "pair_style eam/alloy\n",
+            "# Bouhadja et al., J. Chem. Phys. 138, 224510 (2013) \n",
+            "# create groups ###\n",
+            "group Al type 1\n",
+            "group Ca type 2\n",
+            "group O type 3\n",
+            "group Si type 4\n",
+            "### set charges ###\n",
+            "set type 1 charge 1.8\n",
+            "set type 2 charge 1.2\n",
+            "set type 3 charge -1.2\n",
+            "set type 4 charge 2.4\n",
+            "### Bouhadja Born-Mayer-Huggins + Coulomb Potential Parameters ###\n",
+            "pair_style born/coul/dsf 0.25 8.0\n",
+            "pair_coeff 1 1 0.002900 0.068000 1.570400 14.049800 0.000000\n",
+            "pair_coeff 1 2 0.003200 0.074000 1.957200 17.171000 0.000000\n",
+            "pair_coeff 1 3 0.007500 0.164000 2.606700 34.574700 0.000000\n",
+            "pair_coeff 1 4 0.002500 0.057000 1.505600 18.811600 0.000000\n",
+            "pair_coeff 2 2 0.003500 0.080000 2.344000 20.985600 0.000000\n",
+            "pair_coeff 2 3 0.007700 0.178000 2.993500 42.255600 0.000000\n",
+            "pair_coeff 2 4 0.002700 0.063000 1.892400 22.990700 0.000000\n",
+            "pair_coeff 3 3 0.012000 0.263000 3.643000 85.084000 0.000000\n",
+            "pair_coeff 3 4 0.007000 0.156000 2.541900 46.293000 0.000000\n",
+            "pair_coeff 4 4 0.001200 0.046000 1.440800 25.187300 0.000000\n",
+            "pair_modify shift yes\n",
             "variable dumptime equal 100 \n",
             "dump 1 all custom ${dumptime} dump.out id type xsu ysu zsu fx fy fz vx vy vz\n",
             'dump_modify 1 sort id format line "%d %d %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g"\n',
             "fix ensemble all nve\n",
             "variable thermotime equal 100 \n",
             "timestep 0.001\n",
+            "velocity all create None 80996 dist gaussian\n",
             "thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol\n",
             "thermo_modify format float %20.15g\n",
             "thermo ${thermotime}\n",
@@ -460,3 +521,63 @@ def test_calc_minimize_pressure_3d(self):
         ]
         for line in content_expected:
             self.assertIn(line, content)
+
+
+class TestGlassPotential(unittest.TestCase):
+    def setUp(self):
+        self.static_path = os.path.abspath(
+            os.path.join("..", os.path.dirname(__file__), "static")
+        )
+
+    def test_bouhadja(self):
+        potential = [
+            "# Bouhadja et al., J. Chem. Phys. 138, 224510 (2013) \n",
+            "units metal\n",
+            "dimension 3\n",
+            "atom_style charge\n",
+            "\n",
+            "# create groups ###\n",
+            "group Al type 1\n",
+            "group Ca type 2\n",
+            "group O type 3\n",
+            "group Si type 4\n",
+            "\n### set charges ###\n",
+            "set type 1 charge 1.8\n",
+            "set type 2 charge 1.2\n",
+            "set type 3 charge -1.2\n",
+            "set type 4 charge 2.4\n",
+            "\n### Bouhadja Born-Mayer-Huggins + Coulomb Potential Parameters ###\n",
+            "pair_style born/coul/dsf 0.25 8.0\n",
+            "pair_coeff 1 1 0.002900 0.068000 1.570400 14.049800 0.000000\n",
+            "pair_coeff 1 2 0.003200 0.074000 1.957200 17.171000 0.000000\n",
+            "pair_coeff 1 3 0.007500 0.164000 2.606700 34.574700 0.000000\n",
+            "pair_coeff 1 4 0.002500 0.057000 1.505600 18.811600 0.000000\n",
+            "pair_coeff 2 2 0.003500 0.080000 2.344000 20.985600 0.000000\n",
+            "pair_coeff 2 3 0.007700 0.178000 2.993500 42.255600 0.000000\n",
+            "pair_coeff 2 4 0.002700 0.063000 1.892400 22.990700 0.000000\n",
+            "pair_coeff 3 3 0.012000 0.263000 3.643000 85.084000 0.000000\n",
+            "pair_coeff 3 4 0.007000 0.156000 2.541900 46.293000 0.000000\n",
+            "pair_coeff 4 4 0.001200 0.046000 1.440800 25.187300 0.000000\n",
+            "\npair_modify shift yes\n",
+        ]
+        element_lst = ["Al", "Ca", "O", "Si"]
+        potential_lst, potential_replace, species = _get_potential(
+            potential=pandas.DataFrame(
+                {"Config": [potential], "Species": [element_lst]}
+            ),
+            resource_path=os.path.join(self.static_path, "potential"),
+        )
+        for i, l in enumerate(potential):
+            if i in [1, 2, 3]:
+                self.assertFalse(l in potential_lst)
+            else:
+                self.assertTrue(l in potential_lst)
+
+        for k, v in {
+            "units": "units metal\n",
+            "dimension": "dimension 3\n",
+            "atom_style": "atom_style charge\n",
+        }.items():
+            self.assertEqual(potential_replace[k], v)
+
+        self.assertEqual(species, element_lst)