Skip to content

Commit

Permalink
Merge 1c2cb94 into a68a509
Browse files Browse the repository at this point in the history
  • Loading branch information
Michael Ashton committed May 22, 2020
2 parents a68a509 + 1c2cb94 commit d6afd39
Show file tree
Hide file tree
Showing 6 changed files with 277 additions and 8 deletions.
82 changes: 74 additions & 8 deletions pyiron/sphinx/base.py
Expand Up @@ -24,6 +24,7 @@
from pyiron.sphinx.potential import SphinxJTHPotentialFile
from pyiron.sphinx.potential import find_potential_file \
as find_potential_file_jth
from pyiron.sphinx.volumetric_data import SphinxVolumetricData
from pyiron.base.settings.generic import Settings
from pyiron.base.generic.parameters import GenericParameters

Expand Down Expand Up @@ -438,7 +439,7 @@ def load_basis_group(self):
"""
self.input.sphinx.basis.setdefault("eCut", self.input["EnCut"]/RYDBERG_TO_EV)
self.input.sphinx.basis.setdefault("kPoint", Group())
if "KpointCoords" is not None:
if self.input["KpointCoords"] is not None:
self.input.sphinx.basis.kPoint.setdefault(
"coords", np.array(self.input["KpointCoords"])
)
Expand Down Expand Up @@ -1242,6 +1243,42 @@ def write_input(self):
f.write("include <parameters.sx>;\n\n")
f.write(self.input.sphinx.to_sphinx(indent=0))

def get_charge_density(self):
"""
Gets the charge density from the hdf5 file. This value is normalized by the volume
Returns:
pyiron.atomistics.volumetric.generic.VolumetricData
"""
if self.status not in [
"finished", "warning", "not_converged"
]:
return
else:
with self.project_hdf5.open("output") as ho:
cd_obj = SphinxVolumetricData()
cd_obj.from_hdf(ho, "charge_density")
cd_obj.atoms = self.get_structure(-1)
return cd_obj

def get_electrostatic_potential(self):
"""
Gets the electrostatic potential from the hdf5 file.
Returns:
pyiron.atomistics.volumetric.generic.VolumetricData
"""
if self.status not in [
"finished", "warning", "not_converged"
]:
return
else:
with self.project_hdf5.open("output") as ho:
es_obj = SphinxVolumetricData()
es_obj.from_hdf(ho, "electrostatic_potential")
es_obj.atoms = self.get_structure(-1)
return es_obj

def collect_output(self, force_update=False):
"""
Collects the outputs and stores them to the hdf file
Expand Down Expand Up @@ -1399,14 +1436,8 @@ def validate_ready_to_run(self):
}

if np.any([len(all_groups[group]) == 0 for group in all_groups]):
# warnings.warn("The following input groups have not been loaded:\n"
# + ", ".join([
# g for g in all_groups if len(all_groups[g]) == 0
# ])
# + "\n"
# + "They are set to default values."
# )
self.load_default_groups()

if self.structure is None:
raise AssertionError(
"Structure not set; set it via job.structure = "
Expand Down Expand Up @@ -1825,6 +1856,8 @@ class Output(object):
def __init__(self, job):
self._job = job
self._parse_dict = defaultdict(list)
self.charge_density = SphinxVolumetricData()
self.electrostatic_potential = SphinxVolumetricData()

@staticmethod
def splitter(arr, counter):
Expand Down Expand Up @@ -2241,6 +2274,26 @@ def collect_relaxed_hist(self, file_name="relaxHist.sx", cwd=None):
* BOHR_TO_ANGSTROM
)

def collect_charge_density(self, file_name, cwd):
if (
file_name in os.listdir(cwd)
and os.stat(posixpath.join(cwd, file_name)).st_size != 0
):
self.charge_density.from_file(
filename=posixpath.join(cwd, file_name),
normalize=True
)

def collect_electrostatic_potential(self, file_name, cwd):
if (
file_name in os.listdir(cwd) and
os.stat(posixpath.join(cwd, file_name)).st_size != 0
):
self.electrostatic_potential.from_file(
filename=posixpath.join(cwd, file_name),
normalize=False
)

