#!/usr/bin/env python3 import math import os import glob ########################################### version="0.8" 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 calc_gaussian(dateiname): class Daten: Wellenzahl=[] Frequenz=[] VibTemp=[] Sv=[] Inertia=[] AvInertia=[] Sr=[] Gewichtung=[] Entropie=[] pass freq1=[] freq2=[] print("Log-File",dateiname,"wird verwendet...") dateieingabe[len(dateieingabe)-4:len(dateieingabe)-1] ausgabedatei=dateiname[0:len(dateiname)-4]+".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] test3=save2.find("\\") save3=save2[0:test3] Energy=float(save1+save3) else: test2=ausschneiden.find("\\") if test2==-1: print("kein \\ gefunden") Energy=float(string[position+4:position+15]) else: print("\\ befindet sich an Stelle",test2) Energy=float(string[position+4:position+3+test2]) #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)): #global freq2 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("Okay") ########################################### 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: datei="" 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? (Default: Gaussian)") #GAUSSIAN if TurbomolOderGaussian=="G" or TurbomolOderGaussian=="": while not os.path.isfile(datei): print("Im aktuellen Verzeichnis sind folgende log-Dateien vorhanden: ") dateien=glob.glob("*.log") print(dateien) dateieingabe=str(input("Log-Datei? [*.log, [A]lle]")) if not dateieingabe[len(dateieingabe)-4:len(dateieingabe)-1]=="*.log": datei=dateieingabe+".log" else: datei=dateieingabe if dateieingabe=="A": for datei in dateien: calc_gaussian(datei) else: if not os.path.isfile(datei): print("Datei",datei,"existiert nicht...") calc_gaussian(datei) #TURBOMOL elif TurbomolOderGaussian=="T" or TurbomolOderGaussian=="": class Daten: Wellenzahl=[] Frequenz=[] VibTemp=[] Sv=[] Inertia=[] AvInertia=[] Sr=[] Gewichtung=[] Entropie=[] pass freq1=[] freq2=[] translation=[] rotation=[] Gesamtentropie=0 GaussianEntropie=0 #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")