--- /dev/null	2004-02-18 08:26:44.000000000 -0700
+++ /home/sjplimp/pizza/scripts/distance.py	2006-01-18 11:47:41.000000000 -0700
@@ -0,0 +1,75 @@
+#!/usr/bin/python
+
+# 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
+# Example: distance.py 0.95 dump.file1
+# Author:  Paul Crozier (Sandia)
+
+# print out 2 atoms less than maxcut apart (with PBC)
+
+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
+d = dump(files,0)
+d.map(1,"id",2,"type",3,"x",4,"y",5,"z")
+
+while 1:
+  time = d.next()
+  if time < 0: break
+  d.unscale(time)
+
+  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)
+  d.aselect.all(time)
+  id,type,x,y,z = d.vecs(time,"id","type","x","y","z")
+  n = len(x)
+
+  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:
+          delx += xprd
+        else:
+          delx -= xprd
+      if (delx*delx < maxcut_sq):
+          
+        dely = y[j] - y[i]
+        if abs(dely) > 0.5*yprd:
+          if dely < 0.0:
+            dely += yprd
+          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: 
+            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))
+
+  d.tselect.none()
+  d.tselect.one(time)
+  print "timestep = ", time
+  d.delete()
--- /dev/null	2004-02-18 08:26:44.000000000 -0700
+++ /home/sjplimp/pizza/scripts/bond_distribute.py	2006-01-18 11:37:58.000000000 -0700
@@ -0,0 +1,122 @@
+#!/usr/bin/python 
+
+# 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
+#          nbin = # of bins per bond type
+#          rmin = min expected bond length
+#          rmax = max expected bond length
+#          outfile = file to write stats to
+#          files = series of dump files
+# Example: bond_distribute.py pore.data 1000 1.5 1.85 bonds.out pore.dump.1
+# Author:  Paul Crozier (Sandia)
+
+# enable script to run from Python directly w/out Pizza.py
+
+import sys
+from dump import dump
+from math import sqrt
+if not globals().has_key("argv"): argv = sys.argv
+
+# main script
+
+if len(argv) < 7:
+  raise StandardError, \
+   "Syntax: bond_distribute.py datafile nbin rmin rmax outfile files ..."
+
+dt = data(argv[1])	
+nbins = int(argv[2])
+rmin = float(argv[3])
+rmax = float(argv[4])
+outfile = argv[5]
+files = ' '.join(argv[6:])
+
+# get the bonds from the data file
+
+bond = dt.get("Bonds")
+nbonds = len(bond)
+btype = nbonds * [0]
+iatom = nbonds * [0]
+jatom = nbonds * [0]
+for i in xrange(nbonds):
+  btype[i] = int(bond[i][1] - 1)
+  iatom[i] = int(bond[i][2] - 1)
+  jatom[i] = int(bond[i][3] - 1)
+
+ntypes = 0
+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] 
+
+# read snapshots one-at-a-time
+
+d = dump(files,0)
+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] 
+  zprd = box[5] - box[2]
+  
+  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
+      else:
+        delx -= xprd
+    if abs(dely) > 0.5*yprd:
+      if dely < 0.0:
+        dely += yprd
+      else:
+        dely -= yprd
+    if abs(delz) > 0.5*zprd:
+      if delz < 0.0:
+        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)): 
+      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
+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), 
+  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 
+fp.close()
--- /dev/null	2004-02-18 08:26:44.000000000 -0700
+++ /home/sjplimp/pizza/scripts/angle_distribute.py	2006-01-18 11:38:16.000000000 -0700
@@ -0,0 +1,156 @@
+#!/usr/bin/python 
+
+# Script:  angle_distribute.py 
+# Purpose: binned angle distributions by angle type
+# Syntax:  angle_distribute.py datafile nbin theta_min theta_max outfile files ...
+#          datafile = lammps data file
+#          nbin = # of bins per angle type
+#          theta_min = min expected angle
+#          theta_max = max expected angle length
+#          outfile = file to write stats to
+#          files = series of dump files
+# Example: angle_distribute.py pore.data 1000 110. 120. angles.out pore.dump.1
+# Author:  Paul Crozier (Sandia)
+
+# enable script to run from Python directly w/out Pizza.py
+
+import sys
+from dump import dump
+from math import sqrt,acos,atan
+if not globals().has_key("argv"): argv = sys.argv
+
+# main script
+
+if len(argv) < 7:
+  raise StandardError, \
+  "Syntax: angle_distribute.py datafile nbin theta_min theta_max outfile files ..."
+
+dt = data(argv[1])	
+nbins = int(argv[2])
+theta_min = float(argv[3])
+theta_max = float(argv[4])
+outfile = argv[5]
+files = ' '.join(argv[6:])
+
+# get the angles from the data file
+
+angle = dt.get("Angles")
+nangles = len(angle)
+atype = nangles * [0]
+iatom = nangles * [0]
+jatom = nangles * [0]
+katom = nangles * [0]
+for i in xrange(nangles):
+  atype[i] = int(angle[i][1] - 1)
+  iatom[i] = int(angle[i][2] - 1)
+  jatom[i] = int(angle[i][3] - 1)
+  katom[i] = int(angle[i][4] - 1)
+
+ntypes = 0
+for i in xrange(nangles): ntypes = max(angle[i][1],ntypes)
+ntypes = int(ntypes)
+ncount = ntypes * [0]
+bin = nbins * [0]
+for i in xrange(nbins): 
+  bin[i] = ntypes * [0] 
+
+# read snapshots one-at-a-time
+
+d = dump(files,0)
+d.map(1,"id",2,"type",3,"x",4,"y",5,"z")
+
+PI = 4.0*atan(1.0)
+
+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] 
+  zprd = box[5] - box[2]
+  
+  d.unscale() 
+  d.sort()
+  x,y,z = d.vecs(time,"x","y","z")
+   
+  for i in xrange(nangles):
+  
+    delx1 = x[iatom[i]] - x[jatom[i]]
+    dely1 = y[iatom[i]] - y[jatom[i]]
+    delz1 = z[iatom[i]] - z[jatom[i]]
+        
+    if abs(delx1) > 0.5*xprd:
+      if delx1 < 0.0:
+        delx1 += xprd
+      else:
+        delx1 -= xprd
+    if abs(dely1) > 0.5*yprd:
+      if dely1 < 0.0:
+        dely1 += yprd
+      else:
+        dely1 -= yprd
+    if abs(delz1) > 0.5*zprd:
+      if delz1 < 0.0:
+        delz1 += zprd
+      else:
+        delz1 -= zprd
+        
+    r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1)
+    
+    delx2 = x[katom[i]] - x[jatom[i]]
+    dely2 = y[katom[i]] - y[jatom[i]]
+    delz2 = z[katom[i]] - z[jatom[i]]
+        
+    if abs(delx2) > 0.5*xprd:
+      if delx2 < 0.0:
+        delx2 += xprd
+      else:
+        delx2 -= xprd
+    if abs(dely2) > 0.5*yprd:
+      if dely2 < 0.0:
+        dely2 += yprd
+      else:
+        dely2 -= yprd
+    if abs(delz2) > 0.5*zprd:
+      if delz2 < 0.0:
+        delz2 += zprd
+      else:
+        delz2 -= zprd
+        
+    r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2)
+    
+    c = delx1*delx2 + dely1*dely2 + delz1*delz2
+    c /= r1*r2
+        
+    if (c > 1.0): c = 1.0
+    if (c < -1.0): c = -1.0
+        
+    theta = 180.0*acos(c)/PI
+    
+    ibin = int(nbins*(theta - theta_min)/(theta_max - theta_min) + 0.5)
+    if ((ibin >= 0) and (ibin <= nbins-1)): 
+      bin[ibin][atype[i]] += nbins
+      ncount[atype[i]] += 1
+    else:
+      print "Warning: angle outside specified range"
+      print "angle type:", atype[i]+1
+      print "angle number:", i
+  print time,    
+      
+print
+print "Printing normalized angle distributions to",outfile
+    
+fp = open(outfile,"w")
+theta_range = theta_max - theta_min
+for i in xrange(nbins):
+  print >>fp, theta_min + theta_range*float(i)/float(nbins), 
+  for j in xrange(ntypes):
+    if (ncount[j] > 0):
+      print >>fp, float(bin[i][j])/float(ncount[j])/theta_range,
+    else:
+      print >>fp, 0.0,    
+  print >>fp 
+fp.close()