From 147335bcf8a695f4134304905d91a5c9cee1a1c6 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 12:23:24 -0700 Subject: [PATCH 01/22] lammps interactive: move functionality to separate class --- pyiron_atomistics/lammps/interactive.py | 283 +++++------------- pyiron_atomistics/lammps/wrapper.py | 375 ++++++++++++++++++++++++ 2 files changed, 448 insertions(+), 210 deletions(-) create mode 100644 pyiron_atomistics/lammps/wrapper.py diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 88b01db7d..b858c1df3 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -14,6 +14,22 @@ from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.control import LammpsControl from pyiron_atomistics.atomistics.job.interactive import GenericInteractive +from pyiron_atomistics.lammps.wrapper import ( + interactive_lib_command, + interactive_positions_setter, + interactive_positions_getter, + interactive_cells_getter, + interactive_cells_setter, + interactive_volume_getter, + interactive_forces_getter, + interactive_structure_setter, + interactive_energy_pot_getter, + interactive_energy_tot_getter, + interactive_steps_getter, + interactive_temperatures_getter, + interactive_indices_getter, + interactive_pressures_getter +) try: # mpi4py is only supported on Linux and Mac Os X @@ -75,105 +91,67 @@ def interactive_mpi_communicator(self, comm): self._interactive_mpi_communicator = comm def _interactive_lib_command(self, command): - self._logger.debug("Lammps library: " + command) - self._interactive_library.command(command) + interactive_lib_command(lmp=self._interactive_library, logger=self._logger, command=command) def interactive_positions_getter(self): uc = UnitConverter(units=self.units) - positions = np.reshape( - np.array(self._interactive_library.gather_atoms("x", 1, 3)), - (len(self.structure), 3), + positions = interactive_positions_getter( + lmp=self._interactive_library, + number_of_atoms=len(self.structure), + prism=self._prism ) - if _check_ortho_prism(prism=self._prism): - positions = np.matmul(positions, self._prism.R.T) positions = uc.convert_array_to_pyiron_units(positions, label="positions") return positions.tolist() def interactive_positions_setter(self, positions): - if _check_ortho_prism(prism=self._prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, self._prism.R) - positions = np.array(positions).flatten() - if self.server.run_mode.interactive and self.server.cores == 1: - self._interactive_library.scatter_atoms( - "x", 1, 3, (len(positions) * c_double)(*positions) - ) - else: - self._interactive_library.scatter_atoms("x", positions) - self._interactive_lib_command("change_box all remap") + interactive_positions_setter( + lmp=self._interactive_library, + logger=self._logger, + positions=positions, + prism=self._prism, + cores=self.server.cores, + interactive=self.server.run_mode.interactive + ) def interactive_cells_getter(self): uc = UnitConverter(units=self.units) - cc = np.array( - [ - [self._interactive_library.get_thermo("lx"), 0, 0], - [ - self._interactive_library.get_thermo("xy"), - self._interactive_library.get_thermo("ly"), - 0, - ], - [ - self._interactive_library.get_thermo("xz"), - self._interactive_library.get_thermo("yz"), - self._interactive_library.get_thermo("lz"), - ], - ] - ) + cc = interactive_cells_getter(lmp=self._interactive_library) return uc.convert_array_to_pyiron_units( self._prism.unfold_cell(cc), label="cells" ) def interactive_cells_setter(self, cell): - self._prism = UnfoldingPrism(cell) - lx, ly, lz, xy, xz, yz = self._prism.get_lammps_prism() - if _check_ortho_prism(prism=self._prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" - ) - - is_skewed = self._structure_current.is_skewed(tolerance=1.0e-8) - was_skewed = self._structure_previous.is_skewed(tolerance=1.0e-8) - - if is_skewed: - if not was_skewed: - self._interactive_lib_command("change_box all triclinic") - self._interactive_lib_command( - "change_box all x final 0 %f y final 0 %f z final 0 %f \ - xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, xy, xz, yz) - ) - elif was_skewed: - self._interactive_lib_command( - "change_box all x final 0 %f y final 0 %f z final 0 %f \ - xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, 0.0, 0.0, 0.0) - ) - self._interactive_lib_command("change_box all ortho") - else: - self._interactive_lib_command( - "change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" - % (lx, ly, lz) - ) + self._prism = interactive_cells_setter( + lmp=self._interactive_library, + logger=self._logger, + cell=cell, + structure_current=self._structure_current, + structure_previous=self._structure_previous + ) def interactive_volume_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.get_thermo("vol"), label="volume" + interactive_volume_getter(lmp=self._interactive_library), + label="volume" ) def interactive_forces_getter(self): uc = UnitConverter(units=self.units) - ff = np.reshape( - np.array(self._interactive_library.gather_atoms("f", 1, 3)), - (len(self.structure), 3), + ff = interactive_forces_getter( + lmp=self._interactive_library, + prism=self._prism, + number_of_atoms=len(self.structure) ) - if _check_ortho_prism(prism=self._prism): - ff = np.matmul(ff, self._prism.R.T) ff = uc.convert_array_to_pyiron_units(ff, label="forces") return ff.tolist() def interactive_execute(self): - self._interactive_lib_command(self._interactive_run_command) + interactive_lib_command( + lmp=self._interactive_library, + logger=self._logger, + command=self._interactive_run_command + ) def _interactive_lammps_input(self): del self.input.control["dump___1"] @@ -480,120 +458,22 @@ def interactive_fetch(self): self._logger.debug("interactive run - done") def interactive_structure_setter(self, structure): - old_symbols = self.structure.get_species_symbols() - new_symbols = structure.get_species_symbols() - if any(old_symbols != new_symbols): - raise ValueError( - f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" - ) - self._interactive_lib_command("clear") - self._set_selective_dynamics() - self._interactive_lib_command("units " + self.input.control["units"]) - self._interactive_lib_command( - "dimension " + str(self.input.control["dimension"]) + self._prism, control_dict = interactive_structure_setter( + lmp=self._interactive_library, + logger=self._logger, + structure_current=structure, + structure_previous=self.structure, + units=self.input.control["units"], + dimension=self.input.control["dimension"], + boundary=self.input.control["boundary"], + atom_style=self.input.control["atom_style"], + el_eam_lst=self.input.potential.get_element_lst(), + calc_md=self._generic_input["calc_mode"] == "md", + interactive=self.server.run_mode.interactive, + cores=self.server.cores ) - self._interactive_lib_command("boundary " + self.input.control["boundary"]) - self._interactive_lib_command("atom_style " + self.input.control["atom_style"]) - - self._interactive_lib_command("atom_modify map array") - self._prism = UnfoldingPrism(structure.cell) - if _check_ortho_prism(prism=self._prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" - ) - xhi, yhi, zhi, xy, xz, yz = self._prism.get_lammps_prism() - if self._prism.is_skewed(): - self._interactive_lib_command( - "region 1 prism" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " " - + str(xy) - + " " - + str(xz) - + " " - + str(yz) - + " units box" - ) - else: - self._interactive_lib_command( - "region 1 block" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " units box" - ) - el_struct_lst = self.structure.get_species_symbols() - el_obj_lst = self.structure.get_species_objects() - el_eam_lst = self.input.potential.get_element_lst() - if self.input.control["atom_style"] == "full": - self._interactive_lib_command( - "create_box " - + str(len(el_eam_lst)) - + " 1 " - + "bond/types 1 " - + "angle/types 1 " - + "extra/bond/per/atom 2 " - + "extra/angle/per/atom 2 " - ) - else: - self._interactive_lib_command("create_box " + str(len(el_eam_lst)) + " 1") - el_dict = {} - for id_eam, el_eam in enumerate(el_eam_lst): - if el_eam in el_struct_lst: - id_el = list(el_struct_lst).index(el_eam) - el = el_obj_lst[id_el] - el_dict[el] = id_eam + 1 - self._interactive_lib_command( - "mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass) - ) - else: - self._interactive_lib_command( - "mass {0:3d} {1:f}".format(id_eam + 1, 1.00) - ) - positions = structure.positions.flatten() - if _check_ortho_prism(prism=self._prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, self._prism.R) - positions = positions.flatten() - try: - elem_all = np.array( - [el_dict[el] for el in structure.get_chemical_elements()] - ) - except KeyError: - missing = set(structure.get_chemical_elements()).difference(el_dict.keys()) - missing = ", ".join([el.Abbreviation for el in missing]) - raise ValueError( - f"Structure contains elements [{missing}], that are not present in the potential!" - ) - if self.server.run_mode.interactive and self.server.cores == 1: - self._interactive_library.create_atoms( - n=len(structure), - id=None, - type=(len(elem_all) * c_int)(*elem_all), - x=(len(positions) * c_double)(*positions), - v=None, - image=None, - shrinkexceed=False, - ) - else: - self._interactive_library.create_atoms( - n=len(structure), - id=None, - type=elem_all, - x=positions, - v=None, - image=None, - shrinkexceed=False, - ) - self._interactive_lib_command("change_box all remap") + for key, value in control_dict.items(): + self.input.control[key.replace(" ", "___")] = value self._interactive_lammps_input() self._interactive_set_potential() @@ -680,7 +560,7 @@ def update_potential(self): def interactive_indices_getter(self): uc = UnitConverter(units=self.units) - lammps_indices = np.array(self._interactive_library.gather_atoms("type", 0, 1)) + lammps_indices = interactive_indices_getter(lmp=self._interactive_library) indices = uc.convert_array_to_pyiron_units( self.remap_indices(lammps_indices), label="indices" ) @@ -709,25 +589,29 @@ def interactive_indices_setter(self, indices): def interactive_energy_pot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.get_thermo("pe"), label="energy_pot" + interactive_energy_pot_getter(lmp=self._interactive_library), + label="energy_pot" ) def interactive_energy_tot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.get_thermo("etotal"), label="energy_tot" + interactive_energy_tot_getter(lmp=self._interactive_library), + label="energy_tot" ) def interactive_steps_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.get_thermo("step"), label="steps" + interactive_steps_getter(lmp=self._interactive_library), + label="steps" ) def interactive_temperatures_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.get_thermo("temp"), label="temperature" + interactive_temperatures_getter(lmp=self._interactive_library), \ + label="temperature" ) def interactive_stress_getter(self): @@ -763,28 +647,7 @@ def interactive_stress_getter(self): def interactive_pressures_getter(self): uc = UnitConverter(units=self.units) - pp = np.array( - [ - [ - self._interactive_library.get_thermo("pxx"), - self._interactive_library.get_thermo("pxy"), - self._interactive_library.get_thermo("pxz"), - ], - [ - self._interactive_library.get_thermo("pxy"), - self._interactive_library.get_thermo("pyy"), - self._interactive_library.get_thermo("pyz"), - ], - [ - self._interactive_library.get_thermo("pxz"), - self._interactive_library.get_thermo("pyz"), - self._interactive_library.get_thermo("pzz"), - ], - ] - ) - if _check_ortho_prism(prism=self._prism): - rotation_matrix = self._prism.R.T - pp = rotation_matrix.T @ pp @ rotation_matrix + pp = interactive_pressures_getter(lmp=self._interactive_library, prism=self._prism) return uc.convert_array_to_pyiron_units(pp, label="pressure") def interactive_close(self): diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py new file mode 100644 index 000000000..bffb82e87 --- /dev/null +++ b/pyiron_atomistics/lammps/wrapper.py @@ -0,0 +1,375 @@ +from ctypes import c_double, c_int +import numpy as np +import warnings +from pyiron_atomistics.lammps.structure import UnfoldingPrism +from pyiron_atomistics.lammps.base import _check_ortho_prism + + +def interactive_positions_getter(lmp, number_of_atoms, prism): + positions = np.reshape( + np.array(lmp.gather_atoms("x", 1, 3)), + (number_of_atoms, 3), + ) + if _check_ortho_prism(prism=prism): + positions = np.matmul(positions, prism.R.T) + return positions + + +def interactive_positions_setter(lmp, logger, positions, prism, cores, interactive): + if _check_ortho_prism(prism=prism): + positions = np.array(positions).reshape(-1, 3) + positions = np.matmul(positions, prism) + positions = np.array(positions).flatten() + if interactive and cores == 1: + lmp.scatter_atoms( + "x", 1, 3, (len(positions) * c_double)(*positions) + ) + else: + lmp.scatter_atoms("x", positions) + interactive_lib_command(lmp=lmp, logger=logger, command="change_box all remap") + + +def interactive_lib_command(lmp, logger, command): + logger.debug("Lammps library: " + command) + lmp.command(command) + + +def interactive_cells_getter(lmp): + return np.array( + [ + [lmp.get_thermo("lx"), 0, 0], + [ + lmp.get_thermo("xy"), + lmp.get_thermo("ly"), + 0, + ], + [ + lmp.get_thermo("xz"), + lmp.get_thermo("yz"), + lmp.get_thermo("lz"), + ], + ] + ) + + +def interactive_cells_setter(lmp, logger, cell, structure_current, structure_previous): + prism = UnfoldingPrism(cell) + lx, ly, lz, xy, xz, yz = prism.get_lammps_prism() + if _check_ortho_prism(prism=prism): + warnings.warn( + "Warning: setting upper trangular matrix might slow down the calculation" + ) + + is_skewed = structure_current.is_skewed(tolerance=1.0e-8) + was_skewed = structure_previous.is_skewed(tolerance=1.0e-8) + + if is_skewed: + if not was_skewed: + interactive_lib_command( + lmp=lmp, + logger=logger, + command="change_box all triclinic" + ) + interactive_lib_command( + lmp=lmp, + logger=logger, + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" % (lx, ly, lz, xy, xz, yz) + ) + elif was_skewed: + interactive_lib_command( + lmp=lmp, + logger=logger, + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" % (lx, ly, lz, 0.0, 0.0, 0.0) + ) + interactive_lib_command( + lmp=lmp, + logger=logger, + command="change_box all ortho" + ) + else: + interactive_lib_command( + lmp=lmp, + logger=logger, + command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" % (lx, ly, lz) + ) + return prism + + +def interactive_volume_getter(lmp): + return lmp.get_thermo("vol") + + +def interactive_forces_getter(lmp, prism, number_of_atoms): + ff = np.reshape( + np.array(lmp.gather_atoms("f", 1, 3)), + (number_of_atoms, 3), + ) + if _check_ortho_prism(prism=prism): + ff = np.matmul(ff, prism.R.T) + return ff + + +def interactive_structure_setter(lmp, logger, structure_current, structure_previous, units, dimension, boundary, atom_style, el_eam_lst, calc_md, interactive, cores): + old_symbols = structure_previous.get_species_symbols() + new_symbols = structure_current.get_species_symbols() + if any(old_symbols != new_symbols): + raise ValueError( + f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" + ) + interactive_lib_command(lmp=lmp, logger=logger, command="clear") + control_dict = set_selective_dynamics(structure=structure_current, calc_md=calc_md) + interactive_lib_command(lmp=lmp, logger=logger, command="units " + units) + interactive_lib_command(lmp=lmp, logger=logger, command="dimension " + str(dimension)) + interactive_lib_command(lmp=lmp, logger=logger, command="boundary " + boundary) + interactive_lib_command(lmp=lmp, logger=logger, command="atom_style " + atom_style) + + interactive_lib_command(lmp=lmp, logger=logger, command="atom_modify map array") + prism = UnfoldingPrism(structure_current.cell) + if _check_ortho_prism(prism=prism): + warnings.warn( + "Warning: setting upper trangular matrix might slow down the calculation" + ) + xhi, yhi, zhi, xy, xz, yz = prism.get_lammps_prism() + if prism.is_skewed(): + interactive_lib_command(lmp=lmp, logger=logger, command="region 1 prism" + + " 0.0 " + + str(xhi) + + " 0.0 " + + str(yhi) + + " 0.0 " + + str(zhi) + + " " + + str(xy) + + " " + + str(xz) + + " " + + str(yz) + + " units box" + ) + else: + interactive_lib_command(lmp=lmp, logger=logger, command="region 1 block" + + " 0.0 " + + str(xhi) + + " 0.0 " + + str(yhi) + + " 0.0 " + + str(zhi) + + " units box" + ) + el_struct_lst = structure_current.get_species_symbols() + el_obj_lst = structure_current.get_species_objects() + if atom_style == "full": + interactive_lib_command(lmp=lmp, logger=logger, command="create_box " + + str(len(el_eam_lst)) + + " 1 " + + "bond/types 1 " + + "angle/types 1 " + + "extra/bond/per/atom 2 " + + "extra/angle/per/atom 2 " + ) + else: + interactive_lib_command(lmp=lmp, logger=logger, command="create_box " + str(len(el_eam_lst)) + " 1") + el_dict = {} + for id_eam, el_eam in enumerate(el_eam_lst): + if el_eam in el_struct_lst: + id_el = list(el_struct_lst).index(el_eam) + el = el_obj_lst[id_el] + el_dict[el] = id_eam + 1 + interactive_lib_command(lmp=lmp, logger=logger, command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass)) + else: + interactive_lib_command(lmp=lmp, logger=logger, command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00)) + positions = structure_current.positions.flatten() + if _check_ortho_prism(prism=prism): + positions = np.array(positions).reshape(-1, 3) + positions = np.matmul(positions, prism.R) + positions = positions.flatten() + try: + elem_all = np.array( + [el_dict[el] for el in structure_current.get_chemical_elements()] + ) + except KeyError: + missing = set(structure_current.get_chemical_elements()).difference(el_dict.keys()) + missing = ", ".join([el.Abbreviation for el in missing]) + raise ValueError( + f"Structure contains elements [{missing}], that are not present in the potential!" + ) + if interactive and cores == 1: + lmp.create_atoms( + n=len(structure_current), + id=None, + type=(len(elem_all) * c_int)(*elem_all), + x=(len(positions) * c_double)(*positions), + v=None, + image=None, + shrinkexceed=False, + ) + else: + lmp.create_atoms( + n=len(structure_current), + id=None, + type=elem_all, + x=positions, + v=None, + image=None, + shrinkexceed=False, + ) + interactive_lib_command(lmp=lmp, logger=logger, command="change_box all remap") + return prism, control_dict + + +def set_selective_dynamics(structure, calc_md): + control_dict = {} + if "selective_dynamics" in structure._tag_list.keys(): + if structure.selective_dynamics._default is None: + structure.selective_dynamics._default = [True, True, True] + sel_dyn = np.logical_not(structure.selective_dynamics.list()) + # Enter loop only if constraints present + if len(np.argwhere(np.any(sel_dyn, axis=1)).flatten()) != 0: + all_indices = np.arange(len(structure), dtype=int) + constraint_xyz = np.argwhere(np.all(sel_dyn, axis=1)).flatten() + not_constrained_xyz = np.setdiff1d(all_indices, constraint_xyz) + # LAMMPS starts counting from 1 + constraint_xyz += 1 + ind_x = np.argwhere(sel_dyn[not_constrained_xyz, 0]).flatten() + ind_y = np.argwhere(sel_dyn[not_constrained_xyz, 1]).flatten() + ind_z = np.argwhere(sel_dyn[not_constrained_xyz, 2]).flatten() + constraint_xy = not_constrained_xyz[np.intersect1d(ind_x, ind_y)] + 1 + constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 + constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 + constraint_x = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + + 1 + ) + constraint_y = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + + 1 + ) + constraint_z = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + + 1 + ) + control_dict = {} + if len(constraint_xyz) > 0: + control_dict["group constraintxyz"] = "id " + " ".join( + [str(ind) for ind in constraint_xyz] + ) + control_dict[ + "fix constraintxyz" + ] = "constraintxyz setforce 0.0 0.0 0.0" + if calc_md: + control_dict[ + "velocity constraintxyz" + ] = "set 0.0 0.0 0.0" + if len(constraint_xy) > 0: + control_dict["group constraintxy"] = "id " + " ".join( + [str(ind) for ind in constraint_xy] + ) + control_dict[ + "fix constraintxy" + ] = "constraintxy setforce 0.0 0.0 NULL" + if calc_md: + control_dict[ + "velocity constraintxy" + ] = "set 0.0 0.0 NULL" + if len(constraint_yz) > 0: + control_dict["group constraintyz"] = "id " + " ".join( + [str(ind) for ind in constraint_yz] + ) + control_dict[ + "fix constraintyz" + ] = "constraintyz setforce NULL 0.0 0.0" + if calc_md: + control_dict[ + "velocity constraintyz" + ] = "set NULL 0.0 0.0" + if len(constraint_zx) > 0: + control_dict["group constraintxz"] = "id " + " ".join( + [str(ind) for ind in constraint_zx] + ) + control_dict[ + "fix constraintxz" + ] = "constraintxz setforce 0.0 NULL 0.0" + if calc_md: + control_dict[ + "velocity constraintxz" + ] = "set 0.0 NULL 0.0" + if len(constraint_x) > 0: + control_dict["group constraintx"] = "id " + " ".join( + [str(ind) for ind in constraint_x] + ) + control_dict[ + "fix constraintx" + ] = "constraintx setforce 0.0 NULL NULL" + if calc_md: + control_dict[ + "velocity constraintx" + ] = "set 0.0 NULL NULL" + if len(constraint_y) > 0: + control_dict["group constrainty"] = "id " + " ".join( + [str(ind) for ind in constraint_y] + ) + control_dict[ + "fix constrainty" + ] = "constrainty setforce NULL 0.0 NULL" + if calc_md: + control_dict[ + "velocity constrainty" + ] = "set NULL 0.0 NULL" + if len(constraint_z) > 0: + control_dict["group constraintz"] = "id " + " ".join( + [str(ind) for ind in constraint_z] + ) + control_dict[ + "fix constraintz" + ] = "constraintz setforce NULL NULL 0.0" + if calc_md: + control_dict[ + "velocity constraintz" + ] = "set NULL NULL 0.0" + return control_dict + + +def interactive_indices_getter(lmp): + return np.array(lmp.gather_atoms("type", 0, 1)) + + +def interactive_energy_pot_getter(lmp): + return lmp.get_thermo("pe") + + +def interactive_energy_tot_getter(lmp): + return lmp.get_thermo("etotal") + + +def interactive_steps_getter(lmp): + return lmp.get_thermo("step") + + +def interactive_temperatures_getter(lmp): + return lmp.get_thermo("temp") + + +def interactive_pressures_getter(lmp, prism): + pp = np.array( + [ + [ + lmp.get_thermo("pxx"), + lmp.get_thermo("pxy"), + lmp.get_thermo("pxz"), + ], + [ + lmp.get_thermo("pxy"), + lmp.get_thermo("pyy"), + lmp.get_thermo("pyz"), + ], + [ + lmp.get_thermo("pxz"), + lmp.get_thermo("pyz"), + lmp.get_thermo("pzz"), + ], + ] + ) + if _check_ortho_prism(prism=prism): + rotation_matrix = prism.R.T + pp = rotation_matrix.T @ pp @ rotation_matrix + return pp From c8a10884d8a70f7cf05813223abd9c181fd4ada3 Mon Sep 17 00:00:00 2001 From: pyiron-runner Date: Mon, 27 Feb 2023 19:50:53 +0000 Subject: [PATCH 02/22] Format black --- pyiron_atomistics/lammps/interactive.py | 36 +++--- pyiron_atomistics/lammps/wrapper.py | 148 ++++++++++++------------ 2 files changed, 93 insertions(+), 91 deletions(-) diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index b858c1df3..ef8ebac62 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -28,7 +28,7 @@ interactive_steps_getter, interactive_temperatures_getter, interactive_indices_getter, - interactive_pressures_getter + interactive_pressures_getter, ) @@ -91,14 +91,16 @@ def interactive_mpi_communicator(self, comm): self._interactive_mpi_communicator = comm def _interactive_lib_command(self, command): - interactive_lib_command(lmp=self._interactive_library, logger=self._logger, command=command) + interactive_lib_command( + lmp=self._interactive_library, logger=self._logger, command=command + ) def interactive_positions_getter(self): uc = UnitConverter(units=self.units) positions = interactive_positions_getter( lmp=self._interactive_library, number_of_atoms=len(self.structure), - prism=self._prism + prism=self._prism, ) positions = uc.convert_array_to_pyiron_units(positions, label="positions") return positions.tolist() @@ -110,7 +112,7 @@ def interactive_positions_setter(self, positions): positions=positions, prism=self._prism, cores=self.server.cores, - interactive=self.server.run_mode.interactive + interactive=self.server.run_mode.interactive, ) def interactive_cells_getter(self): @@ -126,14 +128,13 @@ def interactive_cells_setter(self, cell): logger=self._logger, cell=cell, structure_current=self._structure_current, - structure_previous=self._structure_previous + structure_previous=self._structure_previous, ) def interactive_volume_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_volume_getter(lmp=self._interactive_library), - label="volume" + interactive_volume_getter(lmp=self._interactive_library), label="volume" ) def interactive_forces_getter(self): @@ -141,7 +142,7 @@ def interactive_forces_getter(self): ff = interactive_forces_getter( lmp=self._interactive_library, prism=self._prism, - number_of_atoms=len(self.structure) + number_of_atoms=len(self.structure), ) ff = uc.convert_array_to_pyiron_units(ff, label="forces") return ff.tolist() @@ -150,7 +151,7 @@ def interactive_execute(self): interactive_lib_command( lmp=self._interactive_library, logger=self._logger, - command=self._interactive_run_command + command=self._interactive_run_command, ) def _interactive_lammps_input(self): @@ -470,7 +471,7 @@ def interactive_structure_setter(self, structure): el_eam_lst=self.input.potential.get_element_lst(), calc_md=self._generic_input["calc_mode"] == "md", interactive=self.server.run_mode.interactive, - cores=self.server.cores + cores=self.server.cores, ) for key, value in control_dict.items(): self.input.control[key.replace(" ", "___")] = value @@ -590,28 +591,27 @@ def interactive_energy_pot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( interactive_energy_pot_getter(lmp=self._interactive_library), - label="energy_pot" + label="energy_pot", ) def interactive_energy_tot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( interactive_energy_tot_getter(lmp=self._interactive_library), - label="energy_tot" + label="energy_tot", ) def interactive_steps_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_steps_getter(lmp=self._interactive_library), - label="steps" + interactive_steps_getter(lmp=self._interactive_library), label="steps" ) def interactive_temperatures_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_temperatures_getter(lmp=self._interactive_library), \ - label="temperature" + interactive_temperatures_getter(lmp=self._interactive_library), + label="temperature", ) def interactive_stress_getter(self): @@ -647,7 +647,9 @@ def interactive_stress_getter(self): def interactive_pressures_getter(self): uc = UnitConverter(units=self.units) - pp = interactive_pressures_getter(lmp=self._interactive_library, prism=self._prism) + pp = interactive_pressures_getter( + lmp=self._interactive_library, prism=self._prism + ) return uc.convert_array_to_pyiron_units(pp, label="pressure") def interactive_close(self): diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index bffb82e87..1e13ca5fa 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -21,9 +21,7 @@ def interactive_positions_setter(lmp, logger, positions, prism, cores, interacti positions = np.matmul(positions, prism) positions = np.array(positions).flatten() if interactive and cores == 1: - lmp.scatter_atoms( - "x", 1, 3, (len(positions) * c_double)(*positions) - ) + lmp.scatter_atoms("x", 1, 3, (len(positions) * c_double)(*positions)) else: lmp.scatter_atoms("x", positions) interactive_lib_command(lmp=lmp, logger=logger, command="change_box all remap") @@ -66,31 +64,28 @@ def interactive_cells_setter(lmp, logger, cell, structure_current, structure_pre if is_skewed: if not was_skewed: interactive_lib_command( - lmp=lmp, - logger=logger, - command="change_box all triclinic" + lmp=lmp, logger=logger, command="change_box all triclinic" ) interactive_lib_command( lmp=lmp, logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" % (lx, ly, lz, xy, xz, yz) + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" + % (lx, ly, lz, xy, xz, yz), ) elif was_skewed: interactive_lib_command( lmp=lmp, logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" % (lx, ly, lz, 0.0, 0.0, 0.0) - ) - interactive_lib_command( - lmp=lmp, - logger=logger, - command="change_box all ortho" + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" + % (lx, ly, lz, 0.0, 0.0, 0.0), ) + interactive_lib_command(lmp=lmp, logger=logger, command="change_box all ortho") else: interactive_lib_command( lmp=lmp, logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" % (lx, ly, lz) + command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" + % (lx, ly, lz), ) return prism @@ -109,7 +104,20 @@ def interactive_forces_getter(lmp, prism, number_of_atoms): return ff -def interactive_structure_setter(lmp, logger, structure_current, structure_previous, units, dimension, boundary, atom_style, el_eam_lst, calc_md, interactive, cores): +def interactive_structure_setter( + lmp, + logger, + structure_current, + structure_previous, + units, + dimension, + boundary, + atom_style, + el_eam_lst, + calc_md, + interactive, + cores, +): old_symbols = structure_previous.get_species_symbols() new_symbols = structure_current.get_species_symbols() if any(old_symbols != new_symbols): @@ -119,7 +127,9 @@ def interactive_structure_setter(lmp, logger, structure_current, structure_previ interactive_lib_command(lmp=lmp, logger=logger, command="clear") control_dict = set_selective_dynamics(structure=structure_current, calc_md=calc_md) interactive_lib_command(lmp=lmp, logger=logger, command="units " + units) - interactive_lib_command(lmp=lmp, logger=logger, command="dimension " + str(dimension)) + interactive_lib_command( + lmp=lmp, logger=logger, command="dimension " + str(dimension) + ) interactive_lib_command(lmp=lmp, logger=logger, command="boundary " + boundary) interactive_lib_command(lmp=lmp, logger=logger, command="atom_style " + atom_style) @@ -131,7 +141,10 @@ def interactive_structure_setter(lmp, logger, structure_current, structure_previ ) xhi, yhi, zhi, xy, xz, yz = prism.get_lammps_prism() if prism.is_skewed(): - interactive_lib_command(lmp=lmp, logger=logger, command="region 1 prism" + interactive_lib_command( + lmp=lmp, + logger=logger, + command="region 1 prism" + " 0.0 " + str(xhi) + " 0.0 " @@ -144,40 +157,56 @@ def interactive_structure_setter(lmp, logger, structure_current, structure_previ + str(xz) + " " + str(yz) - + " units box" + + " units box", ) else: - interactive_lib_command(lmp=lmp, logger=logger, command="region 1 block" + interactive_lib_command( + lmp=lmp, + logger=logger, + command="region 1 block" + " 0.0 " + str(xhi) + " 0.0 " + str(yhi) + " 0.0 " + str(zhi) - + " units box" + + " units box", ) el_struct_lst = structure_current.get_species_symbols() el_obj_lst = structure_current.get_species_objects() if atom_style == "full": - interactive_lib_command(lmp=lmp, logger=logger, command="create_box " + interactive_lib_command( + lmp=lmp, + logger=logger, + command="create_box " + str(len(el_eam_lst)) + " 1 " + "bond/types 1 " + "angle/types 1 " + "extra/bond/per/atom 2 " - + "extra/angle/per/atom 2 " + + "extra/angle/per/atom 2 ", ) else: - interactive_lib_command(lmp=lmp, logger=logger, command="create_box " + str(len(el_eam_lst)) + " 1") + interactive_lib_command( + lmp=lmp, logger=logger, command="create_box " + str(len(el_eam_lst)) + " 1" + ) el_dict = {} for id_eam, el_eam in enumerate(el_eam_lst): if el_eam in el_struct_lst: id_el = list(el_struct_lst).index(el_eam) el = el_obj_lst[id_el] el_dict[el] = id_eam + 1 - interactive_lib_command(lmp=lmp, logger=logger, command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass)) + interactive_lib_command( + lmp=lmp, + logger=logger, + command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass), + ) else: - interactive_lib_command(lmp=lmp, logger=logger, command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00)) + interactive_lib_command( + lmp=lmp, + logger=logger, + command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00), + ) positions = structure_current.positions.flatten() if _check_ortho_prism(prism=prism): positions = np.array(positions).reshape(-1, 3) @@ -188,7 +217,9 @@ def interactive_structure_setter(lmp, logger, structure_current, structure_previ [el_dict[el] for el in structure_current.get_chemical_elements()] ) except KeyError: - missing = set(structure_current.get_chemical_elements()).difference(el_dict.keys()) + missing = set(structure_current.get_chemical_elements()).difference( + el_dict.keys() + ) missing = ", ".join([el.Abbreviation for el in missing]) raise ValueError( f"Structure contains elements [{missing}], that are not present in the potential!" @@ -237,95 +268,64 @@ def set_selective_dynamics(structure, calc_md): constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 constraint_x = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] - + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + 1 ) constraint_y = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] - + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + 1 ) constraint_z = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] - + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + 1 ) control_dict = {} if len(constraint_xyz) > 0: control_dict["group constraintxyz"] = "id " + " ".join( [str(ind) for ind in constraint_xyz] ) - control_dict[ - "fix constraintxyz" - ] = "constraintxyz setforce 0.0 0.0 0.0" + control_dict["fix constraintxyz"] = "constraintxyz setforce 0.0 0.0 0.0" if calc_md: - control_dict[ - "velocity constraintxyz" - ] = "set 0.0 0.0 0.0" + control_dict["velocity constraintxyz"] = "set 0.0 0.0 0.0" if len(constraint_xy) > 0: control_dict["group constraintxy"] = "id " + " ".join( [str(ind) for ind in constraint_xy] ) - control_dict[ - "fix constraintxy" - ] = "constraintxy setforce 0.0 0.0 NULL" + control_dict["fix constraintxy"] = "constraintxy setforce 0.0 0.0 NULL" if calc_md: - control_dict[ - "velocity constraintxy" - ] = "set 0.0 0.0 NULL" + control_dict["velocity constraintxy"] = "set 0.0 0.0 NULL" if len(constraint_yz) > 0: control_dict["group constraintyz"] = "id " + " ".join( [str(ind) for ind in constraint_yz] ) - control_dict[ - "fix constraintyz" - ] = "constraintyz setforce NULL 0.0 0.0" + control_dict["fix constraintyz"] = "constraintyz setforce NULL 0.0 0.0" if calc_md: - control_dict[ - "velocity constraintyz" - ] = "set NULL 0.0 0.0" + control_dict["velocity constraintyz"] = "set NULL 0.0 0.0" if len(constraint_zx) > 0: control_dict["group constraintxz"] = "id " + " ".join( [str(ind) for ind in constraint_zx] ) - control_dict[ - "fix constraintxz" - ] = "constraintxz setforce 0.0 NULL 0.0" + control_dict["fix constraintxz"] = "constraintxz setforce 0.0 NULL 0.0" if calc_md: - control_dict[ - "velocity constraintxz" - ] = "set 0.0 NULL 0.0" + control_dict["velocity constraintxz"] = "set 0.0 NULL 0.0" if len(constraint_x) > 0: control_dict["group constraintx"] = "id " + " ".join( [str(ind) for ind in constraint_x] ) - control_dict[ - "fix constraintx" - ] = "constraintx setforce 0.0 NULL NULL" + control_dict["fix constraintx"] = "constraintx setforce 0.0 NULL NULL" if calc_md: - control_dict[ - "velocity constraintx" - ] = "set 0.0 NULL NULL" + control_dict["velocity constraintx"] = "set 0.0 NULL NULL" if len(constraint_y) > 0: control_dict["group constrainty"] = "id " + " ".join( [str(ind) for ind in constraint_y] ) - control_dict[ - "fix constrainty" - ] = "constrainty setforce NULL 0.0 NULL" + control_dict["fix constrainty"] = "constrainty setforce NULL 0.0 NULL" if calc_md: - control_dict[ - "velocity constrainty" - ] = "set NULL 0.0 NULL" + control_dict["velocity constrainty"] = "set NULL 0.0 NULL" if len(constraint_z) > 0: control_dict["group constraintz"] = "id " + " ".join( [str(ind) for ind in constraint_z] ) - control_dict[ - "fix constraintz" - ] = "constraintz setforce NULL NULL 0.0" + control_dict["fix constraintz"] = "constraintz setforce NULL NULL 0.0" if calc_md: - control_dict[ - "velocity constraintz" - ] = "set NULL NULL 0.0" + control_dict["velocity constraintz"] = "set NULL NULL 0.0" return control_dict From bfef117d7409d48d3858ab88630704a62ca29410 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 12:58:41 -0700 Subject: [PATCH 03/22] fix bug --- pyiron_atomistics/lammps/wrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index bffb82e87..2a02398d7 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -18,7 +18,7 @@ def interactive_positions_getter(lmp, number_of_atoms, prism): def interactive_positions_setter(lmp, logger, positions, prism, cores, interactive): if _check_ortho_prism(prism=prism): positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, prism) + positions = np.matmul(positions, prism.R) positions = np.array(positions).flatten() if interactive and cores == 1: lmp.scatter_atoms( From 2790ee674f230b436426b2a43beed9fd551caa15 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 14:11:02 -0700 Subject: [PATCH 04/22] Introduce lammps wrapper --- pyiron_atomistics/lammps/interactive.py | 144 ++--- pyiron_atomistics/lammps/wrapper.py | 743 +++++++++++++----------- 2 files changed, 431 insertions(+), 456 deletions(-) diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index ef8ebac62..683f8967f 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -3,7 +3,6 @@ # Distributed under the terms of "New BSD License", see the LICENSE file. from ctypes import c_double, c_int -import importlib import numpy as np import os import pandas as pd @@ -14,22 +13,7 @@ from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.control import LammpsControl from pyiron_atomistics.atomistics.job.interactive import GenericInteractive -from pyiron_atomistics.lammps.wrapper import ( - interactive_lib_command, - interactive_positions_setter, - interactive_positions_getter, - interactive_cells_getter, - interactive_cells_setter, - interactive_volume_getter, - interactive_forces_getter, - interactive_structure_setter, - interactive_energy_pot_getter, - interactive_energy_tot_getter, - interactive_steps_getter, - interactive_temperatures_getter, - interactive_indices_getter, - interactive_pressures_getter, -) +from pyiron_atomistics.lammps.wrapper import PyironLammpsLibrary try: # mpi4py is only supported on Linux and Mac Os X @@ -91,66 +75,47 @@ def interactive_mpi_communicator(self, comm): self._interactive_mpi_communicator = comm def _interactive_lib_command(self, command): - interactive_lib_command( - lmp=self._interactive_library, logger=self._logger, command=command - ) + self._interactive_library.interactive_lib_command(command=command) def interactive_positions_getter(self): uc = UnitConverter(units=self.units) - positions = interactive_positions_getter( - lmp=self._interactive_library, - number_of_atoms=len(self.structure), - prism=self._prism, - ) + positions = self._interactive_library.interactive_positions_getter() positions = uc.convert_array_to_pyiron_units(positions, label="positions") return positions.tolist() def interactive_positions_setter(self, positions): - interactive_positions_setter( - lmp=self._interactive_library, - logger=self._logger, + self._interactive_library.interactive_positions_setter( positions=positions, - prism=self._prism, - cores=self.server.cores, - interactive=self.server.run_mode.interactive, ) def interactive_cells_getter(self): uc = UnitConverter(units=self.units) - cc = interactive_cells_getter(lmp=self._interactive_library) + cc = self._interactive_library.interactive_cells_getter() return uc.convert_array_to_pyiron_units( self._prism.unfold_cell(cc), label="cells" ) def interactive_cells_setter(self, cell): - self._prism = interactive_cells_setter( - lmp=self._interactive_library, - logger=self._logger, + self._interactive_library.interactive_cells_setter( cell=cell, - structure_current=self._structure_current, - structure_previous=self._structure_previous, + structure=self._structure_current ) def interactive_volume_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_volume_getter(lmp=self._interactive_library), label="volume" + self._interactive_library.interactive_volume_getter(), + label="volume" ) def interactive_forces_getter(self): uc = UnitConverter(units=self.units) - ff = interactive_forces_getter( - lmp=self._interactive_library, - prism=self._prism, - number_of_atoms=len(self.structure), - ) + ff = self._interactive_library.interactive_forces_getter() ff = uc.convert_array_to_pyiron_units(ff, label="forces") return ff.tolist() def interactive_execute(self): - interactive_lib_command( - lmp=self._interactive_library, - logger=self._logger, + self._interactive_library.interactive_lib_command( command=self._interactive_run_command, ) @@ -216,18 +181,13 @@ def _reset_interactive_run_command(self): def interactive_initialize_interface(self): self._create_working_directory() - if self.server.run_mode.interactive and self.server.cores == 1: - lammps = getattr(importlib.import_module("lammps"), "lammps") - if self._log_file is None: - self._log_file = os.path.join(self.working_directory, "log.lammps") - self._interactive_library = lammps( - cmdargs=["-screen", "none", "-log", self._log_file], - comm=self._interactive_mpi_communicator, - ) - else: - self._interactive_library = LammpsLibrary( - cores=self.server.cores, working_directory=self.working_directory - ) + self._interactive_library = PyironLammpsLibrary( + working_directory=self.working_directory, + cores=self.server.cores, + comm=self._interactive_mpi_communicator, + logger=self._logger, + log_file=self._log_file + ) if not all(self.structure.pbc): self.input.control["boundary"] = " ".join( ["p" if coord else "f" for coord in self.structure.pbc] @@ -410,7 +370,7 @@ def run_if_interactive(self): self._reset_interactive_run_command() if self._user_fix_external is not None: self._interactive_library.set_fix_external_callback( - "fix_external", self._user_fix_external.fix_external + fix_id="fix_external", callback=self._user_fix_external.fix_external, caller=None ) counter = 0 iteration_max = int( @@ -459,19 +419,14 @@ def interactive_fetch(self): self._logger.debug("interactive run - done") def interactive_structure_setter(self, structure): - self._prism, control_dict = interactive_structure_setter( - lmp=self._interactive_library, - logger=self._logger, - structure_current=structure, - structure_previous=self.structure, + self._prism, control_dict = self._interactive_library.interactive_structure_setter( + structure=structure, units=self.input.control["units"], dimension=self.input.control["dimension"], boundary=self.input.control["boundary"], atom_style=self.input.control["atom_style"], el_eam_lst=self.input.potential.get_element_lst(), calc_md=self._generic_input["calc_mode"] == "md", - interactive=self.server.run_mode.interactive, - cores=self.server.cores, ) for key, value in control_dict.items(): self.input.control[key.replace(" ", "___")] = value @@ -561,56 +516,43 @@ def update_potential(self): def interactive_indices_getter(self): uc = UnitConverter(units=self.units) - lammps_indices = interactive_indices_getter(lmp=self._interactive_library) + lammps_indices = self._interactive_library.interactive_indices_getter() indices = uc.convert_array_to_pyiron_units( self.remap_indices(lammps_indices), label="indices" ) return indices.tolist() def interactive_indices_setter(self, indices): - el_struct_lst = self._structure_current.get_species_symbols() - el_obj_lst = self._structure_current.get_species_objects() - el_eam_lst = self.input.potential.get_element_lst() - el_dict = {} - for id_eam, el_eam in enumerate(el_eam_lst): - if el_eam in el_struct_lst: - id_el = list(el_struct_lst).index(el_eam) - el = el_obj_lst[id_el] - el_dict[el] = id_eam + 1 - elem_all = np.array( - [el_dict[self._structure_current.species[el]] for el in indices] + self._interactive_library.interactive_indices_setter( + indices=indices, + el_eam_lst=self.input.potential.get_element_lst() ) - if self.server.run_mode.interactive and self.server.cores == 1: - self._interactive_library.scatter_atoms( - "type", 0, 1, (len(elem_all) * c_int)(*elem_all) - ) - else: - self._interactive_library.scatter_atoms("type", elem_all) def interactive_energy_pot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_energy_pot_getter(lmp=self._interactive_library), + self._interactive_library.interactive_energy_pot_getter(), label="energy_pot", ) def interactive_energy_tot_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_energy_tot_getter(lmp=self._interactive_library), + self._interactive_library.interactive_energy_tot_getter(), label="energy_tot", ) def interactive_steps_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_steps_getter(lmp=self._interactive_library), label="steps" + self._interactive_library.interactive_steps_getter(), + label="steps" ) def interactive_temperatures_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - interactive_temperatures_getter(lmp=self._interactive_library), + self._interactive_library.interactive_temperatures_getter(), label="temperature", ) @@ -623,33 +565,15 @@ def interactive_stress_getter(self): numpy.array: Nx3x3 np array of stress/atom """ if not "stress" in self.interactive_cache.keys(): - self._interactive_lib_command("compute st all stress/atom NULL") - self._interactive_lib_command("run 0") + ss = self._interactive_library.interactive_stress_getter(enable_stress_computation=True) self.interactive_cache["stress"] = [] - id_lst = self._interactive_library.extract_atom("id", 0) - id_lst = np.array([id_lst[i] for i in range(len(self.structure))]) - 1 - id_lst = np.arange(len(id_lst))[np.argsort(id_lst)] - ind = np.array([0, 3, 4, 3, 1, 5, 4, 5, 2]) - ss = self._interactive_library.extract_compute("st", 1, 2) - ss = np.array( - [ss[i][j] for i in range(len(self.structure)) for j in range(6)] - ).reshape(-1, 6)[id_lst] - ss = ( - ss[:, ind].reshape(len(self.structure), 3, 3) - / constants.eV - * constants.bar - * constants.angstrom**3 - ) - if _check_ortho_prism(prism=self._prism): - ss = np.einsum("ij,njk->nik", self._prism.R, ss) - ss = np.einsum("nij,kj->nik", ss, self._prism.R) + else: + ss = self._interactive_library.interactive_stress_getter(enable_stress_computation=False) return ss def interactive_pressures_getter(self): uc = UnitConverter(units=self.units) - pp = interactive_pressures_getter( - lmp=self._interactive_library, prism=self._prism - ) + pp = self._interactive_library.interactive_pressures_getter() return uc.convert_array_to_pyiron_units(pp, label="pressure") def interactive_close(self): diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 63ed69c9d..b70564087 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -1,375 +1,426 @@ from ctypes import c_double, c_int +import importlib import numpy as np +import os +from scipy import constants import warnings + from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.base import _check_ortho_prism +try: # mpi4py is only supported on Linux and Mac Os X + from pylammpsmpi import LammpsLibrary +except ImportError: + pass + + +class PyironLammpsLibrary(object): + def __int__(self, working_directory, cores=1, comm=None, logger=None, log_file=None): + self._logger = logger + self._prism = None + self._structure = None + self._cores = cores + if self._cores == 1: + lammps = getattr(importlib.import_module("lammps"), "lammps") + if log_file is None: + log_file = os.path.join(working_directory, "log.lammps") + self._interactive_library = lammps( + cmdargs=["-screen", "none", "-log", log_file], + comm=comm, + ) + else: + self._interactive_library = LammpsLibrary( + cores=self._cores, working_directory=working_directory + ) -def interactive_positions_getter(lmp, number_of_atoms, prism): - positions = np.reshape( - np.array(lmp.gather_atoms("x", 1, 3)), - (number_of_atoms, 3), - ) - if _check_ortho_prism(prism=prism): - positions = np.matmul(positions, prism.R.T) - return positions - - -def interactive_positions_setter(lmp, logger, positions, prism, cores, interactive): - if _check_ortho_prism(prism=prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, prism.R) - positions = np.array(positions).flatten() - if interactive and cores == 1: - lmp.scatter_atoms("x", 1, 3, (len(positions) * c_double)(*positions)) - else: - lmp.scatter_atoms("x", positions) - interactive_lib_command(lmp=lmp, logger=logger, command="change_box all remap") - - -def interactive_lib_command(lmp, logger, command): - logger.debug("Lammps library: " + command) - lmp.command(command) + def interactive_lib_command(self, command): + if self._logger is not None: + self._logger.debug("Lammps library: " + command) + self._interactive_library.command(command) + def interactive_positions_getter(self): + positions = np.reshape( + np.array(self._interactive_library.gather_atoms("x", 1, 3)), + (len(self._structure), 3), + ) + if _check_ortho_prism(prism=self._prism): + positions = np.matmul(positions, self._prism.R.T) + return positions + + def interactive_positions_setter(self, positions): + if _check_ortho_prism(prism=self._prism): + positions = np.array(positions).reshape(-1, 3) + positions = np.matmul(positions, self._prism.R) + positions = np.array(positions).flatten() + if self._cores == 1: + self._interactive_library.scatter_atoms("x", 1, 3, (len(positions) * c_double)(*positions)) + else: + self._interactive_library.scatter_atoms("x", positions) + self.interactive_lib_command(command="change_box all remap") -def interactive_cells_getter(lmp): - return np.array( - [ - [lmp.get_thermo("lx"), 0, 0], - [ - lmp.get_thermo("xy"), - lmp.get_thermo("ly"), - 0, - ], + def interactive_cells_getter(self): + return np.array( [ - lmp.get_thermo("xz"), - lmp.get_thermo("yz"), - lmp.get_thermo("lz"), - ], - ] - ) - - -def interactive_cells_setter(lmp, logger, cell, structure_current, structure_previous): - prism = UnfoldingPrism(cell) - lx, ly, lz, xy, xz, yz = prism.get_lammps_prism() - if _check_ortho_prism(prism=prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" + [self._interactive_library.get_thermo("lx"), 0, 0], + [ + self._interactive_library.get_thermo("xy"), + self._interactive_library.get_thermo("ly"), + 0, + ], + [ + self._interactive_library.get_thermo("xz"), + self._interactive_library.get_thermo("yz"), + self._interactive_library.get_thermo("lz"), + ], + ] ) - is_skewed = structure_current.is_skewed(tolerance=1.0e-8) - was_skewed = structure_previous.is_skewed(tolerance=1.0e-8) - - if is_skewed: - if not was_skewed: - interactive_lib_command( - lmp=lmp, logger=logger, command="change_box all triclinic" + def interactive_cells_setter(self, cell, structure): + self._prism = UnfoldingPrism(cell) + lx, ly, lz, xy, xz, yz = self._prism.get_lammps_prism() + if _check_ortho_prism(prism=self._prism): + warnings.warn( + "Warning: setting upper trangular matrix might slow down the calculation" ) - interactive_lib_command( - lmp=lmp, - logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, xy, xz, yz), - ) - elif was_skewed: - interactive_lib_command( - lmp=lmp, - logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, 0.0, 0.0, 0.0), - ) - interactive_lib_command(lmp=lmp, logger=logger, command="change_box all ortho") - else: - interactive_lib_command( - lmp=lmp, - logger=logger, - command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" - % (lx, ly, lz), - ) - return prism - -def interactive_volume_getter(lmp): - return lmp.get_thermo("vol") + is_skewed = structure.is_skewed(tolerance=1.0e-8) + was_skewed = self._structure.is_skewed(tolerance=1.0e-8) - -def interactive_forces_getter(lmp, prism, number_of_atoms): - ff = np.reshape( - np.array(lmp.gather_atoms("f", 1, 3)), - (number_of_atoms, 3), - ) - if _check_ortho_prism(prism=prism): - ff = np.matmul(ff, prism.R.T) - return ff - - -def interactive_structure_setter( - lmp, - logger, - structure_current, - structure_previous, - units, - dimension, - boundary, - atom_style, - el_eam_lst, - calc_md, - interactive, - cores, -): - old_symbols = structure_previous.get_species_symbols() - new_symbols = structure_current.get_species_symbols() - if any(old_symbols != new_symbols): - raise ValueError( - f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" - ) - interactive_lib_command(lmp=lmp, logger=logger, command="clear") - control_dict = set_selective_dynamics(structure=structure_current, calc_md=calc_md) - interactive_lib_command(lmp=lmp, logger=logger, command="units " + units) - interactive_lib_command( - lmp=lmp, logger=logger, command="dimension " + str(dimension) - ) - interactive_lib_command(lmp=lmp, logger=logger, command="boundary " + boundary) - interactive_lib_command(lmp=lmp, logger=logger, command="atom_style " + atom_style) - - interactive_lib_command(lmp=lmp, logger=logger, command="atom_modify map array") - prism = UnfoldingPrism(structure_current.cell) - if _check_ortho_prism(prism=prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" - ) - xhi, yhi, zhi, xy, xz, yz = prism.get_lammps_prism() - if prism.is_skewed(): - interactive_lib_command( - lmp=lmp, - logger=logger, - command="region 1 prism" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " " - + str(xy) - + " " - + str(xz) - + " " - + str(yz) - + " units box", - ) - else: - interactive_lib_command( - lmp=lmp, - logger=logger, - command="region 1 block" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " units box", - ) - el_struct_lst = structure_current.get_species_symbols() - el_obj_lst = structure_current.get_species_objects() - if atom_style == "full": - interactive_lib_command( - lmp=lmp, - logger=logger, - command="create_box " - + str(len(el_eam_lst)) - + " 1 " - + "bond/types 1 " - + "angle/types 1 " - + "extra/bond/per/atom 2 " - + "extra/angle/per/atom 2 ", - ) - else: - interactive_lib_command( - lmp=lmp, logger=logger, command="create_box " + str(len(el_eam_lst)) + " 1" - ) - el_dict = {} - for id_eam, el_eam in enumerate(el_eam_lst): - if el_eam in el_struct_lst: - id_el = list(el_struct_lst).index(el_eam) - el = el_obj_lst[id_el] - el_dict[el] = id_eam + 1 - interactive_lib_command( - lmp=lmp, - logger=logger, - command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass), + if is_skewed: + if not was_skewed: + self.interactive_lib_command(command="change_box all triclinic") + self.interactive_lib_command( + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" + % (lx, ly, lz, xy, xz, yz), + ) + elif was_skewed: + self.interactive_lib_command( + command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" + % (lx, ly, lz, 0.0, 0.0, 0.0), ) + self.interactive_lib_command(command="change_box all ortho") else: - interactive_lib_command( - lmp=lmp, - logger=logger, - command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00), + self.interactive_lib_command( + command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" + % (lx, ly, lz), ) - positions = structure_current.positions.flatten() - if _check_ortho_prism(prism=prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, prism.R) - positions = positions.flatten() - try: - elem_all = np.array( - [el_dict[el] for el in structure_current.get_chemical_elements()] - ) - except KeyError: - missing = set(structure_current.get_chemical_elements()).difference( - el_dict.keys() - ) - missing = ", ".join([el.Abbreviation for el in missing]) - raise ValueError( - f"Structure contains elements [{missing}], that are not present in the potential!" - ) - if interactive and cores == 1: - lmp.create_atoms( - n=len(structure_current), - id=None, - type=(len(elem_all) * c_int)(*elem_all), - x=(len(positions) * c_double)(*positions), - v=None, - image=None, - shrinkexceed=False, - ) - else: - lmp.create_atoms( - n=len(structure_current), - id=None, - type=elem_all, - x=positions, - v=None, - image=None, - shrinkexceed=False, - ) - interactive_lib_command(lmp=lmp, logger=logger, command="change_box all remap") - return prism, control_dict + def interactive_volume_getter(self): + return self._interactive_library.get_thermo("vol") -def set_selective_dynamics(structure, calc_md): - control_dict = {} - if "selective_dynamics" in structure._tag_list.keys(): - if structure.selective_dynamics._default is None: - structure.selective_dynamics._default = [True, True, True] - sel_dyn = np.logical_not(structure.selective_dynamics.list()) - # Enter loop only if constraints present - if len(np.argwhere(np.any(sel_dyn, axis=1)).flatten()) != 0: - all_indices = np.arange(len(structure), dtype=int) - constraint_xyz = np.argwhere(np.all(sel_dyn, axis=1)).flatten() - not_constrained_xyz = np.setdiff1d(all_indices, constraint_xyz) - # LAMMPS starts counting from 1 - constraint_xyz += 1 - ind_x = np.argwhere(sel_dyn[not_constrained_xyz, 0]).flatten() - ind_y = np.argwhere(sel_dyn[not_constrained_xyz, 1]).flatten() - ind_z = np.argwhere(sel_dyn[not_constrained_xyz, 2]).flatten() - constraint_xy = not_constrained_xyz[np.intersect1d(ind_x, ind_y)] + 1 - constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 - constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 - constraint_x = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + 1 + def interactive_forces_getter(self): + ff = np.reshape( + np.array(self._interactive_library.gather_atoms("f", 1, 3)), + (len(self._structure), 3), + ) + if _check_ortho_prism(prism=self._prism): + ff = np.matmul(ff, self._prism.R.T) + return ff + + def interactive_structure_setter( + self, + structure, + units, + dimension, + boundary, + atom_style, + el_eam_lst, + calc_md=True, + ): + old_symbols = self._structure.get_species_symbols() + new_symbols = structure.get_species_symbols() + if any(old_symbols != new_symbols): + raise ValueError( + f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" ) - constraint_y = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + 1 + self.interactive_lib_command(command="clear") + control_dict = self._set_selective_dynamics(structure=structure, calc_md=calc_md) + self.interactive_lib_command(command="units " + units) + self.interactive_lib_command(command="dimension " + str(dimension)) + self.interactive_lib_command(command="boundary " + boundary) + self.interactive_lib_command(command="atom_style " + atom_style) + + self.interactive_lib_command(command="atom_modify map array") + self._prism = UnfoldingPrism(structure.cell) + if _check_ortho_prism(prism=self._prism): + warnings.warn( + "Warning: setting upper trangular matrix might slow down the calculation" ) - constraint_z = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + 1 + xhi, yhi, zhi, xy, xz, yz = self._prism.get_lammps_prism() + if self._prism.is_skewed(): + self.interactive_lib_command( + command="region 1 prism" + + " 0.0 " + + str(xhi) + + " 0.0 " + + str(yhi) + + " 0.0 " + + str(zhi) + + " " + + str(xy) + + " " + + str(xz) + + " " + + str(yz) + + " units box", ) - control_dict = {} - if len(constraint_xyz) > 0: - control_dict["group constraintxyz"] = "id " + " ".join( - [str(ind) for ind in constraint_xyz] - ) - control_dict["fix constraintxyz"] = "constraintxyz setforce 0.0 0.0 0.0" - if calc_md: - control_dict["velocity constraintxyz"] = "set 0.0 0.0 0.0" - if len(constraint_xy) > 0: - control_dict["group constraintxy"] = "id " + " ".join( - [str(ind) for ind in constraint_xy] - ) - control_dict["fix constraintxy"] = "constraintxy setforce 0.0 0.0 NULL" - if calc_md: - control_dict["velocity constraintxy"] = "set 0.0 0.0 NULL" - if len(constraint_yz) > 0: - control_dict["group constraintyz"] = "id " + " ".join( - [str(ind) for ind in constraint_yz] + else: + self.interactive_lib_command( + command="region 1 block" + + " 0.0 " + + str(xhi) + + " 0.0 " + + str(yhi) + + " 0.0 " + + str(zhi) + + " units box", + ) + el_struct_lst = structure.get_species_symbols() + el_obj_lst = structure.get_species_objects() + if atom_style == "full": + self.interactive_lib_command( + command="create_box " + + str(len(el_eam_lst)) + + " 1 " + + "bond/types 1 " + + "angle/types 1 " + + "extra/bond/per/atom 2 " + + "extra/angle/per/atom 2 ", + ) + else: + self.interactive_lib_command( + command="create_box " + str(len(el_eam_lst)) + " 1" + ) + el_dict = {} + for id_eam, el_eam in enumerate(el_eam_lst): + if el_eam in el_struct_lst: + id_el = list(el_struct_lst).index(el_eam) + el = el_obj_lst[id_el] + el_dict[el] = id_eam + 1 + self.interactive_lib_command( + command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass), ) - control_dict["fix constraintyz"] = "constraintyz setforce NULL 0.0 0.0" - if calc_md: - control_dict["velocity constraintyz"] = "set NULL 0.0 0.0" - if len(constraint_zx) > 0: - control_dict["group constraintxz"] = "id " + " ".join( - [str(ind) for ind in constraint_zx] + else: + self.interactive_lib_command( + command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00), ) - control_dict["fix constraintxz"] = "constraintxz setforce 0.0 NULL 0.0" - if calc_md: - control_dict["velocity constraintxz"] = "set 0.0 NULL 0.0" - if len(constraint_x) > 0: - control_dict["group constraintx"] = "id " + " ".join( - [str(ind) for ind in constraint_x] + positions = structure.positions.flatten() + if _check_ortho_prism(prism=self._prism): + positions = np.array(positions).reshape(-1, 3) + positions = np.matmul(positions, self._prism.R) + positions = positions.flatten() + try: + elem_all = np.array( + [el_dict[el] for el in structure.get_chemical_elements()] + ) + except KeyError: + missing = set(structure.get_chemical_elements()).difference( + el_dict.keys() + ) + missing = ", ".join([el.Abbreviation for el in missing]) + raise ValueError( + f"Structure contains elements [{missing}], that are not present in the potential!" + ) + if self._cores == 1: + self._interactive_library.create_atoms( + n=len(structure), + id=None, + type=(len(elem_all) * c_int)(*elem_all), + x=(len(positions) * c_double)(*positions), + v=None, + image=None, + shrinkexceed=False, + ) + else: + self._interactive_library.create_atoms( + n=len(structure), + id=None, + type=elem_all, + x=positions, + v=None, + image=None, + shrinkexceed=False, + ) + self.interactive_lib_command(command="change_box all remap") + for key, value in control_dict.items(): + self.interactive_lib_command(command=key + " " + value) + self._structure = structure + + @staticmethod + def _set_selective_dynamics(structure, calc_md): + control_dict = {} + if "selective_dynamics" in structure._tag_list.keys(): + if structure.selective_dynamics._default is None: + structure.selective_dynamics._default = [True, True, True] + sel_dyn = np.logical_not(structure.selective_dynamics.list()) + # Enter loop only if constraints present + if len(np.argwhere(np.any(sel_dyn, axis=1)).flatten()) != 0: + all_indices = np.arange(len(structure), dtype=int) + constraint_xyz = np.argwhere(np.all(sel_dyn, axis=1)).flatten() + not_constrained_xyz = np.setdiff1d(all_indices, constraint_xyz) + # LAMMPS starts counting from 1 + constraint_xyz += 1 + ind_x = np.argwhere(sel_dyn[not_constrained_xyz, 0]).flatten() + ind_y = np.argwhere(sel_dyn[not_constrained_xyz, 1]).flatten() + ind_z = np.argwhere(sel_dyn[not_constrained_xyz, 2]).flatten() + constraint_xy = not_constrained_xyz[np.intersect1d(ind_x, ind_y)] + 1 + constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 + constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 + constraint_x = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + 1 ) - control_dict["fix constraintx"] = "constraintx setforce 0.0 NULL NULL" - if calc_md: - control_dict["velocity constraintx"] = "set 0.0 NULL NULL" - if len(constraint_y) > 0: - control_dict["group constrainty"] = "id " + " ".join( - [str(ind) for ind in constraint_y] + constraint_y = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + 1 ) - control_dict["fix constrainty"] = "constrainty setforce NULL 0.0 NULL" - if calc_md: - control_dict["velocity constrainty"] = "set NULL 0.0 NULL" - if len(constraint_z) > 0: - control_dict["group constraintz"] = "id " + " ".join( - [str(ind) for ind in constraint_z] + constraint_z = ( + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + 1 ) - control_dict["fix constraintz"] = "constraintz setforce NULL NULL 0.0" - if calc_md: - control_dict["velocity constraintz"] = "set NULL NULL 0.0" - return control_dict - - -def interactive_indices_getter(lmp): - return np.array(lmp.gather_atoms("type", 0, 1)) - - -def interactive_energy_pot_getter(lmp): - return lmp.get_thermo("pe") - - -def interactive_energy_tot_getter(lmp): - return lmp.get_thermo("etotal") - - -def interactive_steps_getter(lmp): - return lmp.get_thermo("step") - - -def interactive_temperatures_getter(lmp): - return lmp.get_thermo("temp") - - -def interactive_pressures_getter(lmp, prism): - pp = np.array( - [ - [ - lmp.get_thermo("pxx"), - lmp.get_thermo("pxy"), - lmp.get_thermo("pxz"), - ], + control_dict = {} + if len(constraint_xyz) > 0: + control_dict["group constraintxyz"] = "id " + " ".join( + [str(ind) for ind in constraint_xyz] + ) + control_dict["fix constraintxyz"] = "constraintxyz setforce 0.0 0.0 0.0" + if calc_md: + control_dict["velocity constraintxyz"] = "set 0.0 0.0 0.0" + if len(constraint_xy) > 0: + control_dict["group constraintxy"] = "id " + " ".join( + [str(ind) for ind in constraint_xy] + ) + control_dict["fix constraintxy"] = "constraintxy setforce 0.0 0.0 NULL" + if calc_md: + control_dict["velocity constraintxy"] = "set 0.0 0.0 NULL" + if len(constraint_yz) > 0: + control_dict["group constraintyz"] = "id " + " ".join( + [str(ind) for ind in constraint_yz] + ) + control_dict["fix constraintyz"] = "constraintyz setforce NULL 0.0 0.0" + if calc_md: + control_dict["velocity constraintyz"] = "set NULL 0.0 0.0" + if len(constraint_zx) > 0: + control_dict["group constraintxz"] = "id " + " ".join( + [str(ind) for ind in constraint_zx] + ) + control_dict["fix constraintxz"] = "constraintxz setforce 0.0 NULL 0.0" + if calc_md: + control_dict["velocity constraintxz"] = "set 0.0 NULL 0.0" + if len(constraint_x) > 0: + control_dict["group constraintx"] = "id " + " ".join( + [str(ind) for ind in constraint_x] + ) + control_dict["fix constraintx"] = "constraintx setforce 0.0 NULL NULL" + if calc_md: + control_dict["velocity constraintx"] = "set 0.0 NULL NULL" + if len(constraint_y) > 0: + control_dict["group constrainty"] = "id " + " ".join( + [str(ind) for ind in constraint_y] + ) + control_dict["fix constrainty"] = "constrainty setforce NULL 0.0 NULL" + if calc_md: + control_dict["velocity constrainty"] = "set NULL 0.0 NULL" + if len(constraint_z) > 0: + control_dict["group constraintz"] = "id " + " ".join( + [str(ind) for ind in constraint_z] + ) + control_dict["fix constraintz"] = "constraintz setforce NULL NULL 0.0" + if calc_md: + control_dict["velocity constraintz"] = "set NULL NULL 0.0" + return control_dict + + def interactive_indices_getter(self): + return np.array(self._interactive_library.gather_atoms("type", 0, 1)) + + def interactive_energy_pot_getter(self): + return self._interactive_library.get_thermo("pe") + + def interactive_energy_tot_getter(self): + return self._interactive_library.get_thermo("etotal") + + def interactive_steps_getter(self): + return self._interactive_library.get_thermo("step") + + def interactive_temperatures_getter(self): + return self._interactive_library.get_thermo("temp") + + def interactive_pressures_getter(self): + pp = np.array( [ - lmp.get_thermo("pxy"), - lmp.get_thermo("pyy"), - lmp.get_thermo("pyz"), - ], - [ - lmp.get_thermo("pxz"), - lmp.get_thermo("pyz"), - lmp.get_thermo("pzz"), - ], - ] - ) - if _check_ortho_prism(prism=prism): - rotation_matrix = prism.R.T - pp = rotation_matrix.T @ pp @ rotation_matrix - return pp + [ + self._interactive_library.get_thermo("pxx"), + self._interactive_library.get_thermo("pxy"), + self._interactive_library.get_thermo("pxz"), + ], + [ + self._interactive_library.get_thermo("pxy"), + self._interactive_library.get_thermo("pyy"), + self._interactive_library.get_thermo("pyz"), + ], + [ + self._interactive_library.get_thermo("pxz"), + self._interactive_library.get_thermo("pyz"), + self._interactive_library.get_thermo("pzz"), + ], + ] + ) + if _check_ortho_prism(prism=self._prism): + rotation_matrix = self._prism.R.T + pp = rotation_matrix.T @ pp @ rotation_matrix + return pp + + def interactive_indices_setter(self, indices, el_eam_lst): + el_struct_lst = self._structure.get_species_symbols() + el_obj_lst = self._structure.get_species_objects() + el_dict = {} + for id_eam, el_eam in enumerate(el_eam_lst): + if el_eam in el_struct_lst: + id_el = list(el_struct_lst).index(el_eam) + el = el_obj_lst[id_el] + el_dict[el] = id_eam + 1 + elem_all = np.array( + [el_dict[self._structure.species[el]] for el in indices] + ) + if self._cores == 1: + self._interactive_library.scatter_atoms( + "type", 0, 1, (len(elem_all) * c_int)(*elem_all) + ) + else: + self._interactive_library.scatter_atoms("type", elem_all) + + def interactive_stress_getter(self, enable_stress_computation=True): + """ + This gives back an Nx3x3 array of stress/atom defined in http://lammps.sandia.gov/doc/compute_stress_atom.html + Keep in mind that it is stress*volume in eV. Further discussion can be found on the website above. + + Returns: + numpy.array: Nx3x3 np array of stress/atom + """ + if enable_stress_computation: + self.interactive_lib_command("compute st all stress/atom NULL") + self.interactive_lib_command("run 0") + id_lst = self._interactive_library.extract_atom("id", 0) + id_lst = np.array([id_lst[i] for i in range(len(self._structure))]) - 1 + id_lst = np.arange(len(id_lst))[np.argsort(id_lst)] + ind = np.array([0, 3, 4, 3, 1, 5, 4, 5, 2]) + ss = self._interactive_library.extract_compute("st", 1, 2) + ss = np.array( + [ss[i][j] for i in range(len(self._structure)) for j in range(6)] + ).reshape(-1, 6)[id_lst] + ss = ( + ss[:, ind].reshape(len(self._structure), 3, 3) + / constants.eV + * constants.bar + * constants.angstrom**3 + ) + if _check_ortho_prism(prism=self._prism): + ss = np.einsum("ij,njk->nik", self._prism.R, ss) + ss = np.einsum("nij,kj->nik", ss, self._prism.R) + return ss + + def close(self): + if self._interactive_library is not None: + self._interactive_library.close() + + def set_fix_external_callback(self, fix_id, callback, caller=None): + self._interactive_library.set_fix_external_callback( + fix_id=fix_id, callback=callback, caller=caller + ) From 82ad2ef56ca20a46af3e54a60547cda3a867a8ea Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 14:13:03 -0700 Subject: [PATCH 05/22] remove unused imports --- pyiron_atomistics/lammps/interactive.py | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 683f8967f..c03219aaa 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -2,24 +2,16 @@ # Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department # Distributed under the terms of "New BSD License", see the LICENSE file. -from ctypes import c_double, c_int import numpy as np import os import pandas as pd import warnings -from scipy import constants -from pyiron_atomistics.lammps.base import LammpsBase, _check_ortho_prism +from pyiron_atomistics.lammps.base import LammpsBase from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.control import LammpsControl from pyiron_atomistics.atomistics.job.interactive import GenericInteractive from pyiron_atomistics.lammps.wrapper import PyironLammpsLibrary - - -try: # mpi4py is only supported on Linux and Mac Os X - from pylammpsmpi import LammpsLibrary -except ImportError: - pass from pyiron_atomistics.lammps.units import UnitConverter __author__ = "Osamu Waseda, Jan Janssen" From ad04434b88e68f4ac71cc9d4aded7e0f413712fd Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 14:24:18 -0700 Subject: [PATCH 06/22] minor bugfixes --- pyiron_atomistics/lammps/interactive.py | 9 +++----- pyiron_atomistics/lammps/wrapper.py | 29 +++++++++++++++++++++---- 2 files changed, 28 insertions(+), 10 deletions(-) diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index c03219aaa..20f8b60f3 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -82,15 +82,14 @@ def interactive_positions_setter(self, positions): def interactive_cells_getter(self): uc = UnitConverter(units=self.units) - cc = self._interactive_library.interactive_cells_getter() return uc.convert_array_to_pyiron_units( - self._prism.unfold_cell(cc), label="cells" + self._interactive_library.interactive_cells_getter(), + label="cells" ) def interactive_cells_setter(self, cell): self._interactive_library.interactive_cells_setter( cell=cell, - structure=self._structure_current ) def interactive_volume_getter(self): @@ -411,7 +410,7 @@ def interactive_fetch(self): self._logger.debug("interactive run - done") def interactive_structure_setter(self, structure): - self._prism, control_dict = self._interactive_library.interactive_structure_setter( + self._interactive_library.interactive_structure_setter( structure=structure, units=self.input.control["units"], dimension=self.input.control["dimension"], @@ -420,8 +419,6 @@ def interactive_structure_setter(self, structure): el_eam_lst=self.input.potential.get_element_lst(), calc_md=self._generic_input["calc_mode"] == "md", ) - for key, value in control_dict.items(): - self.input.control[key.replace(" ", "___")] = value self._interactive_lammps_input() self._interactive_set_potential() diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index b70564087..2ce5b3ac4 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -59,7 +59,7 @@ def interactive_positions_setter(self, positions): self.interactive_lib_command(command="change_box all remap") def interactive_cells_getter(self): - return np.array( + cc = np.array( [ [self._interactive_library.get_thermo("lx"), 0, 0], [ @@ -74,8 +74,9 @@ def interactive_cells_getter(self): ], ] ) + return self._prism.unfold_cell(cc) - def interactive_cells_setter(self, cell, structure): + def interactive_cells_setter(self, cell): self._prism = UnfoldingPrism(cell) lx, ly, lz, xy, xz, yz = self._prism.get_lammps_prism() if _check_ortho_prism(prism=self._prism): @@ -83,8 +84,8 @@ def interactive_cells_setter(self, cell, structure): "Warning: setting upper trangular matrix might slow down the calculation" ) - is_skewed = structure.is_skewed(tolerance=1.0e-8) - was_skewed = self._structure.is_skewed(tolerance=1.0e-8) + is_skewed = cell_is_skewed(cell=cell, tolerance=1.0e-8) + was_skewed = cell_is_skewed(cell=self._structure.cell, tolerance=1.0e-8) if is_skewed: if not was_skewed: @@ -424,3 +425,23 @@ def set_fix_external_callback(self, fix_id, callback, caller=None): self._interactive_library.set_fix_external_callback( fix_id=fix_id, callback=callback, caller=caller ) + + +def cell_is_skewed(cell, 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 = np.abs(np.linalg.det(cell)) + prod = np.linalg.norm(cell, axis=-1).prod() + if volume > 0: + if abs(volume - prod) / volume < tolerance: + return False + return True From a6224f7bd60844051581b70935e0fcfdf733a727 Mon Sep 17 00:00:00 2001 From: pyiron-runner Date: Mon, 27 Feb 2023 21:26:52 +0000 Subject: [PATCH 07/22] Format black --- pyiron_atomistics/lammps/interactive.py | 26 +++++------ pyiron_atomistics/lammps/wrapper.py | 57 ++++++++++++++++--------- 2 files changed, 52 insertions(+), 31 deletions(-) diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 20f8b60f3..1f236d929 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -83,8 +83,7 @@ def interactive_positions_setter(self, positions): def interactive_cells_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.interactive_cells_getter(), - label="cells" + self._interactive_library.interactive_cells_getter(), label="cells" ) def interactive_cells_setter(self, cell): @@ -95,8 +94,7 @@ def interactive_cells_setter(self, cell): def interactive_volume_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.interactive_volume_getter(), - label="volume" + self._interactive_library.interactive_volume_getter(), label="volume" ) def interactive_forces_getter(self): @@ -177,7 +175,7 @@ def interactive_initialize_interface(self): cores=self.server.cores, comm=self._interactive_mpi_communicator, logger=self._logger, - log_file=self._log_file + log_file=self._log_file, ) if not all(self.structure.pbc): self.input.control["boundary"] = " ".join( @@ -361,7 +359,9 @@ def run_if_interactive(self): self._reset_interactive_run_command() if self._user_fix_external is not None: self._interactive_library.set_fix_external_callback( - fix_id="fix_external", callback=self._user_fix_external.fix_external, caller=None + fix_id="fix_external", + callback=self._user_fix_external.fix_external, + caller=None, ) counter = 0 iteration_max = int( @@ -513,8 +513,7 @@ def interactive_indices_getter(self): def interactive_indices_setter(self, indices): self._interactive_library.interactive_indices_setter( - indices=indices, - el_eam_lst=self.input.potential.get_element_lst() + indices=indices, el_eam_lst=self.input.potential.get_element_lst() ) def interactive_energy_pot_getter(self): @@ -534,8 +533,7 @@ def interactive_energy_tot_getter(self): def interactive_steps_getter(self): uc = UnitConverter(units=self.units) return uc.convert_array_to_pyiron_units( - self._interactive_library.interactive_steps_getter(), - label="steps" + self._interactive_library.interactive_steps_getter(), label="steps" ) def interactive_temperatures_getter(self): @@ -554,10 +552,14 @@ def interactive_stress_getter(self): numpy.array: Nx3x3 np array of stress/atom """ if not "stress" in self.interactive_cache.keys(): - ss = self._interactive_library.interactive_stress_getter(enable_stress_computation=True) + ss = self._interactive_library.interactive_stress_getter( + enable_stress_computation=True + ) self.interactive_cache["stress"] = [] else: - ss = self._interactive_library.interactive_stress_getter(enable_stress_computation=False) + ss = self._interactive_library.interactive_stress_getter( + enable_stress_computation=False + ) return ss def interactive_pressures_getter(self): diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 2ce5b3ac4..e463dbb7a 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -15,7 +15,9 @@ class PyironLammpsLibrary(object): - def __int__(self, working_directory, cores=1, comm=None, logger=None, log_file=None): + def __int__( + self, working_directory, cores=1, comm=None, logger=None, log_file=None + ): self._logger = logger self._prism = None self._structure = None @@ -53,7 +55,9 @@ def interactive_positions_setter(self, positions): positions = np.matmul(positions, self._prism.R) positions = np.array(positions).flatten() if self._cores == 1: - self._interactive_library.scatter_atoms("x", 1, 3, (len(positions) * c_double)(*positions)) + self._interactive_library.scatter_atoms( + "x", 1, 3, (len(positions) * c_double)(*positions) + ) else: self._interactive_library.scatter_atoms("x", positions) self.interactive_lib_command(command="change_box all remap") @@ -135,7 +139,9 @@ def interactive_structure_setter( f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" ) self.interactive_lib_command(command="clear") - control_dict = self._set_selective_dynamics(structure=structure, calc_md=calc_md) + control_dict = self._set_selective_dynamics( + structure=structure, calc_md=calc_md + ) self.interactive_lib_command(command="units " + units) self.interactive_lib_command(command="dimension " + str(dimension)) self.interactive_lib_command(command="boundary " + boundary) @@ -215,9 +221,7 @@ def interactive_structure_setter( [el_dict[el] for el in structure.get_chemical_elements()] ) except KeyError: - missing = set(structure.get_chemical_elements()).difference( - el_dict.keys() - ) + missing = set(structure.get_chemical_elements()).difference(el_dict.keys()) missing = ", ".join([el.Abbreviation for el in missing]) raise ValueError( f"Structure contains elements [{missing}], that are not present in the potential!" @@ -268,62 +272,79 @@ def _set_selective_dynamics(structure, calc_md): constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 constraint_x = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] + + 1 ) constraint_y = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] + + 1 ) constraint_z = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + 1 + not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] + + 1 ) control_dict = {} if len(constraint_xyz) > 0: control_dict["group constraintxyz"] = "id " + " ".join( [str(ind) for ind in constraint_xyz] ) - control_dict["fix constraintxyz"] = "constraintxyz setforce 0.0 0.0 0.0" + control_dict[ + "fix constraintxyz" + ] = "constraintxyz setforce 0.0 0.0 0.0" if calc_md: control_dict["velocity constraintxyz"] = "set 0.0 0.0 0.0" if len(constraint_xy) > 0: control_dict["group constraintxy"] = "id " + " ".join( [str(ind) for ind in constraint_xy] ) - control_dict["fix constraintxy"] = "constraintxy setforce 0.0 0.0 NULL" + control_dict[ + "fix constraintxy" + ] = "constraintxy setforce 0.0 0.0 NULL" if calc_md: control_dict["velocity constraintxy"] = "set 0.0 0.0 NULL" if len(constraint_yz) > 0: control_dict["group constraintyz"] = "id " + " ".join( [str(ind) for ind in constraint_yz] ) - control_dict["fix constraintyz"] = "constraintyz setforce NULL 0.0 0.0" + control_dict[ + "fix constraintyz" + ] = "constraintyz setforce NULL 0.0 0.0" if calc_md: control_dict["velocity constraintyz"] = "set NULL 0.0 0.0" if len(constraint_zx) > 0: control_dict["group constraintxz"] = "id " + " ".join( [str(ind) for ind in constraint_zx] ) - control_dict["fix constraintxz"] = "constraintxz setforce 0.0 NULL 0.0" + control_dict[ + "fix constraintxz" + ] = "constraintxz setforce 0.0 NULL 0.0" if calc_md: control_dict["velocity constraintxz"] = "set 0.0 NULL 0.0" if len(constraint_x) > 0: control_dict["group constraintx"] = "id " + " ".join( [str(ind) for ind in constraint_x] ) - control_dict["fix constraintx"] = "constraintx setforce 0.0 NULL NULL" + control_dict[ + "fix constraintx" + ] = "constraintx setforce 0.0 NULL NULL" if calc_md: control_dict["velocity constraintx"] = "set 0.0 NULL NULL" if len(constraint_y) > 0: control_dict["group constrainty"] = "id " + " ".join( [str(ind) for ind in constraint_y] ) - control_dict["fix constrainty"] = "constrainty setforce NULL 0.0 NULL" + control_dict[ + "fix constrainty" + ] = "constrainty setforce NULL 0.0 NULL" if calc_md: control_dict["velocity constrainty"] = "set NULL 0.0 NULL" if len(constraint_z) > 0: control_dict["group constraintz"] = "id " + " ".join( [str(ind) for ind in constraint_z] ) - control_dict["fix constraintz"] = "constraintz setforce NULL NULL 0.0" + control_dict[ + "fix constraintz" + ] = "constraintz setforce NULL NULL 0.0" if calc_md: control_dict["velocity constraintz"] = "set NULL NULL 0.0" return control_dict @@ -377,9 +398,7 @@ def interactive_indices_setter(self, indices, el_eam_lst): id_el = list(el_struct_lst).index(el_eam) el = el_obj_lst[id_el] el_dict[el] = id_eam + 1 - elem_all = np.array( - [el_dict[self._structure.species[el]] for el in indices] - ) + elem_all = np.array([el_dict[self._structure.species[el]] for el in indices]) if self._cores == 1: self._interactive_library.scatter_atoms( "type", 0, 1, (len(elem_all) * c_int)(*elem_all) From 9b030e64598fab1c449a13175b3a54f7aa0fd344 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 15:36:45 -0700 Subject: [PATCH 08/22] fix tests --- pyiron_atomistics/lammps/wrapper.py | 6 ++- tests/lammps/test_interactive.py | 71 +++++++++++++++++++---------- 2 files changed, 50 insertions(+), 27 deletions(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 2ce5b3ac4..aaf46d952 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -15,12 +15,14 @@ class PyironLammpsLibrary(object): - def __int__(self, working_directory, cores=1, comm=None, logger=None, log_file=None): + def __init__(self, working_directory, cores=1, comm=None, logger=None, log_file=None, library=None): self._logger = logger self._prism = None self._structure = None self._cores = cores - if self._cores == 1: + if library is not None: + self._interactive_library = library + elif self._cores == 1: lammps = getattr(importlib.import_module("lammps"), "lammps") if log_file is None: log_file = os.path.join(working_directory, "log.lammps") diff --git a/tests/lammps/test_interactive.py b/tests/lammps/test_interactive.py index ac33d5a6e..f76626054 100644 --- a/tests/lammps/test_interactive.py +++ b/tests/lammps/test_interactive.py @@ -8,6 +8,8 @@ from pyiron_base import Project, ProjectHDFio from pyiron_atomistics.atomistics.structure.atoms import Atoms from pyiron_atomistics.lammps.lammps import Lammps +from pyiron_atomistics.lammps.wrapper import PyironLammpsLibrary, cell_is_skewed +from pyiron_atomistics.lammps.structure import UnfoldingPrism class InteractiveLibrary(object): @@ -23,9 +25,30 @@ def scatter_atoms(self, *args): class TestLammpsInteractive(unittest.TestCase): def setUp(self): - self.job._interactive_library = InteractiveLibrary() - self.minimize_job._interactive_library = InteractiveLibrary() - self.minimize_control_job._interactive_library = InteractiveLibrary() + self.job._interactive_library = PyironLammpsLibrary( + working_directory=self.job.working_directory, + cores=1, + comm=None, + logger=None, + log_file=None, + library=InteractiveLibrary() + ) + self.minimize_job._interactive_library = PyironLammpsLibrary( + working_directory=self.job.working_directory, + cores=1, + comm=None, + logger=None, + log_file=None, + library=InteractiveLibrary() + ) + self.minimize_control_job._interactive_library = PyironLammpsLibrary( + working_directory=self.job.working_directory, + cores=1, + comm=None, + logger=None, + log_file=None, + library=InteractiveLibrary() + ) @classmethod def setUpClass(cls): @@ -64,41 +87,39 @@ def tearDownClass(cls): def test_interactive_cells_setter(self): atoms = Atoms("Fe8", positions=np.zeros((8, 3)), cell=np.eye(3), pbc=True) - self.job._structure_previous = atoms.copy() - self.job._structure_current = atoms.copy() - self.job.interactive_cells_setter(self.job._structure_current.cell) + self.job._interactive_library._structure = atoms.copy() + self.job.interactive_cells_setter(atoms.cell) self.assertEqual( - self.job._interactive_library._command[-1], + self.job._interactive_library._interactive_library._command[-1], "change_box all x final 0 1.000000 y final 0 1.000000 z final 0 1.000000 remap units box", ) - self.job._structure_previous = atoms.copy() - self.job._structure_current = atoms.copy() - self.job._structure_previous.cell[1, 0] += 0.01 - self.job.interactive_cells_setter(self.job._structure_current.cell) + self.job._interactive_library._structure.cell[1, 0] += 0.01 + self.job.interactive_cells_setter(atoms.cell) self.assertEqual( - self.job._interactive_library._command[-1], + self.job._interactive_library._interactive_library._command[-1], "change_box all ortho", ) - self.job._structure_previous = atoms.copy() - self.job._structure_current = atoms.copy() - self.job._structure_current.cell[1, 0] += 0.01 - self.job.interactive_cells_setter(self.job._structure_current.cell) + structure_current = atoms.copy() + self.job._interactive_library._structure = atoms.copy() + structure_current.cell[1, 0] += 0.01 + self.job.interactive_cells_setter(structure_current.cell) self.assertEqual( - self.job._interactive_library._command[-2], + self.job._interactive_library._interactive_library._command[-2], "change_box all triclinic", ) def test_interactive_positions_setter(self): + self.job._interactive_library._prism = UnfoldingPrism(cell=self.job.structure.cell) self.job.interactive_positions_setter(np.arange(6).reshape(2, 3)) - self.assertTrue(self.job._interactive_library._command[0].startswith("x 1 3")) + self.assertTrue(self.job._interactive_library._interactive_library._command[0].startswith("x 1 3")) self.assertEqual( - self.job._interactive_library._command[1], "change_box all remap" + self.job._interactive_library._interactive_library._command[1], "change_box all remap" ) def test_interactive_execute(self): self.job._interactive_lammps_input() self.assertEqual( - self.job._interactive_library._command, + self.job._interactive_library._interactive_library._command, [ "fix ensemble all nve", "variable dumptime equal 100", @@ -120,25 +141,25 @@ def test_calc_minimize_input(self): self.minimize_job._interactive_lammps_input() self.assertEqual( - self.minimize_control_job._interactive_library._command, - self.minimize_job._interactive_library._command + self.minimize_control_job._interactive_library._interactive_library._command, + self.minimize_job._interactive_library._interactive_library._command ) # Ensure that pressure inputs are being parsed OK self.minimize_job.calc_minimize(pressure=0) self.minimize_job._interactive_lammps_input() self.assertTrue(("fix ensemble all box/relax iso 0.0" in - self.minimize_job._interactive_library._command)) + self.minimize_job._interactive_library._interactive_library._command)) self.minimize_job.calc_minimize(pressure=[0.0, 0.0, 0.0]) self.minimize_job._interactive_lammps_input() self.assertTrue(("fix ensemble all box/relax x 0.0 y 0.0 z 0.0 couple none" in - self.minimize_job._interactive_library._command)) + self.minimize_job._interactive_library._interactive_library._command)) self.minimize_job.calc_minimize(pressure=[1, 2, None, 0., 0., None]) self.minimize_job._interactive_lammps_input() self.assertTrue(("fix ensemble all box/relax x 10000.0 y 20000.0 xy 0.0 xz 0.0 couple none" in - self.minimize_job._interactive_library._command)) + self.minimize_job._interactive_library._interactive_library._command)) def test_fix_external(self): def f(x, nt, nl): From 294534ccea8754f79168262fb8add02965cbd407 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 15:47:19 -0700 Subject: [PATCH 09/22] Add library parameter --- pyiron_atomistics/lammps/wrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 73e5ba283..42ececd24 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -16,7 +16,7 @@ class PyironLammpsLibrary(object): def __init__( - self, working_directory, cores=1, comm=None, logger=None, log_file=None + self, working_directory, cores=1, comm=None, logger=None, log_file=None, library=None ): self._logger = logger self._prism = None From 802504e4654272480847ad9b30db6638f62910be Mon Sep 17 00:00:00 2001 From: pyiron-runner Date: Mon, 27 Feb 2023 22:49:24 +0000 Subject: [PATCH 10/22] Format black --- pyiron_atomistics/lammps/wrapper.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 42ececd24..e2456b3c8 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -16,7 +16,13 @@ class PyironLammpsLibrary(object): def __init__( - self, working_directory, cores=1, comm=None, logger=None, log_file=None, library=None + self, + working_directory, + cores=1, + comm=None, + logger=None, + log_file=None, + library=None, ): self._logger = logger self._prism = None From a8befe69ece4c4862f399be65b4f35459cfc6440 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 27 Feb 2023 16:21:13 -0700 Subject: [PATCH 11/22] another bug fix --- pyiron_atomistics/lammps/wrapper.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index 42ececd24..31cb744c8 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -134,12 +134,13 @@ def interactive_structure_setter( el_eam_lst, calc_md=True, ): - old_symbols = self._structure.get_species_symbols() - new_symbols = structure.get_species_symbols() - if any(old_symbols != new_symbols): - raise ValueError( - f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" - ) + if self._structure is not None: + old_symbols = self._structure.get_species_symbols() + new_symbols = structure.get_species_symbols() + if any(old_symbols != new_symbols): + raise ValueError( + f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" + ) self.interactive_lib_command(command="clear") control_dict = self._set_selective_dynamics( structure=structure, calc_md=calc_md From 3d59ea39aa64e3f587e655e8489caf2a16125aae Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Fri, 3 Mar 2023 08:48:50 -0700 Subject: [PATCH 12/22] Update wrapper.py --- pyiron_atomistics/lammps/wrapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py index e66bd2107..964966e2f 100644 --- a/pyiron_atomistics/lammps/wrapper.py +++ b/pyiron_atomistics/lammps/wrapper.py @@ -6,7 +6,7 @@ import warnings from pyiron_atomistics.lammps.structure import UnfoldingPrism -from pyiron_atomistics.lammps.base import _check_ortho_prism +from pyiron_atomistics.lammps.output import _check_ortho_prism try: # mpi4py is only supported on Linux and Mac Os X from pylammpsmpi import LammpsLibrary From 651d8061b6b01f15be8282813a7e2a2616f391e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20Jan=C3=9Fen?= Date: Sun, 14 May 2023 11:05:28 -0700 Subject: [PATCH 13/22] Use pyiron_lammps --- .ci_support/environment.yml | 1 + pyiron_atomistics/lammps/interactive.py | 2 +- pyiron_atomistics/lammps/wrapper.py | 475 ------------------------ setup.py | 1 + tests/lammps/test_interactive.py | 2 +- 5 files changed, 4 insertions(+), 477 deletions(-) delete mode 100644 pyiron_atomistics/lammps/wrapper.py diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index ad534ff85..565e902e0 100644 --- a/.ci_support/environment.yml +++ b/.ci_support/environment.yml @@ -16,6 +16,7 @@ dependencies: - phonopy =2.17.1 - pint =0.20.1 - pyiron_base =0.5.33 +- pyiron_lammps =0.2.0 - pymatgen =2022.11.7 - pyscal =2.10.18 - scikit-learn =1.2.1 diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 1f236d929..3edfe7d2c 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -11,7 +11,7 @@ from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.control import LammpsControl from pyiron_atomistics.atomistics.job.interactive import GenericInteractive -from pyiron_atomistics.lammps.wrapper import PyironLammpsLibrary +from pyiron_lammps.wrapper import PyironLammpsLibrary from pyiron_atomistics.lammps.units import UnitConverter __author__ = "Osamu Waseda, Jan Janssen" diff --git a/pyiron_atomistics/lammps/wrapper.py b/pyiron_atomistics/lammps/wrapper.py deleted file mode 100644 index 964966e2f..000000000 --- a/pyiron_atomistics/lammps/wrapper.py +++ /dev/null @@ -1,475 +0,0 @@ -from ctypes import c_double, c_int -import importlib -import numpy as np -import os -from scipy import constants -import warnings - -from pyiron_atomistics.lammps.structure import UnfoldingPrism -from pyiron_atomistics.lammps.output import _check_ortho_prism - -try: # mpi4py is only supported on Linux and Mac Os X - from pylammpsmpi import LammpsLibrary -except ImportError: - pass - - -class PyironLammpsLibrary(object): - def __init__( - self, - working_directory, - cores=1, - comm=None, - logger=None, - log_file=None, - library=None, - ): - self._logger = logger - self._prism = None - self._structure = None - self._cores = cores - if library is not None: - self._interactive_library = library - elif self._cores == 1: - lammps = getattr(importlib.import_module("lammps"), "lammps") - if log_file is None: - log_file = os.path.join(working_directory, "log.lammps") - self._interactive_library = lammps( - cmdargs=["-screen", "none", "-log", log_file], - comm=comm, - ) - else: - self._interactive_library = LammpsLibrary( - cores=self._cores, working_directory=working_directory - ) - - def interactive_lib_command(self, command): - if self._logger is not None: - self._logger.debug("Lammps library: " + command) - self._interactive_library.command(command) - - def interactive_positions_getter(self): - positions = np.reshape( - np.array(self._interactive_library.gather_atoms("x", 1, 3)), - (len(self._structure), 3), - ) - if _check_ortho_prism(prism=self._prism): - positions = np.matmul(positions, self._prism.R.T) - return positions - - def interactive_positions_setter(self, positions): - if _check_ortho_prism(prism=self._prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, self._prism.R) - positions = np.array(positions).flatten() - if self._cores == 1: - self._interactive_library.scatter_atoms( - "x", 1, 3, (len(positions) * c_double)(*positions) - ) - else: - self._interactive_library.scatter_atoms("x", positions) - self.interactive_lib_command(command="change_box all remap") - - def interactive_cells_getter(self): - cc = np.array( - [ - [self._interactive_library.get_thermo("lx"), 0, 0], - [ - self._interactive_library.get_thermo("xy"), - self._interactive_library.get_thermo("ly"), - 0, - ], - [ - self._interactive_library.get_thermo("xz"), - self._interactive_library.get_thermo("yz"), - self._interactive_library.get_thermo("lz"), - ], - ] - ) - return self._prism.unfold_cell(cc) - - def interactive_cells_setter(self, cell): - self._prism = UnfoldingPrism(cell) - lx, ly, lz, xy, xz, yz = self._prism.get_lammps_prism() - if _check_ortho_prism(prism=self._prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" - ) - - is_skewed = cell_is_skewed(cell=cell, tolerance=1.0e-8) - was_skewed = cell_is_skewed(cell=self._structure.cell, tolerance=1.0e-8) - - if is_skewed: - if not was_skewed: - self.interactive_lib_command(command="change_box all triclinic") - self.interactive_lib_command( - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, xy, xz, yz), - ) - elif was_skewed: - self.interactive_lib_command( - command="change_box all x final 0 %f y final 0 %f z final 0 %f xy final %f xz final %f yz final %f remap units box" - % (lx, ly, lz, 0.0, 0.0, 0.0), - ) - self.interactive_lib_command(command="change_box all ortho") - else: - self.interactive_lib_command( - command="change_box all x final 0 %f y final 0 %f z final 0 %f remap units box" - % (lx, ly, lz), - ) - - def interactive_volume_getter(self): - return self._interactive_library.get_thermo("vol") - - def interactive_forces_getter(self): - ff = np.reshape( - np.array(self._interactive_library.gather_atoms("f", 1, 3)), - (len(self._structure), 3), - ) - if _check_ortho_prism(prism=self._prism): - ff = np.matmul(ff, self._prism.R.T) - return ff - - def interactive_structure_setter( - self, - structure, - units, - dimension, - boundary, - atom_style, - el_eam_lst, - calc_md=True, - ): - if self._structure is not None: - old_symbols = self._structure.get_species_symbols() - new_symbols = structure.get_species_symbols() - if any(old_symbols != new_symbols): - raise ValueError( - f"structure has different chemical symbols than old one: {new_symbols} != {old_symbols}" - ) - self.interactive_lib_command(command="clear") - control_dict = self._set_selective_dynamics( - structure=structure, calc_md=calc_md - ) - self.interactive_lib_command(command="units " + units) - self.interactive_lib_command(command="dimension " + str(dimension)) - self.interactive_lib_command(command="boundary " + boundary) - self.interactive_lib_command(command="atom_style " + atom_style) - - self.interactive_lib_command(command="atom_modify map array") - self._prism = UnfoldingPrism(structure.cell) - if _check_ortho_prism(prism=self._prism): - warnings.warn( - "Warning: setting upper trangular matrix might slow down the calculation" - ) - xhi, yhi, zhi, xy, xz, yz = self._prism.get_lammps_prism() - if self._prism.is_skewed(): - self.interactive_lib_command( - command="region 1 prism" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " " - + str(xy) - + " " - + str(xz) - + " " - + str(yz) - + " units box", - ) - else: - self.interactive_lib_command( - command="region 1 block" - + " 0.0 " - + str(xhi) - + " 0.0 " - + str(yhi) - + " 0.0 " - + str(zhi) - + " units box", - ) - el_struct_lst = structure.get_species_symbols() - el_obj_lst = structure.get_species_objects() - if atom_style == "full": - self.interactive_lib_command( - command="create_box " - + str(len(el_eam_lst)) - + " 1 " - + "bond/types 1 " - + "angle/types 1 " - + "extra/bond/per/atom 2 " - + "extra/angle/per/atom 2 ", - ) - else: - self.interactive_lib_command( - command="create_box " + str(len(el_eam_lst)) + " 1" - ) - el_dict = {} - for id_eam, el_eam in enumerate(el_eam_lst): - if el_eam in el_struct_lst: - id_el = list(el_struct_lst).index(el_eam) - el = el_obj_lst[id_el] - el_dict[el] = id_eam + 1 - self.interactive_lib_command( - command="mass {0:3d} {1:f}".format(id_eam + 1, el.AtomicMass), - ) - else: - self.interactive_lib_command( - command="mass {0:3d} {1:f}".format(id_eam + 1, 1.00), - ) - positions = structure.positions.flatten() - if _check_ortho_prism(prism=self._prism): - positions = np.array(positions).reshape(-1, 3) - positions = np.matmul(positions, self._prism.R) - positions = positions.flatten() - try: - elem_all = np.array( - [el_dict[el] for el in structure.get_chemical_elements()] - ) - except KeyError: - missing = set(structure.get_chemical_elements()).difference(el_dict.keys()) - missing = ", ".join([el.Abbreviation for el in missing]) - raise ValueError( - f"Structure contains elements [{missing}], that are not present in the potential!" - ) - if self._cores == 1: - self._interactive_library.create_atoms( - n=len(structure), - id=None, - type=(len(elem_all) * c_int)(*elem_all), - x=(len(positions) * c_double)(*positions), - v=None, - image=None, - shrinkexceed=False, - ) - else: - self._interactive_library.create_atoms( - n=len(structure), - id=None, - type=elem_all, - x=positions, - v=None, - image=None, - shrinkexceed=False, - ) - self.interactive_lib_command(command="change_box all remap") - for key, value in control_dict.items(): - self.interactive_lib_command(command=key + " " + value) - self._structure = structure - - @staticmethod - def _set_selective_dynamics(structure, calc_md): - control_dict = {} - if "selective_dynamics" in structure._tag_list.keys(): - if structure.selective_dynamics._default is None: - structure.selective_dynamics._default = [True, True, True] - sel_dyn = np.logical_not(structure.selective_dynamics.list()) - # Enter loop only if constraints present - if len(np.argwhere(np.any(sel_dyn, axis=1)).flatten()) != 0: - all_indices = np.arange(len(structure), dtype=int) - constraint_xyz = np.argwhere(np.all(sel_dyn, axis=1)).flatten() - not_constrained_xyz = np.setdiff1d(all_indices, constraint_xyz) - # LAMMPS starts counting from 1 - constraint_xyz += 1 - ind_x = np.argwhere(sel_dyn[not_constrained_xyz, 0]).flatten() - ind_y = np.argwhere(sel_dyn[not_constrained_xyz, 1]).flatten() - ind_z = np.argwhere(sel_dyn[not_constrained_xyz, 2]).flatten() - constraint_xy = not_constrained_xyz[np.intersect1d(ind_x, ind_y)] + 1 - constraint_yz = not_constrained_xyz[np.intersect1d(ind_y, ind_z)] + 1 - constraint_zx = not_constrained_xyz[np.intersect1d(ind_z, ind_x)] + 1 - constraint_x = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_x, ind_y), ind_z)] - + 1 - ) - constraint_y = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_y, ind_z), ind_x)] - + 1 - ) - constraint_z = ( - not_constrained_xyz[np.setdiff1d(np.setdiff1d(ind_z, ind_x), ind_y)] - + 1 - ) - control_dict = {} - if len(constraint_xyz) > 0: - control_dict["group constraintxyz"] = "id " + " ".join( - [str(ind) for ind in constraint_xyz] - ) - control_dict[ - "fix constraintxyz" - ] = "constraintxyz setforce 0.0 0.0 0.0" - if calc_md: - control_dict["velocity constraintxyz"] = "set 0.0 0.0 0.0" - if len(constraint_xy) > 0: - control_dict["group constraintxy"] = "id " + " ".join( - [str(ind) for ind in constraint_xy] - ) - control_dict[ - "fix constraintxy" - ] = "constraintxy setforce 0.0 0.0 NULL" - if calc_md: - control_dict["velocity constraintxy"] = "set 0.0 0.0 NULL" - if len(constraint_yz) > 0: - control_dict["group constraintyz"] = "id " + " ".join( - [str(ind) for ind in constraint_yz] - ) - control_dict[ - "fix constraintyz" - ] = "constraintyz setforce NULL 0.0 0.0" - if calc_md: - control_dict["velocity constraintyz"] = "set NULL 0.0 0.0" - if len(constraint_zx) > 0: - control_dict["group constraintxz"] = "id " + " ".join( - [str(ind) for ind in constraint_zx] - ) - control_dict[ - "fix constraintxz" - ] = "constraintxz setforce 0.0 NULL 0.0" - if calc_md: - control_dict["velocity constraintxz"] = "set 0.0 NULL 0.0" - if len(constraint_x) > 0: - control_dict["group constraintx"] = "id " + " ".join( - [str(ind) for ind in constraint_x] - ) - control_dict[ - "fix constraintx" - ] = "constraintx setforce 0.0 NULL NULL" - if calc_md: - control_dict["velocity constraintx"] = "set 0.0 NULL NULL" - if len(constraint_y) > 0: - control_dict["group constrainty"] = "id " + " ".join( - [str(ind) for ind in constraint_y] - ) - control_dict[ - "fix constrainty" - ] = "constrainty setforce NULL 0.0 NULL" - if calc_md: - control_dict["velocity constrainty"] = "set NULL 0.0 NULL" - if len(constraint_z) > 0: - control_dict["group constraintz"] = "id " + " ".join( - [str(ind) for ind in constraint_z] - ) - control_dict[ - "fix constraintz" - ] = "constraintz setforce NULL NULL 0.0" - if calc_md: - control_dict["velocity constraintz"] = "set NULL NULL 0.0" - return control_dict - - def interactive_indices_getter(self): - return np.array(self._interactive_library.gather_atoms("type", 0, 1)) - - def interactive_energy_pot_getter(self): - return self._interactive_library.get_thermo("pe") - - def interactive_energy_tot_getter(self): - return self._interactive_library.get_thermo("etotal") - - def interactive_steps_getter(self): - return self._interactive_library.get_thermo("step") - - def interactive_temperatures_getter(self): - return self._interactive_library.get_thermo("temp") - - def interactive_pressures_getter(self): - pp = np.array( - [ - [ - self._interactive_library.get_thermo("pxx"), - self._interactive_library.get_thermo("pxy"), - self._interactive_library.get_thermo("pxz"), - ], - [ - self._interactive_library.get_thermo("pxy"), - self._interactive_library.get_thermo("pyy"), - self._interactive_library.get_thermo("pyz"), - ], - [ - self._interactive_library.get_thermo("pxz"), - self._interactive_library.get_thermo("pyz"), - self._interactive_library.get_thermo("pzz"), - ], - ] - ) - if _check_ortho_prism(prism=self._prism): - rotation_matrix = self._prism.R.T - pp = rotation_matrix.T @ pp @ rotation_matrix - return pp - - def interactive_indices_setter(self, indices, el_eam_lst): - el_struct_lst = self._structure.get_species_symbols() - el_obj_lst = self._structure.get_species_objects() - el_dict = {} - for id_eam, el_eam in enumerate(el_eam_lst): - if el_eam in el_struct_lst: - id_el = list(el_struct_lst).index(el_eam) - el = el_obj_lst[id_el] - el_dict[el] = id_eam + 1 - elem_all = np.array([el_dict[self._structure.species[el]] for el in indices]) - if self._cores == 1: - self._interactive_library.scatter_atoms( - "type", 0, 1, (len(elem_all) * c_int)(*elem_all) - ) - else: - self._interactive_library.scatter_atoms("type", elem_all) - - def interactive_stress_getter(self, enable_stress_computation=True): - """ - This gives back an Nx3x3 array of stress/atom defined in http://lammps.sandia.gov/doc/compute_stress_atom.html - Keep in mind that it is stress*volume in eV. Further discussion can be found on the website above. - - Returns: - numpy.array: Nx3x3 np array of stress/atom - """ - if enable_stress_computation: - self.interactive_lib_command("compute st all stress/atom NULL") - self.interactive_lib_command("run 0") - id_lst = self._interactive_library.extract_atom("id", 0) - id_lst = np.array([id_lst[i] for i in range(len(self._structure))]) - 1 - id_lst = np.arange(len(id_lst))[np.argsort(id_lst)] - ind = np.array([0, 3, 4, 3, 1, 5, 4, 5, 2]) - ss = self._interactive_library.extract_compute("st", 1, 2) - ss = np.array( - [ss[i][j] for i in range(len(self._structure)) for j in range(6)] - ).reshape(-1, 6)[id_lst] - ss = ( - ss[:, ind].reshape(len(self._structure), 3, 3) - / constants.eV - * constants.bar - * constants.angstrom**3 - ) - if _check_ortho_prism(prism=self._prism): - ss = np.einsum("ij,njk->nik", self._prism.R, ss) - ss = np.einsum("nij,kj->nik", ss, self._prism.R) - return ss - - def close(self): - if self._interactive_library is not None: - self._interactive_library.close() - - def set_fix_external_callback(self, fix_id, callback, caller=None): - self._interactive_library.set_fix_external_callback( - fix_id=fix_id, callback=callback, caller=caller - ) - - -def cell_is_skewed(cell, 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 = np.abs(np.linalg.det(cell)) - prod = np.linalg.norm(cell, axis=-1).prod() - if volume > 0: - if abs(volume - prod) / volume < tolerance: - return False - return True diff --git a/setup.py b/setup.py index ba0a0855d..62f3e49c7 100644 --- a/setup.py +++ b/setup.py @@ -54,6 +54,7 @@ 'phonopy==2.17.1', 'pint==0.20.1', 'pyiron_base==0.5.33', + 'pyiron_lammps==0.2.0' 'pymatgen==2022.11.7', 'scipy==1.10.1', 'seekpath==2.0.1', diff --git a/tests/lammps/test_interactive.py b/tests/lammps/test_interactive.py index f76626054..345faa7bb 100644 --- a/tests/lammps/test_interactive.py +++ b/tests/lammps/test_interactive.py @@ -8,7 +8,7 @@ from pyiron_base import Project, ProjectHDFio from pyiron_atomistics.atomistics.structure.atoms import Atoms from pyiron_atomistics.lammps.lammps import Lammps -from pyiron_atomistics.lammps.wrapper import PyironLammpsLibrary, cell_is_skewed +from pyiron_lammps.wrapper import PyironLammpsLibrary from pyiron_atomistics.lammps.structure import UnfoldingPrism From b86f957809cb8e41121cddaf76cd62c6cffcc16b Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Sun, 14 May 2023 17:23:56 -0600 Subject: [PATCH 14/22] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index f5180dca5..bea1c717b 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,7 @@ 'phonopy==2.18.0', 'pint==0.21', 'pyiron_base==0.5.37', - 'pyiron_lammps==0.2.0' + 'pyiron_lammps==0.2.0', 'pymatgen==2023.3.23', 'scipy==1.10.1', 'seekpath==2.0.1', From 4fa2238617347c77cd62201a1b6e93763225f432 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20Jan=C3=9Fen?= Date: Wed, 16 Aug 2023 22:16:06 -0600 Subject: [PATCH 15/22] Switch from pyiron_lammps to pylammpsmpi --- .ci_support/environment.yml | 2 +- pyiron_atomistics/lammps/interactive.py | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index e0b8c9f66..9a3060f48 100644 --- a/.ci_support/environment.yml +++ b/.ci_support/environment.yml @@ -15,7 +15,7 @@ dependencies: - phonopy =2.18.0 - pint =0.21 - pyiron_base =0.5.37 -- pyiron_lammps =0.2.0 +- pylammpsmpi =0.2.2 - pymatgen =2023.3.23 - pyscal =2.10.18 - scikit-learn =1.2.2 diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 3edfe7d2c..48f5e401a 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -7,13 +7,21 @@ import pandas as pd import warnings +from pyiron_base import ImportAlarm + from pyiron_atomistics.lammps.base import LammpsBase from pyiron_atomistics.lammps.structure import UnfoldingPrism from pyiron_atomistics.lammps.control import LammpsControl from pyiron_atomistics.atomistics.job.interactive import GenericInteractive -from pyiron_lammps.wrapper import PyironLammpsLibrary from pyiron_atomistics.lammps.units import UnitConverter +with ImportAlarm( + "Gpaw relies on the gpaw module but this is unavailable. Please ensure your python environment contains gpaw, " + "e.g. by running `conda install -c conda-forge gpaw`." +) as import_alarm: + from pylammpsmpi import LammpsASELibrary + + __author__ = "Osamu Waseda, Jan Janssen" __copyright__ = ( "Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH - " @@ -168,9 +176,10 @@ def _reset_interactive_run_command(self): df = pd.DataFrame(self.input.control.dataset) self._interactive_run_command = " ".join(df.T[df.index[-1]].values) + @import_alarm def interactive_initialize_interface(self): self._create_working_directory() - self._interactive_library = PyironLammpsLibrary( + self._interactive_library = LammpsASELibrary( working_directory=self.working_directory, cores=self.server.cores, comm=self._interactive_mpi_communicator, From fe9c9259c9230c5acc731e31fd7224881b01f90a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20Jan=C3=9Fen?= Date: Wed, 16 Aug 2023 22:23:05 -0600 Subject: [PATCH 16/22] Fixes --- .ci_support/environment.yml | 1 - pyiron_atomistics/lammps/interactive.py | 4 ++-- setup.py | 1 - 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index 28bdde114..dd6261bfe 100644 --- a/.ci_support/environment.yml +++ b/.ci_support/environment.yml @@ -15,7 +15,6 @@ dependencies: - phonopy =2.20.0 - pint =0.22 - pyiron_base =0.6.3 -- pylammpsmpi =0.2.2 - pymatgen =2023.8.10 - pyscal =2.10.18 - scikit-learn =1.3.0 diff --git a/pyiron_atomistics/lammps/interactive.py b/pyiron_atomistics/lammps/interactive.py index 48f5e401a..abeffc7b8 100644 --- a/pyiron_atomistics/lammps/interactive.py +++ b/pyiron_atomistics/lammps/interactive.py @@ -16,8 +16,8 @@ from pyiron_atomistics.lammps.units import UnitConverter with ImportAlarm( - "Gpaw relies on the gpaw module but this is unavailable. Please ensure your python environment contains gpaw, " - "e.g. by running `conda install -c conda-forge gpaw`." + "Lammps interactive relies on the lammps module but this is unavailable. Please ensure your python environment" + "contains lammps, e.g. by running `conda install -c conda-forge lammps`." ) as import_alarm: from pylammpsmpi import LammpsASELibrary diff --git a/setup.py b/setup.py index fe8ca7606..70836302c 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,6 @@ 'phonopy==2.20.0', 'pint==0.22', 'pyiron_base==0.6.3', - 'pylammpsmpi==0.2.2', 'pymatgen==2023.8.10', 'scipy==1.11.1', 'seekpath==2.1.0', From 04fd1353b1d2d8ac95a6a513fc6ebe2cc1de50ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20Jan=C3=9Fen?= Date: Wed, 16 Aug 2023 22:41:44 -0600 Subject: [PATCH 17/22] update tests but they still require the lammps library --- tests/lammps/test_interactive.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/lammps/test_interactive.py b/tests/lammps/test_interactive.py index 345faa7bb..a1d354c76 100644 --- a/tests/lammps/test_interactive.py +++ b/tests/lammps/test_interactive.py @@ -8,7 +8,7 @@ from pyiron_base import Project, ProjectHDFio from pyiron_atomistics.atomistics.structure.atoms import Atoms from pyiron_atomistics.lammps.lammps import Lammps -from pyiron_lammps.wrapper import PyironLammpsLibrary +from pylammpsmpi import LammpsASELibrary from pyiron_atomistics.lammps.structure import UnfoldingPrism @@ -25,7 +25,7 @@ def scatter_atoms(self, *args): class TestLammpsInteractive(unittest.TestCase): def setUp(self): - self.job._interactive_library = PyironLammpsLibrary( + self.job._interactive_library = LammpsASELibrary( working_directory=self.job.working_directory, cores=1, comm=None, @@ -33,7 +33,7 @@ def setUp(self): log_file=None, library=InteractiveLibrary() ) - self.minimize_job._interactive_library = PyironLammpsLibrary( + self.minimize_job._interactive_library = LammpsASELibrary( working_directory=self.job.working_directory, cores=1, comm=None, @@ -41,7 +41,7 @@ def setUp(self): log_file=None, library=InteractiveLibrary() ) - self.minimize_control_job._interactive_library = PyironLammpsLibrary( + self.minimize_control_job._interactive_library = LammpsASELibrary( working_directory=self.job.working_directory, cores=1, comm=None, From d2b06c3e9f9e42d4f8927b060aae10b19cd34f2f Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 20 Nov 2023 14:27:34 +0100 Subject: [PATCH 18/22] Update environment.yml --- .ci_support/environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index dd6261bfe..34b26611e 100644 --- a/.ci_support/environment.yml +++ b/.ci_support/environment.yml @@ -15,6 +15,7 @@ dependencies: - phonopy =2.20.0 - pint =0.22 - pyiron_base =0.6.3 +- pylammpsmpi =0.2.6 - pymatgen =2023.8.10 - pyscal =2.10.18 - scikit-learn =1.3.0 From d8c6146e8b61b60badbdaf3372257b1aac3870e8 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Mon, 20 Nov 2023 14:27:58 +0100 Subject: [PATCH 19/22] Update setup.py --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index 70836302c..5a7dcd12a 100644 --- a/setup.py +++ b/setup.py @@ -53,6 +53,7 @@ 'phonopy==2.20.0', 'pint==0.22', 'pyiron_base==0.6.3', + 'pylammpsmpi==0.2.6', 'pymatgen==2023.8.10', 'scipy==1.11.1', 'seekpath==2.1.0', From d427080759d5d7f248a67276e77bda2053054aba Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Thu, 23 Nov 2023 08:27:58 +0100 Subject: [PATCH 20/22] Update environment.yml --- .ci_support/environment.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml index a2e6f2d73..c82937199 100644 --- a/.ci_support/environment.yml +++ b/.ci_support/environment.yml @@ -18,7 +18,6 @@ dependencies: - pint =0.22 - pyiron_base =0.6.9 - pylammpsmpi =0.2.6 -- pymatgen =2023.8.10 - pyscal =2.10.18 - scikit-learn =1.3.2 - scipy =1.11.3 From 0c7583730b370046e15c6bffb23c014272fbe046 Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Thu, 23 Nov 2023 08:36:16 +0100 Subject: [PATCH 21/22] Update setup.py --- setup.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/setup.py b/setup.py index 901666b3d..ed6db552b 100644 --- a/setup.py +++ b/setup.py @@ -54,8 +54,6 @@ 'phonopy==2.20.0', 'pint==0.22', 'pyiron_base==0.6.9', - 'pylammpsmpi==0.2.6', - 'pymatgen==2023.8.10', 'scipy==1.11.3', 'seekpath==2.1.0', 'scikit-learn==1.3.2', From 9b6d133c99a6b71e1235991019b9529d1eac8f2a Mon Sep 17 00:00:00 2001 From: Jan Janssen Date: Thu, 23 Nov 2023 08:37:39 +0100 Subject: [PATCH 22/22] Update setup.py --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index ed6db552b..c874a0aa2 100644 --- a/setup.py +++ b/setup.py @@ -54,6 +54,7 @@ 'phonopy==2.20.0', 'pint==0.22', 'pyiron_base==0.6.9', + 'pylammpsmpi==0.2.6', 'scipy==1.11.3', 'seekpath==2.1.0', 'scikit-learn==1.3.2',