Calculate radial distribution function

#########################################
## Description: Python script for calculating rdf by reading lines from an xyz data
## Author: Kevin May 2014
## Usage: python RDF.py [traj.dat] [box_length (A)] [# of bins] [# of frames]
## Input: xyz
## Output: rdf
## Units: Angstron
## Other Notes: Note the unit and you may want to change to nm when plotting in Matlab
## Other Notes: NO pbc is considered
#########################################

import sys
import numpy as np
from operator import sub
import math

# Check Arguments
if len(sys.argv) == 5:
  trajfile      = sys.argv[1]
  boxlength     = float(sys.argv[2])
  n_bins        = int(sys.argv[3])
  NFRAME        = int(sys.argv[4])
else:
  print "Insufficent number of arguments. USAGE: $ python RDF.py [traj.dat] [box_length (A)] [# of bins] [# of frames]"
  sys.exit(1)

f = open(trajfile)
lines = f.readlines()
f.close()
print('This file has ' + str(len(lines)) + ' lines')
print('You are going to use ' + str(NFRAME) + ' frames')

N = 20
#NFRAME = 50000
#dr = 2.0
#n_bins = 1000
#boxlength = 300
print('I guess it is a ' + str(N) + '-atom system\nIf I am wrong, modify accordingly')

maxbin = n_bins
box_half = boxlength
dr = float(box_half/maxbin)
print('Your bin width will be ' + str(dr) + ' (A)')
hist = np.zeros((maxbin,maxbin))


print "Calculating g(r)..."
for frame in range(0, NFRAME):
  for start in range(0, N-1):
    compo1 = [float(x) for x in lines[frame * N * 2 + start].split()]
    compo2 = [float(x) for x in lines[frame * N * 2 + start + 1].split()]
    compo3 = [float(x) for x in lines[frame * N * 2 + start + N].split()]
    compo4 = [float(x) for x in lines[frame * N * 2 + start + N + 1].split()]

    # distance between i and j
    dd = map(sub, compo2, compo1)
    rij = np.linalg.norm(dd[0:3])

    bin = int(rij/dr)
    hist[start][bin] += 1

#-------------------------------------------------------------------#

# Write to file

for start in range(0, N-1):
  rdfout = './P'+repr(start+1)+'_P'+repr(start+2)+'_rdf.dat'
  print "Writing output to " + rdfout
  rdfout=open(rdfout, "w")
  for bin in range(0, maxbin):
    rdfout.write("%13.8g %13.8g\n"%(bin*dr, hist[start][bin] / NFRAME / 2.0))