diff --git a/examples/group_energy.py b/examples/group_energy.py index f5e478b..193b867 100644 --- a/examples/group_energy.py +++ b/examples/group_energy.py @@ -6,11 +6,11 @@ # return distance sq between 2 atoms with PBC def distance(box,x1,y1,z1,x2,y2,z2): - + delx = x2 - x1 dely = y2 - y1 delz = z2 - z1 - + xprd = box[3] - box[0] yprd = box[4] - box[1] zprd = box[5] - box[2] @@ -30,7 +30,7 @@ def distance(box,x1,y1,z1,x2,y2,z2): delz += zprd else: delz -= zprd - + distsq = delx*delx + dely*dely + delz*delz return distsq @@ -39,10 +39,10 @@ def distance(box,x1,y1,z1,x2,y2,z2): if len(argv) < 3: raise StandardError,"group_energy.py data.file dump.file1 dump.file2 ..." -dt = data(argv[1]) # data file +dt = data(argv[1]) # data file q = dt.get("Atoms",4) -files = ' '.join(argv[2:]) # dump files +files = ' '.join(argv[2:]) # dump files d = dump(files,0) d.map(1,"id",2,"type",3,"x",4,"y",5,"z") @@ -97,8 +97,8 @@ while 1: rsq = distance(box,x1[i],y1[i],z1[i],x2[j],y2[j],z2[j]) if rsq < maxcut_sq: eng_coul,eng_vdwl = p.single(rsq,typei,type2[j],qi,q[id2[j]]) - e_coul_sum += eng_coul - e_vdwl_sum += eng_vdwl + e_coul_sum += eng_coul + e_vdwl_sum += eng_vdwl print "eng_coul = %g at timestep %d" % (e_coul_sum,time) print "eng_vdwl = %g at timestep %d" % (e_vdwl_sum,time) diff --git a/examples/test_cdata.py b/examples/test_cdata.py index 3965f82..234b145 100644 --- a/examples/test_cdata.py +++ b/examples/test_cdata.py @@ -30,5 +30,5 @@ c.select("A","B","C","linebox","organelle","nuchalf") g = gl(c) v = vcr(g) - + print "all done ... type CTRL-D to exit Pizza.py" diff --git a/examples/test_data.py b/examples/test_data.py index b5c4bd8..84040b7 100644 --- a/examples/test_data.py +++ b/examples/test_data.py @@ -1,5 +1,5 @@ # simple test of data tool -# requires files/data.micelle +# requires files/data.micelle and dump.micelle # creates tmp.data d = data("files/data.micelle") @@ -24,7 +24,7 @@ while 1: index,time,flag = d.iterator(flag) if flag == -1: break time,box,atoms,bonds,tris,lines= d.viz(index) - + d.write("tmp.data") print "all done ... type CTRL-D to exit Pizza.py" diff --git a/examples/test_vcr.py b/examples/test_vcr.py index 4ffa053..9daabd2 100644 --- a/examples/test_vcr.py +++ b/examples/test_vcr.py @@ -1,5 +1,4 @@ # simple test of vcr tool -# requires files/bucky*png files d = dump("files/dump.micelle") dt = data("files/data.micelle") diff --git a/scripts/bond_distribute.py b/scripts/bond_distribute.py index de2a600..a059cf4 100644 --- a/scripts/bond_distribute.py +++ b/scripts/bond_distribute.py @@ -1,6 +1,6 @@ -#!/usr/bin/python +#!/usr/bin/python -# Script: bond_distribute.py +# Script: bond_distribute.py # Purpose: binned bond length distributions by bond type # Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ... # datafile = lammps data file @@ -25,7 +25,7 @@ if len(argv) < 7: raise StandardError, \ "Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ..." -dt = data(argv[1]) +dt = data(argv[1]) nbins = int(argv[2]) rmin = float(argv[3]) rmax = float(argv[4]) @@ -49,8 +49,8 @@ for i in xrange(nbonds): ntypes = max(bond[i][1],ntypes) ntypes = int(ntypes) ncount = ntypes * [0] bin = nbins * [0] -for i in xrange(nbins): - bin[i] = ntypes * [0] +for i in xrange(nbins): + bin[i] = ntypes * [0] # read snapshots one-at-a-time @@ -60,24 +60,24 @@ d.map(1,"id",2,"type",3,"x",4,"y",5,"z") while 1: time = d.next() if time == -1: break - + box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - + xprd = box[3] - box[0] - yprd = box[4] - box[1] + yprd = box[4] - box[1] zprd = box[5] - box[2] - - d.unscale() + + d.unscale() d.sort() x,y,z = d.vecs(time,"x","y","z") - + for i in xrange(nbonds): - + delx = x[jatom[i]] - x[iatom[i]] dely = y[jatom[i]] - y[iatom[i]] delz = z[jatom[i]] - z[iatom[i]] - + if abs(delx) > 0.5*xprd: if delx < 0.0: delx += xprd @@ -93,30 +93,30 @@ while 1: delz += zprd else: delz -= zprd - + r = sqrt(delx*delx + dely*dely + delz*delz) - + ibin = int(nbins*(r - rmin)/(rmax - rmin) + 0.5) - if ((ibin >= 0) and (ibin <= nbins-1)): + if ((ibin >= 0) and (ibin <= nbins-1)): bin[ibin][btype[i]] += nbins ncount[btype[i]] += 1 else: print "Warning: bond distance outside specified range" print "Bond type:", btype[i]+1 print "Bond number:", i - print time, - + print time, + print print "Printing bond distance normalized distribution to",outfile - + fp = open(outfile,"w") rrange = rmax - rmin for i in xrange(nbins): - print >>fp, rmin + rrange*float(i)/float(nbins), + print >>fp, rmin + rrange*float(i)/float(nbins), for j in xrange(ntypes): if (ncount[j] > 0): print >>fp, float(bin[i][j])/float(ncount[j])/rrange, else: - print >>fp, 0.0, - print >>fp + print >>fp, 0.0, + print >>fp fp.close() diff --git a/scripts/cluster.py b/scripts/cluster.py index c932f8a..80cf2ee 100644 --- a/scripts/cluster.py +++ b/scripts/cluster.py @@ -31,11 +31,11 @@ def distance(atom1,atom2,box): x2 = atom2[2] y2 = atom2[3] z2 = atom2[4] - + delx = x2 - x1 dely = y2 - y1 delz = z2 - z1 - + xprd = box[3] - box[0] yprd = box[4] - box[1] zprd = box[5] - box[2] @@ -55,7 +55,7 @@ def distance(atom1,atom2,box): delz += zprd else: delz -= zprd - + distsq = delx*delx + dely*dely + delz*delz return distsq @@ -92,7 +92,7 @@ while 1: sys.stdout.flush() # loop over all type1 atoms - + n = len(atoms) for i in xrange(n): itype = atoms[i][1] @@ -101,15 +101,15 @@ while 1: # loop over all type2 atoms # increment cluster count if distance is within cutoff - + for j in xrange(n): jtype = atoms[j][1] - if jtype != type2 or i == j: continue + if jtype != type2 or i == j: continue distsq = distance(atoms[i],atoms[j],box) if distsq < cutsq: ncount += 1 # increment histogram count - + if ncount >= nbin: cluster[nbin-1] += 1 else: cluster[ncount] += 1 diff --git a/scripts/density.py b/scripts/density.py index 3ed6724..a0c0c6e 100644 --- a/scripts/density.py +++ b/scripts/density.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python # Script: density.py # Purpose: binned atom density by atom type @@ -43,19 +43,19 @@ while 1: bin = nbins * [0] for i in xrange(nbins): bin[i] = ntypes * [0] first = 0 - + box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) + vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) if direction == "x": type,x = d.vecs(time,"type","x") elif direction == "y": type,x = d.vecs(time,"type","y") elif direction == "z": type,x = d.vecs(time,"type","z") - + type = map(int,type) natoms = len(type) for i in xrange(natoms): type[i] -= 1 - + for i in xrange(natoms): ibin = int(nbins*x[i] + 0.5) if (ibin < 0): ibin += nbins @@ -63,15 +63,15 @@ while 1: bin[ibin][type[i]] += nbins/vol nsnaps += 1 print time, - + print print "Printing ", direction, "-directional density distribution in mol/L to",outfile conversion = 1660.53873 # convert from atoms/Angs^3 to mol/L - + fp = open(outfile,"w") for i in xrange(nbins): - print >>fp, float(i)/float(nbins), + print >>fp, float(i)/float(nbins), for j in xrange(ntypes): print >>fp, conversion*bin[i][j]/nsnaps, - print >>fp + print >>fp fp.close() diff --git a/scripts/density2d.py b/scripts/density2d.py index 8823e68..e87ad31 100644 --- a/scripts/density2d.py +++ b/scripts/density2d.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python # Script: density2d.py # Purpose: binned atom density by atom type @@ -65,12 +65,12 @@ while 1: ntypes = int(ntypes) bin = np.zeros(shape=(nbins,nbins,ntypes)) first = 0 - + box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) + vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) - if direction == "z": + if direction == "z": type,x,y,z = d.vecs(time,"type","x","y","z") bidirect = 'x/y' dx = box[3] - box[0] @@ -90,7 +90,7 @@ while 1: y0 = box[2] + float(dy)/float(nbins)/2.0 zmax = min(zmax,box[4]) zmin = max(zmin,box[1]) - elif direction == "x": + elif direction == "x": type,x,y,z = d.vecs(time,"type","y","z","x") bidirect = 'y/z' dx = box[4] - box[1] @@ -105,7 +105,7 @@ while 1: type = map(int,type) natoms = len(type) for i in xrange(natoms): type[i] -= 1 - + for i in xrange(natoms): ibin = int(nbins*x[i]) jbin = int(nbins*y[i]) @@ -119,10 +119,10 @@ while 1: nsnaps += 1 print time, -print +print print "Printing %s-mapped density distribution for %s-slice [%.2f,%.2f] in mol/L to %s" %(bidirect, direction, zmin, zmax, outfile) conversion = 1660.53873 # convert from atoms/Angs^3 to mol/L - + fp = open(outfile,"w") # '''Uncomment for column headers. Commented for consistency with density.py''' # print >>fp, " %8s %8s " %('ra', 'rb'), @@ -134,5 +134,5 @@ for i in xrange(nbins): print >>fp, " %8.3f %8.3f " %(float(i)/float(nbins)*float(dx)+float(x0), float(j)/float(nbins)*float(dy)+float(y0)), for k in xrange(ntypes): print >>fp, " %8.3f " % (conversion*bin[j][i][k]/nsnaps), - print >>fp + print >>fp fp.close() diff --git a/scripts/density_area.py b/scripts/density_area.py index fd10d89..e65a7f0 100644 --- a/scripts/density_area.py +++ b/scripts/density_area.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python # Script: density_area.py # Purpose: binned atom density by atom type and running area under the curve @@ -46,19 +46,19 @@ while 1: bin = nbins * [0] for i in xrange(nbins): bin[i] = ntypes * [0] first = 0 - + box = (d.snaps[-1].xlo,d.snaps[-1].ylo,d.snaps[-1].zlo, d.snaps[-1].xhi,d.snaps[-1].yhi,d.snaps[-1].zhi) - vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) + vol = (box[3] - box[0]) * (box[4] - box[1]) * (box[5] - box[2]) if direction == "x": type,x = d.vecs(time,"type","x") elif direction == "y": type,x = d.vecs(time,"type","y") elif direction == "z": type,x = d.vecs(time,"type","z") - + type = map(int,type) natoms = len(type) for i in xrange(natoms): type[i] -= 1 - + for i in xrange(natoms): ibin = int(nbins*x[i] + 0.5) if (ibin < 0): ibin += nbins @@ -66,12 +66,12 @@ while 1: bin[ibin][type[i]] += nbins/vol nsnaps += 1 print time, - + print print "Printing ",direction,"-directional density distribution in mol/L to", \ outfile conversion = 1660.53873 # convert from atoms/Angs^3 to mol/L - + # Output as x, density_1, area_1, ... fp = open(outfile,"w") @@ -93,5 +93,5 @@ for i in xrange(nbins): yden[i][j] = conversion*bin[i][j]/nsnaps sum[j] += 0.5 * (xden[i] - xden[i-1]) * (yden[i][j] + yden[i-1][j]) print >>fp, yden[i][j], sum[j], - print >>fp + print >>fp fp.close() diff --git a/scripts/distance.py b/scripts/distance.py index 255b86e..2738847 100644 --- a/scripts/distance.py +++ b/scripts/distance.py @@ -1,6 +1,6 @@ #!/usr/bin/python -# Script: distance.py +# Script: distance.py # Purpose: check if any atom pairs are closer than specified distance # Syntax: distance.py maxcut dump.file1 dump.file2 ... # maxcut = flag atoms which are less than this distance apart @@ -13,11 +13,11 @@ from math import sqrt if len(argv) < 3: raise StandardError,"distance.py maxcut dump.file1 dump.file2 ..." - + maxcut = float(argv[1]) maxcut_sq = maxcut*maxcut - -files = ' '.join(argv[2:]) # dump files + +files = ' '.join(argv[2:]) # dump files d = dump(files,0) d.map(1,"id",2,"type",3,"x",4,"y",5,"z") @@ -35,10 +35,10 @@ while 1: xprd = box[3] - box[0] yprd = box[4] - box[1] zprd = box[5] - box[2] - + for i in xrange(n): for j in xrange(i+1,n): - + delx = x[j] - x[i] if abs(delx) > 0.5*xprd: if delx < 0.0: @@ -46,7 +46,7 @@ while 1: else: delx -= xprd if (delx*delx < maxcut_sq): - + dely = y[j] - y[i] if abs(dely) > 0.5*yprd: if dely < 0.0: @@ -54,17 +54,17 @@ while 1: else: dely -= yprd if ((dely*dely + delx*delx) < maxcut_sq): - + delz = z[j] - z[i] if abs(delz) > 0.5*zprd: if delz < 0.0: delz += zprd else: delz -= zprd - + rsq = delx*delx + dely*dely + delz*delz - - if rsq < maxcut_sq: + + if rsq < maxcut_sq: print "time = %d, id[i] = %d, id[j] = %d," \ " type[i] = %d, type[j] = %d, distance = %g" % \ (time, id[i], id[j], type[i], type[j], sqrt(rsq)) diff --git a/scripts/flux.py b/scripts/flux.py index dec15a8..1a4ba05 100644 --- a/scripts/flux.py +++ b/scripts/flux.py @@ -1,4 +1,4 @@ -#!/usr/bin/python +#!/usr/bin/python # Script: flux.py # Purpose: flux of atoms through a user-defined plane @@ -41,38 +41,38 @@ flag = 0 while 1: which,time,flag = d.iterator(flag) if flag == -1: break - + if direction == "x": id,type,x = d.vecs(time,"id","type","x") lo = d.snaps[which].xlo hi = d.snaps[which].xhi elif direction == "y": id,type,x = d.vecs(time,"id","type","y") - lo = d.snaps[which].ylo - hi = d.snaps[which].yhi - elif direction == "z": + lo = d.snaps[which].ylo + hi = d.snaps[which].yhi + elif direction == "z": id,type,x = d.vecs(time,"id","type","z") - lo = d.snaps[which].zlo - hi = d.snaps[which].zhi - - prd = hi - lo - plane = lo + scaled_plane*prd - + lo = d.snaps[which].zlo + hi = d.snaps[which].zhi + + prd = hi - lo + plane = lo + scaled_plane*prd + print time, sys.stdout.flush() - + natoms = len(x) if jconfig == 0: x_initial = (natoms+1) * [0] jconfig += 1 typeflux = ntypes * [0] - + for i in xrange(natoms): id[i] = int(id[i]) type[i] = int(type[i]) if jconfig == 1: x_initial[id[i]] = x[i] if x_initial[id[i]] < plane and x[i] > plane : - crossings = int((x[i] - plane)/prd) + 1 + crossings = int((x[i] - plane)/prd) + 1 typeflux[type[i]] += crossings elif x_initial[id[i]] > plane and x[i] < plane : crossings = int((plane - x[i])/prd) + 1 @@ -83,5 +83,5 @@ while 1: print >>f,typeflux[j+1], print >>f print - + f.close() diff --git a/scripts/movie.py b/scripts/movie.py index 611077c..a38c943 100644 --- a/scripts/movie.py +++ b/scripts/movie.py @@ -4,7 +4,7 @@ # Purpose: create images from LAMMPS dump snapshots # Syntax: movie.py raster/svg theta phi dump.1 dump.2 ... # raster/svg = style of image to create -# theta/phi = vertical (z) and azimuthal angle to view from +# theta/phi = vertical (z) and azimuthal angle to view from # files = one or more dump files # Example: movie.py svg 60 130 dump.* # Author: Steve Plimpton (Sandia) diff --git a/src/DEFAULTS.py b/src/DEFAULTS.py index 1423d0b..4250f3c 100644 --- a/src/DEFAULTS.py +++ b/src/DEFAULTS.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. # -------------- diff --git a/src/animate.py b/src/animate.py index 0b696db..cb672fd 100644 --- a/src/animate.py +++ b/src/animate.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. # animate tool @@ -12,19 +12,19 @@ oneline = "Animate a series of image files" docstr = """ a = animate("image*.png") create GUI to animate set of image files - + Actions (same as GUI widgets): -a.first() go to first frame +a.first() go to first frame a.prev() go to previous frame a.back() play backwards from current frame to start -a.stop() stop on current frame +a.stop() stop on current frame a.play() play from current frame to end -a.next() go to next frame -a.last() go to last frame +a.next() go to next frame +a.last() go to last frame -a.frame(31) set frame slider -a.delay(0.4) set delay slider +a.frame(31) set frame slider +a.delay(0.4) set delay slider """ # History @@ -63,7 +63,7 @@ class animate: self.delay_msec = 0 # convert filestr into full list of files - + list = str.split(filestr) self.files = [] for file in list: self.files += glob.glob(file) @@ -71,18 +71,18 @@ class animate: if self.nframes == 0: raise StandardError, "No files to load" # load all images - + self.images = [] for i in xrange(self.nframes): self.images.append(PhotoImage(file=self.files[i])) # grab Tk instance from main - + from __main__ import tkroot self.tkroot = tkroot # GUI control window - + win1 = Toplevel(tkroot) win1.title("Pizza.py animate tool") @@ -95,7 +95,7 @@ class animate: button6 = Button(holder1,text=">",command=self.next).pack(side=LEFT) button7 = Button(holder1,text=">>",command=self.last).pack(side=LEFT) holder1.pack(side=TOP) - + holder2 = Frame(win1) self.slider_frame = Scale(holder2,from_=0,to=self.nframes-1, command=self.frame,orient=HORIZONTAL, @@ -106,12 +106,12 @@ class animate: self.slider_frame.pack(side=LEFT) self.slider_delay.pack(side=LEFT) holder2.pack(side=TOP) - + holder3 = Frame(win1) self.label_frame = Label(holder3) self.label_frame.pack(side=LEFT) holder3.pack(side=TOP) - + # image window win2 = Toplevel(tkroot) @@ -120,7 +120,7 @@ class animate: tkroot.update_idletasks() # force window to appear # display 1st image - + self.index = 0 self.display(self.index) @@ -135,19 +135,19 @@ class animate: def last(self): self.index = self.nframes - 1 self.display(self.index) - + # -------------------------------------------------------------------- def previous(self): if self.index > 0: self.index -= 1 self.display(self.index) - + # -------------------------------------------------------------------- def next(self): if self.index < self.nframes - 1: self.index += 1 self.display(self.index) - + # -------------------------------------------------------------------- def back(self): @@ -157,7 +157,7 @@ class animate: self.index = self.nframes - 1 self.display(self.index) self.loop() - + # -------------------------------------------------------------------- def play(self): @@ -167,15 +167,15 @@ class animate: self.index = 0 self.display(self.index) self.loop() - + # -------------------------------------------------------------------- def stop(self): self.loop_flag = 0 - + # -------------------------------------------------------------------- # loop forward or back until end of animation - + def loop(self): if self.loop_flag == 1 and self.index == self.nframes - 1: self.loop_flag = 0 @@ -190,7 +190,7 @@ class animate: # -------------------------------------------------------------------- # display a frame corresponding to iframe - + def display(self,iframe): self.image_pane.configure(image=self.images[iframe]) self.slider_frame.set(iframe) @@ -209,4 +209,4 @@ class animate: self.delay_value = float(value) self.slider_delay.set(self.delay_value) self.delay_msec = int(1000*self.delay_value) - + diff --git a/src/bdump.py b/src/bdump.py index 9e4f4af..420c907 100644 --- a/src/bdump.py +++ b/src/bdump.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. # bdump tool @@ -12,14 +12,14 @@ oneline = "Read dump files with bond info" docstr = """ b = bdump("dump.one") read in one or more dump files -b = bdump("dump.1 dump.2.gz") can be gzipped -b = bdump("dump.*") wildcard expands to multiple files -b = bdump("dump.*",0) two args = store filenames, but don't read +b = bdump("dump.1 dump.2.gz") can be gzipped +b = bdump("dump.*") wildcard expands to multiple files +b = bdump("dump.*",0) two args = store filenames, but don't read incomplete and duplicate snapshots are deleted no column name assignment is performed -time = b.next() read next snapshot from dump files +time = b.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 @@ -93,7 +93,7 @@ class bdump: for word in words: self.flist += glob.glob(word) if len(self.flist) == 0 and len(list) == 1: raise StandardError,"no bdump file specified" - + if len(list) == 1: self.increment = 0 self.read_all() @@ -148,15 +148,15 @@ class bdump: snap = self.read_snapshot(f) if not snap: self.nextfile += 1 - if self.nextfile == len(self.flist): return -1 + if self.nextfile == len(self.flist): return -1 f.close() - self.eof = 0 - continue + self.eof = 0 + continue self.eof = f.tell() f.close() try: self.findtime(snap.time) - continue + continue except: break self.snaps.append(snap) @@ -170,7 +170,7 @@ class bdump: # 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() @@ -202,7 +202,7 @@ class bdump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise StandardError, "bdump map() requires pairs of mappings" @@ -242,7 +242,7 @@ class bdump: del self.snaps[i] else: i += 1 - + # -------------------------------------------------------------------- # return list of bonds to viz for snapshot isnap # if called with flag, then index is timestep, so convert to snapshot index @@ -267,7 +267,7 @@ class bdump: # create line list from id,type,atom1,atom2 # abs() of type since could be negative - + bonds = [] for i in xrange(snap.natoms): atom = snap.atoms[i] diff --git a/src/cdata.py b/src/cdata.py index 3020fc5..d74a9d0 100644 --- a/src/cdata.py +++ b/src/cdata.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. # how to specify and print 2d particles @@ -13,11 +13,11 @@ oneline = "Read, create, manipulate ChemCell data files" docstr = """ -c = cdata() create a datafile object +c = cdata() create a datafile object c = cdata("mem.surf") read in one or more ChemCell data files c = cdata("mem.part.gz mem.surf") can be gzipped -c = cdata("mem.*") wildcard expands to multiple files -c.read("mem.surf") read in one or more data files +c = cdata("mem.*") wildcard expands to multiple files +c.read("mem.surf") read in one or more data files read() has same argument options as constructor files contain the following kinds of entries, each of which becomes an object @@ -30,10 +30,10 @@ c.read("mem.surf") read in one or more data files ID can be any number or string, must be unique c.box(ID,xlo,ylo,zlo,xhi,yhi,zhi) create a box region -c.sphere(ID,x,y,z,r) create a sphere region -c.shell(ID,x,y,z,r,rinner) create a shell region -c.cyl(ID,'x',c1,c2,r,lo,hi) create a axis-aligned cylinder region -c.cap(ID,'x',c1,c2,r,lo,hi) create a axis-aligned capped-cylinder region +c.sphere(ID,x,y,z,r) create a sphere region +c.shell(ID,x,y,z,r,rinner) create a shell region +c.cyl(ID,'x',c1,c2,r,lo,hi) create a axis-aligned cylinder region +c.cap(ID,'x',c1,c2,r,lo,hi) create a axis-aligned capped-cylinder region c.q(ID,q1,q2,...) set region triangulation quality factors box() can create an axis-aligned plane, line, or point if lo=hi @@ -57,14 +57,14 @@ c.bins(ID,nx,ny) set binning parameters for a surf for surftri(), one or more tri indices (1-N) must be listed for surfselect(), test is string like "$x < 2.0 and $y > 0.0" bins are used when particles are created inside/outside a surf - -c.part(ID,n,id_in) create N particles inside object id_in -c.part(ID,n,id_in,id_out) particles are also outside object id_out + +c.part(ID,n,id_in) create N particles inside object id_in +c.part(ID,n,id_in,id_out) particles are also outside object id_out c.part2d(ID,n,id_on) create 2d particles on object id_on c.partarray(ID,nx,nz,nz,x,y,z,dx,dy,dz) create 3d grid of particles c.partring(ID,n,x,y,z,r,'x') create ring of particles c.partsurf(ID,id_on) change surf of existing 2d particle group -c.seed(43284) set random # seed (def = 12345) +c.seed(43284) set random # seed (def = 12345) generate particle positions randomly (unless otherwise noted) for part(), id_in and id_out must be IDs of a surf, region, or union object @@ -88,9 +88,9 @@ c.project(ID,ID2,dx,dy,dz,eps,fg) project particles in ID to surf of obj ID2 particles are converted to 2d assigned to surf ID2 c.center(ID,x,y,z) set center point of object -c.trans(ID,dx,dy,dz) translate an object +c.trans(ID,dx,dy,dz) translate an object c.rotate(ID,'x',1,1,0,'z',-1,1,0) rotate an object -c.scale(ID,sx,sy,sz) scale an object +c.scale(ID,sx,sy,sz) scale an object objects must be surface or particle group, regions cannot be changed for center(), default is middle of bounding box (set when obj is created) @@ -98,11 +98,11 @@ c.scale(ID,sx,sy,sz) scale an object object is rotated so that it's current xyz axes point along new ones rotation and scaling occur relative to center point -c.union(ID,id1,id2,...) create a new union object from id1,id2,etc +c.union(ID,id1,id2,...) create a new union object from id1,id2,etc c.join(ID,id1,id2,...) create a new object by joining id1,id2,etc c.delete(id1,id2,...) delete one or more objects c.rename(ID,IDnew) rename an object -c.copy(ID,IDnew) create a new object as copy of old object +c.copy(ID,IDnew) create a new object as copy of old object for union, all lower-level objects must be of surface, region, or union style for join, all joined objects must be of same style: group, surf, line @@ -114,14 +114,14 @@ c.unselect(id1,id2,...) unselect one or more objects c.unselect() unselect all objects selection applies to write() and viz() - -c.write("file") write all selected objs to ChemCell file -c.write("file",id1,id2,...) write only listed & selected objects to file -c.append("file") append all selected objs to ChemCell file -c.append("file",id1,id2,...) append only listed & selected objects + +c.write("file") write all selected objs to ChemCell file +c.write("file",id1,id2,...) write only listed & selected objects to file +c.append("file") append all selected objs to ChemCell file +c.append("file",id1,id2,...) append only listed & selected objects union objects are skipped, not written to file - + index,time,flag = c.iterator(0/1) loop over single snapshot time,box,atoms,bonds,tris,lines = c.viz(index) return list of viz objects @@ -199,7 +199,7 @@ class cdata: else: f = open(file) # read all entries in file - + while 1: line = f.readline() if not line: break @@ -214,7 +214,7 @@ class cdata: raise StandardError, "unrecognized ChemCell data file" # create a surface object from set of triangles or facets - + if flag == "triangles" or flag == "facets": tmp,id,nvert,ntri = line.split() nvert = int(nvert) @@ -242,7 +242,7 @@ class cdata: connections.append([int(value) for value in list[1:]]) else: connections = connect(nvert,ntri,triangles) - + obj = Surface() obj.select = 1 self.ids[id] = len(self.objs) @@ -255,7 +255,7 @@ class cdata: obj.triangles = triangles obj.connections = connections obj.center() - + print id, sys.stdout.flush() @@ -286,7 +286,7 @@ class cdata: obj.npart = npart obj.xyz = xyz obj.center() - + print id, sys.stdout.flush() @@ -297,7 +297,7 @@ class cdata: id = words[1] style = words[2] args = words[3:] - + if style == "box": obj = Box(*args) obj.substyle = BOX @@ -309,7 +309,7 @@ class cdata: self.objs.append(obj) obj.id = id obj.style = REGION - + print id, sys.stdout.flush() @@ -398,7 +398,7 @@ class cdata: cmd = "obj.q%d = arg" % n exec cmd n += 1 - + # -------------------------------------------------------------------- # create a line object with single line @@ -413,7 +413,7 @@ class cdata: obj.style = LINE obj.nline = 0 obj.pairs = [] - + obj.addline(args) # -------------------------------------------------------------------- @@ -430,7 +430,7 @@ class cdata: obj.style = LINE obj.nline = 0 obj.pairs = [] - + xlo,ylo,zlo,xhi,yhi,zhi = args obj.addline([xlo,ylo,zlo,xhi,ylo,zlo]) obj.addline([xlo,yhi,zlo,xhi,yhi,zlo]) @@ -444,7 +444,7 @@ class cdata: obj.addline([xhi,ylo,zlo,xhi,ylo,zhi]) obj.addline([xhi,yhi,zlo,xhi,yhi,zhi]) obj.addline([xlo,yhi,zlo,xlo,yhi,zhi]) - + # -------------------------------------------------------------------- # create a triangulated surface object from a region object @@ -454,7 +454,7 @@ class cdata: region = self.objs[self.ids[id_region]] region.triangulate() - + obj = Surface() obj.select = 1 self.ids[id] = len(self.objs) @@ -477,7 +477,7 @@ class cdata: raise StandardError,"ID %s is already in use" % id o = self.objs[self.ids[id_surf]] - + obj = Surface() obj.select = 1 self.ids[id] = len(self.objs) @@ -490,7 +490,7 @@ class cdata: obj.triangles = [] # subtract 1 from tri and vert to convert to C indexing from (1-N) - + for i in list: v1 = o.triangles[i-1][0] v2 = o.triangles[i-1][1] @@ -503,7 +503,7 @@ class cdata: obj.ntri += 1 # make any connections in new set of triangles - + obj.connections = connect(obj.nvert,obj.ntri,obj.triangles) obj.center() @@ -516,7 +516,7 @@ class cdata: raise StandardError,"ID %s is already in use" % id o = self.objs[self.ids[id_surf]] - + obj = Surface() obj.select = 1 self.ids[id] = len(self.objs) @@ -547,7 +547,7 @@ class cdata: # loop over triangles in id_surf # 3 vertices must satisfy all 3 tests for tri's inclusion in new surf obj - + for tri in o.triangles: v1 = tri[0] - 1 v2 = tri[1] - 1 @@ -600,7 +600,7 @@ class cdata: if out_id: out_obj = self.objs[self.ids[out_id]] # pre-process SURFACE objects to bin their triangles for faster searching - + if in_obj.style == SURFACE: in_obj.inside_prep() if out_id and out_obj.style == SURFACE: out_obj.inside_prep() @@ -612,7 +612,7 @@ class cdata: zsize = zhi-zlo # generate particles until have enough that satisfy in/out constraints - + count = attempt = 0 while count < npart: attempt += 1 @@ -626,7 +626,7 @@ class cdata: obj.center() print "Created %d particles in %d attempts" % (count,attempt) - + # -------------------------------------------------------------------- # create a group object with npart 2d particles on surface of on_id object @@ -649,12 +649,12 @@ class cdata: on_obj.style != UNION: raise StandardError,"Illegal ID to place particles on" totalarea = on_obj.area() - + for count in xrange(npart): area = self.random() * totalarea pt,norm = on_obj.loc2d(area,self.random) obj.xyz.append(pt) - + obj.center() print "Created %d particles on area of %g" % (npart,totalarea) @@ -663,7 +663,7 @@ class cdata: # array size = Nx by Ny by Nz # lower left corner of array = x,y,z # array spacing = dx,dy,dz - + def partarray(self,id,nx,ny,nz,x,y,z,dx,dy,dz): if self.ids.has_key(id): raise StandardError,"ID %s is already in use" % id @@ -685,7 +685,7 @@ class cdata: for i in xrange(nx): xnew = x + i*dx obj.xyz.append([xnew,ynew,znew]) - + obj.center() print "Created %d particles" % (nx*ny*nz) @@ -693,7 +693,7 @@ class cdata: # create a ring of N particles # ring center and radius = x,y,z,r # ring axis = 'x' or 'y' or 'z' - + def partring(self,id,n,x,y,z,r,axis): if self.ids.has_key(id): raise StandardError,"ID %s is already in use" % id @@ -723,7 +723,7 @@ class cdata: ynew = y + r * sin(i*deltheta) znew = z obj.xyz.append([xnew,ynew,znew]) - + obj.center() print "Created %d particles" % n @@ -779,19 +779,19 @@ class cdata: # move along dir until get within EPS of surf # factor = multiply bracketing distance by this amount each iteration # maxscale = max multiple of dir vector to bracket in each direction - + factor = 2 maxscale = 10.0 - + for i in xrange(obj.npart): x,y,z = obj.xyz[i] if flag: dir = [dx-x,dy-y,dz-z] else: dir = [dx,dy,dz] normalize(dir) - + # start = in/out at starting pt # stop = in/out at bracketing pt - + start = obj_on.inside(x,y,z) if start: stop = 0 else: stop = 1 @@ -800,8 +800,6 @@ class cdata: # bracket pt = xyz +/- scale*dir # multiply scale by factor each iteration - #print "AAA",i,start,stop,x,y,z,dir - scale = EPS bracket = start while scale < maxscale: @@ -820,7 +818,7 @@ class cdata: # bisection search to zoom in to within EPS of surface # separation = distance between 2 points - + delx = xnew-x; dely = ynew-y; delz = znew-z separation = sqrt(delx*delx + dely*dely + delz*delz) while separation > EPS: @@ -838,7 +836,7 @@ class cdata: obj.xyz[i][0] = x obj.xyz[i][1] = y obj.xyz[i][2] = z - + obj.on_id = id2 obj.center() @@ -891,7 +889,7 @@ class cdata: raise StandardError,"Can only use rotate() on a surface or group object" # create xyz old and new - + xnew = ynew = znew = None if axis1 == 'x': xnew = [i1,j1,k1] elif axis1 == 'y': ynew = [i1,j1,k1] @@ -952,7 +950,7 @@ class cdata: obj.xyz[i][0] = obj.xc + sx * (obj.xyz[i][0] - obj.xc) obj.xyz[i][1] = obj.yc + sy * (obj.xyz[i][1] - obj.yc) obj.xyz[i][2] = obj.zc + sz * (obj.xyz[i][2] - obj.zc) - + # -------------------------------------------------------------------- # create union object from list of other objects @@ -969,7 +967,7 @@ class cdata: # -------------------------------------------------------------------- # join objects in list to form a new object # only possible for group, surface, line objects - + def join(self,id,*list): style = self.objs[self.ids[list[0]]].style if style == GROUP: obj = Group() @@ -994,14 +992,14 @@ class cdata: elif style == LINE: obj.nline = 0 obj.pairs = [] - + for id in list: o = self.objs[self.ids[id]] if o.style != style: raise StandardError,"All joined objects must be of same style" # force deep copy of particle coords - + if style == GROUP: if o.on_id != obj.on_id: raise StandardError,"Particle group surfaces do not match" @@ -1010,10 +1008,10 @@ class cdata: obj.xyz.append(xyz) obj.npart += o.npart obj.center() - + # force deep copy of triangle vertices and indices # increment vertex and triangle indices b/c now have previous surfaces - + elif style == SURFACE: for i in xrange(o.nvert): vert = o.vertices[i][:] @@ -1031,7 +1029,7 @@ class cdata: obj.center() # force deep copy of line pt pairs - + elif style == LINE: pairs = o.pairs[:] obj.pairs += (pairs) @@ -1040,7 +1038,7 @@ class cdata: # -------------------------------------------------------------------- # delete each object in list # reset values in ids since some indices are decremented - + def delete(self,*list): for id in list: i = self.ids[id] @@ -1049,7 +1047,7 @@ class cdata: for key in self.ids.keys(): j = self.ids[key] if j > i: self.ids[key] = j-1 - + # -------------------------------------------------------------------- # rename the ID of an object # check that new name doesn't already exist @@ -1061,7 +1059,7 @@ class cdata: self.ids[idnew] = i self.objs[i].id = idnew del self.ids[idold] - + # -------------------------------------------------------------------- # create a deep copy of an object and assign it a new ID # check that new name doesn't already exist @@ -1078,7 +1076,7 @@ class cdata: # -------------------------------------------------------------------- # set selection flag for each object in list # if list is empty, select all - + def select(self,*list): if len(list) == 0: list = self.ids.keys() for id in list: @@ -1088,7 +1086,7 @@ class cdata: # -------------------------------------------------------------------- # unset selection flag for each object in list # if list is empty, unselect all - + def unselect(self,*list): if len(list) == 0: list = self.ids.keys() for id in list: @@ -1098,7 +1096,7 @@ class cdata: # -------------------------------------------------------------------- # write out list of objects to data file via filewrite() # open a new file - + def write(self,file,*list): if not len(list): vlist = range(len(self.objs)) else: @@ -1108,11 +1106,11 @@ class cdata: fp = open(file,'w') self.filewrite(fp,vlist) fp.close() - + # -------------------------------------------------------------------- # append list of objects to data file via filewrite() # open existing file for appending - + def append(self,file,*list): if not len(list): vlist = range(len(self.objs)) else: @@ -1126,7 +1124,7 @@ class cdata: # -------------------------------------------------------------------- # write out list of objects to previously opened data file # for particles, write as 3d or 2d depending on on_id - + def filewrite(self,fp,vlist): for index in vlist: obj = self.objs[index] @@ -1166,12 +1164,12 @@ class cdata: return 0,0,-1 # -------------------------------------------------------------------- - # return list of atoms and triangles to viz for cdata object + # return list of atoms and triangles and lines to viz for cdata object def viz(self,isnap): if isnap: raise StandardError, "cannot call cdata.viz() with isnap != 0" - + # create atom list from sum of all particle groups # id = running count # type = running type of particle group @@ -1187,7 +1185,7 @@ class cdata: atoms.append([id,itype,xyz[0],xyz[1],xyz[2]]) # no bonds - + bonds = [] # create triangle list from sum of all surfaces and regions @@ -1221,7 +1219,7 @@ class cdata: for pair in obj.pairs: id += 1 lines.append([id,itype] + pair) - + return 0,self.bbox(),atoms,bonds,tris,lines # -------------------------------------------------------------------- @@ -1285,7 +1283,7 @@ IR = 2836 class Random: def __init__(self,seed): self.seed = seed - + def __call__(self): k = self.seed/IQ self.seed = IA*(self.seed-k*IQ) - IR*k @@ -1299,9 +1297,9 @@ class Surface: def __init__(self): self.nbinx = self.nbiny = 0 - + # bounding box - + def bbox(self): list = [float(vert[0]) for vert in self.vertices] xlo = min(list) @@ -1316,7 +1314,7 @@ class Surface: return (xlo,ylo,zlo,xhi,yhi,zhi) # set center point explicitly or set to middle of bounding box - + def center(self,*xyz): if len(xyz): self.xc = xyz[0] @@ -1332,7 +1330,7 @@ class Surface: # each bin overlayed by triangle's xy bounding box stores the triangle index # EPSILON insures that bins completely overlay surface bounding box # self.xlo, self.ylo, self.dxinv, self.dyinv, self.bin are used by inside() - + def inside_prep(self): self.xlo,self.ylo,self.zlo,self.xhi,self.yhi,self.zhi = self.bbox() xsize = self.xhi - self.xlo @@ -1348,13 +1346,13 @@ class Surface: print "Binning %d triangles into %d by %d bins ..." % \ (self.ntri,self.nbinx,self.nbiny) - + self.bin = [] for i in xrange(self.nbinx): self.bin.append(self.nbiny*[0]) for j in xrange(self.nbiny): self.bin[i][j] = [] - + for m in xrange(self.ntri): v1 = self.vertices[self.triangles[m][0]-1] v2 = self.vertices[self.triangles[m][1]-1] @@ -1370,7 +1368,7 @@ class Surface: for i in xrange(ilo,ihi+1): for j in xrange(jlo,jhi+1): self.bin[i][j].append(m) - + print "Done with binning" # check for inside assumes that surf is a closed set of triangles @@ -1384,13 +1382,13 @@ class Surface: # is xy pt inside tri in xy plane (including edges and vertices) ? # compute ztri = z value of xy pt on plane of 3d tri # if z of pt is <= ztri, then line segment intersects the tri - + def inside(self,x,y,z): ix = int((x - self.xlo) * self.dxinv) iy = int((y - self.ylo) * self.dyinv) if ix < 0 or ix >= self.nbinx or iy < 0 or iy >= self.nbiny: return 0 n = len(self.bin[ix][iy]) - + hit = 0 for m in xrange(n): itri = self.bin[ix][iy][m] @@ -1399,9 +1397,9 @@ class Surface: v1 = self.vertices[tri[0]-1] v2 = self.vertices[tri[1]-1] v3 = self.vertices[tri[2]-1] - + # is x,y in bounding box of 2d triangle ? - + if x < v1[0] and x < v2[0] and x < v3[0]: continue if x > v1[0] and x > v2[0] and x > v3[0]: continue if y < v1[1] and y < v2[1] and y < v3[1]: continue @@ -1410,7 +1408,7 @@ class Surface: # is x,y inside 2d triangle ? # cross product of each edge with vertex-to-point must have same sign # cross product = 0 is OK, means point is on edge or vertex - + c1 = (v2[0]-v1[0])*(y-v1[1]) - (v2[1]-v1[1])*(x-v1[0]) c2 = (v3[0]-v2[0])*(y-v2[1]) - (v3[1]-v2[1])*(x-v2[0]) c3 = (v1[0]-v3[0])*(y-v3[1]) - (v1[1]-v3[1])*(x-v3[0]) @@ -1429,7 +1427,7 @@ class Surface: # denom (vx - vy * wx/xy) cannot be 0 since would imply vx/vy = wx/wy # and thus e1 would be parallel to e2 # if wy = 0, vy cannot be 0, since e1 would be parallel to e2 - + px = x - v1[0]; py = y - v1[1] vx = v2[0] - v1[0] @@ -1448,7 +1446,7 @@ class Surface: # surface area # areas = cummulative total area of all triangles # triangle area = 1/2 of magnitude of cross product of 2 edge vectors - + def area(self): a = 3*[0] self.areas = [] @@ -1465,7 +1463,7 @@ class Surface: return sum # return a random location on one of triangles - + def loc2d(self,area,random): for i in xrange(self.ntri): if area < self.areas[i]: break @@ -1502,7 +1500,7 @@ class Group: return (xlo,ylo,zlo,xhi,yhi,zhi) # set center point explicitly or set to middle of bounding box - + def center(self,*xyz): if len(xyz): self.xc = xyz[0] @@ -1528,7 +1526,7 @@ class Box: self.q1 = self.q2 = self.q3 = 1 # bounding box around region - + def bbox(self): return (self.xlo,self.ylo,self.zlo,self.xhi,self.yhi,self.zhi) @@ -1544,7 +1542,7 @@ class Box: # set nvert,ntri,vertices,triangles,connections # convert vertices from unit box to lo/hi box # convert tuples returned by box_triangluate (used for dict lookup) to lists - + def triangulate(self): vertices,triangles = box_triangulate(self.q1,self.q2,self.q3) self.nvert = len(vertices) @@ -1559,9 +1557,9 @@ class Box: for i in xrange(self.ntri): self.triangles.append([triangles[i][0],triangles[i][1],triangles[i][2]]) self.connections = connect(self.nvert,self.ntri,self.triangles) - + # surface area of region - + def area(self): xsize = self.xhi - self.xlo ysize = self.yhi - self.ylo @@ -1582,7 +1580,7 @@ class Box: return sum # return a random location on surface of region - + def loc2d(self,area,random): xsize = self.xhi - self.xlo ysize = self.yhi - self.ylo @@ -1627,7 +1625,7 @@ class Sphere: self.x+self.r,self.y+self.r,self.z+self.r) # return 1,0 for xyz inside/outside the region - + def inside(self,x,y,z): dx = x - self.x dy = y - self.y @@ -1661,15 +1659,15 @@ class Sphere: for i in xrange(self.ntri): self.triangles.append([triangles[i][0],triangles[i][1],triangles[i][2]]) self.connections = connect(self.nvert,self.ntri,self.triangles) - + # surface area of region - + def area(self): value = 4 * pi * self.r*self.r return value # return a random location on surface of region - + def loc2d(self,area,random): while 1: x = random() - 0.5 @@ -1681,7 +1679,7 @@ class Sphere: return [self.x + self.r*c[0], self.y + self.r*c[1], self.z + self.r*c[2]],c # ChemCell text to create the region - + def command(self): return "%s sphere %g %g %g %g" % (self.id,self.x,self.y,self.z,self.r) @@ -1694,7 +1692,7 @@ class Shell(Sphere): def __init__(self,*list): self.rinner = float(list[4]) self.innersq = self.rinner*self.rinner - + def inside(self,x,y,z): dx = x - self.x dy = y - self.y @@ -1877,7 +1875,7 @@ class Capped: elif self.axis == 'z': return (self.c1-self.r,self.c2-self.r,self.lo-self.r, \ self.c1+self.r,self.c2+self.r,self.hi+self.r) - + # return 1,0 for xyz inside/outside the region def inside(self,x,y,z): @@ -2049,7 +2047,7 @@ class Line: def addline(self,coords): self.pairs.append(list(coords)) - + # -------------------------------------------------------------------- # union object that contains other objects # pass in parent cdata IDs and object list to constructor @@ -2080,7 +2078,7 @@ class Union: # return 1,0 for xyz inside/outside the union # inside union if inside any of its child objects - + def inside(self,x,y,z): for obj in self.objs: if obj.inside(x,y,z): return 1 @@ -2088,7 +2086,7 @@ class Union: # surface area of union # areas = cummulative total area for child objects - + def area(self): self.areas = [] sum = 0.0 @@ -2098,13 +2096,13 @@ class Union: return sum # return a random location on surface of one of child objects - + def loc2d(self,area,random): for i in xrange(len(self.objs)): if area < self.areas[i]: break if i > 0: area -= self.areas[i] return self.objs[i].loc2d(area,random) - + # -------------------------------------------------------------------- # return c = a x b @@ -2161,10 +2159,10 @@ def connect(nvert,ntri,triangles): v2tri = [] for i in xrange(nvert): v2tri.append([]); - + for i in xrange(ntri): for vert in triangles[i]: v2tri[vert-1].append(i) - + # loop over triangles to reset connections connections = [] @@ -2174,17 +2172,17 @@ def connect(nvert,ntri,triangles): connect = 6*[0] connections.append(connect) - + # tri123 = list of tris attached to each vertex of triangle i - + v = triangles[i] tri1 = v2tri[v[0]-1] tri2 = v2tri[v[1]-1] tri3 = v2tri[v[2]-1] - + # loop over attached tris and look for 2nd vertex in each of 3 edges # when find it, set connection triangle and edge values - + for itri in tri1: if itri == i: continue if v[1] in triangles[itri]: diff --git a/src/cfg.py b/src/cfg.py index 8cefd38..bbd930e 100644 --- a/src/cfg.py +++ b/src/cfg.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. # cfg tool @@ -11,14 +11,14 @@ oneline = "Convert LAMMPS snapshots to AtomEye CFG format" docstr = """ -c = cfg(d) d = object containing atom coords (dump, data) +c = cfg(d) d = object containing atom coords (dump, data) c.one() write all snapshots to tmp.cfg c.one("new") write all snapshots to new.cfg c.many() write snapshots to tmp0000.cfg, tmp0001.cfg, etc c.many("new") write snapshots to new0000.cfg, new0001.cfg, etc -c.single(N) write snapshot for timestep N to tmp.cfg -c.single(N,"file") write snapshot for timestep N to file.cfg +c.single(N) write snapshot for timestep N to tmp.cfg +c.single(N,"file") write snapshot for timestep N to file.cfg """ # History @@ -46,7 +46,7 @@ class cfg: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def one(self,*args): @@ -68,16 +68,16 @@ class cfg: print >>f,"Number of particles = %d " % len(atoms) print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time print >>f,"H0(1,1) = %20.10f A " % xlen - print >>f,"H0(1,2) = 0.0 A " - print >>f,"H0(1,3) = 0.0 A " - print >>f,"H0(2,1) = 0.0 A " + print >>f,"H0(1,2) = 0.0 A " + print >>f,"H0(1,3) = 0.0 A " + print >>f,"H0(2,1) = 0.0 A " print >>f,"H0(2,2) = %20.10f A " % ylen - print >>f,"H0(2,3) = 0.0 A " - print >>f,"H0(3,1) = 0.0 A " - print >>f,"H0(3,2) = 0.0 A " + print >>f,"H0(2,3) = 0.0 A " + print >>f,"H0(3,1) = 0.0 A " + print >>f,"H0(3,2) = 0.0 A " print >>f,"H0(3,3) = %20.10f A " % zlen print >>f,"#" - + for atom in atoms: itype = int(atom[1]) xfrac = (atom[2]-box[0])/xlen @@ -85,14 +85,14 @@ class cfg: zfrac = (atom[4]-box[2])/zlen # print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]) print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac) - + print time, sys.stdout.flush() n += 1 - + f.close() print "\nwrote %d snapshots to %s in CFG format" % (n,file) - + # -------------------------------------------------------------------- def many(self,*args): @@ -104,7 +104,7 @@ class cfg: which,time,flag = self.data.iterator(flag) if flag == -1: break time,box,atoms,bonds,tris,lines = self.data.viz(which) - + if n < 10: file = root + "000" + str(n) elif n < 100: @@ -112,7 +112,7 @@ class cfg: elif n < 1000: file = root + "0" + str(n) else: - file = root + str(n) + file = root + str(n) file += ".cfg" f = open(file,"w") @@ -123,16 +123,16 @@ class cfg: print >>f,"Number of particles = %d " % len(atoms) print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time print >>f,"H0(1,1) = %20.10f A " % xlen - print >>f,"H0(1,2) = 0.0 A " - print >>f,"H0(1,3) = 0.0 A " - print >>f,"H0(2,1) = 0.0 A " + print >>f,"H0(1,2) = 0.0 A " + print >>f,"H0(1,3) = 0.0 A " + print >>f,"H0(2,1) = 0.0 A " print >>f,"H0(2,2) = %20.10f A " % ylen - print >>f,"H0(2,3) = 0.0 A " - print >>f,"H0(3,1) = 0.0 A " - print >>f,"H0(3,2) = 0.0 A " + print >>f,"H0(2,3) = 0.0 A " + print >>f,"H0(3,1) = 0.0 A " + print >>f,"H0(3,2) = 0.0 A " print >>f,"H0(3,3) = %20.10f A " % zlen print >>f,"#" - + for atom in atoms: itype = int(atom[1]) xfrac = (atom[2]-box[0])/xlen @@ -140,14 +140,14 @@ class cfg: zfrac = (atom[4]-box[2])/zlen # print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]) print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac) - + print time, sys.stdout.flush() f.close() n += 1 - + print "\nwrote %s snapshots in CFG format" % n - + # -------------------------------------------------------------------- def single(self,time,*args): @@ -166,16 +166,16 @@ class cfg: print >>f,"Number of particles = %d " % len(atoms) print >>f,"# Timestep %d \n#\nA = 1.0 Angstrom" % time print >>f,"H0(1,1) = %20.10f A " % xlen - print >>f,"H0(1,2) = 0.0 A " - print >>f,"H0(1,3) = 0.0 A " - print >>f,"H0(2,1) = 0.0 A " + print >>f,"H0(1,2) = 0.0 A " + print >>f,"H0(1,3) = 0.0 A " + print >>f,"H0(2,1) = 0.0 A " print >>f,"H0(2,2) = %20.10f A " % ylen - print >>f,"H0(2,3) = 0.0 A " - print >>f,"H0(3,1) = 0.0 A " - print >>f,"H0(3,2) = 0.0 A " + print >>f,"H0(2,3) = 0.0 A " + print >>f,"H0(3,1) = 0.0 A " + print >>f,"H0(3,2) = 0.0 A " print >>f,"H0(3,3) = %20.10f A " % zlen print >>f,"#" - + for atom in atoms: itype = int(atom[1]) xfrac = (atom[2]-box[0])/xlen @@ -183,5 +183,5 @@ class cfg: zfrac = (atom[4]-box[2])/zlen # print >>f,"1.0 %d %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f " % (itype,xfrac,yfrac,zfrac,atom[5],atom[6],atom[7]) print >>f,"1.0 %d %15.10f %15.10f %15.10f 0.0 0.0 0.0 " % (itype,xfrac,yfrac,zfrac) - + f.close() diff --git a/src/chain.py b/src/chain.py index cc6a1c2..bd795b6 100644 --- a/src/chain.py +++ b/src/chain.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. # chain tool @@ -12,11 +12,11 @@ oneline = "Create bead-spring chains for LAMMPS input" docstr = """ c = chain(N,rho) setup box with N monomers at reduced density rho -c = chain(N,rho,1,1,2) x,y,z = aspect ratio of box (def = 1,1,1) +c = chain(N,rho,1,1,2) x,y,z = aspect ratio of box (def = 1,1,1) c.seed = 48379 set random # seed (def = 12345) -c.mtype = 2 set type of monomers (def = 1) -c.btype = 1 set type of bonds (def = 1) +c.mtype = 2 set type of monomers (def = 1) +c.btype = 1 set type of bonds (def = 1) c.blen = 0.97 set length of bonds (def = 0.97) c.dmin = 1.02 set min dist from i-1 to i+1 site (def = 1.02) @@ -24,11 +24,11 @@ c.id = "chain" set molecule ID to chain # (default) c.id = "end1" set molecule ID to count from one end of chain c.id = "end2" set molecule ID to count from either end of chain -c.build(100,10) create 100 chains, each of length 10 +c.build(100,10) create 100 chains, each of length 10 can be invoked multiple times interleaved with different settings must fill box with total of N monomers - + c.write("data.file") write out all built chains to LAMMPS data file """ @@ -60,7 +60,7 @@ from data import data # Class definition class chain: - + # -------------------------------------------------------------------- def __init__(self,n,rhostar,*list): @@ -112,7 +112,7 @@ class chain: x = self.xlo + self.random()*self.xprd y = self.ylo + self.random()*self.yprd z = self.zlo + self.random()*self.zprd - ix = iy = iz = 0 + ix = iy = iz = 0 else: restriction = True while restriction: @@ -134,7 +134,7 @@ class chain: dz = z - atoms[-2][5] if math.sqrt(dx*dx + dy*dy + dz*dz) <= self.dmin: restriction = True - + x,y,z,ix,iy,iz = self.pbc(x,y,z,ix,iy,iz) idatom = id_atom_prev + imonomer + 1 if self.id == "chain": @@ -147,12 +147,12 @@ class chain: idmol = nper - imonomer else: raise StandardError,"chain ID is not a valid value" - + atoms.append([idatom,idmol,self.mtype,x,y,z,ix,iy,iz]) if imonomer: - bondid = id_bond_prev + imonomer + bondid = id_bond_prev + imonomer bonds.append([bondid,self.btype,idatom-1,idatom]) - + self.atoms += atoms self.bonds += bonds @@ -186,7 +186,7 @@ class chain: lines = [] for i in range(atypes): lines.append("%d 1.0\n" % (i+1)) d.sections["Masses"] = lines - + lines = [] for atom in self.atoms: line = "%d %d %d %g %g %g %d %d %d\n" % \ @@ -194,7 +194,7 @@ class chain: atom[6], atom[7], atom[8]) lines.append(line) d.sections["Atoms"] = lines - + lines = [] for bond in self.bonds: line = "%d %d %d %d\n" % (bond[0], bond[1], bond[2], bond[3]) diff --git a/src/clog.py b/src/clog.py index 2137b84..432caa5 100644 --- a/src/clog.py +++ b/src/clog.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. # clog tool @@ -21,12 +21,12 @@ c = clog("log.cell","",0) 3rd arg = average all runs if specify 2nd arg, it delimits a time section no 2nd arg or empty string, use default which is ChemCell specific if specify any 3rd arg, average all runs, assume all start at time 0 - + nvec = c.nvec # of vectors of thermo info nlen = c.nlen length of each vectors names = c.names list of vector names a,b,... = c.get("A","B",...) return one or more vectors of values -c.write("file.txt") write all vectors to a file +c.write("file.txt") write all vectors to a file c.write("file.txt","A","B",...) write listed vectors to a file get and write allow abbreviated (uniquely) vector names @@ -78,12 +78,12 @@ class clog: if len(list) > 1 and len(list[1]): self.firststr = list[1] if len(list) == 3: self.ave = 1 - + self.read_all() # -------------------------------------------------------------------- # read all log data from all files - + def read_all(self): self.read_header(self.flist[0]) if self.nvec == 0: raise StandardError,"log file has no values" @@ -117,9 +117,9 @@ class clog: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: @@ -145,9 +145,9 @@ class clog: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: @@ -205,7 +205,7 @@ class clog: self.nlen = nlen self.data = data - + # -------------------------------------------------------------------- def read_header(self,file): @@ -261,38 +261,38 @@ class clog: 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 + last = 1 s1 = 0 elif s1 >= 0 and s2 == -1: # found s1, but no s2 last = 1 s1 = txt.find("\n",s1) + 1 s2 = txt.rfind("\n",s1) + 1 - eof -= len(txt) - s2 + eof -= len(txt) - s2 elif s1 == -1 and s2 == -1: # found neither # could be end-of-file section - # or entire read was one chunk + # 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 + eof -= len(txt) - start # reset eof to "Loop" + break - last = 1 # entire read is a chunk + last = 1 # entire read is a chunk s1 = 0 s2 = txt.rfind("\n",s1) + 1 - eof -= len(txt) - s2 - if s1 == s2: break + 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 - + 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]), diff --git a/src/data.py b/src/data.py index b2c9344..5e2277b 100644 --- a/src/data.py +++ b/src/data.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. # data tool @@ -12,7 +12,7 @@ oneline = "Read, write, manipulate LAMMPS data files" docstr = """ d = data("data.poly") read a LAMMPS data file, can be gzipped -d = data() create an empty data file +d = data() create an empty data file d.map(1,"id",3,"x") assign names to atom columns (1-N) @@ -26,17 +26,17 @@ d.reorder("Atoms",1,3,2,4,5) reorder columns (1-N) in a data file section 1,3,2,4,5 = new order of previous columns, can delete columns this way -d.title = "My LAMMPS data file" set title of the data file +d.title = "My LAMMPS data file" set title of the data file d.headers["atoms"] = 1500 set a header value d.sections["Bonds"] = lines set a section to list of lines (with newlines) -d.delete("bonds") delete a keyword or section of data file +d.delete("bonds") delete a keyword or section of data file d.delete("Bonds") -d.replace("Atoms",5,vec) replace Nth column of section with vector -d.newxyz(dmp,1000) replace xyz in Atoms with xyz of snapshot N +d.replace("Atoms",5,vec) replace Nth column of section with vector +d.newxyz(dmp,1000) replace xyz in Atoms with xyz of snapshot N newxyz assumes id,x,y,z are defined in both data and dump files also replaces ix,iy,iz if they are defined - + index,time,flag = d.iterator(0/1) loop over single data file snapshot time,box,atoms,bonds,tris,lines = d.viz(index) return list of viz objects @@ -53,7 +53,7 @@ time,box,atoms,bonds,tris,lines = d.viz(index) return list of viz objects NULL if bonds do not exist tris = NULL lines = NULL - + d.write("data.new") write a LAMMPS data file """ @@ -65,7 +65,7 @@ d.write("data.new") write a LAMMPS data file # Variables # title = 1st line of data file -# names = dictionary with atom attributes as keys, col #s as values +# names = dictionary with atom attributes as keys, col #s as values # headers = dictionary with header name as key, value or tuple as values # sections = dictionary with section name as key, array of lines as values # nselect = 1 = # of snapshots @@ -85,7 +85,7 @@ class data: def __init__(self,*list): self.nselect = 1 - + if len(list) == 0: self.title = "LAMMPS data file" self.names = {} @@ -99,7 +99,7 @@ class data: self.title = f.readline() self.names = {} - + headers = {} while 1: line = f.readline() @@ -109,16 +109,16 @@ class data: found = 0 for keyword in hkeywords: if line.find(keyword) >= 0: - found = 1 - words = line.split() - if keyword == "xlo xhi" or keyword == "ylo yhi" or \ - keyword == "zlo zhi": - headers[keyword] = (float(words[0]),float(words[1])) - elif keyword == "xy xz yz": - headers[keyword] = \ + found = 1 + words = line.split() + if keyword == "xlo xhi" or keyword == "ylo yhi" or \ + keyword == "zlo zhi": + headers[keyword] = (float(words[0]),float(words[1])) + elif keyword == "xy xz yz": + headers[keyword] = \ (float(words[0]),float(words[1]),float(words[2])) else: - headers[keyword] = int(words[0]) + headers[keyword] = int(words[0]) if not found: break @@ -128,11 +128,11 @@ class data: for pair in skeywords: keyword,length = pair[0],pair[1] if keyword == line: - found = 1 + found = 1 if not headers.has_key(length): raise StandardError, \ "data section %s has no matching header value" % line - f.readline() + f.readline() list = [] for i in xrange(headers[length]): list.append(f.readline()) sections[keyword] = list @@ -143,7 +143,7 @@ class data: if not line: break line = line.strip() - + f.close() self.headers = headers self.sections = sections @@ -225,21 +225,21 @@ class data: if dm.scaled(nsnap): scaleflag = 1 else: scaleflag = 0 dm.sort(ntime) - if scaleflag: dm.unscale(ntime) + x,y,z = dm.vecs(ntime,"x","y","z") if scaleflag: dm.scale(ntime) self.replace("Atoms",self.names['x']+1,x) self.replace("Atoms",self.names['y']+1,y) self.replace("Atoms",self.names['z']+1,z) - + if dm.names.has_key("ix") and self.names.has_key("ix"): ix,iy,iz = dm.vecs(ntime,"ix","iy","iz") self.replace("Atoms",self.names['ix']+1,ix) self.replace("Atoms",self.names['iy']+1,iy) self.replace("Atoms",self.names['iz']+1,iz) - + # -------------------------------------------------------------------- # delete header value or section from data file @@ -259,19 +259,19 @@ class data: if self.headers.has_key(keyword): if keyword == "xlo xhi" or keyword == "ylo yhi" or \ keyword == "zlo zhi": - pair = self.headers[keyword] - print >>f,pair[0],pair[1],keyword + pair = self.headers[keyword] + print >>f,pair[0],pair[1],keyword elif keyword == "xy xz yz": - triple = self.headers[keyword] - print >>f,triple[0],triple[1],triple[2],keyword + triple = self.headers[keyword] + print >>f,triple[0],triple[1],triple[2],keyword else: - print >>f,self.headers[keyword],keyword + print >>f,self.headers[keyword],keyword for pair in skeywords: keyword = pair[0] if self.sections.has_key(keyword): print >>f,"\n%s\n" % keyword for line in self.sections[keyword]: - print >>f,line, + print >>f,line, f.close() # -------------------------------------------------------------------- @@ -287,19 +287,19 @@ class data: def findtime(self,n): if n == 0: return 0 raise StandardError, "no step %d exists" % (n) - + # -------------------------------------------------------------------- # return list of atoms and bonds to viz for data object def viz(self,isnap): if isnap: raise StandardError, "cannot call data.viz() with isnap != 0" - + id = self.names["id"] type = self.names["type"] x = self.names["x"] y = self.names["y"] z = self.names["z"] - + xlohi = self.headers["xlo xhi"] ylohi = self.headers["ylo yhi"] zlohi = self.headers["zlo zhi"] @@ -331,8 +331,8 @@ class data: float(atom1words[z]), float(atom2words[x]),float(atom2words[y]), float(atom2words[z]), - float(atom1words[type]),float(atom2words[type])]) - + float(atom1words[type]),float(atom2words[type])]) + tris = [] lines = [] return 0,box,atoms,bonds,tris,lines @@ -357,18 +357,18 @@ class data: hkeywords = ["atoms","lines","tris", "bonds","angles","dihedrals","impropers", - "atom types","bond types","angle types","dihedral types", - "improper types","xlo xhi","ylo yhi","zlo zhi","xy xz yz"] + "atom types","bond types","angle types","dihedral types", + "improper types","xlo xhi","ylo yhi","zlo zhi","xy xz yz"] skeywords = [["Masses","atom types"], ["Atoms","atoms"],["Lines","lines"],["Triangles","tris"], ["Bonds","bonds"], - ["Angles","angles"],["Dihedrals","dihedrals"], - ["Impropers","impropers"],["Velocities","atoms"], + ["Angles","angles"],["Dihedrals","dihedrals"], + ["Impropers","impropers"],["Velocities","atoms"], ["Pair Coeffs","atom types"], - ["Bond Coeffs","bond types"],["Angle Coeffs","angle types"], - ["Dihedral Coeffs","dihedral types"], - ["Improper Coeffs","improper types"], + ["Bond Coeffs","bond types"],["Angle Coeffs","angle types"], + ["Dihedral Coeffs","dihedral types"], + ["Improper Coeffs","improper types"], ["BondBond Coeffs","angle types"], ["BondAngle Coeffs","angle types"], ["MiddleBondTorsion Coeffs","dihedral types"], diff --git a/src/dump.py b/src/dump.py index f700cb8..5f5f892 100644 --- a/src/dump.py +++ b/src/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. # dump tool @@ -12,15 +12,15 @@ 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 -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 @@ -31,21 +31,21 @@ 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 -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 +56,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 - headd = 0/1 for no/yes snapshot header, app = 0/1 for write vs append + head = 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 +85,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,7 +107,7 @@ 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 @@ -117,7 +117,7 @@ fx,fy,... = d.vecs(1000,"fx","fy",...) return vector(s) for timestep N 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 +d.extra(obj) extract bond/tri/line info from obj iterator() loops over selected timesteps iterator() called with arg = 0 first time, with arg = 1 on subsequent calls @@ -137,7 +137,7 @@ d.extra(obj) extract bond/tri/line info from obj if extra() used to define lines, else NULL 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, + obj can be data object for bonds, cdata object for tris and lines, bdump object for bonds, mdump object for tris, ldump object for lines """ @@ -204,7 +204,7 @@ class dump: # -------------------------------------------------------------------- def __init__(self,*input,**kwargs): - + self.snaps = [] self.nsnaps = self.nselect = 0 self.names = {} @@ -220,16 +220,16 @@ class dump: self.multiprocflag = 0 self.fileNums = [] self.objextra = None - + outputfl = True if isinstance(input[0],dict): # multiprocessing code (the [0] comes from the asteriks in the argumentlist) dictionary = input[0] - + # check whether to output or not if "debugMode" in dictionary: outputfl = dictionary["debugMode"] - + if outputfl: print "number of subprocess:", os.getpid() - + self.flist = dictionary["filelist"] self.multiprocflag = 1 self.increment = 0 @@ -241,7 +241,7 @@ class dump: for word in words: self.flist += glob.glob(word) if len(self.flist) == 0 and len(input) == 1: raise StandardError,"no dump file specified" - + if len(input) == 1: self.increment = 0 self.read_all(output=outputfl) @@ -256,11 +256,11 @@ class dump: # read all snapshots from each file # test for gzipped files - + # check whether to output or not outputfl = True if "output" in kwargs: outputfl = kwargs["output"] - + if outputfl: print "reading dump file..." for i, file in enumerate(self.flist): @@ -292,7 +292,7 @@ class dump: self.tselect.all(output=outputfl) # set default names for atom columns if file wasn't self-describing - + if len(self.snaps) == 0: if outputfl: print "no column assignments made" elif len(self.names): @@ -300,7 +300,7 @@ class dump: else: if outputfl: print "no column assignments made" pass - + # if snapshots are scaled, unscale them if (not self.names.has_key("x")) or \ @@ -331,15 +331,15 @@ class dump: snap = self.read_snapshot(f) if not snap: self.nextfile += 1 - if self.nextfile == len(self.flist): return -1 + if self.nextfile == len(self.flist): return -1 f.close() - self.eof = 0 - continue + self.eof = 0 + continue self.eof = f.tell() f.close() try: self.findtime(snap.time) - continue + continue except: break # select the new snapshot with all its atoms @@ -360,7 +360,7 @@ class dump: # assign column names if not already done and file is self-describing # set scale_original to 0/1/-1 for unscaled/scaled/unknown # convert xs,xu to x - + def read_snapshot(self,f): try: snap = Snap() @@ -407,7 +407,7 @@ class dump: 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() ncol = len(words) @@ -430,7 +430,7 @@ class dump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise StandardError, "dump map() requires pairs of mappings" @@ -507,7 +507,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 @@ -520,7 +520,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: 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 @@ -555,10 +555,10 @@ 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 ..." - + id = self.names["id"] x = self.names["x"] y = self.names["y"] @@ -584,10 +584,10 @@ class dump: # 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 = max(self.names.values()) #len(self.snaps[0].atoms[0]) pairs = self.names.items() @@ -603,11 +603,11 @@ class dump: # if arg = numeric, sort atoms in single step def sort(self,*list, **kwargs): - + # check whether to output or not outputfl = True if "output" in kwargs: outputfl = kwargs["output"] - + if len(list) == 0: if outputfl: print "Sorting selected snapshots ..." id = self.names["id"] @@ -655,7 +655,7 @@ class dump: print >>f,snap.ylo,snap.yhi print >>f,snap.zlo,snap.zhi print >>f,"ITEM: ATOMS",namestr - + atoms = snap.atoms nvalues = len(atoms[0]) for i in xrange(snap.natoms): @@ -679,7 +679,7 @@ class dump: if not snap.tselect: continue print snap.time, sys.stdout.flush() - + file = root + "." + str(snap.time) f = open(file,"w") print >>f,"ITEM: TIMESTEP" @@ -691,7 +691,7 @@ class dump: print >>f,snap.ylo,snap.yhi print >>f,snap.zlo,snap.zhi print >>f,"ITEM: ATOMS",namestr - + atoms = snap.atoms nvalues = len(atoms[0]) for i in xrange(snap.natoms): @@ -733,7 +733,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] @@ -745,7 +745,7 @@ class dump: if not snap.tselect: continue for i in xrange(snap.natoms): if snap.aselect[i]: exec ceq - + # -------------------------------------------------------------------- # set a column value via an input vec for all selected snapshots/atoms @@ -765,7 +765,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 @@ -831,7 +831,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: @@ -847,13 +847,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 StandardError, "no columns specified" columns = [] @@ -908,7 +908,7 @@ class dump: del self.snaps[i] else: i += 1 - + # -------------------------------------------------------------------- # iterate over selected snapshots @@ -920,12 +920,12 @@ 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 # 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: @@ -949,7 +949,7 @@ class dump: # create atom list needed by viz from id,type,x,y,z # need Numeric/Numpy mode here - + atoms = [] for i in xrange(snap.natoms): if not snap.aselect[i]: continue @@ -989,7 +989,7 @@ class dump: 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 @@ -1000,7 +1000,7 @@ class dump: tmp1,tmp2,tmp3,tmp4,tmp5,lines = self.objextra.viz(time,1) return time,box,atoms,bonds,tris,lines - + # -------------------------------------------------------------------- def findtime(self,n): @@ -1045,7 +1045,7 @@ class dump: def extra(self,arg): # data object, grab bonds statically - + if type(arg) is types.InstanceType and ".data" in str(arg.__class__): self.bondflag = 0 try: @@ -1062,7 +1062,7 @@ class dump: raise StandardError,"could not extract bonds from data object" # cdata object, grab tris and lines statically - + elif type(arg) is types.InstanceType and ".cdata" in str(arg.__class__): self.triflag = self.lineflag = 0 try: @@ -1077,32 +1077,32 @@ class dump: raise StandardError,"could not extract tris/lines from cdata object" # mdump object, grab tris dynamically - + 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 tris 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 else: raise StandardError,"unrecognized argument to dump.extra()" - + # -------------------------------------------------------------------- def compare_atom(self,a,b): @@ -1111,7 +1111,7 @@ class dump: elif a[0] > b[0]: return 1 else: - return 0 + return 0 # -------------------------------------------------------------------- # one snapshot @@ -1126,15 +1126,15 @@ class tselect: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def all(self,**kwargs): - + # check whether to output or not outputfl = True if "output" in kwargs: outputfl = kwargs["output"] - + data = self.data for snap in data.snaps: snap.tselect = 1 @@ -1178,7 +1178,7 @@ class tselect: data.nselect -= 1 data.aselect.all() print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) - + # -------------------------------------------------------------------- def test(self,teststr): @@ -1224,7 +1224,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/src/dump2force.py b/src/dump2force.py index c45121c..44bcd08 100644 --- a/src/dump2force.py +++ b/src/dump2force.py @@ -4,10 +4,10 @@ A simple routine to load in a LIGGGHTS hybrid dump file containing contact and contact force data and convert into a .vtk unstructured grid which can be used to visualise the force network. This routine also writes the length of the connection between particles, in order -to be able to filter out incorrect connections (produced by the +to be able to filter out incorrect connections (produced by the "deform" fix) -This routine is based on Mark Bentley's dump2force (Space Research Institute, +This routine is based on Mark Bentley's dump2force (Space Research Institute, Austrian Academy of Sciences, mark.bentley@oeaw.ac.at) contributing author: Stefan Radl, TU Graz (radl@tugraz.at) @@ -39,7 +39,7 @@ import sys, os # Check for command line arguments if len(sys.argv) != 2: sys.exit('Usage: dump2force.py , where filename is a SINGLE filename; typically dump.') - + elif len(sys.argv) == 2: # we have one input param, that should be parsed as a filename filename = str(sys.argv[1]) if not os.path.isfile(filename): @@ -101,7 +101,7 @@ while timestep >= 0: # one datasets has some missing, data for the previous timestep are still displayed - # this means that it is better here to generate "empty" files for these timesteps. - if forcedata.snaps[fileindex].natoms == 0: + if forcedata.snaps[fileindex].natoms == 0: vtufile = fileprefix+'_'+str(timestep)+'.vtu' vtufile = os.path.join(outputdir,vtufile) vtuwrite = file(vtufile,'w') @@ -118,7 +118,7 @@ while timestep >= 0: """) - + else: # ****************************************** # Cell and connection lists @@ -168,7 +168,7 @@ while timestep >= 0: -np.array(forcedata.snaps[fileindex].atoms[:,forcedata.names["y2"]],dtype=np.float64))**2 \ + \ (np.array(forcedata.snaps[fileindex].atoms[:,forcedata.names["z1"]],dtype=np.float64) \ - -np.array(forcedata.snaps[fileindex].atoms[:,forcedata.names["z2"]],dtype=np.float64))**2 + -np.array(forcedata.snaps[fileindex].atoms[:,forcedata.names["z2"]],dtype=np.float64))**2 connectionLength = np.sqrt(connectionLength) @@ -217,7 +217,7 @@ while timestep >= 0: y = np.zeros( npoints, dtype=np.float64) z = np.zeros( npoints, dtype=np.float64) - counter = 0 + counter = 0 for id in ids: if id in id1: index = id1.index(id) @@ -229,10 +229,10 @@ while timestep >= 0: xtemp,ytemp,ztemp = forcedata.snaps[fileindex].atoms[index,forcedata.names["x2"]], \ forcedata.snaps[fileindex].atoms[index,forcedata.names["y2"]], \ forcedata.snaps[fileindex].atoms[index,forcedata.names["z2"]] - + x[counter]=xtemp y[counter]=ytemp - z[counter]=ztemp + z[counter]=ztemp counter += 1 # Now create the connectivity list - this corresponds to pairs of IDs, but referencing @@ -241,7 +241,7 @@ while timestep >= 0: # If the periodic flag is set for a given interactions, DO NOT connect the points # (to avoid lines that cross the simulation domain) - + # Mask out periodic interactions from the cell (connectivity) array # newList = [word for (word, mask) in zip(s,b) if mask] id1_masked = [ident for (ident,mask) in zip(id1,np.invert(periodic)) if mask] diff --git a/src/ensight.py b/src/ensight.py index 7c68ef1..d2cd2a1 100644 --- a/src/ensight.py +++ b/src/ensight.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. # ensight tool @@ -11,7 +11,7 @@ oneline = "Convert LAMMPS snapshots or meshes to Ensight format" docstr = """ -e = ensight(d) d = object with atoms or elements (dump,data,mdump) +e = ensight(d) d = object with atoms or elements (dump,data,mdump) e.change = 1 set to 1 if element nodal xyz change with time (def = 0) e.maxtype = 10 max particle type, set if query to data will be bad @@ -74,12 +74,12 @@ class ensight: self.which = 1 else: raise StandardError,"unrecognized object passed to ensight" - + # -------------------------------------------------------------------- def one(self,*args): if len(args) % 2 == 0: root = "tmp" - else: + else: root = args[0] args = args[1:] @@ -90,7 +90,7 @@ class ensight: if self.which == 0 and self.maxtype == 0: self.maxtype = self.data.maxtype() - + # write Ensight *.case header file f = open("%s.case" % root,"w") @@ -140,7 +140,7 @@ class ensight: else: self.variable_file_elements(vfiles[i],pairs[i][1],etype,values) print >>vfiles[i],"END TIME STEP" - + print time, sys.stdout.flush() n += 1 @@ -156,7 +156,7 @@ class ensight: def increment(self,*args): if len(args) % 2 == 0: root = "tmp" - else: + else: root = args[0] args = args[1:] @@ -214,7 +214,7 @@ class ensight: else: self.variable_file_elements(vfiles[i],pairs[i][1],etype,values) print >>vfiles[i],"END TIME STEP" - + print time, sys.stdout.flush() n += 1 @@ -236,7 +236,7 @@ class ensight: def many(self,*args): if len(args) % 2 == 0: root = "tmp" - else: + else: root = args[0] args = args[1:] @@ -268,20 +268,20 @@ class ensight: files = [] if n < 10: file = root + "000" + str(n) + ".xyz" - for pair in pairs: - files.append(root + "000" + str(n) + "." + pair[0]) + for pair in pairs: + files.append(root + "000" + str(n) + "." + pair[0]) elif n < 100: file = root + "00" + str(n) + ".xyz" - for pair in pairs: - files.append(root + "00" + str(n) + "." + pair[0]) + for pair in pairs: + files.append(root + "00" + str(n) + "." + pair[0]) elif n < 1000: file = root + "0" + str(n) + ".xyz" - for pair in pairs: - files.append(root + "0" + str(n) + "." + pair[0]) + for pair in pairs: + files.append(root + "0" + str(n) + "." + pair[0]) else: file = root + str(n) + ".xyz" - for pair in pairs: - files.append(root + str(n) + "." + pair[0]) + for pair in pairs: + files.append(root + str(n) + "." + pair[0]) if self.which == 0: f = open(file,"w") @@ -309,19 +309,19 @@ class ensight: self.variable_file_atoms(f,pairs[i][1],atoms,values) else: self.variable_file_elements(f,pairs[i][1],etype,values) - f.close() + f.close() print time, sys.stdout.flush() n += 1 - + print "\nwrote %s snapshots in Ensight format" % n # -------------------------------------------------------------------- def single(self,time,*args): if len(args) % 2 == 0: root = "tmp" - else: + else: root = args[0] args = args[1:] @@ -343,7 +343,7 @@ class ensight: which = self.data.findtime(time) etype = 0 - + f = open(root + ".xyz","w") if self.which == 0: time,box,atoms,bonds,tris,lines = self.data.viz(which) @@ -353,7 +353,7 @@ class ensight: self.coord_file_elements(f,box,nodes,elements) etype = len(elements[0]) f.close() - + for i in range(len(pairs)): values = self.data.vecs(time,pairs[i][0]) f = open(root + "." + pairs[i][0],"w") @@ -362,14 +362,14 @@ class ensight: else: self.variable_file_elements(f,pairs[i][1],etype,values) f.close() - + # -------------------------------------------------------------------- # write Ensight case file def case_file(self,f,root,pairs,multifile,nsnaps,times): print >>f,"# Ensight case file\n" print >>f,"FORMAT\ntype: ensight gold\n" - + if self.which == 0: if multifile: # print >>f,"GEOMETRY\nmodel: %s****.xyz change_coords_only\n" % root @@ -411,7 +411,7 @@ class ensight: if i % 10 == 9: print >>f print >>f print >>f - + if not multifile: print >>f,"FILE" print >>f,"file set: 1" diff --git a/src/gl.py b/src/gl.py index b0312cf..467ffce 100644 --- a/src/gl.py +++ b/src/gl.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. # gl tool @@ -16,8 +16,8 @@ g = gl(d) create OpenGL display for data in d d = atom snapshot object (dump, data) g.bg("black") set background color (def = "black") -g.size(N) set image size to NxN -g.size(N,M) set image size to NxM +g.size(N) set image size to NxN +g.size(N,M) set image size to NxM g.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) g.shift(x,y) translate by x,y pixels in view window (def = 0,0) g.zoom(0.5) scale image by factor (def = 1) @@ -27,7 +27,7 @@ g.box(0/1/2,"red",4) set box edge thickness g.file = "image" file prefix for created images (def = "image") g.show(N) show image of snapshot at timestep N - + g.all() make images of all selected snapshots g.all(P) images of all, start file label at P g.all(N,M,P) make M images of snapshot N, start label at P @@ -40,12 +40,12 @@ g.pan() no pan during all() (default) g.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() g.select = "" no extra aselect (default) - + %g varies from 0.0 to 1.0 from beginning to end of all() - -g.acol(2,"green") set atom colors by atom type (1-N) -g.acol([2,4],["red","blue"]) 1st arg = one type or list of types -g.acol(0,"blue") 2nd arg = one color or list of colors + +g.acol(2,"green") set atom colors by atom type (1-N) +g.acol([2,4],["red","blue"]) 1st arg = one type or list of types +g.acol(0,"blue") 2nd arg = one color or list of colors g.acol(range(20),["red","blue"]) if list lengths unequal, interpolate g.acol(range(10),"loop") assign colors in loop, randomly ordered @@ -55,23 +55,23 @@ g.acol(range(10),"loop") assign colors in loop, randomly ordered g.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() -g.bcol() set bond color, same args as acol() -g.brad() set bond thickness, same args as arad() +g.bcol() set bond color, same args as acol() +g.brad() set bond thickness, same args as arad() -g.tcol() set triangle color, same args as acol() -g.tfill() set triangle fill, 0 fill, 1 line, 2 both +g.tcol() set triangle color, same args as acol() +g.tfill() set triangle fill, 0 fill, 1 line, 2 both g.lcol() set line color, same args as acol() g.lrad() set line thickness, same args as arad() g.adef() set atom/bond/tri/line properties to default -g.bdef() default = "loop" for colors, 0.45 for radii -g.tdef() default = 0.25 for bond/line thickness -g.ldef() default = 0 fill +g.bdef() default = "loop" for colors, 0.45 for radii +g.tdef() default = 0.25 for bond/line thickness +g.ldef() default = 0 fill by default 100 types are assigned if atom/bond/tri/line has type > # defined properties, is an error - + from vizinfo import colors access color list print colors list defined color names and RGB values colors["nickname"] = [R,G,B] set new RGB values from 0 to 255 @@ -145,7 +145,7 @@ class gl: self.azphi = 30 self.scale = 1.0 self.xshift = self.yshift = 0 - + self.file = "image" self.boxflag = 0 self.bxcol = [1,1,0] @@ -162,7 +162,7 @@ class gl: self.nsides = 10 self.theta_amplify = 2 self.shiny = 2 - + self.clipflag = 0 self.clipxlo = self.clipylo = self.clipzlo = 0.0 self.clipxhi = self.clipyhi = self.clipzhi = 1.0 @@ -187,7 +187,7 @@ class gl: self.bdef() self.tdef() self.ldef() - + self.center = 3*[0] self.view = 3*[0] self.up = 3*[0] @@ -209,7 +209,7 @@ class gl: if not ynew: self.ypixels = self.xpixels else: self.ypixels = ynew self.create_window() - + # -------------------------------------------------------------------- def axis(self,value): @@ -221,7 +221,7 @@ class gl: def create_window(self): if self.root: self.root.destroy() - + from __main__ import tkroot self.root = Toplevel(tkroot) self.root.title('Pizza.py gl tool') @@ -230,7 +230,7 @@ class gl: double=1,depth=1) self.w.pack(expand=YES) # self.w.pack(expand=YES,fill=BOTH) - + glViewport(0,0,self.xpixels,self.ypixels) glEnable(GL_LIGHTING); glEnable(GL_LIGHT0); @@ -245,7 +245,7 @@ class gl: self.w.parent = self self.w.tkRedraw() tkroot.update_idletasks() # force window to appear - + # -------------------------------------------------------------------- def clip(self,which,value): @@ -312,7 +312,7 @@ class gl: self.up[1] = sin(pi*self.azphi/180) self.up[2] = 0.0 else: - dot = self.view[2] # dot = (0,0,1) . view + dot = self.view[2] # dot = (0,0,1) . view self.up[0] = -dot*self.view[0] # up projected onto v = dot * v self.up[1] = -dot*self.view[1] # up perp to v = up - dot * v self.up[2] = 1.0 - dot*self.view[2] @@ -323,7 +323,7 @@ class gl: # -------------------------------------------------------------------- # reset ztheta,azphi and thus view,up.right # called as function from Pizza.py - + def rotate(self,ztheta,azphi): self.ztheta = ztheta self.azphi = azphi @@ -364,11 +364,11 @@ class gl: # rotate view,up around axis of rotation = old x new # right = up x view # reset ztheta,azphi from view - + def mouse_rotate(self,xnew,ynew,xold,yold): # change y pixels to measure from bottom of window instead of top - + yold = self.ypixels - yold ynew = self.ypixels - ynew @@ -405,7 +405,7 @@ class gl: axis[1] = rot[0]*self.right[1] + rot[1]*self.up[1] + rot[2]*self.view[1] axis[2] = rot[0]*self.right[2] + rot[1]*self.up[2] + rot[2]*self.view[2] axis = vecnorm(axis) - + # view is changed by (axis x view) scaled by theta # up is changed by (axis x up) scaled by theta # force up to be perp to view via up_perp = up - (up . view) view @@ -466,14 +466,14 @@ class gl: # output: eye = distance to view scene from # xto,yto,zto = point to look to # xfrom,yfrom,zfrom = point to look from - + def setview(self): if not self.ready: return # no distance since no scene yet - + self.eye = 3 * self.distance / self.scale xfactor = 0.5*self.eye*self.xshift/self.xpixels yfactor = 0.5*self.eye*self.yshift/self.ypixels - + self.xto = self.center[0] - xfactor*self.right[0] - yfactor*self.up[0] self.yto = self.center[1] - xfactor*self.right[1] - yfactor*self.up[1] self.zto = self.center[2] - xfactor*self.right[2] - yfactor*self.up[2] @@ -484,7 +484,7 @@ class gl: # -------------------------------------------------------------------- # box attributes, also used for triangle lines - + def box(self,*args): self.boxflag = args[0] if len(args) > 1: @@ -498,7 +498,7 @@ class gl: # -------------------------------------------------------------------- # grab all selected snapshots from data object # add GL-specific info to each bond - + def reload(self): print "Loading data into gl tool ..." data = self.data @@ -527,7 +527,7 @@ class gl: self.bondframes.append(bonds) self.triframes.append(tris) self.lineframes.append(lines) - + print time, sys.stdout.flush() print @@ -543,11 +543,11 @@ class gl: def nolabel(self): self.cachelist = -self.cachelist self.labels = [] - + # -------------------------------------------------------------------- # show a single snapshot # distance from snapshot box or max box for all selected steps - + def show(self,ntime): data = self.data which = data.findtime(ntime) @@ -569,7 +569,7 @@ class gl: self.cachelist = -self.cachelist self.w.tkRedraw() self.save() - + # -------------------------------------------------------------------- def pan(self,*list): @@ -582,7 +582,7 @@ class gl: self.ztheta_stop = list[3] self.azphi_stop = list[4] self.scale_stop = list[5] - + # -------------------------------------------------------------------- def all(self,*list): @@ -613,7 +613,7 @@ class gl: if flag == -1: break fraction = float(i) / (ncount-1) - + if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,time) @@ -636,7 +636,7 @@ class gl: fraction*(self.scale_stop - self.scale_start) self.viewupright() - if n == nstart or self.panflag: self.center = compute_center(box) + if n == nstart or self.panflag: self.center = compute_center(box) if bonds: self.bonds_augment(bonds) @@ -651,7 +651,7 @@ class gl: self.cachelist = -self.cachelist self.w.tkRedraw() self.save(file) - + print time, sys.stdout.flush() i += 1 @@ -691,7 +691,7 @@ class gl: fraction*(self.scale_stop - self.scale_start) self.viewupright() - if n == nstart or self.panflag: self.center = compute_center(box) + if n == nstart or self.panflag: self.center = compute_center(box) if bonds: self.bonds_augment(bonds) @@ -729,19 +729,19 @@ class gl: # -------------------------------------------------------------------- # draw the GL scene - + def redraw(self,o): # clear window to background color - + glClearColor(self.bgcol[0],self.bgcol[1],self.bgcol[2],0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) # not ready if no scene yet - + if not self.ready: return # set view from eye, distance, 3 lookat vectors (from,to,up) - + glMatrixMode(GL_PROJECTION) glLoadIdentity() if self.orthoflag: @@ -757,14 +757,14 @@ class gl: # draw scene from display list if caching allowed and list hasn't changed # else redraw and store as new display list if caching allowed - + if self.cache and self.cachelist > 0: glCallList(self.cachelist); else: if self.cache: if self.cachelist < 0: glDeleteLists(-self.cachelist,1) self.cachelist = glGenLists(1) glNewList(self.cachelist,GL_COMPILE_AND_EXECUTE) - + # draw box, clip-box, xyz axes, lines glDisable(GL_LIGHTING) @@ -781,7 +781,7 @@ class gl: red,green,blue = self.vizinfo.lcolor[itype] glColor3f(red,green,blue) thick = self.vizinfo.lrad[itype] - glLineWidth(thick) + glLineWidth(thick) glBegin(GL_LINES) glVertex3f(line[2],line[3],line[4]) glVertex3f(line[5],line[6],line[7]) @@ -840,7 +840,7 @@ class gl: if self.tridraw: fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] - + if fillflag != 1: if fillflag: glEnable(GL_POLYGON_OFFSET_FILL) @@ -919,7 +919,7 @@ class gl: gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) glPopMatrix() - if self.tridraw: + if self.tridraw: fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] if fillflag != 1: @@ -973,7 +973,7 @@ class gl: glEnd() glEnable(GL_LIGHTING) glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) - + if self.cache: glEndList() glFlush() @@ -981,16 +981,16 @@ class gl: # -------------------------------------------------------------------- # make new call list for each atom type # called when atom color/rad/quality is changed - + def make_atom_calllist(self): # extend calllist array if necessary - + if self.vizinfo.nacolor > self.nclist: for i in range(self.vizinfo.nacolor-self.nclist): self.calllist.append(0) self.nclist = self.vizinfo.nacolor # create new calllist for each atom type - + for itype in xrange(1,self.vizinfo.nacolor+1): if self.calllist[itype]: glDeleteLists(self.calllist[itype],1) ilist = glGenLists(1) @@ -999,12 +999,12 @@ class gl: red,green,blue = self.vizinfo.acolor[itype] rad = self.vizinfo.arad[itype] glColor3f(red,green,blue); - + # glPointSize(10.0*rad) # glBegin(GL_POINTS) # glVertex3f(0.0,0.0,0.0) # glEnd() - + glMaterialfv(GL_FRONT,GL_EMISSION,[red,green,blue,1.0]); glMaterialf(GL_FRONT,GL_SHININESS,self.shiny); glutSolidSphere(rad,self.nslices,self.nstacks) @@ -1013,7 +1013,7 @@ class gl: # -------------------------------------------------------------------- # augment bond info returned by viz() with info needed for GL draw # info = length, theta, -dy, dx for bond orientation - + def bonds_augment(self,bonds): for bond in bonds: dx = bond[5] - bond[2] @@ -1044,7 +1044,7 @@ class gl: glLineWidth(self.bxthick) glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) - + glBegin(GL_LINE_LOOP) glVertex3f(xlo,ylo,zlo) glVertex3f(xhi,ylo,zlo) @@ -1079,7 +1079,7 @@ class gl: if yhi-ylo > delta: delta = yhi-ylo if zhi-zlo > delta: delta = zhi-zlo delta *= 0.1 - + glLineWidth(self.bxthick) glBegin(GL_LINES) @@ -1098,7 +1098,7 @@ class gl: def save(self,file=None): self.w.update() # force image on screen to be current before saving it - + pstring = glReadPixels(0,0,self.xpixels,self.ypixels, GL_RGBA,GL_UNSIGNED_BYTE) snapshot = Image.fromstring("RGBA",(self.xpixels,self.ypixels),pstring) @@ -1108,14 +1108,14 @@ class gl: snapshot.save(file + ".png") # -------------------------------------------------------------------- - + def adef(self): self.vizinfo.setcolors("atom",range(100),"loop") self.vizinfo.setradii("atom",range(100),0.45) self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() - + # -------------------------------------------------------------------- def bdef(self): @@ -1128,14 +1128,14 @@ class gl: def tdef(self): self.vizinfo.setcolors("tri",range(100),"loop") - self.vizinfo.setfills("tri",range(100),0) + self.vizinfo.setfills("tri",range(100),0) self.cachelist = -self.cachelist self.w.tkRedraw() # -------------------------------------------------------------------- def ldef(self): - self.vizinfo.setcolors("line",range(100),"loop") + self.vizinfo.setcolors("line",range(100),"loop") self.vizinfo.setradii("line",range(100),0.25) self.cachelist = -self.cachelist self.w.tkRedraw() @@ -1147,29 +1147,29 @@ class gl: self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() - + # -------------------------------------------------------------------- def arad(self,atypes,radii): - self.vizinfo.setradii("atom",atypes,radii) + self.vizinfo.setradii("atom",atypes,radii) self.make_atom_calllist() self.cachelist = -self.cachelist self.w.tkRedraw() - + # -------------------------------------------------------------------- def bcol(self,btypes,colors): self.vizinfo.setcolors("bond",btypes,colors) self.cachelist = -self.cachelist self.w.tkRedraw() - + # -------------------------------------------------------------------- def brad(self,btypes,radii): self.vizinfo.setradii("bond",btypes,radii) self.cachelist = -self.cachelist self.w.tkRedraw() - + # -------------------------------------------------------------------- def tcol(self,ttypes,colors): @@ -1208,10 +1208,10 @@ class MyOpengl(Opengl): args = (self,master,cnf) Opengl.__init__(*args,**kw) Opengl.autospin_allowed = 0 - + # redraw Opengl scene # call parent redraw() method - + def tkRedraw(self,*dummy): if not self.initialised: return self.tk.call(self._w,'makecurrent') @@ -1220,7 +1220,7 @@ class MyOpengl(Opengl): # left button translate # access parent xshift/yshift and call parent trans() method - + def tkTranslate(self,event): dx = event.x - self.xmouse dy = event.y - self.ymouse @@ -1240,7 +1240,7 @@ class MyOpengl(Opengl): # right button zoom # access parent scale and call parent zoom() method - + def tkScale(self,event): scale = 1 - 0.01 * (event.y - self.ymouse) if scale < 0.001: scale = 0.001 diff --git a/src/gnu.py b/src/gnu.py index f6f0167..d99ab38 100644 --- a/src/gnu.py +++ b/src/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. # gnu tool @@ -11,12 +11,12 @@ 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 +29,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: @@ -94,7 +94,7 @@ except: PIZZA_GNUTERM = "x11" # Class definition class gnu: - + # -------------------------------------------------------------------- def __init__(self): @@ -102,7 +102,7 @@ class gnu: self.file = "tmp.gnu" self.figures = [] self.select(1) - + # -------------------------------------------------------------------- def stop(self): @@ -114,7 +114,7 @@ class gnu: def __call__(self,command): self.GNUPLOT.write(command + '\n') self.GNUPLOT.flush() - + # -------------------------------------------------------------------- def enter(self): @@ -152,7 +152,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) @@ -160,7 +160,7 @@ class gnu: self.save(newfile) n += 1 - + # -------------------------------------------------------------------- # write list of equal-length vectors to filename @@ -201,7 +201,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 @@ -212,7 +212,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 @@ -221,7 +221,7 @@ class gnu: fig.ncurves = self.figures[self.current-1].ncurves self.figures[self.current-1] = fig self.draw() - + # -------------------------------------------------------------------- def aspect(self,value): @@ -245,12 +245,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() # -------------------------------------------------------------------- @@ -259,7 +259,7 @@ class gnu: self.figures[self.current-1].nlabel = 0 self.figures[self.current-1].labels = [] self.draw() - + # -------------------------------------------------------------------- def title(self,*strings): @@ -276,13 +276,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): @@ -291,7 +291,7 @@ class gnu: else: self.figures[self.current-1].xlog = 1 self.draw() - + # -------------------------------------------------------------------- def ylog(self): @@ -300,7 +300,7 @@ class gnu: else: self.figures[self.current-1].ylog = 1 self.draw() - + # -------------------------------------------------------------------- def curve(self,num,color): @@ -316,10 +316,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 + '"' @@ -331,11 +331,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 [*:*]") @@ -365,7 +365,7 @@ class figure: def __init__(self): self.ncurves = 0 - self.colors = [] + self.colors = [] self.title = "" self.xtitle = "" self.ytitle = "" diff --git a/src/histo.py b/src/histo.py index f1bd621..f9ce04c 100644 --- a/src/histo.py +++ b/src/histo.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. # histo tool @@ -36,7 +36,7 @@ class histo: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def compute(self,dim,nbins,lo=None,hi=None): @@ -46,7 +46,7 @@ class histo: else: raise StandardError,"illegal dim value" y = nbins*[0] - + count = 0 n = flag = 0 while 1: @@ -67,19 +67,19 @@ class histo: delta = (hi-lo) / nbins; invdelta = 1.0/delta - + for atom in atoms: coord = atom[idim] ibin = int((coord-lo) * invdelta) if ibin < 0 or ibin >= nbins: continue y[ibin] += 1 count += 1 - + n += 1 x = nbins*[0] for i in xrange(nbins): x[i] = (i+0.5)*delta - + print "histogram snapshots = ",n print "histogram counts (per snap) = %d (%g)" % (count,float(count)/n) print "histogram bounds = ",lo,hi diff --git a/src/image.py b/src/image.py index f8611b5..3aea788 100644 --- a/src/image.py +++ b/src/image.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. # image tool @@ -13,14 +13,14 @@ oneline = "View and manipulate images" docstr = """ i = image("my1.gif my2.gif") display thumbnails of matching images i = image("*.png *.gif") wildcards allowed -i = image("") blank string matches all image suffixes -i = image() no display window opened if no arg +i = image("") blank string matches all image suffixes +i = image() no display window opened if no arg image suffixes for blank string = *.png, *.bmp, *.gif, *.tiff, *.tif click on a thumbnail to view it full-size click on thumbnail again to remove full-sized version -i.view("*.png *.gif") display thumbnails of matching images +i.view("*.png *.gif") display thumbnails of matching images view arg is same as constructor arg @@ -62,7 +62,7 @@ except: PIZZA_MONTAGE = "montage" # Class definition class image: - + # -------------------------------------------------------------------- def __init__(self,filestr=None): @@ -84,24 +84,24 @@ class image: from __main__ import tkroot # GUI control window - + gui = Toplevel(tkroot) gui.title('Pizza.py image tool') - + scroll = \ Pmw.ScrolledFrame(gui,usehullsize=1,hull_width=420,hull_height=500) pane = scroll.interior() - + ncolumns = 4 for i in xrange(len(files)): - + # create new row frame if 1st in column - + if i % ncolumns == 0: rowframe = Frame(pane) oneframe = Frame(rowframe) - + # create a thumbnail of image - + im = Image.open(files[i]) imt = im.copy() imt.thumbnail((60,60),Image.ANTIALIAS) @@ -109,7 +109,7 @@ class image: imt.save("tmp." + basename) thumbnail = ImageTk.PhotoImage(file = "tmp." + basename) os.remove("tmp." + basename) - + # read in full size image # create a thumbnail object that links to it # create button that calls the thumbnail, label with filename @@ -119,18 +119,18 @@ class image: obj = thumbnails(gui,files[i],big,thumbnail) Button(oneframe,image=thumbnail,command=obj.display).pack(side=TOP) Label(oneframe,text=basename).pack(side=BOTTOM) - + # pack into row frame - + oneframe.pack(side=LEFT) if (i+1) % ncolumns == 0: rowframe.pack(side=TOP) - + if len(files) % ncolumns != 0: rowframe.pack(side=TOP) - scroll.pack(side=LEFT) + scroll.pack(side=LEFT) # -------------------------------------------------------------------- # wrapper on ImageMagick convert command - + def convert(self,file1,file2,switch=""): if file1.find('*') < 0 or file2.find('*') < 0: cmd = "%s %s %s %s" % (PIZZA_CONVERT,switch,file1,file2) @@ -168,7 +168,7 @@ class image: for j in range(nsets): cmd += " %s" % fileargs[j] commands.getoutput(cmd) return - + nfiles = len(glob.glob(fileargs[0])) filesets = [] for i in range(nsets-1): @@ -197,7 +197,7 @@ class image: # -------------------------------------------------------------------- # thumbnail class - + class thumbnails: def __init__(self,root,name,bigimage,thumbimage): @@ -207,19 +207,19 @@ class thumbnails: self.name = name self.bigexist = 0 self.window = None - + def display(self): # destroy the big image window - + if self.bigexist: self.bigexist = 0 if self.window: self.window.destroy() - self.window = None - + self.window = None + # create a new window with the big image - + else: self.bigexist = 1 self.window = Toplevel(self.root) @@ -230,4 +230,4 @@ class thumbnails: # list of file extensions to test for # could add any extensions that PIL recognizes -extensions = ["*.png", "*.bmp", "*.gif", "*.tiff", "*.tif"] +extensions = ["*.png", "*.bmp", "*.gif", "*.tiff", "*.tif"] diff --git a/src/ldump.py b/src/ldump.py index 30bf1bf..4f07456 100644 --- a/src/ldump.py +++ b/src/ldump.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. # ldump tool @@ -12,14 +12,14 @@ oneline = "Read dump files with line segment info" docstr = """ l = ldump("dump.one") read in one or more dump files -l = ldump("dump.1 dump.2.gz") can be gzipped -l = ldump("dump.*") wildcard expands to multiple files -l = ldump("dump.*",0) two args = store filenames, but don't read +l = ldump("dump.1 dump.2.gz") can be gzipped +l = ldump("dump.*") wildcard expands to multiple files +l = ldump("dump.*",0) two args = store filenames, but don't read incomplete and duplicate snapshots are deleted no column name assignment is performed -time = l.next() read next snapshot from dump files +time = l.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 @@ -43,7 +43,7 @@ time,box,atoms,bonds,tris,lines = l.viz(index) return list of viz objects lines = id,type,x1,y1,z1,x2,y2,z2 for each line as 2d array id,type are from associated atom -l.owrap(...) wrap lines to same image as their atoms +l.owrap(...) wrap lines to same image as their atoms owrap() is called by dump tool's owrap() useful for wrapping all molecule's atoms/lines the same so it is contiguous @@ -100,7 +100,7 @@ class ldump: for word in words: self.flist += glob.glob(word) if len(self.flist) == 0 and len(list) == 1: raise StandardError,"no ldump file specified" - + if len(list) == 1: self.increment = 0 self.read_all() @@ -155,15 +155,15 @@ class ldump: snap = self.read_snapshot(f) if not snap: self.nextfile += 1 - if self.nextfile == len(self.flist): return -1 + if self.nextfile == len(self.flist): return -1 f.close() - self.eof = 0 - continue + self.eof = 0 + continue self.eof = f.tell() f.close() try: self.findtime(snap.time) - continue + continue except: break self.snaps.append(snap) @@ -175,7 +175,7 @@ class ldump: # -------------------------------------------------------------------- # read a single snapshot from file f # return snapshot or 0 if failed - + def read_snapshot(self,f): try: snap = Snap() @@ -216,7 +216,7 @@ class ldump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise StandardError, "ldump map() requires pairs of mappings" @@ -263,7 +263,7 @@ class ldump: del self.snaps[i] else: i += 1 - + # -------------------------------------------------------------------- # return list of lines to viz for snapshot isnap # if called with flag, then index is timestep, so convert to snapshot index @@ -291,7 +291,7 @@ class ldump: # create line list from id,type,end1x,end1y,end2x,end2y # don't add line if all 4 values are 0 since not a line - + lines = [] for i in xrange(snap.natoms): atom = snap.atoms[i] @@ -322,7 +322,7 @@ class ldump: # idump = index of my line I in dump's atoms # jdump = atom J in dump's atoms that atom I was owrapped on # delx,dely = offset applied to atom I and thus to line I - + for i in xrange(snap.natoms): tag = atoms[i][id] idump = idsdump[tag] diff --git a/src/log.py b/src/log.py index 7020c21..73c35ea 100644 --- a/src/log.py +++ b/src/log.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. # log tool @@ -28,7 +28,7 @@ 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") 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 @@ -89,7 +89,7 @@ class log: # -------------------------------------------------------------------- # 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" @@ -100,7 +100,7 @@ class log: print # sort entries by timestep, cull duplicates - + self.data.sort(self.compare) self.cull() self.nlen = len(self.data) @@ -133,9 +133,9 @@ class log: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: @@ -161,9 +161,9 @@ class log: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: @@ -224,7 +224,7 @@ class log: keywords.insert(0,"Step") i = 0 for keyword in keywords: - self.names.append(keyword) + self.names.append(keyword) self.ptr[keyword] = i i += 1 @@ -234,7 +234,7 @@ class log: line = txt[s1:s2] words = line.split() for i in range(len(words)): - self.names.append(words[i]) + self.names.append(words[i]) self.ptr[words[i]] = i self.nvec = len(self.names) @@ -273,43 +273,43 @@ class log: 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 + 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 + 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 + s1 = txt.find("\n",s1) + 1 s2 = txt.rfind("\n",s1) + 1 - eof -= len(txt) - s2 + eof -= len(txt) - s2 elif s1 == -1 and s2 == -1: # found neither # could be end-of-file section - # or entire read was one chunk + # 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 + eof -= len(txt) - start # reset eof to "Loop" + break - last = 1 # entire read is a chunk + 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 + 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") diff --git a/src/lpp.py b/src/lpp.py index 1f1001d..9840022 100755 --- a/src/lpp.py +++ b/src/lpp.py @@ -18,7 +18,7 @@ import exceptions import getopt class lpp: - + #============================================================================= # creates a filelist, seperates it to sublists # creates multiple processes @@ -27,7 +27,7 @@ class lpp: # calls dump, vtk and manyGran for the given list of files # returns 0 #============================================================================= - + def __init__(self, *list, **kwargs): # do argument parsing, raise errors if non-integers were given # this can be changed if one wants less overhead but use more memory: @@ -36,7 +36,7 @@ class lpp: self.chunksize = 8 self.overwrite = True - if "--chunksize" in kwargs: + if "--chunksize" in kwargs: try: if int(kwargs["--chunksize"]) > 0: self.chunksize = int(kwargs["--chunksize"]) @@ -44,7 +44,7 @@ class lpp: except ValueError: raise ValueError, "Invalid or no argument given for chunksize" - if "--cpunum" in kwargs: + if "--cpunum" in kwargs: try: if int(kwargs["--cpunum"]) > 0 and int(kwargs["--cpunum"]) <= self.cpunum: self.cpunum = int(kwargs["--cpunum"]) @@ -55,11 +55,11 @@ class lpp: # do not overwrite existing files if "--no-overwrite" in kwargs: self.overwrite = False - + # suppress output with 'False' if "--debug" in kwargs: self.debugMode = True else: self.debugMode = False - + if "--quiet" in kwargs: self.output = False self.debugMode = False @@ -72,7 +72,7 @@ class lpp: starttime = time.time() if self.debugMode: print "number of process:", os.getpid() - + # check whether file-list is nonempty self.flist = [] # get file list for windows @@ -89,15 +89,15 @@ class lpp: raise StandardError, "no dump file specified" if listlen == 1 and self.overwrite == False: raise StandardError, "Cannot process single dump files with --no-overwrite." - + if self.output: print "Working with", self.cpunum, "processes..." - + # seperate list in pieces+rest self.slices = [] - + residualPresent = int(bool(listlen-floor(listlen/self.chunksize)*self.chunksize)) - + for i in xrange(int(floor(listlen/self.chunksize))+residualPresent): slice = self.flist[i*self.chunksize:(i+1)*self.chunksize] self.slices.append(slice) @@ -111,31 +111,31 @@ class lpp: "output":output,\ "overwrite":self.overwrite} \ for i in xrange(len(self.slices))] - + if self.debugMode: print "dumpInput:",dumpInput - + numberOfRuns = len(dumpInput) i = 0 while i < len(dumpInput): - if self.output: + if self.output: print "calculating chunks",i+1,"-",min(i+self.cpunum,numberOfRuns),"of",numberOfRuns - + if self.debugMode: print "input of this \"map\": ",dumpInput[i:i+self.cpunum] - + # create job_server job_server = multiprocessing.Pool(processes = self.cpunum) - + # map lppWorker on all inputs via job_server (parallelly) job_server.map_async(lppWorker,dumpInput[i:i+self.cpunum]).get(9999999) - + # close jobserver job_server.close() job_server.join() i += self.cpunum - + endtime = time.time() - if self.output: + if self.output: print "wrote", listlen,"granular snapshots in VTK format" print "time needed:",endtime-starttime,"sec" @@ -144,7 +144,7 @@ def lppWorker(input): debugMode = input["debugMode"] outfileName = input["output"] overwrite = input["overwrite"] - + flistlen = len(flist) # generate name of manyGran splitfname = flist[0].rsplit(".") @@ -154,7 +154,7 @@ def lppWorker(input): granName = outfileName + splitfname[len(splitfname)-1] else: granName = outfileName - + # if no-overwrite: read timestamp in first line of file # if corresponding dump-file does not already exists: add it to 'shortFlist' # shortFlist ... list of files to finally be processed by dump, and vtk. @@ -174,13 +174,13 @@ def lppWorker(input): ff.close() except: continue - + # generate filename from time like in vtk, # check if file exists; if yes: do not add to list filename,file_bb,file_walls = vtk.generateFilename(granName,[time],0) if not os.path.isfile(filename): shortFlist.append(f) - + # call dump, vtk, manyGran on shortFlist try: d = dump({"filelist":shortFlist, "debugMode":debugMode}) @@ -189,7 +189,7 @@ def lppWorker(input): v.manyGran(granName,fileNos=d.fileNums,output=debugMode) except KeyboardInterrupt: raise - + return 0 def printHelp(): @@ -223,10 +223,10 @@ if __name__ == "__main__": # except: # if sys.exc_type == exceptions.SystemExit: # pass - # else: + # else: # print sys.exc_info() #=========================================================================== else: printHelp() - + diff --git a/src/matlab.py b/src/matlab.py index 9cde9b0..2c55ea2 100644 --- a/src/matlab.py +++ b/src/matlab.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. # matlab tool @@ -11,12 +11,12 @@ oneline = "Create plots via MatLab numerical analysis program" docstr = """ -m = matlab() start up MatLab -m.stop() shut down MatLab process - +m = matlab() start up MatLab +m.stop() shut down MatLab process + m.plot(a) plot vector A against linear index -m.plot(a,b) plot B against A -m.plot(a,b,c,d,...) plot B against A, D against C, etc +m.plot(a,b) plot B against A +m.plot(a,b,c,d,...) plot B against A, D against C, etc m.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,10 +29,10 @@ m.mplot(M,N,S,"file",a,b,...) multiple plots saved to file0000.eps, etc m("c = a + b") execute string in MatLab -m.enter() enter MatLab shell +m.enter() enter MatLab shell matlab> c = a + b type commands directly to MatLab -matlab> exit, quit exit MatLab shell - +matlab> exit, quit exit MatLab shell + m.export("data",range(100),a,...) create file with columns of numbers all vectors must be of equal length @@ -40,12 +40,12 @@ m.export("data",range(100),a,...) create file with columns of numbers cols = importdata('data') plot(cols(:,1),cols(:,2)) -m.select(N) figure N becomes the current plot - +m.select(N) figure N becomes the current plot + subsequent commands apply to this plot -m.hide(N) delete window for figure N -m.save("file") save current plot as file.eps +m.hide(N) delete window for figure N +m.save("file") save current plot as file.eps Set attributes for current plot: @@ -101,7 +101,7 @@ except: PIZZA_MATLAB = "matlab -nosplash -nodesktop -nojvm" # Class definition class matlab: - + # -------------------------------------------------------------------- def __init__(self): @@ -109,7 +109,7 @@ class matlab: self.file = "tmp.matlab" self.figures = [] self.select(1) - + # -------------------------------------------------------------------- def stop(self): @@ -121,7 +121,7 @@ class matlab: def __call__(self,command): self.MATLAB.write(command + '\n') self.MATLAB.flush() - + # -------------------------------------------------------------------- def enter(self): @@ -146,11 +146,11 @@ class matlab: self.export(file,vectors[i],vectors[i+1]) self.figures[self.current-1].ncurves = len(vectors)/2 self.draw() - + # -------------------------------------------------------------------- # create multiple plots from growing vectors, save to numbered files # don't plot empty vector, create a [0] instead - + def mplot(self,start,stop,skip,file,*vectors): n = 0 for i in range(start,stop,skip): @@ -212,7 +212,7 @@ class matlab: self.__call__(cmd) self.__call__("!touch tmp.done") while not os.path.exists("tmp.done"): continue - + # -------------------------------------------------------------------- # restore default attributes by creating a new fig object @@ -221,7 +221,7 @@ class matlab: fig.ncurves = self.figures[self.current-1].ncurves self.figures[self.current-1] = fig self.draw() - + # -------------------------------------------------------------------- def aspect(self,value): @@ -245,12 +245,12 @@ class matlab: 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() # -------------------------------------------------------------------- @@ -259,7 +259,7 @@ class matlab: self.figures[self.current-1].nlabel = 0 self.figures[self.current-1].labels = [] self.draw() - + # -------------------------------------------------------------------- def title(self,*strings): @@ -276,13 +276,13 @@ class matlab: 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): @@ -291,7 +291,7 @@ class matlab: else: self.figures[self.current-1].xlog = 1 self.draw() - + # -------------------------------------------------------------------- def ylog(self): @@ -300,7 +300,7 @@ class matlab: else: self.figures[self.current-1].ylog = 1 self.draw() - + # -------------------------------------------------------------------- def curve(self,num,*settings): @@ -348,16 +348,16 @@ class matlab: self.__call__("ylabel('%s','FontSize',16)" % fig.ytitle) if fig.xlimit == 0 or fig.ylimit == 0: self.__call__("axis auto") - if fig.xlimit: + if fig.xlimit: self.__call__("xlim([%g,%g])" % (fig.xlimit[0],fig.xlimit[1])) - if fig.ylimit: + if fig.ylimit: self.__call__("ylim([%g,%g])" % (fig.ylimit[0],fig.ylimit[1])) for i in range(fig.nlabels): x = fig.labels[i][0] y = fig.labels[i][1] text = fig.labels[i][2] # kludge to set label font size - self.__call__("text(%g,%g,'%s','FontSize',16)" % (x,y,text)) + self.__call__("text(%g,%g,'%s','FontSize',16)" % (x,y,text)) # -------------------------------------------------------------------- # class to store settings for a single plot diff --git a/src/mdump.py b/src/mdump.py index 82c5dc6..089b4d3 100644 --- a/src/mdump.py +++ b/src/mdump.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. # mdump tool @@ -12,13 +12,13 @@ oneline = "Read, write, manipulate mesh dump files" docstr = """ m = mdump("mesh.one") read in one or more mesh dump files -m = mdump("mesh.1 mesh.2.gz") can be gzipped -m = mdump("mesh.*") wildcard expands to multiple files -m = mdump("mesh.*",0) two args = store filenames, but don't read +m = mdump("mesh.1 mesh.2.gz") can be gzipped +m = mdump("mesh.*") wildcard expands to multiple files +m = mdump("mesh.*",0) two args = store filenames, but don't read incomplete and duplicate snapshots are deleted -time = m.next() read next snapshot from dump files +time = m.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,20 +28,20 @@ time = m.next() read next snapshot from dump files m.map(2,"temperature") assign names to element value columns (1-N) -m.tselect.all() select all timesteps -m.tselect.one(N) select only timestep N -m.tselect.none() deselect all timesteps -m.tselect.skip(M) select every Mth step +m.tselect.all() select all timesteps +m.tselect.one(N) select only timestep N +m.tselect.none() deselect all timesteps +m.tselect.skip(M) select every Mth step m.tselect.test("$t >= 100 and $t < 10000") select matching timesteps -m.delete() delete non-selected timesteps +m.delete() delete non-selected timesteps selecting a timestep also selects all elements 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 -m.eselect.all() select all elems in all steps -m.eselect.all(N) select all elems in one step +m.eselect.all() select all elems in all steps +m.eselect.all(N) select all elems in one step m.eselect.test("$id > 100 and $type == 2") select match elems in all steps m.eselect.test("$id > 100 and $type == 2",N) select matching elems in one step @@ -52,7 +52,7 @@ m.eselect.test("$id > 100 and $type == 2",N) select matching elems in one step Python comparison syntax: == != < > <= >= and or $name must end with a space -t = m.time() return vector of selected timestep values +t = m.time() return vector of selected timestep values fx,fy,... = m.vecs(1000,"fx","fy",...) return vector(s) for timestep N vecs() returns vectors with one value for each selected elem in the timestep @@ -161,7 +161,7 @@ class mdump: 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() @@ -228,7 +228,7 @@ class mdump: # reference definitions of nodes and elements in previous timesteps self.reference() - + self.nsnaps = len(self.snaps) print "read %d snapshots" % self.nsnaps @@ -253,15 +253,15 @@ class mdump: snap = self.read_snapshot(f) if not snap: self.nextfile += 1 - if self.nextfile == len(self.flist): return -1 + if self.nextfile == len(self.flist): return -1 f.close() - self.eof = 0 - continue + self.eof = 0 + continue self.eof = f.tell() f.close() try: self.findtime(snap.time) - continue + continue except: break # select the new snapshot with all its elements @@ -308,7 +308,7 @@ class mdump: 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 n: words = f.readline().split() @@ -340,7 +340,7 @@ class mdump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise StandardError, "mdump map() requires pairs of mappings" @@ -407,7 +407,7 @@ class mdump: snap = self.snaps[self.findtime(n)] if not snap.evalues: raise StandardError, "snapshot has no element values" - + if len(list) == 0: raise StandardError, "no columns specified" columns = [] @@ -441,7 +441,7 @@ class mdump: # -------------------------------------------------------------------- # delete successive snapshots with duplicate time stamp # if have same timestamp, combine them if internal flags are different - + def cull(self): i = 1 while i < len(self.snaps): @@ -480,7 +480,7 @@ class mdump: # -------------------------------------------------------------------- # insure every snapshot has node and element connectivity info # if not, point it at most recent shapshot that does - + def reference(self): for i in xrange(len(self.snaps)): if not self.snaps[i].nflag: @@ -520,11 +520,11 @@ class mdump: self.iterate = i return i,self.snaps[i].time,1 return 0,0,-1 - + # -------------------------------------------------------------------- # return list of triangles to viz for snapshot isnap # if called with flag, then index is timestep, so convert to snapshot index - + def viz(self,index,flag=0): if not flag: isnap = index else: @@ -536,7 +536,7 @@ class mdump: i += 1 isnap = i - 1 snap = self.snaps[isnap] - + time = snap.time box = [snap.xlo,snap.ylo,snap.zlo,snap.xhi,snap.yhi,snap.zhi] if self.etype == "": type = -1 @@ -547,7 +547,7 @@ class mdump: # create triangle list from all elements # for type, either use element type (-1) or user-defined column in evalues - + tris = [] nodes = snap.nodes for i in xrange(snap.nelements): @@ -557,7 +557,7 @@ class mdump: else: evalues = [] # single tri, normal = up - + if snap.eflag == 1: v1 = nodes[int(element[2])-1][2:5].tolist() v2 = nodes[int(element[3])-1][2:5].tolist() @@ -568,7 +568,7 @@ class mdump: else: tris.append([element[0],evalue[type]] + list + n) # single tet, convert to 4 tris, normals = out - + elif snap.eflag == 2: v1 = nodes[int(element[2])-1][2:5].tolist() v2 = nodes[int(element[3])-1][2:5].tolist() @@ -592,7 +592,7 @@ class mdump: else: tris.append([element[0],evalue[type]] + list + n) # single square, convert to 2 tris, normals = up - + elif snap.eflag == 3: v1 = nodes[int(element[2])-1][2:5].tolist() v2 = nodes[int(element[3])-1][2:5].tolist() @@ -608,7 +608,7 @@ class mdump: else: tris.append([element[0],evalue[type]] + list + n) # single cube, convert to 12 tris, normals = out - + elif snap.eflag == 4: v1 = nodes[int(element[2])-1][2:5].tolist() v2 = nodes[int(element[3])-1][2:5].tolist() @@ -666,7 +666,7 @@ class mdump: n = normal(list[0:3],list[3:6],list[6:9]) if type == -1: tris.append([element[0],element[1]] + list + n) else: tris.append([element[0],evalue[type]] + list + n) - + lines = [] return time,box,atoms,bonds,tris,lines @@ -674,7 +674,7 @@ class mdump: # -------------------------------------------------------------------- # return lists of node/element info for snapshot isnap # if called with flag, then index is timestep, so convert to snapshot index - + def mviz(self,index,flag=0): if not flag: isnap = index else: @@ -686,14 +686,14 @@ class mdump: i += 1 isnap = i - 1 snap = self.snaps[isnap] - + time = snap.time box = [snap.xlo,snap.ylo,snap.zlo,snap.xhi,snap.yhi,snap.zhi] nvalues = [] if snap.nvalueflag: nvalues = snap.nvalues evalues = [] if snap.nvalueflag: evalues = snap.evalues - + return time,box,snap.nodes,snap.elements,nvalues,evalues # -------------------------------------------------------------------- @@ -718,7 +718,7 @@ class mdump: 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] - + # -------------------------------------------------------------------- def compare_atom(self,a,b): @@ -727,7 +727,7 @@ class mdump: elif a[0] > b[0]: return 1 else: - return 0 + return 0 # -------------------------------------------------------------------- # one snapshot @@ -742,7 +742,7 @@ class tselect: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def all(self): @@ -789,7 +789,7 @@ class tselect: data.nselect -= 1 data.eselect.all() print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) - + # -------------------------------------------------------------------- def test(self,teststr): @@ -835,7 +835,7 @@ class eselect: 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/src/pair.py b/src/pair.py index ad64d4b..80240f7 100644 --- a/src/pair.py +++ b/src/pair.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. # pair tool @@ -15,15 +15,15 @@ p = pair("lj/charmm/coul/charmm") create pair object for specific pair style available styles: lj/cut, lj/cut/coul/cut, lj/charmm/coul/charmm -p.coeff(d) extract pairwise coeffs from data object -p.init(cut1,cut2,...) setup based on coeffs and cutoffs +p.coeff(d) extract pairwise coeffs from data object +p.init(cut1,cut2,...) setup based on coeffs and cutoffs init args are specific to pair style: lj/cut = cutlj lj/cut/coul/cut = cutlj,cut_coul (cut_coul optional) lj/charmm/coul/charmm = cutlj_inner,cutlj,cutcoul_inner,cut_coul (last 2 optional) - + e_vdwl,e_coul = p.single(rsq,itype,jtype,q1,q2,...) compute LJ/Coul energy pairwise energy between 2 atoms at distance rsq with their attributes @@ -96,9 +96,9 @@ class pair: epsilon = data.get("Pair Coeffs",2) sigma = data.get("Pair Coeffs",3) ntypes = len(epsilon) - + self.lj3 = [] - self.lj4 = [] + self.lj4 = [] for i in xrange(ntypes): self.lj3.append(ntypes * [0]) self.lj4.append(ntypes * [0]) @@ -107,7 +107,7 @@ class pair: sigma_ij = sqrt(sigma[i]*sigma[j]) self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0); self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0); - + # -------------------------------------------------------------------- # args = cutlj @@ -117,20 +117,20 @@ class pair: # -------------------------------------------------------------------- # args = rsq,itype,jtype - + def single_lj_cut(self,list): rsq = list[0] itype = list[1] jtype = list[2] - + r2inv = 1.0/rsq - + if rsq < self.cut_ljsq: r6inv = r2inv*r2inv*r2inv eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype]) else: eng_vdwl = 0.0 - - return eng_vdwl + + return eng_vdwl # -------------------------------------------------------------------- # -------------------------------------------------------------------- @@ -140,9 +140,9 @@ class pair: epsilon = data.get("Pair Coeffs",2) sigma = data.get("Pair Coeffs",3) ntypes = len(epsilon) - + self.lj3 = [] - self.lj4 = [] + self.lj4 = [] for i in xrange(ntypes): self.lj3.append(ntypes * [0]) self.lj4.append(ntypes * [0]) @@ -151,41 +151,41 @@ class pair: sigma_ij = sqrt(sigma[i]*sigma[j]) self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0); self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0); - + # -------------------------------------------------------------------- # args = cutlj, cut_coul (cut_coul optional) def init_lj_cut_coul_cut(self,list): - self.qqr2e = 332.0636 # convert energy to kcal/mol + self.qqr2e = 332.0636 # convert energy to kcal/mol cut_lj = list[0] self.cut_ljsq = cut_lj*cut_lj - + if len(list) == 1: cut_coul = cut_lj else: cut_coul = list[1] self.cut_coulsq = cut_coul*cut_coul # -------------------------------------------------------------------- # args = rsq,itype,jtype,q1,q2 - + def single_lj_cut_coul_cut(self,list): rsq = list[0] itype = list[1] jtype = list[2] q1 = list[3] q2 = list[4] - + r2inv = 1.0/rsq - + if rsq < self.cut_coulsq: eng_coul = self.qqr2e * q1*q2*sqrt(r2inv) else: eng_coul = 0.0 - + if rsq < self.cut_ljsq: r6inv = r2inv*r2inv*r2inv eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype]) else: eng_vdwl = 0.0 - - return eng_coul,eng_vdwl - + + return eng_coul,eng_vdwl + # -------------------------------------------------------------------- # -------------------------------------------------------------------- # lj/charmm/coul/charmm methods @@ -194,9 +194,9 @@ class pair: epsilon = data.get("Pair Coeffs",2) sigma = data.get("Pair Coeffs",3) ntypes = len(epsilon) - + self.lj3 = [] - self.lj4 = [] + self.lj4 = [] for i in xrange(ntypes): self.lj3.append(ntypes * [0]) self.lj4.append(ntypes * [0]) @@ -205,18 +205,18 @@ class pair: sigma_ij = 0.5 * (sigma[i] + sigma[j]) self.lj3[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,12.0); self.lj4[i][j] = 4.0 * epsilon_ij * pow(sigma_ij,6.0); - + # -------------------------------------------------------------------- # args = cutlj_inner,cutlj,cutcoul_inner,cut_coul (last 2 optional) def init_lj_charmm_coul_charmm(self,list): - self.qqr2e = 332.0636 # convert energy to kcal/mol + self.qqr2e = 332.0636 # convert energy to kcal/mol cut_lj_inner = list[0] cut_lj = list[1] self.cut_lj_innersq = cut_lj_inner*cut_lj_inner self.cut_ljsq = cut_lj*cut_lj - + if len(list) == 2: cut_coul_inner = cut_lj_inner cut_coul = cut_lj @@ -235,16 +235,16 @@ class pair: # -------------------------------------------------------------------- # args = rsq,itype,jtype,q1,q2 - + def single_lj_charmm_coul_charmm(self,list): rsq = list[0] itype = list[1] jtype = list[2] q1 = list[3] q2 = list[4] - + r2inv = 1.0/rsq - + if rsq < self.cut_coulsq: eng_coul = self.qqr2e * q1*q2*sqrt(r2inv) if rsq > self.cut_coul_innersq: @@ -253,7 +253,7 @@ class pair: self.denom_coul eng_coul *= switch1 else: eng_coul = 0.0 - + if rsq < self.cut_ljsq: r6inv = r2inv*r2inv*r2inv eng_vdwl = r6inv*(self.lj3[itype][jtype]*r6inv-self.lj4[itype][jtype]) @@ -263,5 +263,5 @@ class pair: self.denom_lj eng_vdwl *= switch1 else: eng_vdwl = 0.0 - + return eng_coul,eng_vdwl diff --git a/src/patch.py b/src/patch.py index 2bf0cd6..6877d11 100644 --- a/src/patch.py +++ b/src/patch.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. # patch tool @@ -13,10 +13,10 @@ oneline = "Create patchy Lennard-Jones particles for LAMMPS input" docstr = """ p = patch(vfrac) setup box with a specified volume fraction p = patch(vfrac,1,1,2) x,y,z = aspect ratio of box (def = 1,1,1) - -p.seed = 48379 set random # seed (def = 12345) -p.randomized = 0 1 = choose next particle randomly, 0 = as generated -p.dim = 2 set dimension of created box (def = 3) + +p.seed = 48379 set random # seed (def = 12345) +p.randomized = 0 1 = choose next particle randomly, 0 = as generated +p.dim = 2 set dimension of created box (def = 3) p.blen = 0.97 set length of tether bonds (def = 0.97) p.dmin = 1.02 set min r from i-1 to i+1 tether site (def = 1.02) p.lattice = [Nx,Ny,Nz] generate Nx by Ny by Nz lattice of particles @@ -25,7 +25,7 @@ p.lattice = [Nx,Ny,Nz] generate Nx by Ny by Nz lattice of particles lattice = [0,0,0] = generate N particles randomly, default p.build(100,"hex2",1,2,3) create 100 "hex2" particles with types 1,2,3 - + can be invoked multiple times keywords: c60hex2: diam,1,2,3 = C-60 with 2 hex patches and ctr part, types 1,2,3 @@ -46,7 +46,7 @@ p.build(100,"hex2",1,2,3) create 100 "hex2" particles with types 1,2,3 from Alo to Ahi and Blo to Bhi, line type m linetri: Alo,Ahi,Blo,Bhi,m = 3-line 2d triangle with random base from Alo to Ahi and height Blo to Bhi, type m - + p.write("data.patch") write out system to LAMMPS data file """ @@ -69,7 +69,7 @@ from data import data # Class definition class patch: - + # -------------------------------------------------------------------- def __init__(self,vfrac,*list): @@ -179,7 +179,7 @@ class patch: #xp[0] = 1; xp[1] = 0; xp[2] = 0 #yp[0] = 0; yp[1] = 1; yp[2] = 0 #zp[0] = 0; zp[1] = 0; zp[2] = 1 - + # random origin or lattice site for new particle if latflag == 0: @@ -195,15 +195,15 @@ class patch: zorig = self.zlo + iz*self.zprd/self.lattice[2] #xorig = 0; yorig = 0; zorig = 0 - + # unpack bonds in molecule before atoms so idatom = all previous atoms - + for bond in molecule[1]: idbond += 1 bonds.append([idbond,bond[0],bond[1]+idatom+1,bond[2]+idatom+1]) # unpack triples in molecule as displacements from associated atom - + for triple in molecule[2]: triples.append(triple) # unpack atoms in molecule @@ -223,7 +223,7 @@ class patch: ix = iy = iz = 0 x,y,z,ix,iy,iz = self.pbc(x,y,z,ix,iy,iz) atoms.append([idatom,idmol,atom[0],x,y,z,ix,iy,iz]) - + elif self.style == "tri": for i,atom in enumerate(molecule[0]): idatom += 1 @@ -239,7 +239,7 @@ class patch: if not triples: triflag = 0 else: triflag = 1 atoms.append([idatom,idmol,atom[0],triflag,mass,x,y,z,ix,iy,iz]) - + if triflag: triple = triples[i] xtri = triple[0] @@ -272,7 +272,7 @@ class patch: list = [atom[2] for atom in atoms] atypes = max(list) - + d = data() d.title = "LAMMPS data file for Nanoparticles" d.headers["atoms"] = len(atoms) @@ -345,7 +345,7 @@ class patch: if self.lattice[0]*self.lattice[1] != len(self.molecules): raise StandardError,"lattice inconsistent with # of molecules" else: latflag = 0 - + idatom = idbond = idmol = 0 atoms = [] bonds = [] @@ -357,10 +357,10 @@ class patch: if self.randomized: i = int(self.random()*len(self.molecules)) else: i = 0 molecule = self.molecules.pop(i) - + idmol += 1 segments = [] - + # xp[2],yp[2] = randomly oriented, normalized basis vectors # xp is in random direction # yp is (0,0,1) crossed into xp @@ -387,15 +387,15 @@ class patch: xorig = self.xlo + ix*self.xprd/self.lattice[0] yorig = self.ylo + iy*self.yprd/self.lattice[1] zorig = 0.0 - + # unpack bonds in molecule before atoms so idatom = all previous atoms - + for bond in molecule[1]: idbond += 1 bonds.append([idbond,bond[0],bond[1]+idatom+1,bond[2]+idatom+1]) # unpack segments in molecule as displacements from associated atom - + for segment in molecule[3]: segments.append(segment) # unpack atoms in molecule @@ -445,7 +445,7 @@ class patch: segment[2],segment[3],tmp = \ self.pbc_near(segment[2],segment[3],0,x,y,z) lines.append([idatom] + segment) - + # create the data file list = [atom[2] for atom in atoms] @@ -464,7 +464,7 @@ class patch: d.headers["zlo zhi"] = (self.zlo,self.zhi) # atoms section of data file - + records = [] if self.style == "molecular": for atom in atoms: @@ -502,7 +502,7 @@ class patch: # -------------------------------------------------------------------- # adjust x,y,z to be inside periodic box - + def pbc(self,x,y,z,ix,iy,iz): if x < self.xlo: x += self.xprd @@ -526,7 +526,7 @@ class patch: # -------------------------------------------------------------------- # adjust xnew,ynew,znew to be near x,y,z in periodic sense - + def pbc_near(self,xnew,ynew,znew,x,y,z): if x-xnew > 0.5*self.xprd: xnew += self.xprd elif xnew-x > 0.5*self.xprd: xnew -= self.xprd @@ -540,7 +540,7 @@ class patch: # params = diam,type1,type2,type3 # type1 = type of non-patch atoms, type2 = type of patch atoms # type3 = type of center-of-sphere atom - + def c60hex2(self,*params): template = BUCKY_60 diam = params[0] @@ -549,19 +549,19 @@ class patch: atoms = make_sphere(template,diam,params[1],patches) volume = 4.0/3.0 * pi * diam*diam*diam/8 return atoms,[],[],[],volume - + # -------------------------------------------------------------------- # params = diam,type1,type2 # type1 = type of large center atom, type2 = type of hex patch atoms - + def hex2(self,*params): diam = params[0] type1 = params[1] type2 = params[2] - + atoms = [] atoms.append([type1,0.0,0.0,0.0]) - + atoms.append(atom_on_sphere(diam,type2,0.5*diam,0.0,0.0)) atoms.append(atom_on_sphere(diam,type2,0.5*diam,1.0,0.0)) atoms.append(atom_on_sphere(diam,type2,0.5*diam,-1.0,0.0)) @@ -580,19 +580,19 @@ class patch: volume = 4.0/3.0 * pi * diam*diam*diam/8 return atoms,[],[],[],volume - + # -------------------------------------------------------------------- # params = diam,type1,type2 # type1 = type of large center atom, type2 = type of hex patch atoms - + def hex4(self,*params): diam = params[0] type1 = params[1] type2 = params[2] - + atoms = [] atoms.append([type1,0.0,0.0,0.0]) - + atoms.append(atom_on_sphere(diam,type2,0.5*diam,0.0,0.0)) atoms.append(atom_on_sphere(diam,type2,0.5*diam,1.0,0.0)) atoms.append(atom_on_sphere(diam,type2,0.5*diam,-1.0,0.0)) @@ -632,13 +632,13 @@ class patch: # params = diam,nring,type1,type2 # nring = # of particles in ring # type1 = type of large center atom, type2 = type of ring atoms - + def ring(self,*params): diam = params[0] nring = params[1] type1 = params[2] type2 = params[3] - + atoms = [] atoms.append([type1,0.0,0.0,0.0]) @@ -655,7 +655,7 @@ class patch: # m12,m12type = length of tethers on each side of ball (m12 = 0 = no tether) # set three types of bonds: # 1 = big to small, 2 = small to small, 3 = across two tethers - + def ball(self,*params): diam = params[0] m1 = params[1] @@ -684,14 +684,14 @@ class patch: if i == 0: bonds.append([1,0,1+m1]) else: bonds.append([2,1+m1+i-1,1+m1+i]) if m1 and m2: bonds.append([3,1,m1+1]) - + volume = 4.0/3.0 * pi * diam*diam*diam/8 + (m1+m2)*pi/6.0 return atoms,bonds,[],[],volume # -------------------------------------------------------------------- # params = type1,type2 # type1 = type of 1/3 layers, type2 = type of middle layer - + def tri5(self,*params): template = TRI5_HOLLOW nlayer = 3 @@ -724,12 +724,12 @@ class patch: ntype = params[3] m1type = params[4] m2type = params[5] - + atoms = [] for i in range(n): x,y,z = i*self.blen,0.0,0.0 atoms.append([ntype,x,y,z]) - + if m1: atoms += tether(m1,m1type,self.blen,self.dmin, atoms[n-2],atoms[n-1],self.random) if m2: atoms += tether(m2,m2type,self.blen,self.dmin, @@ -742,7 +742,7 @@ class patch: for i in range(m2): if i == 0: bonds.append([1,0,n+m1]) else: bonds.append([1,n+m1+i-1,n+m1+i]) - + volume = (n+m1+m2) * pi / 6.0 return atoms,bonds,[],[],volume @@ -750,7 +750,7 @@ class patch: # params = nsize,m1,m2,m3,ntype,m1type,m2type,m3type # nsize,ntype = size,type of triangle center # m123,m123type = length of tethers on each corner (m123 = 0 = no tether) - + def tri(self,*params): nsize = params[0] m1 = params[1] @@ -760,7 +760,7 @@ class patch: m1type = params[5] m2type = params[6] m3type = params[7] - + atoms = [] for i in range(nsize): n = nsize - i @@ -794,7 +794,7 @@ class patch: # params = m1,m2,m3,m4,m5,m6,ntype,m1type,m2type,m3type,m4type,m5type,m6type # ntype = type of hex center # m123456,m123456type = length of tethers on each corner (m = 0 = no tether) - + def hex(self,*params): m1 = params[0] m2 = params[1] @@ -809,7 +809,7 @@ class patch: m4type = params[10] m5type = params[11] m6type = params[12] - + atoms = [] atoms.append([ntype,0.0,0.0,0.0]) atoms.append([ntype,self.blen,0.0,0.0]) @@ -818,7 +818,7 @@ class patch: atoms.append([ntype,-self.blen/2.0,self.blen*sqrt(3.0)/2.0,0.0]) atoms.append([ntype,self.blen/2.0,-self.blen*sqrt(3.0)/2.0,0.0]) atoms.append([ntype,-self.blen/2.0,-self.blen*sqrt(3.0)/2.0,0.0]) - + n = len(atoms) if m1: atoms += tether(m1,m1type,self.blen,self.dmin, atoms[0],atoms[1],self.random) @@ -864,7 +864,7 @@ class patch: def dimer(self,*params): sep = params[0] type = params[1] - + atoms = [] x,y,z = 0.0,0.0,0.0 atoms.append([type,x,y,z]) @@ -887,7 +887,7 @@ class patch: if n % 2 == 0: raise StandardError, "N in patch::star2d is not odd" middle = n/2 - + atoms = [] x,y,z = 0.0,0.0,0.0 atoms.append([type,x,y,z]) @@ -919,7 +919,7 @@ class patch: height = (m-1) * sep width = (n-1) * sep - + atoms = [] for i in range(n): x,y,z = i*sep,0.0,0.0 @@ -941,7 +941,7 @@ class patch: # params = a,type # a = edge length of tet # type = type of each vertex in tet - + def tritet(self,*params): a = params[0] type = params[1] @@ -985,7 +985,7 @@ class patch: # Blo to Bhi = bounds of edge length in y of box # Clo to Chi = bounds of edge length in y of box # type = type of each vertex in rectangular box - + def tribox(self,*params): alo = params[0] ahi = params[1] @@ -994,7 +994,7 @@ class patch: clo = params[4] chi = params[5] type = params[6] - + a = alo + self.random()*(ahi-alo) b = blo + self.random()*(bhi-blo) c = clo + self.random()*(chi-clo) @@ -1063,7 +1063,7 @@ class patch: # Alo to Ahi = bounds of edge length in x of rectangle # Blo to Bhi = bounds of edge length in y of rectangle # type = type of each line segment in rectangle - + def linebox(self,*params): alo = float(params[0]) ahi = float(params[1]) @@ -1073,7 +1073,7 @@ class patch: a = alo + self.random()*(ahi-alo) b = blo + self.random()*(bhi-blo) - + atoms = [] segments = [] @@ -1099,7 +1099,7 @@ class patch: # Alo to Ahi = bounds of base length in x of triangle # Blo to Bhi = bounds of heigth in y of isosceles triangle # type = type of each line segment in triangle - + def linetri(self,*params): alo = float(params[0]) ahi = float(params[1]) @@ -1109,7 +1109,7 @@ class patch: base = alo + self.random()*(ahi-alo) ht = blo + self.random()*(bhi-blo) - + atoms = [] segments = [] @@ -1208,7 +1208,7 @@ TRI5_HOLLOW = ((0,0,0),(1,0,0),(2,0,0),(3,0,0),(4,0,0), (1.0,2*sqrt(3)/2,0),(3.0,2*sqrt(3)/2,0), (1.5,3*sqrt(3)/2,0),(2.5,3*sqrt(3)/2,0), (2.0,4*sqrt(3)/2,0)) - + SIMPLE_7 = ((0,0,0),(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)) # C60 with added center point at end diff --git a/src/pdbfile.py b/src/pdbfile.py index 1713ada..9b2238c 100644 --- a/src/pdbfile.py +++ b/src/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. # pdb tool @@ -22,7 +22,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, ... @@ -36,7 +36,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) @@ -87,7 +87,7 @@ class pdbfile: # flist = full list of all PDB input file names # append .pdb if needed - + if filestr: list = filestr.split() flist = [] @@ -107,7 +107,7 @@ class pdbfile: raise StandardError, "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() @@ -117,7 +117,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: @@ -135,7 +135,7 @@ class pdbfile: f = open(file,'w') # use template PDB file with each snapshot - + if self.data: n = flag = 0 while 1: @@ -153,7 +153,7 @@ class pdbfile: print >>f,"END" print file, sys.stdout.flush() - + f.close() print "\nwrote %d datasets to %s in PDB format" % (n,file) @@ -189,7 +189,7 @@ class pdbfile: f = open(file,'w') self.convert(f,which) f.close() - + print time, sys.stdout.flush() n += 1 @@ -206,13 +206,13 @@ class pdbfile: else: file = root + str(n) file += ".pdb" - + f = open(file,'w') f.write(open(infile,'r').read()) f.close() print file, sys.stdout.flush() - + n += 1 print "\nwrote %d datasets to %s*.pdb in PDB format" % (n,root) @@ -239,7 +239,7 @@ class pdbfile: self.convert(f,which) else: f.write(open(self.files[time],'r').read()) - + f.close() # -------------------------------------------------------------------- @@ -258,8 +258,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/src/pizza.py b/src/pizza.py index 138c1c4..e92027c 100755 --- a/src/pizza.py +++ b/src/pizza.py @@ -5,7 +5,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. # Change log: @@ -26,12 +26,12 @@ type ? for help, CTRL-D to quit help = """ pizza.py switch arg(s) switch arg(s) ... - -s silent (else print start-up help) - -t log dump raster load only these tools - -x raster rasmol load all tools except these + -s silent (else print start-up help) + -t log dump raster load only these tools + -x raster rasmol load all tools except these -f mine.py arg1 arg2 run script file with args - -c "vec = range(100)" run Python command - -q quit (else interactive) + -c "vec = range(100)" run Python command + -q quit (else interactive) Everything typed at the ">" prompt is a Python command @@ -39,13 +39,13 @@ Additional commands available at ">" prompt: ? print help message ?? one-line for each tool and script ? raster list tool commands or script syntax - ?? energy.py full documentation of tool or script + ?? energy.py full documentation of tool or script !ls -l shell command @cd .. cd to a new directory @log tmp.log log all commands typed so far to file @run block.py arg1 arg2 run script file with args @time d = dump("*.dump") time a command - + Tools: """ @@ -86,13 +86,13 @@ def trap(type,value,tback): global argv # only check SyntaxErrors - + if not isinstance(value,exceptions.SyntaxError): sys.__excepthook__(type,value,tback) return - + # special commands at top level only, not in indented text entry - + if value.text[0].isspace(): sys.__excepthook__(type,value,tback) return @@ -105,19 +105,19 @@ def trap(type,value,tback): if value.text[0] == "?": words = value.text.split() - + if len(words) == 1 and words[0] == "?": print intro[1:] % version print help[1:]," ", for tool in tools: print tool, print - + elif len(words) == 1 and words[0] == "??": for tool in tools: exec "oneline = oneline_%s" % tool print "%-11s%s" % (tool,oneline) print - + scripts = [] for dir in PIZZA_SCRIPTS[1:]: list = glob.glob("%s/*.py" % dir) @@ -134,7 +134,7 @@ def trap(type,value,tback): if flag: doc = line[line.find("Purpose:")+8:] else: doc = " not available\n" print "%-20s%s" % (filename,doc), - + elif len(words) == 2 and words[0] == "?": if words[1][-3:] == ".py": fileflag = 0 @@ -154,7 +154,7 @@ def trap(type,value,tback): break if not fileflag: print "%s is not a recognized script" % words[1] - + else: if words[1] in tools: exec "txt = docstr_%s" % words[1] @@ -164,7 +164,7 @@ def trap(type,value,tback): print txt else: print "%s is not a recognized tool" % words[1] - + elif len(words) == 2 and words[0] == "??": if words[1][-3:] == ".py": fileflag = 0 @@ -187,7 +187,7 @@ def trap(type,value,tback): exec "print docstr_%s" % words[1] else: print "%s is not a recognized class" % words[1] - + return # shell command like !ls, !ls -l @@ -200,7 +200,7 @@ def trap(type,value,tback): # for run and time, use namespace in execfile and exec commands # else variables defined in script/command # won't be set in top-level Pizza.py - + if value.text[0] == "@": words = value.text.split() if words[0][1:] == "cd": @@ -237,9 +237,9 @@ def trap(type,value,tback): t2 = clock() print "CPU time = ",t2-t1 return - + # unrecognized command, let system handle error - + sys.__excepthook__(type,value,tback) # ------------------------------------------------------------------------- @@ -306,7 +306,7 @@ if not silent: print intro[1:] % version, if len(yes_tools) > 0 and len(no_tools) > 0: print "ERROR: cannot use -t and -x switches together" sys.exit() - + # ------------------------------------------------------------------------- # tools = list of tool names to import # if -t switch was used, tools = just those files @@ -424,7 +424,7 @@ for task in tasks: cmd = "" for arg in argv: cmd += arg + " " exec cmd - + # ------------------------------------------------------------------------- # store global namespace # swap in a new exception handler diff --git a/src/plotview.py b/src/plotview.py index 05abc3d..7cc83c5 100644 --- a/src/plotview.py +++ b/src/plotview.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. # plotview tool @@ -62,27 +62,27 @@ class plotview: self.plot = plot # create GUI - + from __main__ import tkroot root = Toplevel(tkroot) root.title('Pizza.py plotview tool') - + self.frame1 = Frame(root) self.frame2 = Frame(root) self.frame3 = Frame(root) - + Button(self.frame1,text="Print As:",command=self.save).pack(side=TOP) self.entry = Entry(self.frame1,width=16) self.entry.insert(0,"tmp") self.entry.pack(side=TOP) - + Label(self.frame2,text="Select").pack(side=LEFT) Label(self.frame2,text = "Display").pack(side=RIGHT) self.nplots = source.nvec self.names = source.names self.x = self.names[0] - + self.radiovar = IntVar() self.checkbuttons = [] self.checkvars = [] @@ -91,18 +91,18 @@ class plotview: # for each vector (not including 1st) # create a plot and title it # create a line in GUI with selection and check button - + for i in range(self.nplots): self.plot.select(i+1) self.plot.xtitle(self.x) self.plot.ytitle(self.names[i]) self.plot.title(self.names[i]) - + b = BooleanVar() b.set(0) self.checkvars.append(b) self.checkold.append(0) - + line = Frame(self.frame3) rtitle = "%d %s" % (i+1,self.names[i]) Radiobutton(line, text=rtitle, value=i+1, variable=self.radiovar, @@ -119,7 +119,7 @@ class plotview: # -------------------------------------------------------------------- # set radio button and checkbox - + def select(self,n): self.plot.select(n) self.radiovar.set(n) @@ -130,7 +130,7 @@ class plotview: def yes(self,n): if not self.checkvars[n-1].get(): self.checkbuttons[n-1].invoke() - + # -------------------------------------------------------------------- # only invoke if currently set @@ -143,7 +143,7 @@ class plotview: oldtext = self.entry.get() self.entry.delete(0,len(oldtext)) self.entry.insert(0,newtext) - + # -------------------------------------------------------------------- def save(self): @@ -154,7 +154,7 @@ class plotview: # -------------------------------------------------------------------- # called when any radio selection button is clicked - + def radioselect(self): self.select(self.radiovar.get()) @@ -163,7 +163,7 @@ class plotview: # draws or hides plot # loop is to find which checkbox changed status # grab x,y data to plot out of source object - + def check(self): for i in range(self.nplots): if int(self.checkvars[i].get()) != self.checkold[i]: diff --git a/src/rasmol.py b/src/rasmol.py index 513eae3..5427688 100644 --- a/src/rasmol.py +++ b/src/rasmol.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. # rasmol tool @@ -20,7 +20,7 @@ r.show(N,"my.rasmol") use file as RasMol script r.all() make images of all selected snapshots with def script r.all("my.rasmol") use file as RasMol script - + r.run(N) run RasMol interactivly on snapshot N r.run(N,"new.rasmol") adjust via mouse or RasMol commands r.run(N,"new.rasmol","old.rasmol") type quit to save RasMol script file @@ -57,7 +57,7 @@ class rasmol: def __init__(self,pdb): self.pdb = pdb self.file = "image" - + # -------------------------------------------------------------------- def start(self): @@ -88,7 +88,7 @@ class rasmol: def show(self,*list): # create tmp.pdb with atom data - + n = list[0] self.pdb.single(n,"tmp.pdb") @@ -103,21 +103,21 @@ class rasmol: rasmol_text = rasmol_template # write rasmol_text to tmp.rasmol, substituting tmp.pdb for filename - + f = open("tmp.rasmol","w") text = rasmol_text % "tmp.pdb" print >>f,text f.close() # run RasMol to create image in tmp.gif - + self.start() self.__call__("source tmp.rasmol") self.__call__("write tmp.gif") self.stop() # display the image - + cmd = "%s tmp.gif" % (PIZZA_DISPLAY) commands.getoutput(cmd) @@ -134,7 +134,7 @@ class rasmol: rasmol_text = re.sub('load pdb ".*"','load pdb "%s"',rasmol_text) else: rasmol_text = rasmol_template - + # iterate over all timesteps # write snapshot to tmpN.pdb # write RasMol input script to tmpN.rasmol @@ -193,7 +193,7 @@ class rasmol: commands.getoutput("rm tmp*.pdb") commands.getoutput("rm tmp*.rasmol") - + # -------------------------------------------------------------------- def run(self,*list): @@ -214,18 +214,18 @@ class rasmol: rasmol_text = rasmol_template # write rasmol_text to tmp.rasmol - + f = open("tmp.rasmol","w") text = rasmol_template % "tmp.pdb" print >>f,text f.close() # run RasMol to create image in tmp.gif - + self.start() self.__call__("source tmp.rasmol") self.enter() - + if len(list) > 1: newfile = list[1] else: newfile = "tmp.rasmol" self.__call__("write script %s" % newfile) diff --git a/src/raster.py b/src/raster.py index 8d6d03d..abdfb69 100644 --- a/src/raster.py +++ b/src/raster.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. # raster tool @@ -16,8 +16,8 @@ r = raster(d) create Raster3d wrapper for data in d d = atom snapshot object (dump, data) r.bg("black") set background color (def = "black") -r.size(N) set image size to NxN -r.size(N,M) set image size to NxM +r.size(N) set image size to NxN +r.size(N,M) set image size to NxM r.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) r.shift(x,y) translate by x,y pixels in view window (def = 0,0) r.zoom(0.5) scale image by factor (def = 1) @@ -27,7 +27,7 @@ r.box(0/1/2,"red",4) set box edge thickness r.file = "image" file prefix for created images (def = "image") r.show(N) show image of snapshot at timestep N - + r.all() make images of all selected snapshots r.all(P) images of all, start file label at P r.all(N,M,P) make M images of snapshot N, start label at P @@ -40,18 +40,18 @@ r.pan() no pan during all() (default) r.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() r.select = "" no extra aselect (default) - + %g varies from 0.0 to 1.0 from beginning to end of all() r.label(x,y,"h",size,"red","This is a label") add label to each image r.nolabel() delete all labels - + x,y coords = -0.5 to 0.5, "h" or "t" for Helvetica or Times font size = fontsize (e.g. 10), "red" = color of text - -r.acol(2,"green") set atom colors by atom type (1-N) -r.acol([2,4],["red","blue"]) 1st arg = one type or list of types -r.acol(0,"blue") 2nd arg = one color or list of colors + +r.acol(2,"green") set atom colors by atom type (1-N) +r.acol([2,4],["red","blue"]) 1st arg = one type or list of types +r.acol(0,"blue") 2nd arg = one color or list of colors r.acol(range(20),["red","blue"]) if list lengths unequal, interpolate r.acol(range(10),"loop") assign colors in loop, randomly ordered @@ -61,23 +61,23 @@ r.acol(range(10),"loop") assign colors in loop, randomly ordered r.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() -r.bcol() set bond color, same args as acol() -r.brad() set bond thickness, same args as arad() +r.bcol() set bond color, same args as acol() +r.brad() set bond thickness, same args as arad() -r.tcol() set triangle color, same args as acol() -r.tfill() set triangle fill, 0 fill, 1 line, 2 both +r.tcol() set triangle color, same args as acol() +r.tfill() set triangle fill, 0 fill, 1 line, 2 both r.lcol() set line color, same args as acol() r.lrad() set line thickness, same args as arad() r.adef() set atom/bond/tri/line properties to default -r.bdef() default = "loop" for colors, 0.45 for radii -r.tdef() default = 0.25 for bond/line thickness -r.ldef() default = 0 fill +r.bdef() default = "loop" for colors, 0.45 for radii +r.tdef() default = 0.25 for bond/line thickness +r.ldef() default = 0 fill by default 100 types are assigned if atom/bond/tri/line has type > # defined properties, is an error - + from vizinfo import colors access color list print colors list defined color names and RGB values colors["nickname"] = [R,G,B] set new RGB values from 0 to 255 @@ -136,7 +136,7 @@ class raster: self.scale = 1.0 self.xshift = self.yshift = 0 self.eye = 50.0 - + self.file = "image" self.boxflag = 0 self.bxcol = [1,1,0] @@ -158,14 +158,14 @@ class raster: from vizinfo import colors self.bgcol = [colors[color][0]/255.0,colors[color][1]/255.0, colors[color][2]/255.0] - + # -------------------------------------------------------------------- def size(self,xnew,ynew=None): self.xpixels = xnew if not ynew: self.ypixels = self.xpixels else: self.ypixels = ynew - + # -------------------------------------------------------------------- def rotate(self,ztheta,azphi): @@ -182,7 +182,7 @@ class raster: def zoom(self,scale): self.scale = scale; - + # -------------------------------------------------------------------- def box(self,*args): @@ -195,7 +195,7 @@ class raster: # -------------------------------------------------------------------- # scale down point-size by 3x - + def label(self,x,y,font,point,color,text): from vizinfo import colors scaledcolor = [colors[color][0]/255.0,colors[color][1]/255.0, @@ -207,12 +207,12 @@ class raster: def nolabel(self): self.labels = [] - + # -------------------------------------------------------------------- # show a single snapshot # distance from snapshot box or max box for all selected steps # always pre-call single() to re-center simulation data - + def show(self,ntime): data = self.data which = data.findtime(ntime) @@ -227,7 +227,7 @@ class raster: self.xtrans = float(nums[0][0]) self.ytrans = float(nums[0][1]) self.ztrans = float(nums[0][2]) - + self.single(0,self.file,box,atoms,bonds,tris,lines) cmd = "%s %s.png" % (PIZZA_DISPLAY,self.file) commands.getoutput(cmd) @@ -244,7 +244,7 @@ class raster: self.ztheta_stop = list[3] self.azphi_stop = list[4] self.scale_stop = list[5] - + # -------------------------------------------------------------------- def all(self,*list): @@ -275,7 +275,7 @@ class raster: if flag == -1: break fraction = float(i) / (ncount-1) - + if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,time) @@ -297,9 +297,9 @@ class raster: self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) - if n == nstart or self.panflag: + if n == nstart or self.panflag: self.xtrans = self.ytrans = self.ztrans = 0.0 - output = self.single(1,file,box,atoms,bonds,tris,lines) + output = self.single(1,file,box,atoms,bonds,tris,lines) nums = re.findall("translation to:\s*(\S*)\s*(\S*)\s*(\S*)\s",output) self.xtrans = float(nums[0][0]) self.ytrans = float(nums[0][1]) @@ -344,14 +344,14 @@ class raster: self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) - if n == nstart or self.panflag: + if n == nstart or self.panflag: self.xtrans = self.ytrans = self.ztrans = 0.0 - output = self.single(1,file,box,atoms,bonds,tris,lines) + output = self.single(1,file,box,atoms,bonds,tris,lines) nums = re.findall("translation to:\s*(\S*)\s*(\S*)\s*(\S*)\s",output) self.xtrans = float(nums[0][0]) self.ytrans = float(nums[0][1]) self.ztrans = float(nums[0][2]) - + self.single(0,file,box,atoms,bonds,tris,lines) print n, sys.stdout.flush() @@ -379,7 +379,7 @@ class raster: # draw box if boxflag or flag is set # flag = 1 is a pre-call of single to set frame size correctly # this keeps the view fixed even if atoms move around - + if self.boxflag or flag: box_write(f,box,self.bxcol,self.bxthick) ncolor = self.vizinfo.nacolor @@ -421,7 +421,7 @@ class raster: if itype > ncolor: raise StandardError,"line type too big" color = self.vizinfo.lcolor[itype] thick = self.vizinfo.lrad[itype] - print >>f,3 + print >>f,3 print >>f,line[2],line[3],line[4],thick, \ line[5],line[6],line[7],thick,color[0],color[1],color[2] @@ -432,7 +432,7 @@ class raster: print >>f,11 print >>f,label[0],label[1],0.0,label[5][0],label[5][1],label[5][2] print >>f,label[6] - + f.close() if len(self.labels) == 0: @@ -448,45 +448,45 @@ class raster: def adef(self): self.vizinfo.setcolors("atom",range(100),"loop") self.vizinfo.setradii("atom",range(100),0.45) - + # -------------------------------------------------------------------- def bdef(self): self.vizinfo.setcolors("bond",range(100),"loop") self.vizinfo.setradii("bond",range(100),0.25) - + # -------------------------------------------------------------------- def tdef(self): - self.vizinfo.setcolors("tri",range(100),"loop") - self.vizinfo.setfills("tri",range(100),0) + self.vizinfo.setcolors("tri",range(100),"loop") + self.vizinfo.setfills("tri",range(100),0) # -------------------------------------------------------------------- def ldef(self): - self.vizinfo.setcolors("line",range(100),"loop") + self.vizinfo.setcolors("line",range(100),"loop") self.vizinfo.setradii("line",range(100),0.25) - + # -------------------------------------------------------------------- def acol(self,atypes,colors): self.vizinfo.setcolors("atom",atypes,colors) - + # -------------------------------------------------------------------- def arad(self,atypes,radii): - self.vizinfo.setradii("atom",atypes,radii) - + self.vizinfo.setradii("atom",atypes,radii) + # -------------------------------------------------------------------- def bcol(self,btypes,colors): self.vizinfo.setcolors("bond",btypes,colors) - + # -------------------------------------------------------------------- def brad(self,btypes,radii): self.vizinfo.setradii("bond",btypes,radii) - + # -------------------------------------------------------------------- def tcol(self,ttypes,colors): @@ -525,34 +525,34 @@ def box_write(f,box,color,thick): red = color[0] green = color[1] blue = color[2] - - print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,ylo,zlo,thick,xhi,ylo,zlo,thick,red,green,blue) - print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,yhi,zlo,thick,xhi,yhi,zlo,thick,red,green,blue) - print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,ylo,zhi,thick,xhi,ylo,zhi,thick,red,green,blue) - print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,yhi,zhi,thick,xhi,yhi,zhi,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,ylo,zlo,thick,xlo,yhi,zlo,thick,red,green,blue) + (xlo,ylo,zlo,thick,xhi,ylo,zlo,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xhi,ylo,zlo,thick,xhi,yhi,zlo,thick,red,green,blue) + (xlo,yhi,zlo,thick,xhi,yhi,zlo,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,ylo,zhi,thick,xlo,yhi,zhi,thick,red,green,blue) + (xlo,ylo,zhi,thick,xhi,ylo,zhi,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xhi,ylo,zhi,thick,xhi,yhi,zhi,thick,red,green,blue) + (xlo,yhi,zhi,thick,xhi,yhi,zhi,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,ylo,zlo,thick,xlo,ylo,zhi,thick,red,green,blue) + (xlo,ylo,zlo,thick,xlo,yhi,zlo,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xhi,ylo,zlo,thick,xhi,ylo,zhi,thick,red,green,blue) + (xhi,ylo,zlo,thick,xhi,yhi,zlo,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xlo,yhi,zlo,thick,xlo,yhi,zhi,thick,red,green,blue) + (xlo,ylo,zhi,thick,xlo,yhi,zhi,thick,red,green,blue) print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ - (xhi,yhi,zlo,thick,xhi,yhi,zhi,thick,red,green,blue) - + (xhi,ylo,zhi,thick,xhi,yhi,zhi,thick,red,green,blue) + + print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ + (xlo,ylo,zlo,thick,xlo,ylo,zhi,thick,red,green,blue) + print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ + (xhi,ylo,zlo,thick,xhi,ylo,zhi,thick,red,green,blue) + print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ + (xlo,yhi,zlo,thick,xlo,yhi,zhi,thick,red,green,blue) + print >>f,"3\n%g %g %g %g %g %g %g %g %g %g %g" % \ + (xhi,yhi,zlo,thick,xhi,yhi,zhi,thick,red,green,blue) + # -------------------------------------------------------------------- # compute 3x3 rotation matrix for viewing angle @@ -599,9 +599,9 @@ def rotation_matrix(coord1,angle1,coord2,angle2): a33 = 1.0 a12 = a22 = c1 a12 = s1; a21 = -s1 - + # 2nd rotation matrix - + b11 = b12 = b13 = b21 = b22 = b23 = b31 = b32 = b33 = 0.0 if coord2 == 'x': b11 = 1.0 @@ -615,26 +615,26 @@ def rotation_matrix(coord1,angle1,coord2,angle2): b33 = 1.0 b11 = b22 = c2 b12 = s2; b21 = -s2 - + # full matrix c = b*a - + c11 = b11*a11 + b12*a21 + b13*a31 c12 = b11*a12 + b12*a22 + b13*a32 c13 = b11*a13 + b12*a23 + b13*a33 - + c21 = b21*a11 + b22*a21 + b23*a31 c22 = b21*a12 + b22*a22 + b23*a32 c23 = b21*a13 + b22*a23 + b23*a33 - + c31 = b31*a11 + b32*a21 + b33*a31 c32 = b31*a12 + b32*a22 + b33*a32 c33 = b31*a13 + b32*a23 + b33*a33 - + # form rotation matrix # each line padded with 0.0 for 4x4 raster3d matrix matrix = "%g %g %g 0.0\n%g %g %g 0.0\n%g %g %g 0.0" % \ - (c11,c12,c13,c21,c22,c23,c31,c32,c33) + (c11,c12,c13,c21,c22,c23,c31,c32,c33) return matrix # -------------------------------------------------------------------- @@ -653,7 +653,7 @@ F no, shadowed rods look funny 0.25 specular reflection component %g eye position 1 1 1 main light source postion -%s +%s %g %g %g %g 3 mixed object types * diff --git a/src/svg.py b/src/svg.py index 747db72..db1c099 100644 --- a/src/svg.py +++ b/src/svg.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. # svg tool @@ -16,8 +16,8 @@ s = svg(d) create SVG object for data in d d = atom snapshot object (dump, data) s.bg("black") set background color (def = "black") -s.size(N) set image size to NxN -s.size(N,M) set image size to NxM +s.size(N) set image size to NxN +s.size(N,M) set image size to NxM s.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) s.shift(x,y) translate by x,y pixels in view window (def = 0,0) s.zoom(0.5) scale image by factor (def = 1) @@ -40,18 +40,18 @@ s.pan() no pan during all() (default) s.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() s.select = "" no extra aselect (default) - + %g varies from 0.0 to 1.0 from beginning to end of all() s.label(x,y,"h",size,"red","This is a label") add label to each image s.nolabel() delete all labels - + x,y coords = -0.5 to 0.5, "h" or "t" for Helvetica or Times font size = fontsize (e.g. 10), "red" = color of text - -s.acol(2,"green") set atom colors by atom type (1-N) -s.acol([2,4],["red","blue"]) 1st arg = one type or list of types -s.acol(0,"blue") 2nd arg = one color or list of colors + +s.acol(2,"green") set atom colors by atom type (1-N) +s.acol([2,4],["red","blue"]) 1st arg = one type or list of types +s.acol(0,"blue") 2nd arg = one color or list of colors s.acol(range(20),["red","blue"]) if list lengths unequal, interpolate s.acol(range(10),"loop") assign colors in loop, randomly ordered @@ -61,19 +61,19 @@ s.acol(range(10),"loop") assign colors in loop, randomly ordered s.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() -s.bcol() set bond color, same args as acol() -s.brad() set bond thickness, same args as arad() +s.bcol() set bond color, same args as acol() +s.brad() set bond thickness, same args as arad() -s.tcol() set triangle color, same args as acol() -s.tfill() set triangle fill, 0 fill, 1 line, 2 both +s.tcol() set triangle color, same args as acol() +s.tfill() set triangle fill, 0 fill, 1 line, 2 both s.lcol() set line color, same args as acol() s.lrad() set line thickness, same args as arad() s.adef() set atom/bond/tri/line properties to default -s.bdef() default = "loop" for colors, 0.45 for radii -s.tdef() default = 0.25 for bond/line thickness -s.ldef() default = 0 fill +s.bdef() default = "loop" for colors, 0.45 for radii +s.tdef() default = 0.25 for bond/line thickness +s.ldef() default = 0 fill by default 100 types are assigned if atom/bond/tri/line has type > # defined properties, is an error @@ -149,21 +149,21 @@ class svg: self.bdef() self.tdef() self.ldef() - + # -------------------------------------------------------------------- def bg(self,color): from vizinfo import colors self.bgcol = [colors[color][0]/255.0,colors[color][1]/255.0, colors[color][2]/255.0] - + # -------------------------------------------------------------------- def size(self,newx,newy=None): self.xpixels = xnew if not ynew: self.ypixels = self.xpixels else: self.ypixels = ynew - + # -------------------------------------------------------------------- def rotate(self,ztheta,azphi): @@ -175,7 +175,7 @@ class svg: def shift(self,x,y): self.xshift = x; self.yshift = y; - + # -------------------------------------------------------------------- def box(self,*args): @@ -185,12 +185,12 @@ class svg: self.bxcol = [colors[args[1]][0]/255.0,colors[args[1]][1]/255.0, colors[args[1]][2]/255.0] if len(args) > 2: self.bxthick = args[2] - + # -------------------------------------------------------------------- def zoom(self,factor): self.scale = factor - + # -------------------------------------------------------------------- def show(self,ntime): @@ -199,11 +199,11 @@ class svg: time,box,atoms,bonds,tris,lines = data.viz(which) if self.boxflag == 2: box = data.maxbox() self.distance = compute_distance(box) - + self.single(self.file,box,atoms,bonds,tris,lines,1) cmd = "%s %s.svg" % (PIZZA_DISPLAY,self.file) commands.getoutput(cmd) - + # -------------------------------------------------------------------- def pan(self,*list): @@ -215,8 +215,8 @@ class svg: self.scale_start = list[2] self.ztheta_stop = list[3] self.azphi_stop = list[4] - self.scale_stop = list[5] - + self.scale_stop = list[5] + # -------------------------------------------------------------------- def all(self,*list): @@ -247,11 +247,11 @@ class svg: if flag == -1: break fraction = float(i) / (ncount-1) - - if self.select != "": + + if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,time) - time,boxone,atoms,bonds,tris,lines = data.viz(which) + time,boxone,atoms,bonds,tris,lines = data.viz(which) if self.boxflag < 2: box = boxone if n == nstart: self.distance = compute_distance(box) @@ -260,19 +260,19 @@ class svg: elif n < 100: file = self.file + "00" + str(n) elif n < 1000: file = self.file + "0" + str(n) else: file = self.file + str(n) - - if self.panflag: + + if self.panflag: self.ztheta = self.ztheta_start + \ fraction*(self.ztheta_stop - self.ztheta_start) self.azphi = self.azphi_start + \ fraction*(self.azphi_stop - self.azphi_start) self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) - - scaleflag = 0 - if n == nstart or self.panflag: scaleflag = 1 - self.single(file,box,atoms,bonds,tris,lines,scaleflag) + scaleflag = 0 + if n == nstart or self.panflag: scaleflag = 1 + + self.single(file,box,atoms,bonds,tris,lines,scaleflag) print time, sys.stdout.flush() i += 1 @@ -289,8 +289,8 @@ class svg: n = nstart for i in range(ncount): fraction = float(i) / (ncount-1) - - if self.select != "": + + if self.select != "": newstr = self.select % fraction data.aselect.test(newstr,ntime) time,boxone,atoms,bonds,tris,lines = data.viz(which) @@ -302,25 +302,25 @@ class svg: elif n < 100: file = self.file + "00" + str(n) elif n < 1000: file = self.file + "0" + str(n) else: file = self.file + str(n) - - if self.panflag: + + if self.panflag: self.ztheta = self.ztheta_start + \ fraction*(self.ztheta_stop - self.ztheta_start) self.azphi = self.azphi_start + \ fraction*(self.azphi_stop - self.azphi_start) - self.scale = self.scale_start + \ + self.scale = self.scale_start + \ fraction*(self.scale_stop - self.scale_start) scaleflag = 0 - if n == nstart or self.panflag: scaleflag = 1 + if n == nstart or self.panflag: scaleflag = 1 - self.single(file,box,atoms,bonds,tris,lines,scaleflag) + self.single(file,box,atoms,bonds,tris,lines,scaleflag) print n, sys.stdout.flush() n += 1 - - print "\n%d images" % ncount - + + print "\n%d images" % ncount + # -------------------------------------------------------------------- def label(self,x,y,font,point,color,text): @@ -334,11 +334,11 @@ class svg: def nolabel(self): self.labels = [] - + # -------------------------------------------------------------------- def single(self,file,box,atoms,bonds,tris,lines,scaleflag): - + matrix = rotation_matrix('x',-self.ztheta,'z',270.0-self.azphi) if scaleflag: self.factor = self.xpixels*self.scale / (1.6*self.distance) @@ -359,31 +359,31 @@ class svg: tri[0] = 1 newtri = self.transform(tri,matrix) olist.append(newtri) - + bound = 0.25 * self.distance for bond in bonds: newbond = [2,bond[1]] dx = bond[5] - bond[2] dy = bond[6] - bond[3] - dz = bond[7] - bond[4] + dz = bond[7] - bond[4] r = sqrt(dx*dx+dy*dy+dz*dz) if not r: r = 1 rad = self.vizinfo.arad[int(bond[9])] newbond.append(bond[2] + (r/r - rad/r) * dx) newbond.append(bond[3] + (r/r - rad/r) * dy) newbond.append(bond[4] + (r/r - rad/r) * dz) - + # cut off second side of bond - + dx = bond[2] - bond[5] dy = bond[3] - bond[6] - dz = bond[4] - bond[7] + dz = bond[4] - bond[7] r = sqrt(dx*dx+dy*dy+dz*dz) if not r: r = 1 rad = self.vizinfo.arad[int(bond[8])] newbond.append(bond[5] + (r/r - rad/r) * dx) newbond.append(bond[6] + (r/r - rad/r) * dy) - newbond.append(bond[7] + (r/r - rad/r) * dz) + newbond.append(bond[7] + (r/r - rad/r) * dz) if fabs(newbond[2]-newbond[5]) > bound or \ fabs(newbond[3]-newbond[6]) > bound: continue @@ -391,7 +391,7 @@ class svg: newbond = self.transform(newbond,matrix) if newbond[4] < newbond[7]: newbond[4] = newbond[7] olist.append(newbond) - + for line in lines: line[0] = 3 newline = self.transform(line,matrix) @@ -425,8 +425,8 @@ class svg: # write SVG file file += ".svg" - f = open(file,"w") - + f = open(file,"w") + header = ' ' % \ (self.ypixels,self.xpixels) header += '' @@ -437,15 +437,15 @@ class svg: color += 'fill="rgb(%s,%s,%s)"/>' % \ (self.bgcol[0]*255,self.bgcol[1]*255,self.bgcol[2]*255) print >>f,color - + for element in olist: self.write(f,0,element) for label in self.labels: self.write(f,1,label) footer = "" - print >> f,footer - - f.close() - + print >> f,footer + + f.close() + # -------------------------------------------------------------------- # rotate with matrix @@ -476,9 +476,9 @@ class svg: onew.append(matrix[0]*obj[5] + matrix[3]*obj[6] + matrix[6]*obj[7]) onew.append(matrix[1]*obj[5] + matrix[4]*obj[6] + matrix[7]*obj[7]) onew.append(matrix[2]*obj[5] + matrix[5]*obj[6] + matrix[8]*obj[7]) - + return onew - + # -------------------------------------------------------------------- def convert(self,objlist): @@ -504,26 +504,26 @@ class svg: obj[3] = yctr - factor*(obj[3] - offsety) obj[5] = factor*(obj[5] - offsetx) + xctr obj[6] = yctr - factor*(obj[6] - offsety) - + # -------------------------------------------------------------------- def write(self,f,flag,*args): if len(args): obj = args[0] - + if flag == 0: if obj[0] == 0: # atom with its color and radius - itype = int(obj[1]) - if itype > self.vizinfo.nacolor: + itype = int(obj[1]) + if itype > self.vizinfo.nacolor: raise StandardError,"atom type too big" color = self.vizinfo.acolor[itype] rad = self.vizinfo.arad[itype] - print >>f,'' % \ + print >>f,'' % \ (obj[2],obj[3],rad*self.factor, color[0]*255,color[1]*255,color[2]*255,self.thick) - + elif obj[0] == 1: # tri with its color (need to add fill type) - itype = int(obj[1]) - if itype > self.vizinfo.ntcolor: + itype = int(obj[1]) + if itype > self.vizinfo.ntcolor: raise StandardError,"tri type too big" color = self.vizinfo.tcolor[itype] print >>f,'' % \ @@ -532,17 +532,17 @@ class svg: elif obj[0] == 2: # bond with its color and thickness itype = int(obj[1]) - if itype > self.vizinfo.nbcolor: + if itype > self.vizinfo.nbcolor: raise StandardError,"bond type too big" color = self.vizinfo.bcolor[itype] thick = self.vizinfo.brad[itype] - print >>f,'' % \ + print >>f,'' % \ (obj[2],obj[3],obj[5],obj[6], color[0]*255,color[1]*255,color[2]*255,thick*self.factor) - + elif obj[0] == 3: # line with its color and thickness itype = int(obj[1]) - if itype > self.vizinfo.nlcolor: + if itype > self.vizinfo.nlcolor: raise StandardError,"line type too big" color = self.vizinfo.lcolor[itype] thick = self.vizinfo.lrad[itype] @@ -564,56 +564,56 @@ class svg: print >>f,' %s ' % \ (x,y,obj[3],obj[2],color[0]*255,color[1]*255,color[2]*255, color[0]*255,color[1]*255,color[2]*255,obj[5]) - + # -------------------------------------------------------------------- def adef(self): self.vizinfo.setcolors("atom",range(100),"loop") self.vizinfo.setradii("atom",range(100),0.45) - + # -------------------------------------------------------------------- def bdef(self): self.vizinfo.setcolors("bond",range(100),"loop") self.vizinfo.setradii("bond",range(100),0.25) - + # -------------------------------------------------------------------- def tdef(self): - self.vizinfo.setcolors("tri",range(100),"loop") - self.vizinfo.setfills("tri",range(100),0) + self.vizinfo.setcolors("tri",range(100),"loop") + self.vizinfo.setfills("tri",range(100),0) # -------------------------------------------------------------------- def ldef(self): - self.vizinfo.setcolors("line",range(100),"loop") + self.vizinfo.setcolors("line",range(100),"loop") self.vizinfo.setradii("line",range(100),0.25) # -------------------------------------------------------------------- def acol(self,atypes,colors): self.vizinfo.setcolors("atom",atypes,colors) - + # -------------------------------------------------------------------- def arad(self,atypes,radii): - self.vizinfo.setradii("atom",atypes,radii) - + self.vizinfo.setradii("atom",atypes,radii) + # -------------------------------------------------------------------- def bcol(self,btypes,colors): self.vizinfo.setcolors("bond",btypes,colors) - + # -------------------------------------------------------------------- def brad(self,btypes,radii): self.vizinfo.setradii("bond",btypes,radii) - + # -------------------------------------------------------------------- def tcol(self,ttypes,colors): self.vizinfo.setcolors("tri",ttypes,colors) - + # -------------------------------------------------------------------- def tfill(self,ttypes,flags): @@ -696,9 +696,9 @@ def rotation_matrix(coord1,angle1,coord2,angle2): a33 = 1.0 a12 = a22 = c1 a12 = s1; a21 = -s1 - + # 2nd rotation matrix - + b11 = b12 = b13 = b21 = b22 = b23 = b31 = b32 = b33 = 0.0 if coord2 == 'x': b11 = 1.0 @@ -712,23 +712,23 @@ def rotation_matrix(coord1,angle1,coord2,angle2): b33 = 1.0 b11 = b22 = c2 b12 = s2; b21 = -s2 - + # full matrix c = b*a - + c11 = b11*a11 + b12*a21 + b13*a31 c12 = b11*a12 + b12*a22 + b13*a32 c13 = b11*a13 + b12*a23 + b13*a33 - + c21 = b21*a11 + b22*a21 + b23*a31 c22 = b21*a12 + b22*a22 + b23*a32 c23 = b21*a13 + b22*a23 + b23*a33 - + c31 = b31*a11 + b32*a21 + b33*a31 c32 = b31*a12 + b32*a22 + b33*a32 c33 = b31*a13 + b32*a23 + b33*a33 - + # form rotation matrix - + matrix = (c11,c12,c13,c21,c22,c23,c31,c32,c33) return matrix diff --git a/src/tdump.py b/src/tdump.py index bcfd904..01f0464 100644 --- a/src/tdump.py +++ b/src/tdump.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. # tdump tool @@ -12,14 +12,14 @@ oneline = "Read dump files with triangle info" docstr = """ t = tdump("dump.one") read in one or more dump files -t = tdump("dump.1 dump.2.gz") can be gzipped -t = tdump("dump.*") wildcard expands to multiple files -t = tdump("dump.*",0) two args = store filenames, but don't read +t = tdump("dump.1 dump.2.gz") can be gzipped +t = tdump("dump.*") wildcard expands to multiple files +t = tdump("dump.*",0) two args = store filenames, but don't read incomplete and duplicate snapshots are deleted no column name assignment is performed -time = t.next() read next snapshot from dump files +time = t.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 @@ -43,7 +43,7 @@ time,box,atoms,bonds,tris,lines = t.viz(index) return list of viz objects id,type are from associated atom lines = NULL -t.owrap(...) wrap tris to same image as their atoms +t.owrap(...) wrap tris to same image as their atoms owrap() is called by dump tool's owrap() useful for wrapping all molecule's atoms/tris the same so it is contiguous @@ -101,7 +101,7 @@ class tdump: for word in words: self.flist += glob.glob(word) if len(self.flist) == 0 and len(list) == 1: raise StandardError,"no ldump file specified" - + if len(list) == 1: self.increment = 0 self.read_all() @@ -156,15 +156,15 @@ class tdump: snap = self.read_snapshot(f) if not snap: self.nextfile += 1 - if self.nextfile == len(self.flist): return -1 + if self.nextfile == len(self.flist): return -1 f.close() - self.eof = 0 - continue + self.eof = 0 + continue self.eof = f.tell() f.close() try: self.findtime(snap.time) - continue + continue except: break self.snaps.append(snap) @@ -176,7 +176,7 @@ class tdump: # -------------------------------------------------------------------- # read a single snapshot from file f # return snapshot or 0 if failed - + def read_snapshot(self,f): try: snap = Snap() @@ -217,7 +217,7 @@ class tdump: # -------------------------------------------------------------------- # map atom column names - + def map(self,*pairs): if len(pairs) % 2 != 0: raise StandardError, "tdump map() requires pairs of mappings" @@ -248,7 +248,7 @@ class tdump: return 0 # -------------------------------------------------------------------- - + def findtime(self,n): for i in xrange(self.nsnaps): if self.snaps[i].time == n: return i @@ -264,7 +264,7 @@ class tdump: del self.snaps[i] else: i += 1 - + # -------------------------------------------------------------------- # return list of lines to viz for snapshot isnap # if called with flag, then index is timestep, so convert to snapshot index @@ -295,9 +295,9 @@ class tdump: corner3y = self.names["corner3y"] corner3z = self.names["corner3z"] - # create line list from id,type,corner1x,...corner3z - # don't add line if all 4 values are 0 since not a line - + # create tris list from id,type,corner1x,...corner3z + # don't add tri if all 4 values are 0 since not a line + tris = [] for i in xrange(snap.natoms): atom = snap.atoms[i] @@ -334,7 +334,7 @@ class tdump: # idump = index of my line I in dump's atoms # jdump = atom J in dump's atoms that atom I was owrapped on # delx,dely = offset applied to atom I and thus to line I - + for i in xrange(snap.natoms): tag = atoms[i][id] idump = idsdump[tag] @@ -366,20 +366,20 @@ def normal(x,y,z): v1[0] = y[0] - x[0] v1[1] = y[1] - x[1] v1[2] = y[2] - x[2] - + v2 = 3*[0] v2[0] = z[0] - y[0] v2[1] = z[1] - y[1] v2[2] = z[2] - y[2] - + n = 3*[0] n[0] = v1[1]*v2[2] - v1[2]*v2[1] n[1] = v1[2]*v2[0] - v1[0]*v2[2] n[2] = v1[0]*v2[1] - v1[1]*v2[0] - + length = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]) n[0] /= length n[1] /= length n[2] /= length - + return n diff --git a/src/vcr.py b/src/vcr.py index 31ba42a..fd5c992 100644 --- a/src/vcr.py +++ b/src/vcr.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. # vcr tool @@ -16,26 +16,26 @@ v.add(gl) add a gl window to vcr GUI Actions (same as GUI widgets): -v.first() go to first frame +v.first() go to first frame v.prev() go to previous frame v.back() play backwards from current frame to start -v.stop() stop on current frame +v.stop() stop on current frame v.play() play from current frame to end -v.next() go to next frame -v.last() go to last frame +v.next() go to next frame +v.last() go to last frame -v.frame(31) set frame slider -v.delay(0.4) set delay slider -v.q(5) set quality slider +v.frame(31) set frame slider +v.delay(0.4) set delay slider +v.q(5) set quality slider -v.xaxis() view scene from x axis -v.yaxis() view scene from y axis -v.zaxis() view scene from z axis -v.box() toggle bounding box -v.axis() toggle display of xyz axes -v.norm() recenter and resize the view -v.ortho() toggle ortho/perspective button -v.reload() reload all frames from gl viewer data files +v.xaxis() view scene from x axis +v.yaxis() view scene from y axis +v.zaxis() view scene from z axis +v.box() toggle bounding box +v.axis() toggle display of xyz axes +v.norm() recenter and resize the view +v.ortho() toggle ortho/perspective button +v.reload() reload all frames from gl viewer data files v.clipxlo(0.2) clip scene at x lo fraction of box v.clipxhi(1.0) clip at x hi @@ -44,9 +44,9 @@ v.clipyhi(1.0) v.clipzlo(0.2) clip in z v.clipzhi(1.0) -v.save() save current scene to file.png -v.file("image") set filename -v.saveall() toggle save-all checkbox +v.save() save current scene to file.png +v.file("image") set filename +v.saveall() toggle save-all checkbox """ # History @@ -108,7 +108,7 @@ class vcr: Button(frame1,text=">",command=self.next).pack(side=LEFT) Button(frame1,text=">>",command=self.last).pack(side=LEFT) frame1.pack() - + frame2 = Frame(root) self.slider_frame = Scale(frame2,from_=0,to=self.nframes-1, command=self.frame,orient=HORIZONTAL, @@ -121,7 +121,7 @@ class vcr: command=self.delay,orient=HORIZONTAL, label=" Delay") self.slider_frame.pack(side=LEFT) - self.slider_quality.pack(side=LEFT) + self.slider_quality.pack(side=LEFT) self.slider_delay.pack(side=LEFT) frame2.pack() @@ -177,7 +177,7 @@ class vcr: Button(frame6,text="Save As:",command=self.save).pack(side=LEFT) self.entry_file = Entry(frame6,width = 16) self.entry_file.insert(0,"image") - self.entry_file.pack(side=LEFT) + self.entry_file.pack(side=LEFT) self.button_save = Checkbutton(frame6,text="SaveAll",command=self.saveall) self.button_save.pack(side=LEFT) frame6.pack() @@ -197,7 +197,7 @@ class vcr: frame8.pack() # display 1st image - + self.index = 0 self.display() @@ -219,25 +219,25 @@ class vcr: def last(self): self.index = self.nframes - 1 self.slider_frame.set(self.index) - + # -------------------------------------------------------------------- def next(self): if self.index < self.nframes - 1: self.index += 1 self.slider_frame.set(self.index) - + # -------------------------------------------------------------------- def previous(self): if self.index > 0: self.index -= 1 self.slider_frame.set(self.index) - + # -------------------------------------------------------------------- # play backward loop # disable GL caching while animating - + def back(self): if self.loop_flag != 0: return self.loop_flag = -1 @@ -248,7 +248,7 @@ class vcr: self.tkroot.update() if self.saveflag: self.saveloop(0) self.tkroot.after(self.delay_msec,self.loop) - + # -------------------------------------------------------------------- # play forward loop # disable GL caching while animating @@ -271,11 +271,11 @@ class vcr: def stop(self): self.loop_flag = 0 for view in self.viewlist: view.cache = 1 - + # -------------------------------------------------------------------- # loop forward or back until end of animation # if save flag is set, change file name and save snapshots - + def loop(self): if self.loop_flag == 1 and self.index == self.nframes - 1: self.loop_flag = 0 @@ -292,7 +292,7 @@ class vcr: # since slider_frame already called it # but seems to be necessary before GL save() saves window to file # else get previous image saved into file at each frame - + if self.saveflag: for view in self.viewlist: time,natoms = view.display(self.index) self.saveloop(1) @@ -309,7 +309,7 @@ class vcr: def frame(self,value): self.index = int(value) self.display() - + # -------------------------------------------------------------------- def delay(self,value): @@ -357,7 +357,7 @@ class vcr: else: self.boxflag = 1 for view in self.viewlist: view.box(1) - + # -------------------------------------------------------------------- def axis(self): @@ -367,7 +367,7 @@ class vcr: else: self.axisflag = 1 for view in self.viewlist: view.axis(1) - + # -------------------------------------------------------------------- def recenter(self): @@ -386,7 +386,7 @@ class vcr: self.orthoflag = 1 self.button_ortho.config(text="Persp") for view in self.viewlist: view.ortho(1) - + # -------------------------------------------------------------------- def reload(self): @@ -449,7 +449,7 @@ class vcr: # -------------------------------------------------------------------- # set filename for saving - + def file(self,newtext): oldtext = self.entry_file.get() self.entry_file.delete(0,len(oldtext)) @@ -457,7 +457,7 @@ class vcr: # -------------------------------------------------------------------- # toggle save all checkbox - + def saveall(self): if self.saveflag: self.saveflag = 0 @@ -465,11 +465,11 @@ class vcr: else: self.saveflag = 1 self.button_save.select() - + # -------------------------------------------------------------------- # save current image to file # if multiple windows change filenames to file.0.png, file.1.png, etc - + def save(self): file = self.entry_file.get() if len(self.viewlist) == 1: @@ -485,7 +485,7 @@ class vcr: # -------------------------------------------------------------------- # save images when in a play/back loop # flag 0 = first save, flag 1 = continuing save, flag -1 = stop - + def saveloop(self,flag): if flag == -1: self.file(self.fileroot) @@ -502,7 +502,7 @@ class vcr: # -------------------------------------------------------------------- # display index frame and set status strings - + def display(self): for view in self.viewlist: time,natoms = view.display(self.index) self.label_frame.config(text="Frame: %d" % self.index) diff --git a/src/vec.py b/src/vec.py index d9e1639..dfe8f04 100644 --- a/src/vec.py +++ b/src/vec.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. # vec tool @@ -17,12 +17,12 @@ v = vec(array) array = list of numeric vectors skip blank lines and lines that start with non-numeric characters example array with 2 vecs = [[1,2,3,4,5], [10,20,30,40,50]] assigns names = "col1", "col2", etc - + nvec = v.nvec # of vectors -nlen = v.nlen lengths of vectors -names = v.names list of vector names +nlen = v.nlen lengths of vectors +names = v.names list of vector names x,y,... = l.get(1,"col2",...) return one or more vectors of values -l.write("file.txt") write all vectors to a file +l.write("file.txt") write all vectors to a file l.write("file.txt","col1",7,...) write listed vectors to a file get and write allow abbreviated (uniquely) vector names or digits (1-Nvec) @@ -52,7 +52,7 @@ class vec: def __init__(self,data): self.data = [] - + if type(data) == types.StringType: lines = open(data,'r').readlines() for line in lines: @@ -69,7 +69,7 @@ class vec: self.data.append(map(float,values)) else: raise StandardError,"invalid argument to vec" - + if len(self.data) == 0: self.nlen = self.nvec = 0 else: @@ -85,7 +85,7 @@ class vec: self.ptr[self.names[i]] = i print "read %d vectors of length %d" % (self.nvec,self.nlen) - + # -------------------------------------------------------------------- def get(self,*keys): @@ -99,9 +99,9 @@ class vec: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: @@ -127,9 +127,9 @@ class vec: else: count = 0 for i in range(self.nvec): - if self.names[i].find(key) == 0: - count += 1 - index = i + if self.names[i].find(key) == 0: + count += 1 + index = i if count == 1: map.append(index) else: diff --git a/src/vizinfo.py b/src/vizinfo.py index acc421a..c5410d7 100644 --- a/src/vizinfo.py +++ b/src/vizinfo.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. # vizinfo class, not a top-level Pizza.py tool @@ -25,7 +25,7 @@ import types class vizinfo: """ Information holder for Pizza.py visualization tools - + acolor,bcolor,tcolor,lcolor = RGB values for each atom/bond/tri/line type arad = radius of each atom type brad,lrad = thickness of each bond/line type @@ -41,7 +41,7 @@ class vizinfo: setfill() = set triangle fill factor extend() = grow an array """ - + # -------------------------------------------------------------------- def __init__(self): @@ -57,15 +57,15 @@ class vizinfo: self.nbcolor = self.nbrad = 0 self.ntcolor = self.ntfill = 0 self.nlcolor = self.nlrad = 0 - + # -------------------------------------------------------------------- # set color RGB for which = atoms, bonds, triangles - + def setcolors(self,which,ids,rgbs): # convert args into lists if single values # if arg = 0, convert to full-range list - + if type(ids) is types.IntType and ids == 0: if which == "atom": ids = range(self.nacolor) if which == "bond": ids = range(self.nbcolor) @@ -101,11 +101,11 @@ class vizinfo: if max(ids) > self.nlcolor: self.nlcolor = self.extend(self.lcolor,max(ids)) self.nlcolor = self.extend(self.lrad,max(ids)) - + # set color for each type # if list lengths match, set directly, else interpolate # convert final color from 0-255 to 0.0-1.0 - + ntypes = len(ids) nrgbs = len(rgbs) @@ -135,7 +135,7 @@ class vizinfo: if which == "bond": self.bcolor[id] = color if which == "tri": self.tcolor[id] = color if which == "line": self.lcolor[id] = color - + # -------------------------------------------------------------------- # set radii for which = atoms, bonds, lines @@ -143,7 +143,7 @@ class vizinfo: # convert args into lists if single values # if arg = 0, convert to full-range list - + if type(ids) is types.IntType and ids == 0: if which == "atom": ids = range(self.narad) if which == "bond": ids = range(self.nbrad) @@ -199,16 +199,16 @@ class vizinfo: if which == "atom": self.arad[id] = rad if which == "bond": self.brad[id] = rad if which == "line": self.lrad[id] = rad - + # -------------------------------------------------------------------- # set triangle fill style # 0 = fill only, 1 = line only, 2 = fill and line - + def setfills(self,which,ids,fills): # convert args into lists if single values # if arg = 0, convert to full-range list - + if type(ids) is types.IntType and ids == 0: ids = range(self.ntfill) if type(ids) is not types.ListType and type(ids) is not types.TupleType: @@ -237,7 +237,7 @@ class vizinfo: for i in xrange(len(ids)): self.tfill[ids[i]] = int(fills[i]) else: for id in ids: self.tfill[id] = int(fills[0]) - + # -------------------------------------------------------------------- def extend(self,array,n): diff --git a/src/vmd.py b/src/vmd.py index 8b61ee2..1223451 100644 --- a/src/vmd.py +++ b/src/vmd.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. # vmd tool @@ -17,24 +17,24 @@ oneline = "Control VMD from python" docstr = """ -v = vmd() start up VMD -v.stop() shut down VMD instance -v.clear() delete all visualizations +v = vmd() start up VMD +v.stop() shut down VMD instance +v.clear() delete all visualizations -v.rep(style) set default representation style. One of - (Lines|VDW|Licorice|DynamicBonds|Points|CPK) -v.new(file[,type]) load new file (default file type 'lammpstrj') +v.rep(style) set default representation style. One of + (Lines|VDW|Licorice|DynamicBonds|Points|CPK) +v.new(file[,type]) load new file (default file type 'lammpstrj') v.data(file[,atomstyle]) load new data file (default atom style 'full') -v.replace(file[,type]) replace current frames with new file -v.append(file[,type]) append file to current frame(s) +v.replace(file[,type]) replace current frames with new file +v.append(file[,type]) append file to current frame(s) v.set(snap,x,y,z,(True|False)) set coordinates from a pizza.py snapshot to new or current frame -v.frame(frame) set current frame -v.flush() flush pending input to VMD and update GUI -v.read(file) read Tcl script file (e.g. saved state) - -v.enter() enter interactive shell -v.debug([True|False]) display generated VMD script commands? +v.frame(frame) set current frame +v.flush() flush pending input to VMD and update GUI +v.read(file) read Tcl script file (e.g. saved state) + +v.enter() enter interactive shell +v.debug([True|False]) display generated VMD script commands? """ # History @@ -55,7 +55,7 @@ try: from DEFAULTS import PIZZA_VMDARCH except: PIZZA_VMDARCH = "LINUX" try: import pexpect -except: +except: print "pexpect from http://pypi.python.org/pypi/pexpect", \ "is required for vmd tool" raise @@ -63,7 +63,7 @@ except: # Class definition class vmd: - + # -------------------------------------------------------------------- def __init__(self): @@ -95,7 +95,7 @@ class vmd: # open pipe to vmd and wait until we have a prompt self.VMD = pexpect.spawn(self.vmdexe) self.VMD.expect('vmd >') - + # -------------------------------------------------------------------- # post command to vmd and wait until the prompt returns. def __call__(self,command): @@ -105,7 +105,7 @@ class vmd: if self.debugme: print "call+result:"+self.VMD.before return - + # -------------------------------------------------------------------- # exit VMD def stop(self): @@ -190,7 +190,7 @@ class vmd: self.__call__('mol addfile ' + filename + ' mol $tmol type ' + filetype + ' waitfor all') self.__call__('foreach mol [molinfo list] { molinfo $mol set {center_matrix rotate_matrix scale_matrix global_matrix} $viewpoints($mol)}') self.flush() - + # -------------------------------------------------------------------- # replace all frames of a molecule with those from a given file def update(self,filename,filetype='lammpstrj'): @@ -201,7 +201,7 @@ class vmd: self.__call__('mol addfile ' + filename + ' mol $tmol type ' + filetype + ' waitfor all') self.__call__('foreach mol [molinfo list] {molinfo $mol set {center_matrix rotate_matrix scale_matrix global_matrix} $viewpoints($mol)}') self.flush() - + # -------------------------------------------------------------------- # add or overwrite coordinates with coordinates in a snapshot def set(self,snap,x,y,z,append=True): diff --git a/src/vtk.py b/src/vtk.py index 63e4f20..eb18050 100644 --- a/src/vtk.py +++ b/src/vtk.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. # vtk tool @@ -11,13 +11,13 @@ oneline = "Convert LAMMPS snapshots to VTK format" docstr = """ -v = vtk(d) d = object containing atom coords (dump, data) +v = vtk(d) d = object containing atom coords (dump, data) v.one() write all snapshots to tmp.vtk v.one("new") write all snapshots to new.vtk v.many() write snapshots to tmp0000.vtk, tmp0001.vtk, etc v.many("new") write snapshots to new0000.vtk, new0001.vtk, etc -v.single(N) write snapshot for timestep N to tmp.vtk +v.single(N) write snapshot for timestep N to tmp.vtk v.single(N,"file") write snapshot for timestep N to file.vtk surfaces in snapshot will be written to SURF1.vtk, SURF2.vtk, etc @@ -39,12 +39,12 @@ import sys, re # Class definition class vtk: - + # -------------------------------------------------------------------- def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def one(self,*args): @@ -59,16 +59,16 @@ class vtk: sys.stdout.flush() if len(tris): surface(tris) - + allatoms = [] for atom in atoms: allatoms.append(atom) - + while 1: which,time,flag = self.data.iterator(flag) if flag == -1: break time,box,atoms,bonds,tris,lines = self.data.viz(which) - + for atom in atoms: allatoms.append(atom) print time, sys.stdout.flush() @@ -93,7 +93,7 @@ class vtk: if surfflag == 0 and len(tris): surfflag = 1 surface(tris) - + if n < 10: file = root + "000" + str(n) + ".vtk" elif n < 100: @@ -104,19 +104,19 @@ class vtk: file = root + str(n) + ".vtk" particle(file,atoms) - + print time, sys.stdout.flush() n += 1 - + print "\nwrote %s snapshots in VTK format" % n # -------------------------------------------------------------------- def manyGran(self,*args,**kwargs): - + # check whether to output or not outputfl = True if "output" in kwargs: outputfl = kwargs["output"] - + # read startIndex (offset for filename due to parallel processing) startIndex = 0 fileNos = [] @@ -124,14 +124,14 @@ class vtk: fileNos = kwargs["fileNos"] else: fileNos = range(len(self.data.snaps)) - + # output name if len(args) == 0: root = "tmp" else: root = args[0] surfflag = 0 n = flag = 0 - + # iterate over snaps while 1: which,time,flag = self.data.iterator(flag) @@ -143,14 +143,14 @@ class vtk: yhi=self.data.snaps[n].yhi zlo=self.data.snaps[n].zlo zhi=self.data.snaps[n].zhi - + atoms=self.data.snaps[n].atoms names=self.data.names if surfflag == 0 and len(tris): surfflag = 1 surface(tris) - + file, file_bb, file_walls = generateFilename(root,fileNos,n) boundingBox(file_bb,xlo,xhi,ylo,yhi,zlo,zhi) @@ -158,11 +158,11 @@ class vtk: try: nvalues = len(self.data.snaps[0].atoms[0]) except: nvalues = 0 particleGran(file,atoms,names,nvalues) - + if outputfl: print time, if outputfl: sys.stdout.flush() n += 1 - + if outputfl: print "\nwrote %s granular snapshots in VTK format" % n # -------------------------------------------------------------------- @@ -209,7 +209,7 @@ def generateFilename(root,fileNos,n): def surface(tris): ntypes = tris[-1][1] - + for i in xrange(ntypes): itype = i+1 v = {} @@ -234,10 +234,10 @@ def surface(tris): vinverse = {} for key in keys: vinverse[v[key]] = key - + filename = "SURF" + str(itype) + ".vtk" f = open(filename,"w") - + print >>f,"# vtk DataFile Version 3.0" print >>f,"Generated by pizza.py" print >>f,"ASCII" @@ -267,7 +267,7 @@ def surface(tris): def particle(file,atoms): f = open(file,"w") - + print >>f,"# vtk DataFile Version 2.0" print >>f,"Generated by pizza.py" print >>f,"ASCII" @@ -285,13 +285,13 @@ def particle(file,atoms): itype = int(atom[1]) print >>f,itype, print >>f - + f.close() def boundingBox(file,xlo,xhi,ylo,yhi,zlo,zhi): f = open(file,"w") - + print >>f,"# vtk DataFile Version 2.0" print >>f,"Generated by pizza.py" print >>f,"ASCII" @@ -311,14 +311,14 @@ def typestr(o): def particleGran(file,atoms,names,n_values): f = open(file,"w") - + # if no atoms are present if atoms is None: atoms = [] - + # find indices of scalars and vectors scalars, vectors = findScalarsAndVectors(names) - + # print head print >>f,"# vtk DataFile Version 2.0" print >>f,"Generated by lpp.py" @@ -331,7 +331,7 @@ def particleGran(file,atoms,names,n_values): for i in xrange(len(atoms)): print >>f,1,i print >>f,"POINT_DATA",len(atoms) - + if len(atoms) == 0: print >> f f.close() @@ -339,11 +339,11 @@ def particleGran(file,atoms,names,n_values): # print VECTORS for key in vectors.keys(): - + # don't print coodinates again if key == 'x': continue - + vectortype = 'float' if atoms != []: vectortype = typestr(atoms[0][vectors[key]]) @@ -352,11 +352,11 @@ def particleGran(file,atoms,names,n_values): else: vectortype = 'float' else: # if no atoms are present pass - + print >>f,"VECTORS",key,vectortype for atom in atoms: print >>f, atom[vectors[key]], atom[vectors[key]+1], atom[vectors[key]+2] - + # print SCALARS for key in scalars.keys(): scalartype ='' @@ -367,30 +367,30 @@ def particleGran(file,atoms,names,n_values): else: scalartype = 'int' else: # if no atoms are present pass - + print >>f,"SCALARS",key,scalartype,1 print >>f,"LOOKUP_TABLE default" for atom in atoms: print >>f, atom[scalars[key]] - + print >>f f.close() def findScalarsAndVectors(names): - + vectors={} scalars={} - + # create reversed dictionary {position:name} indices = {} for name in names: indices[names[name]]=name - + # fill missing indices (occurrs e.g. if output is like vx vy vz fx fy fz vx vy vz) for i in xrange(max(indices)): if i not in indices: indices[i]="" - + # compile regexes to find vectors regvx = re.compile(".*x") regvy = re.compile(".*y") @@ -398,53 +398,53 @@ def findScalarsAndVectors(names): regf = re.compile("f_.*\[[0-9]+\]") regc = re.compile("c_.*\[[0-9]+\]") regv = re.compile("v_.*\[[0-9]+\]") - + # loop over all indices and look if their names represent a vector (if not: it's a scalar) i = 0 while i<= max(indices): - + if i+2 <= max(indices) and regvx.match(indices[i]) != None and regvy.match(indices[i+1]) != None and regvz.match(indices[i+2]) != None: - + newname='' if len(indices[i]) == 1: newname=indices[i] else: newname=indices[i][:-1] - + vectors[newname]=i i+=3 continue - + if regf.match(indices[i]) != None or regc.match(indices[i]) != None or regv.match(indices[i]) != None: - + name = indices[i] number = int( name.split('[')[1].split(']')[0] ) nextName = name.split('[')[0]+'['+str(number+1)+']' nextButOneName = name.split('[')[0]+'['+str(number+2)+']' - + newname = name[2:-(len(name.split('[')[1])+1)] - + if i+2 <= max(indices) and indices[i+1] == nextName and indices[i+2] == nextButOneName: vectors[newname]=i i+=3 continue - + else: scalars[newname]=i i+=1 continue - + # program only here if not a vector if indices[i] != '': newname = indices[i] scalars[newname]=i i+=1 - + if 'x' not in vectors.keys(): print "vector x y z has to be contained in dump file. please change liggghts input script accordingly." exit() - + return scalars, vectors - - - + + + diff --git a/src/xyz.py b/src/xyz.py index 66699ab..92b6815 100644 --- a/src/xyz.py +++ b/src/xyz.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. # xyz tool @@ -11,14 +11,14 @@ oneline = "Convert LAMMPS snapshots to XYZ format" docstr = """ -x = xyz(d) d = object containing atom coords (dump, data) +x = xyz(d) d = object containing atom coords (dump, data) x.one() write all snapshots to tmp.xyz x.one("new") write all snapshots to new.xyz x.many() write snapshots to tmp0000.xyz, tmp0001.xyz, etc x.many("new") write snapshots to new0000.xyz, new0001.xyz, etc -x.single(N) write snapshot for timestep N to tmp.xyz -x.single(N,"file") write snapshot for timestep N to file.xyz +x.single(N) write snapshot for timestep N to tmp.xyz +x.single(N,"file") write snapshot for timestep N to file.xyz """ # History @@ -41,7 +41,7 @@ class xyz: def __init__(self,data): self.data = data - + # -------------------------------------------------------------------- def one(self,*args): @@ -61,14 +61,14 @@ class xyz: for atom in atoms: itype = int(atom[1]) print >>f,itype,atom[2],atom[3],atom[4] - + print time, sys.stdout.flush() n += 1 - + f.close() print "\nwrote %d snapshots to %s in XYZ format" % (n,file) - + # -------------------------------------------------------------------- def many(self,*args): @@ -80,7 +80,7 @@ class xyz: which,time,flag = self.data.iterator(flag) if flag == -1: break time,box,atoms,bonds,tris,lines = self.data.viz(which) - + if n < 10: file = root + "000" + str(n) elif n < 100: @@ -88,7 +88,7 @@ class xyz: elif n < 1000: file = root + "0" + str(n) else: - file = root + str(n) + file = root + str(n) file += ".xyz" f = open(file,"w") print >>f,len(atoms) @@ -100,9 +100,9 @@ class xyz: sys.stdout.flush() f.close() n += 1 - + print "\nwrote %s snapshots in XYZ format" % n - + # -------------------------------------------------------------------- def single(self,time,*args):