An analytical gradient theory for single-state N-electron valence state perturbation theory (NEVPT2), using both strongly contracted (SC) and partially contracted (PC) internal contraction schemes, is developed. We demonstrate the utility of the developed algorithm in the optimization of the single-state molecular geometry of acrolein, benzyne, benzene, the retinal chromophore PSB3, the GFP chromophore pHBI, and porphine, with the cc-pVTZ basis sets. The SC-NEVPT2 analytical gradients exhibit numerical instability due to the lack of invariance with respect to the rotations among the inactive orbitals. On the other hand, PC-NEVPT2 gives molecular geometries comparable to the second-order complete active space perturbation theory in any tested cases. We discuss possible future developments that will make the NEVPT2 gradient algorithm a powerful tool for optimizing molecular geometries and conducting molecular dynamics simulations of correlated systems.