added read_trajectory for VASP 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.

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.utils import PropertyDict
 from sisl._internal import set_module
 
-
 __all__ = ["stdoutSileVASP", "outSileVASP"]
 
 
@@ -146,6 +148,39 @@ def readE(itt):
         return None
 
 
+    @SileBinder()
+    @sile_fh_open()
+    def read_trajectory(self):
+        f = self.step_to("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()
+            # 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)
 
 add_sile("OUTCAR", stdoutSileVASP, gzip=True)

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]
@@ -35,3 +35,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