Calculate radial distribution function
From Computational Biophysics and Materials Science Group
######################################### ## 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))