VASP forces from OUTCAR

CHANGELOG.md

@@ -8,6 +8,7 @@ we hit release version 1.0.0.
 ## [X.Y.Z] - YYYY-MM-DD
 
 ### Added
+- added `read_trajectory` to read cell vectors, atomic positions, and forces from VASP OUTCAR
 - slicing io files multiple output (still WIP), see #584 for details
   Intention is to have all methods use this method for returning
   multiple values, it should streamline the API.
@@ -43,6 +44,8 @@ we hit release version 1.0.0.
 - `BrillouinZone.merge` allows simple merging of several objects, #537
 
 ### Changed
+- `stdoutSileVASP` will not accept `all=` arguments
+- `stdoutSileVASP.read_energy` returns as default the next item (no longer the last)
 - `txtSileOrca` will not accept `all=` arguments, see #584
 - `stdoutSileOrca` will not accept `all=` arguments, see #584
 - `xyzSile` out from sisl will now default to the extended xyz file-format

src/sisl/io/__init__.py

@@ -141,7 +141,7 @@
    eigenvalSileVASP
    chgSileVASP
    locpotSileVASP
-   outSileVASP
+   stdoutSileVASP
 
 
 Wannier90

src/sisl/io/vasp/stdout.py

@@ -1,14 +1,16 @@
 # This Source Code Form is subject to the terms of the Mozilla Public
 # License, v. 2.0. If a copy of the MPL was not distributed with this
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
+import numpy as np
+
 from .sile import SileVASP
 from ..sile import add_sile, sile_fh_open
+from .._multiple import SileBinder
 
-from sisl.messages import deprecation
+from sisl.messages import deprecation, deprecate_argument
 from sisl.utils import PropertyDict
 from sisl._internal import set_module
 
-
 __all__ = ["stdoutSileVASP", "outSileVASP"]
 
 
@@ -55,9 +57,14 @@ def accuracy_reached(self):
         """ True if the line "reached required accuracy" was found. """
         return self.step_to("reached required accuracy")[0]
 
+    @SileBinder()
     @sile_fh_open()
-    def read_energy(self, all=False):
-        """ Reads the energy specification from OUTCAR and returns energy dictionary in units of eV
+    @deprecate_argument("all", None, "use read_energy[:]() instead to get all entries", from_version="0.14")
+    @deprecation("WARNING: direct calls to stdoutSileVASP.read_energy() no longer returns the last entry! Now the next block on file is returned.", from_version="0.14")
+    def read_energy(self):
+        """ Reads an energy specification block from OUTCAR
+
+        The function steps to the next occurrence of the "Free energy of the ion-electron system" segment
 
         Notes
         -----
@@ -79,15 +86,11 @@ def read_energy(self, all=False):
             energy without entropy= total
             energy(sigma->0)      = sigma0
 
-        Parameters
-        ----------
-        all: bool, optional
-            return a list of energy dictionaries from each step
-
         Returns
         -------
