From 249a2a6783f286716160f7fd3e8bf633e6b0596c Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Wed, 2 Jun 2021 13:09:52 -0400 Subject: [PATCH] Sync copies of pizza --- python/examples/pizza/dump.py | 76 ++-- python/examples/pizza/gnu.py | 48 +-- python/examples/pizza/pdbfile.py | 28 +- tools/python/pizza/dump.py | 684 +++++++++++++------------------ tools/python/pizza/gnu.py | 47 ++- tools/python/pizza/pdbfile.py | 46 ++- 6 files changed, 424 insertions(+), 505 deletions(-) diff --git a/python/examples/pizza/dump.py b/python/examples/pizza/dump.py index 443ced1949..5c7fab33ae 100644 --- a/python/examples/pizza/dump.py +++ b/python/examples/pizza/dump.py @@ -3,7 +3,7 @@ # # 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 +# certain rights in this software. This software is distributed under # the GNU General Public License. # for python3 compatibility @@ -35,7 +35,7 @@ time = d.next() read next snapshot from dump files 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 @@ -227,7 +227,7 @@ class dump: for word in words: self.flist += glob.glob(word) if len(self.flist) == 0 and len(list) == 1: raise Exception("no dump file specified") - + if len(list) == 1: self.increment = 0 self.read_all() @@ -270,7 +270,7 @@ class dump: 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): @@ -341,7 +341,7 @@ class dump: # 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() @@ -414,7 +414,7 @@ class dump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise Exception("dump map() requires pairs of mappings") @@ -492,7 +492,7 @@ class dump: 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 @@ -505,7 +505,7 @@ class dump: 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 @@ -527,7 +527,7 @@ class dump: 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 @@ -542,7 +542,7 @@ class dump: def owrap(self,other): print("Wrapping to other ...") - + id = self.names["id"] x = self.names["x"] y = self.names["y"] @@ -551,7 +551,7 @@ class dump: 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 @@ -568,7 +568,7 @@ class dump: # -------------------------------------------------------------------- # convert column names assignment to a string, in column order - + def names2str(self): ncol = len(self.snaps[0].atoms[0]) pairs = self.names.items() @@ -631,7 +631,7 @@ class dump: print(snap.ylo,snap.yhi,file=f) print(snap.zlo,snap.zhi,file=f) print("ITEM: ATOMS",namestr,file=f) - + atoms = snap.atoms nvalues = len(atoms[0]) for i in range(snap.natoms): @@ -655,7 +655,7 @@ class dump: if not snap.tselect: continue print(snap.time,end='') sys.stdout.flush() - + file = root + "." + str(snap.time) f = open(file,"w") print("ITEM: TIMESTEP",file=f) @@ -667,7 +667,7 @@ class dump: print(snap.ylo,snap.yhi,file=f) print(snap.zlo,snap.zhi,file=f) print("ITEM: ATOMS",namestr,file=f) - + atoms = snap.atoms nvalues = len(atoms[0]) for i in range(snap.natoms): @@ -709,7 +709,7 @@ class dump: 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] @@ -721,7 +721,7 @@ class dump: if not snap.tselect: continue for i in range(snap.natoms): if snap.aselect[i]: exec(ceq) - + # -------------------------------------------------------------------- # set a column value via an input vec for all selected snapshots/atoms @@ -741,7 +741,7 @@ class dump: if snap.aselect[i]: atoms[i][icol] = vec[m] m += 1 - + # -------------------------------------------------------------------- # clone value in col across selected timesteps for atoms with same ID @@ -807,7 +807,7 @@ class dump: columns.append(self.names[name]) values.append(self.nselect * [0]) ncol = len(columns) - + id = self.names["id"] m = 0 for snap in self.snaps: @@ -823,13 +823,13 @@ class dump: 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 Exception("no columns specified") columns = [] @@ -884,7 +884,7 @@ class dump: del self.snaps[i] else: i += 1 - + # -------------------------------------------------------------------- # iterate over selected snapshots @@ -896,11 +896,11 @@ class dump: 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] @@ -914,7 +914,7 @@ class dump: # create atom list needed by viz from id,type,x,y,z # need Numeric/Numpy mode here - + atoms = [] for i in range(snap.natoms): if not snap.aselect[i]: continue @@ -948,12 +948,12 @@ class dump: 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): @@ -997,7 +997,7 @@ class dump: def extra(self,arg): # read bonds from bond dump file - + if type(arg) is types.StringType: try: f = open(arg,'r') @@ -1017,7 +1017,7 @@ class dump: f.close() # convert values to int and absolute value since can be negative types - + if oldnumeric: bondlist = np.zeros((nbonds,4),np.Int) else: bondlist = np.zeros((nbonds,4),np.int) ints = [abs(int(value)) for value in words] @@ -1032,9 +1032,9 @@ class dump: self.bondlist = bondlist except: raise Exception("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 = [] @@ -1050,7 +1050,7 @@ class dump: raise Exception("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) @@ -1064,7 +1064,7 @@ class dump: raise Exception("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 @@ -1074,7 +1074,7 @@ class dump: else: raise Exception("unrecognized argument to dump.extra()") - + # -------------------------------------------------------------------- def compare_atom(self,a,b): @@ -1083,7 +1083,7 @@ class dump: elif a[0] > b[0]: return 1 else: - return 0 + return 0 # -------------------------------------------------------------------- # one snapshot @@ -1098,7 +1098,7 @@ class tselect: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def all(self): @@ -1145,7 +1145,7 @@ class tselect: data.nselect -= 1 data.aselect.all() print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) - + # -------------------------------------------------------------------- def test(self,teststr): @@ -1191,7 +1191,7 @@ class aselect: 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: diff --git a/python/examples/pizza/gnu.py b/python/examples/pizza/gnu.py index 26dc0a5d5a..6e0fc1ee0b 100644 --- a/python/examples/pizza/gnu.py +++ b/python/examples/pizza/gnu.py @@ -3,7 +3,7 @@ # # 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 +# certain rights in this software. This software is distributed under # the GNU General Public License. # for python3 compatibility @@ -16,7 +16,7 @@ oneline = "Create plots via GnuPlot plotting program" docstr = """ g = gnu() start up GnuPlot g.stop() shut down GnuPlot process - + g.plot(a) plot vector A against linear index g.plot(a,b) plot B against A g.plot(a,b,c,d,...) plot B against A, D against C, etc @@ -35,14 +35,14 @@ g("plot 'file.dat' using 2:3 with lines") execute string in GnuPlot g.enter() enter GnuPlot shell gnuplot> plot sin(x) with lines type commands directly to GnuPlot gnuplot> exit, quit exit GnuPlot shell - + g.export("data",range(100),a,...) create file with columns of numbers all vectors must be of equal length could plot from file with GnuPlot command: plot 'data' using 1:2 with lines g.select(N) figure N becomes the current plot - + subsequent commands apply to this plot g.hide(N) delete window for figure N @@ -98,7 +98,7 @@ except ImportError: PIZZA_GNUTERM = "x11" # Class definition class gnu: - + # -------------------------------------------------------------------- def __init__(self): @@ -106,7 +106,7 @@ class gnu: self.file = "tmp.gnu" self.figures = [] self.select(1) - + # -------------------------------------------------------------------- def stop(self): @@ -118,7 +118,7 @@ class gnu: def __call__(self,command): self.GNUPLOT.write(command + '\n') self.GNUPLOT.flush() - + # -------------------------------------------------------------------- def enter(self): @@ -159,7 +159,7 @@ class gnu: if i: partial_vecs.append(vec[:i]) else: partial_vecs.append([0]) self.plot(*partial_vecs) - + if n < 10: newfile = file + "000" + str(n) elif n < 100: newfile = file + "00" + str(n) elif n < 1000: newfile = file + "0" + str(n) @@ -167,7 +167,7 @@ class gnu: self.save(newfile) n += 1 - + # -------------------------------------------------------------------- # write list of equal-length vectors to filename @@ -208,7 +208,7 @@ class gnu: # do not continue until plot file is written out # else script could go forward and change data file # use tmp.done as semaphore to indicate plot is finished - + def save(self,file): self.__call__("set terminal postscript enhanced solid lw 2 color portrait") cmd = "set output '%s.eps'" % file @@ -219,7 +219,7 @@ class gnu: while not os.path.exists("tmp.done"): continue self.__call__("set output") self.select(self.current) - + # -------------------------------------------------------------------- # restore default attributes by creating a new fig object @@ -228,7 +228,7 @@ class gnu: fig.ncurves = self.figures[self.current-1].ncurves self.figures[self.current-1] = fig self.draw() - + # -------------------------------------------------------------------- def aspect(self,value): @@ -252,12 +252,12 @@ class gnu: else: self.figures[self.current-1].ylimit = (values[0],values[1]) self.draw() - + # -------------------------------------------------------------------- def label(self,x,y,text): self.figures[self.current-1].labels.append((x,y,text)) - self.figures[self.current-1].nlabels += 1 + self.figures[self.current-1].nlabels += 1 self.draw() # -------------------------------------------------------------------- @@ -266,7 +266,7 @@ class gnu: self.figures[self.current-1].nlabel = 0 self.figures[self.current-1].labels = [] self.draw() - + # -------------------------------------------------------------------- def title(self,*strings): @@ -283,13 +283,13 @@ class gnu: def xtitle(self,label): self.figures[self.current-1].xtitle = label self.draw() - + # -------------------------------------------------------------------- def ytitle(self,label): self.figures[self.current-1].ytitle = label self.draw() - + # -------------------------------------------------------------------- def xlog(self): @@ -298,7 +298,7 @@ class gnu: else: self.figures[self.current-1].xlog = 1 self.draw() - + # -------------------------------------------------------------------- def ylog(self): @@ -307,7 +307,7 @@ class gnu: else: self.figures[self.current-1].ylog = 1 self.draw() - + # -------------------------------------------------------------------- def curve(self,num,color): @@ -323,10 +323,10 @@ class gnu: def draw(self): fig = self.figures[self.current-1] if not fig.ncurves: return - + cmd = 'set size ratio ' + str(1.0/float(fig.aspect)) self.__call__(cmd) - + cmd = 'set title ' + '"' + fig.title + '"' self.__call__(cmd) cmd = 'set xlabel ' + '"' + fig.xtitle + '"' @@ -338,11 +338,11 @@ class gnu: else: self.__call__("unset logscale x") if fig.ylog: self.__call__("set logscale y") else: self.__call__("unset logscale y") - if fig.xlimit: + if fig.xlimit: cmd = 'set xr [' + str(fig.xlimit[0]) + ':' + str(fig.xlimit[1]) + ']' self.__call__(cmd) else: self.__call__("set xr [*:*]") - if fig.ylimit: + if fig.ylimit: cmd = 'set yr [' + str(fig.ylimit[0]) + ':' + str(fig.ylimit[1]) + ']' self.__call__(cmd) else: self.__call__("set yr [*:*]") @@ -372,7 +372,7 @@ class figure: def __init__(self): self.ncurves = 0 - self.colors = [] + self.colors = [] self.title = "" self.xtitle = "" self.ytitle = "" diff --git a/python/examples/pizza/pdbfile.py b/python/examples/pizza/pdbfile.py index 64e77db75f..efdf32fab1 100644 --- a/python/examples/pizza/pdbfile.py +++ b/python/examples/pizza/pdbfile.py @@ -3,7 +3,7 @@ # # 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 +# certain rights in this software. This software is distributed under # the GNU General Public License. # for python3 compatibility @@ -25,7 +25,7 @@ p = pdbfile("3CRO",d) read in single PDB file with snapshot data if only one 4-char file specified and it is not found, it will be downloaded from http://www.rcsb.org as 3CRO.pdb d arg is object with atom coordinates (dump, data) - + p.one() write all output as one big PDB file to tmp.pdb p.one("mine") write to mine.pdb p.many() write one PDB file per snapshot: tmp0000.pdb, ... @@ -39,7 +39,7 @@ p.single(N,"new") write as new.pdb if one file in str arg and d: one new PDB file per snapshot using input PDB file as template multiple input PDB files with a d is not allowed - + index,time,flag = p.iterator(0) index,time,flag = p.iterator(1) @@ -97,7 +97,7 @@ class pdbfile: # flist = full list of all PDB input file names # append .pdb if needed - + if filestr: list = filestr.split() flist = [] @@ -117,7 +117,7 @@ class pdbfile: raise Exception("no input PDB file(s)") # grab PDB file from http://rcsb.org if not a local file - + if len(self.files) == 1 and len(self.files[0]) == 8: try: open(self.files[0],'r').close() @@ -127,7 +127,7 @@ class pdbfile: urllib.urlretrieve(fetchstr,self.files[0]) if self.data and len(self.files): self.read_template(self.files[0]) - + # -------------------------------------------------------------------- # write a single large PDB file for concatenating all input data or files # if data exists: @@ -145,7 +145,7 @@ class pdbfile: f = open(file,'w') # use template PDB file with each snapshot - + if self.data: n = flag = 0 while 1: @@ -163,7 +163,7 @@ class pdbfile: print("END",file=f) print(file,end='') sys.stdout.flush() - + f.close() print("\nwrote %d datasets to %s in PDB format" % (n,file)) @@ -199,7 +199,7 @@ class pdbfile: f = open(file,'w') self.convert(f,which) f.close() - + print(time,end='') sys.stdout.flush() n += 1 @@ -216,13 +216,13 @@ class pdbfile: else: file = root + str(n) file += ".pdb" - + f = open(file,'w') f.write(open(infile,'r').read()) f.close() print(file,end='') sys.stdout.flush() - + n += 1 print("\nwrote %d datasets to %s*.pdb in PDB format" % (n,root)) @@ -249,7 +249,7 @@ class pdbfile: self.convert(f,which) else: f.write(open(self.files[time],'r').read()) - + f.close() # -------------------------------------------------------------------- @@ -268,8 +268,8 @@ class pdbfile: # -------------------------------------------------------------------- # read a PDB file and store ATOM lines - - def read_template(self,file): + + def read_template(self,file): lines = open(file,'r').readlines() self.atomlines = {} for line in lines: diff --git a/tools/python/pizza/dump.py b/tools/python/pizza/dump.py index 1c6eb5edfd..5c7fab33ae 100644 --- a/tools/python/pizza/dump.py +++ b/tools/python/pizza/dump.py @@ -6,21 +6,25 @@ # certain rights in this software. This software is distributed under # the GNU General Public License. +# for python3 compatibility + +from __future__ import print_function + # 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 +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 - self-describing column names assigned -time = d.next() read next snapshot from dump files +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 @@ -28,24 +32,24 @@ time = d.next() read next snapshot from dump files 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 columns (1-N) +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.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 +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.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 @@ -56,24 +60,24 @@ d.aselect.test("$id > 100 and $type == 2",N) select matching atoms in one step Python comparison syntax: == != < > <= >= and or $name must end with a space -d.write("file") write selected steps/atoms to dump file -d.write("file",head,app) write selected steps/atoms to dump file -d.scatter("tmp") write selected steps/atoms to multiple files +d.write("file") write selected steps/atoms to dump file +d.write("file",head,app) write selected steps/atoms to dump file +d.scatter("tmp") write selected steps/atoms to multiple files write() can be specified with 2 additional flags - head = 0/1 for no/yes snapshot header, app = 0/1 for write vs append + headd = 0/1 for no/yes snapshot header, app = 0/1 for write vs append 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 +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 @@ -85,8 +89,8 @@ d.sort(1000) sort atoms in timestep N 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.setv("type",vector) set a column to a vector of values -d.spread("ke",N,"color") 2nd col = N ints spread over 1st col -d.clone(1000,"color") clone timestep N values to other steps +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 @@ -107,17 +111,18 @@ d.clone(1000,"color") clone timestep N values to other steps 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 +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,lines = d.viz(index) return list of viz objects -d.atype = "color" set column returned as "type" by viz -d.extra(obj) extract bond/tri/line info from obj +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 @@ -125,26 +130,24 @@ d.extra(obj) extract bond/tri/line info from obj time = timestep value flag = -1 when iteration is done, 1 otherwise viz() returns info for selected atoms for specified timestep index - can also call as viz(time,1) and will find index of preceding snapshot time = timestep value - box = \[xlo,ylo,zlo,xhi,yhi,zhi\] + 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 extra() used to define bonds, else NULL + 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() used to define tris, else NULL + 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() used to define lines, else NULL + if extra() was used to define lines, else empty list atype is column name viz() will return as atom type (def = "type") - extra() extracts bonds/tris/lines from obj each time viz() is called - obj can be data object for bonds, cdata object for tris and lines, - bdump object for bonds, tdump object for tris, ldump object for lines. - mdump object for tris + extra() stores list of bonds/tris/lines to return each time viz() is called """ # History # 8/05, Steve Plimpton (SNL): original version # 12/09, David Hart (SNL): allow use of NumPy or Numeric +# 03/17, Richard Berger (Temple U): improve Python 3 compatibility, +# simplify read_snapshot by using reshape # ToDo list # try to optimize this line in read_snap: words += f.readline().split() @@ -156,7 +159,6 @@ d.extra(obj) extract bond/tri/line info from obj # 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() -# scale_original = 0/1/-1 if coords were read in as unscaled/scaled/unknown # nsnaps = # of snapshots # nselect = # of selected snapshots # snaps = list of snapshots @@ -165,36 +167,34 @@ d.extra(obj) extract bond/tri/line info from obj # 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, 2 if dynamic -# bondlist = static list of bonds to return w/ viz() for all snapshots +# 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 w/ viz() for all snapshots -# lineflag = 0 if no lines, 1 if they are defined statically, 2 if dynamic -# linelist = static list of lines to return w/ viz() for all snapshots -# objextra = object to get bonds,tris,lines from dynamically +# 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 -# boxstr = format string after BOX BOUNDS, if it exists -# triclinic = 0/1 for orthogonal/triclinic based on BOX BOUNDS fields # nselect = # of selected atoms in this snapshot # aselect[i] = 0/1 for each atom -# xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz = box bounds (float) +# 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 +import sys, re, glob, types from os import popen from math import * # any function could be used by set() try: - import numpy as np - oldnumeric = False + import numpy as np + oldnumeric = False except: - import Numeric as np - oldnumeric = True + import Numeric as np + oldnumeric = True try: from DEFAULTS import PIZZA_GUNZIP except: PIZZA_GUNZIP = "gunzip" @@ -216,9 +216,9 @@ class dump: self.bondlist = [] self.triflag = 0 self.trilist = [] + self.triobj = 0 self.lineflag = 0 self.linelist = [] - self.objextra = None # flist = list of all dump file names @@ -226,7 +226,7 @@ class dump: 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" + raise Exception("no dump file specified") if len(list) == 1: self.increment = 0 @@ -251,48 +251,57 @@ class dump: snap = self.read_snapshot(f) while snap: self.snaps.append(snap) - print snap.time, + print(snap.time,end='') sys.stdout.flush() snap = self.read_snapshot(f) f.close() - print + 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 + print("read %d snapshots" % self.nsnaps) # select all timesteps and atoms self.tselect.all() - # print column assignments + # set default names for atom columns if file wasn't self-describing - if len(self.names): - print "assigned columns:",self.names2str() + if len(self.snaps) == 0: + print("no column assignments made") + elif len(self.names): + print("assigned columns:",self.names2str()) + elif self.snaps[0].atoms is 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:",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:",self.names2str()) else: - print "no column assignments made" + 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 "dump scaling status is unknown" + print("no unscaling could be performed") elif self.nsnaps > 0: - if self.scale_original == 1: self.unscale() - elif self.scale_original == 0: print "dump is already unscaled" - else: print "dump scaling status is unknown" + 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" + if not self.increment: raise Exception("cannot read incrementally") # read next snapshot in current file using eof as pointer # if fail, try next file @@ -321,7 +330,7 @@ class dump: snap = self.snaps[self.nsnaps] snap.tselect = 1 snap.nselect = snap.natoms - for i in xrange(snap.natoms): snap.aselect[i] = 1 + for i in range(snap.natoms): snap.aselect[i] = 1 self.nsnaps += 1 self.nselect += 1 @@ -330,111 +339,93 @@ class dump: # -------------------------------------------------------------------- # read a single snapshot from file f # return snapshot or 0 if failed - # for first snapshot only: - # assign column names (file must be self-describing) - # set scale_original to 0/1/-1 for unscaled/scaled/unknown - # convert xs,xu to x in names + # 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 + snap.time = int(f.readline().decode().split()[0]) # just grab 1st field item = f.readline() - snap.natoms = int(f.readline()) + snap.natoms = int(f.readline().decode()) snap.aselect = np.zeros(snap.natoms) - item = f.readline() - words = item.split("BOUNDS ") - if len(words) == 1: snap.boxstr = "" - else: snap.boxstr = words[1].strip() - if "xy" in snap.boxstr: snap.triclinic = 1 - else: snap.triclinic = 0 - + item = f.readline().decode() words = f.readline().split() - if len(words) == 2: - snap.xlo,snap.xhi,snap.xy = float(words[0]),float(words[1]),0.0 - else: - snap.xlo,snap.xhi,snap.xy = \ - float(words[0]),float(words[1]),float(words[2]) - + snap.xlo,snap.xhi = float(words[0]),float(words[1]) words = f.readline().split() - if len(words) == 2: - snap.ylo,snap.yhi,snap.xz = float(words[0]),float(words[1]),0.0 - else: - snap.ylo,snap.yhi,snap.xz = \ - float(words[0]),float(words[1]),float(words[2]) - + snap.ylo,snap.yhi = float(words[0]),float(words[1]) words = f.readline().split() - if len(words) == 2: - snap.zlo,snap.zhi,snap.yz = float(words[0]),float(words[1]),0.0 - else: - snap.zlo,snap.zhi,snap.yz = \ - float(words[0]),float(words[1]),float(words[2]) + snap.zlo,snap.zhi = float(words[0]),float(words[1]) - item = f.readline() + item = f.readline().decode() if len(self.names) == 0: - self.scale_original = -1 - xflag = yflag = zflag = -1 words = item.split()[2:] if len(words): for i in range(len(words)): - if words[i] == "x" or words[i] == "xu": - xflag = 0 + if words[i] == "xs" or words[i] == "xu": self.names["x"] = i - elif words[i] == "xs" or words[i] == "xsu": - xflag = 1 - self.names["x"] = i - elif words[i] == "y" or words[i] == "yu": - yflag = 0 + elif words[i] == "ys" or words[i] == "yu": self.names["y"] = i - elif words[i] == "ys" or words[i] == "ysu": - yflag = 1 - self.names["y"] = i - elif words[i] == "z" or words[i] == "zu": - zflag = 0 - self.names["z"] = i - elif words[i] == "zs" or words[i] == "zsu": - zflag = 1 + elif words[i] == "zs" or words[i] == "zu": self.names["z"] = i else: self.names[words[i]] = i - if xflag == 0 and yflag == 0 and zflag == 0: self.scale_original = 0 - if xflag == 1 and yflag == 1 and zflag == 1: self.scale_original = 1 if snap.natoms: - words = f.readline().split() + words = f.readline().decode().split() ncol = len(words) - for i in xrange(1,snap.natoms): - words += f.readline().split() + for i in range(1,snap.natoms): + words += f.readline().decode().split() floats = map(float,words) - if oldnumeric: atoms = np.zeros((snap.natoms,ncol),np.Float) - else: atoms = np.zeros((snap.natoms,ncol),np.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 + if oldnumeric: + atom_data = np.array(list(floats),np.Float) + else: + atom_data = np.array(list(floats),np.float) + + snap.atoms = atom_data.reshape((snap.natoms, ncol)) + else: + snap.atoms = None return snap except: - return 0 + return None + + # -------------------------------------------------------------------- + # 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" + raise Exception("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: @@ -443,16 +434,15 @@ class dump: self.nsnaps -= 1 ndel += 1 else: i += 1 - print "%d snapshots deleted" % ndel - print "%d snapshots remaining" % self.nsnaps + print("%d snapshots deleted" % ndel) + print("%d snapshots remaining" % self.nsnaps) # -------------------------------------------------------------------- # scale coords to 0-1 for all snapshots or just one - # use 6 params as h-matrix to treat orthogonal or triclinic boxes def scale(self,*list): if len(list) == 0: - print "Scaling dump ..." + print("Scaling dump ...") x = self.names["x"] y = self.names["y"] z = self.names["z"] @@ -467,56 +457,20 @@ class dump: # -------------------------------------------------------------------- def scale_one(self,snap,x,y,z): - if snap.xy == 0.0 and snap.xz == 0.0 and snap.yz == 0.0: - xprdinv = 1.0 / (snap.xhi - snap.xlo) - yprdinv = 1.0 / (snap.yhi - snap.ylo) - zprdinv = 1.0 / (snap.zhi - snap.zlo) - atoms = snap.atoms - if type(atoms) != types.NoneType: - atoms[:,x] = (atoms[:,x] - snap.xlo) * xprdinv - atoms[:,y] = (atoms[:,y] - snap.ylo) * yprdinv - atoms[:,z] = (atoms[:,z] - snap.zlo) * zprdinv - else: - xlo_bound = snap.xlo; xhi_bound = snap.xhi - ylo_bound = snap.ylo; yhi_bound = snap.yhi - zlo_bound = snap.zlo; zhi_bound = snap.zhi - xy = snap.xy - xz = snap.xz - yz = snap.yz - xlo = xlo_bound - min((0.0,xy,xz,xy+xz)) - xhi = xhi_bound - max((0.0,xy,xz,xy+xz)) - ylo = ylo_bound - min((0.0,yz)) - yhi = yhi_bound - max((0.0,yz)) - zlo = zlo_bound - zhi = zhi_bound - h0 = xhi - xlo - h1 = yhi - ylo - h2 = zhi - zlo - h3 = yz - h4 = xz - h5 = xy - h0inv = 1.0 / h0 - h1inv = 1.0 / h1 - h2inv = 1.0 / h2 - h3inv = yz / (h1*h2) - h4inv = (h3*h5 - h1*h4) / (h0*h1*h2) - h5inv = xy / (h0*h1) - atoms = snap.atoms - if type(atoms) != types.NoneType: - atoms[:,x] = (atoms[:,x] - snap.xlo)*h0inv + \ - (atoms[:,y] - snap.ylo)*h5inv + \ - (atoms[:,z] - snap.zlo)*h4inv - atoms[:,y] = (atoms[:,y] - snap.ylo)*h1inv + \ - (atoms[:,z] - snap.zlo)*h3inv - atoms[:,z] = (atoms[:,z] - snap.zlo)*h2inv + 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 - # use 6 params as h-matrix to treat orthogonal or triclinic boxes def unscale(self,*list): if len(list) == 0: - print "Unscaling dump ..." + print("Unscaling dump ...") x = self.names["x"] y = self.names["y"] z = self.names["z"] @@ -531,45 +485,19 @@ class dump: # -------------------------------------------------------------------- def unscale_one(self,snap,x,y,z): - if snap.xy == 0.0 and snap.xz == 0.0 and snap.yz == 0.0: - xprd = snap.xhi - snap.xlo - yprd = snap.yhi - snap.ylo - zprd = snap.zhi - snap.zlo - atoms = snap.atoms - if type(atoms) != types.NoneType: - atoms[:,x] = snap.xlo + atoms[:,x]*xprd - atoms[:,y] = snap.ylo + atoms[:,y]*yprd - atoms[:,z] = snap.zlo + atoms[:,z]*zprd - else: - xlo_bound = snap.xlo; xhi_bound = snap.xhi - ylo_bound = snap.ylo; yhi_bound = snap.yhi - zlo_bound = snap.zlo; zhi_bound = snap.zhi - xy = snap.xy - xz = snap.xz - yz = snap.yz - xlo = xlo_bound - min((0.0,xy,xz,xy+xz)) - xhi = xhi_bound - max((0.0,xy,xz,xy+xz)) - ylo = ylo_bound - min((0.0,yz)) - yhi = yhi_bound - max((0.0,yz)) - zlo = zlo_bound - zhi = zhi_bound - h0 = xhi - xlo - h1 = yhi - ylo - h2 = zhi - zlo - h3 = yz - h4 = xz - h5 = xy - atoms = snap.atoms - if type(atoms) != types.NoneType: - atoms[:,x] = snap.xlo + atoms[:,x]*h0 + atoms[:,y]*h5 + atoms[:,z]*h4 - atoms[:,y] = snap.ylo + atoms[:,y]*h1 + atoms[:,z]*h3 - atoms[:,z] = snap.zlo + atoms[:,z]*h2 + 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 ..." + print("Wrapping dump ...") x = self.names["x"] y = self.names["y"] @@ -591,7 +519,7 @@ class dump: # unwrap coords from inside box to outside def unwrap(self): - print "Unwrapping dump ..." + print("Unwrapping dump ...") x = self.names["x"] y = self.names["y"] @@ -611,10 +539,9 @@ class dump: # -------------------------------------------------------------------- # wrap coords to same image as atom ID stored in "other" column - # if dynamic extra lines or triangles defined, owrap them as well def owrap(self,other): - print "Wrapping to other ..." + print("Wrapping to other ...") id = self.names["id"] x = self.names["x"] @@ -631,26 +558,23 @@ class dump: zprd = snap.zhi - snap.zlo atoms = snap.atoms ids = {} - for i in xrange(snap.natoms): + for i in range(snap.natoms): ids[atoms[i][id]] = i - for i in xrange(snap.natoms): + for i in range(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 - # should bonds also be owrapped ? - if self.lineflag == 2 or self.triflag == 2: - self.objextra.owrap(snap.time,xprd,yprd,zprd,ids,atoms,iother,ix,iy,iz) # -------------------------------------------------------------------- # 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() - ncol = len(pairs) str = "" - for i in xrange(ncol): + for i in range(ncol): if i in values: str += pairs[values.index(i)][0] + ' ' return str @@ -661,12 +585,12 @@ class dump: def sort(self,*list): if len(list) == 0: - print "Sorting selected snapshots ..." + 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] + 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) @@ -682,7 +606,7 @@ class dump: atoms = snap.atoms ids = atoms[:,id] ordering = np.argsort(ids) - for i in xrange(len(atoms[0])): + for i in range(len(atoms[0])): atoms[:,i] = np.take(atoms[:,i],ordering) # -------------------------------------------------------------------- @@ -692,47 +616,35 @@ class dump: if len(self.snaps): namestr = self.names2str() if not append: f = open(file,"w") else: f = open(file,"a") - - if "id" in self.names: id = self.names["id"] - else: id = -1 - if "type" in self.names: type = self.names["type"] - else: type = -1 - for snap in self.snaps: if not snap.tselect: continue - print snap.time, + print(snap.time,end='') sys.stdout.flush() if header: - print >>f,"ITEM: TIMESTEP" - print >>f,snap.time - print >>f,"ITEM: NUMBER OF ATOMS" - print >>f,snap.nselect - if snap.boxstr: print >>f,"ITEM: BOX BOUNDS",snap.boxstr - else: print >>f,"ITEM: BOX BOUNDS" - if snap.triclinic: - print >>f,snap.xlo,snap.xhi,snap.xy - print >>f,snap.ylo,snap.yhi,snap.xz - print >>f,snap.zlo,snap.zhi,snap.yz - else: - print >>f,snap.xlo,snap.xhi - print >>f,snap.ylo,snap.yhi - print >>f,snap.zlo,snap.zhi - print >>f,"ITEM: ATOMS",namestr + print("ITEM: TIMESTEP",file=f) + print(snap.time,file=f) + print("ITEM: NUMBER OF ATOMS",file=f) + print(snap.nselect,file=f) + print("ITEM: BOX BOUNDS",file=f) + print(snap.xlo,snap.xhi,file=f) + print(snap.ylo,snap.yhi,file=f) + print(snap.zlo,snap.zhi,file=f) + print("ITEM: ATOMS",namestr,file=f) atoms = snap.atoms nvalues = len(atoms[0]) - for i in xrange(snap.natoms): + for i in range(snap.natoms): if not snap.aselect[i]: continue line = "" - for j in xrange(nvalues): - if j == id or j == type: + for j in range(nvalues): + if (j < 2): line += str(int(atoms[i][j])) + " " else: line += str(atoms[i][j]) + " " - print >>f,line + print(line,file=f) f.close() - print "\n%d snapshots" % self.nselect + print("\n%d snapshots" % self.nselect) # -------------------------------------------------------------------- # write one dump file per snapshot from current selection @@ -741,40 +653,34 @@ class dump: if len(self.snaps): namestr = self.names2str() for snap in self.snaps: if not snap.tselect: continue - print snap.time, + print(snap.time,end='') 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 - if snap.boxstr: print >>f,"ITEM: BOX BOUNDS",snap.boxstr - else: print >>f,"ITEM: BOX BOUNDS" - if snap.triclinic: - print >>f,snap.xlo,snap.xhi,snap.xy - print >>f,snap.ylo,snap.yhi,snap.xz - print >>f,snap.zlo,snap.zhi,snap.yz - else: - print >>f,snap.xlo,snap.xhi - print >>f,snap.ylo,snap.yhi - print >>f,snap.zlo,snap.zhi - print >>f,"ITEM: ATOMS",namestr + print("ITEM: TIMESTEP",file=f) + print(snap.time,file=f) + print("ITEM: NUMBER OF ATOMS",file=f) + print(snap.nselect,file=f) + print("ITEM: BOX BOUNDS",file=f) + print(snap.xlo,snap.xhi,file=f) + print(snap.ylo,snap.yhi,file=f) + print(snap.zlo,snap.zhi,file=f) + print("ITEM: ATOMS",namestr,file=f) atoms = snap.atoms nvalues = len(atoms[0]) - for i in xrange(snap.natoms): + for i in range(snap.natoms): if not snap.aselect[i]: continue line = "" - for j in xrange(nvalues): + for j in range(nvalues): if (j < 2): line += str(int(atoms[i][j])) + " " else: line += str(atoms[i][j]) + " " - print >>f,line + print(line,file=f) f.close() - print "\n%d snapshots" % self.nselect + print("\n%d snapshots" % self.nselect) # -------------------------------------------------------------------- # find min/max across all selected snapshots/atoms for a particular column @@ -786,7 +692,7 @@ class dump: for snap in self.snaps: if not snap.tselect: continue atoms = snap.atoms - for i in xrange(snap.natoms): + for i in range(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] @@ -796,7 +702,7 @@ class dump: # set a column value via an equation for all selected snapshots def set(self,eq): - print "Setting ..." + print("Setting ...") pattern = "\$\w*" list = re.findall(pattern,eq) @@ -813,14 +719,14 @@ class dump: for snap in self.snaps: if not snap.tselect: continue - for i in xrange(snap.natoms): - if snap.aselect[i]: exec ceq + for i in range(snap.natoms): + if snap.aselect[i]: exec(ceq) # -------------------------------------------------------------------- # set a column value via an input vec for all selected snapshots/atoms def setv(self,colname,vec): - print "Setting ..." + print("Setting ...") if not self.names.has_key(colname): self.newcolumn(colname) icol = self.names[colname] @@ -828,10 +734,10 @@ class dump: for snap in self.snaps: if not snap.tselect: continue if snap.nselect != len(vec): - raise StandardError,"vec length does not match # of selected atoms" + raise Exception("vec length does not match # of selected atoms") atoms = snap.atoms m = 0 - for i in xrange(snap.natoms): + for i in range(snap.natoms): if snap.aselect[i]: atoms[i][icol] = vec[m] m += 1 @@ -844,12 +750,12 @@ class dump: icol = self.names[col] id = self.names["id"] ids = {} - for i in xrange(self.snaps[istep].natoms): + for i in range(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): + for i in range(snap.natoms): if not snap.aselect[i]: continue j = ids[atoms[i][id]] atoms[i][icol] = self.snaps[istep].atoms[j][icol] @@ -863,14 +769,14 @@ class dump: inew = self.names[new] min,max = self.minmax(old) - print "min/max = ",min,max + 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): + for i in range(snap.natoms): if not snap.aselect[i]: continue ivalue = int((atoms[i][iold] - min) * invdelta) + 1 if ivalue > n: ivalue = n @@ -894,7 +800,7 @@ class dump: def atom(self,n,*list): if len(list) == 0: - raise StandardError, "no columns specified" + raise Exception("no columns specified") columns = [] values = [] for name in list: @@ -907,11 +813,11 @@ class dump: for snap in self.snaps: if not snap.tselect: continue atoms = snap.atoms - for i in xrange(snap.natoms): + for i in range(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): + raise Exception("could not find atom ID in snapshot") + for j in range(ncol): values[j][m] = atoms[i][columns[j]] m += 1 @@ -925,7 +831,7 @@ class dump: snap = self.snaps[self.findtime(n)] if len(list) == 0: - raise StandardError, "no columns specified" + raise Exception("no columns specified") columns = [] values = [] for name in list: @@ -934,9 +840,9 @@ class dump: ncol = len(columns) m = 0 - for i in xrange(snap.natoms): + for i in range(snap.natoms): if not snap.aselect[i]: continue - for j in xrange(ncol): + for j in range(ncol): values[j][m] = snap.atoms[i][columns[j]] m += 1 @@ -985,7 +891,7 @@ class dump: def iterator(self,flag): start = 0 if flag: start = self.iterate + 1 - for i in xrange(start,self.nsnaps): + for i in range(start,self.nsnaps): if self.snaps[i].tselect: self.iterate = i return i,self.snaps[i].time,1 @@ -993,20 +899,9 @@ class dump: # -------------------------------------------------------------------- # return list of atoms to viz for snapshot isnap - # if called with flag, then index is timestep, so convert to snapshot index # augment with bonds, tris, lines if extra() was invoked - def viz(self,index,flag=0): - if not flag: isnap = index - else: - times = self.time() - n = len(times) - i = 0 - while i < n: - if times[i] > index: break - i += 1 - isnap = i - 1 - + def viz(self,isnap): snap = self.snaps[isnap] time = snap.time @@ -1021,27 +916,23 @@ class dump: # need Numeric/Numpy mode here atoms = [] - for i in xrange(snap.natoms): + for i in range(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 bonds from static or dynamic bond list - # then generate bond coords from bondlist + # 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 added to bonds + # any bond with unselected atom is not returned to viz caller # need Numeric/Numpy mode here bonds = [] if self.bondflag: - if self.bondflag == 1: bondlist = self.bondlist - elif self.bondflag == 2: - tmp1,tmp2,tmp3,bondlist,tmp4,tmp5 = self.objextra.viz(time,1) alist = {} - for i in xrange(len(atoms)): alist[int(atoms[i][0])] = i - for bond in bondlist: + for i in range(len(atoms)): alist[int(atoms[i][0])] = i + for bond in self.bondlist: try: i = alist[bond[2]] j = alist[bond[3]] @@ -1051,32 +942,24 @@ class dump: atom2[2],atom2[3],atom2[4],atom1[1],atom2[1]]) except: continue - # create list of tris from static or dynamic tri list - # if dynamic, could eliminate tris for unselected atoms - tris = [] if self.triflag: if self.triflag == 1: tris = self.trilist elif self.triflag == 2: - tmp1,tmp2,tmp3,tmp4,tris,tmp5 = self.objextra.viz(time,1) - - # create list of lines from static or dynamic tri list - # if dynamic, could eliminate lines for unselected atoms + timetmp,boxtmp,atomstmp,bondstmp, \ + tris,linestmp = self.triobj.viz(time,1) lines = [] - if self.lineflag: - if self.lineflag == 1: lines = self.linelist - elif self.lineflag == 2: - tmp1,tmp2,tmp3,tmp4,tmp5,lines = self.objextra.viz(time,1) + 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 + for i, snap in enumerate(self.snaps): + if snap.time == n: return i + raise Exception("no step %d exists" % n) # -------------------------------------------------------------------- # return maximum box size across all selected snapshots @@ -1086,12 +969,12 @@ class dump: 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 + if xlo is None or snap.xlo < xlo: xlo = snap.xlo + if xhi is None or snap.xhi > xhi: xhi = snap.xhi + if ylo is None or snap.ylo < ylo: ylo = snap.ylo + if yhi is None or snap.yhi > yhi: yhi = snap.yhi + if zlo is None or snap.zlo < zlo: zlo = snap.zlo + if zhi is None or snap.zhi > zhi: zhi = snap.zhi return [xlo,ylo,zlo,xhi,yhi,zhi] # -------------------------------------------------------------------- @@ -1103,21 +986,56 @@ class dump: for snap in self.snaps: if not snap.tselect: continue atoms = snap.atoms - for i in xrange(snap.natoms): + for i in range(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 - # if static, grab once, else store obj to grab dynamically def extra(self,arg): - # data object, grab bonds statically + # read bonds from bond dump file - if type(arg) is types.InstanceType and ".data" in str(arg.__class__): - self.bondflag = 0 + 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 Exception("could not read bonds from dump file") + + words = f.readline().split() + ncol = len(words) + for i in range(1,nbonds): + words += f.readline().split() + f.close() + + # convert values to int and absolute value since can be negative types + + if oldnumeric: bondlist = np.zeros((nbonds,4),np.Int) + else: bondlist = np.zeros((nbonds,4),np.int) + ints = [abs(int(value)) for value in words] + start = 0 + stop = 4 + for i in range(nbonds): + bondlist[i] = ints[start:stop] + start += ncol + stop += ncol + if bondlist: + self.bondflag = 1 + self.bondlist = bondlist + except: + raise Exception("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"] @@ -1129,12 +1047,11 @@ class dump: self.bondflag = 1 self.bondlist = bondlist except: - raise StandardError,"could not extract bonds from data object" + raise Exception("could not extract bonds from data object") - # cdata object, grab tris and lines statically + # request tris/lines from cdata object elif type(arg) is types.InstanceType and ".cdata" in str(arg.__class__): - self.triflag = self.lineflag = 0 try: tmp,tmp,tmp,tmp,tris,lines = arg.viz(0) if tris: @@ -1144,34 +1061,19 @@ class dump: self.lineflag = 1 self.linelist = lines except: - raise StandardError,"could not extract tris/lines from cdata object" + raise Exception("could not extract tris/lines from cdata object") - # mdump object, grab tris dynamically + # request tris from mdump object elif type(arg) is types.InstanceType and ".mdump" in str(arg.__class__): - self.triflag = 2 - self.objextra = arg - - # bdump object, grab bonds dynamically - - elif type(arg) is types.InstanceType and ".bdump" in str(arg.__class__): - self.bondflag = 2 - self.objextra = arg - - # ldump object, grab lines dynamically - - elif type(arg) is types.InstanceType and ".ldump" in str(arg.__class__): - self.lineflag = 2 - self.objextra = arg - - # tdump object, grab tris dynamically - - elif type(arg) is types.InstanceType and ".tdump" in str(arg.__class__): - self.triflag = 2 - self.objextra = arg + try: + self.triflag = 2 + self.triobj = arg + except: + raise Exception("could not extract tris from mdump object") else: - raise StandardError,"unrecognized argument to dump.extra()" + raise Exception("unrecognized argument to dump.extra()") # -------------------------------------------------------------------- @@ -1205,7 +1107,7 @@ class tselect: snap.tselect = 1 data.nselect = len(data.snaps) data.aselect.all() - print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) # -------------------------------------------------------------------- @@ -1217,7 +1119,7 @@ class tselect: data.snaps[i].tselect = 1 data.nselect = 1 data.aselect.all() - print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) # -------------------------------------------------------------------- @@ -1226,7 +1128,7 @@ class tselect: for snap in data.snaps: snap.tselect = 0 data.nselect = 0 - print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) # -------------------------------------------------------------------- @@ -1242,7 +1144,7 @@ class tselect: snap.tselect = 0 data.nselect -= 1 data.aselect.all() - print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) # -------------------------------------------------------------------- @@ -1251,14 +1153,14 @@ class tselect: snaps = data.snaps cmd = "flag = " + teststr.replace("$t","snaps[i].time") ccmd = compile(cmd,'','single') - for i in xrange(data.nsnaps): + for i in range(data.nsnaps): if not snaps[i].tselect: continue - exec ccmd + 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) + print("%d snapshots selected out of %d" % (data.nselect,data.nsnaps)) # -------------------------------------------------------------------- # atom selection class @@ -1275,12 +1177,12 @@ class aselect: 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 + for i in range(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 + for i in range(snap.natoms): snap.aselect[i] = 1 snap.nselect = snap.natoms # -------------------------------------------------------------------- @@ -1303,29 +1205,29 @@ class aselect: if len(args) == 0: # all selected timesteps for snap in data.snaps: if not snap.tselect: continue - for i in xrange(snap.natoms): + for i in range(snap.natoms): if not snap.aselect[i]: continue - exec ccmd + exec(ccmd) if not flag: snap.aselect[i] = 0 snap.nselect -= 1 - for i in xrange(data.nsnaps): + for i in range(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) + 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): + for i in range(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) + 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): + for i in range(snap.natoms): if not snap.aselect[i]: continue - exec ccmd + exec(ccmd) if not flag: snap.aselect[i] = 0 snap.nselect -= 1 diff --git a/tools/python/pizza/gnu.py b/tools/python/pizza/gnu.py index d99ab3811d..6e0fc1ee0b 100644 --- a/tools/python/pizza/gnu.py +++ b/tools/python/pizza/gnu.py @@ -6,17 +6,20 @@ # certain rights in this software. This software is distributed under # the GNU General Public License. +# for python3 compatibility +from __future__ import print_function + # gnu tool oneline = "Create plots via GnuPlot plotting program" docstr = """ -g = gnu() start up GnuPlot -g.stop() shut down GnuPlot process +g = gnu() start up GnuPlot +g.stop() shut down GnuPlot process g.plot(a) plot vector A against linear index -g.plot(a,b) plot B against A -g.plot(a,b,c,d,...) plot B against A, D against C, etc +g.plot(a,b) plot B against A +g.plot(a,b,c,d,...) plot B against A, D against C, etc g.mplot(M,N,S,"file",a,b,...) multiple plots saved to file0000.eps, etc each plot argument can be a tuple, list, or Numeric/NumPy vector @@ -29,21 +32,21 @@ g.mplot(M,N,S,"file",a,b,...) multiple plots saved to file0000.eps, etc g("plot 'file.dat' using 2:3 with lines") execute string in GnuPlot -g.enter() enter GnuPlot shell +g.enter() enter GnuPlot shell gnuplot> plot sin(x) with lines type commands directly to GnuPlot -gnuplot> exit, quit exit GnuPlot shell +gnuplot> exit, quit exit GnuPlot shell g.export("data",range(100),a,...) create file with columns of numbers all vectors must be of equal length could plot from file with GnuPlot command: plot 'data' using 1:2 with lines -g.select(N) figure N becomes the current plot +g.select(N) figure N becomes the current plot subsequent commands apply to this plot -g.hide(N) delete window for figure N -g.save("file") save current plot as file.eps +g.hide(N) delete window for figure N +g.save("file") save current plot as file.eps Set attributes for current plot: @@ -84,12 +87,13 @@ g.curve(N,'r') set color of curve N # Imports and external programs -import types, os +import os +import sys try: from DEFAULTS import PIZZA_GNUPLOT -except: PIZZA_GNUPLOT = "gnuplot" +except ImportError: PIZZA_GNUPLOT = "gnuplot -p" try: from DEFAULTS import PIZZA_GNUTERM -except: PIZZA_GNUTERM = "x11" +except ImportError: PIZZA_GNUTERM = "x11" # Class definition @@ -119,7 +123,10 @@ class gnu: def enter(self): while 1: - command = raw_input("gnuplot> ") + if sys.version_info[0] == 3: + command = input("gnuplot> ") + else: + command = raw_input("gnuplot> ") if command == "quit" or command == "exit": return self.__call__(command) @@ -133,7 +140,7 @@ class gnu: self.export(file,linear,vectors[0]) self.figures[self.current-1].ncurves = 1 else: - if len(vectors) % 2: raise StandardError,"vectors must come in pairs" + if len(vectors) % 2: raise Exception("vectors must come in pairs") for i in range(0,len(vectors),2): file = self.file + ".%d.%d" % (self.current,i/2+1) self.export(file,vectors[i],vectors[i+1]) @@ -167,13 +174,13 @@ class gnu: def export(self,filename,*vectors): n = len(vectors[0]) for vector in vectors: - if len(vector) != n: raise StandardError,"vectors must be same length" + if len(vector) != n: raise Exception("vectors must be same length") f = open(filename,'w') nvec = len(vectors) - for i in xrange(n): - for j in xrange(nvec): - print >>f,vectors[j][i], - print >>f + for i in range(n): + for j in range(nvec): + print(str(vectors[j][i])+" ",file=f,end='') + print ("",file=f) f.close() # -------------------------------------------------------------------- @@ -350,7 +357,7 @@ class gnu: self.__call__("set key off") cmd = 'plot ' - for i in range(fig.ncurves): + for i in range(int(fig.ncurves)): file = self.file + ".%d.%d" % (self.current,i+1) if len(fig.colors) > i and fig.colors[i]: cmd += "'" + file + "' using 1:2 with line %d, " % fig.colors[i] diff --git a/tools/python/pizza/pdbfile.py b/tools/python/pizza/pdbfile.py index 9b2238cbd6..efdf32fab1 100644 --- a/tools/python/pizza/pdbfile.py +++ b/tools/python/pizza/pdbfile.py @@ -6,6 +6,9 @@ # certain rights in this software. This software is distributed under # the GNU General Public License. +# for python3 compatibility +from __future__ import print_function + # pdb tool oneline = "Read, write PDB files in combo with LAMMPS snapshots" @@ -51,6 +54,7 @@ index,time,flag = p.iterator(1) # History # 8/05, Steve Plimpton (SNL): original version +# 3/17, Richard Berger (Temple U): improve Python 3 compatibility # ToDo list # for generic PDB file (no template) from a LJ unit system, @@ -64,7 +68,13 @@ index,time,flag = p.iterator(1) # Imports and external programs -import sys, types, glob, urllib +import sys, glob, urllib +PY3 = sys.version_info[0] == 3 + +if PY3: + string_types = str, +else: + string_types = basestring # Class definition @@ -74,7 +84,7 @@ class pdbfile: def __init__(self,*args): if len(args) == 1: - if type(args[0]) is types.StringType: + if type(args[0]) is string_types: filestr = args[0] self.data = None else: @@ -83,7 +93,7 @@ class pdbfile: elif len(args) == 2: filestr = args[0] self.data = args[1] - else: raise StandardError, "invalid args for pdb()" + else: raise Exception("invalid args for pdb()") # flist = full list of all PDB input file names # append .pdb if needed @@ -94,17 +104,17 @@ class pdbfile: for file in list: if '*' in file: flist += glob.glob(file) else: flist.append(file) - for i in xrange(len(flist)): + for i in range(len(flist)): if flist[i][-4:] != ".pdb": flist[i] += ".pdb" if len(flist) == 0: - raise StandardError,"no PDB file specified" + raise Exception("no PDB file specified") self.files = flist else: self.files = [] if len(self.files) > 1 and self.data: - raise StandardError, "cannot use multiple PDB files with data object" + raise Exception("cannot use multiple PDB files with data object") if len(self.files) == 0 and not self.data: - raise StandardError, "no input PDB file(s)" + raise Exception("no input PDB file(s)") # grab PDB file from http://rcsb.org if not a local file @@ -112,7 +122,7 @@ class pdbfile: try: open(self.files[0],'r').close() except: - print "downloading %s from http://rcsb.org" % self.files[0] + print("downloading %s from http://rcsb.org" % self.files[0]) fetchstr = "http://www.rcsb.org/pdb/cgi/export.cgi/%s?format=PDB&pdbId=2cpk&compression=None" % self.files[0] urllib.urlretrieve(fetchstr,self.files[0]) @@ -142,20 +152,20 @@ class pdbfile: which,time,flag = self.data.iterator(flag) if flag == -1: break self.convert(f,which) - print >>f,"END" - print time, + print("END",file=f) + print(time,end='') sys.stdout.flush() n += 1 else: for file in self.files: f.write(open(file,'r').read()) - print >>f,"END" - print file, + print("END",file=f) + print(file,end='') sys.stdout.flush() f.close() - print "\nwrote %d datasets to %s in PDB format" % (n,file) + print("\nwrote %d datasets to %s in PDB format" % (n,file)) # -------------------------------------------------------------------- # write series of numbered PDB files @@ -190,7 +200,7 @@ class pdbfile: self.convert(f,which) f.close() - print time, + print(time,end='') sys.stdout.flush() n += 1 @@ -210,12 +220,12 @@ class pdbfile: f = open(file,'w') f.write(open(infile,'r').read()) f.close() - print file, + print(file,end='') sys.stdout.flush() n += 1 - print "\nwrote %d datasets to %s*.pdb in PDB format" % (n,root) + print("\nwrote %d datasets to %s*.pdb in PDB format" % (n,root)) # -------------------------------------------------------------------- # write a single PDB file @@ -280,10 +290,10 @@ class pdbfile: if self.atomlines.has_key(id): (begin,end) = self.atomlines[id] line = "%s%8.3f%8.3f%8.3f%s" % (begin,atom[2],atom[3],atom[4],end) - print >>f,line, + print(line,file=f,end='') else: for atom in atoms: begin = "ATOM %6d %2d R00 1 " % (atom[0],atom[1]) middle = "%8.3f%8.3f%8.3f" % (atom[2],atom[3],atom[4]) end = " 1.00 0.00 NONE" - print >>f,begin+middle+end + print(begin+middle+end,file=f)