Added larger saved displacement amplitudes

README.md

@@ -93,6 +93,9 @@ graphrc calculation.out --graph
 
 # Save structures for IRC calculations
 graphrc calculation.out --save-displacement
+
+# Save with larger displacement (values >4 extrapolate beyond the trajectory using the normal mode vector)
+graphrc calculation.out --save-displacement -ds 8
 ```
 
 ---
@@ -503,7 +506,8 @@ output options:
   --save-displacement, -sd
                         Save displaced structure pair
   --displacement-scale DISPLACEMENT_SCALE, -ds DISPLACEMENT_SCALE
-                        Displacement level (1-4, ~0.2-0.8 amplitude) (default: 1)
+                        Displacement scale for saved displaced structures (1-4 uses trajectory frames directly, ~0.2-0.8
+                        amplitude; >4 extrapolates beyond the trajectory using the normal mode vector) (default: 1)
   --no-save             Do not save trajectory to disk (keep in memory only)
 ```
 

src/graphrc/api.py

@@ -96,22 +96,23 @@ def load_trajectory(
             print(f"Reading trajectory from {basename}")
         frames = read_xyz_trajectory(input_file)
         trajectory_file = input_file
-        return {"frames": frames, "frequencies": None, "trajectory_file": trajectory_file}
+        return {"frames": frames, "frequencies": None, "trajectory_file": trajectory_file, "mode_displacement": None}
 
     # QM output file - ORCA files use the direct parser; everything else uses cclib
+    mode_displacement = None
     if is_orca_output(input_file):
         try:
             if print_output:
                 print(f"\nParsing {basename} with direct ORCA output parser...")
-            frequencies, trajectory_string = parse_orca_output(input_file, mode)
+            frequencies, trajectory_string, mode_displacement = parse_orca_output(input_file, mode)
         except Exception as e:
             if print_output:
                 print(f"Direct ORCA parsing failed ({e}), falling back to cclib...")
-            frequencies, trajectory_string = parse_cclib_output(input_file, mode)
+            frequencies, trajectory_string, mode_displacement = parse_cclib_output(input_file, mode)
     else:
         if print_output:
             print(f"\nParsing {basename} with cclib...")
-        frequencies, trajectory_string = parse_cclib_output(input_file, mode)
+        frequencies, trajectory_string, mode_displacement = parse_cclib_output(input_file, mode)
 
     # Convert string to frames
     frames = parse_xyz_string_to_frames(trajectory_string)
@@ -123,7 +124,12 @@ def load_trajectory(
         if print_output:
             print(f"Saved trajectory to {os.path.basename(trajectory_file)}")
 
-    return {"frames": frames, "frequencies": frequencies, "trajectory_file": trajectory_file}
+    return {
+        "frames": frames,
+        "frequencies": frequencies,
+        "trajectory_file": trajectory_file,
+        "mode_displacement": mode_displacement,
+    }
 
 
 def run_vib_analysis(
@@ -242,6 +248,7 @@ def run_vib_analysis(
     )
 
     frames = trajectory_data["frames"]
+    mode_displacement = trajectory_data.get("mode_displacement")
 
     # Now print loading info after metadata will be displayed
     if print_output:
@@ -381,6 +388,7 @@ def run_vib_analysis(
             output_prefix=output_prefix,
             scale=displacement_scale,
             max_level=config.MAX_DISPLACEMENT_LEVEL,
+            mode_displacement=mode_displacement,
             print_output=print_output,
         )
 

src/graphrc/characterize.py

@@ -444,7 +444,7 @@ def characterize_vib_mode(
     has_angles = len(angle_changes) > 0
     has_dihedrals = len(dihedral_changes) > 0
 
-    result = {
+    result: Dict[str, Any] = {
         "has_bond_changes": has_bonds,
         "has_angle_changes": has_angles,
         "has_dihedral_changes": has_dihedrals,

src/graphrc/cli.py

@@ -145,7 +145,9 @@ def main():
         "-ds",
         type=int,
         default=config.DEFAULT_DISPLACEMENT_LEVEL,
-        help=f"Displacement level (1-{config.MAX_DISPLACEMENT_LEVEL}, ~0.2-0.8 amplitude) "
+        help=f"Displacement scale for saved displaced structures "
+        f"(1-{config.MAX_DISPLACEMENT_LEVEL} uses trajectory frames directly, ~0.2-0.8 amplitude; "
+        f">{config.MAX_DISPLACEMENT_LEVEL} extrapolates beyond the trajectory) "
         f"(default: {config.DEFAULT_DISPLACEMENT_LEVEL})",
     )
     output_group.add_argument(

src/graphrc/convert.py

@@ -87,7 +87,7 @@ def parse_cclib_output(output_file, mode):
         for sym, coord in zip(atom_symbols, displaced):
             trj_data += f"{sym} {coord[0]:.6f} {coord[1]:.6f} {coord[2]:.6f}\n"
 
-    return freqs, trj_data
+    return freqs, trj_data, displacement
 
 
 def parse_orca_output(orca_file: str, mode: int):
@@ -242,7 +242,7 @@ def parse_orca_output(orca_file: str, mode: int):
         for sym, coord in zip(atom_symbols, displaced):
             trj_data += f"{sym} {coord[0]:.6f} {coord[1]:.6f} {coord[2]:.6f}\n"
 
-    return freqs, trj_data
+    return freqs, trj_data, displacement
 
 
 def parse_xyz_string_to_frames(trj_data_str: str) -> List[Dict[str, Any]]:

src/graphrc/utils.py

@@ -222,36 +222,38 @@ def save_displacement_pair(
     output_prefix: str,
     scale: int = 1,
     max_level: int = 4,
+    mode_displacement: Optional[np.ndarray] = None,
     print_output: bool = False,
 ) -> Optional[Tuple[str, str]]:
     """
     Save symmetric displaced structures (TS ± scale) as XYZ files.
 
     Uses Python-style negative index wrapping (e.g., frame -1 is the last frame).
+    For scale <= max_level, coordinates are taken directly from trajectory frames.
+    For scale > max_level, mode_displacement (the Hessian normal mode vector) is used
+    to extrapolate: ts_coords ± scale * 0.2 * mode_displacement.
 
     Args:
         frames: List of frame dicts with 'symbols' and 'positions' keys
         ts_frame: Index of TS frame
         output_prefix: Prefix for output files
-        scale: Displacement scale (1-4, corresponding to amplitudes ~0.2-0.8)
-        max_level: Maximum allowed scale
+        scale: Displacement scale (1-max_level uses trajectory frames; >max_level extrapolates)
+        max_level: Trajectory frame limit (beyond this, extrapolation is used)
+        mode_displacement: Normal mode displacement vector from the Hessian (shape: n_atoms x 3).
+            Required for scale > max_level; ignored otherwise.
         print_output: Print status messages
 
     Returns
     -------
         (forward_path, reverse_path) if successful, else None
     """
-    n_frames = len(frames)
-
-    if not (1 <= scale <= max_level):
-        logger.warning(f"Invalid scale {scale} (must be 1-{max_level})")
+    if scale < 1:
+        logger.warning(f"Invalid scale {scale} (must be >= 1)")
         if print_output:
-            print(f"Invalid scale {scale} (must be 1-{max_level}).")
+            print(f"Invalid scale {scale} (must be >= 1).")
         return None
 
-    # Calculate indices with wrapping support
-    minus_idx = ts_frame - scale
-    plus_idx = ts_frame + scale
+    n_frames = len(frames)
 
     # Normalize indices (allow Python-style negative indexing)
     def normalize_index(idx: int) -> Optional[int]:
@@ -260,21 +262,48 @@ def normalize_index(idx: int) -> Optional[int]:
             return idx % n_frames
         return None
 
-    norm_minus = normalize_index(minus_idx)
-    norm_plus = normalize_index(plus_idx)
-
-    if norm_minus is None or norm_plus is None:
-        logger.warning(f"Scale {scale} out of range for TS {ts_frame} (total {n_frames} frames)")
-        if print_output:
-            print(f"Scale {scale} out of range for TS {ts_frame} (total {n_frames} frames).")
+    if scale <= max_level:
+        # Use trajectory frames directly
+        norm_minus = normalize_index(ts_frame - scale)
+        norm_plus = normalize_index(ts_frame + scale)
+
+        if norm_minus is None or norm_plus is None:
+            logger.warning(f"Scale {scale} out of range for TS {ts_frame} (total {n_frames} frames)")
+            if print_output:
+                print(f"Scale {scale} out of range for TS {ts_frame} (total {n_frames} frames).")
+            return None
+
+        paths = write_displaced_structures(frames, prefix=output_prefix, indices=[norm_minus, norm_plus])
+        if len(paths) == 2:
+            logger.info(f"Saved displaced pair (±{scale}): {os.path.basename(paths[0])}, {os.path.basename(paths[1])}")
+            if print_output:
+                print(f"Saved displaced pair (±{scale}): {os.path.basename(paths[0])}, {os.path.basename(paths[1])}")
+            return (paths[0], paths[1])
         return None
-
-    paths = write_displaced_structures(frames, prefix=output_prefix, indices=[norm_minus, norm_plus])
-
-    if len(paths) == 2:
-        logger.info(f"Saved displaced pair (±{scale}): {os.path.basename(paths[0])}, {os.path.basename(paths[1])}")
+    else:
+        # Extrapolate beyond trajectory using the Hessian normal mode vector
+        if mode_displacement is None:
+            logger.warning(f"Scale {scale} exceeds trajectory range and no mode displacement vector available")
+            if print_output:
+                print(f"Scale {scale} exceeds trajectory range and no mode displacement vector available.")
+            return None
+
+        ts_coords = frames[ts_frame]["positions"]
+        symbols = frames[ts_frame]["symbols"]
+        amp = scale * 0.2
+        f_coords = ts_coords + amp * mode_displacement
+        r_coords = ts_coords - amp * mode_displacement
+
+        f_path = f"{output_prefix}_F.xyz"
+        r_path = f"{output_prefix}_R.xyz"
+        write_xyz(f_path, symbols, f_coords)
+        write_xyz(r_path, symbols, r_coords)
+
+        logger.info(
+            f"Saved displaced pair (±{scale}, extrapolated): {os.path.basename(f_path)}, {os.path.basename(r_path)}"
+        )
         if print_output:
-            print(f"Saved displaced pair (±{scale}): {os.path.basename(paths[0])}, {os.path.basename(paths[1])}")
-        return (paths[0], paths[1])
-
-    return None
+            print(
+                f"Saved displaced pair (±{scale}, extrapolated): {os.path.basename(f_path)}, {os.path.basename(r_path)}"
+            )
+        return (f_path, r_path)