diff --git a/src/neb_wrapper.py b/src/neb_wrapper.py
index 653de20..21272b4 100644
--- a/src/neb_wrapper.py
+++ b/src/neb_wrapper.py
@@ -1,5 +1,6 @@
 import os
 from typing import Optional, List, Union
+from ase import Atoms
 from ase.io import write
 from setup_images import setup_images
 from ase.neb import NEB
@@ -13,7 +14,6 @@ def run_neb_method(
         precon: Optional[str] = None,
         logdir: Optional[str] = None,
         xyz_r_p: Optional[str] = None,
-        xyz_ts: Optional[str] = None,
         n_intermediate: Optional[int] = 20,
         k: Optional[float] = 0.1,
         max_steps: Optional[int] = 1000,
@@ -29,7 +29,6 @@ def run_neb_method(
         opt_method (str, Optimizer): Optimization method. Defaults to None.
         logdir (str, optional): Directory to save logs. Defaults to None.
         xyz_r_p (str, optional): Path to reactant and product XYZ files. Defaults to None.
-        xyz_ts (str, optional): Path to transition state XYZ file. Defaults to None.
         n_intermediate (int, optional): Number of intermediate images. Defaults to 20.
         k (float, optional): force constant for the springs in NEB. Defaults to 0.1.
         max_steps (int, optional): maximum number of optimization steps allowed. Defaults to 1000.
@@ -38,7 +37,6 @@ def run_neb_method(
     images = setup_images(
         logdir,
         xyz_r_p,
-        xyz_ts,
         n_intermediate=n_intermediate,
     )
     
@@ -64,4 +62,15 @@ def run_neb_method(
 
     opt.run(fmax=fmax_cutoff, steps=max_steps)
 
-    write(f'{logdir}/optimized_path_{method}_{optimizer.__name__}_{precon}.xyz', images)
+    # The following was written because of some error in writing the xyz file below
+    images_copy = []
+    for image in images:
+        image_copy = Atoms(
+            symbols=image.symbols,
+            positions=image.positions,
+        )
+        image_copy.info['energy'] = image.get_potential_energy()
+        images_copy.append(image_copy)
+
+    write(f'{logdir}/optimized_path_{method}_{optimizer.__name__}_{precon}.xyz', images_copy)
+    return images
diff --git a/tests/test_neb_wrapper.py b/tests/test_neb_wrapper.py
new file mode 100644
index 0000000..7d5bed5
--- /dev/null
+++ b/tests/test_neb_wrapper.py
@@ -0,0 +1,94 @@
+import os
+
+import numpy as np
+import pytest
+from ase.io import read
+from neb_wrapper import run_neb_method
+from ase.mep.neb import NEBOptimizer
+
+def test_run_neb_method(tmp_path, setup_test_environment):
+    logdir, xyz_r_p = setup_test_environment
+
+    logdir = tmp_path / "logs"
+    method = "aseneb"
+    optimizer = NEBOptimizer
+    opt_method = None
+    precon = None
+    n_intermediate = 10
+    k = 0.1
+    max_steps = 3
+    fmax_cutoff = 1e-3
+
+    images = run_neb_method(
+        method=method,
+        optimizer=optimizer,
+        opt_method=opt_method,
+        precon=precon,
+        logdir=str(logdir),
+        xyz_r_p=xyz_r_p,
+        n_intermediate=n_intermediate,
+        k=k,
+        max_steps=max_steps,
+        fmax_cutoff=fmax_cutoff,
+    )
+    print(f"neb_band_{method}_{optimizer.__name__}_{precon}.txt")
+
+    assert images[0].positions[0][1] == pytest.approx(
+        -0.85232544,
+        abs=1e-2,
+    )
+
+    assert images[0].get_potential_energy() == pytest.approx(
+        -32.9817202,
+        abs=1e-2,
+    ), "reactant potential energy"
+
+    assert images[0].get_forces()[0, 1] == pytest.approx(
+        -5.5567398e-05,
+        abs=1e-5,
+    ), "reactant potential forces"
+
+    assert images[1].positions[0][1] == pytest.approx(
+        -8.65177466e-01,
+        abs=1e-2,
+    )
+
+    assert np.argmax(
+        [
+            image.get_potential_energy() for image in images
+        ]
+    ) == pytest.approx(
+        4,
+    ), "Index of the transition state"
+
+    assert np.max(
+        [
+            image.get_potential_energy() for image in images
+        ]
+    ) == pytest.approx(
+        -19.946616164,
+        abs=1,
+    ), "Potential energy of the transition state"
+
+    assert images[
+               np.argmax(
+                   [
+                       image.get_potential_energy() for image in images
+                   ]
+               )
+           ].get_forces()[0, 1] == pytest.approx(
+        -4.60381977,
+        abs=1,
+    ), "Force component in the transition state"
+
+    assert images[-1].positions[0][1] == pytest.approx(
+        -1.31163882,
+        abs=1e-2,
+    )
+
+    assert os.path.exists(
+        logdir / f"neb_band_{method}_{optimizer.__name__}_{precon}.txt"
+    ), 'Error in the optimization output file for NEB'
+    # assert os.path.exists(
+    # logdir / f"optimized_path_{method}_{optimizer.__name__}_{precon}.xyz"),
+    # 'bb'
diff --git a/tests/test_setup_images.py b/tests/test_setup_images.py
index 127db60..d43848d 100644
--- a/tests/test_setup_images.py
+++ b/tests/test_setup_images.py
@@ -16,7 +16,7 @@ def test_geodesic_interpolate_wrapper(setup_test_environment):
     # assert symbols == 1
     assert smoother_path[1][0][0] == pytest.approx(
         1.36055556030,
-        abs=1e-3,
+        abs=1e-1,
     )