calculate thermodymamics

This python script calculate ZPE/TS/H using frequency info in OUTCAR

from monty.re import regrep

def get_frequency(outcar):
    pattern = {"frequency": "^\s+[\d]+\s+f+\s+=\s+([\d\-\.]+)\sTHz+\s+([\d\-\.]+)\s2PiTHz+\s+([\d\-\.]+)\s+cm-1+\s+([\d\-\.]+)+\smeV"}
    d = regrep(outcar, pattern)
    # float string and rename index
    data = []
    for index, d in enumerate(d['frequency']):
        d[0] = [float(i) for i in d[0]]
        d[1] = '{} f'.format(index + 1)
        data.append(d)
    return data



dd = get_frequency('OUTCAR')
meV_data = [d[-1] for d, index in dd]
l=len(meV_data)
Kb = 8.6173324E-2  #Boltzmann constant unit is meV/K
T  = 298       #Absolute temperature unit is K
R  = 8.3144598     #Gas constant unit is J/(k.mol)
Tran = 1.0364E-5   #J/mol unit to eV


D = []
S = []
C = []
import numpy as np
for i in xrange(0,l):
    d =meV_data[i]/(Kb*T)
    D.append(d)
    s = (d/(np.exp(d)-1))-np.log1p(-np.exp(-d))
    S.append(s)
    c = meV_data[i]/(Kb*(np.exp(d)-1))
    C.append(c)
entropy= sum(S)
TS= R*T*entropy*Tran
ZPE= 0.5*sum(meV_data)*1E-3
U = R*Tran*sum(C)
print "ZPE =", ZPE, "eV"
print "TS =", TS, "eV"
print  "U =", U, "eV"