Analytical gradients of the Born-Oppenheimer variational energy of a molecular system with respect to nuclear coordinates are derived and implemented to optimize the molecular geometry in a basis of Singer's explicitly correlated Gaussian functions. The 3N - 6 geometry variables are optimized simultaneously with the linear and non-linear variational parameters of the wave function. The method can be used for global minimization, transition structure prediction, and following reaction paths in problems that require very accurate wave functions. Test results on hydrogen clusters with two and three atoms are presented as an illustration of the method.