Sugestoes de programas para manipular PDF1 - PDF mod2 - Gscan2pdf3 - ScribusNG4 - PDF Editore de quebra, uma ferramenta de multi-rename: GPRenameTudo via GUI ok ?!Abracos !!! Extracting data from a PES scan with gaussianThere are a few reasons to like gaussian, and many reasons not to. Gaussian is fast, and their whitepapers are great resources for learning computational techniques.Without going into discussions about the commercial behaviour of Wavefunction inc., the things I don't like about gaussian is the clunky input format (nwchem has a much more readable syntax), the inscrutable error messages, and the unreadable output. Well, it's not unreadable in a literal sense, but it could certainly be clearer. On the other hand, I've having issues with running some of my PES scans in nwchem -- and I can't find a solution (more about that in a later post)Anyway, here's a python script for extracting optimized structures and energies from a relaxed PES scan in Gaussian 09.First, an example of a simple scan:%nprocshared=2%Chk=methanol.chk#P rB3LYP/6-31g 6D 10F Opt=(modredundant) NoSymm Punch=(MO) Pop=() methanol0 1 ! charge and multiplicity C 0.0351714 0.00548884 0.0351714 H -0.617781 -0.634073 0.667983 H 0.667983 -0.634073 -0.617781 H -0.605139 0.646470 -0.605139 O 0.839603 0.818768 0.839603 H 1.38912 0.201564 1.38912 1 5 S 10 0.1And here's the script, pes_parse_g09:#!/usr/bin/pythonimport sysdef getrawdata(infile): f=open(infile,'r') opt=0 geo=0 struct=[] structure=[] energies=[] energy=[] for line in f: if opt==1 and geo==1 and not ("---" in line): structure+=[line.rstrip()] if 'Coordinates (Angstroms)' in line: if opt==0: opt=1 structure=[] if opt==1 and "--------------------------" in line: if geo==0: geo=1 elif geo==1: geo=0 opt=0 if 'SCF Done' in line: energy=filter(None,line.rstrip('\n').split(' ')) if 'Optimization completed' in line and (opt==0 and geo==0): energies+=[float(energy[4])] opt=0 geo=0 struct+=[structure] structure=[] return struct, energiesdef periodictable(elementnumber): ptable={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:'Cn',\ 113:'Uut',114:'Fl',115:'Uup',116:'Lv',117:'Uus',118:'Uuo'} element=ptable[elementnumber] return elementdef genxyzstring(coords,elementnumber): x_str='%10.5f'% coords[0] y_str='%10.5f'% coords[1] z_str='%10.5f'% coords[2] element=periodictable(int(elementnumber)) xyz_string=element+(3-len(element))*' '+10*' '+\ (8-len(x_str))*' '+x_str+10*' '+(8-len(y_str))*' '+y_str+10*' '+(8-len(z_str))*' '+z_str+'\n' return xyz_stringdef getstructures(rawdata): n=0 for structure in rawdata: n=n+1 num="%03d" % (n,) g=open('structure_'+num+'.xyz','w') itson=False cartesian=[] for item in structure: coords=filter(None,item.split(' ')) coordinates=[float(coords[3]),float(coords[4]),float(coords[5])] element=coords[1] cartesian+=[genxyzstring(coordinates,element)] g.write(str(len(cartesian))+'\n') g.write('Structure '+str(n)+'\n') for line in cartesian: g.write(line) g.close() cartesian=[] return 0 if __name__ == "__main__": infile=sys.argv[1] rawdata,energies=getrawdata(infile) structures=getstructures(rawdata) g=open('energies.dat','w') for n in range(0,len(energies)): g.write(str(n)+'\t'+str(energies[n])+'\n') g.close()And here's what we get from the output:g09 methanol.in |tee methanol.outpes_parse_g09 methanol.logcat structure* > meoh_traj.xyzAnd here's a plot of energies.dat:Posted by lindqvist at 19:28 Email ThisBlogThis!Share to TwitterShare to FacebookShare to PinterestLabels: extract, g09, PES, potential energy surface, python, script6 comments:Anonymous11 July, 2014 08:32you are a god among men. thank you!ReplyReplieslindqvist11 July, 2014 08:42I shall direct the head of the department and the dean of the faculty to this Comment, so that my exulted status will be recognised and tenure unconditionally granted.Either way, happy it helped!ReplyAnonymous08 November, 2014 06:32Thanks for posting this, it's really helpful. One mistake, easily corrected: line 54 should terminate with ,\ReplyReplieslindqvist10 November, 2014 13:57Thanks for spotting that and letting me know -- fixed now.ReplyAnonymous16 June, 2015 19:20Thanks, very useful!ReplyAnonymous07 September, 2015 17:22Very useful, thanks a lot for sharing.Reply Esse tutorial mostra