diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
index 63d6113..071b3a5 100644
--- a/.github/workflows/python-package-conda.yml
+++ b/.github/workflows/python-package-conda.yml
@@ -5,30 +5,37 @@ on: [push]
 jobs:
   build-linux:
     runs-on: ubuntu-latest
-    strategy:
-      max-parallel: 5
 
     steps:
     - uses: actions/checkout@v4
-    - name: Set up Python 3.10
-      uses: actions/setup-python@v3
+
+    - name: Set up Miniconda
+      uses: conda-incubator/setup-miniconda@v3
       with:
-        python-version: '3.10'
-    - name: Add conda to system path
-      run: |
-        # $CONDA is an environment variable pointing to the root of the miniconda directory
-        echo $CONDA/bin >> $GITHUB_PATH
-    - name: Install dependencies
-      run: |
-        conda env update --file environment.yml --name base
+        auto-update-conda: true
+        environment-file: environment.yml
+        activate-environment: base
+        auto-activate-base: true
+
     - name: Lint with flake8
       run: |
-        conda install flake8
-        # stop the build if there are Python syntax errors or undefined names
-        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+        conda run -n base python -m pip install flake8
+        conda run -n base flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
+        conda run -n base flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
+
     - name: Test with pytest
       run: |
-        conda install pytest
-        pytest -q -m "not slow"
+        conda run -n base python -m pip install pytest
+        conda run -n base pytest -q -m "not slow"
+
+    - name: Run simulation
+      run: |
+        mkdir -p outputs
+        conda run -n base python iro2_test_1.py
+        ls -R .
+
+    - name: Upload simulation outputs
+      uses: actions/upload-artifact@v4
+      with:
+        name: irO2-simulation-outputs
+        path: outputs/
diff --git a/.gitignore b/.gitignore
index c291fe1..65ebff9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -17,8 +17,8 @@ build/
 .ipynb_checkpoints/
 
 # Common simulation outputs (adjust as needed)
-*.traj
+#*.traj
 *.xyz
 *.cif
-*.log
+#*.log
 *.out
