Repository: Quantum-Accelerators/quacc
Pull request #2176 (2024-05-24 02:50:24)

Added ase Nudged Elastic Band in the NewtonNet recipe
## Summary of Changes

\>>
Target here is to add ase's NEB as the double ended transition state search method for NewtonNet recipe.

- I need to make sure that geodesic code is added as a preprocessor for the double ended transition state search.
- NEB will take the output of Geodesic
- Sella will them run from the top point of NEB.
- Tests will need to be added for all the files keeping the coverage of more than 99% before taking it for review.
- Currently, there is a lot of unstructured code in the ts.py and in the test. This will need to be cleaned up to match with the rest of the recipe's in QuAcc.
- The summary dictionaries will need to be written so that the data can be passed back to the database.
- Also, the requirements will be added in the NewtonNet's recipe, however, it is a general problem where ideally the calculator can be changed. If a decision is made later on to move it out of NewtonNet then things will (hopefully) be still usable a general path optimization process and can be added along with the static, relax, freq, and thermo jobs. But, right now, NewtonNet seems like an OK place for this to be put in.

 <<

### Checklist

- [ ] I have read the ["Guidelines" section](https://quantum-accelerators.github.io/quacc/dev/contributing.html#guidelines) of the contributing guide. Don't lie! 😉
- [ ] My PR is on a custom branch and is _not_ named `main`.
- [ ] I have added relevant, comprehensive [unit tests](https://quantum-accelerators.github.io/quacc/dev/contributing.html#unit-tests).

### Notes

- Your PR will likely not be merged without proper and thorough tests.
- If you are an external contributor, you will see a comment from [@buildbot-princeton](https://github.com/buildbot-princeton). This is solely for the maintainers.
- When your code is ready for review, ping one of the [active maintainers](https://quantum-accelerators.github.io/quacc/about/contributors.html#active-maintainers).


buildbot-princeton: Can one of the admins verify this patch?
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators) Report
Attention: Patch coverage is `76.54321%` with `19 lines` in your changes are missing coverage. Please review.
> Project coverage is 98.50%. Comparing base [(`ee85da2`)](https://app.codecov.io/gh/Quantum-Accelerators/quacc/commit/ee85da235d962bd8286784b7fc650967280129d9?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators) to head [(`3759dbf`)](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?dropdown=coverage&src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators).
> Report is 5 commits behind head on main.

| [Files](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators) | Patch % | Lines |
|---|---|---|
| [src/quacc/recipes/newtonnet/ts.py](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?src=pr&el=tree&filepath=src%2Fquacc%2Frecipes%2Fnewtonnet%2Fts.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators#diff-c3JjL3F1YWNjL3JlY2lwZXMvbmV3dG9ubmV0L3RzLnB5) | 76.54% | [19 Missing :warning: ](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #2176      +/-   ##
==========================================
- Coverage   99.04%   98.50%   -0.54%     
==========================================
  Files          81       81              
  Lines        3335     3416      +81     
==========================================
+ Hits         3303     3365      +62     
- Misses         32       51      +19     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/Quantum-Accelerators/quacc/pull/2176?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Quantum-Accelerators).

kumaranu: Hi @Andrew-S-Rosen
 
I have added a function in runners/ase.py with its test.

My test is running alright locally but it is complaining that the required packages are not installed such as "sella", "NewtonNet", "geodesic" etc.
https://github.com/Quantum-Accelerators/quacc/actions/runs/9449377214/job/26025492513?pr=2176

Please help resolve this issue. Thanks.

Tagging Sam here as I was discussing this issue with him recently.
@samblau 
kumaranu: > Hi @Andrew-S-Rosen
> 
> I have added a function in runners/ase.py with its test.
> 
> My test is running alright locally but it is complaining that the required packages are not installed such as "sella", "NewtonNet", "geodesic" etc. https://github.com/Quantum-Accelerators/quacc/actions/runs/9449377214/job/26025492513?pr=2176
> 
> Maybe the reason for this problem is that requirements-newtonnet.txt is not called for the packages that are present inside runners/ase.py. But I do not know which file will used to get that list of required packages.
> 
> Please help resolve this issue. Thanks.
> 
> Tagging Sam here as I was discussing this issue with him recently. @samblau

I added the following three lines to requirements.txt inside the tests directory to solve that problem:
newtonnet==1.1.1
geodesic-interpolate @ git+https://github.com/virtualzx-nad/geodesic-interpolate.git
sella==2.3.4
kumaranu: > > Hi @Andrew-S-Rosen
> > I have added a function in runners/ase.py with its test.
> > My test is running alright locally but it is complaining that the required packages are not installed such as "sella", "NewtonNet", "geodesic" etc. https://github.com/Quantum-Accelerators/quacc/actions/runs/9449377214/job/26025492513?pr=2176
> > Maybe the reason for this problem is that requirements-newtonnet.txt is not called for the packages that are present inside runners/ase.py. But I do not know which file will used to get that list of required packages.
> > Please help resolve this issue. Thanks.
> > Tagging Sam here as I was discussing this issue with him recently. @samblau
> 
> I added the following three lines to requirements.txt inside the tests directory to solve that problem: newtonnet==1.1.1 geodesic-interpolate @ git+https://github.com/virtualzx-nad/geodesic-interpolate.git sella==2.3.4

However, the original problem is solved, I saw two new errors, tblite and mlp tests both of which I did not touch.

mlp test were just complaining that they could not find mace. I added the following lines in the requirements.txt to get rid of that error.
mace-torch==0.3.4
torch-dftd==0.4.0

However, the error in tblite seems much different as the tests ran fine but the numbers don't match. I have left that as it is:
https://github.com/Quantum-Accelerators/quacc/actions/runs/9456768276/job/26049315963?pr=2176
Andrew-S-Rosen: @kumaranu: The `tests/requirements.txt` file cannot have the various packages you have added to it. It is meant to include a minimum set of packages corresponding to the `tests-core` CI run. If your tests require things like mace, they need to be running with the `tests-mlps` run. The reason the TBLite issues are happening is because TBLite is not compatible with torch (https://github.com/tblite/tblite/issues/110).

I do not quite understand what your goal is here with `tests/requirements.txt`. You have added NewtonNet tests, which should only run if NewtonNet is installed.
kumaranu: > @kumaranu: The `tests/requirements.txt` file cannot have the various packages you have added to it. It is meant to include a minimum set of packages corresponding to the `tests-core` CI run. If your tests require things like mace, they need to be running with the `tests-mlps` run. The reason the TBLite issues are happening is because TBLite is not compatible with torch ([tblite/tblite#110](https://github.com/tblite/tblite/issues/110)).
> 
> I do not quite understand what your goal is here with `tests/requirements.txt`. You have added NewtonNet tests, which should only run if NewtonNet is installed.

I have removed NewtonNet and Sella from the test_ase.py requirements inside runners directory. It is using BFGS and EMT now.
@Andrew-S-Rosen @samblau 
Andrew-S-Rosen: @kumaranu: Thanks. I will review the PR once you undo all the unrelated changes you've made (there are still many) and remove commented out code blocks.
kumaranu: Hi @Andrew-S-Rosen 

I am not sure why one of my tests is failing on Windows all of a sudden. Please let me know if you see anything that is obvious to you. Thanks.
Andrew-S-Rosen: It looks like a tolerance issue and should be not a problem. You can change the tolerance slightly (just enough for the tests to pass). I am not sure the source either. 
kumaranu: Hi @Andrew-S-Rosen 

I am struggling with the path of traj file for neb. The error is only present on GitHub but not when I run the tests locally.
Please let me know if you see any problems over here.
https://github.com/Quantum-Accelerators/quacc/actions/runs/10085853797/job/27887414542?pr=2176
I have not gone through all of your comments so I apologize for asking if you already mentioned it somewhere else. Thanks.
Andrew-S-Rosen: > Hi @Andrew-S-Rosen
> 
> I am struggling with the path of traj file for neb. The error is only present on GitHub but not when I run the tests locally. Please let me know if you see any problems over here. https://github.com/Quantum-Accelerators/quacc/actions/runs/10085853797/job/27887414542?pr=2176 I have not gone through all of your comments so I apologize for asking if you already mentioned it somewhere else. Thanks.

Please see this comment: https://github.com/Quantum-Accelerators/quacc/pull/2176#discussion_r1690015953.
kumaranu: > > Hi @Andrew-S-Rosen
> > I am struggling with the path of traj file for neb. The error is only present on GitHub but not when I run the tests locally. Please let me know if you see any problems over here. https://github.com/Quantum-Accelerators/quacc/actions/runs/10085853797/job/27887414542?pr=2176 I have not gone through all of your comments so I apologize for asking if you already mentioned it somewhere else. Thanks.
> 
> @kumaranu: Please see this comment: [#2176 (comment)](https://github.com/Quantum-Accelerators/quacc/pull/2176#discussion_r1690015953).

Done this. All checks are passing now. Thanks.
kumaranu: Hi @Andrew-S-Rosen,

Thanks a lot for providing the comments above. I have finished going through them all.

As you mention two things where you might need to pitch in:
1. run_neb method in the Runner class
2. summarize_opt_run with sumarize_neb_run.

I have left those two unresolved.

Tagging @samblau here just for the update.
Andrew-S-Rosen: Thanks, @kumaranu! I recognize that this has likely taken longer than you anticipated and involved more iterations than desired. That said, I feel that the end result here is much more robust and will ensure that your work can continue to be well-maintained into the future.

I am starting my new position Aug. 1st so can't wrap this up immediately, but I will keep it on my to-do list for the near future (middle of August most likely). As you suggested, the two things that need to happen before the final merge are: 1) Ensuring that the `run_neb` function becomes a method of the `Runner` class rather than a standalone function; 2) Ensuring that the `summarize_neb_run` function become a method of the `Summarize` class rather than a standalone function. Both scenarios here are meant to minimize code duplication so that changes to one part of the code are properly reflected in your NEB runner and schema. I will ping you when that's done so you can look it over before we wrap this up.
samblau: hey @kumaranu , any progress on this? Thanks!
kumaranu: > hey @kumaranu , any progress on this? Thanks!

Hi @samblau ,
Things have been quite busy lately. Sorry, I couldn't get back to it earlier.

I reduced the timing of the tests by adding max_iter to them and also got rid of writing to the current directory problem. Let's see if it works out OK or not. Thanks!
Andrew-S-Rosen: Alright friends, @kumaranu and @samblau. It is time for the merge! Is it perfect? No. Is it good enough? We'll find out! 

I did a _lot_ of refactoring here, but things seem to be okay. Let's move.

Patch:
diff --git a/docs/user/recipes/recipes_list.md b/docs/user/recipes/recipes_list.md
index fe58ce0e..28b25eb2 100644
--- a/docs/user/recipes/recipes_list.md
+++ b/docs/user/recipes/recipes_list.md
@@ -107,14 +107,16 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp
 
 <center>
 
-| Name                | Decorator       | Documentation                                | Req'd Extras   |
-| ------------------- | --------------- | -------------------------------------------- | -------------- |
-| NewtonNet Static    | `#!Python @job` | [quacc.recipes.newtonnet.core.static_job][]  |                |
-| NewtonNet Relax     | `#!Python @job` | [quacc.recipes.newtonnet.core.relax_job][]   |                |
-| NewtonNet Frequency | `#!Python @job` | [quacc.recipes.newtonnet.core.freq_job][]    |                |
-| NewtonNet TS        | `#!Python @job` | [quacc.recipes.newtonnet.ts.ts_job][]        | `quacc[sella]` |
-| NewtonNet IRC       | `#!Python @job` | [quacc.recipes.newtonnet.ts.irc_job][]       | `quacc[sella]` |
-| NewtonNet Quasi IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.quasi_irc_job][] | `quacc[sella]` |
+| Name                | Decorator       | Documentation                                 | Req'd Extras   |
+|---------------------| --------------- |-----------------------------------------------|                |
+| NewtonNet Static    | `#!Python @job` | [quacc.recipes.newtonnet.core.static_job][]   |                |
+| NewtonNet Relax     | `#!Python @job` | [quacc.recipes.newtonnet.core.relax_job][]    |                |
+| NewtonNet Frequency | `#!Python @job` | [quacc.recipes.newtonnet.core.freq_job][]     |                |
+| NewtonNet TS        | `#!Python @job` | [quacc.recipes.newtonnet.ts.ts_job][]         | `quacc[sella]` |
+| NewtonNet IRC       | `#!Python @job` | [quacc.recipes.newtonnet.ts.irc_job][]        | `quacc[sella]` |
+| NewtonNet Quasi IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.quasi_irc_job][]  | `quacc[sella]` |
+| NewtonNet neb       | `#!Python @job` | [quacc.recipes.newtonnet.ts.neb_job][]        |                |
+| NewtonNet geodesic  | `#!Python @job` | [quacc.recipes.newtonnet.ts.geodesic_job][]   |                |
 
 </center>
 
diff --git a/src/quacc/atoms/ts.py b/src/quacc/atoms/ts.py
new file mode 100644
index 00000000..cd7329eb
--- /dev/null
+++ b/src/quacc/atoms/ts.py
@@ -0,0 +1,115 @@
+from __future__ import annotations
+
+import logging
+from importlib.util import find_spec
+from typing import TYPE_CHECKING
+
+from ase.atoms import Atoms
+from monty.dev import requires
+
+from quacc.atoms.core import copy_atoms
+
+LOGGER = logging.getLogger(__name__)
+
+has_geodesic_interpolate = bool(find_spec("geodesic_interpolate"))
+
+if has_geodesic_interpolate:
+    from geodesic_interpolate.geodesic import Geodesic
+    from geodesic_interpolate.interpolation import redistribute
+
+if TYPE_CHECKING:
+    from typing import Literal
+
+
+@requires(
+    has_geodesic_interpolate,
+    "geodesic-interpolate must be installed. Refer to the quacc documentation.",
+)
+def geodesic_interpolate_wrapper(
+    reactant: Atoms,
+    product: Atoms,
+    n_images: int = 10,
+    perform_sweep: bool | Literal["auto"] = "auto",
+    redistribute_tol: float = 1e-2,
+    smoother_tol: float = 2e-3,
+    max_iterations: int = 15,
+    max_micro_iterations: int = 20,
+    morse_scaling: float = 1.7,
+    geometry_friction: float = 1e-2,
+    distance_cutoff: float = 3.0,
+    sweep_cutoff_size: int = 35,
+) -> list[Atoms]:
+    """
+    Interpolates between two geometries and optimizes the path with the geodesic method.
+
+    Parameters
+    ----------
+    reactant
+        The ASE Atoms object representing the initial geometry.
+    product
+        The ASE Atoms object representing the final geometry.
+    n_images
+        Number of images for interpolation. Default is 10.
+    perform_sweep
+        Whether to sweep across the path optimizing one image at a time.
+        Default is to perform sweeping updates if there are more than 35 atoms.
+    redistribute_tol
+        the value passed to the tol keyword argument of
+         geodesic_interpolate.interpolation.redistribute. Default is 1e-2.
+    smoother_tol
+        the value passed to the tol keyword argument of geodesic_smoother.smooth
+        or geodesic_smoother.sweep. Default is 2e-3.
+    max_iterations
+        Maximum number of minimization iterations. Default is 15.
+    max_micro_iterations
+        Maximum number of micro iterations for the sweeping algorithm. Default is 20.
+    morse_scaling
+        Exponential parameter for the Morse potential. Default is 1.7.
+    geometry_friction
+        Size of friction term used to prevent very large changes in geometry. Default is 1e-2.
+    distance_cutoff
+        Cut-off value for the distance between a pair of atoms to be included in the coordinate system. Default is 3.0.
+    sweep_cutoff_size
+        Cut off system size that above which sweep function will be called instead of smooth
+        in Geodesic.
+
+    Returns
+    -------
+    list[Atoms]
+        A list of ASE Atoms objects representing the smoothed path between the reactant and product geometries.
+    """
+    reactant = copy_atoms(reactant)
+    product = copy_atoms(product)
+
+    # Read the initial geometries.
+    chemical_symbols = reactant.get_chemical_symbols()
+
+    # First redistribute number of images. Perform interpolation if too few and subsampling if too many images are given
+    raw_interpolated_positions = redistribute(
+        chemical_symbols,
+        [reactant.positions, product.positions],
+        n_images,
+        tol=redistribute_tol,
+    )
+
+    # Perform smoothing by minimizing distance in Cartesian coordinates with redundant internal metric
+    # to find the appropriate geodesic curve on the hyperspace.
+    geodesic_smoother = Geodesic(
+        chemical_symbols,
+        raw_interpolated_positions,
+        morse_scaling,
+        threshold=distance_cutoff,
+        friction=geometry_friction,
+    )
+    if perform_sweep == "auto":
+        perform_sweep = len(chemical_symbols) > sweep_cutoff_size
+    if perform_sweep:
+        geodesic_smoother.sweep(
+            tol=smoother_tol, max_iter=max_iterations, micro_iter=max_micro_iterations
+        )
+    else:
+        geodesic_smoother.smooth(tol=smoother_tol, max_iter=max_iterations)
+    return [
+        Atoms(symbols=chemical_symbols, positions=geom)
+        for geom in geodesic_smoother.path
+    ]
diff --git a/src/quacc/recipes/newtonnet/core.py b/src/quacc/recipes/newtonnet/core.py
index ed632313..929dcecf 100644
--- a/src/quacc/recipes/newtonnet/core.py
+++ b/src/quacc/recipes/newtonnet/core.py
@@ -204,7 +204,7 @@ def freq_job(
     )
 
 
-def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]:
+def _add_stdev_and_hess(summary: dict[str, Any], **calc_kwargs) -> dict[str, Any]:
     """
     Calculate and add standard deviation values and Hessians to the summary.
 
@@ -218,6 +218,10 @@ def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]:
     ----------
     summary
         A dictionary containing information about the molecular trajectory.
+    **calc_kwargs
+        Custom kwargs for the NewtonNet calculator. Set a value to
+        `quacc.Remove` to remove a pre-existing key entirely. For a list of available
+        keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator.
 
     Returns
     -------
@@ -226,11 +230,13 @@ def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]:
         Hessian values.
     """
     settings = get_settings()
+    calc_defaults = {
+        "model_path": settings.NEWTONNET_MODEL_PATH,
+        "settings_path": settings.NEWTONNET_CONFIG_PATH,
+    }
+    calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs)
     for i, atoms in enumerate(summary["trajectory"]):
-        calc = NewtonNet(
-            model_path=settings.NEWTONNET_MODEL_PATH,
-            settings_path=settings.NEWTONNET_CONFIG_PATH,
-        )
+        calc = NewtonNet(**calc_flags)
         results = Runner(atoms, calc).run_calc().calc.results
         summary["trajectory_results"][i]["hessian"] = results["hessian"]
         summary["trajectory_results"][i]["energy_std"] = results["energy_disagreement"]
diff --git a/src/quacc/recipes/newtonnet/ts.py b/src/quacc/recipes/newtonnet/ts.py
index a99a5142..b379b4ea 100644
--- a/src/quacc/recipes/newtonnet/ts.py
+++ b/src/quacc/recipes/newtonnet/ts.py
@@ -5,6 +5,8 @@ from __future__ import annotations
 from importlib.util import find_spec
 from typing import TYPE_CHECKING
 
+import numpy as np
+from ase.mep import NEB
 from monty.dev import requires
 
 from quacc import change_settings, get_settings, job, strip_decorator
@@ -13,6 +15,7 @@ from quacc.runners.ase import Runner
 from quacc.schemas.ase import Summarize
 from quacc.utils.dicts import recursive_dict_merge
 
+has_geodesic_interpolate = bool(find_spec("geodesic_interpolate"))
 has_sella = bool(find_spec("sella"))
 has_newtonnet = bool(find_spec("newtonnet"))
 
@@ -20,7 +23,8 @@ if has_sella:
     from sella import IRC, Sella
 if has_newtonnet:
     from newtonnet.utils.ase_interface import MLAseCalculator as NewtonNet
-
+if has_geodesic_interpolate:
+    from quacc.atoms.ts import geodesic_interpolate_wrapper
 
 if TYPE_CHECKING:
     from typing import Any, Literal
@@ -29,6 +33,8 @@ if TYPE_CHECKING:
     from numpy.typing import NDArray
 
     from quacc.types import (
+        GeodesicSchema,
+        NebSchema,
         NewtonNetIRCSchema,
         NewtonNetQuasiIRCSchema,
         NewtonNetTSSchema,
@@ -85,7 +91,6 @@ def ts_job(
     calc_defaults = {
         "model_path": settings.NEWTONNET_MODEL_PATH,
         "settings_path": settings.NEWTONNET_CONFIG_PATH,
-        "hess_method": "autograd",
     }
     opt_defaults = {
         "optimizer": Sella,
@@ -105,7 +110,7 @@ def ts_job(
     # Run the TS optimization
     dyn = Runner(atoms, calc).run_opt(**opt_flags)
     opt_ts_summary = _add_stdev_and_hess(
-        Summarize(additional_fields={"name": "NewtonNet TS"}).opt(dyn)
+        Summarize(additional_fields={"name": "NewtonNet TS"}).opt(dyn), **calc_flags
     )
 
     # Run a frequency calculation
@@ -267,7 +272,193 @@ def quasi_irc_job(
     return relax_summary
 
 
-def _get_hessian(atoms: Atoms) -> NDArray:
+@job
+@requires(
+    has_newtonnet, "NewtonNet must be installed. Refer to the quacc documentation."
+)
+def neb_job(
+    reactant_atoms: Atoms,
+    product_atoms: Atoms,
+    interpolation_method: Literal["linear", "idpp", "geodesic"] = "linear",
+    relax_job_kwargs: dict[str, Any] | None = None,
+    interpolate_kwargs: dict[str, Any] | None = None,
+    neb_kwargs: dict[str, Any] | None = None,
+    **calc_kwargs,
+) -> NebSchema:
+    """
+    Perform a nudged elastic band (NEB) calculation to find the minimum energy path (MEP) between the given reactant and product structures.
+
+    Parameters
+    ----------
+    reactant_atoms
+        The Atoms object representing the reactant structure.
+    product_atoms
+        The Atoms object representing the product structure.
+    interpolation_method
+        The method to initialize the NEB optimization. There are three choices here, "linear", "idpp" and "geodesic".
+        Defaults to linear.
+    relax_job_kwargs
+        Keyword arguments to use for the [quacc.recipes.newtonnet.core.relax_job][] function.
+    interpolate_kwargs
+        Keyword arguments for the interpolate functions (geodesic, linear or idpp).
+    neb_kwargs
+        Keyword arguments for the NEB calculation.
+    **calc_kwargs
+        Dictionary of custom kwargs for the NewtonNet calculator. Set a value to
+        `quacc.Remove` to remove a pre-existing key entirely. For a list of available
+        keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator.
+
+    Returns
+    -------
+    NebSchema
+        A dictionary containing the NEB results
+    """
+    relax_job_kwargs = relax_job_kwargs or {}
+    settings = get_settings()
+
+    calc_defaults = {
+        "model_path": settings.NEWTONNET_MODEL_PATH,
+        "settings_path": settings.NEWTONNET_CONFIG_PATH,
+        "hess_method": None,
+    }
+
+    interpolate_defaults = {"n_images": 20}
+
+    neb_defaults = {"method": "aseneb", "precon": None}
+    calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs)
+    interpolate_flags = recursive_dict_merge(interpolate_defaults, interpolate_kwargs)
+    neb_flags = recursive_dict_merge(neb_defaults, neb_kwargs)
+
+    # Define calculator
+    calc = NewtonNet(**calc_flags)
+
+    # Run relax job
+    relax_summary_r = strip_decorator(relax_job)(reactant_atoms, **relax_job_kwargs)
+    relax_summary_p = strip_decorator(relax_job)(product_atoms, **relax_job_kwargs)
+
+    if interpolation_method == "geodesic":
+        images = geodesic_interpolate_wrapper(
+            relax_summary_r["atoms"], relax_summary_p["atoms"], **interpolate_flags
+        )
+    else:
+        images = [reactant_atoms]
+        images += [
+            reactant_atoms.copy() for _ in range(interpolate_flags["n_images"] - 2)
+        ]
+        images += [product_atoms]
+        neb = NEB(images)
+
+        # Interpolate linearly the positions of the middle images:
+        neb.interpolate(method=interpolation_method)
+        images = neb.images
+
+    max_steps = neb_flags.pop("max_steps", None)
+    dyn = Runner(images, calc).run_neb(max_steps=max_steps, neb_kwargs=neb_flags)
+
+    return {
+        "relax_reactant": relax_summary_r,
+        "relax_product": relax_summary_p,
+        "initial_images": images,
+        "neb_results": Summarize(
+            additional_fields={
+                "neb_flags": neb_flags,
+                "calc_flags": calc_flags,
+                "interpolate_flags": interpolate_flags,
+            }
+        ).neb(dyn, len(images)),
+    }
+
+
+@job
+@requires(
+    has_newtonnet, "NewtonNet must be installed. Refer to the quacc documentation."
+)
+@requires(
+    has_geodesic_interpolate,
+    "geodesic-interpolate must be installed. Refer to the quacc documentation.",
+)
+def geodesic_job(
+    reactant_atoms: Atoms,
+    product_atoms: Atoms,
+    relax_job_kwargs: dict[str, Any] | None = None,
+    geodesic_interpolate_kwargs: dict[str, Any] | None = None,
+    **calc_kwargs,
+) -> GeodesicSchema:
+    """
+    Perform a quasi-IRC job using the given reactant and product atoms objects.
+
+    Parameters
+    ----------
+    reactant_atoms
+        The Atoms object representing the reactant structure.
+    product_atoms
+        The Atoms object representing the product structure.
+    relax_job_kwargs
+        Keyword arguments to use for [quacc.recipes.newtonnet.core.relax_job][] function.
+    geodesic_interpolate_kwargs
+        Keyword arguments for geodesic_interpolate, by default None.
+    **calc_kwargs
+        Dictionary of custom kwargs for the NewtonNet calculator. Set a value to
+        `quacc.Remove` to remove a pre-existing key entirely. For a list of available
+        keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator.
+
+    Returns
+    -------
+    GeodesicSchema
+        A dictionary containing the following keys:
+            - 'relax_reactant': Summary of the relaxed reactant structure.
+            - 'relax_product': Summary of the relaxed product structure.
+            - 'initial_images': The interpolated images between reactant and product.
+            - 'ts_atoms': ASE atoms object for the highest energy structure for the geodesic path
+    """
+    relax_job_kwargs = relax_job_kwargs or {}
+    settings = get_settings()
+
+    calc_defaults = {
+        "model_path": settings.NEWTONNET_MODEL_PATH,
+        "settings_path": settings.NEWTONNET_CONFIG_PATH,
+        "hess_method": None,
+    }
+
+    geodesic_defaults = {"n_images": 20}
+
+    calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs)
+    calc_flags["hess_method"] = None
+    geodesic_interpolate_flags = recursive_dict_merge(
+        geodesic_defaults, geodesic_interpolate_kwargs
+    )
+
+    # Define calculator
+    reactant_atoms.calc = NewtonNet(**calc_flags)
+    product_atoms.calc = NewtonNet(**calc_flags)
+
+    # Run IRC
+    relax_summary_r = strip_decorator(relax_job)(reactant_atoms, **relax_job_kwargs)
+    relax_summary_p = strip_decorator(relax_job)(product_atoms, **relax_job_kwargs)
+
+    images = geodesic_interpolate_wrapper(
+        relax_summary_r["atoms"].copy(),
+        relax_summary_p["atoms"].copy(),
+        **geodesic_interpolate_flags,
+    )
+
+    potential_energies = []
+    for image in images:
+        image.calc = NewtonNet(**calc_flags)
+        potential_energies.append(image.get_potential_energy())
+
+    ts_index = np.argmax(potential_energies)
+    ts_atoms = images[ts_index]
+
+    return {
+        "relax_reactant": relax_summary_r,
+        "relax_product": relax_summary_p,
+        "initial_images": images,
+        "ts_atoms": ts_atoms,
+    }
+
+
+def _get_hessian(atoms: Atoms, **calc_kwargs) -> NDArray:
     """
     Calculate and retrieve the Hessian matrix for the given molecular configuration.
 
@@ -280,6 +471,10 @@ def _get_hessian(atoms: Atoms) -> NDArray:
     ----------
     atoms
         The ASE Atoms object representing the molecular configuration.
+    **calc_kwargs
+        Dictionary of custom kwargs for the NewtonNet calculator. Set a value to
+        `quacc.Remove` to remove a pre-existing key entirely. For a list of available
+        keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator.
 
     Returns
     -------
@@ -287,10 +482,12 @@ def _get_hessian(atoms: Atoms) -> NDArray:
         The calculated Hessian matrix, reshaped into a 2D array.
     """
     settings = get_settings()
-    ml_calculator = NewtonNet(
-        model_path=settings.NEWTONNET_MODEL_PATH,
-        settings_path=settings.NEWTONNET_CONFIG_PATH,
-    )
-    ml_calculator.calculate(atoms)
-
-    return ml_calculator.results["hessian"].reshape((-1, 3 * len(atoms)))
+    calc_defaults = {
+        "model_path": settings.NEWTONNET_MODEL_PATH,
+        "settings_path": settings.NEWTONNET_CONFIG_PATH,
+        "hess_method": "autograd",
+    }
+    calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs)
+    calc = NewtonNet(**calc_flags)
+    calc.calculate(atoms)
+    return calc.results["hessian"].reshape((-1, 3 * len(atoms)))
diff --git a/src/quacc/runners/ase.py b/src/quacc/runners/ase.py
index 4157c870..74be4a1c 100644
--- a/src/quacc/runners/ase.py
+++ b/src/quacc/runners/ase.py
@@ -2,10 +2,13 @@
 
 from __future__ import annotations
 
+from collections.abc import Callable
+from copy import deepcopy
 from importlib.util import find_spec
 from logging import getLogger
+from pathlib import Path
 from shutil import copy, copytree
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
 from ase.calculators import calculator
@@ -17,21 +20,22 @@ from ase.md.velocitydistribution import (
     Stationary,
     ZeroRotation,
 )
-from ase.optimize import BFGS
+from ase.mep import NEB
+from ase.mep.neb import NEBOptimizer
+from ase.optimize import BFGS, BFGSLineSearch
 from ase.optimize.sciopt import SciPyOptimizer
 from ase.vibrations import Vibrations
 from monty.dev import requires
 from monty.os.path import zpath
 
-from quacc.atoms.core import copy_atoms
 from quacc.runners._base import BaseRunner
-from quacc.runners.prep import terminate
+from quacc.runners.prep import calc_cleanup, calc_setup, terminate
 from quacc.utils.dicts import recursive_dict_merge
 
 LOGGER = getLogger(__name__)
 
 has_sella = bool(find_spec("sella"))
-
+has_geodesic_interpolate = bool(find_spec("geodesic_interpolate"))
 
 if TYPE_CHECKING:
     from collections.abc import Callable
@@ -40,7 +44,7 @@ if TYPE_CHECKING:
 
     from ase.atoms import Atoms
     from ase.calculators.calculator import Calculator
-    from ase.optimize.optimize import Dynamics
+    from ase.optimize.optimize import Dynamics, Optimizer
 
     from quacc.types import (
         Filenames,
@@ -58,7 +62,7 @@ class Runner(BaseRunner):
 
     def __init__(
         self,
-        atoms: Atoms,
+        atoms: Atoms | list[Atoms],
         calculator: Calculator,
         copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None,
     ) -> None:
@@ -68,7 +72,7 @@ class Runner(BaseRunner):
         Parameters
         ----------
         atoms
-            The Atoms object to run calculations on.
+            The Atoms object to run calculations on. A list[Atoms] is used for NEB.
         calculator
             The instantiated ASE calculator object to attach to the Atoms object.
         copy_files
@@ -78,10 +82,15 @@ class Runner(BaseRunner):
         -------
         None
         """
-        self.atoms = copy_atoms(atoms)
-        self.atoms.calc = calculator
         self.copy_files = copy_files
-        self.setup()
+        if isinstance(atoms, list):
+            self.atoms = [image.copy() for image in atoms]
+            for image in self.atoms:
+                image.calc = deepcopy(calculator)
+        else:
+            self.atoms = atoms.copy()
+            self.atoms.calc = calculator
+            self.setup()
 
     def run_calc(
         self, geom_file: str | None = None, properties: list[str] | None = None
@@ -347,6 +356,97 @@ class Runner(BaseRunner):
             optimizer_kwargs=dynamics_kwargs,
         )
 
+    def run_neb(
+        self,
+        relax_cell: bool = False,
+        fmax: float = 0.01,
+        max_steps: int | None = 1000,
+        optimizer: Optimizer = NEBOptimizer,
+        optimizer_kwargs: dict[str, Any] | None = None,
+        neb_kwargs: dict[str, Any] | None = None,
+        run_kwargs: dict[str, Any] | None = None,
+    ) -> Dynamics:
+        """
+        Run an NEB calculation.
+
+        Parameters
+        ----------
+        relax_cell
+            Whether to relax the unit cell shape and volume.
+        fmax
+            Tolerance for the force convergence (in eV/A).
+        max_steps
+            Maximum number of steps to take.
+        optimizer
+            Optimizer class to use.
+        optimizer_kwargs
+            Dictionary of kwargs for the optimizer. Takes all valid kwargs for ASE
+            Optimizer classes.
+        neb_kwargs
+            Dictionary of kwargs for the NEB class.
+        run_kwargs
+            Dictionary of kwargs for the `run()` method of the optimizer.
+
+        Returns
+        -------
+        Dynamics
+            The ASE Dynamics object following an NEB calculation.
+        """
+        images = self.atoms
+        run_kwargs = run_kwargs or {}
+        neb_kwargs = neb_kwargs or {}
+        traj_filename = "opt.traj"
+
+        # Create a parent temporary directory for the NEB run
+        neb_tmpdir, neb_results_dir = calc_setup(None)
+
+        # Adjust optimizer_kwargs to use the parent directory
+        optimizer_kwargs = recursive_dict_merge(
+            {
+                "logfile": str(neb_tmpdir / "opt.log"),
+                "restart": str(neb_tmpdir / "opt.json"),
+            },
+            optimizer_kwargs,
+        )
+
+        if "trajectory" in optimizer_kwargs:
+            msg = "Quacc does not support setting the `trajectory` kwarg."
+            raise ValueError(msg)
+
+        if optimizer == BFGSLineSearch:
+            raise ValueError("BFGSLineSearch is not allowed as optimizer with NEB.")
+
+        # Copy atoms so we don't modify it in-place
+        neb = NEB(images, **neb_kwargs)
+
+        # Perform staging operations
+        for i, image in enumerate(images):
+            image_tmpdir = neb_tmpdir / f"image_{i}"
+            image_tmpdir.mkdir()
+            image.calc.directory = image_tmpdir
+
+        # Define the Trajectory object
+        traj_file = neb_tmpdir / traj_filename
+        traj = Trajectory(traj_file, "w", atoms=neb)
+
+        # Set volume relaxation constraints, if relevant
+        if relax_cell:
+            for i in range(len(images)):
+                if images[i].pbc.any():
+                    images[i] = FrechetCellFilter(images[i])
+
+        dyn = optimizer(neb, **optimizer_kwargs)
+        dyn.attach(traj.write)
+        dyn.run(fmax, max_steps)
+        traj.close()
+        dyn.logfile.close()
+
+        calc_cleanup(None, neb_tmpdir, neb_results_dir)
+        traj.filename = zpath(str(neb_results_dir / traj_filename))
+        dyn.trajectory = traj
+
+        return dyn
+
     def _copy_intermediate_files(
         self, step_number: int, files_to_ignore: list[Path] | None = None
     ) -> None:
diff --git a/src/quacc/schemas/ase.py b/src/quacc/schemas/ase.py
index 084fdad8..68103ac9 100644
--- a/src/quacc/schemas/ase.py
+++ b/src/quacc/schemas/ase.py
@@ -23,6 +23,7 @@ if TYPE_CHECKING:
     from typing import Any, Literal
 
     from ase.atoms import Atoms
+    from ase.io.trajectory import TrajectoryWriter
     from ase.md.md import MolecularDynamics
     from ase.optimize.optimize import Optimizer
     from ase.vibrations import Vibrations
@@ -293,6 +294,89 @@ class Summarize:
             store=store,
         )
 
+    def neb(
+        self,
+        dyn: Optimizer,
+        n_images: int,
+        n_iter_return: int = -1,
+        trajectory: TrajectoryWriter | list[Atoms] | None = None,
+        store: Store | None | DefaultSetting = QuaccDefault,
+    ) -> OptSchema:
+        """
+        Summarize the NEB run results and store them in a database-friendly format.
+
+        Parameters
+        ----------
+        dyn
+            ASE Optimizer object used for the NEB run.
+        n_images
+            Number of images in the NEB run.
+        n_iter_return
+            Number of iterations to return. If -1, all iterations are returned.
+        trajectory
+            Trajectory of the NEB run, either as a Trajectory object or a list of Atoms objects.
+        store
+            Maggma Store object to store the results in. Defaults to `QuaccSettings.STORE`.
+
+        Returns
+        -------
+        OptSchema
+            A dictionary containing the summarized NEB run results.
+        """
+        store = self._settings.STORE if store == QuaccDefault else store
+
+        # Get trajectory
+        if trajectory:
+            atoms_trajectory = trajectory
+        else:
+            atoms_trajectory = read(dyn.trajectory.filename, index=":")  # type: ignore[union-attr]
+
+        if n_iter_return == -1:
+            atoms_trajectory = atoms_trajectory[-(n_images):]
+        else:
+            atoms_trajectory = _get_nth_iteration(
+                atoms_trajectory,
+                int(len(atoms_trajectory) / n_images),
+                n_images,
+                n_iter_return,
+            )
+        trajectory_results = [atoms.calc.results for atoms in atoms_trajectory]
+        ts_index = (
+            np.argmax(
+                [
+                    result["energy"]
+                    for result in trajectory_results[-(n_images - 1) : -1]
+                ]
+            )
+            + 1
+        )
+        ts_atoms = atoms_trajectory[ts_index]
+        base_task_doc = atoms_to_metadata(
+            atoms_trajectory[0], charge_and_multiplicity=self.charge_and_multiplicity
+        )
+
+        # Clean up the opt parameters
+        parameters_opt = dyn.todict()
+        parameters_opt.pop("logfile", None)
+        parameters_opt.pop("restart", None)
+
+        opt_fields = {
+            "parameters_opt": parameters_opt,
+            "trajectory": atoms_trajectory,
+            "trajectory_results": trajectory_results,
+            "ts_atoms": ts_atoms,
+        }
+
+        # Create a dictionary of the inputs/outputs
+        unsorted_task_doc = base_task_doc | opt_fields | self.additional_fields
+
+        return finalize_dict(
+            unsorted_task_doc,
+            directory=None,
+            gzip_file=self._settings.GZIP_FILES,
+            store=store,
+        )
+
 
 class VibSummarize:
     """
@@ -516,3 +600,37 @@ class VibSummarize:
             gzip_file=self._settings.GZIP_FILES,
             store=store,
         )
+
+
+def _get_nth_iteration(
+    neb_trajectory: list[Atoms], n_iter: int, n_images: int, interval: int
+) -> list[Atoms]:
+    """
+    Extract every nth iteration from the NEB trajectory.
+
+    Parameters
+    ----------
+    neb_trajectory
+        List of configurations (length: n_iter * n_images).
+    n_iter
+        Total number of iterations.
+    n_images
+        Number of images per iteration.
+    interval
+        Interval to get every nth iteration.
+
+    Returns
+    -------
+    list[Atoms]
+        List of configurations from every nth iteration.
+    """
+    result = []
+    start_idx, end_idx = 0, 0
+    for i in range(0, n_iter, interval):
+        start_idx = i * n_images
+        end_idx = start_idx + n_images
+
+        result.extend(neb_trajectory[start_idx:end_idx])
+    if end_idx < len(neb_trajectory) - 1:
+        result.extend(neb_trajectory[-(n_images):])
+    return result
diff --git a/src/quacc/types.py b/src/quacc/types.py
index f9db8555..56de12a0 100644
--- a/src/quacc/types.py
+++ b/src/quacc/types.py
@@ -685,6 +685,31 @@ if TYPE_CHECKING:
         irc_job: NewtonNetIRCSchema
         freq_job: VibThermoSchema | None
 
+    class NebSchema(TypedDict):
+        relax_reactant: OptSchema
+        relax_product: OptSchema
+        initial_images: list[Atoms]
+        neb_results: dict
+
+    class NebTsSchema(TypedDict):
+        relax_reactant: OptSchema
+        relax_product: OptSchema
+        initial_images: list[Atoms]
+        neb_results: dict
+        ts_results: OptSchema
+
+    class GeodesicSchema(TypedDict):
+        relax_reactant: OptSchema
+        relax_product: OptSchema
+        initial_images: list[Atoms]
+        ts_atoms: Atoms
+
+    class GeodesicTsSchema(TypedDict):
+        relax_reactant: OptSchema
+        relax_product: OptSchema
+        initial_images: list[Atoms]
+        ts_results: NewtonNetTSSchema
+
     # ----------- Recipe (Q-Chem) type hints -----------
 
     class QchemQuasiIRCSchema(OptSchema):


Test patch:
diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml
index a7d4433c..04cf7688 100644
--- a/.github/workflows/tests.yaml
+++ b/.github/workflows/tests.yaml
@@ -344,7 +344,7 @@ jobs:
           echo "COVERAGE_CORE=sysmon" >> $GITHUB_ENV
 
       - name: Run tests with pytest
-        run: pytest -k 'mlp or newtonnet' --durations=10 --cov=quacc --cov-report=xml
+        run: pytest -k 'mlp or newtonnet or geodesic' --durations=10 --cov=quacc --cov-report=xml
 
       - name: Upload code coverage report to Artifact
         uses: actions/upload-artifact@v4
diff --git a/tests/core/atoms/test_ts.py b/tests/core/atoms/test_ts.py
new file mode 100644
index 00000000..61221882
--- /dev/null
+++ b/tests/core/atoms/test_ts.py
@@ -0,0 +1,83 @@
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+
+@pytest.fixture(scope="module", autouse=True)
+def set_seed():
+    np.random.seed(42)  # noqa: NPY002
+
+
+from importlib.util import find_spec
+
+from ase.atoms import Atoms
+from ase.build import molecule
+
+from quacc.atoms.ts import geodesic_interpolate_wrapper
+
+has_geodesic_interpolate = bool(find_spec("geodesic_interpolate"))
+
+
+@pytest.fixture
+def setup_test_environment(tmp_path):
+    reactant = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.4835950817281542, -1.0145410211301968, -0.13209027203235943],
+            [0.8409564131524673, 0.018549610257914483, -0.07338809662321308],
+            [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891],
+            [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892],
+            [1.402180736662139, 0.944112416574632, -0.12179540364365492],
+            [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071],
+            [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239],
+        ],
+    )
+
+    product = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.348003553501624, 0.4819311116778978, 0.2752537177143993],
+            [0.2386618286631742, -0.3433222966734429, 0.37705518940917926],
+            [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955],
+            [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953],
+            [0.3801391040059668, -1.3793340533058087, 0.71035902765307],
+            [1.9296265384257907, 0.622088341468767, 1.0901733942191298],
+            [-1.090815880212625, 1.0965111343610956, -0.23791518420660265],
+        ],
+    )
+    return reactant, product
+
+
+@pytest.mark.skipif(
+    not has_geodesic_interpolate,
+    reason="geodesic_interpolate function is not available",
+)
+def test_geodesic_interpolate_wrapper(setup_test_environment):
+    reactant, product = setup_test_environment
+
+    # Execute the geodesic_interpolate_wrapper function
+    smoother_path = geodesic_interpolate_wrapper(
+        reactant,
+        product,
+        n_images=20,
+        redistribute_tol=5e-4,
+        smoother_tol=1e-4,
+        max_iterations=10,
+        max_micro_iterations=10,
+        morse_scaling=1.5,
+        geometry_friction=1e-2,
+        distance_cutoff=2.5,
+    )
+    assert smoother_path[1].positions[0][0] == pytest.approx(1.378384900, abs=1e-5)
+    assert smoother_path[5].positions[0][2] == pytest.approx(-0.512075394, abs=1e-5)
+
+
+@pytest.mark.skipif(
+    not has_geodesic_interpolate,
+    reason="geodesic_interpolate function is not available",
+)
+def test_geodesic_interpolate_wrapper_large_system(setup_test_environment):
+    # Test with large system to trigger sweeping updates
+    smoother_path = geodesic_interpolate_wrapper(molecule("C60"), molecule("C60"))
+    assert len(smoother_path) == 10
diff --git a/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py b/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py
index ebb0cb6e..00bd3a70 100644
--- a/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py
+++ b/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py
@@ -2,16 +2,29 @@ from __future__ import annotations
 
 import pytest
 
