diff --git a/pyiron_lammps/structure.py b/pyiron_lammps/structure.py index f87ae15c..d6fd43e4 100644 --- a/pyiron_lammps/structure.py +++ b/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) diff --git a/tests/test_structure.py b/tests/test_structure.py index a78dd113..7b7377a5 100644 --- a/tests/test_structure.py +++ b/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", + ], + )