The python sub-directory contains several Python scripts +that perform common LAMMPS post-processing tasks, like +
+These are simple scripts built on Pizza.py modules. See the +README for more info on Pizza.py and how to use these scripts. +
+The file restart2data.cpp converts a binary LAMMPS restart file into diff -Naur lammps-1Apr08/doc/Section_tools.txt lammps-2Apr08/doc/Section_tools.txt --- lammps-1Apr08/doc/Section_tools.txt 2008-03-28 09:03:22.000000000 -0600 +++ lammps-2Apr08/doc/Section_tools.txt 2008-03-31 16:44:23.000000000 -0600 @@ -53,6 +53,7 @@ "micelle2d"_#micelle "msi2lmp"_#msi "pymol_asphere"_#pymol +"python"_#pythontools "restart2data"_#restart "thermo_extract"_#thermo_extract "vim"_#vim @@ -270,6 +271,21 @@ :line +python tool :h4,link(pythontools) + +The python sub-directory contains several Python scripts +that perform common LAMMPS post-processing tasks, like + +extract thermodynamic info from a log file as columns of numbers +plot two columns of thermodynamic info from a log file +sort the snapshots in a dump file by atom ID +convert dump files into XYZ, CFG, or PDB format for viz by other packages :ul + +These are simple scripts built on "Pizza.py"_pizza modules. See the +README for more info on Pizza.py and how to use these scripts. + +:line + restart2data tool :h4,link(restart) The file restart2data.cpp converts a binary LAMMPS restart file into diff -Naur lammps-1Apr08/src/MAKE/Makefile.mac_mpi lammps-2Apr08/src/MAKE/Makefile.mac_mpi --- lammps-1Apr08/src/MAKE/Makefile.mac_mpi 2007-10-22 15:05:37.000000000 -0600 +++ lammps-2Apr08/src/MAKE/Makefile.mac_mpi 2008-03-28 11:50:59.000000000 -0600 @@ -7,7 +7,7 @@ FFTW = /sw CC = mpic++ -CCFLAGS = -O -MMD -DFFT_FFTW -I${FFTW}/include -DOMPI_SKIP_MPICXX +CCFLAGS = -O -MMD -MG -DFFT_FFTW -I${FFTW}/include -DOMPI_SKIP_MPICXX LINK = mpic++ LINKFLAGS = -O -L${FFTW}/lib diff -Naur lammps-1Apr08/tools/python/README lammps-2Apr08/tools/python/README --- lammps-1Apr08/tools/python/README 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/README 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,111 @@ +This directory contains Python scripts for common LAMMPS +post-processing tasks. + +If you have suggestions or contributions for additional scripts or +functionality that could be added, built on the Pizza.py modules (as +explained below), send email to Steve Plimpton (sjplimp at +sandia.gov). + +log2txt.py convert thermo info in a LAMMPS log file to columns of numbers +logplot.py plot 2 columns of thermo info from a LAMMPS log file +dumpsort.py sort the snapshots of atoms in a LAMMPS dump file by atom ID +dump2cfg.py convert a native LAMMPS dump file to CFG format +dump2xyz.py convert a native LAMMPS dump file to XYZ format +dump2pdb.py convert a native LAMMPS dump file to PDB format + +See the top of each script file for syntax, or just run it with no +arguments to get a syntax message. + +------------------------------ +General info: + +These scripts are very simple. They load Python modules in the pizza +sub-directory that are part of the Pizza.py toolkit, which do the +heavy lifting. + +The modules themselves have a lot more functionality than these +scripts expose, so if you want to do something beyond what these +scripts perform, you should learn about Pizza.py. See this WWW page +for details and download info: + +www.sandia.gov/~sjplimp/pizza.html + +The tools in the Pizza.py src directory are identical to those in the +pizza sub-directory of lammps/tools/python. The header section of the +tool files lists all the functionality that tool supports. + +To use all the features of Pizza.py modules, you need to be familiar +with Python syntax. You can then modify the scripts to invoke +additional Pizza.py functionality or use the Python interpreter itself +to drive the Pizza.py modules. + +------------------------------ +Before you run the scripts: + +To use these scripts you must set the environment variable +LAMMPS_TOOLS_PYTHON in your shell to point to the Pizza.py modules +that the scripts use. This can either be a) the pizza sub-directory +under lammps/tools/python, or b) the src directory in the Pizza.py +package if you have installed Pizza.py on your box. + +For example, on my box, either of these lines in my .cshrc works: + +setenv LAMMPS_TOOLS_PYTHON /home/sjplimp/lammps/tools/python/pizza +setenv LAMMPS_TOOLS_PYTHON /home/sjplimp/pizza/src + +------------------------------ +Running the scripts: + +As with any Python script, you can run these scripts in one of two +ways. You may want to setup aliases so that you can run them from the +directory where your data files are. + +% python log2txt.py args ... +% log2txt.py args ... + +The latter requires 2 things: + +1) that the script be made "executable", e.g. type "chmod +x log2txt.py" +2) that the 1st line of the script is the path of the Python installed + on your box, e.g. /usr/local/bin/python + +IMPORTANT NOTE: If you run the logplot.py script using the 1st method +above, you should type + +% python -i logplot.py args ... + +so that the plot produced by GnuPlot stays visible before Python +exits. + +------------------------------ +Dependencies: + +To use the logplot.py script you need to have GnuPlot installed on +your system and its executable "gnuplot" in your path. + +To use any of the scripts which load the dump module to read LAMMPS +dump files, you must have the Python package Numeric installed in your +Python. See http://numeric.scipy.org. + +Note that the Pizza.py modules use the older (but still popular) +Numeric package, not the newer numarray package. + +If Numeric is already installed in your Python, you should be able to +type the following without getting an error: + +>>> import Numeric + +Numeric can be downloaded from SourceForge at this WWW site: + +http://sourceforge.net/project/showfiles.php?group_id=1369&package_id=1351 + +As of July 2007, Version 24.2 was fine. Once unpacked, you can type +the following from the Numeric directory to install it in your Python. + +sudo python setup.py install + +On my Linux box, when Numeric installed itself under the Python lib in +/usr/local, it did not set all file permsissions correctly to allow a +user to import it. So I also needed to do this: + +sudo chmod -R 755 /usr/local/lib/python2.5/site-packages/Numeric diff -Naur lammps-1Apr08/tools/python/dump2cfg.py lammps-2Apr08/tools/python/dump2cfg.py --- lammps-1Apr08/tools/python/dump2cfg.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/dump2cfg.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,32 @@ +#!/usr/local/bin/python + +# Script: dump2cfg.py +# Purpose: convert a LAMMPS dump file to CFG format +# Syntax: dump2cfg.py dumpfile Nid Ntype Nx Ny Nz cfgfile +# dumpfile = LAMMPS dump file in native LAMMPS format +# Nid,Ntype,Nx,Ny,Nz = columns #s for ID,type,x,y,z +# (usually 1,2,3,4,5) +# cfgfile = new CFG file +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from dump import dump +from cfg import cfg + +if len(sys.argv) != 8: + raise StandardError, "Syntax: dump2cfg.py dumpfile Nid Ntype Nx Ny Nz cfgfile" + +dumpfile = sys.argv[1] +nid = int(sys.argv[2]) +ntype = int(sys.argv[3]) +nx = int(sys.argv[4]) +ny = int(sys.argv[5]) +nz = int(sys.argv[6]) +cfgfile = sys.argv[7] + +d = dump(dumpfile) +d.map(nid,"id",ntype,"type",nx,"x",ny,"y",nz,"z") +c = cfg(d) +c.one(cfgfile) diff -Naur lammps-1Apr08/tools/python/dump2pdb.py lammps-2Apr08/tools/python/dump2pdb.py --- lammps-1Apr08/tools/python/dump2pdb.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/dump2pdb.py 2008-03-31 16:39:55.000000000 -0600 @@ -0,0 +1,37 @@ +#!/usr/local/bin/python + +# Script: dump2pdb.py +# Purpose: convert a LAMMPS dump file to PDB format +# Syntax: dump2pdb.py dumpfile Nid Ntype Nx Ny Nz pdbfile template +# dumpfile = LAMMPS dump file in native LAMMPS format +# Nid,Ntype,Nx,Ny,Nz = columns #s for ID,type,x,y,z +# (usually 1,2,3,4,5) +# pdbfile = new PDB file +# template = PDB file to use as template for creating new PDB file +# this arg is optional, if not used a generic PDB file is created +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from dump import dump +from pdbfile import pdbfile + +if len(sys.argv) != 8 and len(sys.argv) != 9: + raise StandardError, "Syntax: dump2pdb.py dumpfile Nid Ntype Nx Ny Nz pdbfile template" + +dumpfile = sys.argv[1] +nid = int(sys.argv[2]) +ntype = int(sys.argv[3]) +nx = int(sys.argv[4]) +ny = int(sys.argv[5]) +nz = int(sys.argv[6]) +pfile = sys.argv[7] +if len(sys.argv) == 9: template = sys.argv[8] +else: template = "" + +d = dump(dumpfile) +d.map(nid,"id",ntype,"type",nx,"x",ny,"y",nz,"z") +if template: p = pdbfile(template,d) +else: p = pdbfile(d) +p.one(pfile) diff -Naur lammps-1Apr08/tools/python/dump2xyz.py lammps-2Apr08/tools/python/dump2xyz.py --- lammps-1Apr08/tools/python/dump2xyz.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/dump2xyz.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,32 @@ +#!/usr/local/bin/python + +# Script: dump2xyz.py +# Purpose: convert a LAMMPS dump file to XYZ format +# Syntax: dump2xyz.py dumpfile Nid Ntype Nx Ny Nz xyzfile +# dumpfile = LAMMPS dump file in native LAMMPS format +# Nid,Ntype,Nx,Ny,Nz = columns #s for ID,type,x,y,z +# (usually 1,2,3,4,5) +# xyzfile = new XYZ file +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from dump import dump +from xyz import xyz + +if len(sys.argv) != 8: + raise StandardError, "Syntax: dump2xyz.py dumpfile Nid Ntype Nx Ny Nz xyzfile" + +dumpfile = sys.argv[1] +nid = int(sys.argv[2]) +ntype = int(sys.argv[3]) +nx = int(sys.argv[4]) +ny = int(sys.argv[5]) +nz = int(sys.argv[6]) +xyzfile = sys.argv[7] + +d = dump(dumpfile) +d.map(nid,"id",ntype,"type",nx,"x",ny,"y",nz,"z") +x = xyz(d) +x.one(xyzfile) diff -Naur lammps-1Apr08/tools/python/dumpsort.py lammps-2Apr08/tools/python/dumpsort.py --- lammps-1Apr08/tools/python/dumpsort.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/dumpsort.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,26 @@ +#!/usr/local/bin/python + +# Script: dumpsort.py +# Purpose: sort the snapshots in a LAMMPS dump file by atom ID +# Syntax: dumpsort.py oldfile N newfile +# oldfile = old LAMMPS dump file in native LAMMPS format +# N = column # for atom ID (usually 1) +# newfile = new sorted LAMMPS dump file +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from dump import dump + +if len(sys.argv) != 4: + raise StandardError, "Syntax: dumpsort.py oldfile N newfile" + +oldfile = sys.argv[1] +ncolumn = int(sys.argv[2]) +newfile = sys.argv[3] + +d = dump(oldfile) +d.map(ncolumn,"id") +d.sort() +d.write(newfile) diff -Naur lammps-1Apr08/tools/python/log2txt.py lammps-2Apr08/tools/python/log2txt.py --- lammps-1Apr08/tools/python/log2txt.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/log2txt.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,32 @@ +#!/usr/local/bin/python + +# Script: log2txt.py +# Purpose: extract thermo info from LAMMPS log file +# create a text file of numbers in columns, suitable for plotting +# Syntax: log2txt.py log.lammps data.txt X Y ... +# log.lammps = LAMMPS log file +# data.txt = text file to create +# X Y ... = columns to include (optional), X,Y are thermo keywords +# if no columns listed, all columns are included +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from log import log + +if len(sys.argv) < 3: + raise StandardError, "Syntax: log2txt.py log.lammps data.txt X Y ..." + +logfile = sys.argv[1] +datafile = sys.argv[2] +columns = sys.argv[3:] + +lg = log(logfile) +if columns == []: + lg.write(datafile) +else: + str = "lg.write(datafile," + for word in columns: str += '"' + word + '",' + str = str[:-1] + ')' + eval(str) diff -Naur lammps-1Apr08/tools/python/logplot.py lammps-2Apr08/tools/python/logplot.py --- lammps-1Apr08/tools/python/logplot.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/logplot.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,28 @@ +#!/usr/local/bin/python -i + +# Script: logplot.py +# Purpose: use GnuPlot to plot two columns from a LAMMPS log file +# Syntax: logplot.py log.lammps X Y +# log.lammps = LAMMPS log file +# X,Y = plot Y versus X where X,Y are thermo keywords +# once plot appears, you are in Python interpreter, type C-D to exit +# Author: Steve Plimpton (Sandia), sjplimp at sandia.gov + +import sys,os +path = os.environ["LAMMPS_PYTHON_TOOLS"] +sys.path.append(path) +from log import log +from gnu import gnu + +if len(sys.argv) != 4: + raise StandardError, "Syntax: logplot.py log.lammps X Y" + +logfile = sys.argv[1] +xlabel = sys.argv[2] +ylabel = sys.argv[3] + +lg = log(logfile) +x,y = lg.get(xlabel,ylabel) +g = gnu() +g.plot(x,y) +print "Type Ctrl-D to exit Python" diff -Naur lammps-1Apr08/tools/python/pizza/cfg.py lammps-2Apr08/tools/python/pizza/cfg.py --- lammps-1Apr08/tools/python/pizza/cfg.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/cfg.py 2008-03-31 16:30:20.000000000 -0600 @@ -0,0 +1,187 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# cfg tool + +oneline = "Convert LAMMPS snapshots to AtomEye CFG format" + +docstr = """ +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 +""" + +# History +# 11/06, Aidan Thompson (SNL): original version + +# ToDo list +# should decide if dump is scaled or not, since CFG prints in scaled coords +# this creates a simple AtomEye CFG format +# there is more complex format we could write out +# which allows for extra atom info, e.g. to do atom coloring on +# how to dump for a triclinic box, since AtomEye accepts this + +# Variables +# data = data file to read from + +# Imports and external programs + +import sys + +# Class definition + +class cfg: + + # -------------------------------------------------------------------- + + def __init__(self,data): + self.data = data + + # -------------------------------------------------------------------- + + def one(self,*args): + if len(args) == 0: file = "tmp.cfg" + elif args[0][-4:] == ".cfg": file = args[0] + else: file = args[0] + ".cfg" + + f = open(file,"w") + n = flag = 0 + while 1: + which,time,flag = self.data.iterator(flag) + if flag == -1: break + time,box,atoms,bonds,tris,lines = self.data.viz(which) + + xlen = box[3]-box[0] + ylen = box[4]-box[1] + zlen = box[5]-box[2] + + 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(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(3,3) = %20.10f A " % zlen + print >>f,"#" + + for atom in atoms: + itype = int(atom[1]) + xfrac = (atom[2]-box[0])/xlen + yfrac = (atom[3]-box[1])/ylen + 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): + if len(args) == 0: root = "tmp" + else: root = args[0] + + n = flag = 0 + while 1: + 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: + file = root + "00" + str(n) + elif n < 1000: + file = root + "0" + str(n) + else: + file = root + str(n) + file += ".cfg" + f = open(file,"w") + + xlen = box[3]-box[0] + ylen = box[4]-box[1] + zlen = box[5]-box[2] + + 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(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(3,3) = %20.10f A " % zlen + print >>f,"#" + + for atom in atoms: + itype = int(atom[1]) + xfrac = (atom[2]-box[0])/xlen + yfrac = (atom[3]-box[1])/ylen + 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): + if len(args) == 0: file = "tmp.cfg" + elif args[0][-4:] == ".cfg": file = args[0] + else: file = args[0] + ".cfg" + + which = self.data.findtime(time) + time,box,atoms,bonds,tris,lines = self.data.viz(which) + f = open(file,"w") + + xlen = box[3]-box[0] + ylen = box[4]-box[1] + zlen = box[5]-box[2] + + 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(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(3,3) = %20.10f A " % zlen + print >>f,"#" + + for atom in atoms: + itype = int(atom[1]) + xfrac = (atom[2]-box[0])/xlen + yfrac = (atom[3]-box[1])/ylen + 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 -Naur lammps-1Apr08/tools/python/pizza/dump.py lammps-2Apr08/tools/python/pizza/dump.py --- lammps-1Apr08/tools/python/pizza/dump.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/dump.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,1155 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# dump tool + +oneline = "Read, write, manipulate dump files and particle attributes" + +docstr = """ +d = dump("dump.one") read in one or more dump files +d = dump("dump.1 dump.2.gz") can be gzipped +d = dump("dump.*") wildcard expands to multiple files +d = dump("dump.*",0) two args = store filenames, but don't read + + incomplete and duplicate snapshots are deleted + if atoms have 5 or 8 columns, assign id,type,x,y,z (ix,iy,iz) + atoms will be unscaled if stored in files as scaled + +time = d.next() read next snapshot from dump files + + used with 2-argument constructor to allow reading snapshots one-at-a-time + snapshot will be skipped only if another snapshot has same time stamp + return time stamp of snapshot read + return -1 if no snapshots left or last snapshot is incomplete + no column name assignment or unscaling is performed + +d.map(1,"id",3,"x") assign names to atom columns (1-N) + +d.tselect.all() select all timesteps +d.tselect.one(N) select only timestep N +d.tselect.none() deselect all timesteps +d.tselect.skip(M) select every Mth step +d.tselect.test("$t >= 100 and $t < 10000") select matching timesteps +d.delete() delete non-selected timesteps + + selecting a timestep also selects all atoms in the timestep + skip() and test() only select from currently selected timesteps + test() uses a Python Boolean expression with $t for timestep value + Python comparison syntax: == != < > <= >= and or + +d.aselect.all() select all atoms in all steps +d.aselect.all(N) select all atoms in one step +d.aselect.test("$id > 100 and $type == 2") select match atoms in all steps +d.aselect.test("$id > 100 and $type == 2",N) select matching atoms in one step + + all() with no args selects atoms from currently selected timesteps + test() with one arg selects atoms from currently selected timesteps + test() sub-selects from currently selected atoms + test() uses a Python Boolean expression with $ for atom attributes + Python comparison syntax: == != < > <= >= and or + $name must end with a space + +d.write("file") write selected steps/atoms to dump file +d.scatter("tmp") write selected steps/atoms to mutiple files + + scatter() files are given timestep suffix: e.g. tmp.0, tmp.100, etc + +d.scale() scale x,y,z to 0-1 for all timesteps +d.scale(100) scale atom coords for timestep N +d.unscale() unscale x,y,z to box size to all timesteps +d.unscale(1000) unscale atom coords for timestep N +d.wrap() wrap x,y,z into periodic box via ix,iy,iz +d.unwrap() unwrap x,y,z out of box via ix,iy,iz +d.owrap("other") wrap x,y,z to same image as another atom +d.sort() sort atoms by atom ID in all selected steps +d.sort("x") sort atoms by column value in all steps +d.sort(1000) sort atoms in timestep N + + scale(), unscale(), wrap(), unwrap(), owrap() operate on all steps and atoms + wrap(), unwrap(), owrap() require ix,iy,iz be defined + owrap() requires a column be defined which contains an atom ID + name of that column is the argument to owrap() + x,y,z for each atom is wrapped to same image as the associated atom ID + useful for wrapping all molecule's atoms the same so it is contiguous + +m1,m2 = d.minmax("type") find min/max values for a column +d.set("$ke = $vx * $vx + $vy * $vy") set a column to a computed value +d.spread("ke",N,"color") 2nd col = N ints spread over 1st col +d.clone(1000,"color") clone timestep N values to other steps + + minmax() operates on selected timesteps and atoms + set() operates on selected timesteps and atoms + left hand side column is created if necessary + left-hand side column is unset or unchanged for non-selected atoms + equation is in Python syntax + use $ for column names, $name must end with a space + spread() operates on selected timesteps and atoms + min and max are found for 1st specified column across all selected atoms + atom's value is linear mapping (1-N) between min and max + that is stored in 2nd column (created if needed) + useful for creating a color map + clone() operates on selected timesteps and atoms + values at every timestep are set to value at timestep N for that atom ID + useful for propagating a color map + +t = d.time() return vector of selected timestep values +fx,fy,... = d.atom(100,"fx","fy",...) return vector(s) for atom ID N +fx,fy,... = d.vecs(1000,"fx","fy",...) return vector(s) for timestep N + + atom() returns vectors with one value for each selected timestep + vecs() returns vectors with one value for each selected atom in the timestep + +index,time,flag = d.iterator(0/1) loop over dump snapshots +time,box,atoms,bonds,tris = d.viz(index) return list of viz objects +d.atype = "color" set column returned as "type" by viz +d.extra("dump.bond") read bond list from dump file +d.extra(data) extract bond/tri/line list from data + + iterator() loops over selected timesteps + iterator() called with arg = 0 first time, with arg = 1 on subsequent calls + index = index within dump object (0 to # of snapshots) + time = timestep value + flag = -1 when iteration is done, 1 otherwise + viz() returns info for selected atoms for specified timestep index + time = timestep value + box = [xlo,ylo,zlo,xhi,yhi,zhi] + atoms = id,type,x,y,z for each atom as 2d array + bonds = id,type,x1,y1,z1,x2,y2,z2,t1,t2 for each bond as 2d array + if bonds() was used to define bonds, else empty list + tris = id,type,x1,y1,z1,x2,y2,z2,x3,y3,z3,nx,ny,nz for each tri as 2d array + if extra() was used to define tris, else empty list + lines = id,type,x1,y1,z1,x2,y2,z2 for each line as 2d array + if extra() was used to define lines, else empty list + atype is column name viz() will return as atom type (def = "type") + extra() stores list of bonds/tris/lines to return each time viz() is called +""" + +# History +# 8/05, Steve Plimpton (SNL): original version + +# ToDo list +# try to optimize this line in read_snap: words += f.readline().split() +# allow $name in aselect.test() and set() to end with non-space +# should next() snapshot be auto-unscaled ? + +# Variables +# flist = list of dump file names +# increment = 1 if reading snapshots one-at-a-time +# nextfile = which file to read from via next() +# eof = ptr into current file for where to read via next() +# nsnaps = # of snapshots +# nselect = # of selected snapshots +# snaps = list of snapshots +# names = dictionary of column names: +# key = "id", value = column # (0 to M-1) +# tselect = class for time selection +# aselect = class for atom selection +# atype = name of vector used as atom type by viz extract +# bondflag = 0 if no bonds, 1 if they are defined statically +# bondlist = static list of bonds to viz() return for all snapshots +# only a list of atom pairs, coords have to be created for each snapshot +# triflag = 0 if no tris, 1 if they are defined statically, 2 if dynamic +# trilist = static list of tris to return via viz() for all snapshots +# lineflag = 0 if no lines, 1 if they are defined statically +# linelist = static list of lines to return via viz() for all snapshots +# Snap = one snapshot +# time = time stamp +# tselect = 0/1 if this snapshot selected +# natoms = # of atoms +# nselect = # of selected atoms in this snapshot +# aselect[i] = 0/1 for each atom +# xlo,xhi,ylo,yhi,zlo,zhi = box bounds (float) +# atoms[i][j] = 2d array of floats, i = 0 to natoms-1, j = 0 to ncols-1 + +# Imports and external programs + +import sys, commands, re, glob, types +from os import popen +from math import * # any function could be used by set() +import Numeric + +try: from DEFAULTS import PIZZA_GUNZIP +except: PIZZA_GUNZIP = "gunzip" + +# Class definition + +class dump: + + # -------------------------------------------------------------------- + + def __init__(self,*list): + self.snaps = [] + self.nsnaps = self.nselect = 0 + self.names = {} + self.tselect = tselect(self) + self.aselect = aselect(self) + self.atype = "type" + self.bondflag = 0 + self.bondlist = [] + self.triflag = 0 + self.trilist = [] + self.triobj = 0 + self.lineflag = 0 + self.linelist = [] + + # flist = list of all dump file names + + words = list[0].split() + self.flist = [] + for word in words: self.flist += glob.glob(word) + if len(self.flist) == 0 and len(list) == 1: + raise StandardError,"no dump file specified" + + if len(list) == 1: + self.increment = 0 + self.read_all() + else: + self.increment = 1 + self.nextfile = 0 + self.eof = 0 + + # -------------------------------------------------------------------- + + def read_all(self): + + # read all snapshots from each file + # test for gzipped files + + for file in self.flist: + if file[-3:] == ".gz": + f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r') + else: f = open(file) + + snap = self.read_snapshot(f) + while snap: + self.snaps.append(snap) + print snap.time, + sys.stdout.flush() + snap = self.read_snapshot(f) + + f.close() + print + + # sort entries by timestep, cull duplicates + + self.snaps.sort(self.compare_time) + self.cull() + self.nsnaps = len(self.snaps) + print "read %d snapshots" % self.nsnaps + + # select all timesteps and atoms + + self.tselect.all() + + # set default names for atom columns + + if len(self.snaps) == 0: + print "no column assignments made" + elif self.snaps[0].atoms == None: + print "no column assignments made" + elif len(self.snaps[0].atoms[0]) == 5: + self.map(1,"id",2,"type",3,"x",4,"y",5,"z") + print "assigned columns id,type,x,y,z" + elif len(self.snaps[0].atoms[0]) == 8: + self.map(1,"id",2,"type",3,"x",4,"y",5,"z",6,"ix",7,"iy",8,"iz") + print "assigned columns id,type,x,y,z,ix,iy,iz" + else: + print "no column assignments made" + + # if snapshots are scaled, unscale them + + if not self.names.has_key("x"): + print "no unscaling could be performed" + elif self.nsnaps > 0: + if self.scaled(self.nsnaps-1): self.unscale() + else: print "dump is already unscaled" + + # -------------------------------------------------------------------- + # read next snapshot from list of files + + def next(self): + + if not self.increment: raise StandardError,"cannot read incrementally" + + # read next snapshot in current file using eof as pointer + # if fail, try next file + # if new snapshot time stamp already exists, read next snapshot + + while 1: + f = open(self.flist[self.nextfile],'r') + f.seek(self.eof) + snap = self.read_snapshot(f) + if not snap: + self.nextfile += 1 + if self.nextfile == len(self.flist): return -1 + f.close() + self.eof = 0 + continue + self.eof = f.tell() + f.close() + try: + self.findtime(snap.time) + continue + except: break + + # select the new snapshot with all its atoms + + self.snaps.append(snap) + snap = self.snaps[self.nsnaps] + snap.tselect = 1 + snap.nselect = snap.natoms + for i in xrange(snap.natoms): snap.aselect[i] = 1 + self.nsnaps += 1 + self.nselect += 1 + + return snap.time + + # -------------------------------------------------------------------- + # read a single snapshot from file f + # return snapshot or 0 if failed + + def read_snapshot(self,f): + try: + snap = Snap() + item = f.readline() + snap.time = int(f.readline()) + item = f.readline() + snap.natoms = int(f.readline()) + + snap.aselect = Numeric.zeros(snap.natoms) + + item = f.readline() + words = f.readline().split() + snap.xlo,snap.xhi = float(words[0]),float(words[1]) + words = f.readline().split() + snap.ylo,snap.yhi = float(words[0]),float(words[1]) + words = f.readline().split() + snap.zlo,snap.zhi = float(words[0]),float(words[1]) + + item = f.readline() + if snap.natoms: + words = f.readline().split() + ncol = len(words) + for i in xrange(1,snap.natoms): + words += f.readline().split() + floats = map(float,words) + atoms = Numeric.zeros((snap.natoms,ncol),Numeric.Float) + start = 0 + stop = ncol + for i in xrange(snap.natoms): + atoms[i] = floats[start:stop] + start = stop + stop += ncol + else: atoms = None + snap.atoms = atoms + return snap + except: + return 0 + + # -------------------------------------------------------------------- + # decide if snapshot i is scaled/unscaled from coords of first and last atom + + def scaled(self,i): + ix = self.names["x"] + iy = self.names["y"] + iz = self.names["z"] + natoms = self.snaps[i].natoms + if natoms == 0: return 0 + x1 = self.snaps[i].atoms[0][ix] + y1 = self.snaps[i].atoms[0][iy] + z1 = self.snaps[i].atoms[0][iz] + x2 = self.snaps[i].atoms[natoms-1][ix] + y2 = self.snaps[i].atoms[natoms-1][iy] + z2 = self.snaps[i].atoms[natoms-1][iz] + if x1 >= -0.1 and x1 <= 1.1 and y1 >= -0.1 and y1 <= 1.1 and \ + z1 >= -0.1 and z1 <= 1.1 and x2 >= -0.1 and x2 <= 1.1 and \ + y2 >= -0.1 and y2 <= 1.1 and z2 >= -0.1 and z2 <= 1.1: + return 1 + else: return 0 + + # -------------------------------------------------------------------- + # map atom column names + + def map(self,*pairs): + if len(pairs) % 2 != 0: + raise StandardError, "dump map() requires pairs of mappings" + for i in range(0,len(pairs),2): + j = i + 1 + self.names[pairs[j]] = pairs[i]-1 + + # delete unselected snapshots + + # -------------------------------------------------------------------- + + def delete(self): + ndel = i = 0 + while i < self.nsnaps: + if not self.snaps[i].tselect: + del self.snaps[i] + self.nsnaps -= 1 + ndel += 1 + else: i += 1 + print "%d snapshots deleted" % ndel + print "%d snapshots remaining" % self.nsnaps + + # -------------------------------------------------------------------- + # scale coords to 0-1 for all snapshots or just one + + def scale(self,*list): + if len(list) == 0: + print "Scaling dump ..." + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + for snap in self.snaps: self.scale_one(snap,x,y,z) + else: + i = self.findtime(list[0]) + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + self.scale_one(self.snaps[i],x,y,z) + + # -------------------------------------------------------------------- + + def scale_one(self,snap,x,y,z): + xprdinv = 1.0 / (snap.xhi - snap.xlo) + yprdinv = 1.0 / (snap.yhi - snap.ylo) + zprdinv = 1.0 / (snap.zhi - snap.zlo) + atoms = snap.atoms + atoms[:,x] = (atoms[:,x] - snap.xlo) * xprdinv + atoms[:,y] = (atoms[:,y] - snap.ylo) * yprdinv + atoms[:,z] = (atoms[:,z] - snap.zlo) * zprdinv + + # -------------------------------------------------------------------- + # unscale coords from 0-1 to box size for all snapshots or just one + + def unscale(self,*list): + if len(list) == 0: + print "Unscaling dump ..." + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + for snap in self.snaps: self.unscale_one(snap,x,y,z) + else: + i = self.findtime(list[0]) + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + self.unscale_one(self.snaps[i],x,y,z) + + # -------------------------------------------------------------------- + + def unscale_one(self,snap,x,y,z): + xprd = snap.xhi - snap.xlo + yprd = snap.yhi - snap.ylo + zprd = snap.zhi - snap.zlo + atoms = snap.atoms + atoms[:,x] = snap.xlo + atoms[:,x]*xprd + atoms[:,y] = snap.ylo + atoms[:,y]*yprd + atoms[:,z] = snap.zlo + atoms[:,z]*zprd + + # -------------------------------------------------------------------- + # wrap coords from outside box to inside + + def wrap(self): + print "Wrapping dump ..." + + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + ix = self.names["ix"] + iy = self.names["iy"] + iz = self.names["iz"] + + for snap in self.snaps: + xprd = snap.xhi - snap.xlo + yprd = snap.yhi - snap.ylo + zprd = snap.zhi - snap.zlo + atoms = snap.atoms + atoms[:,x] -= atoms[:,ix]*xprd + atoms[:,y] -= atoms[:,iy]*yprd + atoms[:,z] -= atoms[:,iz]*zprd + + # -------------------------------------------------------------------- + # unwrap coords from inside box to outside + + def unwrap(self): + print "Unwrapping dump ..." + + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + ix = self.names["ix"] + iy = self.names["iy"] + iz = self.names["iz"] + + for snap in self.snaps: + xprd = snap.xhi - snap.xlo + yprd = snap.yhi - snap.ylo + zprd = snap.zhi - snap.zlo + atoms = snap.atoms + atoms[:,x] += atoms[:,ix]*xprd + atoms[:,y] += atoms[:,iy]*yprd + atoms[:,z] += atoms[:,iz]*zprd + + # -------------------------------------------------------------------- + # wrap coords to same image as atom ID stored in "other" column + + def owrap(self,other): + print "Wrapping to other ..." + + id = self.names["id"] + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + ix = self.names["ix"] + iy = self.names["iy"] + iz = self.names["iz"] + iother = self.names[other] + + for snap in self.snaps: + xprd = snap.xhi - snap.xlo + yprd = snap.yhi - snap.ylo + zprd = snap.zhi - snap.zlo + atoms = snap.atoms + ids = {} + for i in xrange(snap.natoms): + ids[atoms[i][id]] = i + for i in xrange(snap.natoms): + j = ids[atoms[i][iother]] + atoms[i][x] += (atoms[i][ix]-atoms[j][ix])*xprd + atoms[i][y] += (atoms[i][iy]-atoms[j][iy])*yprd + atoms[i][z] += (atoms[i][iz]-atoms[j][iz])*zprd + + # -------------------------------------------------------------------- + # sort atoms by atom ID in all selected timesteps by default + # if arg = string, sort all steps by that column + # if arg = numeric, sort atoms in single step + + def sort(self,*list): + if len(list) == 0: + print "Sorting selected snapshots ..." + id = self.names["id"] + for snap in self.snaps: + if snap.tselect: self.sort_one(snap,id) + elif type(list[0]) is types.StringType: + print "Sorting selected snapshots by %s ..." % list[0] + id = self.names[list[0]] + for snap in self.snaps: + if snap.tselect: self.sort_one(snap,id) + else: + i = self.findtime(list[0]) + id = self.names["id"] + self.sort_one(self.snaps[i],id) + + # -------------------------------------------------------------------- + # sort a single snapshot by ID column + + def sort_one(self,snap,id): + atoms = snap.atoms + ids = atoms[:,id] + ordering = Numeric.argsort(ids) + for i in xrange(len(atoms[0])): + atoms[:,i] = Numeric.take(atoms[:,i],ordering) + + # -------------------------------------------------------------------- + # write a single dump file from current selection + + def write(self,file): + f = open(file,"w") + for snap in self.snaps: + if not snap.tselect: continue + print snap.time, + sys.stdout.flush() + + print >>f,"ITEM: TIMESTEP" + print >>f,snap.time + print >>f,"ITEM: NUMBER OF ATOMS" + print >>f,snap.nselect + print >>f,"ITEM: BOX BOUNDS" + print >>f,snap.xlo,snap.xhi + print >>f,snap.ylo,snap.yhi + print >>f,snap.zlo,snap.zhi + print >>f,"ITEM: ATOMS" + + atoms = snap.atoms + nvalues = len(atoms[0]) + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + line = "" + for j in xrange(nvalues): + if (j < 2): + line += str(int(atoms[i][j])) + " " + else: + line += str(atoms[i][j]) + " " + print >>f,line + f.close() + print "\n%d snapshots" % self.nselect + + # -------------------------------------------------------------------- + # write one dump file per snapshot from current selection + + def scatter(self,root): + for snap in self.snaps: + if not snap.tselect: continue + print snap.time, + sys.stdout.flush() + + file = root + "." + str(snap.time) + f = open(file,"w") + print >>f,"ITEM: TIMESTEP" + print >>f,snap.time + print >>f,"ITEM: NUMBER OF ATOMS" + print >>f,snap.nselect + print >>f,"ITEM: BOX BOUNDS" + print >>f,snap.xlo,snap.xhi + print >>f,snap.ylo,snap.yhi + print >>f,snap.zlo,snap.zhi + print >>f,"ITEM: ATOMS" + + atoms = snap.atoms + nvalues = len(atoms[0]) + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + line = "" + for j in xrange(nvalues): + if (j < 2): + line += str(int(atoms[i][j])) + " " + else: + line += str(atoms[i][j]) + " " + print >>f,line + f.close() + print "\n%d snapshots" % self.nselect + + # -------------------------------------------------------------------- + # find min/max across all selected snapshots/atoms for a particular column + + def minmax(self,colname): + icol = self.names[colname] + min = 1.0e20 + max = -min + for snap in self.snaps: + if not snap.tselect: continue + atoms = snap.atoms + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + if atoms[i][icol] < min: min = atoms[i][icol] + if atoms[i][icol] > max: max = atoms[i][icol] + return (min,max) + + # -------------------------------------------------------------------- + # set a column value via an equation + + def set(self,eq): + print "Setting ..." + pattern = "\$\w*" + list = re.findall(pattern,eq) + + lhs = list[0][1:] + if not self.names.has_key(lhs): + self.newcolumn(lhs) + + for item in list: + name = item[1:] + column = self.names[name] + insert = "snap.atoms[i][%d]" % (column) + eq = eq.replace(item,insert) + ceq = compile(eq,'','single') + + for snap in self.snaps: + if not snap.tselect: continue + for i in xrange(snap.natoms): + if snap.aselect[i]: exec ceq + + # -------------------------------------------------------------------- + # clone value in col across selected timesteps for atoms with same ID + + def clone(self,nstep,col): + istep = self.findtime(nstep) + icol = self.names[col] + id = self.names["id"] + ids = {} + for i in xrange(self.snaps[istep].natoms): + ids[self.snaps[istep].atoms[i][id]] = i + for snap in self.snaps: + if not snap.tselect: continue + atoms = snap.atoms + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + j = ids[atoms[i][id]] + atoms[i][icol] = self.snaps[istep].atoms[j][icol] + + # -------------------------------------------------------------------- + # values in old column are spread as ints from 1-N and assigned to new column + + def spread(self,old,n,new): + iold = self.names[old] + if not self.names.has_key(new): self.newcolumn(new) + inew = self.names[new] + + min,max = self.minmax(old) + print "min/max = ",min,max + + gap = max - min + invdelta = n/gap + for snap in self.snaps: + if not snap.tselect: continue + atoms = snap.atoms + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + ivalue = int((atoms[i][iold] - min) * invdelta) + 1 + if ivalue > n: ivalue = n + if ivalue < 1: ivalue = 1 + atoms[i][inew] = ivalue + + # -------------------------------------------------------------------- + # return vector of selected snapshot time stamps + + def time(self): + vec = self.nselect * [0] + i = 0 + for snap in self.snaps: + if not snap.tselect: continue + vec[i] = snap.time + i += 1 + return vec + + # -------------------------------------------------------------------- + # extract vector(s) of values for atom ID n at each selected timestep + + def atom(self,n,*list): + if len(list) == 0: + raise StandardError, "no columns specified" + columns = [] + values = [] + for name in list: + columns.append(self.names[name]) + values.append(self.nselect * [0]) + ncol = len(columns) + + id = self.names["id"] + m = 0 + for snap in self.snaps: + if not snap.tselect: continue + atoms = snap.atoms + for i in xrange(snap.natoms): + if atoms[i][id] == n: break + if atoms[i][id] != n: + raise StandardError, "could not find atom ID in snapshot" + for j in xrange(ncol): + values[j][m] = atoms[i][columns[j]] + m += 1 + + if len(list) == 1: return values[0] + else: return values + + # -------------------------------------------------------------------- + # extract vector(s) of values for selected atoms at chosen timestep + + def vecs(self,n,*list): + snap = self.snaps[self.findtime(n)] + + if len(list) == 0: + raise StandardError, "no columns specified" + columns = [] + values = [] + for name in list: + columns.append(self.names[name]) + values.append(snap.nselect * [0]) + ncol = len(columns) + + m = 0 + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + for j in xrange(ncol): + values[j][m] = snap.atoms[i][columns[j]] + m += 1 + + if len(list) == 1: return values[0] + else: return values + + # -------------------------------------------------------------------- + # add a new column to every snapshot and set value to 0 + # set the name of the column to str + + def newcolumn(self,str): + ncol = len(self.snaps[0].atoms[0]) + self.map(ncol+1,str) + for snap in self.snaps: + atoms = snap.atoms + newatoms = Numeric.zeros((snap.natoms,ncol+1),Numeric.Float) + newatoms[:,0:ncol] = snap.atoms + snap.atoms = newatoms + + # -------------------------------------------------------------------- + # sort snapshots on time stamp + + def compare_time(self,a,b): + if a.time < b.time: + return -1 + elif a.time > b.time: + return 1 + else: + return 0 + + # -------------------------------------------------------------------- + # delete successive snapshots with duplicate time stamp + + def cull(self): + i = 1 + while i < len(self.snaps): + if self.snaps[i].time == self.snaps[i-1].time: + del self.snaps[i] + else: + i += 1 + + # -------------------------------------------------------------------- + # iterate over selected snapshots + + def iterator(self,flag): + start = 0 + if flag: start = self.iterate + 1 + for i in xrange(start,self.nsnaps): + if self.snaps[i].tselect: + self.iterate = i + return i,self.snaps[i].time,1 + return 0,0,-1 + + # -------------------------------------------------------------------- + # return list of atoms to viz for snapshot isnap + # augment with bonds, tris, lines if extra() was invoked + + def viz(self,isnap): + snap = self.snaps[isnap] + + time = snap.time + box = [snap.xlo,snap.ylo,snap.zlo,snap.xhi,snap.yhi,snap.zhi] + id = self.names["id"] + type = self.names[self.atype] + x = self.names["x"] + y = self.names["y"] + z = self.names["z"] + + # create atom list needed by viz from id,type,x,y,z + # need Numeric mode here + + atoms = [] + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + atom = snap.atoms[i] + atoms.append([atom[id],atom[type],atom[x],atom[y],atom[z]]) + + # create list of current bond coords from static bondlist + # alist = dictionary of atom IDs for atoms list + # lookup bond atom IDs in alist and grab their coords + # try is used since some atoms may be unselected + # any bond with unselected atom is not returned to viz caller + # need Numeric mode here + + bonds = [] + if self.bondflag: + alist = {} + for i in xrange(len(atoms)): alist[int(atoms[i][0])] = i + for bond in self.bondlist: + try: + i = alist[bond[2]] + j = alist[bond[3]] + atom1 = atoms[i] + atom2 = atoms[j] + bonds.append([bond[0],bond[1],atom1[2],atom1[3],atom1[4], + atom2[2],atom2[3],atom2[4],atom1[1],atom2[1]]) + except: continue + + tris = [] + if self.triflag: + if self.triflag == 1: tris = self.trilist + elif self.triflag == 2: + timetmp,boxtmp,atomstmp,bondstmp, \ + tris,linestmp = self.triobj.viz(time,1) + + lines = [] + if self.lineflag: lines = self.linelist + + return time,box,atoms,bonds,tris,lines + + # -------------------------------------------------------------------- + + def findtime(self,n): + for i in xrange(self.nsnaps): + if self.snaps[i].time == n: return i + raise StandardError, "no step %d exists" % n + + # -------------------------------------------------------------------- + # return maximum box size across all selected snapshots + + def maxbox(self): + xlo = ylo = zlo = None + xhi = yhi = zhi = None + for snap in self.snaps: + if not snap.tselect: continue + if xlo == None or snap.xlo < xlo: xlo = snap.xlo + if xhi == None or snap.xhi > xhi: xhi = snap.xhi + if ylo == None or snap.ylo < ylo: ylo = snap.ylo + if yhi == None or snap.yhi > yhi: yhi = snap.yhi + if zlo == None or snap.zlo < zlo: zlo = snap.zlo + if zhi == None or snap.zhi > zhi: zhi = snap.zhi + return [xlo,ylo,zlo,xhi,yhi,zhi] + + # -------------------------------------------------------------------- + # return maximum atom type across all selected snapshots and atoms + + def maxtype(self): + icol = self.names["type"] + max = 0 + for snap in self.snaps: + if not snap.tselect: continue + atoms = snap.atoms + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + if atoms[i][icol] > max: max = atoms[i][icol] + return int(max) + + # -------------------------------------------------------------------- + # grab bonds/tris/lines from another object + + def extra(self,arg): + + # read bonds from bond dump file + + if type(arg) is types.StringType: + try: + f = open(arg,'r') + + item = f.readline() + time = int(f.readline()) + item = f.readline() + nbonds = int(f.readline()) + item = f.readline() + if not re.search("BONDS",item): + raise StandardError, "could not read bonds from dump file" + + words = f.readline().split() + ncol = len(words) + for i in xrange(1,nbonds): + words += f.readline().split() + f.close() + + # convert values to int and absolute value since can be negative types + + bondlist = Numeric.zeros((nbonds,4),Numeric.Int) + ints = [abs(int(value)) for value in words] + start = 0 + stop = 4 + for i in xrange(nbonds): + bondlist[i] = ints[start:stop] + start += ncol + stop += ncol + if bondlist: + self.bondflag = 1 + self.bondlist = bondlist + except: + raise StandardError,"could not read from bond dump file" + + # request bonds from data object + + elif type(arg) is types.InstanceType and ".data" in str(arg.__class__): + try: + bondlist = [] + bondlines = arg.sections["Bonds"] + for line in bondlines: + words = line.split() + bondlist.append([int(words[0]),int(words[1]), + int(words[2]),int(words[3])]) + if bondlist: + self.bondflag = 1 + self.bondlist = bondlist + except: + raise StandardError,"could not extract bonds from data object" + + # request tris/lines from cdata object + + elif type(arg) is types.InstanceType and ".cdata" in str(arg.__class__): + try: + tmp,tmp,tmp,tmp,tris,lines = arg.viz(0) + if tris: + self.triflag = 1 + self.trilist = tris + if lines: + self.lineflag = 1 + self.linelist = lines + except: + raise StandardError,"could not extract tris/lines from cdata object" + + # request tris from mdump object + + elif type(arg) is types.InstanceType and ".mdump" in str(arg.__class__): + try: + self.triflag = 2 + self.triobj = arg + except: + raise StandardError,"could not extract tris from mdump object" + + else: + raise StandardError,"unrecognized argument to dump.extra()" + + # -------------------------------------------------------------------- + + def compare_atom(self,a,b): + if a[0] < b[0]: + return -1 + elif a[0] > b[0]: + return 1 + else: + return 0 + +# -------------------------------------------------------------------- +# one snapshot + +class Snap: + pass + +# -------------------------------------------------------------------- +# time selection class + +class tselect: + + def __init__(self,data): + self.data = data + + # -------------------------------------------------------------------- + + def all(self): + data = self.data + for snap in data.snaps: + snap.tselect = 1 + data.nselect = len(data.snaps) + data.aselect.all() + print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + + # -------------------------------------------------------------------- + + def one(self,n): + data = self.data + for snap in data.snaps: + snap.tselect = 0 + i = data.findtime(n) + data.snaps[i].tselect = 1 + data.nselect = 1 + data.aselect.all() + print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + + # -------------------------------------------------------------------- + + def none(self): + data = self.data + for snap in data.snaps: + snap.tselect = 0 + data.nselect = 0 + print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + + # -------------------------------------------------------------------- + + def skip(self,n): + data = self.data + count = n-1 + for snap in data.snaps: + if not snap.tselect: continue + count += 1 + if count == n: + count = 0 + continue + snap.tselect = 0 + data.nselect -= 1 + data.aselect.all() + print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + + # -------------------------------------------------------------------- + + def test(self,teststr): + data = self.data + snaps = data.snaps + cmd = "flag = " + teststr.replace("$t","snaps[i].time") + ccmd = compile(cmd,'','single') + for i in xrange(data.nsnaps): + if not snaps[i].tselect: continue + exec ccmd + if not flag: + snaps[i].tselect = 0 + data.nselect -= 1 + data.aselect.all() + print "%d snapshots selected out of %d" % (data.nselect,data.nsnaps) + +# -------------------------------------------------------------------- +# atom selection class + +class aselect: + + def __init__(self,data): + self.data = data + + # -------------------------------------------------------------------- + + def all(self,*args): + data = self.data + if len(args) == 0: # all selected timesteps + for snap in data.snaps: + if not snap.tselect: continue + for i in xrange(snap.natoms): snap.aselect[i] = 1 + snap.nselect = snap.natoms + else: # one timestep + n = data.findtime(args[0]) + snap = data.snaps[n] + for i in xrange(snap.natoms): snap.aselect[i] = 1 + snap.nselect = snap.natoms + + # -------------------------------------------------------------------- + + def test(self,teststr,*args): + data = self.data + + # replace all $var with snap.atoms references and compile test string + + pattern = "\$\w*" + list = re.findall(pattern,teststr) + for item in list: + name = item[1:] + column = data.names[name] + insert = "snap.atoms[i][%d]" % column + teststr = teststr.replace(item,insert) + cmd = "flag = " + teststr + ccmd = compile(cmd,'','single') + + if len(args) == 0: # all selected timesteps + for snap in data.snaps: + if not snap.tselect: continue + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + exec ccmd + if not flag: + snap.aselect[i] = 0 + snap.nselect -= 1 + for i in xrange(data.nsnaps): + if data.snaps[i].tselect: + print "%d atoms of %d selected in first step %d" % \ + (data.snaps[i].nselect,data.snaps[i].natoms,data.snaps[i].time) + break + for i in xrange(data.nsnaps-1,-1,-1): + if data.snaps[i].tselect: + print "%d atoms of %d selected in last step %d" % \ + (data.snaps[i].nselect,data.snaps[i].natoms,data.snaps[i].time) + break + + else: # one timestep + n = data.findtime(args[0]) + snap = data.snaps[n] + for i in xrange(snap.natoms): + if not snap.aselect[i]: continue + exec ccmd + if not flag: + snap.aselect[i] = 0 + snap.nselect -= 1 diff -Naur lammps-1Apr08/tools/python/pizza/gnu.py lammps-2Apr08/tools/python/pizza/gnu.py --- lammps-1Apr08/tools/python/pizza/gnu.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/gnu.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,383 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# gnu tool + +oneline = "Create plots via GnuPlot plotting program" + +docstr = """ +g = gnu() start up GnuPlot +g.stop() shut down GnuPlot process + +g.plot(a) plot vector A against linear index +g.plot(a,b) plot B against A +g.plot(a,b,c,d,...) plot B against A, D against C, etc +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 vector + mplot loops over range(M,N,S) and create one plot per iteration + last args are same as list of vectors for plot(), e.g. 1, 2, 4 vectors + each plot is made from a portion of the vectors, depending on loop index i + Ith plot is of b[0:i] vs a[0:i], etc + series of plots saved as file0000.eps, file0001.eps, etc + if use xrange(),yrange() then plot axes will be same for all plots + +g("plot 'file.dat' using 2:3 with lines") execute string in GnuPlot + +g.enter() enter GnuPlot shell +gnuplot> plot sin(x) with lines type commands directly to GnuPlot +gnuplot> exit, quit exit GnuPlot shell + +g.export("data",range(100),a,...) create file with columns of numbers + + all vectors must be of equal length + could plot from file with GnuPlot command: plot 'data' using 1:2 with lines + +g.select(N) figure N becomes the current plot + + subsequent commands apply to this plot + +g.hide(N) delete window for figure N +g.save("file") save current plot as file.eps + +Set attributes for current plot: + +g.erase() reset all attributes to default values +g.aspect(1.3) aspect ratio +g.xtitle("Time") x axis text +g.ytitle("Energy") y axis text +g.title("My Plot") title text +g.title("title","x","y") title, x axis, y axis text +g.xrange(xmin,xmax) x axis range +g.xrange() default x axis range +g.yrange(ymin,ymax) y axis range +g.yrange() default y axis range +g.xlog() toggle x axis between linear and log +g.ylog() toggle y axis between linear and log +g.label(x,y,"text") place label at x,y coords +g.curve(N,'r') set color of curve N + + colors: 'k' = black, 'r' = red, 'g' = green, 'b' = blue + 'm' = magenta, 'c' = cyan, 'y' = yellow +""" + +# History +# 8/05, Matt Jones (BYU): original version +# 9/05, Steve Plimpton: added mplot() method + +# ToDo list +# allow choice of JPG or PNG or GIF when saving ? +# can this be done from GnuPlot or have to do via ImageMagick convert ? +# way to trim EPS plot that is created ? +# hide does not work on Mac aqua +# select does not pop window to front on Mac aqua + +# Variables +# current = index of current figure (1-N) +# figures = list of figure objects with each plot's attributes +# so they aren't lost between replots + +# Imports and external programs + +import types, os + +try: from DEFAULTS import PIZZA_GNUPLOT +except: PIZZA_GNUPLOT = "gnuplot" +try: from DEFAULTS import PIZZA_GNUTERM +except: PIZZA_GNUTERM = "x11" + +# Class definition + +class gnu: + + # -------------------------------------------------------------------- + + def __init__(self): + self.GNUPLOT = os.popen(PIZZA_GNUPLOT,'w') + self.file = "tmp.gnu" + self.figures = [] + self.select(1) + + # -------------------------------------------------------------------- + + def stop(self): + self.__call__("quit") + del self.GNUPLOT + + # -------------------------------------------------------------------- + + def __call__(self,command): + self.GNUPLOT.write(command + '\n') + self.GNUPLOT.flush() + + # -------------------------------------------------------------------- + + def enter(self): + while 1: + command = raw_input("gnuplot> ") + if command == "quit" or command == "exit": return + self.__call__(command) + + # -------------------------------------------------------------------- + # write plot vectors to files and plot them + + def plot(self,*vectors): + if len(vectors) == 1: + file = self.file + ".%d.1" % self.current + linear = range(len(vectors[0])) + self.export(file,linear,vectors[0]) + self.figures[self.current-1].ncurves = 1 + else: + if len(vectors) % 2: raise StandardError,"vectors must come in pairs" + for i in range(0,len(vectors),2): + file = self.file + ".%d.%d" % (self.current,i/2+1) + self.export(file,vectors[i],vectors[i+1]) + 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): + partial_vecs = [] + for vec in vectors: + 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) + else: newfile = file + str(n) + + self.save(newfile) + n += 1 + + # -------------------------------------------------------------------- + # write list of equal-length vectors to filename + + def export(self,filename,*vectors): + n = len(vectors[0]) + for vector in vectors: + if len(vector) != n: raise StandardError,"vectors must be same length" + f = open(filename,'w') + nvec = len(vectors) + for i in xrange(n): + for j in xrange(nvec): + print >>f,vectors[j][i], + print >>f + f.close() + + # -------------------------------------------------------------------- + # select plot N as current plot + + def select(self,n): + self.current = n + if len(self.figures) < n: + for i in range(n - len(self.figures)): + self.figures.append(figure()) + cmd = "set term " + PIZZA_GNUTERM + ' ' + str(n) + self.__call__(cmd) + if self.figures[n-1].ncurves: self.draw() + + # -------------------------------------------------------------------- + # delete window for plot N + + def hide(self,n): + cmd = "set term %s close %d" % (PIZZA_GNUTERM,n) + self.__call__(cmd) + + # -------------------------------------------------------------------- + # save plot to file.eps + # final re-select will reset terminal + # 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 + self.__call__(cmd) + if os.path.exists("tmp.done"): os.remove("tmp.done") + self.draw() + self.__call__("!touch tmp.done") + 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 + + def erase(self): + fig = figure() + fig.ncurves = self.figures[self.current-1].ncurves + self.figures[self.current-1] = fig + self.draw() + + # -------------------------------------------------------------------- + + def aspect(self,value): + self.figures[self.current-1].aspect = value + self.draw() + + # -------------------------------------------------------------------- + + def xrange(self,*values): + if len(values) == 0: + self.figures[self.current-1].xlimit = 0 + else: + self.figures[self.current-1].xlimit = (values[0],values[1]) + self.draw() + + # -------------------------------------------------------------------- + + def yrange(self,*values): + if len(values) == 0: + self.figures[self.current-1].ylimit = 0 + 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.draw() + + # -------------------------------------------------------------------- + + def nolabels(self): + self.figures[self.current-1].nlabel = 0 + self.figures[self.current-1].labels = [] + self.draw() + + # -------------------------------------------------------------------- + + def title(self,*strings): + if len(strings) == 1: + self.figures[self.current-1].title = strings[0] + else: + self.figures[self.current-1].title = strings[0] + self.figures[self.current-1].xtitle = strings[1] + self.figures[self.current-1].ytitle = strings[2] + self.draw() + + # -------------------------------------------------------------------- + + 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): + if self.figures[self.current-1].xlog: + self.figures[self.current-1].xlog = 0 + else: + self.figures[self.current-1].xlog = 1 + self.draw() + + # -------------------------------------------------------------------- + + def ylog(self): + if self.figures[self.current-1].ylog: + self.figures[self.current-1].ylog = 0 + else: + self.figures[self.current-1].ylog = 1 + self.draw() + + # -------------------------------------------------------------------- + + def curve(self,num,color): + fig = self.figures[self.current-1] + while len(fig.colors) < num: fig.colors.append(0) + fig.colors[num-1] = colormap[color] + self.draw() + + # -------------------------------------------------------------------- + # draw a plot with all its settings + # just return if no files of vectors defined yet + + 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 + '"' + self.__call__(cmd) + cmd = 'set ylabel ' + '"' + fig.ytitle + '"' + self.__call__(cmd) + + if fig.xlog: self.__call__("set logscale x") + else: self.__call__("unset logscale x") + if fig.ylog: self.__call__("set logscale y") + else: self.__call__("unset logscale y") + 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: + cmd = 'set yr [' + str(fig.ylimit[0]) + ':' + str(fig.ylimit[1]) + ']' + self.__call__(cmd) + else: self.__call__("set yr [*:*]") + + self.__call__("set nolabel") + for i in range(fig.nlabels): + x = fig.labels[i][0] + y = fig.labels[i][1] + text = fig.labels[i][2] + cmd = 'set label ' + '\"' + text + '\" at ' + str(x) + ',' + str(y) + self.__call__(cmd) + + self.__call__("set key off") + cmd = 'plot ' + for i in range(fig.ncurves): + file = self.file + ".%d.%d" % (self.current,i+1) + if len(fig.colors) > i and fig.colors[i]: + cmd += "'" + file + "' using 1:2 with line %d, " % fig.colors[i] + else: + cmd += "'" + file + "' using 1:2 with lines, " + self.__call__(cmd[:-2]) + +# -------------------------------------------------------------------- +# class to store settings for a single plot + +class figure: + + def __init__(self): + self.ncurves = 0 + self.colors = [] + self.title = "" + self.xtitle = "" + self.ytitle = "" + self.aspect = 1.3 + self.xlimit = 0 + self.ylimit = 0 + self.xlog = 0 + self.ylog = 0 + self.nlabels = 0 + self.labels = [] + +# -------------------------------------------------------------------- +# line color settings + +colormap = {'k':-1, 'r':1, 'g':2, 'b':3, 'm':4, 'c':5, 'y':7} diff -Naur lammps-1Apr08/tools/python/pizza/log.py lammps-2Apr08/tools/python/pizza/log.py --- lammps-1Apr08/tools/python/pizza/log.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/log.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,334 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# log tool + +oneline = "Read LAMMPS log files and extract thermodynamic data" + +docstr = """ +l = log("file1") read in one or more log files +l = log("log1 log2.gz") can be gzipped +l = log("file*") wildcard expands to multiple files +l = log("log.lammps",0) two args = store filename, but don't read + + incomplete and duplicate thermo entries are deleted + +time = l.next() read new thermo info from file + + used with 2-argument constructor to allow reading thermo incrementally + return time stamp of last thermo read + return -1 if no new thermo since last read + +nvec = l.nvec # of vectors of thermo info +nlen = l.nlen length of each vectors +names = l.names list of vector names +t,pe,... = l.get("Time","KE",...) return one or more vectors of values +l.write("file.txt") write all vectors to a file +l.write("file.txt","Time","PE",...) write listed vectors to a file + + get and write allow abbreviated (uniquely) vector names +""" + +# History +# 8/05, Steve Plimpton (SNL): original version + +# ToDo list + +# Variables +# nvec = # of vectors +# nlen = length of each vector +# names = list of vector names +# ptr = dictionary, key = name, value = index into data for which column +# data[i][j] = 2d array of floats, i = 0 to # of entries, j = 0 to nvecs-1 +# style = style of LAMMPS log file, 1 = multi, 2 = one, 3 = gran +# firststr = string that begins a thermo section in log file +# increment = 1 if log file being read incrementally +# eof = ptr into incremental file for where to start next read + +# Imports and external programs + +import sys, re, glob +from os import popen + +try: tmp = PIZZA_GUNZIP +except: PIZZA_GUNZIP = "gunzip" + +# Class definition + +class log: + + # -------------------------------------------------------------------- + + def __init__(self,*list): + self.nvec = 0 + self.names = [] + self.ptr = {} + self.data = [] + + # flist = list of all log file names + + words = list[0].split() + self.flist = [] + for word in words: self.flist += glob.glob(word) + if len(self.flist) == 0 and len(list) == 1: + raise StandardError,"no log file specified" + + if len(list) == 1: + self.increment = 0 + self.read_all() + else: + if len(self.flist) > 1: + raise StandardError,"can only incrementally read one log file" + self.increment = 1 + self.eof = 0 + + # -------------------------------------------------------------------- + # read all thermo from all files + + def read_all(self): + self.read_header(self.flist[0]) + if self.nvec == 0: raise StandardError,"log file has no values" + + # read all files + + for file in self.flist: self.read_one(file) + print + + # sort entries by timestep, cull duplicates + + self.data.sort(self.compare) + self.cull() + self.nlen = len(self.data) + print "read %d log entries" % self.nlen + + # -------------------------------------------------------------------- + + def next(self): + if not self.increment: raise StandardError,"cannot read incrementally" + + if self.nvec == 0: + try: open(self.flist[0],'r') + except: return -1 + self.read_header(self.flist[0]) + if self.nvec == 0: return -1 + + self.eof = self.read_one(self.flist[0],self.eof) + return int(self.data[-1][0]) + + # -------------------------------------------------------------------- + + def get(self,*keys): + if len(keys) == 0: + raise StandardError, "no log vectors specified" + + map = [] + for key in keys: + if self.ptr.has_key(key): + map.append(self.ptr[key]) + else: + count = 0 + for i in range(self.nvec): + if self.names[i].find(key) == 0: + count += 1 + index = i + if count == 1: + map.append(index) + else: + raise StandardError, "unique log vector %s not found" % key + + vecs = [] + for i in range(len(keys)): + vecs.append(self.nlen * [0]) + for j in xrange(self.nlen): + vecs[i][j] = self.data[j][map[i]] + + if len(keys) == 1: return vecs[0] + else: return vecs + + # -------------------------------------------------------------------- + + def write(self,filename,*keys): + if len(keys): + map = [] + for key in keys: + if self.ptr.has_key(key): + map.append(self.ptr[key]) + else: + count = 0 + for i in range(self.nvec): + if self.names[i].find(key) == 0: + count += 1 + index = i + if count == 1: + map.append(index) + else: + raise StandardError, "unique log vector %s not found" % key + else: + map = range(self.nvec) + + f = open(filename,"w") + for i in xrange(self.nlen): + for j in xrange(len(map)): + print >>f,self.data[i][map[j]], + print >>f + f.close() + + # -------------------------------------------------------------------- + + def compare(self,a,b): + if a[0] < b[0]: + return -1 + elif a[0] > b[0]: + return 1 + else: + return 0 + + # -------------------------------------------------------------------- + + def cull(self): + i = 1 + while i < len(self.data): + if self.data[i][0] == self.data[i-1][0]: del self.data[i] + else: i += 1 + + # -------------------------------------------------------------------- + + def read_header(self,file): + str_multi = "----- Step" + str_one = "Step " + + if file[-3:] == ".gz": + txt = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r').read() + else: + txt = open(file).read() + + if txt.find(str_multi) >= 0: + self.firststr = str_multi + self.style = 1 + elif txt.find(str_one) >= 0: + self.firststr = str_one + self.style = 2 + else: + return + + if self.style == 1: + s1 = txt.find(self.firststr) + s2 = txt.find("\n--",s1) + pattern = "\s(\S*)\s*=" + keywords = re.findall(pattern,txt[s1:s2]) + keywords.insert(0,"Step") + i = 0 + for keyword in keywords: + self.names.append(keyword) + self.ptr[keyword] = i + i += 1 + + else: + s1 = txt.find(self.firststr) + s2 = txt.find("\n",s1) + line = txt[s1:s2] + words = line.split() + for i in range(len(words)): + self.names.append(words[i]) + self.ptr[words[i]] = i + + self.nvec = len(self.names) + + # -------------------------------------------------------------------- + + def read_one(self,*list): + + # if 2nd arg exists set file ptr to that value + # read entire (rest of) file into txt + + file = list[0] + if file[-3:] == ".gz": + f = popen("%s -c %s" % (PIZZA_GUNZIP,file),'r') + else: + f = open(file) + + if len(list) == 2: f.seek(list[1]) + txt = f.read() + if file[-3:] == ".gz": eof = 0 + else: eof = f.tell() + f.close() + + start = last = 0 + while not last: + + # chunk = contiguous set of thermo entries (line or multi-line) + # s1 = 1st char on 1st line of chunk + # s2 = 1st char on line after chunk + # set last = 1 if this is last chunk in file, leave 0 otherwise + # set start = position in file to start looking for next chunk + # rewind eof if final entry is incomplete + + s1 = txt.find(self.firststr,start) + s2 = txt.find("Loop time of",start+1) + + if s1 >= 0 and s2 >= 0 and s1 < s2: # found s1,s2 with s1 before s2 + if self.style == 2: + s1 = txt.find("\n",s1) + 1 + elif s1 >= 0 and s2 >= 0 and s2 < s1: # found s1,s2 with s2 before s1 + s1 = 0 + elif s1 == -1 and s2 >= 0: # found s2, but no s1 + last = 1 + s1 = 0 + elif s1 >= 0 and s2 == -1: # found s1, but no s2 + last = 1 + if self.style == 1: + s2 = txt.rfind("\n--",s1) + 1 + else: + s1 = txt.find("\n",s1) + 1 + s2 = txt.rfind("\n",s1) + 1 + eof -= len(txt) - s2 + elif s1 == -1 and s2 == -1: # found neither + # could be end-of-file section + # or entire read was one chunk + + if txt.find("Loop time of",start) == start: # end of file, so exit + eof -= len(txt) - start # reset eof to "Loop" + break + + last = 1 # entire read is a chunk + s1 = 0 + if self.style == 1: + s2 = txt.rfind("\n--",s1) + 1 + else: + s2 = txt.rfind("\n",s1) + 1 + eof -= len(txt) - s2 + if s1 == s2: break + + chunk = txt[s1:s2-1] + start = s2 + + # split chunk into entries + # parse each entry for numeric fields, append to data + + if self.style == 1: + sections = chunk.split("\n--") + pat1 = re.compile("Step\s*(\S*)\s") + pat2 = re.compile("=\s*(\S*)") + for section in sections: + word1 = [re.search(pat1,section).group(1)] + word2 = re.findall(pat2,section) + words = word1 + word2 + self.data.append(map(float,words)) + + else: + lines = chunk.split("\n") + for line in lines: + words = line.split() + self.data.append(map(float,words)) + + # print last timestep of chunk + + print int(self.data[len(self.data)-1][0]), + sys.stdout.flush() + + return eof diff -Naur lammps-1Apr08/tools/python/pizza/pdbfile.py lammps-2Apr08/tools/python/pizza/pdbfile.py --- lammps-1Apr08/tools/python/pizza/pdbfile.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/pdbfile.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,289 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# pdb tool + +oneline = "Read, write PDB files in combo with LAMMPS snapshots" + +docstr = """ +p = pdbfile("3CRO") create pdb object from PDB file or WWW +p = pdbfile("pep1 pep2") read in multiple PDB files +p = pdbfile("pep*") can use wildcards +p = pdbfile(d) read in snapshot data with no PDB file +p = pdbfile("3CRO",d) read in single PDB file with snapshot data + + string arg contains one or more PDB files + don't need .pdb suffix except wildcard must expand to file.pdb + 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, ... +p.many("mine") write as mine0000.pdb, mine0001.pdb, ... +p.single(N) write timestamp N as tmp.pdb +p.single(N,"new") write as new.pdb + + how new PDB files are created depends on constructor inputs: + if no d: one new PDB file for each file in string arg (just a copy) + if only d specified: one new PDB file per snapshot in generic format + 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) + + iterator = loop over number of PDB files + call first time with arg = 0, thereafter with arg = 1 + N = length = # of snapshots or # of input PDB files + index = index of snapshot or input PDB file (0 to N-1) + time = timestep value (time stamp for snapshot, index for multiple PDB) + flag = -1 when iteration is done, 1 otherwise + typically call p.single(time) in iterated loop to write out one PDB file +""" + +# History +# 8/05, Steve Plimpton (SNL): original version + +# ToDo list +# for generic PDB file (no template) from a LJ unit system, +# the atoms in PDB file are too close together + +# Variables +# files = list of input PDB files +# data = data object (ccell,data,dump) to read snapshots from +# atomlines = dict of ATOM lines in original PDB file +# key = atom id, value = tuple of (beginning,end) of line + +# Imports and external programs + +import sys, types, glob, urllib + +# Class definition + +class pdbfile: + + # -------------------------------------------------------------------- + + def __init__(self,*args): + if len(args) == 1: + if type(args[0]) is types.StringType: + filestr = args[0] + self.data = None + else: + filestr = None + self.data = args[0] + elif len(args) == 2: + filestr = args[0] + self.data = args[1] + else: raise StandardError, "invalid args for pdb()" + + # flist = full list of all PDB input file names + # append .pdb if needed + + if filestr: + list = filestr.split() + flist = [] + for file in list: + if '*' in file: flist += glob.glob(file) + else: flist.append(file) + for i in xrange(len(flist)): + if flist[i][-4:] != ".pdb": flist[i] += ".pdb" + if len(flist) == 0: + raise StandardError,"no PDB file specified" + self.files = flist + else: self.files = [] + + if len(self.files) > 1 and self.data: + raise StandardError, "cannot use multiple PDB files with data object" + if len(self.files) == 0 and not self.data: + 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() + except: + print "downloading %s from http://rcsb.org" % self.files[0] + fetchstr = "http://www.rcsb.org/pdb/cgi/export.cgi/%s?format=PDB&pdbId=2cpk&compression=None" % self.files[0] + urllib.urlretrieve(fetchstr,self.files[0]) + + 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: + # only selected atoms returned by extract + # atoms written in order they appear in snapshot + # atom only written if its tag is in PDB template file + # if no data: + # concatenate all input files to one output file + + def one(self,*args): + if len(args) == 0: file = "tmp.pdb" + elif args[0][-4:] == ".pdb": file = args[0] + else: file = args[0] + ".pdb" + + f = open(file,'w') + + # use template PDB file with each snapshot + + if self.data: + n = flag = 0 + while 1: + which,time,flag = self.data.iterator(flag) + if flag == -1: break + self.convert(f,which) + print >>f,"END" + print time, + sys.stdout.flush() + n += 1 + + else: + for file in self.files: + f.write(open(file,'r').read()) + print >>f,"END" + print file, + sys.stdout.flush() + + f.close() + print "\nwrote %d datasets to %s in PDB format" % (n,file) + + # -------------------------------------------------------------------- + # write series of numbered PDB files + # if data exists: + # only selected atoms returned by extract + # atoms written in order they appear in snapshot + # atom only written if its tag is in PDB template file + # if no data: + # just copy all input files to output files + + def many(self,*args): + if len(args) == 0: root = "tmp" + else: root = args[0] + + if self.data: + n = flag = 0 + while 1: + which,time,flag = self.data.iterator(flag) + if flag == -1: break + + if n < 10: + file = root + "000" + str(n) + elif n < 100: + file = root + "00" + str(n) + elif n < 1000: + file = root + "0" + str(n) + else: + file = root + str(n) + file += ".pdb" + + f = open(file,'w') + self.convert(f,which) + f.close() + + print time, + sys.stdout.flush() + n += 1 + + else: + n = 0 + for infile in self.files: + if n < 10: + file = root + "000" + str(n) + elif n < 100: + file = root + "00" + str(n) + elif n < 1000: + file = root + "0" + str(n) + 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) + + # -------------------------------------------------------------------- + # write a single PDB file + # if data exists: + # time is timestamp in snapshot + # only selected atoms returned by extract + # atoms written in order they appear in snapshot + # atom only written if its tag is in PDB template file + # if no data: + # time is index into list of input PDB files + # just copy one input file to output file + + def single(self,time,*args): + if len(args) == 0: file = "tmp.pdb" + elif args[0][-4:] == ".pdb": file = args[0] + else: file = args[0] + ".pdb" + f = open(file,'w') + + if self.data: + which = self.data.findtime(time) + self.convert(f,which) + else: + f.write(open(self.files[time],'r').read()) + + f.close() + + # -------------------------------------------------------------------- + # iterate over list of input files or selected snapshots + # latter is done via data objects iterator + + def iterator(self,flag): + if not self.data: + if not flag: self.iterate = 0 + else: + self.iterate += 1 + if self.iterate > len(self.files): return 0,0,-1 + return self.iterate,self.iterate,1 + + return self.data.iterator(flag) + + # -------------------------------------------------------------------- + # read a PDB file and store ATOM lines + + def read_template(self,file): + lines = open(file,'r').readlines() + self.atomlines = {} + for line in lines: + if line.find("ATOM") == 0: + tag = int(line[4:11]) + begin = line[:30] + end = line[54:] + self.atomlines[tag] = (begin,end) + + # -------------------------------------------------------------------- + # convert one set of atoms to PDB format and write to f + + def convert(self,f,which): + time,box,atoms,bonds,tris,lines = self.data.viz(which) + if len(self.files): + for atom in atoms: + id = atom[0] + if self.atomlines.has_key(id): + (begin,end) = self.atomlines[id] + line = "%s%8.3f%8.3f%8.3f%s" % (begin,atom[2],atom[3],atom[4],end) + print >>f,line, + else: + for atom in atoms: + begin = "ATOM %6d %2d R00 1 " % (atom[0],atom[1]) + middle = "%8.3f%8.3f%8.3f" % (atom[2],atom[3],atom[4]) + end = " 1.00 0.00 NONE" + print >>f,begin+middle+end diff -Naur lammps-1Apr08/tools/python/pizza/xyz.py lammps-2Apr08/tools/python/pizza/xyz.py --- lammps-1Apr08/tools/python/pizza/xyz.py 1969-12-31 17:00:00.000000000 -0700 +++ lammps-2Apr08/tools/python/pizza/xyz.py 2008-03-31 16:28:11.000000000 -0600 @@ -0,0 +1,121 @@ +# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html +# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories +# +# Copyright (2005) Sandia Corporation. Under the terms of Contract +# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains +# certain rights in this software. This software is distributed under +# the GNU General Public License. + +# xyz tool + +oneline = "Convert LAMMPS snapshots to XYZ format" + +docstr = """ +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 +""" + +# History +# 8/05, Steve Plimpton (SNL): original version + +# ToDo list + +# Variables +# data = data file to read from + +# Imports and external programs + +import sys + +# Class definition + +class xyz: + + # -------------------------------------------------------------------- + + def __init__(self,data): + self.data = data + + # -------------------------------------------------------------------- + + def one(self,*args): + if len(args) == 0: file = "tmp.xyz" + elif args[0][-4:] == ".xyz": file = args[0] + else: file = args[0] + ".xyz" + + f = open(file,"w") + n = flag = 0 + while 1: + which,time,flag = self.data.iterator(flag) + if flag == -1: break + time,box,atoms,bonds,tris,lines = self.data.viz(which) + + print >>f,len(atoms) + print >>f,"Atoms" + 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): + if len(args) == 0: root = "tmp" + else: root = args[0] + + n = flag = 0 + while 1: + 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: + file = root + "00" + str(n) + elif n < 1000: + file = root + "0" + str(n) + else: + file = root + str(n) + file += ".xyz" + f = open(file,"w") + print >>f,len(atoms) + print >>f,"Atoms" + for atom in atoms: + itype = int(atom[1]) + print >>f,itype,atom[2],atom[3],atom[4] + print time, + sys.stdout.flush() + f.close() + n += 1 + + print "\nwrote %s snapshots in XYZ format" % n + + # -------------------------------------------------------------------- + + def single(self,time,*args): + if len(args) == 0: file = "tmp.xyz" + elif args[0][-4:] == ".xyz": file = args[0] + else: file = args[0] + ".xyz" + + which = self.data.findtime(time) + time,box,atoms,bonds,tris,lines = self.data.viz(which) + f = open(file,"w") + print >>f,len(atoms) + print >>f,"Atoms" + for atom in atoms: + itype = int(atom[1]) + print >>f,itype,atom[2],atom[3],atom[4] + f.close()