+
+@pytest.fixture(scope="module", autouse=True)
+def set_seed():
+    np.random.seed(42)  # noqa: NPY002
+
+
 pytest.importorskip("sella")
 pytest.importorskip("newtonnet")
 from pathlib import Path
 
 import numpy as np
+from ase import Atoms
 from ase.build import molecule
 
 from quacc import _internally_set_settings
 from quacc.recipes.newtonnet.core import freq_job, relax_job, static_job
-from quacc.recipes.newtonnet.ts import irc_job, quasi_irc_job, ts_job
+from quacc.recipes.newtonnet.ts import (
+    geodesic_job,
+    irc_job,
+    neb_job,
+    quasi_irc_job,
+    ts_job,
+)
 
 current_file_path = Path(__file__).parent
 
@@ -106,7 +119,7 @@ def test_ts_job_with_default_args(tmp_path, monkeypatch):
     assert "freq_job" in output
     assert output["results"]["energy"] == pytest.approx(-6.796914263061945)
     assert output["freq_job"]["results"]["imag_vib_freqs"][0] == pytest.approx(
-        -2426.7398321816004
+        -2426.73983218162, abs=1e-6
     )
 
 
@@ -282,3 +295,103 @@ def test_quasi_irc_job_with_custom_irc_swaps(tmp_path, monkeypatch):
     assert output["irc_job"]["results"]["energy"] == pytest.approx(-9.517354965639784)
     assert output["results"]["energy"] == pytest.approx(-9.517354965639784)
     assert output["freq_job"]["results"]["energy"] == pytest.approx(-9.517354965639784)
