anamo.py

Tuesday, 08 December 2009 00:00
anamo.py Version:2.1
This script is intended to pars NAMD output files and plot energy, temperature, and whatever information NAMD writes down. So far you can choose the following keyworkds for a classical MD: 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 For accelerated MD the following keywords are supported: aTS aTime adV adVAVG aBOND aANGLE aDIHED aIMPRP aELECT aVDW aPOTENTIAL For every selection min/max (with the corresponding frame number) , mean, standard deviation and sum is calculated.

GNU/GPL 2011-05-10 English Linux 20.3 KB 190
#!/usr/bin/python
# -*- coding: utf-8 -*-
#
#--------------------------------------------------------------------------------
#anamo.py v2.1, 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.
#Same applies for aMD. use aTS for the timesteps.
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
#
#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         POTENTIAL         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
usage = "usage: %prog [options] arg"
version = "2.1"
opa = OptionParser(usage=usage, version=version)
 
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 keyworks which will occure on either X or Y Axis (keywords prefixed with \"a\" are for accelerated MDs)",
"TS Time PTS PTime BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP POTENTIAL 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) \
aTS, aTime, adV, adVAVG, aBOND, aANGLE, aDIHED, aIMPRP, aELECT, aVDW, aPOTENTIAL (use aTS or aTIME instead of TS or TIME when you choose keywords for accelerated MDs")
 
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:
	opa.print_help()
	opa.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
 
# Function for natural sorting!
def sort_nicely( l ):
	""" Sort the given list in the way that humans expect.
	http://www.codinghorror.com/blog/2007/12/sorting-for-humans-natural-sort-order.html"""
	convert = lambda text: int(text) if text.isdigit() else text
	alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ]
	l.sort( key=alphanum_key )
	return l
 
#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()
 
# Sort the list naturally (1 20 30 100 200 300)
chmout = sort_nicely(chmout)
 
#Defining the lists
TIMES = []
PTS = []
TS = []; BOND = []; ANGLE = []; DIHED = []; IMPRP = []; ELECT = []; VDW = []; BOUNDARY = []; MISC = []; KINETIC = []; TOTAL = []; TEMP = []; POTENTIAL = []; 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 = []
# Accelerated MD:
aTS = []; adV = []; adVAVG = []; aBOND = []; aANGLE = []; aDIHED = []; aIMPRP = []; 	aELECT = []; 	aVDW = []; aPOTENTIAL = []
 
 
#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         POTENTIAL         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]))
	POTENTIAL.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()))
	#POTENTIAL.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]))
 
#ACCELERATED MD: STEP 500800 dV 133.039 dVAVG 138.397 BOND 639.477 ANGLE 1739.73 DIHED 1773.13 IMPRP 140.725 ELECT -4768.01 VDW -542.019 POTENTIAL -1016.95
def aMD(line):
	aTS.append(int(line.split()[3]))
	adV.append(float(line.split()[5]))
	adVAVG.append(float(line.split()[7]))
	aBOND.append(float(line.split()[9]))
	aANGLE.append(float(line.split()[11]))
	aDIHED.append(float(line.split()[13]))
	aIMPRP.append(float(line.split()[15]))
	aELECT.append(float(line.split()[17]))
	aVDW.append(float(line.split()[19]))
	aPOTENTIAL.append(float(line.split()[21]))
 
LineFunc = {
	'Info': Info,
	'PRESSURE': Pressure,
	'GPRESSURE': GPressure,
	'PRESSAVG': PresAv,
	'GPRESSAVG': GPresAv,
	'ENERGY': Energy,
	'ACCELERATED MD': aMD,
}
 
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)
POTENTIAL = scipy.array(POTENTIAL)
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)
aTS = scipy.array(aTS)
aTime = aTS
aTime = aTime * TIMES[0] * 1e-6
adV = scipy.array(adV)
adVAVG = scipy.array(adVAVG)
aBOND = scipy.array(aBOND)
aANGLE = scipy.array(aANGLE)
aDIHED = scipy.array(aDIHED)
aIMPRP = scipy.array(aIMPRP)
aELECT = scipy.array(aELECT)
aVDW = scipy.array(aVDW)
aPOTENTIAL = scipy.array(aPOTENTIAL)
#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,
	'POTENTIAL': POTENTIAL,
	'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,
	'aTS': aTS,
	'adV': adV,
	'adVAVG': adVAVG,
	'aBOND': aBOND,
	'aANGLE': aANGLE,
	'aDIHED': aDIHED,
	'aIMPRP': aIMPRP,
	'aELECT': aELECT,
	'aVDW': aVDW,
	'aPOTENTIAL': aPOTENTIAL,
}
 
# Print the average for the selected values
print "Some stats for the selected values:"
for XSELECTION,YSELECTION in itertools.izip(xax,yax):
	XDATA = PlotDict[XSELECTION]
	YDATA = PlotDict[YSELECTION]	
	print "%(YSELECTION)s" %{"YSELECTION": YSELECTION}
	print "MIN: %(MIN)s@%(FMIN)s Max: %(MAX)s@%(FMAX)s MEAN: %(MEAN)s STD: %(STD)s | SUM: %(SUM)s" %{"MIN": YDATA.min(), "FMIN": XDATA[YDATA.argmin()], "MAX": YDATA.max(), "FMAX": XDATA[YDATA.argmax()], "MEAN": YDATA.mean(), "STD": YDATA.std(), "SUM": YDATA.sum()}
 
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_png = Xname+"-"+Yname+"_"+saveplot+".png"
		saveimg_svg = Xname+"-"+Yname+"_"+saveplot+".svg"
		pylab.savefig(saveimg_png)
		pylab.savefig(saveimg_svg)
 
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()