In order to predict the form of a stable crystalline polymorph, one compares the free energy of each form, and selects the one with the lowest free energy.  Although this ignores kinetic influences in crystallization, this procedure provides important information for polymorph prediction. However, current methods take only the potential-energy contribution to the free energy into account, and ignore the entropic contribution.  This is done because calculating the entropy is difficult to do rigorously; approximations using lattice dynamics are sometimes applied but the accuracy of these approaches is uncertain.  When there is little difference between the potential contributions to the energy, it is hard to distinguish which of the polymorphs is the most stable.    To calculate a free energy, given a reference system of known or easily calculated energy, we calculate the difference in free energy between the system of interest and the reference system, which in effect gives the free energy of the system of interest.  For this work, we use a harmonic reference system with spring constants given to match configurational correlations measured in the target system.    We consider two approaches to compute the free energy difference between the target and reference systems.  Direct perturbation is not effective, so we examine the performance of overlap sampling approaches, and Bennett's method in particular.  Second, we examine the accuracy of the Normal Mode Monte Carlo (NMMC) method, which, is an approximate treatment that assumes that normal mode coordinates are independent not only in the harmonic system, but also in the reference.  This technique provides much better sampling accuracy than direct or staged-perturbation  methods, but the approximations inherent in its formulation have not been well tested.   We study these approaches as applied to model molecular crystals for which the free energy has been determined by more computationally demanding methods.