+
+
+@pytest.fixture
+def setup_test_environment():
+    reactant = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.4835950817281542, -1.0145410211301968, -0.13209027203235943],
+            [0.8409564131524673, 0.018549610257914483, -0.07338809662321308],
+            [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891],
+            [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892],
+            [1.402180736662139, 0.944112416574632, -0.12179540364365492],
+            [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071],
+            [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239],
+        ],
+    )
+
+    product = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.348003553501624, 0.4819311116778978, 0.2752537177143993],
+            [0.2386618286631742, -0.3433222966734429, 0.37705518940917926],
+            [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955],
+            [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953],
+            [0.3801391040059668, -1.3793340533058087, 0.71035902765307],
+            [1.9296265384257907, 0.622088341468767, 1.0901733942191298],
+            [-1.090815880212625, 1.0965111343610956, -0.23791518420660265],
+        ],
+    )
+
+    return reactant, product
+
+
+def test_neb_job_linear(setup_test_environment):
+    reactant, product = setup_test_environment
+
+    neb_summary = neb_job(
+        reactant, product, interpolation_method="linear", neb_kwargs={"max_steps": 10}
+    )
+
+    assert len(neb_summary["neb_results"]["trajectory_results"]) == 20
+    assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx(
+        0.8815, abs=1e-3
+    )
+    assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx(
+        1.117689, abs=1e-3
+    )
+
+    assert neb_summary["neb_results"]["trajectory_results"][1][
+        "energy"
+    ] == pytest.approx(-21.99570976499786, abs=1e-5)
+
+
+def test_neb_job_idpp(setup_test_environment):
+    reactant, product = setup_test_environment
+
+    neb_summary = neb_job(
+        reactant, product, interpolation_method="idpp", neb_kwargs={"max_steps": 10}
+    )
+
+    assert len(neb_summary["neb_results"]["trajectory_results"]) == 20
+    assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx(
+        0.8815, abs=1e-3
+    )
+    assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx(
+        1.117689, abs=1e-3
+    )
+
+    assert neb_summary["neb_results"]["trajectory_results"][1][
+        "energy"
+    ] == pytest.approx(-23.58905916097854, abs=1e-5)
+
+
+def test_neb_job_geodesic(setup_test_environment):
+    reactant, product = setup_test_environment
+
+    neb_summary = neb_job(
+        reactant, product, interpolation_method="geodesic", neb_kwargs={"max_steps": 10}
+    )
+
+    assert len(neb_summary["neb_results"]["trajectory_results"]) == 20
+    assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx(
+        0.8815, abs=1e-3
+    )
+    assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx(
+        1.117689, abs=1e-3
+    )
+
+    assert neb_summary["neb_results"]["trajectory_results"][1][
+        "energy"
+    ] == pytest.approx(-24.895280838012695, abs=0.05)
+
+
+def test_geodesic_job(setup_test_environment):
+    reactant, product = setup_test_environment
+
+    geodesic_summary = geodesic_job(reactant, product)
+    assert geodesic_summary["ts_atoms"].get_potential_energy() == pytest.approx(
+        -22.597125398584318, abs=1e-6
+    )
diff --git a/tests/core/runners/test_ase.py b/tests/core/runners/test_ase.py
index 3901cb00..cc38617a 100644
--- a/tests/core/runners/test_ase.py
+++ b/tests/core/runners/test_ase.py
@@ -2,16 +2,26 @@ from __future__ import annotations
 
 import glob
 import os