-        PropertyDict : all energies from the "Free energy of the ion-electron system" segment of VASP output
+        PropertyDict : all energies from a single "Free energy of the ion-electron system" segment
         """
+
         name_conv = {
             "alpha": "Z",
             "Ewald": "Ewald",
@@ -101,49 +104,72 @@ def read_energy(self, all=False):
             "Solvation": "solvation",
         }
 
-        def readE(itt):
-            nonlocal name_conv
-            # read the energy tables
-            f = self.step_to("Free energy of the ion-electron system", allow_reread=False)[0]
-            if not f:
-                return None
-            next(itt) # -----
-            line = next(itt)
-            E = PropertyDict()
-            while "----" not in line:
-                key, *v = line.split()
-                if key == "PAW":
-                    E[f"{name_conv[key]}1"] = float(v[-2])
-                    E[f"{name_conv[key]}2"] = float(v[-1])
-                else:
-                    E[name_conv[key]] = float(v[-1])
-                line = next(itt)
-            line = next(itt)
-            E["free"] = float(line.split()[-2])
-            next(itt)
-            line = next(itt)
+        # read the energy tables
+        f = self.step_to("Free energy of the ion-electron system", allow_reread=False)[0]
+        if not f:
+            return None
+        self.readline() # -----
+        line = self.readline()
+        E = PropertyDict()
+        while "----" not in line:
+            key, *v = line.split()
+            if key == "PAW":
+                E[f"{name_conv[key]}1"] = float(v[-2])
+                E[f"{name_conv[key]}2"] = float(v[-1])
+            else:
+                E[name_conv[key]] = float(v[-1])
+            line = self.readline()
+        line = self.readline()
+        E.free = float(line.split()[-2])
+        self.readline()
+        line = self.readline()
+        v = line.split()
+        E.total = float(v[4])
+        E.sigma0 = float(v[-1])
+        return E
+
+    @SileBinder()
+    @sile_fh_open()
+    def read_trajectory(self):
+        """ Reads cell+position+force data from OUTCAR for an ionic trajectory step
+
+        The function steps to the block defined by the "VOLUME and BASIS-vectors are now :"
+        line to first read the cell vectors, then it steps to the "TOTAL-FORCE (eV/Angst)" segment
+        to read the atom positions and forces.
+
+        Returns
+        -------
+        PropertyDict : Trajectory step defined by cell vectors (`.cell`), atom positions (`.xyz`), and forces (`.force`)
+        """
+
+        f = self.step_to("VOLUME and BASIS-vectors are now :", allow_reread=False)[0]
+        if not f:
+            return None
+        for i in range(4):
+            self.readline() # skip 4 lines
+        C = []
+        for i in range(3):
+            line = self.readline()
+            v = line.split()
+            C.append(v[:3]) # direct lattice vectors
+        # read a position-force table
+        f = self.step_to("TOTAL-FORCE (eV/Angst)", allow_reread=False)[0]
+        if not f:
+            return None
+        self.readline() # -----
+        P, F = [], []
+        line = self.readline()
+        while "----" not in line:
             v = line.split()
-            E["total"] = float(v[4])
-            E["sigma0"] = float(v[-1])
-            return E
-
-        itt = iter(self)
-        E = []
-        e = readE(itt)
-        while e is not None:
-            E.append(e)
-            e = readE(itt)
-        try:
-            # this just puts the job_completed flag. But otherwise not used
-            self.cpu_time()
-        except Exception:
-            pass
-
-        if all:
-            return E
-        if len(E) > 0:
-            return E[-1]
-        return None
+            # positions and forces
+            P.append(v[:3])
+            F.append(v[3:6])
+            line = self.readline()
+        step = PropertyDict()
+        step.cell = np.array(C, dtype=np.float64)
+        step.xyz = np.array(P, dtype=np.float64)
+        step.force = np.array(F, dtype=np.float64)
+        return step
 
 
 outSileVASP = deprecation("outSileVASP has been deprecated in favor of stdoutSileVASP.", "0.15")(stdoutSileVASP)

src/sisl/io/vasp/tests/test_stdout.py

@@ -3,7 +3,7 @@
 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
 import pytest
 import os.path as osp
-from sisl.io.vasp.stdout import *
+from sisl.io.vasp.stdout import stdoutSileVASP
 import numpy as np
 
 pytestmark = [pytest.mark.io, pytest.mark.vasp]
@@ -14,9 +14,14 @@ def test_diamond_outcar_energies(sisl_files):
     f = sisl_files(_dir, 'diamond', 'OUTCAR')
     f = stdoutSileVASP(f)
 
-    E = f.read_energy()
-    Eall = f.read_energy(all=True)
+    E0 = f.read_energy()
+    E = f.read_energy[-1]()
+    Eall = f.read_energy[:]()
 
+    assert E0.sigma0 == 0.8569373 # first block
+    assert E.sigma0 == -18.18677613 # last block
+
+    assert E0 == Eall[0]
     assert E == Eall[-1]
     assert len(Eall) > 1
     assert f.completed()
@@ -35,3 +40,60 @@ def test_diamond_outcar_completed(sisl_files):
     f = stdoutSileVASP(f)
 
     assert f.completed()
+
+
+def test_diamond_outcar_trajectory(sisl_files):
+    f = sisl_files(_dir, 'diamond', 'OUTCAR')
+    f = stdoutSileVASP(f)
+
+    step = f.read_trajectory()
+
+    assert step.xyz[0, 1] == 0.0
+    assert step.xyz[1, 1] == 0.89250
+    assert step.force[0, 1] == 0.0
+    assert step.force[1, 0] == 0.0
+
+    traj = f.read_trajectory[:]()
+    assert len(traj) == 1
+
+
+def test_graphene_relax_outcar_trajectory(sisl_files):
+    f = sisl_files(_dir, 'graphene_relax', 'OUTCAR')
+    f = stdoutSileVASP(f)
+
+    step = f.read_trajectory[9]()
+    assert step.cell[0, 0] == 2.462060590
+    assert step.cell[1, 0] == -1.231030295
+    assert step.cell[2, 2] == 9.804915686
+    assert step.force[0, 2] == -0.006138
+    assert step.force[1, 2] == 0.006138
+
+    traj = f.read_trajectory[:]()
+    assert len(traj) == 10
+    assert traj[0].cell[0, 0] == 2.441046239
+    assert traj[0].xyz[0, 2] == 0.00000
+    assert traj[0].xyz[1, 2] == 0.20000
+    assert traj[0].force[0, 2] == 3.448038
+    assert traj[0].force[1, 2] == -3.448038
+    assert traj[-1].cell[2, 2] == 9.804915686
+    assert traj[-1].xyz[1, 2] == -0.00037
+    assert traj[-1].force[0, 2] == -0.006138
+    assert traj[-1].force[1, 2] == 0.006138
+
+
+def test_graphene_md_outcar_trajectory(sisl_files):
+    f = sisl_files(_dir, 'graphene_md', 'OUTCAR')
+    f = stdoutSileVASP(f)
+
+    step = f.read_trajectory[99]()
+    assert step.xyz[0, 0] == 0.09703
+    assert step.force[0, 2] == -0.278082
+
+    traj = f.read_trajectory[:]()
+    assert len(traj) == 100
+    assert traj[0].xyz[0, 0] == 0.12205
+    assert traj[0].xyz[1, 0] == 0.09520
+    assert traj[0].force[0, 1] == -0.017766
+    assert traj[0].force[1, 1] == 0.017766
+    assert traj[-1].xyz[0, 0] == 0.09703
+    assert traj[-1].force[0, 2] == -0.278082