Gaussian09 optimization trajectory: python script

• Note
• • code

Note

Generally, the optimizing trajectory could be viewed in Gaussview by loading the log file when you checked the box “Read Intermediate Geometries”. However, there maybe some noisy jump between two adjacent frames. This happens because Gaussview read coordinates represented by “Standard Orientation” by default. In order to avoid such problem, we could extract the “Input Oritentation” coordinates in the log file, and then load it into the VMD program. I wrote a piece of python script to make it work.

code

高亮的 代碼.

// An highlighted block import sys import subprocessdef get_atomlabel(atomic_number):element_dict = {"1" : "H", "2" : "He", "3" : "Li", "4" : "Be", "5" : "B", \"6" : "C", "7" : "N", "8" : "O", "9" : "F", "10" : "Ne", \"11" : "Na" , "12" : "Mg" , "13" : "Al" , "14" : "Si" , "15" : "P", \"16" : "S" , "17" : "Cl" , "18" : "Ar" , "19" : "K" , "20" : "Ca", \"21" : "Sc" , "22" : "Ti" , "23" : "V" , "24" : "Cr" , "25" : "Mn", \"26" : "Fe" , "27" : "Co" , "28" : "Ni" , "29" : "Cu" , "30" : "Zn", \"31" : "Ga" , "32" : "Ge" , "33" : "As" , "34" : "Se" , "35" : "Br", \"36" : "Kr" , "37" : "Rb" , "38" : "Sr" , "39" : "Y" , "40" : "Zr", \"41" : "Nb" , "42" : "Mo" , "43" : "Tc" , "44" : "Ru" , "45" : "Rh", \"46" : "Pd" , "47" : "Ag" , "48" : "Cd" , "49" : "In" , "50" : "Sn", \"51" : "Sb" , "52" : "Te" , "53" : "I" , "54" : "Xe" , "55" : "Cs", \"56" : "Ba" , "57" : "La" , "58" : "Ce" , "59" : "Pr" , "60" : "Nd", \"61" : "Pm" , "62" : "Sm" , "63" : "Eu" , "64" : "Gd" , "65" : "Tb", \"66" : "Dy" , "67" : "Ho" , "68" : "Er" , "69" : "Tm" , "70" : "Yb", \"71" : "Lu" , "72" : "Hf" , "73" : "Ta" , "74" : "W" , "75" : "Re", \"76" : "Os" , "77" : "Ir" , "78" : "Pt" , "79" : "Au" , "80" : "Hg", \"81" : "Tl" , "82" : "Pb" , "83" : "Bi" , "84" : "Po" , "85" : "At", \"86" : "Rn" , "87" : "Fr" , "88" : "Ra" , "89" : "Ac" , "90" : "Th", \"91" : "Pa" , "92" : "U" , "93" : "Np" , "94" : "Pu" , "95" : "Am", \"96" : "Cm" , "97" : "Bk" , "98" : "Cf" , "99" : "Es" ,"100" : "Fm", \"101": "Md" ,"102" : "No" ,"103" : "Lr" ,"104" : "Rf" ,"105" : "Db", \"106": "Sg" ,"107" : "Bh" ,"108" : "Hs" ,"109" : "Mt" ,"110" : "Ds", \"111": "Rg" ,"112" : "Uub","113" : "Uut","114" : "Uuq","115" : "Uup", \"116": "Uuh","117" : "Uus","118" : "Uuo"}return element_dict[str(atomic_number)]def get_coordindex(logfile):indata=open(logfile,'r')lineindex=0coordindex_begin=[]coordindex_end=[]for line in indata:lineindex=lineindex+1if "Input orientation" in line:coordindex_begin.append(lineindex+5)if "Distance matrix" in line:coordindex_end.append(lineindex-2)coordindex_begin.sort()coordindex_end.sort()#print coordindex_begin#print coordindex_endindata.close()return coordindex_begin,coordindex_enddef get_oneframe(logfile,index_begin,index_end):index0=index_beginindex1=index_endindata=open(logfile,'r')outdata=open(logfile.split('.')[0]+'_traj.xyz','a')outdata.write(str(index1-index0+1)+'\n'+'\n')lineindex=0for line in indata:lineindex=lineindex+1words=line.strip().split()if lineindex >= index0 and lineindex <= index1:outdata.write(get_atomlabel(words[1])+" "+words[3]+" "+words[4]+" "+words[5]+'\n')else :continueindata.close()outdata.close()def get_scfenergy(logfile):indata=open(logfile,'r')outdata=open(logfile.split('.')[0]+'_scfEnergy.dat','w')for line in indata:words=line.strip().split()if "SCF Done" in line:outdata.write(words[4]+'\n')else :continueindata.close()outdata.close()def get_frequencies(logfile):indata=open(logfile,'r')outdata=open(logfile.split('.')[0]+'_freq.dat','w')for line in indata:words=line.strip().split()if "Frequencies" in line:outdata.write(words[2]+'\n'+words[3]+'\n'+words[4]+'\n')else :#print "No Frequencies Found"continueindata.close()outdata.close()def get_maxforce(logfile):indata=open(logfile,'r')outdata=open(logfile.split('.')[0]+'_maxforce.dat','w')for line in indata:words=line.strip().split()if "Cartesian Forces" in line:outdata.write(words[3]+'\n')else :#print "No Frequencies Found"continueindata.close()outdata.close()def main():logfile=sys.argv[1]logfile_short=logfile.split('.')[0]subprocess.call(['rm',logfile_short+'_traj.xyz'])coordindex_begin,coordindex_end=get_coordindex(logfile)for i in range(len(coordindex_begin)):get_oneframe(logfile,coordindex_begin[i],coordindex_end[i])#get_scfenergy(logfile)#get_maxforce(logfile)#get_frequencies(logfile)if __name__ == '__main__': main()

總結(jié)

以上是生活随笔為你收集整理的Gaussian09 optimization trajectory: python script的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。