+from importlib.util import find_spec
 from logging import WARNING, getLogger
 from pathlib import Path
 from shutil import rmtree
 
 import numpy as np
 import pytest
+
+
+@pytest.fixture(scope="module", autouse=True)
+def set_seed():
+    np.random.seed(42)  # noqa: NPY002
+
+
+from ase import Atoms
 from ase.build import bulk, molecule
 from ase.calculators.emt import EMT
 from ase.calculators.lj import LennardJones
 from ase.io import read
+from ase.mep.neb import NEBOptimizer
 from ase.optimize import BFGS, BFGSLineSearch
 from ase.optimize.sciopt import SciPyFminBFGS
 
@@ -19,6 +29,9 @@ from quacc import JobFailure, change_settings, get_settings
 from quacc.runners._base import BaseRunner
 from quacc.runners.ase import Runner
 
+has_geodesic_interpolate = bool(find_spec("geodesic_interpolate"))
+test_files_path = Path(__file__).parent / "test_files"
+
 LOGGER = getLogger(__name__)
 LOGGER.propagate = True
 
@@ -48,6 +61,36 @@ def teardown_function():
             os.remove(os.path.join(results_dir, f))
 
 
+@pytest.fixture
+def setup_test_environment(tmp_path):
+    reactant = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.4835950817281542, -1.0145410211301968, -0.13209027203235943],
+            [0.8409564131524673, 0.018549610257914483, -0.07338809662321308],
+            [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891],
+            [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892],
+            [1.402180736662139, 0.944112416574632, -0.12179540364365492],
+            [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071],
+            [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239],
+        ],
+    )
+
+    product = Atoms(
+        symbols="CCHHCHH",
+        positions=[
+            [1.348003553501624, 0.4819311116778978, 0.2752537177143993],
+            [0.2386618286631742, -0.3433222966734429, 0.37705518940917926],
+            [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955],
+            [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953],
+            [0.3801391040059668, -1.3793340533058087, 0.71035902765307],
+            [1.9296265384257907, 0.622088341468767, 1.0901733942191298],
+            [-1.090815880212625, 1.0965111343610956, -0.23791518420660265],
+        ],
+    )
+    return reactant, product
+
+
 def test_base_runner(tmp_path, monkeypatch):
     monkeypatch.chdir(tmp_path)
     atoms = bulk("Cu")