def collect(self, directory=os.getcwd()):
"""
The collect function, collects all the output from a Sphinx simulation.
Expand All @@ -2256,6 +2309,10 @@ def collect(self, directory=os.getcwd()):
self.collect_energy_struct(file_name="energy-structOpt.dat",
cwd=directory)
self.collect_relaxed_hist(file_name="relaxHist.sx", cwd=directory)
self.collect_electrostatic_potential(file_name="vElStat-eV.sxb",
cwd=directory)
self.collect_charge_density(file_name="rho.sxb",
cwd=directory)
self._job.compress()

def to_hdf(self, hdf, force_update=False):
Expand Down Expand Up @@ -2298,6 +2355,14 @@ def to_hdf(self, hdf, force_update=False):
hdf5_sphinx["bands_eigen_values_initial"] = self._parse_dict[
"bands_eigen_values_initial"
]
if self.electrostatic_potential.total_data is not None:
self.electrostatic_potential.to_hdf(
hdf5_output, group_name="electrostatic_potential"
)
if self.charge_density.total_data is not None:
self.charge_density.to_hdf(
hdf5_output, group_name="charge_density"
)
with hdf5_output.open("generic") as hdf5_generic:
if "dft" not in hdf5_generic.list_groups():
hdf5_generic.create_group("dft")
Expand Down Expand Up @@ -2364,6 +2429,7 @@ def to_hdf(self, hdf, force_update=False):
hdf5_generic["cells"] = np.array(
[self._job.structure.cell])


def from_hdf(self, hdf):
"""
Load output from an HDF5 file
Expand Down
189 changes: 189 additions & 0 deletions pyiron/sphinx/volumetric_data.py
@@ -0,0 +1,189 @@
# 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.

import numpy as np
import scipy
from scipy.io.netcdf import netcdf_file
import os
from pyiron.base.settings.generic import Settings
from pyiron.atomistics.volumetric.generic import VolumetricData

