Intrinsic reaction coordinate

Basics

The Intrinsic Reaction Coordinate (IRC) method aims to trace a minimum energy pathway on a potential energy surface (PES) starting with an optimized transition state (TS) structure. This optimized TS structure sits at a first-order saddle point on the PES where the structure has only one imaginary vibrational frequency mode. To begin, geomeTRIC calculates the Hessian and vibrational frequencies to confirm this imaginary mode. The calculation of Hessian can be carried in parallel using BigChem or Work Queue.

Once the Hessian calculation is completed, the first step is taken in the positive direction of the corresponding eigenvector of the imaginary mode using mass-weighted internal coordinates. The Hessian is updated using the BFGS method and the succeeding steps are guided by the instantaneous acceleration vectors. The IRC method continues taking the steps until it meets the same convergence criteria as geometry optimization.

geomeTRIC repeats the same procedure for the negative direction after the positive direction IRC converges. Once both directions have converged, one of the paths is reversed to concatenated with the other, using the TS structure as the connecting point.

Usage

The IRC method can be used by passing --irc yes on the command line with geometric-optimize. The input geometry should be an optimize TS structure. geomeTRIC adjusts the IRC step size based on the step quality, and the maximum step size will be set equal to the initial trust radius --trust [0.3]. The direction of IRC can be specified using --irc_direction [both]. The forward direction follows the imaginary mode from the frequency analysis as is, while the backward direction follows its opposite. Once convergence is achieved, geomeTRIC will generate an output xyz file containing the IRC pathway.

Example

See examples/1-simple-examples/hcn_hnc_irc directory.

Theory

The IRC method uses a large portion of the same code as the optimization algorithm. The two main differences are in obtaining steps and adjusting the trust radius.

Obtaining the IRC step

geomeTRIC implemented Gonzalez and Schlegel’s mass-weighted internal coordinates IRC method. The mass-weighted Wilson B-matrix (\(\mathbf{B}\)) has elements of \(dq_i / (dx_j \sqrt{m_j})\), where \(m_j\) is the atomic mass, and the G-matrix is calculated as \(\mathbf{G} = \mathbf{B}\mathbf{B}^T\) (See internal coordinate setup).

At the beginning of the first iteration, the mass-weighted step size (\(\mathbf{s}\)) is calculated from the trust radius (\(R_{\mathrm{trust}}\)) as follows:

\[\begin{split}\begin{aligned} & A = \sqrt{\sum_{i=1}^{N_{\mathrm{atoms}}} (\Delta x_i^2 + \Delta y_i^2 + \Delta z_i^2) \times m_i} \\ & \mathrm{s} = R_{\mathrm{trust}} \times A \end{aligned}\end{split}\]

where \(\Delta x\), \(\Delta y\), and \(\Delta z\) are the normalized Cartesian displacements along the imaginary mode on the saddle-point. The conversion factor \(A\) is used in every iteration to convert the trust radius.

Each IRC step starts with taking a half-step towards a pivot point following the instantaneous accelerations (\(-\mathbf{G} \cdot \mathbf{g}\)). The pivot point (\(\mathbf{\mathrm{q}}^*\)) is obtained by:

\[\mathbf{q}^* = \mathbf{q} - \frac{\mathrm{s}}{2} \cdot \frac{\mathbf{G} \cdot \mathbf{g}}{(\mathbf{g}^T \cdot \mathbf{G} \cdot \mathbf{g})^{1/2}}\]

where \(\mathbf{g}\) is the gradients of internal coordinates \(\mathbf{q}\). To reach the next point on the reaction path, another half-step needs to be taken from the pivot point along a vector that is 1) parallel to the acceleration vector of the next point and 2) has a scalar value equal to half of the mass-weighted step size (\(\mathrm{s}/2\)).

First, a guessed point (\(\mathbf{q}_1^\prime\)) is obtained by taking the same step as the initial half-step starting from the pivot point. The guessed point is guided to the next guessed point (\(\mathbf{q}_2^\prime\)) until the vector pointing from the pivot point to the guessed point satisfies the two conditions.

