Analytic Nuclear Gradients

The RHF analytic gradient is:

$$\frac{\partial E}{\partial X_A}=\sum _{\mu \nu }{P}_{\mu \nu }\frac{\partial h_{\mu \nu }}{\partial X_A}+\frac{1}{2}\sum _{\mu \nu \lambda \sigma }{P}_{\mu \nu }{P}_{\lambda \sigma }\frac{\partial (\mu \nu |\lambda \sigma )}{\partial X_A}-\sum _{\mu \nu }{W}_{\mu \nu }\frac{\partial S_{\mu \nu }}{\partial X_A}+\frac{\partial E_{\text{nuc}}}{\partial X_A}$$

where $W=-{\mathbf{C}}_{\text{occ}}{\mathbf{\epsilon }}_{\text{occ}}{\mathbf{C}}_{\text{occ}}^{†}$ is the energy-weighted density matrix.

The UHF gradient has separate α and β components for the density and energy-weighted density matrices.

Gradient Components

Symmetry Acceleration

Analytic gradients use symmetry-accelerated ERI derivatives (skips symmetry-equivalent shell quartets). The detected point group is used to skip symmetry-equivalent ERI shell quartets during gradient computations, matching the GRAD2E SymShl approach.

This provides near-proportional speed-ups based on symmetry order (e.g., 2× for C₂, 4× for D₂ₕ) for quadratically-scaling steps.

Implementation

Gradient integrals are evaluated using the Obara–Saika recurrence for derivatives. The derivative of a Gaussian integral with respect to nuclear position is computed by differentiating the recurrence relations.

The gradient is used by:

• BFGS optimizer (Cartesian)
• Berny optimizer (internal coordinates)
• IRC path following
• NEB pathway optimization
• Geometry-space DIIS (GDIIS/GEDIIS)