From 938867d6bc2cfb2b8afc85c73e1d0a8f3facbcf2 Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Fri, 20 Jan 2023 11:55:32 +0100 Subject: [PATCH 01/10] fix random and collect cells bugs --- pyiron_contrib/atomistics/ipi/ipi_core.py | 4 +++- pyiron_contrib/atomistics/ipi/ipi_jobs.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index 0897ba392..fa26d11ee 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -108,6 +108,7 @@ def run_static(self): [self.working_directory + '/./run_ipi.sh', self.working_directory, str(self.server.cores)]) self.collect_output() self.to_hdf() + self.status.finished = True self.compress() def collect_props(self): @@ -137,7 +138,8 @@ def collect_props(self): self.custom_output.energy_tot = self.custom_output.energy_pot + self.custom_output.energy_kin def collect_cells(self): - f = open(self.working_directory + '/ipi_out.pos_0.xyz') + digits = "{0:0" + str(len(str(self.custom_input.n_beads))) + "}" + f = open(self.working_directory + '/ipi_out.pos_' + digits.format(0) + '.xyz') lines = f.readlines() abc = [] ABC = [] diff --git a/pyiron_contrib/atomistics/ipi/ipi_jobs.py b/pyiron_contrib/atomistics/ipi/ipi_jobs.py index 684c45a87..0fabb2b44 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_jobs.py +++ b/pyiron_contrib/atomistics/ipi/ipi_jobs.py @@ -6,7 +6,7 @@ from pyiron_contrib.atomistics.ipi.ipi_core import IPiCore from xml.etree import ElementTree -from numpy import random +import random from shutil import copy __author__ = "Raynol Dsouza" From 8a718867a8deac9c5088b81d6abd1af8e0092c57 Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Fri, 28 Apr 2023 13:10:55 +0200 Subject: [PATCH 02/10] collect only centroid positions and forces --- pyiron_contrib/atomistics/ipi/ipi_core.py | 45 +++++++++++++++++------ pyiron_contrib/atomistics/ipi/ipi_jobs.py | 11 +++--- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index fa26d11ee..b11ee23d4 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -137,24 +137,47 @@ def collect_props(self): self.custom_output.pressure = np.array([float(i) for i in pressure]) self.custom_output.energy_tot = self.custom_output.energy_pot + self.custom_output.energy_kin - def collect_cells(self): - digits = "{0:0" + str(len(str(self.custom_input.n_beads))) + "}" - f = open(self.working_directory + '/ipi_out.pos_' + digits.format(0) + '.xyz') + @staticmethod + def _collect_traj_helper(filename): + f = open(filename) lines = f.readlines() + starts = [] + for i, x in enumerate(lines): + if x.startswith("#"): + starts.append(i) + starts.append(len(lines) + 1) abc = [] ABC = [] - for x in lines: - if x.startswith("#"): - split_line = x.split() - abc.append([float(i) for i in split_line[2:5]]) - ABC.append([float(i) for i in split_line[5:8]]) + traj = [] + start = starts[0] + for stop in starts[1:]: + temp_list = [] + snap_lines = lines[start:stop - 1] + for i, l in enumerate(snap_lines): + split_l = l.split() + if i != 0: + temp_list.append([float(j) for j in split_l[1:]]) + else: + abc.append([float(j) for j in split_l[2:5]]) + ABC.append([float(j) for j in split_l[5:8]]) + start = stop + traj.append(temp_list) f.close() - self.custom_output.cell_abc = np.array(abc) - self.custom_output.cell_ABC = np.array(ABC) + return np.array(abc), np.array(ABC), np.array(traj) + + def collect_trajectory(self): + #digits = "{0:0" + str(len(str(self.custom_input.n_beads))) + "}" + #f = open(self.working_directory + '/ipi_out.pos_' + digits.format(0) + '.xyz') + positions_out = self._collect_traj_helper(self.working_directory + '/ipi_out.pos.xyz') + forces_out = self._collect_traj_helper(self.working_directory + '/ipi_out.for.xyz') + self.custom_output.cell_abc = np.array(positions_out[0]) + self.custom_output.cell_ABC = np.array(positions_out[1]) + self.custom_output.positions = np.array(positions_out[2]) + self.custom_output.forces = np.array(forces_out[2]) def collect_output(self): self.collect_props() - self.collect_cells() + self.collect_trajectory() def collect_rdf(self): f=open(self.working_directory + '/ipi_out.AlAl.rdf.dat', "r") diff --git a/pyiron_contrib/atomistics/ipi/ipi_jobs.py b/pyiron_contrib/atomistics/ipi/ipi_jobs.py index 0fabb2b44..137e81343 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_jobs.py +++ b/pyiron_contrib/atomistics/ipi/ipi_jobs.py @@ -25,7 +25,7 @@ def __init__(self, project, job_name): def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep=1., temperature_damping_timescale=100., pressure_damping_timescale=1000., - n_ionic_steps=100, prop_n_print=1, traj_n_print=1, seed=None, port=31415, *args, **kwargs): + n_ionic_steps=100, n_print=1, seed=None, port=31415, *args, **kwargs): if not seed: random.seed(self.job_name) seed = random.randint(1, 99999) @@ -36,8 +36,7 @@ def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep= self.custom_input.temperature_damping_timescale = temperature_damping_timescale self.custom_input.pressure_damping_timescale = pressure_damping_timescale self.custom_input.n_ionic_steps = n_ionic_steps - self.custom_input.prop_n_print = prop_n_print - self.custom_input.traj_n_print = traj_n_print + self.custom_input.n_print = n_print self.custom_input.seed = seed self.custom_input.port = port @@ -50,9 +49,9 @@ def ipi_write_helper(self, xml_filename): tree = ElementTree.parse(self.working_directory + '/' + xml_filename) root = tree.getroot() filepath = self.working_directory + '/ipi_input.xml' - root[0][0].attrib['stride'] = str(self.custom_input.prop_n_print) - for i in range(1, 4): - root[0][i].attrib['stride'] = str(self.custom_input.traj_n_print) + for i in range(0, 3): + root[0][i].attrib['stride'] = str(self.custom_input.n_print) + root[0][3].attrib['stride'] = str(self.custom_input.n_ionic_steps) root[1].text = str(self.custom_input.n_ionic_steps) root[2][0].text = str(self.custom_input.seed) root[3][0].text = self.job_name From 25e53297956d1c75d0e43e19ebb746fd0ec64d45 Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Fri, 28 Apr 2023 14:23:41 +0200 Subject: [PATCH 03/10] default rdf thermalize should be 0 --- pyiron_contrib/atomistics/ipi/ipi_core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index b11ee23d4..82ec0728f 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -190,7 +190,7 @@ def collect_rdf(self): f.close() return np.array([float(i) for i in rdf_r]), np.array([float(i) for i in rdf_g_r]) - def get_rdf(self, r_min=2., r_max=5., bins=100, thermalize=50): + def get_rdf(self, r_min=2., r_max=5., bins=100, thermalize=0): self.decompress() rdf_list = [self.working_directory + '/./run_rdf.sh', self.working_directory, From eb334243a975e2507b178fb05f2def24e753f155 Mon Sep 17 00:00:00 2001 From: Raynol Dsouza Date: Mon, 8 May 2023 10:25:34 +0200 Subject: [PATCH 04/10] Update pyiron_contrib/atomistics/ipi/ipi_core.py Co-authored-by: Liam Huber --- pyiron_contrib/atomistics/ipi/ipi_core.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index 82ec0728f..d76d6a2fb 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -141,11 +141,9 @@ def collect_props(self): def _collect_traj_helper(filename): f = open(filename) lines = f.readlines() - starts = [] - for i, x in enumerate(lines): - if x.startswith("#"): - starts.append(i) - starts.append(len(lines) + 1) + starts = [ + i for i, x in enumerate(lines) if x.startswith("#") + ] + [len(lines) + 1] abc = [] ABC = [] traj = [] From 511bfe3105a9ac8158c8f3e85cc9c94516c6f0e7 Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Mon, 8 May 2023 11:09:02 +0200 Subject: [PATCH 05/10] make readability changes --- pyiron_contrib/atomistics/ipi/ipi_core.py | 37 ++++++++++++----------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index d76d6a2fb..eea444713 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -139,33 +139,36 @@ def collect_props(self): @staticmethod def _collect_traj_helper(filename): + """ + For each snapshot collect the cell box edges and angles and the corresponding atom positions. + """ f = open(filename) lines = f.readlines() + # the first line for each snapshot is always the number of atoms. + # Identify the line number where each snapshot starts. starts = [ i for i, x in enumerate(lines) if x.startswith("#") - ] + [len(lines) + 1] - abc = [] - ABC = [] - traj = [] + ] + [len(lines) + 1] # and also record at what line number the last snapshot ends + cell_abc = [] # cell box edges + cell_ABC = [] # cell box angles + trajectory = [] start = starts[0] - for stop in starts[1:]: - temp_list = [] + for stop in starts[1:]: # ignore the numer of atoms line snap_lines = lines[start:stop - 1] - for i, l in enumerate(snap_lines): - split_l = l.split() - if i != 0: - temp_list.append([float(j) for j in split_l[1:]]) - else: - abc.append([float(j) for j in split_l[2:5]]) - ABC.append([float(j) for j in split_l[5:8]]) + # second line contains the cell data. Collect it. + zeroth_line_values = [n for n in snap_lines[0].split()] + cell_abc.append(np.array(zeroth_line_values[2:5], dtype=float)) # collect cell box edges + cell_ABC.append(np.array(zeroth_line_values[5:8], dtype=float)) # collect cell box angles + # the remaining lines of the snapshot contain the positions/forces. Collect them. + temp_list = [] + for l in snap_lines[1:]: + temp_list.append([float(n) for n in l.split()[1:]]) # first column is the species. Ignore it. + trajectory.append(temp_list) start = stop - traj.append(temp_list) f.close() - return np.array(abc), np.array(ABC), np.array(traj) + return np.array(cell_abc), np.array(cell_ABC), np.array(trajectory) def collect_trajectory(self): - #digits = "{0:0" + str(len(str(self.custom_input.n_beads))) + "}" - #f = open(self.working_directory + '/ipi_out.pos_' + digits.format(0) + '.xyz') positions_out = self._collect_traj_helper(self.working_directory + '/ipi_out.pos.xyz') forces_out = self._collect_traj_helper(self.working_directory + '/ipi_out.for.xyz') self.custom_output.cell_abc = np.array(positions_out[0]) From 5022cd951cde5f1fc908d6cecd9b6561a4fd831a Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Mon, 8 May 2023 17:59:15 +0200 Subject: [PATCH 06/10] remove start=stop --- pyiron_contrib/atomistics/ipi/ipi_core.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index eea444713..5634370a2 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -152,9 +152,8 @@ def _collect_traj_helper(filename): cell_abc = [] # cell box edges cell_ABC = [] # cell box angles trajectory = [] - start = starts[0] - for stop in starts[1:]: # ignore the numer of atoms line - snap_lines = lines[start:stop - 1] + for i, stop in enumerate(starts[1:]): # ignore the numer of atoms line + snap_lines = lines[starts[i]:stop - 1] # second line contains the cell data. Collect it. zeroth_line_values = [n for n in snap_lines[0].split()] cell_abc.append(np.array(zeroth_line_values[2:5], dtype=float)) # collect cell box edges @@ -164,7 +163,6 @@ def _collect_traj_helper(filename): for l in snap_lines[1:]: temp_list.append([float(n) for n in l.split()[1:]]) # first column is the species. Ignore it. trajectory.append(temp_list) - start = stop f.close() return np.array(cell_abc), np.array(cell_ABC), np.array(trajectory) From 24749b8f97d46505af358b784d01fdd7be9aff81 Mon Sep 17 00:00:00 2001 From: raynol-dsouza Date: Tue, 9 May 2023 09:26:07 +0200 Subject: [PATCH 07/10] increase the time lammps remains open --- pyiron_contrib/atomistics/ipi/ipi_core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index 5634370a2..84966e2da 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -77,7 +77,7 @@ def write_input_lmp(self): "mass \t 1 " + str(mass) + "\n\n" + \ "include potential.inp\n\n" + \ "fix \t 1 all ipi " + self.job_name + " " + str(self.custom_input.port) + " unix\n" + \ - "run \t 5000000" + "run \t 10000000" with open(filepath, 'w') as file: file.writelines(data) From 81abaf51678c97aac353bd633689679e1cdd2803 Mon Sep 17 00:00:00 2001 From: Raynol Dsouza Date: Mon, 22 May 2023 17:44:39 +0200 Subject: [PATCH 08/10] make provision to constrain a list of atoms --- pyiron_contrib/atomistics/ipi/ipi_jobs.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_jobs.py b/pyiron_contrib/atomistics/ipi/ipi_jobs.py index 137e81343..158cbfe5d 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_jobs.py +++ b/pyiron_contrib/atomistics/ipi/ipi_jobs.py @@ -25,7 +25,7 @@ def __init__(self, project, job_name): def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep=1., temperature_damping_timescale=100., pressure_damping_timescale=1000., - n_ionic_steps=100, n_print=1, seed=None, port=31415, *args, **kwargs): + n_ionic_steps=100, n_print=1, seed=None, port=31415, constrain_ids=None, *args, **kwargs): if not seed: random.seed(self.job_name) seed = random.randint(1, 99999) @@ -39,6 +39,7 @@ def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep= self.custom_input.n_print = n_print self.custom_input.seed = seed self.custom_input.port = port + self.custom_input.constrain_ids = constrain_ids def write_template_file(self): copy(self._templates_directory + '/pimd_template.xml', self.working_directory + '/pimd_template.xml') @@ -59,6 +60,11 @@ def ipi_write_helper(self, xml_filename): root[4][0][1].text = str(self.custom_input.temperature) root[4][2][0][0][0].text = str(self.custom_input.pressure_damping_timescale) root[4][2][0][2].text = str(self.custom_input.timestep) + constrain_ids = self.custom_input.constrain_ids.copy() + if constrain_ids is not None: + if not isinstance(constrain_ids, list): + constrain_ids = constrain_ids.tolist() + root[4][2][2].text = str(constrain_ids) root[4][3][0].text = str(self.custom_input.temperature) root[4][3][1].text = str(self.custom_input.pressure) tree.write(filepath) From 7b560330c9b9d9678b99db58fff3e7eb63094c33 Mon Sep 17 00:00:00 2001 From: Raynol Dsouza Date: Tue, 23 May 2023 14:12:12 +0200 Subject: [PATCH 09/10] update --- pyiron_contrib/atomistics/ipi/ipi_jobs.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_jobs.py b/pyiron_contrib/atomistics/ipi/ipi_jobs.py index 158cbfe5d..8ca40a836 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_jobs.py +++ b/pyiron_contrib/atomistics/ipi/ipi_jobs.py @@ -60,10 +60,9 @@ def ipi_write_helper(self, xml_filename): root[4][0][1].text = str(self.custom_input.temperature) root[4][2][0][0][0].text = str(self.custom_input.pressure_damping_timescale) root[4][2][0][2].text = str(self.custom_input.timestep) - constrain_ids = self.custom_input.constrain_ids.copy() - if constrain_ids is not None: - if not isinstance(constrain_ids, list): - constrain_ids = constrain_ids.tolist() + if self.custom_input.constrain_ids is not None: + if not isinstance(self.custom_input.constrain_ids, list): + constrain_ids = self.custom_input.constrain_ids.tolist() root[4][2][2].text = str(constrain_ids) root[4][3][0].text = str(self.custom_input.temperature) root[4][3][1].text = str(self.custom_input.pressure) From 4c448d5c613a962f689c6648de8783cbf97a37f3 Mon Sep 17 00:00:00 2001 From: Raynol Dsouza Date: Tue, 22 Aug 2023 15:49:12 +0200 Subject: [PATCH 10/10] update ipi jobs --- pyiron_contrib/atomistics/ipi/ipi_core.py | 441 +++++++++++----------- pyiron_contrib/atomistics/ipi/ipi_jobs.py | 10 +- 2 files changed, 230 insertions(+), 221 deletions(-) diff --git a/pyiron_contrib/atomistics/ipi/ipi_core.py b/pyiron_contrib/atomistics/ipi/ipi_core.py index 84966e2da..4e3e5d171 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_core.py +++ b/pyiron_contrib/atomistics/ipi/ipi_core.py @@ -1,218 +1,223 @@ -# coding: utf-8 -# 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 pyiron_atomistics.lammps.interactive import LammpsInteractive -from pyiron_atomistics.lammps.base import Input -from pyiron_base.storage.datacontainer import DataContainer - -import numpy as np -import subprocess -from shutil import copy -from os.path import isdir - -__author__ = "Raynol Dsouza" -__copyright__ = "Copyright 2022, Max-Planck-Institut für Eisenforschung GmbH " \ - "- Computational Materials Design (CM) Department" -__version__ = "0.0" -__maintainer__ = "Raynol Dsouza" -__email__ = "dsouza@mpie.de" -__status__ = "development" -__date__ = "Jan 18, 2023" - -class IPiCore(LammpsInteractive): - - def __init__(self, project, job_name): - super().__init__(project, job_name) - self.input = Input() - self.custom_input = DataContainer(table_name='pimd_input') - self.custom_output = DataContainer(table_name='pimd_output') - self._templates_directory = None - - @property - def templates_directory(self): - return self._templates_directory - - @templates_directory.setter - def templates_directory(self, templates_directory): - if not isinstance(templates_directory, str): - raise TypeError('templates_directory must be a str!') - self._templates_directory = templates_directory - - def calc_npt_md(self): - pass - - def write_potential(self): - self.input.potential.write_file(file_name="potential.inp", cwd=self.working_directory) - self.input.potential.copy_pot_files(self.working_directory) - - def write_init_xyz(self): - filepath = self.working_directory + '/init.xyz' - self.structure.write(filename=filepath, format='xyz') - cell = self.structure.cell.array.diagonal() - angle = self.structure.cell.angles() - with open(filepath, 'r') as file: - data = file.readlines() - data[1] = "# CELL(abcABC): " \ - + str(cell[0]) + " " + str(cell[1]) + " " + str(cell[2]) + " " \ - + str(angle[0]) + " " + str(angle[1]) + " " + str(angle[2]) + \ - " positions{angstrom} cell{angstrom}\n" - with open(filepath, 'w') as file: - file.writelines(data) - - def write_data_lmp(self): - filepath = self.working_directory + '/data.lmp' - self.structure.write(filename=filepath, format='lammps-data') - - def write_input_lmp(self): - filepath = self.working_directory + '/input.lmp' - mass = self.structure.get_masses()[0] - data = "# LAMMPS input file\n\n" + \ - "atom_style \t atomic\n" + \ - "units \t metal\n" + \ - "dimension \t 3\n" + \ - "boundary \t p p p\n" + \ - "\n" + \ - "read_data \t data.lmp\n" + \ - "mass \t 1 " + str(mass) + "\n\n" + \ - "include potential.inp\n\n" + \ - "fix \t 1 all ipi " + self.job_name + " " + str(self.custom_input.port) + " unix\n" + \ - "run \t 10000000" - with open(filepath, 'w') as file: - file.writelines(data) - - def write_ipi_xml(self): - pass - - def write_shell_scripts(self): - copy(self._templates_directory + '/run_ipi.sh', self.working_directory + '/run_ipi.sh') - copy(self._templates_directory + '/run_rdf.sh', self.working_directory + '/run_rdf.sh') - - def write_template_file(self): - pass - - def write_input(self): - if not isdir(self.working_directory): - self.project_hdf5.create_working_directory() - super(IPiCore, self).write_input() - self.write_potential() - self.write_init_xyz() - self.write_data_lmp() - self.write_input_lmp() - self.write_shell_scripts() - self.write_template_file() - self.write_ipi_xml() - - def run_static(self): - subprocess.check_call( - [self.working_directory + '/./run_ipi.sh', self.working_directory, str(self.server.cores)]) - self.collect_output() - self.to_hdf() - self.status.finished = True - self.compress() - - def collect_props(self): - f = open(self.working_directory + '/ipi_out.out', "r") - lines = f.readlines() - time = [] - temperature = [] - energy_kin = [] - energy_pot = [] - volume = [] - pressure = [] - for x in lines: - if not x.startswith('#'): - time.append(x.split()[1]) - temperature.append(x.split()[2]) - energy_kin.append(x.split()[3]) - energy_pot.append(x.split()[4]) - volume.append(x.split()[5]) - pressure.append(x.split()[6]) - f.close() - self.custom_output.time = np.array([float(i) for i in time]) - self.custom_output.temperature = np.array([float(i) for i in temperature]) - self.custom_output.energy_kin = np.array([float(i) for i in energy_kin]) - self.custom_output.energy_pot = np.array([float(i) for i in energy_pot]) - self.custom_output.volume = np.array([float(i) for i in volume]) - self.custom_output.pressure = np.array([float(i) for i in pressure]) - self.custom_output.energy_tot = self.custom_output.energy_pot + self.custom_output.energy_kin - - @staticmethod - def _collect_traj_helper(filename): - """ - For each snapshot collect the cell box edges and angles and the corresponding atom positions. - """ - f = open(filename) - lines = f.readlines() - # the first line for each snapshot is always the number of atoms. - # Identify the line number where each snapshot starts. - starts = [ - i for i, x in enumerate(lines) if x.startswith("#") - ] + [len(lines) + 1] # and also record at what line number the last snapshot ends - cell_abc = [] # cell box edges - cell_ABC = [] # cell box angles - trajectory = [] - for i, stop in enumerate(starts[1:]): # ignore the numer of atoms line - snap_lines = lines[starts[i]:stop - 1] - # second line contains the cell data. Collect it. - zeroth_line_values = [n for n in snap_lines[0].split()] - cell_abc.append(np.array(zeroth_line_values[2:5], dtype=float)) # collect cell box edges - cell_ABC.append(np.array(zeroth_line_values[5:8], dtype=float)) # collect cell box angles - # the remaining lines of the snapshot contain the positions/forces. Collect them. - temp_list = [] - for l in snap_lines[1:]: - temp_list.append([float(n) for n in l.split()[1:]]) # first column is the species. Ignore it. - trajectory.append(temp_list) - f.close() - return np.array(cell_abc), np.array(cell_ABC), np.array(trajectory) - - def collect_trajectory(self): - positions_out = self._collect_traj_helper(self.working_directory + '/ipi_out.pos.xyz') - forces_out = self._collect_traj_helper(self.working_directory + '/ipi_out.for.xyz') - self.custom_output.cell_abc = np.array(positions_out[0]) - self.custom_output.cell_ABC = np.array(positions_out[1]) - self.custom_output.positions = np.array(positions_out[2]) - self.custom_output.forces = np.array(forces_out[2]) - - def collect_output(self): - self.collect_props() - self.collect_trajectory() - - def collect_rdf(self): - f=open(self.working_directory + '/ipi_out.AlAl.rdf.dat', "r") - lines=f.readlines() - rdf_r = [] - rdf_g_r = [] - for x in lines: - rdf_r.append(x.split()[0]) - rdf_g_r.append(x.split()[1]) - f.close() - return np.array([float(i) for i in rdf_r]), np.array([float(i) for i in rdf_g_r]) - - def get_rdf(self, r_min=2., r_max=5., bins=100, thermalize=0): - self.decompress() - rdf_list = [self.working_directory + '/./run_rdf.sh', - self.working_directory, - str(self.custom_input.temperature), - self.structure.get_chemical_symbols()[0], self.structure.get_chemical_symbols()[0], - str(bins), - str(r_min), str(r_max), - str(thermalize)] - subprocess.check_call(rdf_list) - rdf_r, rdf_g_r = self.collect_rdf() - self.compress() - return rdf_r, rdf_g_r - - def to_hdf(self, hdf=None, group_name=None): - super(IPiCore, self).to_hdf(hdf=hdf, group_name=group_name) - self._structure_to_hdf() - self.custom_input.templates_directory = self._templates_directory - self.custom_input.to_hdf(self._hdf5) - self.custom_output.to_hdf(self._hdf5) - - def from_hdf(self, hdf=None, group_name=None): - super(IPiCore, self).from_hdf(hdf=hdf, group_name=group_name) - self._structure_from_hdf() - self.custom_input.from_hdf(self._hdf5) - self._templates_directory = self.custom_input.templates_directory - self.custom_output.from_hdf(self._hdf5) +# coding: utf-8 +# 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 pyiron_atomistics.lammps.interactive import LammpsInteractive +from pyiron_atomistics.lammps.base import Input +from pyiron_base.storage.datacontainer import DataContainer + +import numpy as np +import subprocess +from shutil import copy +from os.path import isdir + +__author__ = "Raynol Dsouza" +__copyright__ = "Copyright 2022, Max-Planck-Institut für Eisenforschung GmbH " \ + "- Computational Materials Design (CM) Department" +__version__ = "0.0" +__maintainer__ = "Raynol Dsouza" +__email__ = "dsouza@mpie.de" +__status__ = "development" +__date__ = "Jan 18, 2023" + +class IPiCore(LammpsInteractive): + + def __init__(self, project, job_name): + super().__init__(project, job_name) + self.input = Input() + self.custom_input = DataContainer(table_name='pimd_input') + self.custom_output = DataContainer(table_name='pimd_output') + self._templates_directory = None + + @property + def templates_directory(self): + return self._templates_directory + + @templates_directory.setter + def templates_directory(self, templates_directory): + if not isinstance(templates_directory, str): + raise TypeError('templates_directory must be a str!') + self._templates_directory = templates_directory + + def calc_npt_md(self): + pass + + def write_potential(self): + self.input.potential.write_file(file_name="potential.inp", cwd=self.working_directory) + self.input.potential.copy_pot_files(self.working_directory) + + def write_init_xyz(self): + filepath = self.working_directory + '/init.xyz' + self.structure.write(filename=filepath, format='xyz') + cell = self.structure.cell.array.diagonal() + angle = self.structure.cell.angles() + with open(filepath, 'r') as file: + data = file.readlines() + data[1] = "# CELL(abcABC): " \ + + str(cell[0]) + " " + str(cell[1]) + " " + str(cell[2]) + " " \ + + str(angle[0]) + " " + str(angle[1]) + " " + str(angle[2]) + \ + " positions{angstrom} cell{angstrom}\n" + with open(filepath, 'w') as file: + file.writelines(data) + + def write_data_lmp(self): + filepath = self.working_directory + '/data.lmp' + self.structure.write(filename=filepath, format='lammps-data') + + def write_input_lmp(self): + filepath = self.working_directory + '/input.lmp' + mass = self.structure.get_masses()[0] + data = "# LAMMPS input file\n\n" + \ + "atom_style \t atomic\n" + \ + "units \t metal\n" + \ + "dimension \t 3\n" + \ + "boundary \t p p p\n" + \ + "\n" + \ + "read_data \t data.lmp\n" + \ + "mass \t 1 " + str(mass) + "\n\n" + \ + "include potential.inp\n\n" + \ + "fix \t 1 all ipi " + self.job_name + " " + str(self.custom_input.port) + " unix\n" + \ + "run \t 10000000" + with open(filepath, 'w') as file: + file.writelines(data) + + def write_ipi_xml(self): + pass + + def write_shell_scripts(self): + copy(self._templates_directory + '/run_ipi.sh', self.working_directory + '/run_ipi.sh') + copy(self._templates_directory + '/run_rdf.sh', self.working_directory + '/run_rdf.sh') + + def write_template_file(self): + pass + + def write_input(self): + if not isdir(self.working_directory): + self.project_hdf5.create_working_directory() + super(IPiCore, self).write_input() + self.write_potential() + self.write_init_xyz() + self.write_data_lmp() + self.write_input_lmp() + self.write_shell_scripts() + self.write_template_file() + self.write_ipi_xml() + + def run_static(self): + subprocess.check_call( + [self.working_directory + '/./run_ipi.sh', self.working_directory, str(self.server.cores)]) + self.collect_output() + self.to_hdf() + self.status.finished = True + self.compress() + + def collect_props(self): + f = open(self.working_directory + '/ipi_out.out', "r") + lines = f.readlines() + time = [] + temperature = [] + energy_kin = [] + energy_pot = [] + volume = [] + pressure = [] + for x in lines: + if not x.startswith('#'): + time.append(x.split()[1]) + temperature.append(x.split()[2]) + energy_kin.append(x.split()[3]) + energy_pot.append(x.split()[4]) + volume.append(x.split()[5]) + pressure.append(x.split()[6]) + f.close() + self.custom_output.time = np.array([float(i) for i in time]) + self.custom_output.temperature = np.array([float(i) for i in temperature]) + self.custom_output.energy_kin = np.array([float(i) for i in energy_kin]) + self.custom_output.energy_pot = np.array([float(i) for i in energy_pot]) + self.custom_output.volume = np.array([float(i) for i in volume]) + self.custom_output.pressure = np.array([float(i) for i in pressure]) + self.custom_output.energy_tot = self.custom_output.energy_pot + self.custom_output.energy_kin + + @staticmethod + def _collect_traj_helper(filename): + """ + For each snapshot collect the cell box edges and angles and the corresponding atom positions. + """ + f = open(filename) + lines = f.readlines() + # the first line for each snapshot is always the number of atoms. + # Identify the line number where each snapshot starts. + starts = [ + i for i, x in enumerate(lines) if x.startswith("#") + ] + [len(lines) + 1] # and also record at what line number the last snapshot ends + cell_abc = [] # cell box edges + cell_ABC = [] # cell box angles + trajectory = [] + for i, stop in enumerate(starts[1:]): # ignore the numer of atoms line + snap_lines = lines[starts[i]:stop - 1] + # second line contains the cell data. Collect it. + zeroth_line_values = [n for n in snap_lines[0].split()] + cell_abc.append(np.array(zeroth_line_values[2:5], dtype=float)) # collect cell box edges + cell_ABC.append(np.array(zeroth_line_values[5:8], dtype=float)) # collect cell box angles + # the remaining lines of the snapshot contain the positions/forces. Collect them. + temp_list = [] + for l in snap_lines[1:]: + temp_list.append([float(n) for n in l.split()[1:]]) # first column is the species. Ignore it. + trajectory.append(temp_list) + f.close() + return np.array(cell_abc), np.array(cell_ABC), np.array(trajectory) + + def collect_trajectory(self): + positions_out = self._collect_traj_helper(self.working_directory + '/ipi_out.x_pos.xyz') + forces_out = self._collect_traj_helper(self.working_directory + '/ipi_out.x_for.xyz') + self.custom_output.cell_abc = np.array(positions_out[0]) + self.custom_output.cell_ABC = np.array(positions_out[1]) + self.custom_output.positions = np.array(positions_out[2]) + self.custom_output.forces = np.array(forces_out[2]) + + def collect_output(self): + self.collect_props() + self.collect_trajectory() + + def collect_rdf(self): + element = self.structure.get_chemical_symbols()[0] + try: + f=open(self.working_directory + '/ipi_out.' + element + element + '.rdf.dat', "r") + except FileNotFoundError: + raise FileNotFoundError('No compiled fortran module for fast calculations have been found. ' + 'Proceeding the calculations is not possible.') + lines=f.readlines() + rdf_r = [] + rdf_g_r = [] + for x in lines: + rdf_r.append(x.split()[0]) + rdf_g_r.append(x.split()[1]) + f.close() + return np.array([float(i) for i in rdf_r]), np.array([float(i) for i in rdf_g_r]) + + def get_rdf(self, r_min=2., r_max=5., bins=100, thermalize=0): + self.decompress() + rdf_list = [self.working_directory + '/./run_rdf.sh', + self.working_directory, + str(self.custom_input.temperature), + self.structure.get_chemical_symbols()[0], self.structure.get_chemical_symbols()[0], + str(bins), + str(r_min), str(r_max), + str(thermalize)] + subprocess.check_call(rdf_list) + rdf_r, rdf_g_r = self.collect_rdf() + self.compress() + return rdf_r, rdf_g_r + + def to_hdf(self, hdf=None, group_name=None): + super(IPiCore, self).to_hdf(hdf=hdf, group_name=group_name) + self._structure_to_hdf() + self.custom_input.templates_directory = self._templates_directory + self.custom_input.to_hdf(self._hdf5) + self.custom_output.to_hdf(self._hdf5) + + def from_hdf(self, hdf=None, group_name=None): + super(IPiCore, self).from_hdf(hdf=hdf, group_name=group_name) + self._structure_from_hdf() + self.custom_input.from_hdf(self._hdf5) + self._templates_directory = self.custom_input.templates_directory + self.custom_output.from_hdf(self._hdf5) diff --git a/pyiron_contrib/atomistics/ipi/ipi_jobs.py b/pyiron_contrib/atomistics/ipi/ipi_jobs.py index 8ca40a836..fe469b7f7 100644 --- a/pyiron_contrib/atomistics/ipi/ipi_jobs.py +++ b/pyiron_contrib/atomistics/ipi/ipi_jobs.py @@ -25,7 +25,7 @@ def __init__(self, project, job_name): def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep=1., temperature_damping_timescale=100., pressure_damping_timescale=1000., - n_ionic_steps=100, n_print=1, seed=None, port=31415, constrain_ids=None, *args, **kwargs): + n_ionic_steps=100, n_print=1, seed=None, port=31415, constrain_ids=None, bead_trajectory=False, *args, **kwargs): if not seed: random.seed(self.job_name) seed = random.randint(1, 99999) @@ -40,6 +40,7 @@ def calc_npt_md(self, temperature=300., pressure=101325e-9, n_beads=4, timestep= self.custom_input.seed = seed self.custom_input.port = port self.custom_input.constrain_ids = constrain_ids + self.custom_input.bead_trajectory = bead_trajectory def write_template_file(self): copy(self._templates_directory + '/pimd_template.xml', self.working_directory + '/pimd_template.xml') @@ -50,9 +51,12 @@ def ipi_write_helper(self, xml_filename): tree = ElementTree.parse(self.working_directory + '/' + xml_filename) root = tree.getroot() filepath = self.working_directory + '/ipi_input.xml' - for i in range(0, 3): + for i in range(0, 5): root[0][i].attrib['stride'] = str(self.custom_input.n_print) - root[0][3].attrib['stride'] = str(self.custom_input.n_ionic_steps) + root[0][5].attrib['stride'] = str(self.custom_input.n_ionic_steps) + if self.custom_input.bead_trajectory: + root[0][3].text = 'positions{angstrom}' + root[0][4].text = 'forces{ev/ang}' root[1].text = str(self.custom_input.n_ionic_steps) root[2][0].text = str(self.custom_input.seed) root[3][0].text = self.job_name