Test write structure

pyiron_lammps/structure.py

@@ -10,6 +10,7 @@
 from collections import OrderedDict
 
 import numpy as np
+from ase.data import atomic_masses, atomic_numbers
 
 from pyiron_lammps.units import UnitConverter
 
@@ -187,7 +188,7 @@ class LammpsStructure(object):
         input_file_name:
     """
 
-    def __init__(self, bond_dict=None, job=None):
+    def __init__(self, bond_dict=None, units="metal"):
         self._string_input = ""
         self._structure = None
         self._potential = None
@@ -197,7 +198,7 @@ def __init__(self, bond_dict=None, job=None):
         self.digits = 10
         self._bond_dict = bond_dict
         self._force_skewed = False
-        self._job = job
+        self._units = units
         self._molecule_ids = []
 
     @property
@@ -237,16 +238,18 @@ def structure(self, structure):
         else:  # self.atom_type == 'atomic'
             input_str = self.structure_atomic()
 
-        if self._structure.velocities is not None:
-            uc = UnitConverter(self._job.units)
-            self._structure.velocities *= uc.pyiron_to_lammps("velocity")
+        if self._structure.get_velocities() is not None:
+            uc = UnitConverter(self._units)
+            self._structure.set_velocities(
+                self._structure.get_velocities() * uc.pyiron_to_lammps("velocity")
+            )
             vels = self.rotate_velocities(self._structure)
             input_str += "Velocities\n\n"
-            if self._structure.dimension == 3:
+            if len(self._structure.positions[0]) == 3:
                 format_str = "{0:d} {1:f} {2:f} {3:f}\n"
                 for id_atom, (x, y, z) in enumerate(vels, start=1):
                     input_str += format_str.format(id_atom, x, y, z)
-            if self._structure.dimension == 2:
+            elif len(self._structure.positions[0]) == 2:
                 format_str = "{0:d} {1:f} {2:f}\n"
                 for id_atom, (x, y) in enumerate(vels, start=1):
                     input_str += format_str.format(id_atom, x, y)
@@ -321,7 +324,7 @@ def lammps_header(
 
         masses = "Masses\n\n"
         for el, idx in species_lammps_id_dict.items():
-            mass = structure._pse[el].AtomicMass
+            mass = atomic_masses[atomic_numbers[el]]
             masses += "{0:3d} {1:f}  # ({2}) \n".format(idx, mass, el)
 
         return atomtypes + "\n" + cell_dimensions + "\n" + masses + "\n"
@@ -343,7 +346,7 @@ def simulation_cell(self):
             + "0. {} zlo zhi\n".format(zhi)
         )
 
-        if self.structure.is_skewed() or self._force_skewed:
+        if is_skewed(self.structure) or self._force_skewed:
             simulation_cell += "{0} {1} {2} xy xz yz\n".format(xy, xz, yz)
 
         return simulation_cell
@@ -652,7 +655,7 @@ def structure_atomic(self):
 
         el_lst = self._structure.get_chemical_symbols()
         for id_atom, (el, coord) in enumerate(zip(el_lst, coords)):
-            dim = self._structure.dimension
+            dim = len(self._structure.positions[0])
             c = np.zeros(3)
             c[:dim] = coord
             atoms += (
@@ -696,7 +699,7 @@ def rotate_velocities(self, structure):
             (list): List of rotated velocities
         """
         prism = UnfoldingPrism(self._structure.cell)
-        vels = [prism.pos_to_lammps(vel) for vel in structure.velocities]
+        vels = [prism.pos_to_lammps(vel) for vel in structure.get_velocities()]
         return vels
 
     def write_file(self, file_name, cwd=None):
@@ -715,9 +718,36 @@ def write_file(self, file_name, cwd=None):
                 f.write(line)
 
 
-def write_lammps_datafile(structure, file_name="lammps.data", cwd=None):
-    lammps_str = LammpsStructure()
-    lammps_str.el_eam_lst = structure.get_species_symbols()
+def is_skewed(structure, tolerance=1.0e-8):
+    """
+    Check whether the simulation box is skewed/sheared. The algorithm compares the box volume
+    and the product of the box length in each direction. If these numbers do not match, the box
+    is considered to be skewed and the function returns True
+
+    Args:
+        tolerance (float): Relative tolerance above which the structure is considered as skewed
+
+    Returns:
+        (bool): Whether the box is skewed or not.
+    """
+    volume = structure.get_volume()
+    prod = np.linalg.norm(structure.cell, axis=-1).prod()
+    if volume > 0:
+        if abs(volume - prod) / volume < tolerance:
+            return False
+    return True
+
+
+def write_lammps_datafile(
+    structure,
+    el_eam_lst,
+    bond_dict=None,
+    units="metal",
+    file_name="lammps.data",
+    cwd=None,
+):
+    lammps_str = LammpsStructure(bond_dict=bond_dict, units=units)
+    lammps_str.el_eam_lst = el_eam_lst
     lammps_str.structure = structure
     lammps_str.write_file(file_name=file_name, cwd=cwd)
 

