Note that the approach developed in this study cannot be applied to molecules that have degenerate active orbitals, for instance, benzene. In the course of determining κ act,γ,αβ tu , we divide non-vanishing contributions by the difference of orbital energies (with quasi-canonical orbitals), which can be singular (Eq. ( )). A similar difficulty has been found with analytic derivatives of the strongly contracted variant of NEVPT2 which is not invariant with respect to rotations among inactive and secondary orbitals. In both cases, the origin of the problem is the arbitrariness of the perturbation energy. Hence, if the target system has a high spatial symmetry or can have accidental degeneracies in the active orbital region, optimizing it with IPEA shift is not recommended.