Gaussian09 optimization trajectory: python script

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 subprocess

def 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=0
    coordindex_begin=[]
    coordindex_end=[]
    for line in indata:
        lineindex=lineindex+1
        if "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_end
    indata.close()
    return coordindex_begin,coordindex_end

def get_oneframe(logfile,index_begin,index_end):
    index0=index_begin
    index1=index_end
    indata=open(logfile,'r')
    outdata=open(logfile.split('.')[0]+'_traj.xyz','a')
    outdata.write(str(index1-index0+1)+'\n'+'\n')
    lineindex=0
    for line in indata:
        lineindex=lineindex+1
        words=line.strip().split()
        if lineindex >= index0 and lineindex <= index1:
            outdata.write(get_atomlabel(words[1])+"  "+words[3]+"  "+words[4]+"  "+words[5]+'\n')
        else :
            continue
    indata.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 :
            continue
    indata.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"
            continue
    indata.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"
            continue
    indata.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()