Script to calculate beta from phonon output

This combines a simple .phonon parser with
code to evaluate equation (2) of
Blanchard et al. 2009 (GCA 73:6565-6578,
doi:10.1016/j.gca.2009.07.034)
Handle different q-pt weights

If using symmetry to reduce the Brillouin zone
sampling, e.q. 2 of Blanchard should deal with
this. This change (1) extends the parser to
report q-point weights (which means we know the
number of q-points as well as the number of
frequencies) and (2) changes the equation for
Beta to use this. We change from this:

\beta = [\Pi_3N \Pi_{q} ...]^{1/(N_qN)}

to this:

\beta = [\Pi_{q} (\Pi_3N ...)^{wq}]^{1/N}

where wq is the weight on each q-pt. If we
don't use symmetry wq is always 1/N_q, but
not if we do.