diff --git a/iro2_test_1.py b/iro2_test_1.py
index 7fa1b1b..d2f73e7 100644
--- a/iro2_test_1.py
+++ b/iro2_test_1.py
@@ -1,119 +1,119 @@
-# pip install ase
-# conda config --add channels conda-forge
-# conda install xtb-python
-# conda search xtb-python --channel conda-forge
-# conda install qcelemental ase
-from pathlib import Path
-from ase.io import read
-import numpy as np
-from ase.constraints import FixAtoms
-from ase import Atom
-from ase.visualize import view
-from ase.io import write
-from ase.visualize import view
-from ase.data import covalent_radii
-from ase.constraints import FixAtoms
-from ase.neighborlist import neighbor_list
-from xtb.ase.calculator import XTB
-from ase.optimize import BFGS
-
-Path("outputs").mkdir(exist_ok=True)
-
-
-def run_relaxation(
-    qe_input: str = "slab_clean_2x2.in",
-    traj: str = "relax_H.traj",
-    log: str = "relax_H.log",
-    fmax: float = 0.05,
-    cut: float = 1.3,
-):
-
-    covalent_radii[1] = 0.6
-
-    slab = read("slab_clean_2x2.in", format="espresso-in")
-    # xTB does not support PBC multipoles in this mode → treat slab as cluster
-    slab.set_pbc((False, False, False))
-    print(slab)
-    print("Atoms:", len(slab))
-    print("Cell (Å):\n", slab.cell)
-    print("PBC:", slab.pbc)
-
-    # Get atomic symbols and positions
-    symbols = slab.get_chemical_symbols()
-    positions = slab.get_positions()
-
-    # Find top of slab
-    z_positions = positions[:, 2]
-    z_top = z_positions.max()
-
-    print(f"Top of slab z = {z_top:.3f} Å")
-
-    # Find surface O atoms within 1.5 Å of top
-    surface_O_indices = [
-        i for i, s in enumerate(symbols)
-        if s == "O" and (z_top - positions[i, 2]) < 1.5
-    ]
-
-    print("Candidate surface O atoms:")
-    for i in surface_O_indices:
-        x, y, z = positions[i]
-        print(f"Index {i}: position = ({x:.3f}, {y:.3f}, {z:.3f})")
-
-    #Adding H* on top of Oxygen 20
-
-    o_idx = 20
-    O_pos = slab.positions[o_idx]
-
-    # Typical O–H bond length
-    OH_distance = 1.0  # Å
-
-    # Place H above O along +z
-    H_pos = O_pos + np.array([0.0, 0.0, OH_distance])
-
-    # Add H atom
-    slab.append(Atom("H", position=H_pos))
-
-    #for atom in slab:
-        #if atom.symbol == "H":
-            #atom.tag = 1
+# # pip install ase
+# # conda config --add channels conda-forge
+# # conda install xtb-python
+# # conda search xtb-python --channel conda-forge
+# # conda install qcelemental ase
+# #from pathlib import Path
+# from ase.io import read
+# import numpy as np
+# from ase.constraints import FixAtoms
+# from ase import Atom
+# from ase.visualize import view
+# from ase.io import write
+# from ase.visualize import view
+# from ase.data import covalent_radii
+# from ase.constraints import FixAtoms
+# from ase.neighborlist import neighbor_list
+# from xtb.ase.calculator import XTB
+# from ase.optimize import BFGS
+
+# Path("outputs").mkdir(exist_ok=True)
+
+
+# def run_relaxation(
+#     qe_input: str = "slab_clean_2x2.in",
+#     traj: str = "relax_H.traj",
+#     log: str = "relax_H.log",
+#     fmax: float = 0.05,
+#     cut: float = 1.3,
+# ):
+
+#     covalent_radii[1] = 0.6
+
+#     slab = read("slab_clean_2x2.in", format="espresso-in")
+#     # xTB does not support PBC multipoles in this mode → treat slab as cluster
+#     slab.set_pbc((False, False, False))
+#     print(slab)
+#     print("Atoms:", len(slab))
+#     print("Cell (Å):\n", slab.cell)
+#     print("PBC:", slab.pbc)
+
+#     # Get atomic symbols and positions
+#     symbols = slab.get_chemical_symbols()
+#     positions = slab.get_positions()
+
+#     # Find top of slab
+#     z_positions = positions[:, 2]
+#     z_top = z_positions.max()
+
+#     print(f"Top of slab z = {z_top:.3f} Å")
+
+#     # Find surface O atoms within 1.5 Å of top
+#     surface_O_indices = [
+#         i for i, s in enumerate(symbols)
+#         if s == "O" and (z_top - positions[i, 2]) < 1.5
+#     ]
+
+#     print("Candidate surface O atoms:")
+#     for i in surface_O_indices:
+#         x, y, z = positions[i]
+#         print(f"Index {i}: position = ({x:.3f}, {y:.3f}, {z:.3f})")
+
+#     #Adding H* on top of Oxygen 20
+
+#     o_idx = 20
+#     O_pos = slab.positions[o_idx]
+
+#     # Typical O–H bond length
+#     OH_distance = 1.0  # Å
+
+#     # Place H above O along +z
+#     H_pos = O_pos + np.array([0.0, 0.0, OH_distance])
+
+#     # Add H atom
+#     slab.append(Atom("H", position=H_pos))
+
+#     #for atom in slab:
+#         #if atom.symbol == "H":
+#             #atom.tag = 1
         
-    #set_highlight_atoms([i for i,a in enumerate(slab) if a.symbol=="H"])
+#     #set_highlight_atoms([i for i,a in enumerate(slab) if a.symbol=="H"])
 
-    print("Total atoms:", len(slab))
-    print("Last atom symbol:", slab[-1].symbol)
-    print("Last atom position:", slab[-1].position)
+#     print("Total atoms:", len(slab))
+#     print("Last atom symbol:", slab[-1].symbol)
+#     print("Last atom position:", slab[-1].position)
 
-    print(f"Added H at position {H_pos}")
-    print("Total atoms now:", len(slab))
+#     print(f"Added H at position {H_pos}")
+#     print("Total atoms now:", len(slab))
 
-    write("slab_with_H.xyz", slab)
+#     write("slab_with_H.xyz", slab)
 
-    view(slab)
+#     view(slab)
 
-    z = slab.positions[:, 2]
-    z_freeze = 20.0   # Å (tune if needed)
+#     z = slab.positions[:, 2]
+#     z_freeze = 20.0   # Å (tune if needed)
 
-    freeze_mask = z < z_freeze
-    slab.set_constraint(FixAtoms(mask=freeze_mask))
+#     freeze_mask = z < z_freeze
+#     slab.set_constraint(FixAtoms(mask=freeze_mask))
 
-    print("Frozen atoms:", int(freeze_mask.sum()), "/", len(slab))
+#     print("Frozen atoms:", int(freeze_mask.sum()), "/", len(slab))
 
