geomeTRIC Geometry Optimization

GeomeTRIC uses a trust radius quasi-Newton method (similar to Levenberg-Marquardt) to determine optimization steps for energy minimization.

1. Internal Coordinate (IC) Step: A proposed step $\boldsymbol{\delta}^{(q)}$ is calculated using the IC Hessian $\mathbf{H}^{(q)}$, the IC gradient $\mathbf{g}^{(q)}$, and a regularization parameter $\lambda$: $$\boldsymbol{\delta}^{(q)} = -(\mathbf{H}^{(q)} + \lambda I)^{-1} \mathbf{g}^{(q)}$$
   • As $\lambda \rightarrow 0$, it recovers the full quasi-Newton step.
   • As $\lambda \rightarrow \infty$, it becomes an infinitely small steepest descent step.
2. Gradient Conversion: The IC gradient is derived from the Cartesian gradient $\mathbf{g}^{(x)}$ via $\mathbf{g}^{(q)} = \mathbf{G}^{\dagger} \mathbf{B} \mathbf{g}^{(x)}$, where $\mathbf{G}$ is inverted using singular value decomposition (SVD), ignoring singular values smaller than $10^{-6}$.
3. Hessian Update: The approximate IC Hessian is updated using the BFGS method at the end of each optimization step.

geomeTRIC uses the Restricted-step Partitioned Rational Function Optimization (RS-P-RFO) method for TS optimizations.

How it works:

1. Approximation: The Potential Energy Surface (PES) is approximated using a [2/2] Padé approximant.
2. Partitioning: The optimization space is partitioned into two subspaces based on the Hessian eigenvalues:
   • Maximization Subspace: Corresponds to the lowest eigenvalue/eigenvector pair (the 'transition vector' or tv). The algorithm seeks the highest eigenvalue in this subspace to move along the reaction path.
   • Minimization Subspace: Corresponds to all other eigenvalues (ot). The algorithm seeks to minimize the energy in these directions.
3. Restriction (RS): The 'Restricted-step' aspect ensures the proposed step stays within the internal coordinate (IC) trust radius. If an unrestricted step ($\alpha = 1$) exceeds the trust radius $R$, the parameter $\alpha$ is iteratively optimized until the step size matches the radius within a tolerance of 0.001.
4. Hessian Updates: Unlike energy minimization (which uses BFGS), TS optimization uses a combination of Murtagh-Sargent (symmetric rank-1) and Powell symmetric Broyden (PSB) updates. This specific combination is used because it does not require the Hessian to remain positive-definite, which is necessary for finding saddle points.

Because internal coordinates (ICs) are nonlinear functions of Cartesian coordinates, an IC step $\boldsymbol{\delta}^{(q)}$ must be iteratively converted into a Cartesian step $\boldsymbol{\delta}^{(x)}$.

The Iterative Process:

1. An initial approximate Cartesian step $\boldsymbol{\delta}^{(x)}_1$ is calculated from the current coordinates $\mathbf{x}_0$ and the IC step $\boldsymbol{\delta}^{(q)}$.
2. The new coordinates are updated: $\mathbf{x}_1 = \mathbf{x}_0 + \boldsymbol{\delta}^{(x)}_1$.
3. Because the step is approximate, a 'gap' is calculated between the target IC step and the actual change in ICs: $\boldsymbol{\delta}^{(q)}_1 = \boldsymbol{\delta}^{(q)} - (\mathbf{q}(\mathbf{x}_1) - \mathbf{q}(\mathbf{x}_0))$.
4. This gap is treated as the new target IC step, and the process repeats using $\mathbf{x}_1$ as the input.

Convergence: The iterations continue until the difference between the target IC step and the actual IC change is less than $10^{-6}$. The final Cartesian step is $\boldsymbol{\delta}^{(x)} = \mathbf{x}_n - \mathbf{x}_0$.

Note: Periodic coordinates are handled specially to account for $2\pi$ jumps during IC difference calculations.

Transition State (TS) optimizations in geomeTRIC use the TRIC (Transition-state Internal Coordinate) system. This system relies on dividing the molecular system into one or more fragments.

To automate this, geomeTRIC uses the Mayer bond order matrix (if the selected engine supports it) to identify fragments. Two atoms are considered bonded if their bond order exceeds a specific threshold.

Key behaviors:

• Thresholding: The threshold is controlled by the --bothre flag. If one or both atoms are transition metals, the threshold is halved.
• Fragment Creation: If the process results in only one fragment, the threshold is automatically increased until at least two fragments (each containing >1 atom) are found, up to a maximum of 0.75.
• Isolated Atoms: Any atoms that remain isolated after fragment identification are automatically connected to their nearest neighbor before the internal coordinate system is constructed.

A Transition State (TS) is an energy maximum along a reaction pathway. In terms of the Potential Energy Surface (PES), a TS is characterized by exactly one negative eigenvalue in the Hessian matrix, representing a vibrational mode where the PES is concave down.

Key differences from energy minimization:

• Goal: While minimization seeks an energy minimum, TS optimization seeks to maximize the energy along the lowest eigenvalue mode of the Hessian.
• Initial Guess: TS optimization is highly sensitive to the initial structure. Performance is best when the starting structure already possesses a negative eigenvalue with an eigenvector pointing along the expected reaction direction. Common ways to generate these guesses include coordinate driving (scanning), internal coordinate interpolation, or reaction path methods like Nudged Elastic Band (NEB).
• Hessian Requirement: TS optimization requires a Hessian calculation at the start. geomeTRIC uses central finite difference of the gradient, which requires $6 \times N_{\text{atoms}}$ gradient calculations.

Internal coordinates (ICs) are generated based on atomic connectivity and the chosen coordinate system type.

Primitive Internal Coordinates

Primitive ICs are a list of individual distances, angles, and dihedral angles.

• Distances: Added for each connected atom pair.
• Angles: Added for every triplet of connected atoms.
• Dihedrals: Added for every "linear" connected quartet.

Each IC provides value(), derivative(), and second_derivative() methods. The derivatives of all primitive ICs can be stacked to form the Wilson B-matrix ($dq_i / dx_j$). The product $G = BB^T$ forms the G-matrix (an overlap matrix).

Delocalization

To improve numerical stability and reduce cost, geomeTRIC can perform a delocalization procedure. This involves diagonalizing the G-matrix and using the resulting eigenvectors (representing non-redundant linear combinations of primitive ICs) as the new coordinate system. This is used in --tric, --dlc, and --hdlc modes.