showPDF1d.py

# GPL v 3.0 
# script name: showPDF1d.py
# written by Naoyuki Miyashita on 12/18/2013 
# Naoyuki Miyashita (c) 2013
#
#-----------------------------------
# Use lProam, pylab, numpy
#-----------------------------------
# ussage:
# python showPDF1d.py \
# "end2end.1.d" 1 \  #filename, columun
# "Angstrom" \       # X axis
# "PDF" \            # Y axis
# 3 20 10 \    # x_min, x_max, #of bins
# "test.out" "pdf"  298  # output, "pdf" or "hist" or "pmf"
# in the case of "pmf", please write temperature at last.
#-----------------------------------
# Please use prody selection rules in the atom selections.
#
from prody import *
from pylab import *
import numpy as np
import lProam

def argv_read():
    inputfile     = sys.argv[1]
    nelement      = int(sys.argv[2])
    x_lab         = sys.argv[3]
    y_lab         = sys.argv[4]
    x_min         = float(sys.argv[5])
    x_max         = float(sys.argv[6])
    nx_bin        = int(sys.argv[7])
    outfile       = sys.argv[8]
    PDForHIST     = sys.argv[9]  # "pdf" or "hist" or "pmf"
    temp = 300
    if PDForHIST == "pmf":
        temp = float(sys.argv[10])

    argv_array=(inputfile,nelement,x_lab,y_lab,x_min,x_max,nx_bin,outfile,PDForHIST,temp)
    return(argv_array)

def main_calc():
    #input
    argv_array=argv_read()
    (inputfile,nelement,x_lab,y_lab,x_min,x_max,nx_bin,outfile,swch,temp)=argv_array
    #
    aa=np.array([])
    f=open(inputfile,'r')
    for line in f.readlines():
        elem=line.split()
        aa=np.append(aa,float(elem[nelement]))
    #print aa
    param_set=(aa,x_min,x_max,nx_bin,swch)
    (hist,bins)=lProam.dist1D(param_set)
    if swch == "pmf":
        (kboltz,a1,b1,c1,d1)=lProam.Constants()
        hist=-1.0*temp*kboltz*np.log(hist)
    mini=np.amin(hist)
    hist=hist-mini
    print hist
    param_set=(hist,bins,x_lab,y_lab)
    lProam.showPlot1(param_set)
    param_set=outfile
    lProam.outputPNG(param_set)
    return()


if __name__ == '__main__':
    main_calc()