diff --git a/mofa/scoring/geometry.py b/mofa/scoring/geometry.py
index 05e74f0b..15dc0f3c 100644
--- a/mofa/scoring/geometry.py
+++ b/mofa/scoring/geometry.py
@@ -34,9 +34,11 @@ def score_mof(self, record: MOFRecord) -> float:
 
         # Get the initial and final structures
         init_strc = read_vasp(StringIO(traj[0]))
-        final_strc: ase.Atoms = read_vasp(StringIO(traj[0]))
+        final_strc: ase.Atoms = read_vasp(StringIO(traj[1]))
 
         # Compute the maximum principal strain
-        strain = np.linalg.solve(final_strc.cell.array, init_strc.cell.array) - np.eye(3)
+        #  Following: https://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-strain
+        strain = np.matmul(init_strc.cell.array, np.linalg.inv(final_strc.cell.array)) - np.eye(3)
+        strain = 0.5 * (strain + strain.T)
         strains = np.linalg.eigvals(strain)
         return np.abs(strains).max()
diff --git a/tests/scoring/test_geometry.py b/tests/scoring/test_geometry.py
index 9a67aa01..73156d4a 100644
--- a/tests/scoring/test_geometry.py
+++ b/tests/scoring/test_geometry.py
@@ -34,4 +34,5 @@ def test_strain(example_record):
     write_vasp(sio, final_atoms)
     example_record.md_trajectory['uff'][1] = sio.getvalue()
 
-    assert scorer.score_mof(example_record) > 0
+    max_strain = scorer.score_mof(example_record)
+    assert np.isclose(max_strain, 0.09647)  # Checked against https://www.cryst.ehu.es/cryst/strain.html