When no IPEA shift is employed, the errors from the numerical gradient with both "Analytic0" and "Analytic" are of the order of 10 -7 hartree/bohr, which should be close to the limit of machine precision. This indicates that the original implementation gives accurate gradients, and the additional constraint conditions do not contribute to the computed gradients, because XMS-CASPT2 (no IPEA) is an invariant theory. Note that the implementation of XMS-CASPT2 in OpenMolcas has small non-invariance owing to the non-uniform first-order interacting space. On the other hand, when IPEA shift is employed, only the full analytic gradient ("Analytic" entry) reproduces the numerical gradient. Other approximate gradients ("Analytic0", "Analytic1", and "Analytic2" entries) can have errors of more than 3.0 × 10 -5 hartrees/bohr for relatively simple molecular and electronic structures, so they do not give analytic gradients. Moreover, considering that the error for S 1 is large, the error may be even larger for more complex electronic structures. The magnitude of the error of the "Analytic" gradients is similar to that in the case without IPEA shift.