cleaned geodesic tests. Added in memory geodesic wrapper too.

src/neb_wrapper.py

@@ -35,7 +35,12 @@ def run_neb_method(
         max_steps (int, optional): maximum number of optimization steps allowed. Defaults to 1000.
         fmax_cutoff (float: optional): convergence cut-off criteria for the NEB optimization. Defaults to 1e-2.
     """
-    images = setup_images(logdir, xyz_r_p, xyz_ts, n_intermediate)
+    images = setup_images(
+        logdir,
+        xyz_r_p,
+        xyz_ts,
+        n_intermediate=n_intermediate,
+    )
 
     mep = NEB(
         images,

src/setup_images.py

@@ -1,8 +1,73 @@
 import os
 from calcs import calc
+from ase import Atoms
 from ase.io import read, write
 from sella_wrapper import sella_wrapper
 
+import logging
+import numpy as np
+from geodesic_interpolate.fileio import read_xyz, write_xyz
+from geodesic_interpolate.interpolation import redistribute
+from geodesic_interpolate.geodesic import Geodesic
+
+
+def geodesic_interpolate_wrapper(
+    r_p_atoms: Atoms,
+    nimages: int = 17,
+    sweep: bool = None,
+    output: str = "interpolated.xyz",
+    tol: float = 2e-3,
+    maxiter: int = 15,
+    microiter: int = 20,
+    scaling: float = 1.7,
+    friction: float = 1e-2,
+    dist_cutoff: float = 3,
+    save_raw: str = None):
+    """
+    Interpolates between two geometries and optimizes the path.
+
+    Parameters:
+    filename (str): XYZ file containing geometries.
+    nimages (int): Number of images. Default is 17.
+    sweep (bool): Sweep across the path optimizing one image at a time.
+                  Default is to perform sweeping updates if there are more than 35 atoms.
+    output (str): Output filename. Default is "interpolated.xyz".
+    tol (float): Convergence tolerance. Default is 2e-3.
+    maxiter (int): Maximum number of minimization iterations. Default is 15.
+    microiter (int): Maximum number of micro iterations for sweeping algorithm. Default is 20.
+    scaling (float): Exponential parameter for Morse potential. Default is 1.7.
+    friction (float): Size of friction term used to prevent very large change of geometry. Default is 1e-2.
+    dist_cutoff (float): Cut-off value for the distance between a pair of atoms to be included in the coordinate system. Default is 3.
+    save_raw (str): When specified, save the raw path after bisections but before smoothing. Default is None.
+    """
+    # Read the initial geometries.
+    symbols = r_p_atoms[0].get_chemical_symbols()
+
+    X = [conf.get_positions() for conf in r_p_atoms]
+
+    if len(X) < 2:
+        raise ValueError("Need at least two initial geometries.")
+
+    # First redistribute number of images. Perform interpolation if too few and subsampling if too many images are given
+    raw = redistribute(symbols, X, nimages, tol=tol * 5)
+    if save_raw is not None:
+        write_xyz(save_raw, symbols, raw)
+
+    # Perform smoothing by minimizing distance in Cartesian coordinates with redundant internal metric
+    # to find the appropriate geodesic curve on the hyperspace.
+    smoother = Geodesic(symbols, raw, scaling, threshold=dist_cutoff, friction=friction)
+    if sweep is None:
+        sweep = len(symbols) > 35
+    try:
+        if sweep:
+            smoother.sweep(tol=tol, max_iter=maxiter, micro_iter=microiter)
+        else:
+            smoother.smooth(tol=tol, max_iter=maxiter)
+    finally:
+        # Save the smoothed path to output file. try block is to ensure output is saved if one ^C the process, or there is an error
+        write_xyz(output, symbols, smoother.path)
+    return symbols, smoother.path
+
 
 def setup_images(
         logdir: str,
@@ -39,9 +104,9 @@ def setup_images(
         write(r_p_path, [reactant.copy(), product.copy()])
 
         # Generate intermediate images using geodesic interpolation
-        output_xyz = os.path.join(logdir, 'output.xyz')
-        os.system(f'geodesic_interpolate {r_p_path} --output {output_xyz} --nimages {n_intermediate}')
-        images = read(output_xyz, index=':')
+        symbols, smoother_path =\
+            geodesic_interpolate_wrapper([reactant.copy(), product.copy()])
+        images = [Atoms(symbols=symbols, positions=conf) for conf in smoother_path]
 
         # Calculate energies and forces for each intermediate image
         for image in images:

tests/test_setup_images.py

@@ -1,8 +1,23 @@
 import os
 import pytest
 from ase import Atoms
-from ase.io import write
+from ase.io import read, write
 from setup_images import setup_images
+from setup_images import geodesic_interpolate_wrapper
+
+
+def test_geodesic_interpolate_wrapper(setup_test_environment):
+    logdir, xyz_r_p = setup_test_environment
+    atoms_object = read(xyz_r_p, index=':')
+    symbols, smoother_path = geodesic_interpolate_wrapper(
+        atoms_object
+    )
+    # assert output == 1
+    # assert symbols == 1
+    assert smoother_path[1][0][0] == pytest.approx(
+        1.36055556030,
+        abs=1e-3,
+    )
 
 
 def test_setup_images(setup_test_environment):
@@ -22,14 +37,35 @@ def test_setup_images(setup_test_environment):
     assert os.path.isfile(logdir / 'reactant_opt.traj'), "Reactant optimization file not found"
     assert os.path.isfile(logdir / 'product_opt.traj'), "Product optimization file not found"
     assert os.path.isfile(logdir / 'r_p.xyz'), "Reactant-Product file not found"
-    assert os.path.isfile(logdir / 'output.xyz'), "Intermediate images file not found"
-    assert os.path.isfile(logdir / 'geodesic_path.xyz'), "Geodesic path file not found"
+
+    assert images[1].get_positions()[0][0] == pytest.approx(
+        1.439202,
+        abs=1e-2,
+    )
 
     # Check energies and forces
     for image in images:
         assert 'energy' in image.info, "Energy not found in image info"
         assert 'forces' in image.arrays, "Forces not found in image arrays"
 
+    assert images[1].get_potential_energy() == pytest.approx(
+        -29.9266674,
+        abs=1,
+    )
+    "Error in first intermediate image's energy for the geodesic path"
+
+    assert images[0].get_potential_energy() == pytest.approx(
+        -32.98172029,
+        abs=1,
+    )
+    "Error in reactant energy prediction for geodesic path."
+
+    assert images[-1].get_potential_energy() == pytest.approx(
+        -25.1172894,
+        abs=1,
+    )
+    "Error in product energy prediction for geodesic path."
+
 
 if __name__ == "__main__":
     pytest.main()