source: util/src/CalculateXRayFromPDBstep.py.in@ a3ffb44

#!@PYTHON@
#
# calculates xray spectra for each time step with coordinates from PDB

import sys, os

exec_prefix="@prefix@"
DATABASE=exec_prefix+"/bin/elements.db" # contains atomic numbers
XRAY="/media/WORK/tmp/CODICE/XR_Intensity"
PYXPLOT="/home/heber/build/pyxplot/pyxplot"
NAME="ettringite_6Gpa"

# get command line parameters
if (len(sys.argv)) < 5:
        print "Usage: "+sys.argv[0]+" <dir> <prefix> <offsetElement> <offsetXYZ> [start from step]"
        sys.exit(1)
else:
        dir=sys.argv[1]
        prefix=sys.argv[2]
        offsetElement=int(sys.argv[3])
        offsetXYZ=int(sys.argv[4])

# generate hash map: element -> atomic number
atomic={}
input = open(DATABASE,"r")
for line in input:
        if "#" in line:
                continue
        entries=line.split()
        element=entries[1].upper()
        number=int(entries[5])
        print "Putting "+element+" -> "+str(number)+" into hash map."
        atomic[element]=number
input.close()

# parse each PDB
if (len(sys.argv))>=5:
  step=int(sys.argv[5])
else:
  step=0
file=dir+"/"+prefix+("%04d" % (step))+".pdb"
while 1:
        try:
                input=open(file, "r")
        except IOError:
                print file+" is not found, breaking."
                break
        print "Current step is "+str(step)+"."
        # recode pdb to "Symbol atomicnumber X Y Z" structure
        output = open(NAME,"w")
        n=0
        for line in input:
                if "#" in line or "CONECT" in line or "MASTER" in line or "AUTHOR" in line or "COMPND" in line:
                        continue
                if "END" in line:
                        break
                entries=line.split()
                try:
                        element=entries[offsetElement]
                except IndexError:
                        print "Could not read element nr."+str(offsetElement)+" in line: "+line+"."
                        continue
                number=atomic.get(element.upper(), None)
                if number != None:
                        output.write("%s\t%d\t%s\t%s\t%s\n" % (element, int(number), entries[offsetXYZ], entries[offsetXYZ+1], entries[offsetXYZ+2]))
                        n=n+1
                else:
                        print "WARNING: atomic number for element "+element+" is not known."
        output.close()
        # call Xray
        outputname=dir+"/xray_intensity"+("%04d" % (step))
        os.system(XRAY+' '+str(n)+' >'+outputname+'.csv')
        # make plotfile and move eps into place
        os.system('/bin/cp -f '+outputname+'.csv C3S_xray_spectrum.csv')
        os.system(PYXPLOT+' xray_spectrum.pyx')
        os.system('/bin/mv C3S_xray_spectrum.eps '+outputname+'.eps')
        # next file
        step+=1
        file=dir+"/"+prefix+("%04d" % (step))+".pdb"

sys.exit(1)