tests/test_structure.py

@@ -1,10 +1,23 @@
 import unittest
 import numpy as np
+import os
+from shutil import rmtree
 from ase.build import bulk
-from pyiron_lammps.structure import structure_to_lammps
+from pyiron_lammps.structure import (
+    structure_to_lammps,
+    write_lammps_datafile,
+)
 
 
 class TestLammpsStructure(unittest.TestCase):
+    def setUp(self):
+        self.output_folder = os.path.abspath(os.path.join(__file__, "..", "structure"))
+        os.makedirs(self.output_folder, exist_ok=True)
+
+    @classmethod
+    def tearDownClass(cls):
+        rmtree(os.path.abspath(os.path.join(__file__, "..", "structure")))
+
     def test_structure_to_lammps_with_velocity(self):
         structure = bulk("Al", a=4.05)
         structure.set_velocities([[1.0, 1.0, 1.0]])
@@ -51,3 +64,87 @@ def test_structure_to_lammps_without_velocity(self):
                 )
             )
         )
+
+    def test_structure_atomic_non_cubic(self):
+        structure = bulk("Al", a=4.05)
+        write_lammps_datafile(
+            structure=structure,
+            el_eam_lst=["Ni", "Al", "H"],
+            file_name="lammps.data",
+            units="metal",
+            cwd=self.output_folder,
+        )
+        with open(os.path.join(self.output_folder, "lammps.data"), "r") as f:
+            self.assertEqual(
+                [
+                    l
+                    for l in f.readlines()
+                    if "xlo xhi" not in l
+                    and "ylo yhi" not in l
+                    and "zlo zhi" not in l
+                    and "xy xz yz" not in l
+                ],
+                [
+                    "Start File for LAMMPS \n",
+                    "1 atoms \n",
+                    "3 atom types \n",
+                    "\n",
+                    "\n",
+                    "Masses\n",
+                    "\n",
+                    "  1 58.693400  # (Ni) \n",
+                    "  2 26.981538  # (Al) \n",
+                    "  3 1.008000  # (H) \n",
+                    "\n",
+                    "Atoms\n",
+                    "\n",
+                    "1 2 0.000000000000000 0.000000000000000 0.000000000000000\n",
+                    "\n",
+                    "Velocities\n",
+                    "\n",
+                    "1 0.000000 0.000000 0.000000\n",
+                ],
+            )
+
+    def test_structure_atomic_cubic(self):
+        structure = bulk("Al", a=4.0, cubic=True)
+        write_lammps_datafile(
+            structure=structure,
+            el_eam_lst=["Ni", "Al", "H"],
+            file_name="lammps_cubic.data",
+            units="metal",
+            cwd=self.output_folder,
+        )
+        with open(os.path.join(self.output_folder, "lammps_cubic.data"), "r") as f:
+            self.assertEqual(
+                f.readlines(),
+                [
+                    "Start File for LAMMPS \n",
+                    "4 atoms \n",
+                    "3 atom types \n",
+                    "\n",
+                    "0. 4.000000000000000 xlo xhi\n",
+                    "0. 4.000000000000000 ylo yhi\n",
+                    "0. 4.000000000000000 zlo zhi\n",
+                    "\n",
+                    "Masses\n",
+                    "\n",
+                    "  1 58.693400  # (Ni) \n",
+                    "  2 26.981538  # (Al) \n",
+                    "  3 1.008000  # (H) \n",
+                    "\n",
+                    "Atoms\n",
+                    "\n",
+                    "1 2 0.000000000000000 0.000000000000000 0.000000000000000\n",
+                    "2 2 0.000000000000000 2.000000000000000 2.000000000000000\n",
+                    "3 2 2.000000000000000 0.000000000000000 2.000000000000000\n",
+                    "4 2 2.000000000000000 2.000000000000000 0.000000000000000\n",
+                    "\n",
+                    "Velocities\n",
+                    "\n",
+                    "1 0.000000 0.000000 0.000000\n",
+                    "2 0.000000 0.000000 0.000000\n",
+                    "3 0.000000 0.000000 0.000000\n",
+                    "4 0.000000 0.000000 0.000000\n",
+                ],
+            )