The natural auxiliary function (NAF) approach is an approximation to decrease the size of the auxiliary basis set required for quantum chemical calculations utilizing the density fitting technique. It has been proven efficient to speed up various correlation models, such as secondorder Møller−Plesset (MP2) theory and coupled-cluster methods. Here, for the first time, we discuss the theory of analytic derivatives for correlation methods employing the NAF approximation on the example of MP2. A detailed algorithm for the gradient calculation with the NAF approximation is proposed in the framework of the method of Lagrange multipliers. To assess the effect of the NAF approximation on gradients and optimized geometric parameters, a series of benchmark calculations on small and medium-sized systems was performed. Our results demonstrate that, for MP2, sufficiently accurate gradients and geometries can be achieved with a moderate time reduction of 15−20% for both small and medium-sized molecules.