diff --git a/tools/python/pizza/dump.py b/tools/python/pizza/dump.py
index ed0f430f4..153ab1a18 100644
--- a/tools/python/pizza/dump.py
+++ b/tools/python/pizza/dump.py
@@ -1,1168 +1,1189 @@
# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
#
# Copyright (2005) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
# dump tool
oneline = "Read, write, manipulate dump files and particle attributes"
docstr = """
d = dump("dump.one") read in one or more dump files
d = dump("dump.1 dump.2.gz") can be gzipped
d = dump("dump.*") wildcard expands to multiple files
d = dump("dump.*",0) two args = store filenames, but don't read
incomplete and duplicate snapshots are deleted
if atoms have 5 or 8 columns, assign id,type,x,y,z (ix,iy,iz)
atoms will be unscaled if stored in files as scaled
time = d.next() read next snapshot from dump files
used with 2-argument constructor to allow reading snapshots one-at-a-time
snapshot will be skipped only if another snapshot has same time stamp
return time stamp of snapshot read
return -1 if no snapshots left or last snapshot is incomplete
no column name assignment or unscaling is performed
d.map(1,"id",3,"x") assign names to atom columns (1-N)
not needed if dump file is self-describing
d.tselect.all() select all timesteps
d.tselect.one(N) select only timestep N
d.tselect.none() deselect all timesteps
d.tselect.skip(M) select every Mth step
d.tselect.test("$t >= 100 and $t < 10000") select matching timesteps
d.delete() delete non-selected timesteps
selecting a timestep also selects all atoms in the timestep
skip() and test() only select from currently selected timesteps
test() uses a Python Boolean expression with $t for timestep value
Python comparison syntax: == != < > <= >= and or
d.aselect.all() select all atoms in all steps
d.aselect.all(N) select all atoms in one step
d.aselect.test("$id > 100 and $type == 2") select match atoms in all steps
d.aselect.test("$id > 100 and $type == 2",N) select matching atoms in one step
all() with no args selects atoms from currently selected timesteps
test() with one arg selects atoms from currently selected timesteps
test() sub-selects from currently selected atoms
test() uses a Python Boolean expression with $ for atom attributes
Python comparison syntax: == != < > <= >= and or
$name must end with a space
d.write("file") write selected steps/atoms to dump file
d.scatter("tmp") write selected steps/atoms to mutiple files
scatter() files are given timestep suffix: e.g. tmp.0, tmp.100, etc
d.scale() scale x,y,z to 0-1 for all timesteps
d.scale(100) scale atom coords for timestep N
d.unscale() unscale x,y,z to box size to all timesteps
d.unscale(1000) unscale atom coords for timestep N
d.wrap() wrap x,y,z into periodic box via ix,iy,iz
d.unwrap() unwrap x,y,z out of box via ix,iy,iz
d.owrap("other") wrap x,y,z to same image as another atom
d.sort() sort atoms by atom ID in all selected steps
d.sort("x") sort atoms by column value in all steps
d.sort(1000) sort atoms in timestep N
scale(), unscale(), wrap(), unwrap(), owrap() operate on all steps and atoms
wrap(), unwrap(), owrap() require ix,iy,iz be defined
owrap() requires a column be defined which contains an atom ID
name of that column is the argument to owrap()
x,y,z for each atom is wrapped to same image as the associated atom ID
useful for wrapping all molecule's atoms the same so it is contiguous
m1,m2 = d.minmax("type") find min/max values for a column
d.set("$ke = $vx * $vx + $vy * $vy") set a column to a computed value
d.spread("ke",N,"color") 2nd col = N ints spread over 1st col
d.clone(1000,"color") clone timestep N values to other steps
minmax() operates on selected timesteps and atoms
set() operates on selected timesteps and atoms
left hand side column is created if necessary
left-hand side column is unset or unchanged for non-selected atoms
equation is in Python syntax
use $ for column names, $name must end with a space
spread() operates on selected timesteps and atoms
min and max are found for 1st specified column across all selected atoms
atom's value is linear mapping (1-N) between min and max
that is stored in 2nd column (created if needed)
useful for creating a color map
clone() operates on selected timesteps and atoms
values at every timestep are set to value at timestep N for that atom ID
useful for propagating a color map
t = d.time() return vector of selected timestep values
fx,fy,... = d.atom(100,"fx","fy",...) return vector(s) for atom ID N
fx,fy,... = d.vecs(1000,"fx","fy",...) return vector(s) for timestep N
atom() returns vectors with one value for each selected timestep
vecs() returns vectors with one value for each selected atom in the timestep
index,time,flag = d.iterator(0/1) loop over dump snapshots
time,box,atoms,bonds,tris = d.viz(index) return list of viz objects
d.atype = "color" set column returned as "type" by viz
d.extra("dump.bond") read bond list from dump file
d.extra(data) extract bond/tri/line list from data
iterator() loops over selected timesteps
iterator() called with arg = 0 first time, with arg = 1 on subsequent calls
index = index within dump object (0 to # of snapshots)
time = timestep value
flag = -1 when iteration is done, 1 otherwise
viz() returns info for selected atoms for specified timestep index
time = timestep value
box = [xlo,ylo,zlo,xhi,yhi,zhi]
atoms = id,type,x,y,z for each atom as 2d array
bonds = id,type,x1,y1,z1,x2,y2,z2,t1,t2 for each bond as 2d array
if bonds() was used to define bonds, else empty list
tris = id,type,x1,y1,z1,x2,y2,z2,x3,y3,z3,nx,ny,nz for each tri as 2d array
if extra() was used to define tris, else empty list
lines = id,type,x1,y1,z1,x2,y2,z2 for each line as 2d array
if extra() was used to define lines, else empty list
atype is column name viz() will return as atom type (def = "type")
extra() stores list of bonds/tris/lines to return each time viz() is called
"""
# History
# 8/05, Steve Plimpton (SNL): original version
# ToDo list
# try to optimize this line in read_snap: words += f.readline().split()
# allow $name in aselect.test() and set() to end with non-space
# should next() snapshot be auto-unscaled ?
# Variables
# flist = list of dump file names
# increment = 1 if reading snapshots one-at-a-time
# nextfile = which file to read from via next()
# eof = ptr into current file for where to read via next()
# nsnaps = # of snapshots
# nselect = # of selected snapshots
# snaps = list of snapshots
# names = dictionary of column names:
# key = "id", value = column # (0 to M-1)
# tselect = class for time selection
# aselect = class for atom selection
# atype = name of vector used as atom type by viz extract
# bondflag = 0 if no bonds, 1 if they are defined statically
# bondlist = static list of bonds to viz() return for all snapshots
# only a list of atom pairs, coords have to be created for each snapshot
# triflag = 0 if no tris, 1 if they are defined statically, 2 if dynamic
# trilist = static list of tris to return via viz() for all snapshots
# lineflag = 0 if no lines, 1 if they are defined statically
# linelist = static list of lines to return via viz() for all snapshots
# Snap = one snapshot
# time = time stamp
# tselect = 0/1 if this snapshot selected
# natoms = # of atoms
# nselect = # of selected atoms in this snapshot
# aselect[i] = 0/1 for each atom
# xlo,xhi,ylo,yhi,zlo,zhi = box bounds (float)
# atoms[i][j] = 2d array of floats, i = 0 to natoms-1, j = 0 to ncols-1
# Imports and external programs
import sys, commands, re, glob, types
from os import popen
from math import * # any function could be used by set()
import Numeric
try: from DEFAULTS import PIZZA_GUNZIP
except: PIZZA_GUNZIP = "gunzip"
# Class definition
class dump:
# --------------------------------------------------------------------
def __init__(self,*list):
self.snaps = []
self.nsnaps = self.nselect = 0
self.names = {}
self.tselect = tselect(self)
self.aselect = aselect(self)
self.atype = "type"
self.bondflag = 0
self.bondlist = []
self.triflag = 0
self.trilist = []
self.triobj = 0
self.lineflag = 0
self.linelist = []
# flist = list of all dump file names
words = list[0].split()
self.flist = []
for word in words: self.flist += glob.glob(word)
if len(self.flist) == 0 and len(list) == 1:
raise StandardError,"no dump file specified"
if len(list) == 1:
self.increment = 0
self.read_all()
else:
self.increment = 1
self.nextfile = 0
self.eof = 0
# --------------------------------------------------------------------
def read_all(self):
# read all snapshots from each file
# test for gzipped files
for file in self.flist:
if file[-3:] == ".gz":
f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r')
else: f = open(file)
snap = self.read_snapshot(f)
while snap:
self.snaps.append(snap)
print snap.time,
sys.stdout.flush()
snap = self.read_snapshot(f)
f.close()
print
# sort entries by timestep, cull duplicates
self.snaps.sort(self.compare_time)
self.cull()
self.nsnaps = len(self.snaps)
print "read %d snapshots" % self.nsnaps
# select all timesteps and atoms
self.tselect.all()
# set default names for atom columns if file wasn't self-describing
if len(self.snaps) == 0:
print "no column assignments made"
elif len(self.names):
- print "column assignments made from self-describing file"
+ print "assigned columns:",self.names2str()
elif self.snaps[0].atoms == None:
print "no column assignments made"
elif len(self.snaps[0].atoms[0]) == 5:
self.map(1,"id",2,"type",3,"x",4,"y",5,"z")
- print "assigned columns id,type,x,y,z"
+ print "assigned columns:",self.names2str()
elif len(self.snaps[0].atoms[0]) == 8:
self.map(1,"id",2,"type",3,"x",4,"y",5,"z",6,"ix",7,"iy",8,"iz")
- print "assigned columns id,type,x,y,z,ix,iy,iz"
+ print "assigned columns:",self.names2str()
else:
print "no column assignments made"
# if snapshots are scaled, unscale them
if (not self.names.has_key("x")) or \
(not self.names.has_key("y")) or \
(not self.names.has_key("z")):
print "no unscaling could be performed"
elif self.nsnaps > 0:
if self.scaled(self.nsnaps-1): self.unscale()
else: print "dump is already unscaled"
# --------------------------------------------------------------------
# read next snapshot from list of files
def next(self):
if not self.increment: raise StandardError,"cannot read incrementally"
# read next snapshot in current file using eof as pointer
# if fail, try next file
# if new snapshot time stamp already exists, read next snapshot
while 1:
- f = open(self.flist[self.nextfile],'r')
+ f = open(self.flist[self.nextfile],'rb')
f.seek(self.eof)
snap = self.read_snapshot(f)
if not snap:
self.nextfile += 1
if self.nextfile == len(self.flist): return -1
f.close()
self.eof = 0
continue
self.eof = f.tell()
f.close()
try:
self.findtime(snap.time)
continue
except: break
# select the new snapshot with all its atoms
self.snaps.append(snap)
snap = self.snaps[self.nsnaps]
snap.tselect = 1
snap.nselect = snap.natoms
for i in xrange(snap.natoms): snap.aselect[i] = 1
self.nsnaps += 1
self.nselect += 1
return snap.time
# --------------------------------------------------------------------
# read a single snapshot from file f
# return snapshot or 0 if failed
# assign column names if not already done and file is self-describing
+ # convert xs,xu to x
def read_snapshot(self,f):
try:
snap = Snap()
item = f.readline()
snap.time = int(f.readline().split()[0]) # just grab 1st field
item = f.readline()
snap.natoms = int(f.readline())
snap.aselect = Numeric.zeros(snap.natoms)
item = f.readline()
words = f.readline().split()
snap.xlo,snap.xhi = float(words[0]),float(words[1])
words = f.readline().split()
snap.ylo,snap.yhi = float(words[0]),float(words[1])
words = f.readline().split()
snap.zlo,snap.zhi = float(words[0]),float(words[1])
item = f.readline()
if len(self.names) == 0:
words = item.split()[2:]
if len(words):
for i in range(len(words)):
- self.names[words[i]] = i
+ if words[i] == "xs" or words[i] == "xu":
+ self.names["x"] = i
+ elif words[i] == "ys" or words[i] == "yu":
+ self.names["y"] = i
+ elif words[i] == "zs" or words[i] == "zu":
+ self.names["z"] = i
+ else: self.names[words[i]] = i
if snap.natoms:
words = f.readline().split()
ncol = len(words)
for i in xrange(1,snap.natoms):
words += f.readline().split()
floats = map(float,words)
atoms = Numeric.zeros((snap.natoms,ncol),Numeric.Float)
start = 0
stop = ncol
for i in xrange(snap.natoms):
atoms[i] = floats[start:stop]
start = stop
stop += ncol
else: atoms = None
snap.atoms = atoms
return snap
except:
return 0
# --------------------------------------------------------------------
# decide if snapshot i is scaled/unscaled from coords of first and last atom
def scaled(self,i):
ix = self.names["x"]
iy = self.names["y"]
iz = self.names["z"]
natoms = self.snaps[i].natoms
if natoms == 0: return 0
x1 = self.snaps[i].atoms[0][ix]
y1 = self.snaps[i].atoms[0][iy]
z1 = self.snaps[i].atoms[0][iz]
x2 = self.snaps[i].atoms[natoms-1][ix]
y2 = self.snaps[i].atoms[natoms-1][iy]
z2 = self.snaps[i].atoms[natoms-1][iz]
if x1 >= -0.1 and x1 <= 1.1 and y1 >= -0.1 and y1 <= 1.1 and \
z1 >= -0.1 and z1 <= 1.1 and x2 >= -0.1 and x2 <= 1.1 and \
y2 >= -0.1 and y2 <= 1.1 and z2 >= -0.1 and z2 <= 1.1:
return 1
else: return 0
# --------------------------------------------------------------------
# map atom column names
def map(self,*pairs):
if len(pairs) % 2 != 0:
raise StandardError, "dump map() requires pairs of mappings"
for i in range(0,len(pairs),2):
j = i + 1
self.names[pairs[j]] = pairs[i]-1
# delete unselected snapshots
# --------------------------------------------------------------------
def delete(self):
ndel = i = 0
while i < self.nsnaps:
if not self.snaps[i].tselect:
del self.snaps[i]
self.nsnaps -= 1
ndel += 1
else: i += 1
print "%d snapshots deleted" % ndel
print "%d snapshots remaining" % self.nsnaps
# --------------------------------------------------------------------
# scale coords to 0-1 for all snapshots or just one
def scale(self,*list):
if len(list) == 0:
print "Scaling dump ..."
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
for snap in self.snaps: self.scale_one(snap,x,y,z)
else:
i = self.findtime(list[0])
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
self.scale_one(self.snaps[i],x,y,z)
# --------------------------------------------------------------------
def scale_one(self,snap,x,y,z):
xprdinv = 1.0 / (snap.xhi - snap.xlo)
yprdinv = 1.0 / (snap.yhi - snap.ylo)
zprdinv = 1.0 / (snap.zhi - snap.zlo)
atoms = snap.atoms
atoms[:,x] = (atoms[:,x] - snap.xlo) * xprdinv
atoms[:,y] = (atoms[:,y] - snap.ylo) * yprdinv
atoms[:,z] = (atoms[:,z] - snap.zlo) * zprdinv
# --------------------------------------------------------------------
# unscale coords from 0-1 to box size for all snapshots or just one
def unscale(self,*list):
if len(list) == 0:
print "Unscaling dump ..."
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
for snap in self.snaps: self.unscale_one(snap,x,y,z)
else:
i = self.findtime(list[0])
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
self.unscale_one(self.snaps[i],x,y,z)
# --------------------------------------------------------------------
def unscale_one(self,snap,x,y,z):
xprd = snap.xhi - snap.xlo
yprd = snap.yhi - snap.ylo
zprd = snap.zhi - snap.zlo
atoms = snap.atoms
atoms[:,x] = snap.xlo + atoms[:,x]*xprd
atoms[:,y] = snap.ylo + atoms[:,y]*yprd
atoms[:,z] = snap.zlo + atoms[:,z]*zprd
# --------------------------------------------------------------------
# wrap coords from outside box to inside
def wrap(self):
print "Wrapping dump ..."
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
ix = self.names["ix"]
iy = self.names["iy"]
iz = self.names["iz"]
for snap in self.snaps:
xprd = snap.xhi - snap.xlo
yprd = snap.yhi - snap.ylo
zprd = snap.zhi - snap.zlo
atoms = snap.atoms
atoms[:,x] -= atoms[:,ix]*xprd
atoms[:,y] -= atoms[:,iy]*yprd
atoms[:,z] -= atoms[:,iz]*zprd
# --------------------------------------------------------------------
# unwrap coords from inside box to outside
def unwrap(self):
print "Unwrapping dump ..."
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
ix = self.names["ix"]
iy = self.names["iy"]
iz = self.names["iz"]
for snap in self.snaps:
xprd = snap.xhi - snap.xlo
yprd = snap.yhi - snap.ylo
zprd = snap.zhi - snap.zlo
atoms = snap.atoms
atoms[:,x] += atoms[:,ix]*xprd
atoms[:,y] += atoms[:,iy]*yprd
atoms[:,z] += atoms[:,iz]*zprd
# --------------------------------------------------------------------
# wrap coords to same image as atom ID stored in "other" column
def owrap(self,other):
print "Wrapping to other ..."
id = self.names["id"]
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
ix = self.names["ix"]
iy = self.names["iy"]
iz = self.names["iz"]
iother = self.names[other]
for snap in self.snaps:
xprd = snap.xhi - snap.xlo
yprd = snap.yhi - snap.ylo
zprd = snap.zhi - snap.zlo
atoms = snap.atoms
ids = {}
for i in xrange(snap.natoms):
ids[atoms[i][id]] = i
for i in xrange(snap.natoms):
j = ids[atoms[i][iother]]
atoms[i][x] += (atoms[i][ix]-atoms[j][ix])*xprd
atoms[i][y] += (atoms[i][iy]-atoms[j][iy])*yprd
atoms[i][z] += (atoms[i][iz]-atoms[j][iz])*zprd
+ # --------------------------------------------------------------------
+ # convert column names assignment to a string, in column order
+
+ def names2str(self):
+ ncol = len(self.snaps[0].atoms[0])
+ pairs = self.names.items()
+ values = self.names.values()
+ str = ""
+ for i in xrange(ncol):
+ if i in values: str += pairs[values.index(i)][0] + ' '
+ return str
+
# --------------------------------------------------------------------
# sort atoms by atom ID in all selected timesteps by default
# if arg = string, sort all steps by that column
# if arg = numeric, sort atoms in single step
def sort(self,*list):
if len(list) == 0:
print "Sorting selected snapshots ..."
id = self.names["id"]
for snap in self.snaps:
if snap.tselect: self.sort_one(snap,id)
elif type(list[0]) is types.StringType:
print "Sorting selected snapshots by %s ..." % list[0]
id = self.names[list[0]]
for snap in self.snaps:
if snap.tselect: self.sort_one(snap,id)
else:
i = self.findtime(list[0])
id = self.names["id"]
self.sort_one(self.snaps[i],id)
# --------------------------------------------------------------------
# sort a single snapshot by ID column
def sort_one(self,snap,id):
atoms = snap.atoms
ids = atoms[:,id]
ordering = Numeric.argsort(ids)
for i in xrange(len(atoms[0])):
atoms[:,i] = Numeric.take(atoms[:,i],ordering)
# --------------------------------------------------------------------
# write a single dump file from current selection
def write(self,file):
+ if len(self.snaps): namestr = self.names2str()
f = open(file,"w")
for snap in self.snaps:
if not snap.tselect: continue
print snap.time,
sys.stdout.flush()
print >>f,"ITEM: TIMESTEP"
print >>f,snap.time
print >>f,"ITEM: NUMBER OF ATOMS"
print >>f,snap.nselect
print >>f,"ITEM: BOX BOUNDS"
print >>f,snap.xlo,snap.xhi
print >>f,snap.ylo,snap.yhi
print >>f,snap.zlo,snap.zhi
- print >>f,"ITEM: ATOMS"
+ print >>f,"ITEM: ATOMS",namestr
atoms = snap.atoms
nvalues = len(atoms[0])
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
line = ""
for j in xrange(nvalues):
if (j < 2):
line += str(int(atoms[i][j])) + " "
else:
line += str(atoms[i][j]) + " "
print >>f,line
f.close()
print "\n%d snapshots" % self.nselect
# --------------------------------------------------------------------
# write one dump file per snapshot from current selection
def scatter(self,root):
+ if len(self.snaps): namestr = self.names2str()
for snap in self.snaps:
if not snap.tselect: continue
print snap.time,
sys.stdout.flush()
file = root + "." + str(snap.time)
f = open(file,"w")
print >>f,"ITEM: TIMESTEP"
print >>f,snap.time
print >>f,"ITEM: NUMBER OF ATOMS"
print >>f,snap.nselect
print >>f,"ITEM: BOX BOUNDS"
print >>f,snap.xlo,snap.xhi
print >>f,snap.ylo,snap.yhi
print >>f,snap.zlo,snap.zhi
- print >>f,"ITEM: ATOMS"
+ print >>f,"ITEM: ATOMS",namestr
atoms = snap.atoms
nvalues = len(atoms[0])
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
line = ""
for j in xrange(nvalues):
if (j < 2):
line += str(int(atoms[i][j])) + " "
else:
line += str(atoms[i][j]) + " "
print >>f,line
f.close()
print "\n%d snapshots" % self.nselect
# --------------------------------------------------------------------
# find min/max across all selected snapshots/atoms for a particular column
def minmax(self,colname):
icol = self.names[colname]
min = 1.0e20
max = -min
for snap in self.snaps:
if not snap.tselect: continue
atoms = snap.atoms
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
if atoms[i][icol] < min: min = atoms[i][icol]
if atoms[i][icol] > max: max = atoms[i][icol]
return (min,max)
# --------------------------------------------------------------------
# set a column value via an equation
def set(self,eq):
print "Setting ..."
pattern = "\$\w*"
list = re.findall(pattern,eq)
lhs = list[0][1:]
if not self.names.has_key(lhs):
self.newcolumn(lhs)
for item in list:
name = item[1:]
column = self.names[name]
insert = "snap.atoms[i][%d]" % (column)
eq = eq.replace(item,insert)
ceq = compile(eq,'','single')
for snap in self.snaps:
if not snap.tselect: continue
for i in xrange(snap.natoms):
if snap.aselect[i]: exec ceq
# --------------------------------------------------------------------
# clone value in col across selected timesteps for atoms with same ID
def clone(self,nstep,col):
istep = self.findtime(nstep)
icol = self.names[col]
id = self.names["id"]
ids = {}
for i in xrange(self.snaps[istep].natoms):
ids[self.snaps[istep].atoms[i][id]] = i
for snap in self.snaps:
if not snap.tselect: continue
atoms = snap.atoms
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
j = ids[atoms[i][id]]
atoms[i][icol] = self.snaps[istep].atoms[j][icol]
# --------------------------------------------------------------------
# values in old column are spread as ints from 1-N and assigned to new column
def spread(self,old,n,new):
iold = self.names[old]
if not self.names.has_key(new): self.newcolumn(new)
inew = self.names[new]
min,max = self.minmax(old)
print "min/max = ",min,max
gap = max - min
invdelta = n/gap
for snap in self.snaps:
if not snap.tselect: continue
atoms = snap.atoms
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
ivalue = int((atoms[i][iold] - min) * invdelta) + 1
if ivalue > n: ivalue = n
if ivalue < 1: ivalue = 1
atoms[i][inew] = ivalue
# --------------------------------------------------------------------
# return vector of selected snapshot time stamps
def time(self):
vec = self.nselect * [0]
i = 0
for snap in self.snaps:
if not snap.tselect: continue
vec[i] = snap.time
i += 1
return vec
# --------------------------------------------------------------------
# extract vector(s) of values for atom ID n at each selected timestep
def atom(self,n,*list):
if len(list) == 0:
raise StandardError, "no columns specified"
columns = []
values = []
for name in list:
columns.append(self.names[name])
values.append(self.nselect * [0])
ncol = len(columns)
id = self.names["id"]
m = 0
for snap in self.snaps:
if not snap.tselect: continue
atoms = snap.atoms
for i in xrange(snap.natoms):
if atoms[i][id] == n: break
if atoms[i][id] != n:
raise StandardError, "could not find atom ID in snapshot"
for j in xrange(ncol):
values[j][m] = atoms[i][columns[j]]
m += 1
if len(list) == 1: return values[0]
else: return values
# --------------------------------------------------------------------
# extract vector(s) of values for selected atoms at chosen timestep
def vecs(self,n,*list):
snap = self.snaps[self.findtime(n)]
if len(list) == 0:
raise StandardError, "no columns specified"
columns = []
values = []
for name in list:
columns.append(self.names[name])
values.append(snap.nselect * [0])
ncol = len(columns)
m = 0
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
for j in xrange(ncol):
values[j][m] = snap.atoms[i][columns[j]]
m += 1
if len(list) == 1: return values[0]
else: return values
# --------------------------------------------------------------------
# add a new column to every snapshot and set value to 0
# set the name of the column to str
def newcolumn(self,str):
ncol = len(self.snaps[0].atoms[0])
self.map(ncol+1,str)
for snap in self.snaps:
atoms = snap.atoms
newatoms = Numeric.zeros((snap.natoms,ncol+1),Numeric.Float)
newatoms[:,0:ncol] = snap.atoms
snap.atoms = newatoms
# --------------------------------------------------------------------
# sort snapshots on time stamp
def compare_time(self,a,b):
if a.time < b.time:
return -1
elif a.time > b.time:
return 1
else:
return 0
# --------------------------------------------------------------------
# delete successive snapshots with duplicate time stamp
def cull(self):
i = 1
while i < len(self.snaps):
if self.snaps[i].time == self.snaps[i-1].time:
del self.snaps[i]
else:
i += 1
# --------------------------------------------------------------------
# iterate over selected snapshots
def iterator(self,flag):
start = 0
if flag: start = self.iterate + 1
for i in xrange(start,self.nsnaps):
if self.snaps[i].tselect:
self.iterate = i
return i,self.snaps[i].time,1
return 0,0,-1
# --------------------------------------------------------------------
# return list of atoms to viz for snapshot isnap
# augment with bonds, tris, lines if extra() was invoked
def viz(self,isnap):
snap = self.snaps[isnap]
time = snap.time
box = [snap.xlo,snap.ylo,snap.zlo,snap.xhi,snap.yhi,snap.zhi]
id = self.names["id"]
type = self.names[self.atype]
x = self.names["x"]
y = self.names["y"]
z = self.names["z"]
# create atom list needed by viz from id,type,x,y,z
# need Numeric mode here
atoms = []
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
atom = snap.atoms[i]
atoms.append([atom[id],atom[type],atom[x],atom[y],atom[z]])
# create list of current bond coords from static bondlist
# alist = dictionary of atom IDs for atoms list
# lookup bond atom IDs in alist and grab their coords
# try is used since some atoms may be unselected
# any bond with unselected atom is not returned to viz caller
# need Numeric mode here
bonds = []
if self.bondflag:
alist = {}
for i in xrange(len(atoms)): alist[int(atoms[i][0])] = i
for bond in self.bondlist:
try:
i = alist[bond[2]]
j = alist[bond[3]]
atom1 = atoms[i]
atom2 = atoms[j]
bonds.append([bond[0],bond[1],atom1[2],atom1[3],atom1[4],
atom2[2],atom2[3],atom2[4],atom1[1],atom2[1]])
except: continue
tris = []
if self.triflag:
if self.triflag == 1: tris = self.trilist
elif self.triflag == 2:
timetmp,boxtmp,atomstmp,bondstmp, \
tris,linestmp = self.triobj.viz(time,1)
lines = []
if self.lineflag: lines = self.linelist
return time,box,atoms,bonds,tris,lines
# --------------------------------------------------------------------
def findtime(self,n):
for i in xrange(self.nsnaps):
if self.snaps[i].time == n: return i
raise StandardError, "no step %d exists" % n
# --------------------------------------------------------------------
# return maximum box size across all selected snapshots
def maxbox(self):
xlo = ylo = zlo = None
xhi = yhi = zhi = None
for snap in self.snaps:
if not snap.tselect: continue
if xlo == None or snap.xlo < xlo: xlo = snap.xlo
if xhi == None or snap.xhi > xhi: xhi = snap.xhi
if ylo == None or snap.ylo < ylo: ylo = snap.ylo
if yhi == None or snap.yhi > yhi: yhi = snap.yhi
if zlo == None or snap.zlo < zlo: zlo = snap.zlo
if zhi == None or snap.zhi > zhi: zhi = snap.zhi
return [xlo,ylo,zlo,xhi,yhi,zhi]
# --------------------------------------------------------------------
# return maximum atom type across all selected snapshots and atoms
def maxtype(self):
icol = self.names["type"]
max = 0
for snap in self.snaps:
if not snap.tselect: continue
atoms = snap.atoms
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
if atoms[i][icol] > max: max = atoms[i][icol]
return int(max)
# --------------------------------------------------------------------
# grab bonds/tris/lines from another object
def extra(self,arg):
# read bonds from bond dump file
if type(arg) is types.StringType:
try:
f = open(arg,'r')
item = f.readline()
time = int(f.readline())
item = f.readline()
nbonds = int(f.readline())
item = f.readline()
if not re.search("BONDS",item):
raise StandardError, "could not read bonds from dump file"
words = f.readline().split()
ncol = len(words)
for i in xrange(1,nbonds):
words += f.readline().split()
f.close()
# convert values to int and absolute value since can be negative types
bondlist = Numeric.zeros((nbonds,4),Numeric.Int)
ints = [abs(int(value)) for value in words]
start = 0
stop = 4
for i in xrange(nbonds):
bondlist[i] = ints[start:stop]
start += ncol
stop += ncol
if bondlist:
self.bondflag = 1
self.bondlist = bondlist
except:
raise StandardError,"could not read from bond dump file"
# request bonds from data object
elif type(arg) is types.InstanceType and ".data" in str(arg.__class__):
try:
bondlist = []
bondlines = arg.sections["Bonds"]
for line in bondlines:
words = line.split()
bondlist.append([int(words[0]),int(words[1]),
int(words[2]),int(words[3])])
if bondlist:
self.bondflag = 1
self.bondlist = bondlist
except:
raise StandardError,"could not extract bonds from data object"
# request tris/lines from cdata object
elif type(arg) is types.InstanceType and ".cdata" in str(arg.__class__):
try:
tmp,tmp,tmp,tmp,tris,lines = arg.viz(0)
if tris:
self.triflag = 1
self.trilist = tris
if lines:
self.lineflag = 1
self.linelist = lines
except:
raise StandardError,"could not extract tris/lines from cdata object"
# request tris from mdump object
elif type(arg) is types.InstanceType and ".mdump" in str(arg.__class__):
try:
self.triflag = 2
self.triobj = arg
except:
raise StandardError,"could not extract tris from mdump object"
else:
raise StandardError,"unrecognized argument to dump.extra()"
# --------------------------------------------------------------------
def compare_atom(self,a,b):
if a[0] < b[0]:
return -1
elif a[0] > b[0]:
return 1
else:
return 0
# --------------------------------------------------------------------
# one snapshot
class Snap:
pass
# --------------------------------------------------------------------
# time selection class
class tselect:
def __init__(self,data):
self.data = data
# --------------------------------------------------------------------
def all(self):
data = self.data
for snap in data.snaps:
snap.tselect = 1
data.nselect = len(data.snaps)
data.aselect.all()
print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps)
# --------------------------------------------------------------------
def one(self,n):
data = self.data
for snap in data.snaps:
snap.tselect = 0
i = data.findtime(n)
data.snaps[i].tselect = 1
data.nselect = 1
data.aselect.all()
print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps)
# --------------------------------------------------------------------
def none(self):
data = self.data
for snap in data.snaps:
snap.tselect = 0
data.nselect = 0
print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps)
# --------------------------------------------------------------------
def skip(self,n):
data = self.data
count = n-1
for snap in data.snaps:
if not snap.tselect: continue
count += 1
if count == n:
count = 0
continue
snap.tselect = 0
data.nselect -= 1
data.aselect.all()
print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps)
# --------------------------------------------------------------------
def test(self,teststr):
data = self.data
snaps = data.snaps
cmd = "flag = " + teststr.replace("$t","snaps[i].time")
ccmd = compile(cmd,'','single')
for i in xrange(data.nsnaps):
if not snaps[i].tselect: continue
exec ccmd
if not flag:
snaps[i].tselect = 0
data.nselect -= 1
data.aselect.all()
print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps)
# --------------------------------------------------------------------
# atom selection class
class aselect:
def __init__(self,data):
self.data = data
# --------------------------------------------------------------------
def all(self,*args):
data = self.data
if len(args) == 0: # all selected timesteps
for snap in data.snaps:
if not snap.tselect: continue
for i in xrange(snap.natoms): snap.aselect[i] = 1
snap.nselect = snap.natoms
else: # one timestep
n = data.findtime(args[0])
snap = data.snaps[n]
for i in xrange(snap.natoms): snap.aselect[i] = 1
snap.nselect = snap.natoms
# --------------------------------------------------------------------
def test(self,teststr,*args):
data = self.data
# replace all $var with snap.atoms references and compile test string
pattern = "\$\w*"
list = re.findall(pattern,teststr)
for item in list:
name = item[1:]
column = data.names[name]
insert = "snap.atoms[i][%d]" % column
teststr = teststr.replace(item,insert)
cmd = "flag = " + teststr
ccmd = compile(cmd,'','single')
if len(args) == 0: # all selected timesteps
for snap in data.snaps:
if not snap.tselect: continue
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
exec ccmd
if not flag:
snap.aselect[i] = 0
snap.nselect -= 1
for i in xrange(data.nsnaps):
if data.snaps[i].tselect:
print "%d atoms of %d selected in first step %d" % \
(data.snaps[i].nselect,data.snaps[i].natoms,data.snaps[i].time)
break
for i in xrange(data.nsnaps-1,-1,-1):
if data.snaps[i].tselect:
print "%d atoms of %d selected in last step %d" % \
(data.snaps[i].nselect,data.snaps[i].natoms,data.snaps[i].time)
break
else: # one timestep
n = data.findtime(args[0])
snap = data.snaps[n]
for i in xrange(snap.natoms):
if not snap.aselect[i]: continue
exec ccmd
if not flag:
snap.aselect[i] = 0
snap.nselect -= 1
diff --git a/tools/python/pizza/log.py b/tools/python/pizza/log.py
index 6ff694240..7020c218a 100644
--- a/tools/python/pizza/log.py
+++ b/tools/python/pizza/log.py
@@ -1,334 +1,334 @@
# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
#
# Copyright (2005) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
# log tool
oneline = "Read LAMMPS log files and extract thermodynamic data"
docstr = """
l = log("file1") read in one or more log files
l = log("log1 log2.gz") can be gzipped
l = log("file*") wildcard expands to multiple files
l = log("log.lammps",0) two args = store filename, but don't read
incomplete and duplicate thermo entries are deleted
time = l.next() read new thermo info from file
used with 2-argument constructor to allow reading thermo incrementally
return time stamp of last thermo read
return -1 if no new thermo since last read
nvec = l.nvec # of vectors of thermo info
nlen = l.nlen length of each vectors
names = l.names list of vector names
t,pe,... = l.get("Time","KE",...) return one or more vectors of values
l.write("file.txt") write all vectors to a file
l.write("file.txt","Time","PE",...) write listed vectors to a file
get and write allow abbreviated (uniquely) vector names
"""
# History
# 8/05, Steve Plimpton (SNL): original version
# ToDo list
# Variables
# nvec = # of vectors
# nlen = length of each vector
# names = list of vector names
# ptr = dictionary, key = name, value = index into data for which column
# data[i][j] = 2d array of floats, i = 0 to # of entries, j = 0 to nvecs-1
# style = style of LAMMPS log file, 1 = multi, 2 = one, 3 = gran
# firststr = string that begins a thermo section in log file
# increment = 1 if log file being read incrementally
# eof = ptr into incremental file for where to start next read
# Imports and external programs
import sys, re, glob
from os import popen
try: tmp = PIZZA_GUNZIP
except: PIZZA_GUNZIP = "gunzip"
# Class definition
class log:
# --------------------------------------------------------------------
def __init__(self,*list):
self.nvec = 0
self.names = []
self.ptr = {}
self.data = []
# flist = list of all log file names
words = list[0].split()
self.flist = []
for word in words: self.flist += glob.glob(word)
if len(self.flist) == 0 and len(list) == 1:
raise StandardError,"no log file specified"
if len(list) == 1:
self.increment = 0
self.read_all()
else:
if len(self.flist) > 1:
raise StandardError,"can only incrementally read one log file"
self.increment = 1
self.eof = 0
# --------------------------------------------------------------------
# read all thermo from all files
def read_all(self):
self.read_header(self.flist[0])
if self.nvec == 0: raise StandardError,"log file has no values"
# read all files
for file in self.flist: self.read_one(file)
print
# sort entries by timestep, cull duplicates
self.data.sort(self.compare)
self.cull()
self.nlen = len(self.data)
print "read %d log entries" % self.nlen
# --------------------------------------------------------------------
def next(self):
if not self.increment: raise StandardError,"cannot read incrementally"
if self.nvec == 0:
try: open(self.flist[0],'r')
except: return -1
self.read_header(self.flist[0])
if self.nvec == 0: return -1
self.eof = self.read_one(self.flist[0],self.eof)
return int(self.data[-1][0])
# --------------------------------------------------------------------
def get(self,*keys):
if len(keys) == 0:
raise StandardError, "no log vectors specified"
map = []
for key in keys:
if self.ptr.has_key(key):
map.append(self.ptr[key])
else:
count = 0
for i in range(self.nvec):
if self.names[i].find(key) == 0:
count += 1
index = i
if count == 1:
map.append(index)
else:
raise StandardError, "unique log vector %s not found" % key
vecs = []
for i in range(len(keys)):
vecs.append(self.nlen * [0])
for j in xrange(self.nlen):
vecs[i][j] = self.data[j][map[i]]
if len(keys) == 1: return vecs[0]
else: return vecs
# --------------------------------------------------------------------
def write(self,filename,*keys):
if len(keys):
map = []
for key in keys:
if self.ptr.has_key(key):
map.append(self.ptr[key])
else:
count = 0
for i in range(self.nvec):
if self.names[i].find(key) == 0:
count += 1
index = i
if count == 1:
map.append(index)
else:
raise StandardError, "unique log vector %s not found" % key
else:
map = range(self.nvec)
f = open(filename,"w")
for i in xrange(self.nlen):
for j in xrange(len(map)):
print >>f,self.data[i][map[j]],
print >>f
f.close()
# --------------------------------------------------------------------
def compare(self,a,b):
if a[0] < b[0]:
return -1
elif a[0] > b[0]:
return 1
else:
return 0
# --------------------------------------------------------------------
def cull(self):
i = 1
while i < len(self.data):
if self.data[i][0] == self.data[i-1][0]: del self.data[i]
else: i += 1
# --------------------------------------------------------------------
def read_header(self,file):
str_multi = "----- Step"
str_one = "Step "
if file[-3:] == ".gz":
txt = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r').read()
else:
txt = open(file).read()
if txt.find(str_multi) >= 0:
self.firststr = str_multi
self.style = 1
elif txt.find(str_one) >= 0:
self.firststr = str_one
self.style = 2
else:
return
if self.style == 1:
s1 = txt.find(self.firststr)
s2 = txt.find("\n--",s1)
pattern = "\s(\S*)\s*="
keywords = re.findall(pattern,txt[s1:s2])
keywords.insert(0,"Step")
i = 0
for keyword in keywords:
self.names.append(keyword)
self.ptr[keyword] = i
i += 1
else:
s1 = txt.find(self.firststr)
s2 = txt.find("\n",s1)
line = txt[s1:s2]
words = line.split()
for i in range(len(words)):
self.names.append(words[i])
self.ptr[words[i]] = i
self.nvec = len(self.names)
# --------------------------------------------------------------------
def read_one(self,*list):
# if 2nd arg exists set file ptr to that value
# read entire (rest of) file into txt
file = list[0]
if file[-3:] == ".gz":
- f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r')
+ f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'rb')
else:
- f = open(file)
+ f = open(file,'rb')
if len(list) == 2: f.seek(list[1])
txt = f.read()
if file[-3:] == ".gz": eof = 0
else: eof = f.tell()
f.close()
start = last = 0
while not last:
# chunk = contiguous set of thermo entries (line or multi-line)
# s1 = 1st char on 1st line of chunk
# s2 = 1st char on line after chunk
# set last = 1 if this is last chunk in file, leave 0 otherwise
# set start = position in file to start looking for next chunk
# rewind eof if final entry is incomplete
s1 = txt.find(self.firststr,start)
s2 = txt.find("Loop time of",start+1)
if s1 >= 0 and s2 >= 0 and s1 < s2: # found s1,s2 with s1 before s2
if self.style == 2:
s1 = txt.find("\n",s1) + 1
elif s1 >= 0 and s2 >= 0 and s2 < s1: # found s1,s2 with s2 before s1
s1 = 0
elif s1 == -1 and s2 >= 0: # found s2, but no s1
last = 1
s1 = 0
elif s1 >= 0 and s2 == -1: # found s1, but no s2
last = 1
if self.style == 1:
s2 = txt.rfind("\n--",s1) + 1
else:
s1 = txt.find("\n",s1) + 1
s2 = txt.rfind("\n",s1) + 1
eof -= len(txt) - s2
elif s1 == -1 and s2 == -1: # found neither
# could be end-of-file section
# or entire read was one chunk
if txt.find("Loop time of",start) == start: # end of file, so exit
eof -= len(txt) - start # reset eof to "Loop"
break
last = 1 # entire read is a chunk
s1 = 0
if self.style == 1:
s2 = txt.rfind("\n--",s1) + 1
else:
s2 = txt.rfind("\n",s1) + 1
eof -= len(txt) - s2
if s1 == s2: break
chunk = txt[s1:s2-1]
start = s2
# split chunk into entries
# parse each entry for numeric fields, append to data
if self.style == 1:
sections = chunk.split("\n--")
pat1 = re.compile("Step\s*(\S*)\s")
pat2 = re.compile("=\s*(\S*)")
for section in sections:
word1 = [re.search(pat1,section).group(1)]
word2 = re.findall(pat2,section)
words = word1 + word2
self.data.append(map(float,words))
else:
lines = chunk.split("\n")
for line in lines:
words = line.split()
self.data.append(map(float,words))
# print last timestep of chunk
print int(self.data[len(self.data)-1][0]),
sys.stdout.flush()
return eof