|
|
#!/usr/bin/env python3
|
|
|
import math
|
|
|
import os
|
|
|
import glob
|
|
|
###########################################
|
|
|
version="0.8vs"
|
|
|
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
|
|
|
###########################################
|
|
|
|
|
|
while wiederholen:
|
|
|
|
|
|
class Daten:
|
|
|
Wellenzahl=[]
|
|
|
Frequenz=[]
|
|
|
VibTemp=[]
|
|
|
Sv=[]
|
|
|
Inertia=[]
|
|
|
AvInertia=[]
|
|
|
Sr=[]
|
|
|
Gewichtung=[]
|
|
|
Entropie=[]
|
|
|
pass
|
|
|
dateiname=""
|
|
|
freq1=[]
|
|
|
freq2=[]
|
|
|
translation=[]
|
|
|
rotation=[]
|
|
|
Gesamtentropie=0
|
|
|
GaussianEntropie=0
|
|
|
|
|
|
Tempeingabe=input("Temperatur? [298.15 K]")
|
|
|
if not Tempeingabe=="":
|
|
|
Temp=float(Tempeingabe)
|
|
|
else:
|
|
|
Temp=298.15
|
|
|
print("T =",Temp)
|
|
|
|
|
|
TurbomolOderGaussian=input("[T]urbomol oder [G]aussian? ")
|
|
|
|
|
|
#GAUSSIAN
|
|
|
if TurbomolOderGaussian=="G":
|
|
|
while not os.path.isfile(dateiname):
|
|
|
print("Im aktuellen Verzeichnis sind folgende log-Dateien vorhanden: ")
|
|
|
dateien=glob.glob("*.log")
|
|
|
print(dateien)
|
|
|
dateieingabe=str(input("Log-Datei? [*.log]"))
|
|
|
print(dateieingabe[len(dateieingabe)-4:len(dateieingabe)-1])
|
|
|
if not dateieingabe[len(dateieingabe)-4:len(dateieingabe)-1]=="*.log":
|
|
|
dateiname=dateieingabe+".log"
|
|
|
else:
|
|
|
dateiname=dateieingabe
|
|
|
if not dateieingabe:
|
|
|
dateieingabe=dateien
|
|
|
if not os.path.isfile(dateiname):
|
|
|
print("Datei",dateiname,"existiert nicht...")
|
|
|
|
|
|
print("Log-File",dateiname,"wird verwendet...")
|
|
|
ausgabedatei=dateieingabe+".gibbs"
|
|
|
ausgabe = open(ausgabedatei,"w")
|
|
|
print("Ausgabe",ausgabedatei,"wird verwendet...")
|
|
|
|
|
|
print("Lese Datei",dateiname, "...")
|
|
|
datei=open(dateiname)
|
|
|
string = datei.read()
|
|
|
#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 open(dateiname):
|
|
|
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<EFBFBD>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)
|
|
|
|
|
|
#TURBOMOL
|
|
|
elif TurbomolOderGaussian=="T" or TurbomolOderGaussian=="":
|
|
|
#Dateiabfrage
|
|
|
dateiname=str(input("Ordner? [vibspectrum]"))
|
|
|
if not dateiname:
|
|
|
dateiname="vibspectrum"
|
|
|
else:
|
|
|
dateiname=dateiname+"/vibspectrum"
|
|
|
if not os.path.isfile(dateiname):
|
|
|
print("vibspectrum existiert nicht...")
|
|
|
else:
|
|
|
ausgabedatei=dateiname+".gibbs"
|
|
|
print("Ausgabe",ausgabedatei,"wird verwendet...")
|
|
|
print("Lese Datei",dateiname,"...")
|
|
|
datei=open(dateiname)
|
|
|
string = datei.read()
|
|
|
position=allindices(string, "a ")
|
|
|
j=0
|
|
|
for i in range(0,len(position)):
|
|
|
zahl=float(string[position[i]+10:position[i]+25])
|
|
|
if zahl>=0:
|
|
|
Daten.Wellenzahl.append(zahl)
|
|
|
Daten.Frequenz.append(Daten.Wellenzahl[j]*Lichtgeschwindigkeit*100)
|
|
|
Daten.VibTemp.append(Daten.Frequenz[j]*Planck/Boltzmann)
|
|
|
Daten.Sv.append(R*((Daten.VibTemp[j]/Temp)/((math.exp(Daten.VibTemp[j]/Temp))-1)-math.log(1-(math.exp((-1)*Daten.VibTemp[j]/Temp)))))
|
|
|
Daten.Inertia.append(Planck/(8*math.pi**2*Daten.Frequenz[j]))
|
|
|
Daten.AvInertia.append((Daten.Inertia[j]*Bav)/(Daten.Inertia[j]+Bav))
|
|
|
Daten.Sr.append(R*(0.5+math.log(((8*math.pi**3*Daten.AvInertia[j]*Boltzmann*Temp)/(Planck**2))**(0.5))))
|
|
|
Daten.Gewichtung.append(1/(1+(Grenze/Daten.Frequenz[j])**4))
|
|
|
Daten.Entropie.append(Daten.Gewichtung[j]*Daten.Sv[j]+((1-Daten.Gewichtung[j])*Daten.Sr[j]))
|
|
|
j=j+1
|
|
|
else:
|
|
|
print("Imaginäre Frequenz", zahl, "cm^-1 wird ignoriert...")
|
|
|
|
|
|
#Berechnung
|
|
|
Sv=sum(Daten.Sv)/CalInJoule/AuInKcal/1000
|
|
|
Sr=sum(Daten.Entropie)/CalInJoule/AuInKcal/1000
|
|
|
dS=Sr-Sv
|
|
|
|
|
|
#Anzeige
|
|
|
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+="S(harmonic):\t"+str(Sv)+" au/K\t"+str(Sv*AuInKcal*1000)+" cal/mol*K\t"+str(Sv*AuInKcal*CalInJoule*1000)+" J/mol*K\n"
|
|
|
schreiben+="S(Grimme):\t"+str(Sr)+ " au/K\t"+ str(Sr*AuInKcal*1000)+ " cal/mol*K\t"+ str(Sr*AuInKcal*CalInJoule*1000)+ " J/mol*K\n"
|
|
|
schreiben+="S(G)-S(h):\t"+ str(dS)+ " au/K\t"+ str(dS*AuInKcal*1000)+ " cal/mol*K\t"+ str(dS*AuInKcal*CalInJoule*1000)+ " J/mol*K\n"
|
|
|
schreiben+="************************************************************************************************************\n"
|
|
|
|
|
|
print(schreiben)
|
|
|
|
|
|
#Ausgabe
|
|
|
ausgabe=open(ausgabedatei, "w")
|
|
|
ausgabe.write(schreiben)
|
|
|
ausgabe.close()
|
|
|
else:
|
|
|
print("Was soll ich tun?")
|
|
|
|
|
|
#Weitermachen
|
|
|
abfrage=input("Noch eine Datei? [j/N]")
|
|
|
if not abfrage=="j": wiederholen=False
|
|
|
print("Copyright 2015 SAW")
|