  We report the analytical nuclear gradient theory for complete active space
second-order perturbation theory (CASPT2) with imaginary shift, which is
commonly used to avoid divergence of the perturbation expression. Our
formulation is based on the Lagrangian approach and is an extension of the
algorithm for CASPT2 nuclear gradients with real shift. The working equations
are derived and implemented into an efficient parallel program. Numerical
examples are presented for the ground- and excited-state geometries and conical
intersections of a green fluorescent protein model chromophore, $p$-HBDI$^-$.
We also report timing benchmarks with adenine, $p$-HBDI$^-$, and iron
porphyrin. It is demonstrated that the energies and geometries obtained with
the imaginary shift improve accuracy at a minor additional cost which is mainly
associated with evaluating the effective density matrix elements for the
imaginary shift term.
