The first analytic derivative with respect to nuclear coordinates is derived for the fragment molecular orbital method used with electrostatic approximations. It is shown how response contributions due to the coupling between the electronic state of fragments and the embedding potential can be accurately and efficiently calculated. The accuracy of the analytic gradient is shown in comparison to the numerical gradient for a polypeptide and water clusters. For the largest system the overall speedup of about 4.9 is observed when approximations are used in FMO. (c) 2012 Elsevier B.V. All rights reserved.