Gibbs/gibbs_gaussian.py

#!/usr/bin/env python3
import math
import os
import glob
###########################################
version="0.8g"
print ("Gibbs",version)
###########################################
Lichtgeschwindigkeit=299792458
Planck=6.626176E-034
Boltzmann=1.380662E-023
R=8.31441
Bav=1E-043
Grenzwellen=100
Grenze=Grenzwellen*Lichtgeschwindigkeit*100
CalInJoule=4.184
AuInKcal=627.506
wiederholen=True

###########################################
def is_number(s):
    try:
        float(s)
        return True
    except ValueError:
        return False
###########################################
def allindices(string2, sub, listindex=[], offset=0):
    i = string2.find(sub, offset)
    while i >= 0:
        listindex.append(i)
        i = string2.find(sub, i + 1)
    return listindex
###########################################

for dateiname in glob.glob('*.log'):

	class Daten:
		Wellenzahl=[]
		Frequenz=[]
		VibTemp=[]
		Sv=[]
		Inertia=[]
		AvInertia=[]
		Sr=[]
		Gewichtung=[]
		Entropie=[]
	pass

	freq1=[]
	freq2=[]
	translation=[]
	rotation=[]
	Gesamtentropie=0
	GaussianEntropie=0

	Temp=298.15
    
	print("Log-File",dateiname,"wird verwendet...")
	ausgabedatei=dateiname+".gibbs"
	ausgabe = open(ausgabedatei,"w")
	print("Ausgabe",ausgabedatei,"wird verwendet...")
	print("Lese Datei",dateiname, "...")
	datei=open(dateiname)
	string = datei.read()
		
	warten=input("Weiter?")

	#E
	position=string.find("MP2=")
	ausschneiden=string[position+4:position+16]
	test=ausschneiden.find("\n")
	if not test==-1:
	    print("MP2-Energie ist dösig angezeigt...")
	    save1=string[position+4:position+4+test]
	    save2=string[position+4+test+2:position+18]
	    Energy=float(save1+save2)
	else:    
	    Energy=float(ausschneiden)
	#H
	position=string.find("Thermal correction to Enthalpy=")
	EnthalpyCorrection=float(string[position+49:position+59])
	Enthalpy=Energy+EnthalpyCorrection
	#S
	position=string.find("Total                  ")
	EntropyGaussian_cal=float(string[position+61:position+70])    
	EntropyGaussian=EntropyGaussian_cal/AuInKcal/1000
	#G
	position=string.find("Sum of electronic and thermal Free Energies=")
	FreeEnergy=float(string[position+52:position+66])
	#Frequenzen
	#for line in string:
	#    if "Frequencies" in line:
	#    freq1.append(line)
	#for i in range(0,len(freq1)):
	#    freq2+=freq1[i].split(" ")
	#for i in range(0,len(freq2)):
	#    if is_number(freq2[i]):
	#	ausgabe.write(freq2[i])
	#	ausgabe.write("\n")
	#ausgabe.close()
	#for line in open(ausgabedatei):
	#    if is_number(line):
	#	if float(line)>0:
	#	    Daten.Wellenzahl.append(line)
	#	else:
	#	    print("Imaginäre Frequenz", float(line), "wird ignoriert...")
	#print("Anzahl der Schwingungen:",len(Daten.Wellenzahl))

	#Array initialisieren
	Daten.Frequenz=[0]*len(Daten.Wellenzahl)
	Daten.VibTemp=[0]*len(Daten.Wellenzahl)
	Daten.Sv=[0]*len(Daten.Wellenzahl)
	Daten.Inertia=[0]*len(Daten.Wellenzahl)
	Daten.AvInertia=[0]*len(Daten.Wellenzahl)
	Daten.Sr=[0]*len(Daten.Wellenzahl)
	Daten.Gewichtung=[0]*len(Daten.Wellenzahl)
	Daten.Entropie=[0]*len(Daten.Wellenzahl)

	for i in range(0,len(Daten.Wellenzahl)):
	    tmp=Daten.Wellenzahl[i]
	    tmp2=float(tmp.replace("\n","0"))
	    Daten.Wellenzahl[i]=tmp2
	    tmp3=tmp2*Lichtgeschwindigkeit*100
	    Daten.Frequenz[i]=tmp3
	    tmp4=tmp3*Planck/Boltzmann
	    Daten.VibTemp[i]=tmp4
	    tmp5=R*((tmp4/Temp)/((math.exp(tmp4/Temp))-1)-math.log(1-(math.exp((-1)*tmp4/Temp))))
	    Daten.Sv[i]=tmp5
	    tmp6=Planck/(8*math.pi**2*tmp3)
	    Daten.Inertia[i]=tmp6
	    tmp7=(tmp6*Bav)/(tmp6+Bav)
	    Daten.AvInertia[i]=tmp7
	    tmp8=R*(0.5+math.log(((8*math.pi**3*tmp7*Boltzmann*Temp)/(Planck**2))**(0.5)))
	    Daten.Sr[i]=tmp8
	    tmp9=1/(1+(Grenze/tmp3)**4)
	    Daten.Gewichtung[i]=tmp9
	    tmp10=tmp9*tmp5+((1-tmp9)*tmp8)
	    Daten.Entropie[i]=tmp10

	#Berechnung
	Sv=sum(Daten.Sv)/CalInJoule/AuInKcal/1000
	Sr=sum(Daten.Entropie)/CalInJoule/AuInKcal/1000
	S=EntropyGaussian-Sv+Sr
	dG=EnthalpyCorrection-Temp*S+Energy
	#dG=Enthalpy-Temp*S
	#dG=FreeEnergy-EntropyGaussian+S

	#Ausgabe
	schreiben="Wellenzahl [1/cm]; Frequenz [1/s]; Vibrational Temperature [K]; Entropie Sv [J /mol*K]; Moment of Inertia; Average Moment of Inertia; Entropie Sr [J /mol*K]; Gewichtung; Entropie S [J /mol*K]\n"
	for i in range(0,len(Daten.Wellenzahl)):
	    schreiben+=str(Daten.Wellenzahl[i])+";"+str(Daten.Frequenz[i])+";"+str(Daten.VibTemp[i])+";"+str(Daten.Sv[i])+";"+str(Daten.Inertia[i])+";"+str(Daten.AvInertia[i])+";"+str(Daten.Sr[i])+";"+str(Daten.Gewichtung[i])+";"+str(Daten.Entropie[i])+"\n"
	schreiben+="#############################################################################################################\n"
	schreiben+="E(Gaussian) = "+str(Energy)+" a.u.; "+str(Energy*AuInKcal)+" kcal/mol; "+str(Energy*AuInKcal*CalInJoule)+" kJ/mol\n"
	schreiben+="H(Gaussian) = "+str(Enthalpy)+" a.u.; "+str(Enthalpy*AuInKcal)+" kcal/mol; "+str(Enthalpy*AuInKcal*CalInJoule)+" kJ/mol\n"
	schreiben+="S(Gaussian) = "+str(EntropyGaussian)+" a.u.; "+str(EntropyGaussian*AuInKcal)+" cal/mol*K; "+str(EntropyGaussian*AuInKcal*CalInJoule)+" J/mol*K\n"
	schreiben+="G(Gaussian) = "+str(FreeEnergy)+" a.u.; "+str(FreeEnergy*AuInKcal)+" kcal/mol; "+str(FreeEnergy*AuInKcal*CalInJoule)+" kJ/mol\n"
	schreiben+="S(korrigiert) = "+str(S)+" a.u.; "+str(S*AuInKcal)+" cal/mol*K; "+str(S*CalInJoule)+" J/mol*K\n"
	schreiben+="G(korrigiert) = "+str(dG)+" a.u.; "+str(dG*AuInKcal)+" kcal/mol; "+str(dG*AuInKcal*CalInJoule)+" kJ/mol\n"
	schreiben+="#############################################################################################################\n"
	schreiben+="G(Gaussian)-G(korrigiert) = "+str(FreeEnergy-dG)+" a.u.; "+str(FreeEnergy*AuInKcal-dG*AuInKcal)+" kcal/mol; "+str(FreeEnergy*AuInKcal*CalInJoule-dG*AuInKcal*CalInJoule)+" kJ/mol\n"
	schreiben+="#############################################################################################################\n"
	ausgabe = open(ausgabedatei,"w")
	ausgabe.write(schreiben)
	ausgabe.close()

	#Anzeige
	print(schreiben)
        
print("Copyright 2015 SAW")