-    iH = len(slab) - 1  # last atom is H
-    cut = 1.3           # Å check radius around H
+#     iH = len(slab) - 1  # last atom is H
+#     cut = 1.3           # Å check radius around H
 
-    i, j, d = neighbor_list("ijd", slab, cutoff=[cut]*len(slab))
-    close = [(jj, dd) for ii, jj, dd in zip(i, j, d) if ii == iH]
+#     i, j, d = neighbor_list("ijd", slab, cutoff=[cut]*len(slab))
+#     close = [(jj, dd) for ii, jj, dd in zip(i, j, d) if ii == iH]
 
 
 
-    print("Neighbors within 1.3 Å of H:", close)
+#     print("Neighbors within 1.3 Å of H:", close)
 
-    # Optimize structure with xTB
-    slab.calc = XTB(method="GFN2-xTB")
-    opt = BFGS(slab, trajectory=traj, logfile=log)
-    opt.run(fmax=fmax)
+#     # Optimize structure with xTB
+#     slab.calc = XTB(method="GFN2-xTB")
+#     opt = BFGS(slab, trajectory=traj, logfile=log)
+#     opt.run(fmax=fmax)
 
-    return slab, close
+#     return slab, close
 
-if __name__ == "__main__":
-    run_relaxation()
+# if __name__ == "__main__":
+#     run_relaxation()
diff --git a/scripts/iro2_slab_setup.py b/scripts/iro2_slab_setup.py
new file mode 100644
index 0000000..b5c204e
--- /dev/null
+++ b/scripts/iro2_slab_setup.py
@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+from ase import Atom
+from ase.constraints import FixAtoms
+from ase.io import read, write
+from ase.neighborlist import neighbor_list
+
+
+def setup_structure(
+    input_file: str = "inputs/slab_clean_2x2.in",
+    o_index: int = 20,
+    oh_distance: float = 1.0,
+    z_freeze: float = 20.0,
+    neighbor_cutoff: float = 1.5,
+    outputs_dir: str = "outputs",
+):
+    """
+    Read slab, place H above a chosen O, freeze atoms below z_freeze,
+    and write outputs for the next step.
+    """
+    input_path = Path(input_file)
+    if not input_path.exists():
+        raise FileNotFoundError(f"Input file not found: {input_path}")
+
+    slab = read(str(input_path))
+
+    n0 = len(slab)
+    if not (0 <= o_index < n0):
+        raise IndexError(f"o_index={o_index} out of range for {n0} atoms (0..{n0-1})")
+
+    # Slab z-extent (for diagnostics)
+    zmax = slab.positions[:, 2].max()
+
+    # Position of selected O
+    o_pos = slab.positions[o_index].copy()
+
+    # Place H above O along +z
+    h_pos = o_pos + np.array([0.0, 0.0, float(oh_distance)])
+
+    # Sanity check: ensure H is above the slab top
+    if h_pos[2] <= zmax:
+        print(
+            f"[warn] H z={h_pos[2]:.3f} Å is not above slab top zmax={zmax:.3f} Å. "
+            "Raising H to be above slab."
+        )
+    h_pos[2] = zmax + float(oh_distance)
+
+    slab.append(Atom("H", position=h_pos))
+    h_index = len(slab) - 1
+
+    slab.append(Atom("H", position=h_pos))
+    h_index = len(slab) - 1
+
+    # Freeze atoms below z_freeze
+    freeze_mask = slab.positions[:, 2] < float(z_freeze)
+    slab.set_constraint(FixAtoms(mask=freeze_mask))
+
+    # Neighbors of H (for sanity check)
+    i, j, d = neighbor_list("ijd", slab, cutoff=float(neighbor_cutoff))
+    neigh = [(int(jj), float(dd)) for ii, jj, dd in zip(i, j, d) if int(ii) == h_index]
+    neigh_sorted = sorted(neigh, key=lambda x: x[1])
+
+    outdir = Path(outputs_dir)
+    (outdir / "results").mkdir(parents=True, exist_ok=True)
+
+    # Write structures
+    write(str(outdir / "slab_with_H.traj"), slab)
+    write(str(outdir / "slab_with_H.xyz"), slab)
+
+    meta = {
+        "input_file": str(input_path),
+        "num_atoms_initial": int(n0),
+        "num_atoms_total": int(len(slab)),
+        "o_index": int(o_index),
+        "h_index": int(h_index),
+        "oh_distance_A": float(oh_distance),
+        "z_freeze_A": float(z_freeze),
+        "neighbor_cutoff_A": float(neighbor_cutoff),
+        "o_position_A": [float(x) for x in o_pos],
+        "h_position_A": [float(x) for x in h_pos],
+        "num_frozen": int(np.count_nonzero(freeze_mask)),
+        "neighbors_of_H": [{"index": idx, "distance_A": dist} for idx, dist in neigh_sorted],
+    }
+    (outdir / "results" / "slab_setup_metadata.json").write_text(json.dumps(meta, indent=2))
+
+    print(f"[slab_setup] Loaded {input_path} (atoms={n0})")
+    print(f"[slab_setup] Placed H at index {h_index} above O index {o_index} by {oh_distance} Å")
+    print(f"[slab_setup] Frozen atoms: {meta['num_frozen']} (z < {z_freeze} Å)")
+    print(f"[slab_setup] H neighbors within {neighbor_cutoff} Å: {neigh_sorted}")
+
+    return slab, meta
+
+
+def main():
+    p = argparse.ArgumentParser(description="IrO2 slab setup: place H above O, freeze bottom layers, write outputs/")
+    p.add_argument("--input", default="inputs/slab_clean_2x2.in")
+    p.add_argument("--o-index", type=int, default=20, help="0-based index of O atom to adsorb H onto")
+    p.add_argument("--oh-dist", type=float, default=1.0)
+    p.add_argument("--z-freeze", type=float, default=20.0)
+    p.add_argument("--neighbor-cutoff", type=float, default=1.5)
+    p.add_argument("--outputs", default="outputs")
+    args = p.parse_args()
+
+    setup_structure(
+        input_file=args.input,
+        o_index=args.o_index,
+        oh_distance=args.oh_dist,
+        z_freeze=args.z_freeze,
+        neighbor_cutoff=args.neighbor_cutoff,
+        outputs_dir=args.outputs,
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/optimization_iro2_3.py b/scripts/optimization_iro2_3.py
new file mode 100644
index 0000000..20fa13e
--- /dev/null
+++ b/scripts/optimization_iro2_3.py
@@ -0,0 +1,106 @@
+import json
+from pathlib import Path
+from ase.io import read
+import numpy as np
+
+def analyze_result(result_dir: str = "results"):
+    #Analyze optimization results and print summary
+    result_path = Path(result_dir)
+    
+    # Find all result JSON files
+    result_files = list(result_path.glob("*_results.json"))
+    
+    if not result_files:
+        print(f"No results found in {result_dir}/")
+        return
+    
+    print(f"\n{'='*70}")
+    print(f"OPTIMIZATION RESULTS SUMMARY")
+    print(f"{'='*70}\n")
+    
+    for result_file in sorted(result_files):
+        with open(result_file) as f:
+            data = json.load(f)
+        
+        print(f"File: {result_file.name}")
+        print(f"  Structure:     {Path(data['structure_file']).name}")
+        print(f"  Method:        {data['method']}")
+        print(f"  Steps:         {data['n_steps']}/{data['max_steps']}")
+        print(f"  Time:          {data['time_minutes']:.1f} min")
+        print(f"  Energy:        {data['e_initial']:.4f} → {data['e_final']:.4f} eV")
+        print(f"  ΔE:            {data['delta_e']:.4f} eV")
+        print(f"  fmax:          {data['fmax_initial']:.4f} → {data['fmax_final']:.4f} eV/Å")
+        print(f"  Converged:     {'YES' if data['converged'] else 'NO'}")
+        print()
+
+
+def compare_structures(initial_xyz: str, final_xyz: str):
+    """
+    Compare initial and final structures
+    """
+    initial = read(initial_xyz)
+    final = read(final_xyz)
+    
+    # Find H atom (last atom)
+    h_idx = len(initial) - 1
+    
+    h_init = initial[h_idx].position
+    h_final = final[h_idx].position
+    h_displacement = np.linalg.norm(h_final - h_init)
+    
+    print(f"\n{'='*70}")
+    print(f"STRUCTURE COMPARISON")
+    print(f"{'='*70}")
+    print(f"Initial: {initial_xyz}")
+    print(f"Final:   {final_xyz}")
+    print(f"\nH atom displacement:")
+    print(f"  Initial position: ({h_init[0]:.3f}, {h_init[1]:.3f}, {h_init[2]:.3f})")
+    print(f"  Final position:   ({h_final[0]:.3f}, {h_final[1]:.3f}, {h_final[2]:.3f})")
+    print(f"  Total movement:   {h_displacement:.3f} Å")
+    print(f"{'='*70}\n")
+    
+    # Calculate RMSD for all atoms
+    rmsd = np.sqrt(((initial.positions - final.positions)**2).sum(axis=1).mean())
+    print(f"Overall RMSD: {rmsd:.4f} Å")
+
+
+def visualize_trajectory(traj_file: str):
+    """
+    Quick stats from trajectory
+    """
+    from ase.io import read as read_traj
+    
+    traj = read_traj(traj_file, index=":")
+    
+    energies = [atoms.get_potential_energy() for atoms in traj]
+    
+    print(f"\n{'='*70}")
+    print(f"TRAJECTORY ANALYSIS")
+    print(f"{'='*70}")
+    print(f"File: {traj_file}")
+    print(f"Steps: {len(traj)}")
+    print(f"Energy range: {min(energies):.4f} to {max(energies):.4f} eV")
+    print(f"Total ΔE: {energies[-1] - energies[0]:.4f} eV")
+    print(f"{'='*70}\n")
+
+
+if __name__ == "__main__":
+    import argparse
+    
+    parser = argparse.ArgumentParser(description="Analyze optimization results")
+    parser.add_argument("--summary", action="store_true", help="Show all results summary")
+    parser.add_argument("--compare", nargs=2, metavar=("INITIAL", "FINAL"), 
+                        help="Compare two structures")
+    parser.add_argument("--traj", help="Analyze trajectory file")
+    parser.add_argument("--dir", default="results", help="Results directory")
+    
+    args = parser.parse_args()
+    
+    if args.summary or (not args.compare and not args.traj):
+        analyze_result(args.dir)
+    
+    if args.compare:
+        compare_structures(args.compare[0], args.compare[1])
+    
+    if args.traj:
+        visualize_trajectory(args.traj)
diff --git a/scripts/optimization_iro2_test.py b/scripts/optimization_iro2_test.py
new file mode 100644
index 0000000..9aa5957
--- /dev/null
+++ b/scripts/optimization_iro2_test.py
@@ -0,0 +1,189 @@
+from pathlib import Path
+from ase.io import read, write
+from xtb.ase.calculator import XTB
+from ase.optimize import BFGS
+import numpy as np
+import json
+import time
+
+def run_optimization(
+    structure_file: str,
+    fmax: float = 0.05,
+    max_steps: int = 200,
+    method: str = "GFN2-xTB",
+    output_dir: str = "results",
+):
+    # Create output directory
+    Path(output_dir).mkdir(exist_ok=True)
+    
+    # Start timing
+    start_time = time.time()
+    
+    print(f"{'='*70}")
+    print(f"XTB GEOMETRY OPTIMIZATION")
+    print(f"{'='*70}")
+    print(f"Structure: {structure_file}")
+    print(f"Method: {method}")
+    print(f"Convergence: fmax < {fmax} eV/Å")
+    print(f"Max steps: {max_steps}")
+    print(f"{'='*70}\n")
+    
+    # Load structure
+    print("Loading structure...")
+    slab = read(structure_file)
+    
+    # Check constraints
+    if slab.constraints:
+        n_constrained = sum(c.get_indices().size for c in slab.constraints)
+        print(f" Constraints loaded: {n_constrained} frozen atoms")
+    else:
+        print("WARNING: No constraints found!")
+    
+    print(f"Total atoms: {len(slab)}")
+    print(f"Mobile atoms: {len(slab) - n_constrained}")
+    
+    # Setup xTB calculator
+    print(f"\nSetting up {method} calculator...")
+    slab.calc = XTB(method=method)
+    
+    # Initial energy and forces
+    print("\nCalculating initial state...")
+    e_initial = slab.get_potential_energy()
+    forces = slab.get_forces()
+    fmax_initial = np.sqrt((forces**2).sum(axis=1)).max()
+    
+    print(f"  Initial energy: {e_initial:.6f} eV")
+    print(f"  Initial fmax:   {fmax_initial:.6f} eV/Å")
+    
+    if fmax_initial < fmax:
+        print(f"\n Already converged! No optimization needed.")
+        final_file = f"{output_dir}/{Path(structure_file).stem}_final.xyz"
+        write(final_file, slab)
+        return slab, {
+            "converged": True,
+            "initial": True,
+            "e_initial": e_initial,
+            "e_final": e_initial,
+            "fmax_initial": fmax_initial,
+            "fmax_final": fmax_initial,
+            "n_steps": 0,
+            "time_seconds": time.time() - start_time,
+        }
+    
+    # Setup output files
+    base_name = Path(structure_file).stem
+    traj_file = f"{output_dir}/{base_name}_opt.traj"
+    log_file = f"{output_dir}/{base_name}_opt.log"
+    
+    print(f"\n{'='*70}")
+    print(f"Starting optimization...")
+    print(f"{'='*70}")
+    print(f"Trajectory: {traj_file}")
+    print(f"Log file:   {log_file}")
+    print(f"{'='*70}\n")
+    
+    # Run optimization
+    opt = BFGS(slab, trajectory=traj_file, logfile=log_file)
+    
+    try:
+        opt.run(fmax=fmax, steps=max_steps)
+        converged = True
+    except Exception as e:
+        print(f"\n  Optimization stopped: {e}")
+        converged = False
+    
+    # Get final results
+    e_final = slab.get_potential_energy()
+    forces_final = slab.get_forces()
+    fmax_final = np.sqrt((forces_final**2).sum(axis=1)).max()
+    n_steps = opt.get_number_of_steps()
+    elapsed_time = time.time() - start_time
+    
+    # Print results
+    print(f"\n{'='*70}")
+    print(f"OPTIMIZATION COMPLETE")
+    print(f"{'='*70}")
+    print(f"Steps taken:        {n_steps}")
+    print(f"Time elapsed:       {elapsed_time:.1f} seconds ({elapsed_time/60:.1f} min)")
+    print(f"Initial energy:     {e_initial:.6f} eV")
+    print(f"Final energy:       {e_final:.6f} eV")
+    print(f"Energy change:      {e_final - e_initial:.6f} eV")
+    print(f"Initial fmax:       {fmax_initial:.6f} eV/Å")
+    print(f"Final fmax:         {fmax_final:.6f} eV/Å")
+    print(f"Converged:          {'YES' if fmax_final < fmax else 'NO'}")
+    print(f"{'='*70}")
+    
+    # Save final structure
+    final_file = f"{output_dir}/{base_name}_final.xyz"
+    write(final_file, slab)
+    print(f"\nFinal structure saved to: {final_file}")
+    
+    # Save results summary
+    results = {
+        "structure_file": structure_file,
+        "method": method,
+        "fmax_target": fmax,
+        "max_steps": max_steps,
+        "n_steps": n_steps,
+        "converged": fmax_final < fmax,
+        "e_initial": float(e_initial),
+        "e_final": float(e_final),
+        "delta_e": float(e_final - e_initial),
+        "fmax_initial": float(fmax_initial),
+        "fmax_final": float(fmax_final),
+        "time_seconds": elapsed_time,
+        "time_minutes": elapsed_time / 60,
+    }
+    
+    results_file = f"{output_dir}/{base_name}_results.json"
+    with open(results_file, 'w') as f:
+        json.dump(results, f, indent=2)
+    
+    print(f"Results summary saved to: {results_file}")
+    print(f"{'='*70}\n")
+    
+    return slab, results
+
+
+#Quick test with limited steps
+
+def quick_test(structure_file: str, steps: int = 5):
+ 
+    print(f"\n{'='*70}")
+    print(f"QUICK TEST MODE - {steps} STEPS ONLY")
+    print(f"{'='*70}\n")
+    
+    return run_optimization(
+        structure_file,
+        fmax=0.05,
+        max_steps=steps,
+        output_dir="test_results"
+    )
+
+
+if __name__ == "__main__":
+    import argparse
+    
+    parser = argparse.ArgumentParser(description="Run xTB optimization")
+    parser.add_argument("structure", help="Prepared structure file (.xyz)")
+    parser.add_argument("--fmax", type=float, default=0.05, help="Force convergence (eV/Å)")
+    parser.add_argument("--steps", type=int, default=200, help="Maximum steps")
+    parser.add_argument("--method", default="GFN2-xTB", choices=["GFN1-xTB", "GFN2-xTB"], 
+                        help="xTB method (GFN1 is faster)")
+    parser.add_argument("--test", action="store_true", help="Quick test mode (5 steps)")
+    parser.add_argument("--output", default="results", help="Output directory")
+    
+    args = parser.parse_args()
+    
+    if args.test:
+        slab, results = quick_test(args.structure, steps=5)
+    else:
+        slab, results = run_optimization(
+            args.structure,
+            fmax=args.fmax,
+            max_steps=args.steps,
+            method=args.method,
+            output_dir=args.output,
+        )
+    
+    print("\n Done!")