source: util/src/CalculatePairCorrelation.py.in@ aa8542

#! @PYTHON@
#
# determines all distance between each atomic pair and prints minimum and respective ids

import sys, string, re, math
wrerr = sys.stderr.write
wrout = sys.stdout.write
CELLLENGTH = 5

# check for parameters
if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <datafile> <Atom1> <Atom2>\n")
  sys.exit(1)

#Get box boundaries (open parameter file)
input = open(sys.argv[1], "r")
filetype=sys.argv[1].split(".")
line=input.readline()
line=input.readline()

line=input.readline()
atoms=line.split()
max=[float(atoms[1]), float(atoms[2]), float(atoms[3])]
min=[float(atoms[1]), float(atoms[2]), float(atoms[3])]
Listatoms=[]
Center=[0, 0, 0]
n=0
for line in input:
  if "#" in line:
    continue
  atom=line.split()
  if (atom[0] == sys.argv[2]) or (atom[0] == sys.argv[3]):
    Listatoms.append([str(atom[0]), float(atom[1]), float(atom[2]), float(atom[3])])
  for i in range(3):
    Center[i]+=float(atom[i+1])
    if float(atom[i+1]) < min[i]:
      min[i] = float(atom[i+1])
    if float(atom[i+1]) > max[i]:
      max[i] = float(atom[i+1])
print "Parsed "+str(len(Listatoms))+" atoms of type "+sys.argv[2]+" or "+sys.argv[3]+" parsed."
print "Found Box bounds [%f,%f]x[%f,%f]x[%f,%f]." % (min[0],max[0],min[1],max[1],min[2],max[2])

# find radius of cell
RMAX=0
for i in range(3):
  Center[i]=Center[i]/float(len(Listatoms))
  if RMAX < (max[i]-min[i]):
    RMAX=(max[i]-min[i])
RMAX=RMAX/2
R=(RMAX*RMAX)/4
print "Radius of set is "+str(RMAX)+" and of inner circle is "+str(math.sqrt(R))+"."

# find reference set of points
cell=[]
for i in range(len(Listatoms)):
  dist=0
  for j in range(3):
    dist+=(Listatoms[i][j+1]-Center[j])*(Listatoms[i][j+1]-Center[j])
  if dist < R:
    cell.append(i)
if len(cell) < 5:
  print "WARNING: We have found only "+str(len(cell))+" reference points."

# now, put each distance into bins
NumBins=len(Listatoms)
BinWidth=(RMAX*2)/float(NumBins)
Bins=[]
for i in range(NumBins+1):
  Bins.append(0)
print "I have created "+str(NumBins)+" bins with width "+str(BinWidth)+"."
# for each reference point
for ref in range(len(cell)):
  # for each other point in cluster
  for other in range(len(Listatoms)):
    if (cell[ref] != other):
      dist=0
      for j in range(3):
        dist+=(Listatoms[other][j+1]-Listatoms[ cell[ref] ][j+1])*(Listatoms[other][j+1]-Listatoms[ cell[ref] ][j+1])
      dist=math.sqrt(dist)
      #print "Distance between "+str( cell[ref] )+" and "+str(other)+" is "+str(dist)+", putting into Bin "+str(math.floor(dist/BinWidth))+"."
      Bins[int(math.floor(dist/BinWidth))]+=1

output=open("histogramm_"+sys.argv[2]+"-"+sys.argv[3]+".csv","w")
output.write("# BinNr\tDistStart\tDistCenter\tCount\n")
for i in range(1,len(Bins)):
  output.write("%d\t%lg\t%lg\t%lg\n" % (i,BinWidth*i, BinWidth*(i+1/2), Bins[i]/(4.*math.pi*BinWidth*float(i)*BinWidth*float(i))))