diff --git a/.ci_support/environment-old.yml b/.ci_support/environment-old.yml index 3657537..08a44b8 100644 --- a/.ci_support/environment-old.yml +++ b/.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 diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index 6433b7c..e4c4b99 100644 --- a/.ci_support/environment.yml +++ b/.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 diff --git a/binder/environment.yml b/binder/environment.yml index 6f6903d..14dc185 100644 --- a/binder/environment.yml +++ b/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 diff --git a/pyiron_lammps/compatibility/file.py b/pyiron_lammps/compatibility/file.py index 709b873..0d21256 100644 --- a/pyiron_lammps/compatibility/file.py +++ b/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( diff --git a/pyiron_lammps/compatibility/structure.py b/pyiron_lammps/compatibility/structure.py index 75302df..bc92dfe 100644 --- a/pyiron_lammps/compatibility/structure.py +++ b/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 diff --git a/tests/test_water_compatibility_file.py b/tests/test_water_compatibility_file.py new file mode 100644 index 0000000..643e447 --- /dev/null +++ b/tests/test_water_compatibility_file.py @@ -0,0 +1,272 @@ +import unittest +import os +import shutil +import numpy as np +import pandas +from ase.atoms import Atoms +import ase.units as units + +from pyiron_lammps.compatibility.file import lammps_file_interface_function +from pyiron_lammps.compatibility.structure import ( + LammpsStructureCompatibility, + write_lammps_datafile, + get_charge, + _find_line_by_prefix, +) + + +class TestCompatibilityFile(unittest.TestCase): + def setUp(self): + self.working_dir = os.path.abspath(os.path.join(__file__, "..", "lmp")) + self.static_path = os.path.abspath( + os.path.join("..", os.path.dirname(__file__), "static") + ) + self.keys = [ + "steps", + "natoms", + "cells", + "indices", + "forces", + "velocities", + "unwrapped_positions", + "positions", + "temperature", + "energy_pot", + "energy_tot", + "volume", + "pressures", + ] + + def tearDown(self): + if os.path.exists(self.working_dir): + shutil.rmtree(self.working_dir) + + def test_water_calculation(self): + density = 1.0e-24 # g/A^3 + n_mols = 27 + mol_mass_water = 18.015 # g/mol + + # Determining the supercell size size + mass = mol_mass_water * n_mols / units.mol # g + vol_h2o = mass / density # in A^3 + a = vol_h2o ** (1.0 / 3.0) # A + + # Constructing the unitcell + n = int(round(n_mols ** (1.0 / 3.0))) + + dx = 0.7 + r_O = [0, 0, 0] + r_H1 = [dx, dx, 0] + r_H2 = [-dx, dx, 0] + unit_cell = (a / n) * np.eye(3) + water_potential = pandas.DataFrame( + { + "Name": ["H2O_tip3p"], + "Filename": [[]], + "Model": ["TIP3P"], + "Species": [["H", "O"]], + "Config": [ + [ + "# @potential_species H_O ### species in potential\n", + "# W.L. Jorgensen et.al., The Journal of Chemical Physics 79, 926 (1983); https://doi.org/10.1063/1.445869\n", + "#\n", + "\n", + "units real\n", + "dimension 3\n", + "atom_style full\n", + "\n", + "# create groups ###\n", + "group O type 2\n", + "group H type 1\n", + "\n", + "## set charges - beside manually ###\n", + "set group O charge -0.830\n", + "set group H charge 0.415\n", + "\n", + "### TIP3P Potential Parameters ###\n", + "pair_style lj/cut/coul/long 10.0\n", + "pair_coeff * * 0.0 0.0 \n", + "pair_coeff 2 2 0.102 3.188 \n", + "bond_style harmonic\n", + "bond_coeff 1 450 0.9572\n", + "angle_style harmonic\n", + "angle_coeff 1 55 104.52\n", + "kspace_style pppm 1.0e-5\n", + "\n", + ] + ], + } + ) + water_base = Atoms( + symbols=["H", "H", "O"], + positions=[r_H1, r_H2, r_O], + cell=unit_cell, + pbc=True, + ) + water = water_base.repeat([n, n, n]) + _, parsed_output, job_crashed = lammps_file_interface_function( + working_directory="lmp", + structure=water, + potential=water_potential, + calc_mode="md", + calc_kwargs={"temperature": 300, "n_ionic_steps": 1e4, "time_step": 0.01}, + units="real", + lmp_command="cp " + + str(os.path.join(self.static_path, "compatibility_output")) + + "/* .", + resource_path=None, + ) + self.assertFalse(job_crashed) + for key in self.keys: + self.assertIn(key, parsed_output["generic"]) + + +class TestLammpsStructureCompatibility(unittest.TestCase): + def setUp(self): + self.output_folder = os.path.abspath( + os.path.join(__file__, "..", "structure_comp") + ) + os.makedirs(self.output_folder, exist_ok=True) + + def tearDown(self): + if os.path.exists(self.output_folder): + shutil.rmtree(self.output_folder) + + def test_structure_full(self): + bond_dict = { + "O": { + "element_list": ["H"], + "cutoff_list": [1.2], + "max_bond_list": [2], + "bond_type_list": [1], + "angle_type_list": [1], + } + } + q_dict = {"H": 0.41, "O": -0.82} + ls = LammpsStructureCompatibility( + atom_type="full", bond_dict=bond_dict, q_dict=q_dict + ) + ls.el_eam_lst = ["H", "O"] + + atoms = Atoms( + "H2O", + positions=[(0, 0.757, 0.586), (0, -0.757, 0.586), (0, 0, 0)], + cell=np.eye(3) * 15, + ) + ls.structure = atoms + output_lines = ls._string_input.split("\n") + + atoms_part_found = False + bonds_part_found = False + angles_part_found = False + + for i, line in enumerate(output_lines): + if line.strip() == "Atoms": + atoms_part_found = True + self.assertIn("1 1 1 0.410000", output_lines[i + 2]) + self.assertIn("2 1 1 0.410000", output_lines[i + 3]) + self.assertIn("3 1 2 -0.820000", output_lines[i + 4]) + if line.strip() == "Bonds": + bonds_part_found = True + line1 = output_lines[i + 2].strip() + line2 = output_lines[i + 3].strip() + self.assertTrue( + ("1 1 3 1" == line1 and "2 1 3 2" == line2) + or ("1 1 3 2" == line1 and "2 1 3 1" == line2) + ) + + if line.strip() == "Angles": + angles_part_found = True + self.assertTrue( + "1 1 1 3 2" in output_lines[i + 2] + or "1 1 2 3 1" in output_lines[i + 2] + ) + + self.assertTrue(atoms_part_found) + self.assertTrue(bonds_part_found) + self.assertTrue(angles_part_found) + + def test_structure_full_no_angles(self): + bond_dict = { + "O": { + "element_list": ["H"], + "cutoff_list": [1.2], + "max_bond_list": [2], + "bond_type_list": [1], + "angle_type_list": [None], + } + } + q_dict = {"H": 0.41, "O": -0.82} + ls = LammpsStructureCompatibility( + atom_type="full", bond_dict=bond_dict, q_dict=q_dict + ) + ls.el_eam_lst = ["H", "O"] + + atoms = Atoms( + "H2O", + positions=[(0, 0.757, 0.586), (0, -0.757, 0.586), (0, 0, 0)], + cell=np.eye(3) * 15, + ) + ls.structure = atoms + + self.assertNotIn("Angles", ls._string_input) + + def test_write_lammps_datafile_full(self): + atoms = Atoms( + "H2O", + positions=[(0, 0.757, 0.586), (0, -0.757, 0.586), (0, 0, 0)], + cell=np.eye(3) * 15, + ) + potential_lst = ["set group H charge 0.41", "set group O charge -0.82"] + + write_lammps_datafile( + structure=atoms, + potential_elements=["H", "O"], + file_name="lammps.data.full", + working_directory=self.output_folder, + atom_type="full", + potential_lst=potential_lst, + ) + + with open(os.path.join(self.output_folder, "lammps.data.full"), "r") as f: + content = f.read() + + self.assertIn("Atoms", content) + self.assertIn("Bonds", content) + self.assertIn("Angles", content) + self.assertIn("1 1 1 0.410000", content) + self.assertIn("2 1 1 0.410000", content) + self.assertIn("3 1 2 -0.820000", content) + self.assertTrue( + ("1 1 3 1" in content and "2 1 3 2" in content) + or ("1 1 3 2" in content and "2 1 3 1" in content) + ) + self.assertTrue("1 1 1 3 2" in content or "1 1 2 3 1" in content) + + def test_get_charge(self): + potential_lines = [ + "set group Fe charge 2.0", + "set group O charge -2.0", + " set group Si charge 4.0 # comment", + ] + self.assertEqual(get_charge(potential_lines, "Fe"), 2.0) + self.assertEqual(get_charge(potential_lines, "O"), -2.0) + self.assertEqual(get_charge(potential_lines, "Si"), 4.0) + with self.assertRaises(NameError): + get_charge(potential_lines, "C") + + def test_find_line_by_prefix(self): + lines = ["some line", "prefix then something", "another line"] + self.assertEqual( + _find_line_by_prefix(lines, "prefix"), ["prefix", "then", "something"] + ) + with self.assertRaises(ValueError): + _find_line_by_prefix(lines, "nonexistent") + + def test_structure_setter_charge(self): + ls = LammpsStructureCompatibility(atom_type="charge") + ls.el_eam_lst = ["Fe"] + atoms = Atoms("Fe", positions=[[0, 0, 0]], cell=np.eye(3)) + atoms.set_initial_charges([1.5]) + ls.structure = atoms + self.assertIn("1 1 1.500000 0.000000", ls._string_input)