In order to overcome the difficulty of optimizing molecular geometry using quantum Monte Carlo methods, we introduce various approximations to the exact force expectation value. We follow Pulay's suggestion [Mol. Phys. 17, 153 (1969)] to correct the Hellmann-Feynman estimator by introducing the contributions due to the changes in the wave function with respect to the nuclear positions. When used in conjunction with energy-optimized explicitly correlated trial wave functions for H-2 and LiH, these approximations appear to yield accurate forces using both the variational and diffusion Monte Carlo methods. Also, the accuracy of the second-order estimate of the Hellmann-Feynman force estimator was investigated employing our energy-optimized trial wave functions, and an erratic behavior was uncovered for some of the studied bond lengths. The additional computational cost required to compute the corrections to the Hellmann-Feynman estimator was found to be only a small fraction of the cost for a simple mean energy calculation. The same approach could be exploited also in computing the derivative of other energy-dependent quantum-mechanical observables. (C) 2003 American Institute of Physics.