The analytic gradient is derived for the frozen domain formulation of the fragment molecular orbital (FMO) method combined with the polarizable continuum model. The accuracy is tested in comparison to full FMO calculations for a representative set of systems in terms of the gradient accuracy, protein-ligand binding energies, and optimized structures. The frozen domain method reproduced geometries optimized with full FMO within 0.03-0.09 angstrom in terms of reduced mean square deviations, whereas a single-point gradient calculation is accelerated by the factor of 38 (Trp-cage protein in explicit solvent, PDB: 1L2Y) and 12 (crambin, PDB: 1CRN). The method is applied to a geometry optimization of the K-Ras protein-ligand complex (4Q03) using two domain definitions, and the optimized structures are consistent with experiment. Pair interaction analysis is used to identify residues important in binding the ligand.