A method is presented for efficient calculation of free energies for molecular crystals. The method is based on a generalisation of the local harmonic approximation. These calculations are rapid enough to allow optimisation of the free energy with respect to all atomic and crystallographic coordinates at finite temperatures for complex molecular crystals. The procedure has been illustrated for methane hydrate, which has about 150 atoms in the unit cell, and ice Ih. The resulting free energies of the fully optimised structures have been used to predict the conditions for three-phase equilibrium between methane hydrate, ice and methane gas, and shown to give results in excellent agreement with experiment.