BOLZ  = 1.380658E-23;
HBAR  = 1.054589E-34;

UENERGY_IN_JOULE = 1.602192E-19;
ULENGTH_IN_METER = 1.E-10;
UMASS_IN_KG = 1.66053E-27;
UTIME_IN_SECOND = ULENGTH_IN_METER*sqrt(UMASS_IN_KG/UENERGY_IN_JOULE);
UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/pi;
UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12;

load freq.out;
a = size(freq);

fprintf(1,'Zero-pt temperature = %f\n', 0.5*HBAR*sum(freq)*UFREQ_IN_HZ/a(1)/BOLZ*2*pi);

