IRC: Gonzalez–Schlegel (GS) Method

Overview

The second-order method of Gonzalez–Schlegel (GS) has been implemented in MAPLE. It operates in mass-weighted Cartesian coordinates and uses a two-phase pivot-step plus constrained-optimization approach at each IRC point:

Pivot step — A steepest-descent step is taken along the negative gradient direction in mass-weighted coordinates to produce a trial point on a hypersphere centered at the current IRC point.
Constrained optimization — The trial point is refined by minimizing the energy on the surface of the hypersphere, using the perpendicular component of the gradient, until the corrector converges.

The Hessian is updated along the path using BFGS/Bofill update formulas, avoiding full Hessian recomputation at every step. Adaptive step-size control adjusts the arc length based on path curvature, and an optional parabolic fit detects when the path has descended into a minimum basin.

The GS method integrates both the forward and backward IRC paths, starting from the transition state.

Parameters

All GS parameters can be specified in the #irc(...) header line. The complete GSParams dataclass:

Parameter	Type	Default	Description
`target_mode`	int	1	Index of the imaginary vibrational mode to follow (1-based).
`step_length_bohr`	float	0.10	Step length along the IRC path in Bohr.
`max_steps`	int	50	Maximum number of IRC steps in each direction.
`max_micro_cycles`	int	20	Maximum micro-iterations per IRC macro-step (constrained optimization on hypersphere).
`micro_step_thresh`	float	1e-3	Convergence threshold for micro-iterations.
`hessian_recalc`	int	None	Recalculate the full Hessian every N steps. None means never recalculate.
`hessian_update`	str	`bofill`	Hessian update method: `"bfgs"` or `"bofill"`.
`f_max_th`	float	2e-3	Max force convergence threshold (Hartree/Bohr).
`f_rms_th`	float	5e-4	RMS force convergence threshold.
`print_each`	bool	True	Print step details for each IRC point.
`write_traj`	bool	True	Write trajectory XYZ files.

Example: SN2 reaction

In this example, we will perform the calculations for the F- + CH3F -> FCH3 + F- reaction involved in the original GS reference. The TS structure is derived via the P-RFO method. Note that the system carries a net charge of -1 in this case.

Input Example

#model=aimnet2
#irc(method=gs)
#device=gpu0

XYZ  -1  1  /path/to/transition_state.xyz

Output Visualization

After a brief wait for the calculation to complete, we can visualize the reaction path and the energy profile along the reaction coordinate.

GS IRC animation — **Trajectory:** Simplified view of the reaction path

GS IRC graph — **Graph:** Energy trend along the IRC path

Ref Structure

Transition State

6
Iteration 26  Energy = -239.7111960742
H   0.0064725364  1.0714480485 -0.8400337833
H   0.9377518596 -0.5437668524 -0.8426007651
H  -0.9267066183 -0.5426693610 -0.8429760387
C   0.0058392639 -0.0049960268 -0.8418701156
F   0.0054541606 -0.0082444202  1.0624162248
F   0.0062246908 -0.0017475723 -2.7461566766