We have applied the electrostatically embedded many-body (EE-MB) method truncated at the two-body level (also called the pairwise additive EE-MB method or the EE-PA approximation) and the three-body level (called EE-3B) to calculate the gradient of the potential energy for a simulation box containing 64 water molecules. We employed the B3LYP density functional with the 6-31+G-(d,p) basis set for this test case. We found that the EE-PA method is able to reproduce the magnitude of the gradient from a B3LYP/6-31+G(d,p) calculation on the entire system to within 1.0% with a 1.3% error for the maximum component of the gradient. Furthermore, the EE-3B method is able to reproduce the magnitude of the gradient to within 0.1% with a 0.2% error for the maximum component of the gradient. The good performance of the EE-MB methods for calculating forces and the highly parallel nature of these methods make them well suited for use in molecular dynamics simulations. Furthermore, since the methods can be used for efficient and accurate calculations of forces with any level of electronic structure theory that has analytic gradients and with any electronic structure package that allows for the presence of a field of point charges, these methods can readily be used with a wide variety of density functional theory and wave function theory methods.