#!/usr/bin/python # -*- coding: utf-8 -*- # #-------------------------------------------------------------------------------- #anamo.py v1.7, Copyright Bjoern Olausson #-------------------------------------------------------------------------------- #This program is free software; you can redistribute it and/or modify #it under the terms of the GNU General Public License as published by #the Free Software Foundation; either version 2 of the License, or #(at your option) any later version. # #This program is distributed in the hope that it will be useful, #but WITHOUT ANY WARRANTY; without even the implied warranty of #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #GNU General Public License for more details. # #To view the license visit #http://www.gnu.org/licenses/old-licenses/gpl-2.0.html #or write to #Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA #-------------------------------------------------------------------------------- #-------------------------------------------------------------------------------- # #This script is intended to pars NAMD output files and plot energy, temperature, and #whatever information NAMD writes down. #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #If you have different Values for "outputEnergies" and "outputPressure" use Time/TS if you want to plot #information affected by this value. Use PTime/PTS if you want to plot Information which is affcted by #this keyword. I didn't test this so far. If you use it double check the resulting plot. #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # #PRESSURE: PTS PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ #GPRESSURE: PTS PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ #PRESSAVG: PTS PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ #GPRESSAVG: PTS PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ #ENERGY: TS BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP TOTAL2 TOTAL3 TEMPAVG PRESSURE GPRESSURE VOLUME PRESSAVG GPRESSAVG ##---------- --------- --------- --------- --------- --------- # # PTS PXX *PXY PXZ *PYX PYY +PYZ PZX +PZY PZZ # PRESSURE: 29101000 113.817 270.048 -425.735 270.048 -577.887 116.238 -425.735 116.238 -25.2852 #GPRESSURE: 29101000 -131.017 120.797 -384.488 167.638 -452.296 173.325 -159.551 -48.9657 64.5387 # PRESSAVG: 29101000 -16.186 5.86015 4.60358 5.86015 -12.8187 -85.7577 4.60358 -85.7577 47.1837 #GPRESSAVG: 29101000 -15.7132 9.52987 5.07437 -4.25737 -14.0235 -91.353 3.07567 -79.6831 48.6288 # ENERGY: 29101000 2805.2103 8458.6888 4626.4405 53.0736 -39150.1769 -1947.9501 0.0000 0.0000 13354.4343 -11800.2794 301.4109 -11674.5603 -11671.2349 300.4348 -163.1183 -172.9249 182937.7460 6.0596 6.2974 #---------- --------- --------- --------- --------- --------- import sys, os, shutil, string, glob, subprocess, itertools, pylab, scipy, re, time from optparse import OptionParser, OptionGroup opa = OptionParser() opa.add_option("-f", action="append", dest="chmout", metavar="FILE", help="NAMD out file") opa.add_option("-o", action="store", dest="saveplot", metavar="FILE", help="Filename to which we save the plot (.png) and RAW data (.dat)") opa.add_option("-s", action="store_true", dest="singleplot", help="Plot to single image") group = OptionGroup(opa, "You can choose freely from the following Values which will occure on either X or Y Axis", "TS Time PTS PTime BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP TOTAL2 TOTAL3 TEMPAVG PRESSURE GPRESSURE VOLUME PRES SAVG GPRESSAVG PTS PXX PXY PXZ PYX PYY PYZ PZX PZY PZZ GPXX GPXY GPXZ GPYX GPYY GPYZ GPZX GPZY GPZZ PXXAV PXYAV PXZAV PYXAV PYYAV PYZAV PZXAV PZYAV PZZAV GPXXAV GPXYAV GPXZAV GPYXAV GPYYAV GPYZAV GPZXAV GPZYAV GPZZAV (If you use one of Pressure Tensors, use for TS PTS or for Time PTime)") group.add_option("-m", action="store_true", dest="Mini", help="Use this option to only analys the minimization part of the output file. Default is, to skip the minimization steps and only plot the MD part") group.add_option("-n", action="store_true", dest="MDMIN", help="Use this option to analys the MD and minimization part of the output file.") group.add_option("-x", action="append", dest="xax", help="Choose Value from above you want to appear on the X-axis.") group.add_option("-y", action="append", dest="yax", help="Choose Value from above you want to appear on the Y-axis. If you want to plot multible Y values, you have to asign every Y value a X value and vice versa. For example: anaout.py -x TS -x TS -y PRESSURE -y TEMPAVG") opa.add_option_group(group) # instruct optparse to parse the program's command line (options, args) = opa.parse_args() # Checking for required otions if not options.chmout: parser.error("-f must be set") #chmout = options.chmout.split(" ") chmout = options.chmout Mini = options.Mini MDMIN = options.MDMIN xax = options.xax yax = options.yax saveplot = options.saveplot singleplot = options.singleplot #If "*" is found in a parameter, glob it. for element in chmout : if '*' in element : chmout.extend(glob.glob(element)) chmout.remove(element) chmout.sort() #Defining the lists TIMES = [] PTS = [] TS = []; BOND = []; ANGLE = []; DIHED = []; IMPRP = []; ELECT = []; VDW = []; BOUNDARY = []; MISC = []; KINETIC = []; TOTAL = []; TEMP = []; TOTAL2 = []; TOTAL3 = []; TEMPAVG = []; PRESSURE = []; GPRESSURE = []; VOLUME = []; PRESSAVG = []; GPRESSAVG = [] PXX = []; PXY = []; PXZ = []; PYX = []; PYY = []; PYZ = []; PZX = []; PZY = []; PZZ = [] GPXX = []; GPXY = []; GPXZ = []; GPYX = []; GPYY = []; GPYZ = []; GPZX = []; GPZY = []; GPZZ = [] PXXAV = []; PXYAV = []; PXZAV = []; PYXAV = []; PYYAV = []; PYZAV = []; PZXAV = []; PZYAV = []; PZZAV = [] GPXXAV = []; GPXYAV = []; GPXZAV = []; GPYXAV = []; GPYYAV = []; GPYZAV = []; GPZXAV = []; GPZYAV = []; GPZZAV = [] #Regular expression to find the start of the MD simulation #MDp = re.compile('Info: Finished startup') MDstart = re.compile('TCL: Running for [0-9]+ steps') #Regular expression to find the Timestep used TSre = re.compile('Info: +TIMESTEP +[1-9]+') #Regular expression to see if "minimization on" was used in the config file MINonly = re.compile('CONJUGATE GRADIENT MINIMIZATION ACTIVE') #Regular expression to see if "mini" was run MINstart = re.compile('TCL: Minimizing for [0-9]+ steps') #Info: TIMESTEP 2 def Info(line) : if TSre.match(line): MyTS = re.sub("\D", "", line) if MyTS.isdigit(): TIMES.append(int(MyTS)) else: raise Exception('Timestep is not a digit') #PRESSURE: 28601100 -13.0939 160.713 -191.914 160.713 -540.017 -233.812 -191.914 -233.812 1224.65 def Pressure(line) : PTS.append(int(line.split()[1])) PXX.append(float(line.split()[2])) PXY.append(float(line.split()[3])) PXZ.append(float(line.split()[4])) PYX.append(float(line.split()[5])) PYY.append(float(line.split()[6])) PYZ.append(float(line.split()[7])) PZX.append(float(line.split()[8])) PZY.append(float(line.split()[9])) PZZ.append(float(line.split()[10])) #GPRESSURE: 28601100 35.4952 117.701 -181.483 248.367 -208.371 -129.89 -404.508 -349.542 913.535 def GPressure(line) : GPXX.append(float(line.split()[2])) GPXY.append(float(line.split()[3])) GPXZ.append(float(line.split()[4])) GPYX.append(float(line.split()[5])) GPYY.append(float(line.split()[6])) GPYZ.append(float(line.split()[7])) GPZX.append(float(line.split()[8])) GPZY.append(float(line.split()[9])) GPZZ.append(float(line.split()[10])) #PRESSAVG: 28601100 11.1302 78.4552 38.3246 78.4552 2.34895 -59.1307 38.3246 -59.1307 2.71294 def PresAv(line) : PXXAV.append(float(line.split()[2])) PXYAV.append(float(line.split()[3])) PXZAV.append(float(line.split()[4])) PYXAV.append(float(line.split()[5])) PYYAV.append(float(line.split()[6])) PYZAV.append(float(line.split()[7])) PZXAV.append(float(line.split()[8])) PZYAV.append(float(line.split()[9])) PZZAV.append(float(line.split()[10])) #GPRESSAVG: 28601100 6.71473 68.2579 37.7828 76.8936 1.54956 -64.4147 42.1503 -46.2259 3.81618 def GPresAv(line) : GPXXAV.append(float(line.split()[2])) GPXYAV.append(float(line.split()[3])) GPXZAV.append(float(line.split()[4])) GPYXAV.append(float(line.split()[5])) GPYYAV.append(float(line.split()[6])) GPYZAV.append(float(line.split()[7])) GPZXAV.append(float(line.split()[8])) GPZYAV.append(float(line.split()[9])) GPZZAV.append(float(line.split()[10])) #ENERGY: 29101000 2805.2103 8458.6888 4626.4405 53.0736 -39150.1769 -1947.9501 0.0000 0.0000 13354.4343 -11800.2794 301.4109 -11674.5603 -11671.2349 300.4348 -163.1183 -172.9249 182937.7460 6.0596 6.2974 #ENERGY: TS BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP TOTAL2 TOTAL3 TEMPAVG PRESSURE GPRESSURE VOLUME PRESSAVG GPRESSAVG def Energy(line) : TS.append(int(line.split()[1])) BOND.append(float(line.split()[2])) ANGLE.append(float(line.split()[3])) DIHED.append(float(line.split()[4])) IMPRP.append(float(line.split()[5])) ELECT.append(float(line.split()[6])) VDW.append(float(line.split()[7])) BOUNDARY.append(float(line.split()[8])) MISC.append(float(line.split()[9])) KINETIC.append(float(line.split()[10])) TOTAL.append(float(line.split()[11])) TEMP.append(float(line.split()[12])) TOTAL2.append(float(line.split()[13])) TOTAL3.append(float(line.split()[14])) TEMPAVG.append(float(line.split()[15])) #TS.append(int(line[7:16].strip())) #BOND.append(float(line[16:31].strip())) #ANGLE.append(float(line[31:46].strip())) #DIHED.append(float(line[46:61].strip())) #IMPRP.append(float(line[61:76].strip())) #ELECT.append(float(line[76:96].strip())) #VDW.append(float(line[96:112].strip())) #BOUNDARY.append(float(line[112:126].strip())) #MISC.append(float(line[126:141].strip())) #KINETIC.append(float(line[141:156].strip())) #TOTAL.append(float(line[156:176].strip())) #TEMP.append(float(line[176:191].strip())) #TOTAL2.append(float(line[191:206].strip())) #TOTAL3.append(float(line[206:221].strip())) #TEMPAVG.append(float(line[221:236].strip())) # If "minimization on" was used in the config file, these values are not populated if MINs or MINo : PRESSURE.append(float(0)) GPRESSURE.append(float(0)) VOLUME.append(float(0)) PRESSAVG.append(float(0)) GPRESSAVG.append(float(0)) else: PRESSURE.append(float(line.split()[16])) GPRESSURE.append(float(line.split()[17])) VOLUME.append(float(line.split()[18])) PRESSAVG.append(float(line.split()[19])) GPRESSAVG.append(float(line.split()[20])) LineFunc = { 'Info': Info, 'PRESSURE': Pressure, 'GPRESSURE': GPressure, 'PRESSAVG': PresAv, 'GPRESSAVG': GPresAv, 'ENERGY': Energy } fnum = 0 nMDs = 0 nMINs = 0 nMINo = 0 for outfile in chmout : MINo = "" MINs = "" MDs = "" #Add Zeros for the first (Group-)Pressure averages cause there are no #averages for the first timestep in the outfile when you use the keyword #printGroupPressure PresAv('0 0 0 0 0 0 0 0 0 0 0') GPresAv('0 0 0 0 0 0 0 0 0 0 0') mm = False for lines in open(outfile, "r") : #Find TIMESTEP and extract it if TSre.search(lines): LineFunc[LineStart[0]](lines) #Check if this is a "minimization only" run if MINonly.search(lines): MINo = True MINs = False MDs = False nMINo += 1 print "%s mini only part(s) found" %(nMINo) #Check if a minimization was run prior to MD if MINstart.search(lines): MINs = True MDs = False nMINs += 1 print "%s mini part(s) found" %(nMINs) #Check for "MD" start if MDstart.search(lines): MDs = True MINs = False nMDs += 1 print "%s md part(s) found" %(nMDs) #Get the Line description to invoke the corresponding function LineStart = lines.split(":") #Per default we analyse only the MD run part of the logfile #Check if the MD startpoint is reached and start grabbing the lines if MDs and not MINs and not Mini and not MDMIN : if LineStart[0] in LineFunc : LineFunc[LineStart[0]](lines) #Only analyse the minimizations bevor/between/after a MD run if Mini: #Check if the MINIMIZATION startpoint is reached if MINs and not MDs: if LineStart[0] in LineFunc : LineFunc[LineStart[0]](lines) #Analyse both MINIMIZATION and MD part if MDMIN: if MINs or MDs: if LineStart[0] in LineFunc : LineFunc[LineStart[0]](lines) #If this is a MINIMIZATION only run... we default to this instead to the MD part if MINo: if LineStart[0] in LineFunc : LineFunc[LineStart[0]](lines) fnum = fnum + 1 print print "---------------------DONE---------------------" if nMDs == 0 and nMINo == 0 and not Mini: print "Only mini part(s) (%s) found, but there was no option passed to analyse minimization" %(nMINs) raise Exception("Add \"-m\" to your options and run again") else: print "%s mini part(s) found overall" %(nMINs) print "%s mini only part(s) found overall" %(nMINo) print "%s md part(s) found overall" %(nMDs) print "---------------------DONE---------------------" #Convert Lists to scipy arrays TS = scipy.array(TS) PTS = scipy.array(PTS) #Calculate time from Timesteps in ns Time = TS Time = Time * TIMES[0] * 1e-6 PTime = PTS PTime = PTime * TIMES[0] * 1e-6 BOND = scipy.array(BOND) ANGLE = scipy.array(ANGLE) DIHED = scipy.array(DIHED) IMPRP = scipy.array(IMPRP) ELECT = scipy.array(ELECT) VDW = scipy.array(VDW) BOUNDARY = scipy.array(BOUNDARY) MISC = scipy.array(MISC) KINETIC = scipy.array(KINETIC) TOTAL = scipy.array(TOTAL) TEMP = scipy.array(TEMP) TOTAL2 = scipy.array(TOTAL2) TOTAL3 = scipy.array(TOTAL3) TEMPAVG = scipy.array(TEMPAVG) PRESSURE = scipy.array(PRESSURE) GPRESSURE = scipy.array(GPRESSURE) VOLUME = scipy.array(VOLUME) PRESSAVG = scipy.array(PRESSAVG) PXX = scipy.array(PXX) PXY = scipy.array(PYY) PXZ = scipy.array(PZZ) PYX = scipy.array(PYX) PYY = scipy.array(PYY) PYZ = scipy.array(PYZ) PZX = scipy.array(PZX) PZY = scipy.array(PZY) PZZ = scipy.array(PZZ) GPXX = scipy.array(GPXX) GPXY = scipy.array(GPXY) GPXZ = scipy.array(GPXZ) GPYX = scipy.array(GPYX) GPYY = scipy.array(GPYY) GPYZ = scipy.array(GPYZ) GPZX = scipy.array(GPZX) GPZY = scipy.array(GPZY) GPZZ = scipy.array(GPZZ) PXXAV = scipy.array(PXXAV) PXYAV = scipy.array(PXYAV) PXZAV = scipy.array(PXZAV) PYXAV = scipy.array(PYXAV) PYYAV = scipy.array(PYYAV) PYZAV = scipy.array(PYZAV) PZXAV = scipy.array(PZXAV) PZYAV = scipy.array(PZYAV) PZZAV = scipy.array(PZZAV) GPXXAV = scipy.array(GPXXAV) GPXYAV = scipy.array(GPXYAV) GPXZAV = scipy.array(GPXZAV) GPYXAV = scipy.array(GPYXAV) GPYYAV = scipy.array(GPYYAV) GPYZAV = scipy.array(GPYZAV) GPZXAV = scipy.array(GPZXAV) GPZYAV = scipy.array(GPZYAV) GPZZAV = scipy.array(GPZZAV) #COMBINED = scipy.array([[TS], [Time], [TOTEner], [TOTKe], [ENERgy], [TEMPerature], [GRMS], [HFCTote], [HFCKe], [EHFCor], [VIRKe], [BONDs], [ANGLes], [UREYb], [DIHEdrals], [IMPRopers], [CMAPs], [VDWaals], [ELEC], [HBONds], [ASP], [USER], [IMNBvdw], [IMELec], [IMHBnd], [RXNField], [EXTElec], [EWKSum], [EWSElf], [EWEXcl], [EWQCor], [EWUTil], [HARMonic], [CDIHedral], [CIC], [RESDistance], [NOE], [GEO], [MDIP], [SSBP], [SHEL], [DROFfa], [VIRE], [VIRI], [PRESSE], [PRESSI], [VOLUme]]) PlotDict = { 'TS': TS, 'PTS': PTS, 'Time': Time, 'PTime': PTime, 'BOND': BOND, 'ANGLE': ANGLE, 'DIHED': DIHED, 'IMPRP': IMPRP, 'ELECT': ELECT, 'VDW': VDW, 'BOUNDARY': BOUNDARY, 'MISC': MISC, 'KINETIC': KINETIC, 'TOTAL': TOTAL, 'TEMP': TEMP, 'TOTAL2': TOTAL2, 'TOTAL3': TOTAL3, 'TEMPAVG': TEMPAVG, 'PRESSURE': PRESSURE, 'GPRESSURE': GPRESSURE, 'VOLUME': VOLUME, 'PRESSAVG': PRESSAVG, 'GPRESSAVG': GPRESSAVG, 'PXX': PXX, 'PXY': PXY, 'PXZ': PXZ, 'PYX': PYX, 'PYY': PYY, 'PYZ': PYZ, 'PZX': PZX, 'PZY': PZY, 'PZZ': PZZ, 'GPXX': GPXX, 'GPXY': GPXY, 'GPXZ': GPXZ, 'GPYX': GPYX, 'GPYY': GPYY, 'GPYZ': GPYZ, 'GPZX': GPZX, 'GPZY': GPZY, 'GPZZ': GPZZ, 'PXXAV': PXXAV, 'PXYAV': PXYAV, 'PXZAV': PXZAV, 'PYXAV': PYXAV, 'PYYAV': PYYAV, 'PYZAV': PYZAV, 'PZXAV': PZXAV, 'PZYAV': PZYAV, 'PZZAV': PZZAV, 'GPXXAV': GPXXAV, 'GPXYAV': GPXYAV, 'GPXZAV': GPXZAV, 'GPYXAV': GPYXAV, 'GPYYAV': GPYYAV, 'GPYZAV': GPYZAV, 'GPZXAV': GPZXAV, 'GPZYAV': GPZYAV, 'GPZZAV': GPZZAV, } def SaveData(Xvalues,Yvalues, Xname, Yname, saveplot) : if saveplot : savedata = Xname+"-"+Yname+"_"+saveplot+".dat" fout = open(savedata, "w") for i,j in itertools.izip(Xvalues,Yvalues) : #print i fout.writelines("%s %s\n" % (i,j)) fout.close() def SaveImage(Xname, Yname, saveplot) : if saveplot : saveimg = Xname+"-"+Yname+"_"+saveplot+".png" pylab.savefig(saveimg) def DynaPltSingle(Xvalues,Yvalues,Xname,Yname) : #Do some nice plots later SaveData(Xvalues, Yvalues, Xname, Yname, saveplot) pylab.plot(Xvalues, Yvalues) Xlegend.append(Xname) Ylegend.append(Yname) def DynaPltMulti(Xvalues,Yvalues,Xname,Yname, figure) : pylab.figure(figure) pylab.plot(Xvalues, Yvalues) pylab.xlabel(Xname) pylab.ylabel(Yname) SaveImage(Xname, Yname, saveplot) SaveData(Xvalues,Yvalues, Xname, Yname, saveplot) fign = 0 Xlegend = [] Xlabel = [] Ylegend = [] Ylabel = [] if singleplot : for xAxis,yAxis in itertools.izip(xax,yax) : DynaPltSingle(PlotDict[xAxis], PlotDict[yAxis], xAxis, yAxis) pylab.legend((Ylegend), shadow = True) pylab.ylabel(Ylegend) pylab.xlabel(Xlegend) SaveImage(xAxis, 'MULTI', saveplot) else : figure = 0 for xAxis,yAxis in itertools.izip(xax,yax) : DynaPltMulti(PlotDict[xAxis], PlotDict[yAxis], xAxis, yAxis, figure) figure = figure + 1 pylab.show()