The energy gradient of the second-order Moller-Plesset perturbation theory (MP2) has been implemented in conjunction with the fragment molecular orbital-based molecular dynamics (FMO-MD) method including up to three-body correction (FMO3). A hybrid integral-direct approach of both atomic and molecular orbital indices was utilized with parallelism for the gradient calculations. A droplet model consisting of 64 water molecules was then simulated with the 6-31G* basis set. The first peak position of O-O radial distribution function was evaluated to be 2.75 angstrom at the FMO3-MP2 level, whereas the corresponding Hartree-Fock (FMO3-HF) value was 2.89 angstrom. Comparison with an X-ray value of 2.73 angstrom indicated better reliability of the MP2 gradient for FMO-MD. (C) 2011 Elsevier B. V. All rights reserved.