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

#! @PYTHON@
#
# calculates the density via elements.db

exec_prefix="@prefix@"
ELEMENTSDB=exec_prefix+"/bin/elements.db"
import sys, string, re, math
wrerr = sys.stderr.write
wrout = sys.stdout.write

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

offset=int(sys.argv[2])
offsetXYZ=int(sys.argv[3])


# parse elements symbol and mass into hash list
print "Parsing elements db at "+ELEMENTSDB
elementfile=open(ELEMENTSDB,"r")
elements={}
for line in elementfile:
  if "#" in line:
    continue
  entries=line.split()
  if len(entries) > 6:
    elements[entries[1]]=float(entries[6]) 
  else:
    print "Not parsed: "+line
elementfile.close()

# go through given data file and sum up
inputfile=open(sys.argv[1],"r")
weight=0.
min=[0., 0., 0.]
max=[0., 0., 0.]
for line in inputfile:
  if "#" in line:
    continue
  entries=line.split()
  if weight == 0.:
    for i in range(3):
      min[i]=float(entries[offsetXYZ+i])
      max[i]=float(entries[offsetXYZ+i])
  else:
    for i in range(3):
      if (min[i] > float(entries[offsetXYZ+i])):
        min[i]=float(entries[offsetXYZ+i])
      if (max[i] < float(entries[offsetXYZ+i])):
        max[i]=float(entries[offsetXYZ+i])
  weight += elements[entries[offset]]
inputfile.close()

volume=1.
wrout("Bounds are ")
for i in range(3):
  wrout("[%g, %g]\t" % (min[i], max[i]))
  volume=volume*(max[i]-min[i])
wrout("\n")
print "Total Volume of bounding box is "+str(volume)+" angstrom^3."
print "Total weight is "+str(weight)+" atomicmassunit."
print "Total density is "+str(weight/volume)+" atomicmassunit/angstrom^3."
print "Total density is "+str(1.6605387*weight/volume)+" g/cm^3."