If we define the following two vectors:

(1)\[\begin{split}\begin{aligned} & \mathbf{p} = \mathbf{q}_n^\prime - \mathbf{q}^* \\ & \Delta\mathbf{q} = \mathbf{q}_{n+1}^\prime - \mathbf{q}_n^\prime \end{aligned}\end{split}\]

\(\Delta\mathbf{q}\) can be updated while keeping the scaler value of \(\mathbf{p}\) equal to \(\mathrm{s}/2\) until \(\mathbf{q}_n^\prime\) reaches the point where \(\mathbf{p} + \Delta\mathbf{q}\) is parallel to the acceleration vectors at point \(\mathbf{q}_{n+1}^\prime\). The vectors, Hessian, and gradients in mass-weighted internal coordinates can be expressed as

(2)\[\begin{split}\begin{aligned} & \Delta\mathbf{q}_\mathrm{M} = \mathbf{G}^{-1/2} \Delta\mathbf{q}\\ & \mathbf{p}_\mathrm{M} = \mathbf{G}^{-1/2} \mathbf{p}\\ & \mathbf{g}_\mathrm{M} = \mathbf{G}^{1/2} \mathbf{g}^{\prime}\\ & \mathbf{H}_\mathrm{M} = \mathbf{G}^{1/2} \mathbf{H} \mathbf{G}^{1/2}\\ \end{aligned}\end{split}\]

where \(\mathbf{g}^{\prime}\) represents the estimated gradients at the point \(\mathbf{q}_n^\prime\), using a quadratic expansion. \(\mathbf{G}\) is calculated at \(\mathbf{q}_n^\prime\) as well.

The step size constraint can be expressed as:

(3)\[(\mathbf{p}_\mathrm{M} + \Delta\mathbf{q}_\mathrm{M})^{T}(\mathbf{p}_\mathrm{M} + \Delta\mathbf{q}_\mathrm{M}) = (\frac{\mathrm{s}}{2})^2\]

The other condition is satisfied at the convergence point (the next point), when the following equation holds true:

(4)\[(\mathbf{g}_\mathrm{M} - \lambda \mathbf{p}_\mathrm{M}) + (\mathbf{H}_\mathrm{M} - \lambda \mathbf{I})\Delta\mathbf{q}_\mathrm{M} = 0\]

where \(\lambda\) is the Lagrangian multiplier and \(\mathbf{I}\) is the identity matrix.

Eq. (4) can be rearranged as follows:

(5)\[\Delta\mathbf{q}_\mathrm{M} = -(\mathbf{H}_\mathrm{M} - \lambda \mathbf{I})^{-1}(\mathbf{g}_\mathrm{M} - \lambda \mathbf{p}_\mathrm{M})\]

\(\lambda\) is calculated iteratively after introducing Eq. (5) to (3). \(\Delta\mathbf{q}_\mathrm{M}\) is then used to move \(\mathbf{q}_n^\prime\) to \(\mathbf{q}_{n+1}^\prime\) and new Eq. (1) and Eq. (2) are defined to calculate the next \(\Delta\mathbf{q}_\mathrm{M}\). This process repeats until the norm of \(\Delta\mathbf{q}\) falls below 1e-6. It then takes the rest of the half-step along \(\mathbf{p} + \Delta\mathbf{q}\) from the pivot point, which completes an iteration.

Trust radius adjustment

The step quality (\(Q\)) is calculated in the same way as the energy minimization step quality. The trust radius is adjusted as follows:

• \(Q \geq 0.75\) : “Good” step, trust radius is increased by a factor of \(\sqrt{2}\), but not greater than the maximum.

• \(0.75 > Q \geq 0.50\) : “Okay” step, trust radius is unchanged.

• \(Q < 0.50\) : Step is rejected, trust radius is decreased by setting it to \(0.5 \times \mathrm{min}(R_{\mathrm{trust}}, \mathrm{RMSD})\), but not lower than the minimum