Merge pull request #1249 from pyiron/energy_free

Add AtomisticsMaintenance

pyiron_atomistics/project.py

@@ -17,6 +17,8 @@
     Creator as CreatorCore,
     deprecate,
 )
+from pyiron_base.state.logger import logger
+from pyiron_base.project.maintenance import Maintenance, LocalMaintenance
 
 try:
     from pyiron_base import ProjectGUI
@@ -29,11 +31,13 @@
 from pyiron_atomistics.atomistics.structure.periodic_table import PeriodicTable
 from pyiron_atomistics.lammps.potential import LammpsPotentialFile
 from pyiron_atomistics.vasp.potential import VaspPotential
+from pyiron_atomistics.vasp.base import _vasp_generic_energy_free_affected
 import pyiron_atomistics.atomistics.structure.pyironase as ase
 from pyiron_atomistics.atomistics.structure.atoms import Atoms
 from pyiron_atomistics.atomistics.structure.factory import StructureFactory
 from pyiron_atomistics.atomistics.master.parallel import pipe
 
+import numpy as np
 
 __author__ = "Joerg Neugebauer, Jan Janssen"
 __copyright__ = (
@@ -50,6 +54,95 @@
     raise AssertionError()
 
 
+def _vasp_energy_kin_affected(job):
+    e_kin = job.project_hdf5["output/generic/dft/scf_energy_kin"]
+    return e_kin is not None and not isinstance(e_kin, np.ndarray)
+
+
+class AtomisticsLocalMaintenance(LocalMaintenance):
+    def vasp_energy_pot_as_free_energy(
+        self, recursive: bool = True, progress: bool = True, **kwargs
+    ):
+        """
+        Ensure generic potential energy is the electronic free energy.
+
+        This ensures that the energy is consistent with the forces and stresses.
+        In version 0.1.0 of the Vasp job (pyiron_atomistics<=0.3.10) a combination of bugs in Vasp and pyiron caused the
+        potential energy reported to be the internal energy of electronic system extrapolated to zero smearing instead.
+
+        Args:
+            recursive (bool): search subprojects [True/False] - True by default
+            progress (bool): if True (default), add an interactive progress bar to the iteration
+            **kwargs (dict): Optional arguments for filtering with keys matching the project database column name
+                            (eg. status="finished"). Asterisk can be used to denote a wildcard, for zero or more
+                            instances of any character
+        """
+        kwargs["hamilton"] = "Vasp"
+        kwargs["status"] = "finished"
+
+        found_energy_kin_job = False
+
+        for job in self._project.iter_jobs(
+            recursive=recursive, progress=progress, convert_to_object=False, **kwargs
+        ):
+            if _vasp_generic_energy_free_affected(job):
+                job.project_hdf5["output/generic/energy_pot"] = np.array(
+                    [
+                        e[-1]
+                        for e in job.project_hdf5["output/generic/dft/scf_energy_free"]
+                    ]
+                )
+                found_energy_kin_job |= _vasp_energy_kin_affected(job)
+
+        if found_energy_kin_job:
+            logger.warn(
+                "Found at least one Vasp MD job with wrong kinetic energy.  Apply vasp_correct_energy_kin to fix!"
+            )
+
+    def vasp_correct_energy_kin(
+        self, recursive: bool = True, progress: bool = True, **kwargs
+    ):
+        """
+        Ensure kinetic and potential energy are correctly parsed for AIMD Vasp jobs.
+
+        Version 0.1.0 of the Vasp job (pyiron_atomistics<=0.3.10) incorrectly parsed the kinetic energy during MD runs,
+        such that it only reported the kinetic energy of the final ionic step and subtracted it from the electronic
+        energy instead of adding it.
+
+        Args:
+            recursive (bool): search subprojects [True/False] - True by default
+            progress (bool): if True (default), add an interactive progress bar to the iteration
+            **kwargs (dict): Optional arguments for filtering with keys matching the project database column name
+                            (eg. status="finished"). Asterisk can be used to denote a wildcard, for zero or more
+                            instances of any character
+        """
+        kwargs["hamilton"] = "Vasp"
+        kwargs["status"] = "finished"
+        for job in self._project.iter_jobs(
+            recursive=recursive, progress=progress, convert_to_object=False, **kwargs
+        ):
+            # only Vasp jobs of version 0.1.0 were affected
+            if job["HDF_VERSION"] != "0.1.0":
+                continue
+            if not _vasp_energy_kin_affected(job):
+                continue
+
+            job.decompress()
+            job = job.to_object()
+            job.status.collect = True
+            job.run()
+
+
+class AtomisticsMaintenance(Maintenance):
+    def __init__(self, project):
+        """
+        Args:
+            (project): pyiron project to do maintenance on
+        """
+        super().__init__(project=project)
+        self._local = AtomisticsLocalMaintenance(project)
+
+
 class Project(ProjectCore):
     """
     Welcome to pyiron! The project is the central class in pyiron, all other objects can be
@@ -140,6 +233,12 @@ def __init__(
         # TODO: instead of re-initialzing, auto-update pyiron_base creator with factories, like we update job class
         #  creation
 
+    @property
+    def maintenance(self):
+        if self._maintenance is None:
+            self._maintenance = AtomisticsMaintenance(self)
+        return self._maintenance
+
     def create_job(self, job_type, job_name, delete_existing_job=False):
         """
         Create one of the following jobs:

pyiron_atomistics/vasp/base.py

@@ -54,6 +54,20 @@
 __date__ = "Sep 1, 2017"
 
 
+def _vasp_generic_energy_free_affected(job):
+    """
+    Checks whether the value saved in output/generic/energy_pot matches the electronic free energy.
+    """
+    if job.project_hdf5.get("HDF_VERSION", "0.1.0") == "0.1.0":
+        energy_free = np.array(
+            [e[-1] for e in job.project_hdf5["output/generic/dft/scf_energy_free"]]
+        )
+        energy_pot = job.project_hdf5["output/generic/energy_pot"]
+        return not np.allclose(energy_free, energy_pot)
+    else:
+        return False
+
+
 class VaspBase(GenericDFTJob):
     """
     Class to setup and run and analyze VASP simulations which is a derivative of pyiron_atomistics.objects.job.generic.GenericJob.
@@ -98,6 +112,7 @@ def __init__(self, project, job_name):
         self._compress_by_default = True
         self.get_enmax_among_species = get_enmax_among_potentials
         state.publications.add(self.publication)
+        self.__hdf_version__ = "0.2.0"
 
     @property
     def structure(self):
@@ -767,6 +782,13 @@ def to_hdf(self, hdf=None, group_name=None):
         self._structure_to_hdf()
         self.input.to_hdf(self._hdf5)
         self._output_parser.to_hdf(self._hdf5)
+        if _vasp_generic_energy_free_affected(self):
+            self.logger.warn(
+                "Generic energy_pot does not match electronic free energy! "
+                "Generic energies is not consistent to generic forces and stress, "
+                "call project.maintenance.local.vasp_energy_pot_as_free_energy() "
+                "to correct generic energy!"
+            )
 
     def from_hdf(self, hdf=None, group_name=None):
         """
@@ -2051,15 +2073,33 @@ def collect(self, directory=os.getcwd(), sorted_indices=None):
             log_dict["forces"] = self.vp_new.vasprun_dict["forces"]
             log_dict["cells"] = self.vp_new.vasprun_dict["cells"]
             log_dict["volume"] = np.linalg.det(self.vp_new.vasprun_dict["cells"])
-            # log_dict["total_energies"] = self.vp_new.vasprun_dict["total_energies"]
-            log_dict["energy_tot"] = self.vp_new.vasprun_dict["total_energies"]
+            # The vasprun parser also returns the energies printed again after the final SCF cycle under the key
+            # "total_energies", but due to a bug in the VASP output, the energies reported there are wrong in Vasp 5.*;
+            # instead use the last energy from the scf cycle energies
+            # BUG link: https://ww.vasp.at/forum/viewtopic.php?p=19242
+            try:
+                # bug report is not specific to which Vasp5 versions are affected; be safe and workaround for all of
+                # them
+                is_vasp5 = self.vp_new.vasprun_dict["generator"]["version"].startswith(
+                    "5."
+                )
+            except KeyError:  # in case the parser didn't read the version info
+                is_vasp5 = True
+            if is_vasp5:
+                log_dict["energy_pot"] = np.array(
+                    [e[-1] for e in self.vp_new.vasprun_dict["scf_fr_energies"]]
+                )
+            else:
+                # total energies refers here to the total energy of the electronic system, not the total system of
+                # electrons plus (potentially) moving ions; hence this is the energy_pot
+                log_dict["energy_pot"] = self.vp_new.vasprun_dict["total_fr_energies"]
             if "kinetic_energies" in self.vp_new.vasprun_dict.keys():
-                log_dict["energy_pot"] = (
-                    log_dict["energy_tot"]
-                    - self.vp_new.vasprun_dict["kinetic_energies"]
+                log_dict["energy_tot"] = (
+                    log_dict["energy_pot"]
+                    + self.vp_new.vasprun_dict["kinetic_energies"]
                 )
             else:
-                log_dict["energy_pot"] = log_dict["energy_tot"]
+                log_dict["energy_tot"] = log_dict["energy_pot"]
             log_dict["steps"] = np.arange(len(log_dict["energy_tot"]))
             log_dict["positions"] = self.vp_new.vasprun_dict["positions"]
             log_dict["forces"][:, sorted_indices] = log_dict["forces"].copy()
@@ -2209,9 +2249,13 @@ def collect(self, directory=os.getcwd(), sorted_indices=None):
                 self.vp_new.vasprun_dict["parameters"]["electronic"]["NELECT"]
             )
             if "kinetic_energies" in self.vp_new.vasprun_dict.keys():
+                # scf_energy_kin is for backwards compatibility
                 self.generic_output.dft_log_dict[
                     "scf_energy_kin"
                 ] = self.vp_new.vasprun_dict["kinetic_energies"]
+                self.generic_output.dft_log_dict[
+                    "energy_kin"
+                ] = self.vp_new.vasprun_dict["kinetic_energies"]
 
         if (
             "LOCPOT" in files_present

pyiron_atomistics/vasp/vasprun.py

@@ -79,6 +79,7 @@ def parse_root_to_dict(self, filename):
         d["total_energies"] = list()
         d["total_fr_energies"] = list()
         d["total_0_energies"] = list()
+        d["kinetic_energies"] = list()
         d["stress_tensors"] = list()
         for _, leaf in ETree.iterparse(filename):
             if leaf.tag in ["generator", "incar"]:
@@ -115,6 +116,10 @@ def parse_root_to_dict(self, filename):
         d["total_energies"] = np.array(d["total_energies"])
         d["total_fr_energies"] = np.array(d["total_fr_energies"])
         d["total_0_energies"] = np.array(d["total_0_energies"])
+        if len(d["kinetic_energies"]) > 0:
+            d["kinetic_energies"] = np.array(d["kinetic_energies"])
+        else:
+            del d["kinetic_energies"]
         d["scf_energies"] = d["scf_energies"]
         d["scf_dipole_moments"] = d["scf_dipole_moments"]
         d["scf_fr_energies"] = d["scf_fr_energies"]
@@ -441,7 +446,7 @@ def parse_calc_to_dict(self, node, d):
                     if i.attrib["name"] == "e_0_energy":
                         d["total_0_energies"].append(float(i.text))
                     if i.attrib["name"] == "kinetic":
-                        d["kinetic_energies"] = float(i.text)
+                        d["kinetic_energies"].append(float(i.text))
             if item.tag == "eigenvalues":
                 self.parse_eigenvalues_to_dict(item, d)
 

tests/table/test_datamining.py

@@ -67,14 +67,14 @@ def filter_job_type_db(job_table):
         self.assertEqual(df["majority_element"].values[0], "Fe")
         self.assertEqual(df["minority_element_list"].values[0], [])
         self.assertEqual(df["job_name"].values[0], "full_job_sample")
-        self.assertEqual(df["energy_tot"].values[0], -17.7331698)
-        self.assertEqual(df["energy_free"].values[0], -17.7379867884)
-        self.assertEqual(df["energy_int"].values[0], -17.72353582)
-        self.assertEqual(df["alat"].values[0], 0.0)
+        self.assertAlmostEqual(df["energy_tot"].values[0], -17.73798679)
+        self.assertAlmostEqual(df["energy_free"].values[0], -17.73798679)
+        self.assertAlmostEqual(df["energy_int"].values[0], -17.72353582)
+        self.assertAlmostEqual(df["alat"].values[0], 0.0)
         self.assertEqual(df["magnetic_structure"].values[0], "ferro-magnetic")
-        self.assertEqual(df["avg. plane waves"].values[0], 196.375)
-        self.assertEqual(df["energy_tot_wo_kin_corr"].values[0], -17.6003698)
-        self.assertTrue(np.isclose(df["volume"].values[0], 21.95199999999999))
+        self.assertAlmostEqual(df["avg. plane waves"].values[0], 196.375)
+        self.assertAlmostEqual(df["energy_tot_wo_kin_corr"].values[0], -17.60518679)
+        self.assertAlmostEqual(df["volume"].values[0], 21.95199999999999)
 
 
 def get_alat(job):