1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
#!/usr/local/bin/python-2.5/bin/python
Info="""
Module name: lmp2data.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 x y z spin radius ...
NOTE: The radius must be the i'th 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 lmp2data.py test.lammpstrj\n"
return
def makeradii(infile,outfile,column,flag_all):
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; grep -m 1 \"ITEM: ATOMS\" %s > params"%(infile,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()
tmp=open("params",'r')
ids=tmp.readline().split()[2:]
os.system("rm -rf frames atoms params")
arry=numpy.zeros((atoms,frames),dtype=str)
framecnt=0
header=9
ecount=0
if flag_all==True: atom_type="nuclei and electron"
else: atom_type="electron"
print "Extracting %s %s per frame from %s ... "%(atom_type,ids[column],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[column-1])
r=lparse[column]
# if (float(r)!=0):
arry[id-1][framecnt]=r
print arry[id-1][framecnt],r,raw_input()
if (float(r)!=0) and (framecnt==0): ecount+=1
# else: arry[id-1][framecnt]=r
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
if outfile=="":
outfile=infile+'.%s'%(ids[column])
fout=open(outfile,'w')
else: fout=open(outfile,'w')
print "Writing %s/frame table to %s ... "%(ids[column],outfile),
sys.stdout.flush()
for i in range(frames):
fout.writelines('\tF'+str(i))
fout.writelines("\n")
e=1
for a in range(atoms):
if flag_all==True:
sys.stdout.write("%d/%d%s"%(a+1,atoms,(int(log10(a+1))+int(log10(atoms))+3)*"\b"))
sys.stdout.flush()
fout.writelines("%d\t"%(a+1))
for f in range(frames):
fout.writelines("%s\t"%(arry[a][f]))
fout.writelines("\n")
else:
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("%s\t"%(arry[a][f]))
fout.writelines("\n")
print
print "DONE .... GOODBYE !!"
fout.close()
fin.close()
if __name__ == '__main__':
# set defaults
outfile = ""
flag_all = False
column=6 # default = radius
# check for input:
opts, argv = getopt(sys.argv[1:], 'c:o:ha')
# 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()
if opt == '-o': # output file name
outfile=arg
if opt == '-a': # all nuclii+electrons
flag_all=True
if opt == '-c': # select column from lammpstrj file to tabulate
column=int(arg)
makeradii(infile,outfile,column,flag_all)
|
