forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1691 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
8f8456fd3a
commit
93ed818975
|
|
@ -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
|
||||
|
|
@ -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)
|
||||
|
|
@ -0,0 +1,32 @@
|
|||
#!/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
|
||||
# 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
|
||||
# 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:
|
||||
raise StandardError, "Syntax: dump2pdb.py dumpfile Nid Ntype Nx Ny Nz pdbfile"
|
||||
|
||||
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]
|
||||
|
||||
d = dump(dumpfile)
|
||||
d.map(nid,"id",ntype,"type",nx,"x",ny,"y",nz,"z")
|
||||
p = pdbfile(d)
|
||||
p.one(pfile)
|
||||
|
|
@ -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)
|
||||
|
|
@ -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)
|
||||
|
|
@ -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)
|
||||
|
|
@ -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"
|
||||
|
|
@ -0,0 +1,185 @@
|
|||
# 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)
|
||||
|
||||
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)
|
||||
|
||||
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()
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -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}
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -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()
|
||||
Loading…
Reference in New Issue