#!/usr/bin/python3.6
from numpy import *

from math import *
import sys

'''
program vibpartitionClassic_eV.py
	- this program does the same as vibpartition_eV.py (and also the usage is identical) but
	it uses semi-classical instead of quantum mechanical expressions
'''

class Thermodynamics:
    def __init__(self, sigma, T, p, nfile='vib', mysmall=1.):
        AVOGADRO = 6.02214e23
        self.PLANCK = 6.62608e-34
        self.BOLTZMANN = 1.38065e-23
        self.CL = 2.997925e08
        self.eVTOJ = 1.60218e-19

        qvib, Uvib, Hvib, Svib, Gvib, zpe = self.part_vib(T, nfile, mysmall)
        # print 'qvib:',qvib

        Uvib /= self.eVTOJ
        Hvib /= self.eVTOJ
        Svib /= self.eVTOJ
        Gvib /= self.eVTOJ
        zpe /= self.eVTOJ

        print('')
        # print 'T [K]','Hvib [J/mol]','ZPE [J/mol]','Hvib-ZPE [J/mol]'
        # print T,Hvib,zpe,Hvib-zpe
        # print 'T [K]','Hvib [J/mol]','Svib [J/K/mol]','Gvib [J/mol]'
        # print T,Hvib,Svib,Hvib-T*Svib
        print('T [K]', 'ZPE [eV]', 'Hvib-ZPE [eV]', 'Hvib [eV]', 'Svib [eV/K]', 'Gvib [eV]')
        print(T, zpe, Hvib-zpe, Hvib, Svib, Gvib)

    def part_vib(self, T, nfile, mysmall):
        EXP = 2.7182818
        modes = open(nfile)
        vibs = []
        vib_count = 0
        vib_total = 0
        for line in modes.readlines():
            # vib_total+=1
            row = line.split()
            # read in wavenumbers in cm-1
            if len(row) == 11:
                if row[1] == "f" and row[2] == "=":
                    vib_total += 1
                    vm = float(row[3])
                    if vm > mysmall:
                        vib_count += 1
                        # convert to energy in J
                        vm = vm*self.PLANCK*self.CL*1e2
                        vibs.append(vm)
                elif row[1] == "f/i":
                    vib_total += 1
        modes.close()
        print("nDOF:", vib_total)
        print("nVIB:", vib_count)
        zpe = 0.
        # print "ZPE (J):", zpe/1.60218e-19

        qvib = 1.
        # Uvib=0.5*vib_count*self.BOLTZMANN*T
        # Hvib=0.5*vib_count*self.BOLTZMANN*T
        Uvib = vib_count*self.BOLTZMANN*T
        Hvib = vib_count*self.BOLTZMANN*T
        Svib = 0.
        Gvib = 0.
        for i in range(len(vibs)):
            part = EXP**(-vibs[i]/(self.BOLTZMANN*T))
            qvib = qvib*part
            # Gvib+=-self.BOLTZMANN*T*log(part)
            Gvib += log(self.BOLTZMANN*T/vibs[i])
        Gvib *= -self.BOLTZMANN*T
        Svib = (Hvib-Gvib)/T
        return qvib, Uvib, Hvib, Svib, Gvib, zpe


mysmall = 1.
T = float(sys.argv[1])
nfile = sys.argv[2]
print(nfile)
Thermodynamics(2, T, 1e5, nfile, mysmall)
