The authors present a method based on a linear response theory that allows one to optimize the geometries of quantum mechanical/molecular mechanical (QM/MM) systems on the free energy surfaces. Two different forms of linear response free energy functionals are introduced, and electronic wave functions of the QM region, as well as the responses of electrostatic and Lennard-Jones potentials between QM and MM regions, are self-consistently determined. The covariant matrix relating the QM charge distribution to the MM response is evaluated by molecular dynamics (MD) simulation of the MM system. The free energy gradients with respect to the QM atomic coordinates are also calculated using the MD trajectory results. They apply the present method to calculate the free energy profiles of Menshutkin-type reaction of NH3 with CH3Cl and Claisen rearrangement of allyl vinyl ether in aqueous solution. For the Menshutkin reaction, the free energy profile calculated with the modified linear response free energy functional is in good agreement with that by the free energy perturbation calculations. They examine the nonequilibrium solvation effect on the transmission coefficient and the kinetic isotope effect for the Claisen rearrangement.