[e1a46d] | 1 | #!@PYTHON@
|
---|
| 2 | #
|
---|
| 3 | # calculates xray spectra for each time step with coordinates from PDB
|
---|
| 4 |
|
---|
| 5 | import sys, os
|
---|
| 6 |
|
---|
| 7 | exec_prefix="@prefix@"
|
---|
| 8 | DATABASE=exec_prefix+"/bin/elements.db" # contains atomic numbers
|
---|
| 9 | XRAY="/media/WORK/tmp/CODICE/XR_Intensity"
|
---|
| 10 | PYXPLOT="/home/heber/build/pyxplot/pyxplot"
|
---|
| 11 | NAME="ettringite_6Gpa"
|
---|
| 12 |
|
---|
| 13 | # get command line parameters
|
---|
| 14 | if (len(sys.argv)) < 5:
|
---|
| 15 | print "Usage: "+sys.argv[0]+" <dir> <prefix> <offsetElement> <offsetXYZ> [start from step]"
|
---|
| 16 | sys.exit(1)
|
---|
| 17 | else:
|
---|
| 18 | dir=sys.argv[1]
|
---|
| 19 | prefix=sys.argv[2]
|
---|
| 20 | offsetElement=int(sys.argv[3])
|
---|
| 21 | offsetXYZ=int(sys.argv[4])
|
---|
| 22 |
|
---|
| 23 | # generate hash map: element -> atomic number
|
---|
| 24 | atomic={}
|
---|
| 25 | input = open(DATABASE,"r")
|
---|
| 26 | for line in input:
|
---|
| 27 | if "#" in line:
|
---|
| 28 | continue
|
---|
| 29 | entries=line.split()
|
---|
| 30 | element=entries[1].upper()
|
---|
| 31 | number=int(entries[5])
|
---|
| 32 | print "Putting "+element+" -> "+str(number)+" into hash map."
|
---|
| 33 | atomic[element]=number
|
---|
| 34 | input.close()
|
---|
| 35 |
|
---|
| 36 | # parse each PDB
|
---|
| 37 | if (len(sys.argv))>=5:
|
---|
| 38 | step=int(sys.argv[5])
|
---|
| 39 | else:
|
---|
| 40 | step=0
|
---|
| 41 | file=dir+"/"+prefix+("%04d" % (step))+".pdb"
|
---|
| 42 | while 1:
|
---|
| 43 | try:
|
---|
| 44 | input=open(file, "r")
|
---|
| 45 | except IOError:
|
---|
| 46 | print file+" is not found, breaking."
|
---|
| 47 | break
|
---|
| 48 | print "Current step is "+str(step)+"."
|
---|
| 49 | # recode pdb to "Symbol atomicnumber X Y Z" structure
|
---|
| 50 | output = open(NAME,"w")
|
---|
| 51 | n=0
|
---|
| 52 | for line in input:
|
---|
| 53 | if "#" in line or "CONECT" in line or "MASTER" in line or "AUTHOR" in line or "COMPND" in line:
|
---|
| 54 | continue
|
---|
| 55 | if "END" in line:
|
---|
| 56 | break
|
---|
| 57 | entries=line.split()
|
---|
| 58 | try:
|
---|
| 59 | element=entries[offsetElement]
|
---|
| 60 | except IndexError:
|
---|
| 61 | print "Could not read element nr."+str(offsetElement)+" in line: "+line+"."
|
---|
| 62 | continue
|
---|
| 63 | number=atomic.get(element.upper(), None)
|
---|
| 64 | if number != None:
|
---|
| 65 | output.write("%s\t%d\t%s\t%s\t%s\n" % (element, int(number), entries[offsetXYZ], entries[offsetXYZ+1], entries[offsetXYZ+2]))
|
---|
| 66 | n=n+1
|
---|
| 67 | else:
|
---|
| 68 | print "WARNING: atomic number for element "+element+" is not known."
|
---|
| 69 | output.close()
|
---|
| 70 | # call Xray
|
---|
| 71 | outputname=dir+"/xray_intensity"+("%04d" % (step))
|
---|
| 72 | os.system(XRAY+' '+str(n)+' >'+outputname+'.csv')
|
---|
| 73 | # make plotfile and move eps into place
|
---|
| 74 | os.system('/bin/cp -f '+outputname+'.csv C3S_xray_spectrum.csv')
|
---|
| 75 | os.system(PYXPLOT+' xray_spectrum.pyx')
|
---|
| 76 | os.system('/bin/mv C3S_xray_spectrum.eps '+outputname+'.eps')
|
---|
| 77 | # next file
|
---|
| 78 | step+=1
|
---|
| 79 | file=dir+"/"+prefix+("%04d" % (step))+".pdb"
|
---|
| 80 |
|
---|
| 81 | sys.exit(1)
|
---|