Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F66014983
lmp2radii.pyx
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Jun 7, 17:39
Size
2 KB
Mime Type
text/x-python
Expires
Sun, Jun 9, 17:39 (2 d)
Engine
blob
Format
Raw Data
Handle
18158528
Attached To
rLAMMPS lammps
lmp2radii.pyx
View Options
#!/usr/local/bin/python-2.5/bin/python
Info="""
Module name: lmp2radii.py
Author: (c) Andres Jaramillo-Botero
California Institute of Technology
ajaramil@wag.caltech.edu
Project: pEFF
Version: August 2009
Extracts the electron radii from a lammps trajectory dump of style custom:
dump 1 all custom period dump_file id type q spin eradius x y z...
NOTE: The radius must be the 5th column per trajectory entry in the dump file
"""
# import essentials:
import sys, os
from math import log10
from shutil import rmtree
from getopt import gnu_getopt as getopt
import numpy
def printHelp():
print Info
print "Usage: python lmp2radii.pyx test.lammpstrj\n"
return
def makeradii(infile):
print "Reading %s ... [WAIT]"%infile,
fin = open(infile,'r')
lines = fin.xreadlines()
print 7*"\b"+"[DONE]"
frame=0
radii=[]
# grep the number of frames and atoms/frame
os.system("grep TIMESTEP %s | wc -l > frames; grep -m 1 -A 1 ATOMS %s > atoms"%(infile,infile))
tmp=open("frames",'r')
frames=int(tmp.readline().split()[0])
tmp.close()
tmp=open("atoms",'r')
atoms=int(tmp.readlines()[1].split()[0])
tmp.close()
os.system("rm -rf frames atoms")
arry=numpy.zeros((atoms,frames),dtype=float)
framecnt=0
header=9
ecount=0
print "Extracting electron radii per frame from %s ... "%(infile),
for i,line in enumerate(lines):
lo=(atoms+header)*framecnt+header
hi=lo+atoms
if (i<lo):
continue
elif (i >= lo) and (i < hi):
lparse=line.split()
id=int(lparse[0])
r=float(lparse[4])
if (r!=0):
arry[id-1][framecnt]=r
if (framecnt==0): ecount+=1
if (i==lo+1):
sys.stdout.write("%d/%d%s"%(framecnt+1,frames,(int(log10(framecnt+1))+3+int(log10(frames)))*"\b"))
sys.stdout.flush()
if (i == hi+1):
framecnt+=1
print
print "Writing radii/frame table to %s ... "%(infile+'.out'),
sys.stdout.flush()
fout=open(infile+'.out','w')
for i in range(frames):
fout.writelines('\tF'+str(i))
fout.writelines("\n")
e=1
for a in range(atoms):
if arry[a][0] == 0.0: continue
else:
sys.stdout.write("%d/%d%s"%(e,ecount,(int(log10(e))+int(log10(ecount))+3)*"\b"))
sys.stdout.flush()
e+=1
fout.writelines("%d\t"%(a+1))
for f in range(frames):
fout.writelines("%f\t"%(arry[a][f]))
fout.writelines("\n")
print
print "Done !! (generated radii/frame table) \n"
fout.close()
fin.close()
if __name__ == '__main__':
# set defaults
# check for input:
opts, argv = getopt(sys.argv[1:], 'h')
# if no input, print help and exit
if len(argv) != 1:
printHelp()
sys.exit(1)
else:
infile=argv[0]
# read options
for opt, arg in opts:
if opt == '-h': # -h: print help
printHelp()
makeradii(infile)
Event Timeline
Log In to Comment