The initial implementation of (full) CASPT2 goes back to 1992, 8 but the complexity of this MRPT has hindered the development of its analytic derivative. The first analytic gradient implementation for CASPT2 was achieved in 2003, and that for GVVPT2 was developed in the same year by evaluating the derivative of wavefunction parameters explic-itly (solving coupled-perturbed equations). The former method employs a partial internal contraction scheme; this avoids the evaluation of higher-order reduced density matrices (RDMs), and therefore the computational demand may be high for large active spaces. After more than 10 years, analytic gradients for the fully internally contracted CASPT2 were implemented using an automatic code-generation technique. Despite the complexity of the formalism, development of analytic derivatives is ongoing, and they have been proposed for many MRPTs: CASPT2, GVVPT2, 25 NEVPT2, (X)MCQDPT2, and (SA-)DSRG-MRPT2. Some of the above implementations can be used for locating minimum energy conical intersections (MECIs) and performing MD simulations. However, it is well recognized that CASPT2 using the standard Fock operator tends to underestimate excitation energies. The origin of this underestimation has been attributed to an overestimation of open-shell (typically excited) states, for which the ionization potentialelectron affinity (IPEA) shift has been suggested as a remedy. Some studies in the literature have statistically proved that vertical excitation energies are indeed improved with IPEA shift. It has also been reported that geometrical parameters and molecular properties may be improved. In spite of its potential utility, no analytic derivatives using IPEA shift have yet been developed; this may be partially due to the requirement for introducing an empirical parameter into "ab initio" methodology rather than technical difficulties.