# -*- coding: utf-8 -*-
#
#--------------------------------------------------------------------------------
#anamo.py v1.6, 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
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:
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) :
fout.writelines("%s %s
" % (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()
Home