  We present a novel approach to efficiently implement thermal expansion in the
quasi-harmonic approximation (QHA) for both isotropic and more importantly,
anisotropic expansion. In this approach, we rapidly determine a crystal's
equilibrium volume and shape at a given temperature by integrating along the
gradient of expansion from zero Kelvin up to the desired temperature. We
compare our approach to previous isotropic methods that rely on a brute-force
grid search to determine the free energy minimum, which is infeasible to carry
out for anisotropic expansion, as well as quasi-anisotropic approaches that
take into account the contributions to anisotropic expansion from the lattice
energy. We compare these methods for experimentally known polymorphs of
piracetam and resorcinol and show that both isotropic methods agree to within
error up to 300 K. Using the Gr\"{u}neisen parameter causes up to 0.04 kcal/mol
deviation in the Gibbs free energy, but for polymorph free energy differences
there is a cancellation in error with all isotropic methods within 0.025
kcal/mol at 300 K.
  Anisotropic expansion allows the crystals to relax into lattice geometries
0.01-0.23 kcal/mol lower in energy at 300 K relative to isotropic expansion.
For polymorph free energy differences all QHA methods produced results within
0.02 kcal/mol of each other for resorcinol and 0.12 kcal/mol for piracetam, the
two molecules tested here, demonstrating a cancellation of error for isotropic
methods.
  We also find that when expanding in more than a single volume variable, there
is a non-negligible rate of failure of the basic approximations of QHA.
Specifically, while expanding into new harmonic modes as the box vectors are
increased, the system often falls into alternate, structurally distinct
harmonic modes unrelated by continuous deformation from the original harmonic
mode.