__author__ = "Sudarsan Surendralal"
__copyright__ = (
"Copyright 2020, Max-Planck-Institut für Eisenforschung GmbH - "
"Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Sudarsan Surendralal"
__email__ = "surendralal@mpie.de"
__status__ = "production"
__date__ = "Sep 1, 2017"


BOHR_TO_ANGSTROM = (
scipy.constants.physical_constants["Bohr radius"][0] /
scipy.constants.angstrom
)


class SphinxVolumetricData(VolumetricData):
"""
General class for parsing and manipulating volumetric static within Sphinx.
The basic idea of the Base class is adapted from the pymatgen vasp
VolumetricData class:
http://pymatgen.org/_modules/pymatgen/io/vasp/outputs.html#VolumetricData
"""


def __init__(self):
super(SphinxVolumetricData, self).__init__()
self.atoms = None
self._diff_data = None
self._total_data = None


def from_file(self, filename, normalize=True):
"""
Parses volumetric data from a sphinx binary (.sxb) file.
Args:
filename (str): Path of file to parse
normalize (boolean): Flag to normalize by the volume of the cell
"""
try:
vol_data_list = self._read_vol_data(
filename=filename,
normalize=normalize
)
except (ValueError, IndexError, TypeError):
raise ValueError("Unable to parse file: {}".format(filename))
self._total_data = vol_data_list[0]
if len(vol_data_list) == 2:
self._total_data = vol_data_list[0] + vol_data_list[1]
self._diff_data = vol_data_list[0] - vol_data_list[1]


@staticmethod
def _read_vol_data(filename, normalize=True):
"""
Parses the Sphinx volumetric data files (rho.sxb and vElStat-eV.sxb).
Args:
filename (str): File to be parsed
normalize (bool): Normalize the data with respect to the volume
(probably sensible for rho)
Returns:
list: A list of the volumetric data (length >1 for density
files with spin)
"""
if not os.path.getsize(filename) > 0:
s = Settings()
s.logger.warning("File:" + filename + "seems to be empty! ")
return None, None

with netcdf_file(filename, mmap=False) as f:
dim = [int(d) for d in f.variables["dim"]]
volume = 1.0
if normalize:
cell = f.variables["cell"].data * BOHR_TO_ANGSTROM
volume = np.abs(np.linalg.det(cell))
if "mesh" in f.variables:
# non-spin polarized
total_data_list = [
np.array(f.variables["mesh"][:]).reshape(dim) / volume
]
elif "mesh-0" in f.variables and "mesh-1" in f.variables:
# spin-polarized
total_data_list = [
np.array(f.variables["mesh-0"][:]).reshape(dim) / volume,
np.array(f.variables["mesh-1"][:]).reshape(dim) / volume
]
else:
raise ValueError(
"Unexpected keys in the netcdf file's variables: neither "
f"'mesh' nor 'mesh-0' and 'mesh-1' found in {f.variables}."
)

if len(total_data_list) == 0:
s = Settings()
s.logger.warning(
"File:"
+ filename
+ "seems to be corrupted/empty even after parsing!"
)
return None

return total_data_list


@property
def total_data(self):
"""
numpy.ndarray: Total volumtric data (3D)
"""
return self._total_data


@total_data.setter
def total_data(self, val):
self._total_data = val


@property
def diff_data(self):
"""
numpy.ndarray: Volumtric difference data (3D)
"""
return self._diff_data


@diff_data.setter
def diff_data(self, val):
self._diff_data = val


def to_hdf(self, hdf5, group_name="volumetric_data"):
"""
Writes the data as a group to a HDF5 file
Args:
hdf5 (pyiron.base.generic.hdfio.ProjectHDFio): The
HDF file/path to write the data
group_name (str): The name of the group under which
the data must be stored
"""
with hdf5.open(group_name) as hdf_vd:
hdf_vd["TYPE"] = str(type(self))
hdf_vd["total"] = self.total_data
if self.diff_data is not None:
hdf_vd["diff"] = self.diff_data


def from_hdf(self, hdf5, group_name="volumetric_data"):
"""
Extract a VolumetricData instance from an HDF5 file.
Args:
hdf5 (pyiron.base.generic.hdfio.ProjectHDFio): The HDF
file/path to read the data
group_name (str): The name of the group under which
the data have been stored
Returns:
pyiron.atomistics.volumetric.generic.VolumetricData: The
VolumetricData instance
"""
with hdf5.open(group_name) as hdf_vd:
self._total_data = hdf_vd["total"]
if len(self._total_data) == 2:
self.is_spin_polarized = True
if "diff" in hdf_vd.list_nodes():
self._diff_data = hdf_vd["diff"]
5 changes: 5 additions & 0 deletions tests/sphinx/test_base.py
Expand Up @@ -527,6 +527,11 @@ def test_collect_2_5(self):
len(output._parse_dict[list_one]), len(output._parse_dict[list_two])
)

rho = self.sphinx_2_5._output_parser.charge_density
vel = self.sphinx_2_5._output_parser.electrostatic_potential
self.assertIsNotNone(rho.total_data)
self.assertIsNotNone(vel.total_data)

def test_check_band_occupancy(self):
self.sphinx_2_5.collect_output()
self.assertTrue(self.sphinx_2_5.output.check_band_occupancy())
Expand Down
@@ -0,0 +1,9 @@
cell = [[5.351704387296569, 0.0, 0.0], [3.2769738239427953e-16, 5.351704387296569, 0.0], [3.2769738239427953e-16, 3.2769738239427953e-16, 5.351704387296569]];
species {
element = "Fe";
atom { label = "spin_2"; coords = [0.03779452250915656, 0.0, 0.0]; }
}
species {
element = "Ni";
atom { label = "spin_2"; coords = [2.6758521936482844, 2.6758521936482844, 2.6758521936482844]; }
}
Binary file not shown.
Binary file not shown.

0 comments on commit d6afd39

Please sign in to comment.