diff --git a/pyiron/sphinx/base.py b/pyiron/sphinx/base.py
index 6a435ec68..d852f5437 100644
--- a/pyiron/sphinx/base.py
+++ b/pyiron/sphinx/base.py
@@ -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
 
@@ -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"])
                 )
@@ -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
@@ -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 = "
@@ -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):
@@ -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.
@@ -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):
@@ -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")
@@ -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
diff --git a/pyiron/sphinx/volumetric_data.py b/pyiron/sphinx/volumetric_data.py
new file mode 100644
index 000000000..1d28a076a
--- /dev/null
+++ b/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"]
diff --git a/tests/sphinx/test_base.py b/tests/sphinx/test_base.py
index 91dd5e6ea..1e92d3ece 100644
--- a/tests/sphinx/test_base.py
+++ b/tests/sphinx/test_base.py
@@ -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())
diff --git a/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/relaxedStr.sx b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/relaxedStr.sx
new file mode 100644
index 000000000..74d2eb966
--- /dev/null
+++ b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/relaxedStr.sx
@@ -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]; }
+}
diff --git a/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/rho.sxb b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/rho.sxb
new file mode 100644
index 000000000..f3c6c211e
Binary files /dev/null and b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/rho.sxb differ
diff --git a/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/vElStat-eV.sxb b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/vElStat-eV.sxb
new file mode 100644
index 000000000..8d48cbb6a
Binary files /dev/null and b/tests/static/sphinx/sphinx_test_2_5_hdf5/sphinx_test_2_5/vElStat-eV.sxb differ