The equations for the response terms for the fragment molecular orbital (FMO) method interfaced with the density functional theory (DFT) gradient are derived and implemented. Compared to the previous FMO–DFT gradient, which lacks response terms, the FMO–DFT analytic gradient has improved accuracy for a variety of functionals, when compared to numerical gradients. The FMO–DFT gradient agrees with the fully ab initio DFT gradient in which no fragmentation is performed, while reducing the nonlinear scaling associated with standard DFT. Solving for the response terms requires the solution of the coupled perturbed Kohn–Sham (CPKS) equations, where the CPKS equations are solved through a decoupled Z-vector procedure called the self-consistent Z-vector method. FMO–DFT is a nonvariational method and the FMO–DFT gradient is unique compared to standard DFT gradients in that the FMO–DFT gradient requires terms from both DFT and time-dependent density functional theory (TDDFT) theories.