We report the development of the theory and computer program for analytical nuclear energy gradients for (extended) multi-state complete active space perturbation theory (CASPT2) with full internal contraction. The vertical shifts are also considered in this work. This is an extension of the fully internally contracted CASPT2 nuclear gradient program, recently developed for a state-specific variant by us [MacLeod and Shiozaki, J. Chem. Phys. 142, 051103 (2015)]; in this extension, the so-called λ equation is solved to account for the variation of the multi-state CASPT2 energies with respect to the change in the amplitudes obtained in the preceding statespecific CASPT2 calculations, and the Z-vector equations are modified accordingly. The program is parallelized using the MPI3 remote memory access protocol that allows us to perform efficient one-sided communication.The optimized geometries of the ground and excited states of a copper corrole and benzophenone are presented as numerical examples. The code is publicly available under the GNU General Public License.