@@ -249,3 +292,40 @@ def test_fn_hook(tmp_path, monkeypatch):
         Runner(bulk("Cu"), EMT()).run_opt(fn_hook=fn_hook)
     with pytest.raises(ValueError, match="Test error"):
         raise err.value.parent_error
+
+
+def test_run_neb(monkeypatch, tmp_path):
+    monkeypatch.chdir(tmp_path)
+    geodesic_path = test_files_path / "geodesic_path.xyz"
+    images = read(geodesic_path, index=":")
+
+    neb_kwargs = {"method": "aseneb", "precon": None}
+    dyn = Runner(images, EMT()).run_neb(optimizer=NEBOptimizer, neb_kwargs=neb_kwargs)
+    traj = read(dyn.trajectory.filename, index=":")
+
+    assert traj[-1].calc.results is not None
+    assert not os.path.exists(tmp_path / "opt.log")
+
+
+def test_run_neb2():
+    geodesic_path = test_files_path / "geodesic_path.xyz"
+
+    images = read(geodesic_path, index=":")
+
+    with pytest.raises(
+        ValueError, match="BFGSLineSearch is not allowed as optimizer with NEB."
+    ):
+        Runner(images, EMT()).run_neb(
+            optimizer=BFGSLineSearch, neb_kwargs={"method": "aseneb", "precon": None}
+        )
+
+
+def test_run_neb_raises_value_error_for_trajectory_kwarg():
+    images = [Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.74]])]
+
+    with pytest.raises(
+        ValueError, match="Quacc does not support setting the `trajectory` kwarg."
+    ):
+        Runner(images, EMT()).run_neb(
+            optimizer=NEBOptimizer, optimizer_kwargs={"trajectory": "some_traj.traj"}
+        )
diff --git a/tests/core/runners/test_files/geodesic_path.xyz b/tests/core/runners/test_files/geodesic_path.xyz
new file mode 100644
index 00000000..73c288b5
--- /dev/null
+++ b/tests/core/runners/test_files/geodesic_path.xyz
@@ -0,0 +1,90 @@
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.60210258      -0.85402158       0.29381286
+C       -0.28783951      -0.53322413      -0.44372747
+H       -0.95727408       0.42884906       0.07350937
+H        0.35481921       0.65326383       0.53310253
+C        0.83497720      -0.13893815      -0.61903514
+H       -0.06453974       0.87896748      -0.88297197
+H       -0.48224567      -0.43489651       1.04530982
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.59223789      -0.81802021       0.31135853
+C       -0.28604018      -0.47152900      -0.42593957
+H       -0.92390973       0.45417612       0.20964610
+H        0.38052911       0.69342674       0.54238992
+C        0.85025969      -0.11855947      -0.60382787
+H       -0.05043198       0.84149609      -1.07238207
+H       -0.56264480      -0.58099027       1.03875497
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.56752580      -0.76523536       0.32900863
+C       -0.29444225      -0.36150756      -0.39642003
+H       -0.88419505       0.46889917       0.41683865
+H        0.42167940       0.75256538       0.57420807
+C        0.85881363      -0.07344971      -0.57888790
+H       -0.01427599       0.74054417      -1.33086389
+H       -0.65510555      -0.76181608       0.98611646
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.53153649      -0.70795790       0.32135314
+C       -0.31196159      -0.20232825      -0.36082476
+H       -0.83211712       0.45957469       0.65330147
+H        0.46707007       0.80230819       0.65159567
+C        0.85547097       0.01215643      -0.55082250
+H        0.02479217       0.53705274      -1.58192589
+H       -0.73479098      -0.90080590       0.86732288
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.49710666      -0.64719683       0.28652142
+C       -0.32455217      -0.02770666      -0.32588575
+H       -0.76903102       0.41846514       0.85063704
+H        0.50991033       0.83413692       0.73959991
+C        0.84936453       0.11862994      -0.53249860
+H        0.03875338       0.26781854      -1.71408411
+H       -0.80155171      -0.96414705       0.69571009
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.46929389      -0.58407604       0.23873741
+C       -0.32819281       0.13632023      -0.29577230
+H       -0.70406684       0.34361913       0.98764228
+H        0.54736483       0.85455598       0.81084680
+C        0.84529882       0.22536273      -0.52705498
+H        0.02725281      -0.01111080      -1.71369497
+H       -0.85695070      -0.96467122       0.49929576
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.45166801      -0.52656333       0.20165987
+C       -0.32191784       0.26695481      -0.27071285
+H       -0.64751808       0.25267214       1.06356077
+H        0.57434755       0.87413911       0.85281977
+C        0.84453688       0.30577859      -0.53229994
+H       -0.00510047      -0.24219076      -1.62060929
+H       -0.89601604      -0.93079055       0.30558168
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.44176149      -0.48136401       0.18034320
+C       -0.31132747       0.36096930      -0.25222570
+H       -0.60331694       0.16140671       1.09768002
+H        0.58999854       0.89461629       0.87334640
+C        0.84542407       0.35773560      -0.54324267
+H       -0.04447183      -0.41002497      -1.49647966
+H       -0.91806787      -0.88333892       0.14057842
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.43559245      -0.44795824       0.16566337
+C       -0.30056794       0.42882307      -0.24077310
+H       -0.56621759       0.07575027       1.10766820
+H        0.59486995       0.91244138       0.88346800
+C        0.84607809       0.39270844      -0.55602468
+H       -0.07914315      -0.53074541      -1.37277954
+H       -0.93061181      -0.83101951       0.01277776
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.43009268      -0.42437008       0.15071523
+C       -0.29152973       0.47863977      -0.23705145
+H       -0.53028723      -0.00274016       1.10374152
+H        0.58906640       0.92480496       0.88885267
+C        0.84668896       0.41932394      -0.56813409
+H       -0.10017589      -0.61881108      -1.26384123
+H       -0.94385519      -0.77684734      -0.07428264
diff --git a/tests/core/schemas/test_ase.py b/tests/core/schemas/test_ase.py
index f32a6cd3..019b20ae 100644
--- a/tests/core/schemas/test_ase.py
+++ b/tests/core/schemas/test_ase.py
@@ -8,6 +8,7 @@ import pytest
 from ase.build import bulk, molecule
 from ase.calculators.emt import EMT
 from ase.io import read
