Lammps function water

.ci_support/environment-old.yml

@@ -8,3 +8,4 @@ dependencies:
 - numpy =1.26.0
 - pandas =2.0.3
 - scipy =1.11.1
+- structuretoolkit =0.0.14

.ci_support/environment.yml

@@ -8,3 +8,4 @@ dependencies:
 - numpy =2.3.5
 - pandas =2.3.3
 - scipy =1.16.3
+- structuretoolkit =0.0.38

binder/environment.yml

@@ -9,3 +9,4 @@ dependencies:
 - lammps =2024.08.29=*_mpi_openmpi_*
 - hatchling =1.27.0
 - hatch-vcs =0.4.0
+- structuretoolkit =0.0.38

pyiron_lammps/compatibility/file.py

@@ -11,9 +11,9 @@
     calc_static,
 )
 from pyiron_lammps.compatibility.constraints import set_selective_dynamics
+from pyiron_lammps.compatibility.structure import write_lammps_datafile
 from pyiron_lammps.output import parse_lammps_output
 from pyiron_lammps.potential import get_potential_by_name
-from pyiron_lammps.structure import write_lammps_datafile
 
 
 def lammps_file_interface_function(
@@ -173,6 +173,7 @@ def lammps_file_interface_function(
         file_name="lammps.data",
         working_directory=working_directory,
         atom_type=atom_type,
+        potential_lst=potential_lst,
     )
 
     shell = subprocess.check_output(

pyiron_lammps/compatibility/structure.py

@@ -1,5 +1,5 @@
 from collections import OrderedDict
-from typing import Dict, Optional
+from typing import Dict, Optional, Union
 
 import numpy as np
 from ase.atoms import Atoms
@@ -15,9 +15,11 @@ def __init__(
         bond_dict: Optional[Dict] = None,
         units: str = "metal",
         atom_type: str = "atomic",
+        q_dict: Optional[Dict] = None,
     ):
         super().__init__(bond_dict=bond_dict, units=units, atom_type=atom_type)
         self._molecule_ids = []
+        self.q_dict = q_dict
 
     @property
     def structure(self) -> Optional[Atoms]:
@@ -39,11 +41,12 @@ def structure(self, structure):
 
         """
         self._structure = structure
-        if self.atom_type == "full":
+        self._molecule_ids = np.ones(len(self.structure), dtype=int)
+        if self._atom_type == "full":
             input_str = self.structure_full()
-        elif self.atom_type == "bond":
+        elif self._atom_type == "bond":
             input_str = self.structure_bond()
-        elif self.atom_type == "charge":
+        elif self._atom_type == "charge":
             input_str = self.structure_charge()
         else:  # self.atom_type == 'atomic'
             input_str = self.structure_atomic()
@@ -57,7 +60,6 @@ def structure_bond(self):
 
         """
         species_lammps_id_dict = self.get_lammps_id_dict(self.el_eam_lst)
-        self.molecule_ids = None
         # analyze structure to get molecule_ids, bonds, angles etc
         coords = self.rotate_positions(self._structure)
 
@@ -70,7 +72,7 @@ def structure_bond(self):
         format_str = "{0:d} {1:d} {2:d} {3:f} {4:f} {5:f} "
         if len(self._structure.positions[0]) == 3:
             for id_atom, (x, y, z) in enumerate(coords):
-                id_mol = self.molecule_ids[id_atom]
+                id_mol = self._molecule_ids[id_atom]
                 atoms += (
                     format_str.format(
                         id_atom + 1,
@@ -84,7 +86,7 @@ def structure_bond(self):
                 )
         elif len(self._structure.positions[0]) == 2:
             for id_atom, (x, y) in enumerate(coords):
-                id_mol = self.molecule_ids[id_atom]
+                id_mol = self._molecule_ids[id_atom]
                 atoms += (
                     format_str.format(
                         id_atom + 1,
@@ -115,26 +117,22 @@ def structure_bond(self):
                 bond_type[i, j] = count
                 bond_type[j, i] = count
 
-        if self.structure.bonds is None:
-            if self.cutoff_radius is None:
-                bonds_lst = get_bonds(structure=self.structure, max_shells=1)
-            else:
-                bonds_lst = get_bonds(
-                    structure=self.structure, radius=self.cutoff_radius
-                )
-            bonds = []
+        if self.cutoff_radius is None:
+            bonds_lst = get_bonds(structure=self.structure, max_shells=1)
+        else:
+            bonds_lst = get_bonds(structure=self.structure, radius=self.cutoff_radius)
+        bonds = []
 
-            for ia, i_bonds in enumerate(bonds_lst):
-                el_i = el_dict[elements[ia]]
-                for el_j, b_lst in i_bonds.items():
-                    b_type = bond_type[el_i][el_dict[el_j]]
-                    for i_shell, ib_shell_lst in enumerate(b_lst):
-                        for ib in np.unique(ib_shell_lst):
-                            if ia < ib:  # avoid double counting of bonds
-                                bonds.append([ia + 1, ib + 1, b_type])
+        for ia, i_bonds in enumerate(bonds_lst):
+            el_i = el_dict[elements[ia]]
+            for el_j, b_lst in i_bonds.items():
+                b_type = bond_type[el_i][el_dict[el_j]]
+                for i_shell, ib_shell_lst in enumerate(b_lst):
+                    for ib in np.unique(ib_shell_lst):
+                        if ia < ib:  # avoid double counting of bonds
+                            bonds.append([ia + 1, ib + 1, b_type])
 
-            self.structure.bonds = np.array(bonds)
-        bonds = self.structure.bonds
+            bonds = np.array(bonds)
 
         bonds_str = "Bonds \n\n"
         for i_bond, (i_a, i_b, b_type) in enumerate(bonds):
@@ -165,15 +163,9 @@ def structure_full(self):
 
         """
         species_lammps_id_dict = self.get_lammps_id_dict(self.el_eam_lst)
-        self.molecule_ids = None
         coords = self.rotate_positions(self._structure)
 
         # extract electric charges from potential file
-        q_dict = {
-            species_name: self.potential.get_charge(species_name)
-            for species_name in set(self.structure.get_chemical_symbols())
-        }
-
         bonds_lst, angles_lst = [], []
         bond_type_lst, angle_type_lst = [], []
         # Using a cutoff distance to draw the bonds instead of the number of neighbors
@@ -250,9 +242,9 @@ def structure_full(self):
             atoms += (
                 format_str.format(
                     id_atom + 1,
-                    self.molecule_ids[id_atom],
+                    self._molecule_ids[id_atom],
                     species_lammps_id_dict[el],
-                    q_dict[el],
+                    self.q_dict[el],
                     coord[0],
                     coord[1],
                     coord[2],
@@ -334,3 +326,98 @@ def get_bonds(
         norm_order=2,
     )
     return neighbors.get_bonds(radius=radius, max_shells=max_shells, prec=prec)
+
+
+def write_lammps_datafile(
+    structure: Atoms,
+    potential_elements: Union[np.ndarray, list],
+    bond_dict: Optional[Dict] = None,
+    units: str = "metal",
+    file_name: str = "lammps.data",
+    working_directory: Optional[str] = None,
+    atom_type: str = "atomic",
+    potential_lst: list[str] = [],
+) -> None:
+    if atom_type == "full":
+        q_dict = {
+            el: get_charge(potential_line_lst=potential_lst, element_symbol=el)
+            for el in set(structure.get_chemical_symbols())
+        }
+        bond_dict = {
+            "O": {
+                "element_list": ["H"],
+                "cutoff_list": [2.0],
+                "max_bond_list": [2],
+                "bond_type_list": [1],
+                "angle_type_list": [1],
+            }
+        }
+    else:
+        q_dict = {}
+    lammps_str = LammpsStructureCompatibility(
+        bond_dict=bond_dict,
+        units=units,
+        atom_type=atom_type,
+        q_dict=q_dict,
+    )
+    lammps_str.el_eam_lst = potential_elements
+    lammps_str.structure = structure
+    lammps_str.write_file(file_name=file_name, cwd=working_directory)
+
+
+def _find_line_by_prefix(potential_line_lst, prefix):
+    """
+    Find a line that starts with the given prefix.  Differences in white
+    space are ignored.  Raises a ValueError if not line matches the prefix.
+
+    Args:
+        prefix (str): line prefix to search for
+
+    Returns:
+        list: words of the matching line
+
+    Raises:
+        ValueError: if not matching line was found
+    """
+
+    def isprefix(prefix, lst):
+        if len(prefix) > len(lst):
+            return False
+        return all(n == l for n, l in zip(prefix, lst))
+
+    # compare the line word by word to also match lines that differ only in
+    # whitespace
+    prefix = prefix.split()
+    for l in potential_line_lst:
+        words = l.strip().split()
+        if isprefix(prefix, words):
+            return words
+
+    raise ValueError('No line with prefix "{}" found.'.format(" ".join(prefix)))
+
+
+def get_charge(potential_line_lst, element_symbol):
+    """
+    Return charge for element. If potential does not specify a charge,
+    raise a :class:NameError.  Only makes sense for potentials
+    with pair_style "full".
+
+    Args:
+        element_symbol (str): short symbol for element
+
+    Returns:
+        float: charge speicified for the given element
+
+    Raises:
+        NameError: if potential does not specify charge for this element
+    """
+
+    try:
+        line = "set group {} charge".format(element_symbol)
+        return float(
+            _find_line_by_prefix(potential_line_lst=potential_line_lst, prefix=line)[4]
+        )
+
+    except ValueError:
+        msg = "potential does not specify charge for element {}".format(element_symbol)
+        raise NameError(msg) from None