+from ase.mep import NEB
 from ase.optimize import BFGS
 from ase.vibrations import Vibrations
 from maggma.stores import MemoryStore
@@ -429,3 +430,26 @@ def test_errors(tmp_path, monkeypatch):
     vib.run()
     with pytest.raises(ValueError, match="Unsupported thermo_method"):
         VibSummarize(vib, directory=tmp_path).vib_and_thermo("bad")
+
+
+def test_summarize_neb(monkeypatch, tmp_path):
+    monkeypatch.chdir(tmp_path)
+
+    # Read initial and final states:
+    images = read(FILE_DIR / "test_files" / "geodesic_path.xyz", index=":")
+    neb = NEB(images)
+    for image in images:
+        image.calc = EMT()
+    optimizer = BFGS(neb, trajectory="opt.traj")
+    optimizer.run(fmax=0.05)
+
+    neb_summary = Summarize().neb(
+        optimizer, len(images), n_iter_return=10, trajectory=read("opt.traj", index=":")
+    )
+    assert neb_summary["trajectory_results"][-2]["energy"] == pytest.approx(
+        1.0875330074934446, abs=1e-4
+    )
+
+    ts_atoms = neb_summary["ts_atoms"]
+    ts_atoms.calc = EMT()
+    assert ts_atoms.get_potential_energy() == pytest.approx(1.1603536513693768, 1e-4)
diff --git a/tests/core/schemas/test_files/geodesic_path.xyz b/tests/core/schemas/test_files/geodesic_path.xyz
new file mode 100644
index 00000000..73c288b5
--- /dev/null
+++ b/tests/core/schemas/test_files/geodesic_path.xyz
@@ -0,0 +1,90 @@
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.60210258      -0.85402158       0.29381286
+C       -0.28783951      -0.53322413      -0.44372747
+H       -0.95727408       0.42884906       0.07350937
+H        0.35481921       0.65326383       0.53310253
+C        0.83497720      -0.13893815      -0.61903514
+H       -0.06453974       0.87896748      -0.88297197
+H       -0.48224567      -0.43489651       1.04530982
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.59223789      -0.81802021       0.31135853
+C       -0.28604018      -0.47152900      -0.42593957
+H       -0.92390973       0.45417612       0.20964610
+H        0.38052911       0.69342674       0.54238992
+C        0.85025969      -0.11855947      -0.60382787
+H       -0.05043198       0.84149609      -1.07238207
+H       -0.56264480      -0.58099027       1.03875497
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.56752580      -0.76523536       0.32900863
+C       -0.29444225      -0.36150756      -0.39642003
+H       -0.88419505       0.46889917       0.41683865
+H        0.42167940       0.75256538       0.57420807
+C        0.85881363      -0.07344971      -0.57888790
+H       -0.01427599       0.74054417      -1.33086389
+H       -0.65510555      -0.76181608       0.98611646
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.53153649      -0.70795790       0.32135314
+C       -0.31196159      -0.20232825      -0.36082476
+H       -0.83211712       0.45957469       0.65330147
+H        0.46707007       0.80230819       0.65159567
+C        0.85547097       0.01215643      -0.55082250
+H        0.02479217       0.53705274      -1.58192589
+H       -0.73479098      -0.90080590       0.86732288
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.49710666      -0.64719683       0.28652142
+C       -0.32455217      -0.02770666      -0.32588575
+H       -0.76903102       0.41846514       0.85063704
+H        0.50991033       0.83413692       0.73959991
+C        0.84936453       0.11862994      -0.53249860
+H        0.03875338       0.26781854      -1.71408411
+H       -0.80155171      -0.96414705       0.69571009
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.46929389      -0.58407604       0.23873741
+C       -0.32819281       0.13632023      -0.29577230
+H       -0.70406684       0.34361913       0.98764228
+H        0.54736483       0.85455598       0.81084680
+C        0.84529882       0.22536273      -0.52705498
+H        0.02725281      -0.01111080      -1.71369497
+H       -0.85695070      -0.96467122       0.49929576
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.45166801      -0.52656333       0.20165987
+C       -0.32191784       0.26695481      -0.27071285
+H       -0.64751808       0.25267214       1.06356077
+H        0.57434755       0.87413911       0.85281977
+C        0.84453688       0.30577859      -0.53229994
+H       -0.00510047      -0.24219076      -1.62060929
+H       -0.89601604      -0.93079055       0.30558168
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.44176149      -0.48136401       0.18034320
+C       -0.31132747       0.36096930      -0.25222570
+H       -0.60331694       0.16140671       1.09768002
+H        0.58999854       0.89461629       0.87334640
+C        0.84542407       0.35773560      -0.54324267
+H       -0.04447183      -0.41002497      -1.49647966
+H       -0.91806787      -0.88333892       0.14057842
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.43559245      -0.44795824       0.16566337
+C       -0.30056794       0.42882307      -0.24077310
+H       -0.56621759       0.07575027       1.10766820
+H        0.59486995       0.91244138       0.88346800
+C        0.84607809       0.39270844      -0.55602468
+H       -0.07914315      -0.53074541      -1.37277954
+H       -0.93061181      -0.83101951       0.01277776
+7
+Properties=species:S:1:pos:R:3 pbc="F F F"
+C        0.43009268      -0.42437008       0.15071523
+C       -0.29152973       0.47863977      -0.23705145
+H       -0.53028723      -0.00274016       1.10374152
+H        0.58906640       0.92480496       0.88885267
+C        0.84668896       0.41932394      -0.56813409
+H       -0.10017589      -0.61881108      -1.26384123
+H       -0.94385519      -0.77684734      -0.07428264
diff --git a/tests/requirements-newtonnet.txt b/tests/requirements-newtonnet.txt
index 932e4e19..5603fb55 100644
--- a/tests/requirements-newtonnet.txt
+++ b/tests/requirements-newtonnet.txt
@@ -1,1 +1,3 @@
 newtonnet==1.1.1
+geodesic-interpolate @ git+https://github.com/virtualzx-nad/geodesic-interpolate.git
+sella==2.3.5